Direct Numerical Simulations of Multiphase Flows

1
Elementary Numerical Analysis Finite Difference
Approximations-II
Instructor Hong G. Im University of Michigan
Fall 2005
Lecture notes were prepared jointly with Prof.
Gretar Tryggvason, Department of Mechanical
Engineering, Worcester Polytechnic Institute.
2
Outline

• The linear advection-diffusion equation
• Finite difference approximations
• Finite difference approximation of the linear
  advection-diffusion equation
• Accuracy
• Stability
• Consistency
• Higher order finite difference approximations in
  space
• Time integration implicit versus explicit and
  higher accuracy

3
Solving Partial Differential Equations
4
Model Equations
We will use the model equation
to demonstrate how to solve a partial
differential equation numerically.
Although this equation is much simpler than the
full Navier Stokes equations, it has both an
advection term and a diffusion term. Before
attempting to solve the equation, it is useful to
understand how the analytical solution behaves.
5
Model Equations
For initial conditions of the form
It can be verified by direct substitution that
the solution is given by
k 1.0 U 1.0 D 0.05 A 1.0
which is a decaying traveling wave
6
Finite Difference Approximations of the
Derivatives
7
Finite Difference Approximations
Derive a numerical approximation to the governing
equation, replacing a relation between the
derivatives by a relation between the discrete
nodal values.
The time derivative is found using a forward
EULER method. The approximation can be found by
using a Taylor series.
? t
f(t?t,x)
? t
8
Finite Difference Approximations
Time derivative
Solving this equation for the time derivative
gives
9
Finite Difference Approximations
Spatial first-order derivative
When using FINITE DIFFERENCE approximations, the
values of f are stored at discrete points.
The spatial derivatives of the function are
approximated using a Taylor series.
10
Finite Difference Approximations
Start by expressing the value of f(xh) and
f(x-h) in terms of f(x)
Subtracting the second equation from the first
11
Finite Difference Approximations
The result is
Rearranging this equation to isolate the first
derivative
? Second order central differencing
12
Finite Difference Approximations
Spatial second-order derivative
Repeat the same Taylor expansion as in the
first-order derivative
13
Finite Difference Approximations
Start by expressing the value of f(xh) and
f(x-h) in terms of f(x)
Adding the second equation to the first
14
Finite Difference Approximations
The result is
Rearranging this equation to isolate the second
derivative
? Second order central differencing
15
Finite Difference Approximations
Solving Partial Differential Equations
16
Finite Difference Approximations
A numerical approximation to
is found by replacing the derivatives by the
following approximations
17
Finite Difference Approximations
Using the shorthand notation
gives
18
Finite Difference Approximations
Substituting these approximations into
gives
Solving for the new value and dropping the error
terms yields
19
Finite Difference Approximations
Thus, given f at one time level, f at the next
time level is given by
The value of every point at level n1 is given
explicitly in terms of the values at the level n.
Forward in Time, Central in Space (FTCS)
20
Example A short MATLAB program
The model equation is computed by the spatial
differencing and time integration. Domain
length 2 Two waves (wavelength 1) Periodic
boundary conditions. Compare the numerical
results with the exact solution.
21
one-dimensional advection-diffusion by the FTCS
scheme n21 nstep100 length2.0
hlength/(n-1) dt0.05 D0.05 fzeros(n,1)
yzeros(n,1) exzeros(n,1) time0.0 for i1n,
f(i)0.5sin(2pih(i-1)) end
initial conditions for m1nstep, m, time for
i1n, ex(i)exp(-4pipiDtime)...
0.5sin(2pi(h(i-1)-time)) end
exact solution hold off
plot(f,'linewidt',2) axis(1 n -2.0, 2.0)
plot solution hold on plot(ex,'r','linewidt',
2)pause plot exact solution
yf
store the solution for
i2n-1, f(i)y(i)-0.5(dt/h)(y(i1)-y(i-
1))... D(dt/h2)(y(i1)-2y(i)y(i-1))
advect by centered differences end
f(n)y(n)-0.5(dt/h)(y(2)-y(n-1))...
D(dt/h2)(y(2)-2y(n)y(n-1)) do
endpoints for f(1)f(n)
periodic
boundaries timetimedt end
22
Evolution for U1 D0.05 k1 N21 ?t0.05
Exact
Numerical
23
Accuracy
24
Although the numerical solution shows a
qualitatively similar trend as the analytical
solution, there are significant quantitative
differences. The derivation of the numerical
approximations for the derivatives showed that
the error depends on the size of h and ?t.
First we test for different ?t. Number of
time steps T/ ?t
25
Accuracy
Evolution for U1 D0.05 k1 N21 ?t0.05
Exact
Numerical
26
Accuracy
Repeat with a smaller time-step ?t0.025
Exact
Numerical
27
Accuracy
Repeat with a smaller time-step ?t0.0125
Exact
Numerical
28
Accuracy
Very fine spatial resolution and a small time step
U1 D0.05 k1 N201 ?t0.0005
Exact
Numerical
29
Quantifying the Error I. Order of Accuracy
30
Order of Accuracy
Exact
time0.50
Numerical
n11 E 0.1633 n21 E 0.0403 n41 E
0.0096 n61 E 0.0041 n81 E
0.0022 n101 E 0.0015 n121 E 0.0011 n161
E 9.2600e-04
31
Order of Accuracy
time0.50
Accuracy. Effect of spatial resolution
dt 0.0005 N 11 to N 161
E
32
Order of Accuracy
If the error is of second order
On a log-log plot, the E versus N (or 1/h) curve
should therefore have a slope ?2
33
Order of Accuracy
time0.50
Accuracy. Effect of spatial resolution
dt 0.0005 N 11 to N 161
E
?
34
Why is does the error deviate from the line for
the highest values of N?
Total
Error
Roundoff
Truncation
Number of Steps, N 1/h
35
Quantifying the Error II. Modified Wave Number
36
Modified Wave Number
Order of Accuracy - How mesh refinement
improves the accuracy Modified Wave Number -
How the difference scheme maintains accuracy
at high wave numbers
37
Modified Wave Number
Consider a pure harmonic of period L
The exact derivative is
where k can take on any of the following values
We now ask how the second order central
difference computes the derivative of f
38
Modified Wave Number
Discretize the x axis with a uniform mesh
Finite-difference approximation
Substituting
39
Modified Wave Number
Modified Wave Number
40
Modified Wave Number
Summary
41
Modified Wave Number
42
Enhancing Accuracy - Padé Scheme
2nd Order Central Difference
4th Order Central Difference
4th Order Padé Scheme
Tridiagonal system of equations for
43
Stability
44
Stability
As long as accuracy is reasonable, integration at
larger time steps is more efficient and
desirable. Can we increase the time step
indefinitely? Lets repeat the 1-D
advection-diffusion equation with large time
steps.
45
Stability
Evolution for U1 D0.05 k1 N21 ?t0.2
Exact
Numerical
46
Stability
Instead of decaying as it should, the amplitude
of the numerical solution keeps increasing.
Indeed, if we continued the calculations, we
would eventually produce numbers larger than the
computer can handle. This results in an
overflow or NaN (not a number).
47
Stability
Stability Analysis of 1st Order ODE (Euler Method)
48
Stability
Consider a differential equation
Stability analysis seeks the conditions (of
stepsize h) for which the numerical solution
remains bounded. Taylor Expansion
49
Stability
A Linearized Model Equation
where
For the Euler Method
Solution at time step n
50
Stability
The numerical solution is stable (bounded) if
Stability boundary
If ? is real and negative,
51
Stability
Stability Analysis of Advection-Diffusion
Equation Von Neumann Method
52
Stability
Stability analysis of the full nonlinear system
of equations is in general too complicated. We
study a model problem that in some way mimics the
full equations. The linear advection-diffusion
equation is one such model equation, and we will
apply von Neumann's method to check the stability
of a simple finite difference approximation to
that equation.
53
Stability
Consider the 1-D advection-diffusion equation
In finite-difference form
54
Stability
Substitute
into the finite difference equation
55
Stability
satisfies the equation, we get
Since
Write the error as a wave
Dropping the subscript
56
Stability
The error at node j is
The error at j1 and j -1 can be written as
Substituting and canceling yields
57
Stability
Substituting
into
yields
58
Stability
Dividing by the error amplitude at n
or
59
Stability
The ratio of the error amplitude at n1 and n is
Stability requires that
Since the amplification factor is a complex
number, and k, the wave number of the error, can
be anything, the determination of the stability
limit is slightly involved. We will look at two
special cases (a) U 0 and (b) D 0
60
Stability
(a) Consider first the case when U 0, so the
problem reduces to a pure diffusion
61
Stability
(b) Consider now the other limit where D 0 and
we have a pure advection problem.
Since the amplification factor has the form 1i()
the absolute value of this complex number is
always larger than unity and the method is
unconditionally unstable for this case.
62
Stability
(c) For the general case (Hirsch v.1, pp. 403-408)
Recall amplification factor
Im(G)
Stability Condition
G
Re(G)
?
0
1
2?
63
Stability
In summary, the stability condition for a general
1-D convection-diffusion equation becomes
and
or
Notice that high velocity and low viscosity lead
to instability according to the second
restriction.
64
Stability
For a two-dimensional problem, assume an error of
the form
A stability analysis gives
and
For a three-dimensional problem we get
and
65
Stability
Stability Now you know!
Convergence the solution to the
finite-difference equation approaches the true
solution to the PDE having the same initial and
boundary conditions as the mesh is refined.
Laxs Equivalence Theorem Given a properly posed
initial value problem and a finite-difference
approximation to it that satisfies the
consistency condition, stability is the necessary
and sufficient condition for convergence.
66
Consistency Does the error always go to zero?
67
Consistency
Modified Equation (Finite Difference
Equation) Using the finite difference
approximation, we are effectively solving for
equations that are slightly different from the
original partial differential equations.
Does the finite difference equation approach the
partial differential equation in the limit of
zero ?t and h?
68
Consistency
Consider the 1-D advection-diffusion equation
and its finite-difference approximation
The discrepancy between the two equations can be
found by deriving the modified equation.
69
Consistency
Substituting
into the finite difference equation
70
Therefore, the finite-difference scheme is
effectively solving the above modified
equation. In this case, the error goes to zero as
h ? 0 and ?t ? 0, so the approximation is said
to be consistent. Although most finite
difference approximations are consistent,
innocent-looking modifications can sometimes lead
to approximations that are not!