Advection Diffusion Problems

Numerical methods for differential equations

Advection-diffusion

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

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?

Peclet number 0.

Numerical methods for differential eqns

Advection-diffusion

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.

Finite element discretization

Discretization with finite elements.

Dirichlet b.c.

Could also have constant term.

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

E.g. discretization of transport term: (3^rd 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

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.).

Oscillations and mesh size

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

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/².

Generalization of the hypothesis

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 + (ϕ(ρe)) ϕ (t) such that lim ϕ (t) = 0 t→+∞

If ϕ(t) = 0 go back to the original problem (centered approx 2^nd order).

ϕ(t) = upwinding, 1^st order.

ϕ(t) = t - t + B(t), B(t) = _tt, t ₁t0 Bernoulli function.

Schapfelter - Gummel (SG) 2^nd order.

For constant facing term solution in the nodes analytical (TO).

Multi-dimensional case

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, but vector 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 = 0

Add μ, 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.

Modified test functions

• 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).

Petrov-Galerkin method

Ex. In 1D:

^_K_ij = ∫^dxdⱯ_iq_jdx

^_T_ij = ∫^dxdⱯ_i P_j^{Ɐ dx}

Galekin approach.

Petrov-Galerkin uses Ɐ_q.

Compute Ɐ_ij, Ɐ_ijǼ and compare them with Ɐ_⟨ij⟩ (should be ⱯK = ⱯK, Ɐ applied by some upwind).

Strongly consistent stabilization

• Strongly Consistent Stabilization-^div(Ƴ_uv) +^{b Ƴ}u +^Ɐ = Ɐ u ≇_FⱯ Ƌ₂ on ₂ could use FEMOR.

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}

Alternative approach

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