Download Solving Navier-Stokes: Finite Difference & Cylindrical-Polar Approach - Prof. George J. Hi and more Study notes Chemistry in PDF only on Docsity!
Chapter 12 Numerical Simulation
The classical exact solutions of the Navier-Stokes equations are valuable for the insight they give and as a exact result that any approximate method should be able to duplicate with to an acceptable precision. However, if we limit ourselves to exact solutions or even approximate perturbation or series solutions, a vast number of important engineering problems will be beyond reach. To meet this need, a number of computational fluid mechanics (CFM) codes have become commercially available. A student could use the CFM code as a virtual experiment to learn about fluid mechanics. While this may be satisfactory for someone who will become a design engineer, this is not sufficient for the researcher who may wish to solve a problem that has not been solved before. The user usually has to treat the CFM code as a "black box" and accept the result of the simulation as correct. The researcher must always be questioning if any computed result is valid and if not, what must be done to make it valid. Therefore I believe Ph.D. students should know how to develop their own numerical simulators rather than be just a user. We will start with a working FORTAN and MATLAB code for solving very simple, generic problems in 2-D. The student should be able to examine each part of the code and understand everything with the exception of the algorithm for solving the linear system of equations. Then the student will change boundary conditions, include transient and nonlinear capabilities, include curvilinear coordinates, and compute pressure and stresses.
The Stream Function – Vorticity Method Two-dimensional flow of an incompressible, Newtonian fluid can be formulated with the stream function and vorticity as dependent variables.
( )
2 2 2 2 2 2 2 2
where
( , , 0) , , 0 , 0, 0,
x y
x y
y (^) x
w x y dw w w w w w v v dt t x y x y
v v y x v (^) v w x y
∂ ∂ ∂ æ^ ∂ ∂ ö = + + = (^) ç + ÷ ∂ ∂ ∂ (^) è ∂ ∂ ø
æ (^) ∂ ∂ ö = = ∇ × = (^) ç − (^) ÷ = è ∂^ ∂ ø æ ∂ ∂ ö = = (^) ç − ÷ è ∂^ ∂ ø
v A A
w
We will begin developing a generic code by assuming steady, creeping (zero Reynolds number) flow with specified velocity on the boundaries in a rectangular domain. The PDE are now,
2 2 2 2
2 2 2 2
, 0 , 0 , Re 0
0
x y
w x y x L y L w w x y
This is a pair of linear, elliptic PDEs. The boundary conditions with specified velocity are
specified specified ( )
normal tangential
BC boundary BC (^) boundary
O U
ds
w
= ∇ ×
ò
v v
v n
v
The boundary conditions at a plane of symmetry are
0 0 =constant 0
n
w
n v v
The PDE and definitions are made dimensionless with respect to reference quantities such that the variables are of the order of unity.
2 2 2 2 2 2
2 2 2 2 2
x x y y
o o
o o x o y o
o o x y y x
x y L x y L L L w w U w
w L w x L w y x y
w w x y
v v U L y U L x
é ù (^) ∂ é^ ù∂ ê ú +^ ê^ ú = − ë û ∂^ êë úû∂ < < < <
∂ ∂
- = ∂ ∂ é ù (^) ∂ é ù∂ = (^) ê ú = −ê ú êë úû ∂^ ë û∂
v v
The unit square is now discretized into JMAX by JMAX grid points where the boundary conditions and dependent variables will be evaluated. The grid
spacing is δ = 1/( JMAX − 1). The first and last row and column are boundary
values. The partial derivatives will be approximated by finite differences. For example, the second derivative of vorticity is discretized by a Taylor's series.
2 2 3 3 4 4 (^1 2 3 )
2 2 3 3 4 4 (^1 2 3 )
2 2 4 4 1 1 2 2 4 4
i i i (^) i i x
i i i (^) i i x
i i i i x
w w w w w w x x x x
w w w w w w x x x x
w w w w w w x x x
δ δ δ δ
δ δ δ δ
δ δ
−
− +
æ ∂ ö^ æ^ ∂ ö^ æ^ ∂ ö^ æ^ ∂ ö = + (^) ç ÷+ (^) ç ÷ + (^) ç ÷ + ç ÷ è ∂^ ø (^) è ∂^ ø è ∂^ ø è ∂ ø
æ ∂^ ö æ^ ∂^ ö^ æ^ ∂^ ö^ æ^ ∂ ö = − (^) ç ÷+ (^) ç ÷ − (^) ç ÷ + ç ÷ è ∂^ ø (^) è ∂^ ø è ∂^ ø è ∂ ø
æ (^) ∂ ö (^) − + æ (^) ∂ ö æ∂ ç ÷ =^ −^ ç ÷ + è ∂^ ø è ∂^ ø è∂
1 1 2 2
x wi wi wi (^) O δ δ
− +
é (^) ö ù ê (^) ç ÷ ú êë (^) ø úû − + = +
The finite difference approximation to the PDE at the interior points results in the following set of equations.
Finite difference grid for discretizing PDE (grid point formulation)
1 2 3 … JM1 JMAX
JMAX
JM
xi
yj
2 2 2 2 1, 1, , 1 , 1 , , 2 2 2 1, 1, , 1 , 1 ,
i j i j i j i j i j i j i j i j i j i j i j
w w w w w w i j JMAX
= L −
The vorticity at the boundary is discretized and expressed in terms of the components of velocity at the boundary, the stream function values on the boundary and a stream function value in the interior grid. (A greater accuracy is possible by using two interior points.) The stream function at the first interior point ( i =2) from the x boundary is written with a Taylor's series as follows.
( )
( )
(^1) ( )
2 2 3 (^2 1 ) (^1 ) 2 2 1 2 2 1 1 2 1 2 1 2 2 1 2 2
BC y BC BC y x
BC x
BC (^) BC BC (^) x y
O
x x
O
x x
v O
v (^) v w x y v x y
v v O y
æ ∂^ ö^ æ^ ∂ ö = ± (^) ç ÷ + (^) ç ÷ + è ∂^ ø (^) è ∂ ø æ ∂ ö − (^) æ ∂ ö ç ÷ =^ ç ÷ + è ∂^ ø è^ ∂ ø − = ± +
æ ∂ ∂ ö ≡ (^) ç − ÷ è ∂^ ∂ ø æ ∂ ö ∂ = − (^) ç ÷− è ∂^ ø ∂ − (^) ∂ = − − + ∂
m
m
The choice of sign depends on whether x is increasing or decreasing at the boundary. Similarly, on a y boundary,
( )
2 2 1 (^1 )
BC BC BC wBC vx^ vy O x
The boundary condition on the stream function is specified by the normal component of velocity at the boundaries. Since we have assumed zero normal component of velocity, the stream function is a constant on the boundary, which we specify to be zero. The stream function at the boundary is calculated from the normal component of velocity by numerical integration using the trapezodial rule, e.g.,
1 (^ )^1 (^ ) / j j (^) 2 x (^) j x j v v
ψ = ψ − + é ë^ −+ ùû α
The coefficients for the interior grid points adjacent to a boundary are modified as a result of substitution the boundary value of stream function or the linear equation for the boundary vorticity into the difference equations. The stream function equation is coupled to the vorticity with the e ij coefficient and the vorticity equation is coupled to the stream function through the boundary conditions. For example, at a x = 0 boundary, the coefficients will be modified as follows.
1, 2 1, 2
BC j BC BC BC (^) x j v
i f i j f i j b i j
v f i j f i j b i j v j y e i j e i j b i j
If the boundary condition is that of zero vorticity as on a surface of symmetry, it is necessary to only set to zero the coefficient that multiplies the boundary value of vorticity, i.e.,
2 ( , , 2, 2) 0, when (1, ) 0
i b i j w j
The linear system of equations and unknowns can be written in matrix form by ordering the equations with the i index changing the fastest, i.e., i = 2, 3, …, JMAX -1. The linear system of equations can now be expressed in matrix notation as follows.
2, 2, 2, 3, 3, 3,
1, 1 1, 1 1, 1
where
JM JM JM JM JM JM
w
w
w
é ù ê ú é ù ê^ ú ê ú ê^ ú ê ú ê^ ú = = (^) ê ú ê ú (^) ê ú ê ú (^) ê ú êë úû ê ú ê ú êë úû
A x F
u u x
u
M
M
The non-zero coefficients of coefficient matrix A for JMAX = 5 is illustrated below. Since the first and last grid points of each row and column of grid points are boundary conditions, the number of grid with pairs of unknowns is only 3×3. Only one grid point is isolated from boundaries. The i,j location of the coefficients become clearer if a box is drawn to enclose each 2×2 coefficient.
Most of the elements of the coefficient matrix are zero. In fact, the coefficient matrix is a block pentadiagonal sparse matrix and only the nonzero coefficients need to be stored and processed for solving the linear system of equations. The non-zero coefficients of the coefficient matrix are stored with the row-indexed sparse storage mode and the linear system of equations is solved by preconditioned biconjugate gradient method (Numerical Recipes, 1992). The row-indexed sparse matrix mode requires storing the nonzero coefficients in the array, sa(k), and the column number of the coefficient in the coefficient matrix in the array, ija(k), where k = 1, 2, …. The indices that are needed to identify the coefficients are the i,j grid location ( i,j = 2, 3, …, JMAX -1), the m equation and dependent variable identifier ( m = 1 for stream function; = 2 for vorticity), and the IA row index of the coefficient matrix ( IA = 1, 2, 3, 4, …, 2( JMAX -2)^2 ). The coefficient matrix has two equations for each grid point, the first for the stream function and the second for vorticity. The total number of equations is NN = 2( JMAX -2)^2. The diagonal elements of the coefficient array are first stored in the sa(k) array. The first element of ija(k), ija(1) = NN +2 and can be used to determine the size of the coefficient matrix. The algorithm then cycles over the pairs of rows in the coefficient matrix while keeping track of the i,j locations of the grid point and cycling over equation m = 1 and 2. The off-diagonal coefficients are stored in sa(k) and the column number in the coefficient matrix stored in ija(k). As each row is completed, the k index for the next row is stored in the first NN elements of ija. A test problem for this code is a square box that has one side sliding as to impart a unit tangential component of velocity along this side. All other walls are
0 5 10 15
0
2
4
6
8
10
12
14
16
18
nz = 83
Assignment 12.1 Use the code in the /MatlabScripts/numer/ directory of the CENG501 web site as a starting point. Add the capability to include arbitrarily specified normal component of velocity at the boundaries. The stream function at the boundary should be numerically calculated from the specified velocity at the boundary. Test the code with Couette and Plane-Poiseuille flows first in the x and then in the y directions. Display vector plot of velocity and contour plots of stream function and vorticity. Show you analysis and code.
Assignment 12.2 Form teams of two or three and do one of the following additions to the code. Find suitable problems to test code.
- Calculate pressure and stress. Test problems: Couette and Plane- Poiseuille flows.
- Calculate transient and finite Reynolds number flows. Test problems: Plate suddenly set in motion and flow in cavity.
- Add option for cylindrical polar coordinates. Use the coordinate transformation, z = ln r. Test problems: line source and annular rotational Couette flow.
Form new teams of two or three with one member that developed each part of the code. Include all features into the code. Test the combined features. Work together as a team to do the following simulation assignments.
0 5 10 15 20 25
0
5
10
15
20
25
Velocity Arrays
Velocity vectors for square with top side sliding
Assignment 12.3 Boundary layer on flat plate. Compute boundary layer velocity profiles and drag coefficient and compare with the Blasius solution. What assumption is made about the value of the Reynolds in the Blasius solution?
Assignment 12.4 Flow around cylinder. Assume potential flow far from the cylinder. Calculate drag coefficient as a function of Reynolds number and compare with literature.
Assignment 12.5 Flow around cylinder. Find the Reynolds number at which boundary layer separation occurs.
Assignment 12.7 Flow around cylinder. Find the Reynolds number for the onset of vortex shedding.
FORTRAN on Owlnet
Modern FORTRAN compilers in the MSWindows are very user friendly. They have built in project workspace routines, documentation, and integrated editor and debugger. If you use the FORTRAN compiler on Owlnet, it is useful to know how to use make files to compile and link FORTRAN code. The following is an example make file. You may give it the name, makefile , with no extension.
f77 -c streamf1.f; f77 -c bc1.f; f77 -c coef1.f; f77 -c asolve.f; f77 -c atimes.f; f77 -c dsprsax.f; f77 -c dsprstx.f; f77 -c linbcg.f; f77 -c snrm.f; f77 -o exe streamf1.o bc1.o coef1.o asolve.o atimes.o dsprsax.o dsprstx.o linbcg.o snrm.o
You have to give yourself permission to execute the make file with the one-time command,
chmod +x makefile
The line with the –c option compiles the source code and makes an object code file with the extension of .o. The line with the – o option links the object files into the executable file called exe. This executable file can be executed by typing its name, exe , on the command line.
Transients and Finite Reynolds Number
The development of numerical simulators should start with a problem for which exact solutions exist so that the numerical algorithm can be verified. Steady, creeping flow has many exact solutions and thus it was used as the starting point for developing a numerical simulator for the Navier-Stokes equations. The student should verify that their code is able to duplicate exact solutions to any desired degree of accuracy. The next stage in the development of the code is to add transients and finite Reynolds numbers. Finite Reynolds number results in nonlinear terms in the Navier-Stokes equation. Algorithms for solving linear systems of equations such as the conjugate gradient methods solves linear systems in one call to the routine. Nonlinear terms need to be iterated or lagged during the transient solution. Thus nonlinear steady state flow may be solved as the evolution of a transient solution from initial conditions to steady state. The transient and nonlinear terms now need to be included when making the vorticity equation dimensionless.
2 2 2 2
2 2
x y
o o o o x y x x x
dw w w w w w v v dt t x y x y t t t w U t w U t w t w w v v t L x L y L x y
μ ρ
μ α α ρ
∂ ∂ ∂ é ∂ ∂ ù = + + = (^) ê + ú ∂ ∂ ∂ (^) ë ∂ ∂ û
=
∂ é^ ù^ ∂ é^ ù^ ∂ é^ ù é ∂ ∂ ù
- (^) ê ú + (^) ê ú = (^) ê ú ê + ú ∂ (^) ë û ∂ (^) ë û ∂ (^) ë ûë ∂ ∂ û
Two of the dimensionless groups can be set to unity by expressing the reference time in terms of the ratio of the characteristic velocity and characteristic length. The remaining dimensionless group is the Reynolds number.
2
Re
o x o x
o x x
U t L t L U
t L U L
é ù ê ú=^ Þ^ = ë û é ù é ù ê ú =^ ê ú= ë û ë û
The dimensionless equation for vorticity with the * dropped is
2 2 2 x y 2 2
w w w w w Re v v t x y x y
é (^) ∂ ∂ ∂ ù ∂ ∂ ê +^ +^ ú=^ + ë ∂^ ∂^ ∂^ û ∂^ ∂
The convective derivative can be written in conservative form by use of the continuity equation.
,^ ,^ ,
,
, (^) ,
2 2 2 2 2
but 0, for incompressible flow thus
or
i j (^) i i i j i j i
i i
i j i i j (^) i
x y
v w v w v w
v
v w v w
w v w v w w w Re t x y x y
é (^) ∂ ∂ ∂ ù ∂ ∂ ê +^ +^ ú=^ + ë ∂^ ∂^ ∂^ û ∂^ ∂
v w v w
We will first illustrate how to compute the transient, linear problem before tackling the nonlinear terms. Stability of the transient finite difference equations is greatly improved if the spatial differences are evaluated at the new, n+1 , time level while the time derivative is evaluated with a backward-difference method.
1 , , ,
1 , , ,
n n i j i j i j
n n i j i j i j
w w^ w^ w O t O t t t t w w w
∂ −^ ∆
Recall that the finite difference form of the steady state equations were expressed as a system of equations with coefficients, a, b , , ect. The vorticity equation for the i,j grid point will be rewritten but now with the transient term
included. It is convenient to solve for ∆ w rather than w n+1^ at each time step.
2 1 1 1 1 1 1 1, 1, , 1 , 1 , ,
1, 1, , 1 , 1 ,
n n n n n Re n i j i j i j i j i j i j n n n n n i j i j i j i j i j
a w b w c w d w e w w t f a w b w c w d w e w
We can see that this equation is of the same form as before with only the coefficients e and f of the vorticity equation modified.
2
1, 1, , 1 , 1 ,
' (^) in j in j i jn^ i jn^ i jn
Re e e t f f a w b w c w d w e w
Since the variable in the linear system of equations is now ∆ w rather than w , the
equations for the stream function also need to be modified.
The choice of the finite difference expression for the convective derivative is very important for the stability of the solution. We will illustrate here treatment of only the x derivative. Suppose i-1 , i , and i+ 1 are the grid points where the stream function and vorticity are evaluated and i – ½ and i+ ½ are the mid-points where the velocity-vorticity products are evaluated.
u w^ (^ u w^ )^ i 1/ 2 (^ u w ) i 1/ 2 x x
The dependent variables are known only at the grid points and interpolation will be needed to evaluate in between location. If the product is evaluated at i – ½ and i+ ½ by using the arithmetic average of the values at the grid points on either side, the numerical solution will tend to oscillate if the convective transport dominates the diffusive (viscous) transport. It is preferable to determine the upstream direction from the sign of u (^) i+1/2 or u (^) i-1/2 and use the vorticity from the upstream grid point. This is illustrated below.
u (^) i-1/2 = 0.5*(u(i,j)+u(i-1,j)) if u (^) i-1/2 >0 then (uw)i-1/2 = u (^) i-1/2 wi-1 ; b'=b+δ Re u (^) i-1/ else (uw)i-1/2 = u (^) i-1/2 wi ; e'=e+δ Re u (^) i-1/ endif
u (^) i+1/2 = 0.5*(u(i,j)+u(i+1,j)) if u (^) I+1/2 >0 then (uw)i+1/2 = u (^) i+1/2 wi ; e'=e - δ Re u (^) i+1/ else (uw)i+1/2 = u (^) i+1/2 wi+1 ; a'=a - δ Re u (^) i+1/ endif
The y-direction will be similar except the aspect ratio, α, must be included. These modifications to the coefficients must be made before the coefficients are updated for the boundary conditions.
Calculation of Pressure
The vorticity-stream function method does not require calculation of pressure to determine the flow field. However, if the force or drag on a body or a conduit is of interest, the pressure must be computed to determine the stress field. Here we will derive the Poisson equation for pressure and determine the boundary conditions using the equations of motion. If the pressure is desired only at the boundary, it may be possible to integrate the pressure gradient determined from the equations of motion. The dimensionless equations of motion for incompressible flow of a Newtonian fluid are as follows.
2
2
where
d dt t t
P Re
U P
t t P L U L L U U L P p g z Re
v v v v v v v
v
v x v x
Here all coordinates were scaled with respect to the same length scale. The aspect ratio must be included in the final equation if the coordinates are scaled with respect to different length scales. We will drop the * and use x and y for the coordinates and u and v as the components of velocity. The following derivation follows that of Hoffmann and Chiang (1993). The equations of motion in 2-D in conservative form is as follows.
2 2 2 2 2 2 2 2 2 2
u u u v P u u t x y x Re x y
v u v v P v v t x y y Re x y
∂ ∂ ∂ ∂ æ ∂ ∂ ö
- = − + (^) ç + ÷ ∂ ∂ ∂ ∂ (^) è ∂ ∂ ø
∂ ∂ ∂ ∂ æ^ ∂ ∂ ö
- = − + (^) ç + ÷ ∂ ∂ ∂ ∂ (^) è ∂ ∂ ø
The Laplacian of pressure is determined by taking the divergence of the equations of motion. We will carry out the derivation step wise by first taking the x -derivative of the x -component of the equations of motion.
2 2 2 2 2 2 2
u u u u u v v v u u P u u v u t x x x x x y x y x y x y x Re x
∂ æ ∂ ö ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ç ÷+^ +^ +^ +^ +^ +^ = −^ +^ ∇ ∂ (^) è ∂ (^) ø ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
Two pair of terms cancels because of the equation of continuity for incompressible flow.
2 2 2
(^2 2 2 ) 2 2 2
Thus 1 Re
u u v x x y v u u v u u u x y x x x y
u u u v u u P u v u t x x x x y x y x x
∂ æ ∂ ∂ ö ç +^ ÷= ∂ (^) è ∂ ∂ ø ∂ ∂ ∂ æ ∂ ∂ ö
- = (^) ç + (^) ÷= ∂ ∂ ∂ ∂ (^) è ∂ ∂ ø
∂ (^) æ ∂ (^) ö æ ∂ (^) ö ∂ ∂ ∂ ∂ ∂ ∂ ç ÷ +^ ç ÷ +^ +^ +^ = −^ +^ ∇ ∂ (^) è ∂ (^) ø è ∂ (^) ø ∂ ∂ ∂ ∂ ∂ ∂ ∂
2 2 2 2 2 2 2 2 2 2 2
P P
x y x y x y
∂ ∂ é^ ∂ ψ ∂ ψ æ ∂ψ öù
- = ê − (^) ç ÷ú ∂ ∂ (^) ê ∂ ∂ (^) è ∂ ∂ øú ë û
This is the equation when the coordinates have the same reference length. If the coordinates are normalized with respect to different lengths, then the equation is as follows.
2 2 2 2 2 2 2 2 2 2 2 2 2
P P
x y x y x y
ψ ψ ψ α α
∂ ∂ é^ ∂ ∂ æ ∂ ö ù
- = ê − (^) ç ÷ ú ∂ ∂ (^) ê ∂ ∂ (^) è ∂ ∂ ø ú ë û
The Poisson equation for pressure needs to have boundary conditions for solution. The equations of motion give an expression for the pressure gradient. The normal derivative of pressure at the boundary can be determined for the Neumann-type boundary condition. The equations of motion usually simplify at boundaries if the boundary has no slip BC or if it has parallel flow in the direction either parallel or normal to the boundary. For example, for no-slip
2 2 2 2 2
, at , 0
, at , 0
P v y c u v x Re y P u x c u v y Re x
α
α
If the flow is parallel and in the direction normal to the boundary
2 2 2 2 2
P u y c v u u y x Re y P v x c u v v x y Re x
α
α
Cylindrical-Polar Coordinates
The code for the numerical solution of the Navier-Stokes equations in Cartesian coordinates can be easily modified to cylindrical-polar coordinates. The coordinate transformation is first illustrated for the Laplacian operator.
2 2 2 2 1 2
r , r r r , (^0) o r r r r
∂ (^) æ ∂ (^) ö ∂ ∇ = (^) ç ÷+ ≤ ≤ ≤ ≤ ∂ (^) è ∂ (^) ø ∂
The independent variables are made dimensionless with respect to the boundary parameters.
o
r r r r r r
The Lapacian operator with the dimensionless coordinates after dropping the * is now,
2 2 2 2 2 2 1
o
r r r r r r r
é (^) ∂ (^) æ ∂ (^) ö ∂ ù ∇ = (^) ê ç ÷+ (^) ú ≤ ≤ ≤ ≤ ë ∂^ è^ ∂^ ø ∂ û
The radial coordinate is transformed to the logarithm of the radius.
( ) ( )
2 1
2 2 2 2
ln 1 ln , 0 1, ln ln / exp /
r z r z r r r z dz r dr z r z
r r r r r z
γ γ β γ γ
ψ γ ψ
∂ (^) æ ∂ (^) ö ∂ ç ÷= ∂ (^) è ∂ (^) ø ∂
The Laplacian operator is now as follows,
2 2 2 2 2 2 2 2 2 2 1
o
z r r z r
é ∂ ∂ ù ∇ = (^) ê + (^) ú ≤ ≤ ≤ ≤ = ë ∂^ ∂ û
The finite difference expression for the Laplacian will be same as that for
Cartesian coordinates with (^) ( z , θ (^) )substituted for ( x,y ) except for γ 2 and r^2 factors.
( ) ( )
i j
i j
z i j JMAX
JMAX
z i j
K
The curl operator is modified from that in Cartesian coordinates.
