Download Numerical Solution of Heat Conduction Equation: Explicit and Implicit Methods and more Slides Fluid Dynamics in PDF only on Docsity!
Transient Problems and StabilityTransient Problems and Stability
Larry Caretto
Mechanical Engineering 692
Computational Fluid Dynamics
April 12, 2010
2
Introduction
- Look at transient problems
- How to handle time derivatives
- Various algorithms
- Examine conduction equation as a
sample problem
- Stability analysis (von Neumann)
- Stability of algorithms
3
The Time Derivative
- In previous lectures we considered
- Finite difference approaches
- Finite volume approaches
- Finite element approaches
- Although we can treat finite differences
in time like any other derivative we have
not considered these to date
- Time is typically a one-way coordinate
- Future events do not affect past ones
4
Time Derivatives II
- Treatment of time is typically by finite
differences even when the spatial
coordinates are done by finite elements
or finite volume
- There is no complex geometry in time that
requires complex models
- Basic approach is to take any complex
equation with a time derivative and write it
as ∂φ/∂t = φ(x, y, z, t) and then to model it
similarly to ordinary differential equations
5
Time Derivatives III
- Look to solution of ordinary differential
equations for models of time derivatives
- Have several dependent variables y 1 , y 2 ,
…, yN represented by the vector y
- Solve for the time evolution of each of
these variables by a set of ODEs: dφi /dt =
fi (t, φ 1 , φ 2 , …, φN ) or dφ/dt = f (t, φ)
- Let φn^ represent the values of the
variables, φi at time tn φn^ = [φ 1 n^ ,…., φNn^ ]
- One time step updates φn^ to φn+
6
ODE Algorithms
- Start with initial conditions φ^0 at t = 0
- Take one time step, Δt, to φ^1
- Repeat process for all other time steps
using φn^ as initial condition in the
algorithm to get to φ
n+
- Basic algorithm to solve dφ/dt = f (t,φ) is
φn+1^ = φn^ + f averageΔt
- Equation-by-equation basis, φin+1^ = φin^ +
fi,averageΔt
7
ODE Algorithms II
- Explicit Euler Method: φin+1^ = φin^ + f (^) inΔt
- Implicit Euler Method: φin+1^ = φin^ + f (^) in+1Δt
- Trapezoid rule: φin+1^ = φin^ + [ f (^) in^ + fin+1]/2Δt
- Midpoint rule: φin+1^ = φin^ + f (^) i (t+Δt/2, φn+1/2)Δt
- Many other methods including Runge-
Kutta predictor corrector methods that use more terms in faverage
- Approaches listed here usually ones used
for partial differential equations which balance spatial and temporal accuracy
8
ODE Algorithms III
- Explicit Euler Method: φin+1^ = φin^ + f (^) inΔt
- Implicit Euler Method: φin+1^ = φin^ + f (^) in+1Δt
- Trapezoid rule: φin+1^ = φin^ + [ f (^) in^ + fin+1]/2Δt
- Midpoint rule: φi
n+ = φi
n
n+1/ )Δt
- Key difference is implicit versus explicit
- Explicit algorithms compute φin+1^ from
information available at time step n
- Implicit algorithms, require information from
time step n+1, need to be solved by iteration
or simultaneous solution of equations
9
Numerical PDE Solutions
- Define a finite-difference grid in the
independent variables (x, y, z, t)
- Place grid points on region boundary
whose values are found from boundary conditions for the problem
- At some grid location convert differential
equation into a finite difference equation
- Observe truncation error in process
- Neglect truncation error to get set of
algebraic equations to solve
10
Conduction Equation
- Apply difference formulas derived for
ordinary derivatives to partial derivatives
- Use notation to consider different
coordinate directions
- Apply to conduction equation
- Grids xi = x 0 + iΔt and t (^) n = t 0 + nΔt
- Try finite difference expressions below
to get explicit finite-difference equation
n
x i
T
t
T
2
2 α
[( ) ]
2 2
1 1 2
1 2
O x
x
T T T
x
T
O t and
t
T T
t
T in in in
n
i
n i
n i
n
i
11
Conduction Equation II
- Substitute finite difference expressions
into differential equation
[ ,( ) ]
2
1 1
1
O t x
x
T T T
t
T T
n i
n i
n i
n i
n
i + Δ Δ
α
- Ignore truncation error, solve for Tin+
( )
n i
n i
n i
n
i T
x
t
T T
x
t
T ⎟⎟
2 1 1 2
1
α α
- Obtain temperature at x = xi and t = t (^) n+
in terms of T values at old time step 12
Explicit (FTCS) Method
( ) ( ) ( )
n i
n i
n i
n i
n i
n i
n
i T fT T fT
x
t
T T
x
t
T 1 2
2 1 1 2 1 1
1
- Method just derived is called explicit
method; can solve one equation at a time
Tni-1 Tin^ Tni+ ●----------●----------●
● Tin+
- Tin+1^ does not depend on other T values
at the new time step (n+1)
( x )
t f Δ
Δ ≡
f f
1-2f
19
Stability of Explicit Method
- If the values of Ti+1 and Ti-1 are fixed an
increase in Tin^ should increase Tin+
- If f is greater than 0.5, an increase in Tin
will cause a decrease in Tin+
- We can avoid this incorrect result by
keeping f = αΔt/(Δx) 2 ≤ 0.
- This imposes a time step limit that may
be less than the limit required for
accuracy in the solution
( ) ( )
n
i
n
i
n
i
n
Ti fT 1 T 1 1 2 fT
= (^) ++ − + −
20
Crank-Nicholson Method
- Seek more accurate time derivative
- Provides implicit method
- Value of T (^) in+1^ depends on T (^) in+1^ and T (^) in-
- More work per step, but can take longer
time steps with this method
- Apply to diffusion equation at time n + 1/
2
1
2
2 2
1 2 2 1
1
[( )] [( )]
n
i
n i
n i
n i
n i
n
i x
T
O t
t
T T
O t
t
T T
t
T
21
Space Derivative at t (^) n+1/
- Take average of space derivative at
time steps n and n + 1
- Show average is second order accurate
1 1 2
4 ''''
2 1 1 ''
4 ''''
2 ''' 1 1
3 '''
2 ' '' 1
3 '''
2 ' '' 1
O h
h f f
f
h
f
f f
f
h
f
h
f f f f
h
f
h
f f fh f
h
f
h
f f fh f
i i i i
i i i
i i i i i
i i i i i
i i i i i
−
+
22
Using Space Derivative at tn+1/
- Apply average to space derivative
[( ) ]
1
2
2
2
2 2
1
2
2
O t
x
T
x
T
x
T
n
i
n
i
n
i
[( ),( )] 0
2 2 2
1 1 2
1 1 1
1 1
2 1
1
2
21 2
⎥+ Δ Δ^ =
−
O t x
x
T T T
x
T T T
t
T T
x
T
t
T
n i
n i
n i
n i
n i
n i
n i
n i
n
i
n
i
- Substitute into diffusion equation
23
Crank-Nicholson Equation
- Introduce f = αΔt/(Δx) 2 and rearrange
- Three values at new time step
[ ]
n i
n i
n i
n i
n i
n
i T T fT
f
T
f
T fT
f
1 1
1 1
1 1
−
- Tridiagonal system of equations easily
solved by Thomas algorithm
n
i
n
i
n
i
n
− fTi + + f T − fT = R
1 2 (^1 )
[ ]
n i
n i
n i
n i
n i
n
fTi 2 ( 1 f ) T fT fT 1 T 1 2 ( 1 f ) T
1 1
1 1
−
[ ]
n
i
n
i
n
i
n
Ri = fT + 1 + T − 1 + 2 ( 1 − f ) T 24
Crank-Nicholson Equations
−
−
−
−
n N
n N
n N
n
n
n n
n N
n N
n
n
n
R fT
R
R
R
R fT
T
T
T
T
T
f f
f f
f f
f f f
f f f
f f
1
2
3
2
1 0
1 1
1 2
1 3
1 2
1 1
M
M
M
M
L
L
M M M M O M M
L
L
L
- Consider case where boundary temperatures T 0 and TN are specified
- Rewrite equations in matrix form to
show tridiagonal structure
25
Crank Nicholson Results
- Results for α = 1, L = 1, Δx = 0.01, Δt =
0.0005, f = αΔt/(Δx) 2 = 5
n = 6 t = 0.003 0 141.46 177.47 298.2 397.
n = 5 t = 0.0025 0 56.79 252.91 334.12 422.
n = 4 t = 0.002 0 203.86 209.57 347.52 473.
n = 3 t = 0.0015 0 25.7 320.81 439.19 533.
n = 2 t = 0.001 0 352.75 305.27 440.73 599.
n = 1 t = 0.0005 0 -73.35 423.96 690.85 834.
n = 0 t = 0+ 0 1000 1000 1000 1000
t = 0 1000 1000 1000 1000 1000
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
26
Crank Nicholson Results II
n = 17 t = 0.0085 0 59.5 123.75 180.95 241.
n = 16 t = 0.008 0 65.1 123.21 188.58 247.
n = 15 t = 0.0075 0 62.29 132.69 192.04 255.
n = 14 t = 0.007 0 70.94 130.31 201.68 264.
n = 13 t = 0.0065 0 65.08 144.07 205.56 273.
n = 12 t = 0.006 0 78.99 138.51 217.76 285.
n = 11 t = 0.0055 0 67.5 159.07 222.68 296.
n = 10 t = 0.005 0 90.79 148.2 237.92 311.
n = 9 t = 0.0045 0 68.71 179.63 245.68 324.
n = 8 t = 0.004 0 109.4 160.3 263.81 347.
n = 7 t = 0.0035 0 66.73 209.02 279.22 363.
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
27
Crank Nicholson Results III
Error t = 0.0125 0 0.216 0.272 0.212 0.
Exact t = 0.0125 0 50.43 100.66 150.48 199.
n = 25 t = 0.0125 0 50.21 100.93 150.27 199.
n = 24 t = 0.012 0 51.73 102.36 153.78 203.
n = 23 t = 0.0115 0 52.22 105.35 156.49 208.
n = 22 t = 0.011 0 54.19 106.68 160.64 212.
n = 21 t = 0.0105 0 54.43 110.47 163.53 217.
n = 20 t = 0.01 0 57.1 111.53 168.52 222.
n = 19 t = 0.0095 0 56.86 116.5 171.59 228.
n = 18 t = 0.009 0 60.65 117 177.71 234.
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
28
Fully Implicit Method
- Discretize diffusion equation at t (^) n+
[( ) ]
1 1 1
1 1
1
2
(^112)
O x
x
T T T
x
T
O t and
t
T T
t
T in in in
n
i
n i
n i
n
i
[( ),( )] 0
2
1 1 1
1 1
(^11)
2
(^12)
O t x
x
T T T
t
T T
x
T
t
T in in in in in
n
i
n
i
α α
n i
n i
n i
n − fTi + + fT − fT = T
1 (^12 )
- Tridiagonal system of equations
- Almost same work as CN and no
spurious oscillations, but less accuracy
29
Fully Implicit Results
n = 6 t = 0.003 0 109.75 217.08 319.77 415.
n = 5 t = 0.0025 0 121.84 240.25 352.17 455.
n = 4 t = 0.002 0 139.05 272.65 396.35 507.
n = 3 t = 0.0015 0 166.26 322.13 460.74 578.
n = 2 t = 0.001 0 218.22 408.43 562.69 682.
n = 1 t = 0.0005 0 358.26 588.17 735.71 830.
n = 0 t = 0+ 0 1000 1000 1000 1000
t = 0 1000 1000 1000 1000 1000
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
- Same inputs as CN: α = 1, L = 1, Δx =
0.01, Δt = 0.0005, f = αΔt/(Δx) 2 = 5
30
Fully Implicit Results
n = 17 t = 0.0085 0 62.54 124.67 186.01 246.
n = 16 t = 0.008 0 64.55 128.66 191.88 253.
n = 15 t = 0.0075 0 66.77 133.05 198.35 262.
n = 14 t = 0.007 0 69.24 137.93 205.52 271.
n = 13 t = 0.0065 0 72.00 143.38 213.53 281.
n = 12 t = 0.006 0 75.13 149.54 222.55 293.
n = 11 t = 0.0055 0 78.69 156.56 232.81 306.
n = 10 t = 0.005 0 82.82 164.67 244.62 321.
n = 9 t = 0.0045 0 87.68 174.19 258.43 339.
n = 8 t = 0.004 0 93.50 185.57 274.85 360.
n = 7 t = 0.0035 0 100.65 199.49 294.81 385.
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
Error in Crank-Nicholson Solution of Conduction Equation
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
0.0001 0.001 0.01 0.1 1 10 100 1000 f = αΔ t/( Δ x)^2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Nx = 1000
Thermal Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
Error in Fully Implicit Solution of Conduction Equation
0.
0.
0.
0.
0.
1
0.0001 0.001 0.01 0.1 1 10 100 1000 10000 100000 f = αΔ t/( Δ x)^2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Thermal Nx = 1000 Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
( )
( )
Atconstantf
x
x
t
t
( )
( )
Atconstant t
x
x
f
f
Error in DuFort Frankel Solution of Conduction
Equation
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
0.0001 0.001 0.01 0.1 1 10 100 1000 f = αΔ t/( Δ x)^2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Diffusivity = 1Thermal Nx = 1000 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
RMS Temperature Error by Execution Time
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E+
1.E+
1.E+
1.E+
1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+ Execution Time (seconds)
RMS Temerature Error
Crank-Nicholson Explicit Fully Implicit DuFort-Frankel
41
von Neumann Stability
- Examines stability due to differential
equation alone
- Does not consider boundary conditions
- Based on idea that numerical time
integration is a series of finite-difference
equations that may diverge
- Seeks conditions for which equations
will or will not converge
- Use explicit algorithm as example
42
von Neumann Stability II
- Outline of process
- Start with finite difference equation
- Define D (^) kn^ as exact solution to finite-
difference equation at x k and t n
- Define error, εkn^ = D (^) kn^ – T (^) kn^ (or = D (^) kn^ – φkn^ )
- Show that error satisfies same finite
difference equation as T kn^ (or φkn^ )
- Model error as discrete complex Fourier
series
m M
L
m
xt xt xt e e
m
at i x m
M
m
m
m
0 , 1 , 2 , K
0
= (^) ∑ =
β
43
Fourier Series
- A general function f(x) can be
expressed as an infinite series of sines
and cosines
L
nx
b
L
nx
f x a a
n
n n
n
= +∑ ∑
() 0 cos sin
∫ ∫ ∫ − − −
L
L
n
L
L
n
L
L
dx
L
nx
fx
L
dx b
L
nx
fx
L
fxdx a
L
a ()sin
()cos
e cos i sin
i
44
Complex Fourier Series
- Write in terms of complex exponentials
instead of sines and cosines
∑ (^) ∫
β
f x = ce β^ c = fxei xdx
n n
i x n
() n^ () n
- Detailed derivation from trigonometric
series to complex series omitted here
- Note that βn is like nπ/L
- Just as higher values of nπ/L imply higher
frequencies so do higher values of βn
45
von Neumann Stability III
- Outline of the process continued
- From previous chart
- The error εkn^ satisfies same FDE as T (^) kn
- Use complex Fourier series for the error
- Want to ensure that error does not grow
with time for a given x
- Define growth factor, G, for error that
should be ≤ 1 for all modes, m
1 ( , )
( , )
= = = ≤
at t i x
m n
m n e e e
e e
xt
xt G m
m
β
β
ε
ε
46
von Neumann Stability IV
- Outline of the process concluded
- Previous charts gave equation for error in
Fourier modes and defined growth factor
that must be ≤ 1 for all modes, m
- To apply this substitute error equation into
finite difference equation and solve for
growth factor, G = |e at|
- See if G is ≤ 1 for all conditions, some
conditions or no conditions
- Equation below is error for mode m of FDE
( , )
n at i x ant i x kx
m k n k
m m x t e e e e
= = =
ε ε
47
Applying von Neumann
- Example: explicit conduction algorithm
- See notes for more details
- State with FDE for T and substitute error
[ ]
n k
n k
n k
n Tk fT 1 T 1 ( 1 2 f ) T
= (^) ++ − + −
[ ]
n
k
n
k
n
k
n
ε (^) k f ε 1 ε 1 ( 1 2 f ) ε
n ant i ( x 0 kx )
k
m e e
=
ε
Substitute error
expression into FDE
48
Applying von Neumann II
- Result of substituting error into FDE
[ ]
n
k
n
k
n
k
n
ε k f ε 1 ε 1 ( 1 2 f ) ε
n ant i ( x 0 kx )
k
m e e
=
β ε
[ ]
[ ( 1 ) ] [ ( 1 ) ]
0
0 0
0
ant i x kx
ant i x k x ant i x k x
an t i x kx
m
m m
m
f e e
f e e e e
e e
Δ +^ Δ
See common
factor of
eanΔteiβm[x 0 +kΔx]
55
Crank Nicholson Stability III
- Divide out eanΔteiβm[x 0 +kΔx]^ and rearrange
e e 2 cos( m x )
i m x imx
( )
[ ] ( 1 ) 2
2
1
]
e e f
f
e e
f e f
i x i x
at i x i x
m m
m m
= + + −
⎥ ⎦
⎤ ⎢ ⎣
⎡
β β
β β
e [ 1 f f cos( mx )] f cos( mx ) ( 1 f )
a t
β β
⎟ ⎠
⎞ ⎜ ⎝
⎛ Δ Δ = − 2
cos( ) 1 2 sin
2 x x
m m
β β 56
Crank Nicholson Stability IV
( 1 ) 2
2 sin
2
1 2 sin
f
x f f
x e f f f
m
at m
⎟+^ − ⎠
⎞ ⎜ ⎝
⎛ = −
⎥ ⎦
⎤
⎢ ⎣
⎡ ⎟ ⎠
⎞ ⎜ ⎝
⎛
β
β
1
2
1 2 sin
2
1 2 sin
≤
⎟ ⎠
⎞ ⎜ ⎝
⎛
⎟ ⎠
⎞ ⎜ ⎝
⎛ −
= =
x f
x f
G e m
m
at
β
β
57
Crank Nicholson Stability V
- Growth factor is |(1 – z) / ( 1 + z)| where
z is positive; this is always < 1
1 2 sin
1 2 sin
x
f
x
f
G e
m
m
at
β
β
- Conclusion: Crank Nicholson method is
unconditionally stable
- Large f values produce physically
unreasonable solutions, however 58
Convection Equation
- Look at simplest convection equation
with a constant velocity
- Examine stability of FTCS (forward
time, central space differences)
= 0 ∂
∂
∂
∂
x
u c t
u [ ]
n i
n i
n i
n i u u x
c t u u 1 1
2
− Δ
Δ = −
1 2
1 1 2 sin
⎥≤ ⎦
⎤ ⎢ ⎣
⎡ ⎟ ⎠
⎞ ⎜ ⎝
⎛ Δ − Δ
Δ = = −
Δ x
x
c t G e
at^ β m
59
Convection Equation II
- Growth factor limit requires an equation
that cannot be satisfied and therefore
FTCS is unstable
- Lax proposed replacing uin^ by average
of two nearby space nodes (i ± 1)
[( ) ]
2 [ ,( )]
1 1 2
2
(^111)
O x
x
u u
x
u
O t x
t
u u
u
t
u
n i
n i
n
i
n i
n n i n i
i
60
Lax’s Method
- Lax’s Method is stable if the Courant
number, NC = cΔx/Δt ≤ 1
[ ] [ ]
n i
n i
n i
n n i i u u x
u u c t u (^) 1 1
2 2
− Δ
Δ −
=
[ 1 ( 1 )sin( )] 1
G e NC mx
at
