Download Stable Numerical Solutions for Parabolic PDEs with Variable Coefficients - Prof. P. Blomgr and more Study notes Mathematics in PDF only on Docsity!
Numerical Solutions to PDEs Lecture Notes #11 — Parabolic PDEsStability, Boundary Conditions;Convection-Diffusion; Variable Coefficients
Peter Blomgren, 〈[email protected]
Department of Mathematics and Statistics
Dynamical Systems GroupComputational Sciences Research Center San Diego State UniversitySan Diego, CA 92182-7720 http://terminus.sdsu.edu/^ Spring 2010
Peter Blomgren,
〈[email protected]
〉^ Stability
— (1/34)
Outline^1 Recap
Parabolic PDEs Schemes: Forward/Backward-Time Central Space Schemes: Crank-Nicolson, Du-Fort Frankel 2 Stability: Lower Order Terms One-step Schemes 3 Dissipation and Smoothness Crank-Nicolson 4 Boundary Conditions 2nd Order One-Sided; Ghost Points 5 Convection-Diffusion Numerics Upwind Differences 6 Variable Coefficients Peter Blomgren,
〈[email protected]
〉^ Stability
— (2/34)
Recap
Parabolic PDEsSchemes: Forward/Backward-Time Central SpaceSchemes: Crank-Nicolson, Du-Fort Frankel
Last Time
1 of 3
A quick introduction to parabolic PDEs: Our model equation isthe one-dimensional heat equation.Exact solutions to the 1D heat equation in infinite space, using theFourier transform.The solution corresponds to a damping of the high-frequencycontent of the initial condition.
⇒^ the parabolic solution operator
is^ dissipative
For^ t^ >
0, the solution of the heat equation is infinitely differentiable.Since parabolic PDEs do not have any characteristics, we needboundary conditions at
every^
boundary. Typically we specify
u
(fixed temperature, “Dirichlet”), the [normal] derivative
ux
(temperature flux, “Neumann”), or a mixture thereof.^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (3/34)
Recap
Parabolic PDEsSchemes: Forward/Backward-Time Central SpaceSchemes: Crank-Nicolson, Du-Fort Frankel
Last Time
2 of 3
Numerical Schemes
for^ ut^
=^ buxx
+^ f^ :
Forward-Time Central-Space
n+1 v − m^
n v m=^ k
nv m+1b
n− 2 v m n+ v m− 1 (^2) h
n+ f m
Explicit; stable when
bμ^ ≤^
1 , where 2
k μ = h ; order-(1,2); dissipa- 2
tive of order 2. Backward-Time Central-Space
n+1 v − m^
n v m=^ k
n+1v m+1b
n+1− 2 v m
n+ v (^) +1 m−^1 (^2) h
n+1+ f m
Implicit; unconditionally stable; order-(1,2); dissipative of order2.^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (4/34)
Recap
Parabolic PDEsSchemes: Forward/Backward-Time Central SpaceSchemes: Crank-Nicolson, Du-Fort Frankel
Last Time
3 of 3
Crank-Nicolsonn^ v^ +1n−^ v^ m m^ k^
"^ nv^ b= 2 +1n−^2 v^ m+^
n+1+1+ v (^) m m− 1 (^2) h
nv −m+1^ +
n 2 v +^ v^ m^
#n^1 m− 1 + (^2) h^2 »n+1^ f^ +f^ m^
Implicit; unconditionally stable; order-(2,2); dissipative of order2, when
μ^ is constant. Du-Fort Frankel
(“fixed leapfrog”) n+1 v − m^
n−^1 v m= 2 k
nv m+1 b
n+1− (v m^
n−^1 + v m^
n) + v m− 1 (^2) h
n+ f m
Explicit; unconditionally stable; order-(2,2); dissipative of order 2,when^ μ
is constant.
It is only consistent if
k/h^ tends to zero
as^ h^ and
k^ go to zero
Peter Blomgren,
〈[email protected]
〉^ Stability
— (5/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
One-step Schemes
Lower Order Terms and Stability
1 of 2
For^ hyperbolic
equations we have the following result: Theorem (Stability of One-Step Schemes) A consistent one-step scheme for the equation
u+^ aut^
+^ bu^ x
is stable
if and only if it is stable for this equation when b^ =^ 0.^
Moreover, when k
=^ λh, and
λ^ is a constant, the stability
condition on g
(hξ,^ k, h)^ is^ |g (θ,^0 ,^ 0)
| ≤^1.
Similar results
do not always apply directly
to^ parabolic
equations.^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (6/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
One-step Schemes
Lower Order Terms and Stability
2 of 2
The problem is that the contribution to the amplification factor from the firstderivative is sometimes (often?)
“^ ”√ Ok which is greater than
O^ (k) as
k^ ց^ 0.
For instance
, the forward-time central-space scheme applied to u=^ but^ xx
−^ au+ x^ cu^ gives the amplification factor^ g^ = 1
−^4 bμ^ sin
„^ «θ 2 −^2 iaλ^ sin(θ
) +^ ck
The term
ck^ does not affect stability, but the term containing
√ λ = k μ^ cannot be
dropped when
μ^ is fixed. In this particular case, we get^ |g
„ (^2) |=^1 −^4 bμ^ sin
„^ ««θ 2 2 22 +^ akμ
(^2) sin(θ)
and since the first derivative term gives an
O^ (k) contribution to
(^2) |g |, it does not
affect stability. (Strikwerda, p.149) This is also true for the backward-timecentral-space, and Crank-Nicolson schemes.^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (7/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Dissipation and Smoothness^ The fact that a dissipative one-step scheme for a parabolic equationgenerates a numerical solution with increased smoothness as
t^ ր
(provided that
μ^ is constant) is a key result, so lets show that it is indeed true...We start with the following theorem Theorem A one-step scheme, consistent with
u=^ but^
,xx
that is dissipative of order 2 with
μ^ constant satisfies n+1^2 ‖v ‖^2
n∑+ ck ν=
ν^ ‖δv ‖+ 20 ≤ ‖v^2
(^2) ‖ 2
for all initial data v
0 and n ≥^0.
Peter Blomgren,
〈[email protected]
〉^ Stability
— (8/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/20,^ dt
μ^ = 20
−1^ −0. 0
0.5^1 T = 0.0000, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1500, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/20, dt = 1/20 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ The Crank-Nicolson scheme applied to the initial condition in panel #1, with zero-flux boundary conditions.
We know that Crank-Nicolson is non-dissipative if
λ^ remains
constant (see next slide).^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (13/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/40,^ dt
μ^ = 40
−1^ −0.^
0 0.^
1 T = 0.0000, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0.^
0 0.^
1 T = 0.1500, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/40, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ The Crank-Nicolson scheme: here we have cut both
h^ and^ k^
in half compared with
the previous slide. On the next slide we show the result of keeping
μ^ =^ k/h
2 constant, in
which case the scheme is dissipative.^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (14/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/40,^ dt
μ^ = 20
−1^ −0. 0
0.5^1 T = 0.0000, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1500, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/40, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ The Crank-Nicolson scheme: here, we finally get some damping in the oscillationsof the solution. — Dissipation is a convergence result!^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (15/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/80,^ dt
μ^ = 80
−1^ −0.^
0 0.^
1 T = 0.0000, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0.^
0 0.^
1 T = 0.1500, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/80, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ Surprisingly(?), refinining in
x^ brings back the over-shoot artefacts.
Peter Blomgren,
〈[email protected]
〉^ Stability
— (16/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/20,^ dt
μ^ = 10
−1^ −0. 0
0.5^1 T = 0.0000, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1500, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/20, dt = 1/40 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ Coarsening in
x^ (dx^ = 1
/20, instead of
dx^ = 1/
40 lessens the “carrying capacity”
of high-frequency content of the grid...^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (17/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
Crank-Nicolson
Example: Crank-Nicolson
dx^ = 1
/20,^ dt
μ^ = 5
−1^ −0.^
0 0.^
1 T = 0.0000, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.0500, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.1000, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0.^
0 0.^
1 T = 0.1500, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2000, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
−1^ −0. 0
0.5^1 T = 0.2500, dx = 1/20, dt = 1/80 (^1) 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
Figure:^ Peter Blomgren,
〈[email protected]
〉^ Stability
— (18/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
2nd Order One-Sided; Ghost Points
Boundary Conditions
(Again)
Since parabolic problems require boundary conditions at everyboundary, there is
less need for numerical boundary conditions
compared with hyperbolic problems.We briefly discuss implementation of the
physical boundary
conditions
: — Implementing the Dirichlet (specified values at the boundary points) boundary conditions is straight-forward.The Neumann (specified flux/derivative) is more of a problem; forinstance,
one-sided differences ∂u(tn,^ x
)v 0 ≈^ ∂x nn −^ v^1 0 h^
∂u(, tn,^ x)M^ ∂x^
nv −M^ ≈ n v M−^1 h
can be used, but these
are however only first-order accurate
and will degrade the accuracy of higher-order schemes.^ Peter Blomgren,
〈[email protected]
〉^ Boundary Conditions
— (19/34)
Stability: Lower Order TermsDissipation and SmoothnessBoundary Conditions
2nd Order One-Sided; Ghost Points
More Accurate Boundary Conditions
1 of 2
Second order one-sided accurate boundary conditions are given by^ ∂u(t,^ n
x)^0 ≈ ∂x
n−v + 4 2 nv −^31 nv 0 2 h^
∂u(, t,^ x)nM^ ∂x^
nv M−≈ −^4 v^2 n+ 3M−^1
nv M 2 h
It is sometimes useful to use second-order central differences andintroduce
“ghost-points”
for the boundary conditions,
e.g.
∂u(t,^ n x)^0 ≈^ ∂x
nv −^ v 1 n −^12 h
How is this useful? — Consider a given flux condition u(t,^ xx^ n
) =^ ϕ 0 (t), thenn nn v −^ v 1 −
1 =^ ϕn 2 h
⇔^
nv=^ −^1
nv −^2 h 1 ϕn
Peter Blomgren,
〈[email protected]
〉^ Boundary Conditions
— (20/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Numerics, 3 of 3
The condition
α^ ≤^ 1, can be re-written
2 b h ≤ a
which is a restriction on the spatial grid-spacing.The quantity
a corresponds to the b^
Reynolds number
in fluid flow,
or the^ Peclet number
in heat flow.
The quantity
ha α = 2 (sometimes 2b
α) is often called the
cell
Reynolds number
or the cell Peclet number
If the grid-spacing
h^ is too large, then the numerical solution cannot resolve the physics and oscillations occur. Theseoscillations are
not^ due to instability (as long as the stability criterion is satisfied, of course) and do not grow; they are only aresult of inadequate resolution.^ Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (25/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Example #
0 2
4 6
8 10 T = 0.0000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.2000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.4000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.6000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.8000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 1.0000, dx = 0.10, dt = 0.0025 (^1) 0.5 (^0) −0.
Figure:^ (Forward-Time Central-Space) Convection-diffusion with
a^ = 10,
b^ = 0.1,
h^ =
0.^1 >^0 .02,
k^ = 0.0025,
μ^ = 1/^4
<^1 /2. We are stable, but have not resolved the physics.
Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (26/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Example #
0 2
4 6
8 10 T = 0.0000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.2000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.4000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.6000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.8000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 1.0000, dx = 0.02, dt = 0.0001 (^1) 0.5 (^0) −0.
Figure:^ (Forward-Time Central-Space) Convection-diffusion with
a^ = 10,
b^ = 0.1,
h^ =
0.^02 ≤^0 .02,^ k^ = 0
.0001,^ μ^
= 1/^4 <
1 /2. We are stable, and have resolved the physics.
Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (27/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Upwind Differences, 1 of 3
In example #2 we had to push the resolution to
h^ = 0.02 (601 points in
[−^1 ,^ 11]) and
k^ = 0. 0001 (10001 time-steps in [
,^ 1]), for a grand total of
6,010,601 space-time grid points. That is a ridiculously high price to payfor such a simple 1D problem!!!One way to avoid the resolution restriction is to use
upwind differencing
of the convection term. This corresponds to a switching betweenbackward differencing when
a^ >^ 0, and forward differencing when
a^ <^ 0,
e.g.^ only differencing in the direction where the (hyperbolic)characteristics come from:n+1^ v^ −^ m^
n v m++^ a k »^ nv^ −^ v^ m^
- nm− 1 +^ a h »^ nv^ −m+
- n− v m = h nv −m+1^ b
n 2 v +^ vm^ n m−^12 h
or
n+1 v = [1 m^ −^2 bμ(1 +
n α)] v +m^ n bμv m+^ +^ bμ(1 + 2
nα)v m−^1
Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (28/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Upwind Differences, 2 of 3
The restriction
2 b h ≤ |a is replaced by| 2 b μ^ +^ |a|
λ^ ≤^1
which is much less restrictive when
b^ is small and
a^ large. If we
want^ μ
i.e. k^
(^2) = h/4, then we must have
(^4) h ≤ a (b^1 −^2
which with
a^ = 10 and
b^ = 0.
1 as in the previous examples is 0.
— 19 times that of the previous restriction.We have, however, also sacrificed the spatial second orderaccuracy, since the first-order upwind difference is first order.^ Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (29/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Example #
0 2
4 6
8 10 T = 0.000000, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.214375, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.428750, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.612500, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.826875, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.980000, dx = 0.35, dt = 0.030625 (^1) 0.5 (^0) −0.
Figure:^ (Upwinding) Convection-diffusion with
a^ = 10,
b^ = 0.1,
h^ = 0.^35
≤^0 .38, k^ =
0 .030625,
μ^ = 1/^4
<^1 /2. We are stable, and have resolved the physics. Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (30/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Example #
0 2
4 6
8 10 T = 0.000000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.200000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.400000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.600000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 0.800000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
0 2
4 6
8 10 T = 1.000000, dx = 0.40, dt = 0.040000 (^1) 0.5 (^0) −0.
Figure:^ (Upwinding) Convection-diffusion with
a^ = 10,^ b^ = 0.1,^
h^ = 0.^40 ≥^0 .38,^
k^ = 0.04,
μ^ = 1/^4
<^1 /2. We are stable, but have not resolved the physics. Peter Blomgren,
〈[email protected]
〉^ Convection-Diffusion
— (31/34)
Convection-DiffusionVariable Coefficients
NumericsUpwind Differences
The Convection-Diffusion Equation
Upwind Differences, 3 of 3
The upwind scheme
n+1 v −^ m^
nv m+^ a k
nnv −^ v^ m^ m −^1 =^ b h
nv −m+^
n 2 v +^ m^ nv m−^12 h
can be rewritten in the formn^ v^ +1n−^ v^ m^ m^ k^
nv m+1+ a
n− v m−^12 h ( =^ b^ +
)^ ahv 2 n −^2 m+^
nnv +^ v^ m^ m −^1 (^2) h
We see that upwinding corresponds to changing the diffusion coefficient,or^ adding artificial viscosity
to suppress oscillations.
There has been much debate regarding the value of theseartificial-viscosity solutions; clearly they may only give qualitativeinformation about the true solution.More details on solving the convection-diffusion equation numerically canbe found in
K.W. Morton
,^ Numerical Solution of Convection-Diffusion
Problems
, Chapman & Hall, London, 1996. Peter Blomgren,^ 〈[email protected]
〉^ Convection-Diffusion
— (32/34)
