(Cfd Ebook) - Solution to Two-Dimensional Incompressible Navier-Stokes Equations
In:
Submitted By 96deepa Words 3456 Pages 14
Solution to two-dimensional Incompressible Navier-Stokes Equations with SIMPLE, SIMPLER and Vorticity-Stream Function Approaches. Driven-Lid Cavity Problem: Solution and Visualization.
Maciej Matyka
Computational Physics Section of Theoretical Physics University of Wroclaw in Poland Department of Physics and Astronomy o Exchange Student at University of Link¨ping in Sweden maq@panoramix.ift.uni.wroc.pl http://panoramix.ift.uni.wroc.pl/∼maq May 8, 2003
Abstract
In that report solution to incompressible Navier - Stokes equations in non - dimensional form will be presented. Standard fundamental methods: SIMPLE, SIMPLER (SIMPLE Revised) and Vorticity-Stream function approach are compared and results of them are analyzed for standard CFD test case - Drived Cavity flow. Different aspect ratios of cavity and different Reynolds numbers are studied.
1
Introduction
The main problem is to solve two-dimensional NavierStokes equations. I will consider two different mathematical formulations of that problem: • u, v, p primitive variables formulation • ζ, ψ vorticity-stream function approach I will provide full solution with both of these methods. First we will consider three standard, primitive component formulations, where fundamental Navier-Stokes equation
will be solved on rectangular, staggered grid. Then, solution on non-staggered grid with vorticity-stream function form of NS equations will be shown.
2
Math background
We will consider two-dimensional Navier-Stokes equations in non-dimensional form1 :
1 We
consider flow without external forces i.e. without gravity.
→ ∂− u − − = −(→ )→ − u u ∂t D=
1 ϕ+ Re
2→ −
u
(1) (2)
Guess (P*) , (U*) n , (V*) n n
Solve Momentum eq. for (U*) n+1 , (V*) n+1
→=0 − u
Where equation (2) is a continuity equation which has to be true for the final result.
Solve Poisson eq. for (P’)
If not coverged
3
Primitive variables formulation
P n+1 = (P*) n + (P’)
U n+1 = (U*) n + (U’) V n+1 = (V*) n + (V’)
If coverged
First we will examine SIMPLE algorithm which is based on primitive variables formulation of NS equations. When we say ”primitive variables” we mean u, v, p where u = (u, v) is a velocity vector, and p is pressure. We can rewrite equation (1) in differential form for both velocity components: ∂u2 ∂uv ∂p 1 ∂2u ∂2v ∂u =− + − + ( + 2) ∂t ∂x ∂y ∂x Re ∂x2 ∂y ∂v 2 ∂uv ∂p 1 ∂2u ∂2v ∂v =− + − + ( + 2) ∂t ∂y ∂x ∂y Re ∂x2 ∂y
Visualize Results
Figure 1: SIMPLE Flow Diagram (3)
3.2
(4) 3.2.1
Numerical Methods in SIMPLE
Staggered Grid
We rewrite continuity equation in the following form: ∂u ∂v + =0 ∂x ∂y (5)
These equations are to be solved with SIMPLE method.
3.1
SIMPLE algorithm
SIMPLE algorithm is one of the fundamental algorithm to solve incompressible NS equations. SIMPLE means: Semi Implicit Method for Pressure Linked Equations. Algorithm used in my calculations is presented in the figure (1). First we have to guess initial values of the pressure field2 (P ∗ )n and set initial value of velocity field - (U ∗ )n , (V ∗ )n . Then equation (3) and (4) is solved to obtain values of (U ∗ )n+1 , (V ∗ )n+1 . Next we have to solve pressure-correction equation:
2
For discretization of differential equations I am using staggered grid. In the figure (2) staggered grid for rectangular area is sketched. Primitive variables are placed in different places. In points i, j on a grid pressure P values, in points i + 0.5, j u x-velocity components and in points i, j + 0.5 v y-velocity components are placed. That simple model of staggered grid gives us possibility to use simple discretization with second order accuracy which will be discussed later. p0,0 v0,0.5 u0.5,0 v i,j-0.5 ui-0.5,j pi,j v i,j+0.5 ui+0.5,j
p =
1 ( ∆t
·V)
(6)
After that we use simple relations to obtain corrected values of pressure and velocity fields. Then we can check if solution is coverged.
2 Subscripts denote computational step, where ”n+1” means current step.
Figure 2: Staggered grid: filled circles P , outline circles U x-velocity, cross V y-velocity component.
2
Differential
∂u ∂t
Discretization un+1 −un ∆t
Type forward, O(h) a1 = −
2 n ˙ i+0.5,j−1 (u2 )n (uv)n ˙ i+0.5,j+1 − (uv)n ˙ i+1.5,j − (u )i−0.5,j − 2 · ∆x 2 · ∆y (11)
∂2u ∂x2
ui+1,j −2∗ui,j +ui−1,j (∆x)2
central, O(h2 ) b1 = −
2 n ˙ i−1,j+0.5 (v 2 )n (v u)n ˙ i+1,j+0.5 − (v u)n ˙ i,j+1.5 − (v )i,j−0.5 − 2 · ∆y 2 · ∆x (12) n n un i+1.5,j − 2 · ui+0.5,j + ui−0.5,j 2 (∆x)
∂u2 ∂x
2 u2 i+1,j −ui−1,j (2·∆x)
central, O(h2 )
(a3 ) =
∂p ∂x pi+1,j −pi,j (∆x)
(13) (14) (15) (16)
forward, O(h) (a4 ) =
n n un i+0.5,j+1 − 2 · ui+0.5,j + ui+0.5,j−1 (∆y)2 n n n vi,j+1.5 − 2 · vi,j+0.5 + vi,j−0.5 (∆y)2
∂p ∂y
pi,j+1 −pi,j (∆x)
forward, O(h) (b3 ) = (b4 ) =
Table 1: Discretizations used in SIMPLE algorithm 3.2.2 Discretization Schemes
n n n vi+1,j+0.5 − 2 · vi,j+0.5 + vi−1,j+0.5 (∆x)2
Let us now examine some numerical methods used in presented solution. For algorithm presented in the figure (1) we have only three equations which have to be discretized on a grid. First we have momentum equations (3) and (4). Discrete schemes used in discretization of momentum equations are presented in a table (1). Using presented discrete form of derivatives I obtain numerical scheme for momentum equations exactly in the form presented in [1]. Equations (3) , (4) discretized on staggered grid can be written3 as follows4 :
−1 un+1 = un (pi+1,j − pi,j )) (7) i+0.5,j + ∆t · (A − (∆x) i+0.5,j n+1 −1 vi,j+0.5 = un (pi,j+1 − pi,j )) (8) i,j+0.5 + ∆t · (B − (∆y)
Now we have defined almost everything. Dotted velocities should be also defined. I use simple expressions for it: u = 0.5 · (ui−0.5,j + ui−0.5,j+1 ) ˙ ˙ ˙ u = 0.5 · (ui+0.5,j + ui+0.5,j+1 ) v = 0.5 · (vi,j+0.5 + vi+1,j+0.5 ) ˙ ˙ v = 0.5 · (vi,j−0.5 + vi+1,j−0.5 ) ˙ 3.2.3 Poisson Equation (17) (18) (19) (20)
where A and B are defined as: A = −a1 + (Re)−1 · (a3 + a4 ) B = −b1 + (Re)−1 · (b3 + b4 ) and respectively we define:
3 Please
For equation I use simple iterative procedure. In the figure (3) points used for calculation of pressure at each (i, j) grid points are marked.
P i,j-1
(9)
(10)
P i-1,j
P i,j
P i+1,j
note than cited [1] reference contains some print mistakes
P i,j+1
there. 4 Generally I show there only an idea how to write discretized equations, they should be rewritten with ”*” and ”’” chars for concrete steps of the algorithm
Figure 3: Points on a grid used in iterative procedure for Poisson equation solving. 3
I use simple 4 points scheme for Laplace operator. Directly from [1] expression for one iterative step of poisson equation solver can be written as follows: pi,j = −a−1 (b · (pi+1,j + pi−1,j ) + c · (pi,j+1 + pi,j−1 ) + d) (21) where a = 2∆t 1 1 + ∆x2 ∆y 2 ∆t ∆x2 ∆t ∆y 2 (22) (23) (24)
In the figure (4) I present SIMPLE Revised algorithm. It is easy to extend existing SIMPLE solution to be SIMPLER one. Treating the boundary conditions and numerical methods used in SIMPLER solution is almost the same as in SIMPLE, so I will not repeat myself.
That iterative procedure is rather simple - we use equation (21) for all interior points on a grid. After that one step of iterative procedure is done. Then we check if solution coverges. We can do it simply to check maximum change of pressure on a grid. If it is bigger than we continue iterative process. Solution should finish when pressure field is exactly coverged ( = 0, but in practice I use different value of for different physical properties of simulated models - it will be discussed later.
Vorticity-Stream Function approach to two-dimensional problem of solving Navier-Stokes equations is rather easy. A different form of equations can be scary at the beginning but, mathematically, we have only two variables which have to be obtained during computations: stream vorticity vector ζ and stream function Ψ. First let us provide some definition which will simplify NS equation. The main goal of that is to remove explicitly Pressure from N-S equations. We can do it with the procedure as follows. First let us define vorticity for 2D case: ζ = |ζ| = | ×V|= ∂u ∂v − ∂x ∂y (26)
And stream function definition is: ∂Ψ =u ∂y ∂Ψ = −v ∂x (27) (28)
Solve Poisson eq. for (P n+1 ) using (U^) n+1 , (V^) n+1
Solve Poisson eq. for (P’ ) using (U*) n+1 , (V*) n+1
We can combine these definitions with equations (3) and (4). It will eliminate pressure from these momentum equations. That combination will give us non-pressure vorticity transport equation which in non-steady form can be written as follows:
If not coverged
Having combined equations (26), (27) and (28) we obtain poisson equation for the Ψ variable:
Un+1 = (U*) n + (U’) V n+1 = (V*) n + (V’)
If coverged
2
Ψ=
∂2Ψ ∂2Ψ + = −ζ ∂x2 ∂y 2
(30)
Visualize Results
Figure 4: Flow chart for SIMPLER algorithm. 4
Now we have all definitions and equations which are needed for vorticity-stream solution. We will solve vorticity transport equation, then new values of ζ will be used to solve equation (30).
Differential
∂ζ ∂t
Discretization ζ n+1 −ζ n ∆t
Type forward, O(h)
4.3
Vorticity-Stream function algorithm
Algorithm of solution for VS function solution is simplier than for SIMPLE method. It is sketched in the figure (6).
∂2ζ ∂x2
ζi+1,j −2∗ζi,j +ζi−1,j (∆x)2
central, O(h2 )
Set initial
Solve Vorticity Transport eq.
∂ζ ∂x
ζi+1,j −ζi−1,j (2·∆x)
central, O(h2 )
Solve Poisson eq. for
∂2Ψ ∂x2
Ψi+1,j −2∗Ψi,j +Ψi−1,j (∆y)2
central, O(h )
2
If not coverged
Table 2: Discretizations used in Vorticity-Stream algorithm
Obtain U,V
If coverged
4.1
Non-Staggered Grid
Visualize Results
Instead of using staggered grid in Vorticity-Stream approach, we will place both ζ and Ψ variables in the same place as it is shown in the figure (5)
Figure 6: Algorithm of Vorticity-Stream solution. First we have to set initial values for ζ and Ψ. I arbitrary set these values to 0. Then Vorticity Transport Equation is solved and new ζ n+1 values are obtained. After that simple iterative procedure is applied to solve the poisson equation. Finally, new values of velocities are easily found from (27) and (28) equation.
5
Figure 5: ζ and Ψ variables in non staggered grid. Discretization is straightforward and easier to implement in a non-staggered grid than in a staggered grid for the SIMPLE algorithm.
Two-dimensional Cavity
Driven-Lid
4.2
Discretization
We will use several schemes to discretize differential equation (26). For Poisson equation we will use the same technique which was presented in the SIMPLE algorithm description, so we will not repeat formulas5 .
5 Formulas for poisson equation will be a little bit different but it is rather easy to obtain it by simple discretization of equation (30).
Let us now provide some examples of practical calculation for implemented methods6 . I will show results of DrivenLid Cavity flow - a standard CFD test case to check the solution. Driven Cavity problem is sketched in the figure (7). Upper lid is moving with u velocity. Main goal is to solve NS equations inside the cavity to obtain velocity field (steady state). First of all we have to decide about boundary conditions for both: SIMPLE and VS approaches which will be quite different.
6 In that section also boundary conditions will be provided, because they are specified especially for the given problem.
The same condition is used for other walls and v velocity components.
5.2
Figure 7: Driven Cavity (lid moving with u constants velocity.
Boundary Stream
Conditions
-
Vorticity
5.1
Boundary Conditions - SIMPLE and SIMPLER
In vorticity-stream formulation I use simple first order expressions for ζ derivatives at the wall. First, we have to set Ψ = 0 at all boundaries. Then for NOSLIP boundary walls we use expression (i.e. for j = ny − 1 row): ζi,0 = 2.0 · Ψi,0 − Ψi,1 ∆y 2 (38)
For SIMPLE(R) method we will use BC as follows: First we have to clear pressure values for all boundaries. We use simple expression: ∂p =0 (31) ∂n where n is normal to the wall. It means that for all i = 0 . . . N X − 1 points of a grid we apply: pi,0 = pi,1 and pi,ny−1 = pi,ny−2 (33) (32)
We apply that procedure for upper and lower wall respectively7 . Then we have to take care of velocities. We would like to apply NOSLIP boundaries for Driven Cavity non-moving walls, so we have to zero values of velocities on every wall. First let us make trivial operation: for every j = 0 . . . N Y − 1 set v0,j = 0 and vnx−1,j = 0 (35) (34)
The same work should be done for u velocities, for i = 0 . . . N X − 1 and for j = N Y − 1. Especially for driven cavity problem we also have to set u velocity equal to 1.0 at j = 0 row, which is done in a straightforward way. One problem is to set boundary conditions at other walls, where no velocity grid points are present. We can do it with a simple linear interpolation of near velocities i.e. for u velocity, for every j = 0 . . . N Y − 1 we set:
7 For
corners simple diagonal values are taken, i.e. p0,0 = p1,1
6
6
Results
In that section some numerical results of calculations with three different techniques will be presented. Since results of calculations are the same I will try to show and compare differences between methods (accuracy, convergence speed). Please note that all comments are under figures.
6.1
Vorticity-Stream, Driven Cavity, Re = 500, Grid: 40x40
Figure 8: Streamlines plot for driven cavity with Re = 500 and 1 : 1 aspect ratio, grid size 40x40. Two vortexes are found in the corners of the Cavity, computed with the Vorticity-Stream approach. Solution visualized with Streamline plot technique.
7
6.2
Different Visualization Techniques, Driven Cavity, Re = 500, Grid: 40x40
Stream Function Distribution
Vorticity Function Distribution
Stream Function Contour Plot
Red - U velocity Green - V velocity
Figure 9: There are presented different types of visualizations generated by my solver. Computations as above Re = 500 and other parameters are the same. (That is only a part of possibility visualizations, more will be available on my web page soon).
8
6.3
SIMPLE, SIMPLER, Driven Cavity, Re = 100, Grid: 21x21
40x40, Aspect Ratio 1:1
30x60, Aspect Ratio 1:2
40x60, Aspect Ratio 2:3
Figure 10: Streamlines plot for SIMPLE (and SIMPLER - because they are the same) calculation of driven cavity with Re = 100 and different grid sizes and aspect ratios.
6.4
Convergence for SIMPLE and Vorticity-Sream algorthms
9
1.5 1.4 1.3 convergence 1.2 1.1 1 0.9 0.8
SIMPLE
0
10000
20000
30000 40000 time step
50000
60000
70000
|v n+1 | |v n |
Figure 11: That figure shows how convergence changes during iteration steps. On y axis we have
variable.
1.5 1.4 1.3 convergence 1.2 1.1 1 0.9 0.8
Vorticity-Stream
0
10000
20000
30000 40000 time step
50000
60000
70000
|v n+1 | |v n |
Figure 12: That figure shows how convergence changes during iteration steps. On y axis we have
Figure 13: Convergence test for three solution algorithms. A lot of problems appeared there. Convergence speed depends on a lot of things, so for different properties of calculation (Reynolds numbers, spatial grid resolution, poisson equation accuracy etc.) different results appears. That results computed for Re = 300 and grid 30x30 shows that Vorticity-Stream function solver converge faster than SIMPLER and SIMPLE. Anyway - more carefully study should be made there to make sure about that results. On the y axis we have |v n−1 − v n | convergence coefficient.
11
7
Calculation For Flows over Obstacles
In that section I present some calculations made to test my SIMPLE solver for problems other than Driven Cavities. There were some problems with boundary conditions and still more work is needed there, but fortunately results are really nice.
Figure 14: Flow of Incompressible fluid over different types of obstacles. Calculations made for Re = 350 and grid size 90x45.
12
8
Conclusion
I have developed three different methods for calculation of incompressible fluid flow. Tests for simple Driven Cavity problem were made. I compared convergence speed for all the methods and it seems that convergence speed depends on problem formulation and physical properties of simulated system. For future I will try to concern more on how to treat boundary conditions for both - pressure based and vorticity-stream function methods. Also some kind of user-friendly software will be released in near future. I would like to thank Grzegorz Juraszek for English language checking.
References
[1] John D. Anderson, Jr. ’Computational Fluid Dynamics: The Basics with Applications’, McGraw-Hill Inc, 1995. Ryszard Grybos, ’Podstawy mechaniki plynow’ (Tom 1 i 2), PWN 1998. David Potter ’Metody obliczeniowe fizyki’, PWN 1982. James D. Bozeman, Charles Dalton, ’Numerical Study of Viscous Flow in Cavity’, Journal of Computational Physics, vol. 12, 1973. J.C. Tannehill, D.A. Anderson, ’Computational Fluid Mechanics and Heat Transfer, Second Edition’, Series in Computational and Physical Processes in Mechanics and Thermal Sciences . C.A.J. Fletcher, ’Computational Techniques for Fluid Dynamics, Volume 2’, Springer . F. H. Harlow, John P. Shannon, ’The Splash of Liquid Drop’, Journal of Applied Physics (vol. 38, n.10 Sept. 1967). J. Welch, F. Harlow, J. Shannon, ’The MAC Method’, Los Alamos Scientific Laboratory of the University of California (1965). N.Foster, D.Metaxas, ’Realistic Animation of Liquids’, Center for Human Modeling and Simulation.
