Computational Fluid Mechanics: Fourth Computer Project¶

This document is online available. The pdf format of this documentation may have some inappropriate image size and not-well organized order of image and its description. You are recommended to see this documentation visiting http://cfm04.readthedocs.org.

Contents:

Code Instruction¶

The present project is aimed to develop a computer program for solving a steady solution with projection correction method. The code being used for answering all the question here is written with Python language. The current version of Python code is an extended version from the previous homework. Therefore, the current Python code has also a feature of pressure correction method as well as artificial compressibility method. This program is to run with simple command:

$ python main.py

Quick instruction for running the simulation¶

The Python code used for this project can be cloned from github.com repository:

$ git clone https://github.com/sayop/CFM04

You can also see the code directly by visiting the website: https://github.com/sayop/CFM04 If you clone the code, you will see the following set of files and directories:

$ sayop@reynolds:~$ ls CFM04/
docs  README.md  src

docs contains the document files set for the current project using Sphinx software. This pdf document is online available at: http://cfm04-gatech.readthedocs.org. The Python script for this simulation is stored in src folder.

Before running the simulation, you need to open the file named input.in using editor for example, VI on unix system:

$ vi input.in

Then, you should be able to see the following set of simulation parameters:

#grid dimension
iDim            20
jDim            20
xmin            0
xmax            10
ymin            0
ymax            10
#boundary conditions
uLeft           0.0
vLeft           0.0
uRight          0.0
vRight          0.0
uBottom         0.0
vBottom         0.0
uUp             10.0
vUp             0.0
# fluid kinematic properties
nu              1.0
pInit           10.0
#simulation setup
pCorr           1
alpha           1.0
pResidual       0.01
maxIter         100000000000
Courant         0.5
dtInit          0.0001
Beta            0.5
residualMin     0.00005
#Post-Process
nIterWrite      100

The parameter’s name above will literally tell you what every single variables indicates in the simulation. For the post-processing as requested in this project, nIterWrite will write a solution plot and CSV file at speicifed interval of time integration number.

The most important feature here is to set the pressure correction method. In order to have this goal, you will need to set pCorr to 1 to switch on its feature. Otherwise, you will run your simulation with artificial compressibility method.

Also, to set the specified Reynolds number, you need to change the uUp that will show a correponding Reynolds number at the beginning of your simulation on screen. Here, 10.0 of uUp will maintain the Reynolds number 100.0. All the input parameters are dimensional quantities. And these variables will be non-dimensionalized when they are transitioned to the main loop of the simulation.

Results¶

Problem1 - a¶

Describe the essential steps of the solution method. Include the discretized equations and implementation of boundary conditions.

In this problem set, we are supposed to solve the Navier-Stokes equations having continuity and momentum conservation equations together. Tensor forms of continuity and momentum equations are given below:

Continuity (incompressible)

\[\frac{\partial u_{i}}{\partial x_{i}} = 0\]
Momentum equation:

\[\frac{\partial u_{i}}{\partial t} + \frac{\partial u_{i}u_{j}}{\partial x_{j}} = -\frac{1}{\rho}\frac{\partial p}{\partial x_{i}} + \nu \frac{\partial}{\partial x_{j}}\left ( \frac{\partial u_{i}}{\partial x_{j}} \right )\]
Non-dimensionalization of the Navier-Stokes equations

In some cases, it is beneficial to non-dimensionalize the given transport equation because it eases the analysis of problem of interest, and also may reduce the number of parameters. The non-dimensionalized form of the Navier-Stokes equation can be achieved by first normalizing the primitive variables as followings:

\[\tilde{u_{i}} = \frac{u_{i}}{U_{\text{ref}}},\;\; \tilde{x_{i}} = \frac{x_{i}}{L_{\text{ref}}},\;\; \tilde{\rho}=\frac{\rho}{\rho_{\text{ref}}},\;\;\tilde{P} = \frac{P}{\rho_{\text{ref}}\, U^{2}_{\text{ref}}},\;\; \tilde{t}=\frac{t}{L/U_{\text{ref}}}\]

For the final form of non-dimensionalized Navier-Stokes equation, tilda, \(\tilde{}\), will be dropped out for brevity and a new non-dimensional physical parameter \(Re\) that represents the flow intertia against the fluid viscosity is introduced. Now we got:

\[\frac{\partial u_{i}}{\partial t} + \frac{\partial u_{i}u_{j}}{\partial x_{j}} = -\frac{1}{\rho}\frac{\partial p}{\partial x_{i}} + \frac{1}{\text{Re}} \frac{\partial}{\partial x_{j}}\left ( \frac{\partial u_{i}}{\partial x_{j}} \right )\]

where the Reynolds number is defined as:

\[\text{Re} = \frac{U_{\text{ref}}L_{\text{ref}}}{\nu}\]
Projection method

Given partial differention equation set is composed of continuity and momentum conservation equations. By contrary to the artificial compressiblity method, the projection method is to solve this equation set with incompressible solution method, so-called, pressure-based method. In this type of method, the continuity equation is not directly solved but used for divergence-free constraint, which is an outcome of incompressible flow continuity feature.

The momentum equation given in the tensor form can be recast in terms of time differential form as shown below. And the convective and diffusion flux terms are combined together to express a simple form of spatially differenciated flow quantities in terms of \(u\) and \(v\). And the pressure derivative term is left alone for a particular purpose of pressure correction.

\[\frac{\partial u_{i}}{\partial t} = H_{i} -\frac{\partial p}{\partial x_{i}}\]

where the convection terms are represented by \(H_{i}\):

\[H_{i} = -\frac{\partial u_{i}u_{j}}{\partial x_{j}} + \frac{1}{\text{Re}}\frac{\partial^{2} u_{i}}{\partial x_{j} \partial x_{j}}\]

Here, the above expression can be again discretized in two separate time steps having different right hand side quantities. Two separate the temporal difference, we employs an intermediate time indicated by \(*\). The intermediate time level of velocity is updated without pressure derivative term. So this time level is incomplete so it is not divergence free.

\[\frac{u^{*}_{i} - u^{n}_{i}}{\Delta t} = H_{i}\]

\[\frac{u^{n+1}_{i} - u^{*}_{i}}{\Delta t} = -\frac{\partial p^{n+1}}{\partial x_{i}}\]

Summing above two equations gives a complete form of temporal difference equation of momentum conservation equation. Then to find the divergence free and next time level of velocity components, we need to first evaluate the intermediate velocity quantities and pressure at the next time level. But, the given form of equation is not well posed because we haven’t yet solved next time level of pressure which also satisfies the divergence free constraint. To do this, we need to take divergence of the equation then it gives a elliptic equation form what is called Poisson’s equation.

\[\frac{1}{\Delta t} \frac{\partial u^{*}}{\partial x_{i}} = \frac{\partial^{2} p^{n+1}}{\partial x_{i} \partial x_{i}}\]

Since the equation is characterized with elliptic equation so it becomes boundary condition problem in the given domain. Because of two dimensional space and non-linearity, it needs to be solved with point-iterative method such as Jacobi method and Gauss-Seidel method. In this project, this solution for pressure is resolved with succesive over-relaxation method to have faster convergence.
Vector form of transport equations

Rewriting the previously drived non-dimensionalized momentum equations in vector form generates a simple format that eases implementation of the numerical method. The above transport equation can be newly formed as shown below:

\[\frac{\partial \vec{U}}{\partial t} + \frac{\partial \vec{E}}{\partial x} + \frac{\partial \vec{F}}{\partial y} + \vec{\triangledown }p = \frac{1}{\text{Re}} \left ( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \right ) \vec{U}\]

where the each of vector elements are summarized below:

\[\begin{split}\vec{U} = \begin{bmatrix}u\\ v \end{bmatrix}, \;\; \vec{E} = \begin{bmatrix} uu \\ uv\end{bmatrix}, \;\; \vec{F} = \begin{bmatrix} uv\\ vv \end{bmatrix}, \;\; \vec{\triangledown} = \begin{bmatrix} \frac{\partial}{\partial x}\\ \frac{\partial}{\partial y} \end{bmatrix}\end{split}\]

Compared to the previous homework that constructs the flux vector including pressure force, the current flux vectors \(E\) and \(F\) contain only convection terms in x and y directions. This is the purpose that the numerical method separates two different steps: evaluating intermediates velocities and projection of this velocities onto the divergence free satisfied velocity field.

Now this is good to go further for descritization because the given task is to solve explicit form of discretization equation. Even though the derived form of transport equation is not linearized, each of vectors above are easily discretized in terms of their elements that are combinations of each primitive variables. Thus, in this project, actual discretization has been doen form the driven transport equation above.
Finding time step algorithm

Contrary to the previous homework, the system of partial differential equation is composed of only momentum conservation equations in x and y direction, that is, this system does not contain continuity equation. Note that the mass conservation only applies as a divergence free constraint. So the system convection velocity is straighforward rather than having to find eigenvalues of coefficient matrices. So the Courant number relation for the two-dimensional system of equation can be defined by:

\[\text{Courant} = \frac{u dt}{dx} + \frac{v dt}{dy}\]

The numerical time step is evaluted as a minimum time step that satifies the given relation at every computational node points.

Problem1 - b, c¶

Consider the case when \(H=W\) (a square cavity). Here, the Reynolds number, \(Re=UW/\nu\), characterizes the flow patters. Compute the steady state solutions for both \(Re=100\) and \(Re=500\). Plot the flow streamlines and centerline profiles (\(u\) vs. \(y\) and \(v\) vs. \(x\) through the center of the domain). For \(Re=100\), valdiate your method by comparing your results to data from given literature.

Re = 100¶

In this test, the lid cavity’s velocity is set to make the Reynolds number set to 100. To see the qualitative effect of different grid spacing, four different grid resolution conditions is employed and compared together in this page.

NxN = 10x10

Observation

Very coarse grid resolution doesn’t produce well-predicted data when it is compared to the reference data.

Even though the streamlines seems to penetrate the wall, it does not necessarily mean it pass through it. It it because the default streamline generation feature of Python does not produce properly when it is resolved on less grid points.

NxN = 20x20

Observation

Denser grid resolution tends to produce better results. The resolved u and v velocities look closer to the reference data.

Compared to 10x10 case, the streamline produced with denser grid resolution looks more physically reasonable.

NxN = 40x40

NxN = 60x60

Observation

Having resolution of 60x60 makes finally the resolved data looks very close to the reference data.

We observed the denser grid size produces the more well-matching data with reference data.

Re = 500¶

In this test, the lid cavity velocity is set to 50 m/s to make the Reynolds number 500. Two different grid spacing are employed to see the qualitative pattern of grid size effect on numerical solution.

NxN = 20x20

NxN = 40x40

Observation

The faster lid cavity velocity makes the distictive vortex at two corners on the bottom where as Re=100 case produces very weak vortext at the same location.

The denser grid resolution makes the vortext look more distinctive and bigger than the coarser grid resolution.

Problem1 - d¶

Examine the method stability with different grids. Determine the maximum time step that leads to a stable solution and compare it to the stability criteria.

Grid spacing test¶

In this project, variable time stepping method is employed such that temporal numerical instability is avoided. Thus, given conditions of Re=100 and Re=500 are numerically stable with relevant choice of Courant number. For this reason, the grid spacing test is performed with much higher Reynolds number condition such that physical viscousity effect is not significantly taken into account. Note that we already know Euler equation having no viscousity is unconditionally unstable with central finite difference method.

On the other hand, having visous terms will attenuate the possible growth of numerical instability and leads to stable solution with properly chosen Courant number. All the preliminary tests with given Re conditions gave stable solutions. This is because viscous terms play a role of smoothing out the stiff solutions so as to result in numerically stable solution.

Here, the test has conducted with Re = 10,000 such that the equation almost goes to Euler equation. We can expect that numerical solution will be unconditionally unstable because our employed finite difference is done with central differencing. The test shown below proves this. On the other hand, the effect of grid spacing on the unstable solution’s temporal evolution is clarified. Even though all the cases go to unstable solution after all, the onset of amplified numerical instability appears at different location in iteration number.

Effect of grid spacing on the temporal evolution of axial velocity residual. (Re = 10,000)

Maximum time step¶

In this code, the variable time step method is used to maintain stable numerically. Therefore, the code does not run with constant time step. The maximum time step test is performed with different set of Courant number condition. The grid spacing is fixed with 20x20 to have fast running of simulation.

As described earlier, the current project is to set the numerical time step variable according to the stability criteria taking Courant number into consideration. The Courant number is essentially based on convective velocity that transport the flow quantities that may propagate the numerical signal waves. Since the current system of equation is based on incompressible flow equation, the convective velocity is simply equivalent to the local velocity components at every node points. As we discovered the numerical instability characteristics with Burger’s equation, the required Courant number is to be less than 1.0. In this regard, the time step experiment in this section seems to well satisfy this requirement. Having Courant number less than 1.0 shows well stabilized solution as we observed from the above plot.

Since the current code is featured with variable time step, single value of time step can be representative of the required time step. Instead, we employed a certain numbers of Courant numbers as listed below. The representative time step value is chosen when the temporal iteration reaches 100.

Courant # dt at 100th iteration

0.5 0.038791

0.8 0.061655

1.0 0.064686

1.1 0.063479

1.2 0.072822

Problem1 - e, f¶

Compare your results with different grid resolutions to evaluate the numerical error and the order of the scheme.

Re = 100

Here, required iteration number for temporal integration is illustrated for the different grid spacing condition. User-defined convergence rate for this test was set to 0.00005 as a normalized residual number.

Observation for iteration number
- As the grid spacing becomes smaller, the required iteration numbers to satisfy the pre-specified convergence rate gets bigger. So 60x60 case takes almost 4,000 iterations to have steady-state solution.
- But it does NOT necessarily mean that it reaches the steady state at later physical time. According to the Courant number definition, the time step for smaller grid spacing must be smaller than the coarser grid resolution.

In this test, the root-mean squared error is evaluted for different grid spacing conditions. The RMS error was resolved by comparing the numerically resolved solution to the Ghia’s data when it reachs the steady-state solution. The steady-state is assumed to be reached when the normalized resdiual number becomes less than the pre-specified convergence rate as mentioned above.

Observation for RMS error
- As expected, the smaller grid spacing, the smaller error is obtained.

Here, the computational times consumed for steady-state solution are compared with different grid resolution. We observe dramatically increased computational time for dense grid spacing. This is because the more grid points consumes more time for resolving Poisson’s equation solution. Even use of successive over-relaxation method for point-iterative method for pressure correction dramatically increases the number of iteration. The massive computational load was made with the pressure correction step at every time step.

Problem1 - g¶

Repeat parts b and c for the case of \(H=1.5W\) (except for validation). How does the flow change in this rectangular cavity as compared to the flow in a sqaure cavity?

Re = 100¶

In this section, two choices of grid spacing conditions were made to see the effect of grid size on the qualitative pattern of flow field. Only the streamline is considered for discussions.

NxN = 20x30

NxN = 40x60

Observation
- Very distinctive flow feature is observed from two different grid spacing choice.
- Higher number of grid points resolves more distinctive vortext at two bottom corners.
- The center of main vortex close to lid is identically positioned with different grid sets.

Re = 500¶

The current test is conducted with higher Reynolds number condition.

NxN = 20x30

NxN = 40x60

Observations
- The higher Reynolds number condition makes two distinctive and opposite direction of vortex vector.