Advection Diffusion Problems

Numerical Methods for Differential Equations

Advection-Diffusion

given vector field

for our case we go to 1D case

diffusive term advection term

b.c: u(0) = 0

u(L) = 1

• We want to figure out an exact solution of the problem.

U(x) = _e^β/μx⁄^eβ/μL − 1

β/μ >> 1

⇒

u(x) ≈ _{e^β/μ(x−L)}⁄^eβ/μL

β/μ = 10

• linear sol. when diffusive term dominates with respect to transport
• when transport dominates with respect to diffusion solution in 2 regions, 1 with strong variations, the other almost constant → high gradient regions = boundary layer → amplitude is of the size of: μ/β

Discretization with Finite Elements

• Dirichlet b.c.

global becomes

1^st and last row adjusted according to b.c.

∀ node:

difference as finite differences

exact solution?

Pecelet number

0

NUMERICAL METHODS FOR DIFFERENTIAL EQNS

ADVECTION-DIFFUSION

given vector field

for our ease we go to 1D case

diffusive term

advection term

• We want to figure out an exact solution of the problem:

• linear sol. when diffusive term dominates with respect to transport

• when transport dominates with respect to diffusion solution with 2 regions: 1 with strong variations, the other almost constant -> high gradient regions = boundary layer

Discretization with Finite Elements

Dirichlet b.c.

could also have constant term

global becomes

1st and last row adjusted according to b.c.

e.g. discretization of transport term: (3rd line of the problem)

Classic 1D F.E. scheme h, equivalent to a centered finite difference scheme

Laplace _i.e. discretization

Central approximation

∀ node:

Finite differences

Peclet number

exact solution?

Pe = _Bh/^2μ

similar meaning of Pe, tells relative importance of advection and diffusion terms

• u = g _i^j, So fundam. solution a substitution we get charact. eqn a b.c. ⇨ Solution

Solution: u _i = _1-^{(1+Pe)/(1-Pe)i}/^{(1+Pe)/(1-Pe)n}

Suppose, for simplicity, n=even, Pe ≈ 1 ⇨ Bh ≈ 2μ (Adv. dominates wrt Diff.)

Where is the problem?

• in the discretization we say something not meaningful from physical point of view!
• u ⁿ⁺¹ - u ^n-1/₂ is sensitive to the direction!!

We don't want to change approximation, only to avoid oscillations work on Pe

Pe = _Bh/^2μ ⇨ can play with h, so that Pe < 1

• h < 2μ/β
• in practice increasing of number of elements not used
• It involves very small elements, big issue in 2D or 3D

Other possibility: Loose something in order of convergence to avoid oscillation

• Upwind finite difference approximation
• only 1^st order accurate
• it adds artificial diffusion, "killing" oscillations

u = _{u ^i-1 + 2u ⁱ - u ⁱ⁺¹}/^h2

• [ β u ⁱ⁺¹ - u ⁱ/_h = 0

u _i - u _i+1 = -u _i-1 + 2u _i - u _i+1/

It is, in practice, another diffusive term centered approximation of 1^st derivative

( μ + _Bh/²) _{u i - u i+1 + 2u i - u i+1}/^h2

equivalent to the original problem, in which physical diffusivity is increased of an artificial diffusion term, τ = _3Bh/²

We can try to generalize the hypothesis:

μ ( 1 + _Bh 1 ) x = μ (1 + Pe) = x

P_e = _Pe B_h 2/h² = P₂ 1 + Pe

< < 1 ALWAYS!Very good for stability

μ = μ (1 + (ϕ(ρe)) ϕ (t) such that lim ϕ (t) = 0t→+∞If ϕ(t) = 0 go back to the original problem ( centered approx 2^nd oder )ϕ(t) = upwinding, 1^st oderϕ(t) = t - t + B(t), B(t) = _tt, t ₁t0 Bernoulli function

Schapfelter - Gummel (SG) 2^nd oderfor constant facing term solution in the nodes analytical (TO)

Original problem in multi-dimensional case:

-μ Δu + (β ∙ ∇) u = 0

discretization parameter, mesh size

(μ - ch(κ)) Δu + (β ∙ ∇) u = 0

stabilization parameter, can be properly chosen

Actually we're adding an isotropic viscous term, same in all directions, butvector fields β affects the result

∫ [(μ + ch ∇) [ (∇u)∇v]] dz +∫ [β ∙ ∇u ∇v] dz = 0

Could be better not to have anisotropic diffusion coefficient (acting more on vector field dir)

-div (∇ ( k ∇v) ), k_ij = ch β_i β_j weakfoundation

∫_Ωch (β ∙ ∇v) (β ∙ ∇v)

projecting viscosity on the direction of vector field β

μ ^d2dy 2 + βd^dydy dx = 0∫ μd^dvdv dx + ∫ βd^dvdv dx = 0add μ, so that∫_Ω (1 + _Bh₂) d^dvdv dx + (2/h) ∫ β d^dvdv dx = 0

μd²d^givend^by x^idxi

Like we have changed test function for diffusivity term = χ

∫ μd^dvdu dx = 0

Test function is giving more weight to upwinding terms and less to downwinding ones

Upwinding done by χ, no more by hand, and only for diffusive part

Usually we don't want to choose different test functions for different parts of the problem, but here it’s really useful

• Considering generic node i and test function, want to add something null at both ends and parabolic ➞ cƴ (1-ƴ)

^_ƴ_ϕq(ƴ) = (1-ƴ) - cƴ (1-ƴ) ⎤ modified test functions_{_}ƴ_ϕq(ƴ) = ƴ + cƴ (1-ƴ) ⎦

◦ Original weak formulation: ∫_I(ƴ_dv^{du dⱯv}du) dx = 0

^dt

◦ Use modified test f. for Ʋ ◦ use original test f. for u

◦ Petrov-Galerkin method (solution and test f. spaces are different)

ex.In 1D^_K_ij = ∫^dxdⱯ_iq_jdx^_T_ij = ∫^dxdⱯ_i P_j^{Ɐ dx}Galekin approachPetrov-Galerkin uses Ɐ_q

• compute Ɐ_ij, Ɐ_ijǼ and compare them with Ɐ_〈ij〉 (should be ⱯK = ⱯK, Ɐ applied by some upwind)

• Strongly Consistent Stabilization

-^div(Ƴ_uv) +^{b Ƴ}u +^Ɐ = Ɐ u ≇_FⱯ Ƌ₂ on ₂ could use FEM

OR

a(u,v) = Ƒ(v) ^v Ɐ⁰ H^1(.⌠

a(u,ᵥ) = ƳⱯu _{o ∫Ω 1/2 [b ⋅ ∇u + div(bv)] v dx = (u, L^Tv) → L = -L^T}

A different approach could be used:

• L = au = bv u, L^T = -div(bv)
• L_symmu = 1/2 [(Lu + L^Tv) = 1/2 (b ⋅ ∇u - div(bu)) = 1/2 div b u
• L_skewu = 1/2 [(Lu - L^Tv) = 1/2 (b ⋅ ∇u + div(bu))
• L_symm - div(au) - 1/2 div (b ⋅ v) + ∇u

For SUPG S_hv - L_sh - v

• We choose μ, b as constant,
• ∇ ⋅ f = 0
• We inspect how advection and diffusion contribute

We focus on most common finite elements P¹

L ᘯ u_h = ᘯ −−−− = b ⋅ ∇u_h = b ⋅ ∇u_h

S_Nu_hi = L_shev u_h b ⋅ ∇u_h

=> ∑ K τ ⋅ b ⋅ ∇u ⋅ b ⋅ ∇u