METU Mechanical Engineering Department ME 485 Computational Fluid Dynamics using Finite Volume Method Spring 2014 (Dr. Sert) Demonstration Demonstration of How SIMPLE Algorithm Works on 1D Co-located C o-located Meshes Problem Definition
3 1 kg/m 0 Pa =2m
Simplify the incompressible flow in the following converging nozzle to be 1D and inviscid. Obtain the f inite volume
1 m/s
solution using 5 cells of equal length. Density of the fluid is and exit pressure is specified as .
= 1 kg/m3
== 01.5mm/s2
. As boundary conditions inlet v elocity is given
==00.1Pma2
Cross sectional area of the nozzle decreases linearly
=0.50.2
/
Exit pressure is set to zero for simplicity. Actually the value of this reference pressure is not important for an incompressible flow, because in INS only the pressure gradient ( ) is important, not the actual pressure values. If the actual exit pressure is , pressure values that we will obtain o btain by setting exit pressure to zero will be gage pressures. To obtain absolute values we can add to all the pressure values.
Analytical Solution
Analytical solution of the problem is governed by the Bernoulli equation.
2 =constant ̇ = =0.5kgs⁄ = ̇ =5ms⁄
Mass flow rate inside the nozzle is constant and its value is known due to the given inlet speed
Using this value we can calculate the speed at any point inside the nozzle. Exit speed is
1
Using the known speed and pressure at the exit the constant of the Bernoulli equation can be calculated as
̇ =0. 5
With the determined.
2 =0 152 =12.5 Pa
equation and the above equation speed and pressure at any point inside the nozzle can be
Nodes and faces of the 5 cell mesh are shown below
∆=0.4 m 1
2
4
3
5
Analytical solution at the faces and nodes is as follows
Face
Node
m2 m/s Pa
m 0.0
1
2
3
4
5
0.50
1.0000
12.0000
0.2
0.46
1.0870
11.9093
0.4
0.42
1.1905
11.7914
0.6
0.38
1.3158
11.6343
0.8
0.34
1.4706
11.4187
1.0
0.30
1.6667
11.1111
1.2
0.26
1.9231
10.6509
1.4
0.22
2.2727
9.9174
1.6
0.18
2.7778
8.6420
1.8
0.14
3.5714
6.1224
2.0
0.10
5.0000
0.0000
Discretization of the x-Momentum Equation
Consider the following cell P with W and E neighbors
W
P
E
For cell P, the discretized x-momentum equation without the viscous and source terms is
= =
= | ∆∀
1
where , are known, calculated using initial guesses or previous iteration values. and of Eqn (1) can be expressed in terms of speeds at nodes using various schemes. Here upwind scheme is used as follows 2
>0 i f ,0 = if <0 → = max,0 max >0 i f ,0 = if <0 → = max,0 max Pressure derivative of Eqn (1) is discretized as
Volume of cell P is
where
| = 2∆ ∆∀=∆
is the cross sectional area at node P.
Substituting these details into Eqn (1) we get
,0 max,0 max,0= 2 max,0 max 2
which can be arranged as
where
Eqn (3) can also be written as
where
= 2 =max,0 = max,0 = max,0 max,0 =0 =̂ 2 and = ̂ = 1 ∑
3
4 5
Although the source term is zero, it is kept in the equations because for boundary cells there may be nonzero contributions to it.
3
Modification of the x-Momentum Equation for Boundary Cells
Cell 1:
P
E 2
1
= = =known
At the west face inlet velocity is known, i.e. in the x-momentum equation flux at the west face is known
This can be taken to the right hand side of the equation to act as a source term. For the pressure derivative one-sided difference can be used instead of central differencing
| = ∆ = =0 = max,0 = max ,0 =
With these, Eqn (3) for cell 1 becomes (modified terms are shown in red)
where
Therefore at an inlet boundary where velocity is given, following changes occur in the x-momentum equation
Coefficient of the ghost neighbor is set to zero. coefficient changes because there is no contribution from the ghost neighbor.
Momentum flux created by the known inlet velocity acts as a source term. Pressure discretization becomes one-sided (not central).
Cell 5: W 4
P
5
= Considering that the right end of the domain is an exit boundary we can use
= | = ∆ = ∆ 2 = 2 2∆
Pressure gradient term can be discretized to make use of the given
value.
4
With these, Eqn (3) for cell 1 becomes (modified terms are shown in red)
= 2 2∆ = max,0 =0 = max,0 max,0
where
Therefore at an exit boundary where pressure is given, following changes occur in the x-momentum equation
Coefficient of the ghost neighbor is set to zero.
Pressure discretization uses the known exit pressure.
Face Velocity Calculation using Rhie-Chow (Momentum) Interpolation and Relaxation
In total, there are 6 faces in the mesh. Consider the following face
L
with neighboring cells L (left) and R (right).
R
Using Rhie-Chow interpolation and velocity under-relaxation, velocity at face is calculated as
where
=̂ 1 ̂ ̂ ̂ = 2 , = 2 ̂
is the velocity under-relaxation factor and
Eqn (6) and the above expressions for
Modification at face 1:
6
is the face velocity of the previous iteration.
and
need to be modified at the boundary faces.
=
̂
Face 1 at the left boundary has a specified inlet velocity, so we do not use Eqn (6) at face 1, i.e. we do not need or .
5
Modification at face 6: 5
4
= ̂
There is no cell on the right of face 6. Instead of central interpolation, one-sided interpolation can be used to calculate and . With the assumption of constant cell size
Also
̂ =̂ ̂ 2̂ , = 2 1/2
term of Eqn (6) can be expressed in terms of the known exit pressure as
which comes from
≈ ∆/−
Pressure Correction (PC) Equation
PC equation is
′ ′ =∗ ∗ ′ =′ =′ ′ and ′ =′ =′ ′ ′ ′ ′ =∗ ∗ = = =
Relating velocity corrections to pressure corrections as follows
PC equation becomes
7
Eqn (7) is modified as follows for boundary cells.
Modification for cell 1:
P 1
E 2
= ′=0 =0 , =
At west face inlet velocity is specified and
. Due to this following changes happen
6
Modification for cell 5: W
5
4
′
P
=
′′ ′ ′ ′ ′ = 1/2
expression need to be modified because there is no
where
′ =0
. Instead of using
we can use half cell as
because exit pressure is fixed. Due to this coefficients of the PC change as follows
=0 , =2
Face Velocity Corrections
L
R
For the above face , velocity correction is done as follows
where
∗
=∗ ′ ′
8
is the velocity calculated previously using Rhie-Chow interpolation. For the boundary faces Eqn (8) is
used as follows Modification for face 1: Inlet velocity is given and no correction is done. Modification for face 6: Similar to the previous step we use
′ ′ ∗ = 1/2
Correct Cell Center Velocities
W
where ′ =0
P
′ ′ ∗ = 2
E
9
Eqn (9) will be modified for the boundary cells. 7
Modification for cell 1: West cell does not exist. Pressure correction difference can be done in a one-sided way.
′ ′ ∗ = 1
Modification for cell 5: East cell does not exist. Pressure correction difference can be done in a one-sided way using the fact that exit pressure is fixed.
′ ′ ∗ = 1/2
Correct Pressures
where
∗
where ′ =0
=∗ ′
is the pressure of the previous iteration (or the initial guess) and
10 is the pressure relaxation factor.
SIMPLE Iterations STEP 1:
As initial guess we can use the inlet velocity and exit pressure at all nodes and faces.
= = = = =1.0 = = = = = =1.0 = = = = =0.0
ITERATION 1:
∗ Knowns: = =1.0 , = =1.0 , =0.5 , =0.42 , = =0.0 =0.0 , =0.42 , =0.0 , =0.5 ̂ = 0.50.0.42 0 =1.1905 , = 0.0.4462 =1.0952 0.42∗ =0.50.460.00.0
STEP 2: Setup x-momentum equation set to solve for
.
Cell 1:
Cell 1 eqn is :
8
Cell 2:
Knowns: = =1. 0 , = =1.0 , =0.42 , =0.34 , = =0.0 =0.42 , =0.34 , =0.0 , =0.0 ̂ = 0.00.0.34421.0 =1.1905 , = 0.0.3384 =1.0952 0.42∗ 0.34∗ =0.380.00.0/2 Cell 2 eqn is :
Cell 3:
Knowns: = =1. 0 , = =1.0 , =0.34 , =0.26 , = =0.0 =0.34 , =0.26 , =0.0 , =0.0 ̂ = 0.00.0.23641.0 =1.3077 , = 0.0.3206 =1.1538 0.34∗ 0.26∗ =0.300.00.0/2 Cell 3 eqn is :
Cell 4:
Knowns: = =1. 0 , = =1.0 , =0.26 , =0.18 , = =0.0 =0.26 , =0.18 , =0.0 , =0.0 ̂ = 0.00.0.12861.0 =1.4444 , = 0.0.2128 =1.2222 0.26∗ 0.18∗ =0.220.00.0/2 Cell 4 eqn is :
Cell 5:
Knowns: = =1.0 , ==1.0 , =0.18 , =0.1 , = = =0.0 =0.18 , =0.1 , =0.0 , =0.0 ̂ = 0.00.0.1181.0 =1.8000 , = 0.0.114 =1.4000 0.18∗ 0.1∗ =0.1420.00.00.0/2 ∗ ∗∗ 1.1.14905706 0.0.4422 0.034 00 00 00 ∗∗ 0.05 00 0.034 0.0.2266 0.018 00 ∗∗ = 00 → ∗∗ = 1.2.97231778 0 0 0 0.18 0.1 ∗ 0 ∗ 0.5000 Cell 5 eqn is :
Discretized x-momentum equation system and the solution for
is
9
STEP 3: Calculate face velocities using Rhie-Chow interpolation and “star” velocities in Step 5.
Face 1: Face 2:
Face 3:
Face 4:
Face 5:
Face 6:
=0.6
. These velocities will be used as
=1.0 ̂ = + =1.2129 , = + =1.1064 → =0.61.21291.10640.00. 010.61.0=1.1277 ̂ = + =1.2715 , = + =1.1357 → =0.61.27151.13570.00. 010.61.0=1.1629 ̂ = + =1.3761 , = + =1.1880 → =0.61.37611.18800.00. 010.61.0=1.2256 ̂ = + =1.6222 , = + =1.3111 → 1. =0.61.62221.31110.00.010. 6 0 =1. 3 733 ̂ =̂ − =1.9778 , = − =1.4889 → =0.61.97781.48890.00.010. 61.0=1.5867
STEP 4: Calculate the coefficients of the pressure correction equation system solve for
Cell 1:
′
values.
=0.0 = =11.10640.42=0.4647 = =0.∗4647 ∗ =0.421.12770.51.0=0.0264 = =11.10640.42=0.4647 = =11. 1 3570. 3 4=0. 3 862 = ∗ =0.∗ 8509 =0.341.16290.421.1277=0.0783 = =11.13570.34=0.3862 = =11. 1 8800. 2 6=0. 3 089 = ∗ =0.∗ 6950 =0.261.22560.341.1629=0.0767 = =11.18800.26=0.3089 = =11. 3 1110. 1 8=0. 2 360 = ∗ =0.∗ 5449 =0.181.37330.261.2256=0.0715 RHS value:
Cell 2:
RHS value:
Cell 3:
RHS value:
Cell 4:
RHS value:
10
= =11.31110.18=0.2360 =0.0 = 2∗ =0.∗ 5338 =0.11.58670.181.3733=0.0885 ′ 0.0.40647467 0.0.84509647 0.30862 00 00 ′′ 0.0.00264783 00 0.30862 0.0.63950089 0.0.53449089 0.20360 ′′ = 0.0.00767715 → [ 0 0 0 0.2360 0.5338 ] ′ 0.0885
Cell 5:
RHS value:
Discretized PC equation system and the solution for
is
′′ 3.3.10321754 ′′ = 2.2.82045175 ′ 1.1463
STEP 5: Correct face velocities using Eqn (8).
Face 1: Face 2: Face 3: Face 4: Face 5: Face 6:
=1.∗0 ′ ′ = =1.12771.10643.07543.1321=1.1905 =∗ ′ ′=1.16291.13572.80453.0754=1.4706 =∗ ′ ′=1.22561.18802.21752.8045=1.9231 =∗ ′ ′=1.37331.31111.14632.2175=2.7778 =∗ (.−) =1.58671.4889 .−.. =5.0000 (No correction for the inlet velocity)
5 plus/minus typos. Results are correct.
Important note: As seen, corrected face velocities satisfy the continuity equation exactly, i.e. mass is conserved exactly in each cell with the corrected face velocities. For example consider cell 3.
Mass balance for cell 3: =10.261.923110.341.4706=0.0000
This is true for all cells. For this 1D problem with specified inlet velocity, these corrected face velocities are the same as the analytical values. For a 2D or 3D problem mass balance in each cell will again be exact, but the face velocities cannot be equal to the analytical values, which are probably not known anyway. This “exact mass balance” at each iteration is an important power of the SIMPLE algorithm. Although in this problem face velocities are exact, cell center velocities and pressures will require a number of iterations to converge and the converged values will not be equal to the analytical values. Accuracy will depend on the used mesh and the selected convergence tolerance.
11
STEP 6: Correct cell center velocities using Eqn (9).
Cell 1: Cell 2: Cell 3: Cell 4: Cell 5:
=∗∗ ′′ ′′⁄1 =1.19051.09523.07543.1321=1.2526 =∗ ′ ′ ⁄2 =1.47061.11762.80453.0754=1.6537 =∗ ′ ′ ⁄2 =1.92311.15382.21752.8045=2.4181 =∗ ′ ⁄′ 2 =2.77781.22221.14632.2175=3.7911 = ⁄1 =5.00001.22221.14632.2175=8.2096 =0. 4 =∗∗ ′′ =0.00.6∗3.1321=1.2529 =∗ ′ =0.00.6∗3.0754=1.2302 =∗ ′ =0.00.6∗2.8045=1.1218 =∗ ′ =0.00.6∗2.2175=0.8870 = =0.00.6∗1.1463=0.4585
STEP 7: Correct pressures using Eqn (10). Use a pressure relaxation value of
Cell 1: Cell 2: Cell 3: Cell 4: Cell 5:
.
STEP 8: Check for convergence. There are many possibilities here. One simple way is to compare velocity and pressure differences of two consecutive iterations. Convergence is declared if the maximum of such differences is less than a certain user specified tolerance value. It is preferred to use not just differences but normalize them in a logical way, e.g. normalize the velocity differences with respect to the given inlet velocity. It is always a good idea to select couple of critical monitoring points in the flow field and watch how variables change at those points before declaring convergence.
If the convergence check fails go to S tep 2 and perform one more iteration. Face and cell center velocities and cell center pressures calculated in this iteration will be used in the next iteration. The solution may also divergence, depending on how far the initial guess is from the exact solution, mesh density and relaxation factors.
ITERATION 2:
You can use the MATLAB code NS_1D_Colocated.m available at the co urse web site to perform a full solution. You can try different meshes, initial guesses, boundary condition implementations, relaxation factors, etc. Results of the second iteration are given below for you to check the correctness of your hand calculations.
STEP 2:
=0.9200 =0.7600 =0.6000 =0.4400 =0.2800
̂ =1.0000 ̂ =1.2526 ̂ =1.6537 ̂ =2.4181 ̂ =3.7911
∗∗ =1.0209 ∗ =1.0707 ∗ =1.1736 ∗ =1.3195 =1.5079 12
STEP 3:
=1.0000 =0.8400 =0.6800 =0.5200 =0.3600 =0.2000 STEP 4:
′′ =0.5818 =0.6141 ′′ =0.5645 ′ =0.2934 =0.5084
̂ =0.8737 ̂ =1.1263 ̂ =1.5043 ̂ =2.0640 ̂ =3.0664 ̂ =4.7967 STEP 5:
=1.0000 =1.1905 =1.4706 =1.9231 =2.7778 =5.0000
∗ =1.0000 ∗ =1.1634 ∗ =1.5043 ∗ =2.0640 ∗ =3.0664 ∗ =4.7957 STEP 6:
=1.0505 =1.0641 =1.0774 =1.0835 =1.7926
STEP 7:
=1.0201 =0.9845 =0.8960 =0.7996 =0.6619
Converged Solution with 5 cells
13
Converged Solution with 100 cells
14
