Computational Fluid Mechanics: Third Computer Project¶

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}\]

• Artificial Compressiblity Method (ACM)

  In the artificial compressibility method (ACM), the continuity equation is modified adding an unsteady term with ariticial compressiblity \(\beta\). To have this new form of continuity equation, an artificial equation of state that relates pressure, \(P\), to artificial density \(\tilde{\rho}\) is emploeyd as following form:

  \[P = \frac{\tilde{\rho}}{\beta}\]

  Finally, the modified continuity equation can then be recast as:

  \[\frac{\partial P}{\partial t} + \frac{1}{\beta} \frac{\partial u_{i}}{\partial x_{i}} = 0\]

• Vector form of transport equations

  Rewriting the previously drived non-dimensionalized continuity and momentum equation 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} = \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}P\\ u\\ v \end{bmatrix}, \;\; \vec{E} = \begin{bmatrix} \frac{u}{\beta}\\ uu + P\\ uv\end{bmatrix}, \;\; \vec{F} = \begin{bmatrix} \frac{v}{\beta}\\ uv\\ vv + P\end{bmatrix}\end{split}\]

  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

  In order to find time step that may stabilize the numerical solution, we need to know system convecting velocity as we pick the coefficient of spatial derivative terms in Burger’s and Euler equations as the convection velocity. The driven system of equation is not a single partial different equation but a set of three different partial different equation. To find the convection speed of numerical information in the time and space domains, we need to first linearize the given system of equations and find the Eigen values. The linearization can be obatained by following process. The driven system of PDE should be reformulated in linearized set of equations:

  \[\frac{\partial \vec{U}}{\partial t} + \left [ A \right ] \frac{\partial \vec{U}}{\partial x} + \left [ B \right ] \frac{\partial \vec{U}}{\partial y} = \frac{1}{\text{Re}} \left ( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \right ) \vec{U}\]

  Now we have found two coefficient matrices of convection terms and the spatial derivatives is now taken with respect to \(\vec{U}\) only. Despite the vector form, the PDE form is a identical with Burger’s equation. The coefficient matrices are below listed:

  \[\begin{split}\left [ A \right ] = \begin{bmatrix} 0 & \frac{1}{\beta} & 0 \\ 1 & 2u & 0\\ 0 & v & u \end{bmatrix}, \;\; \left [ B \right ] = \begin{bmatrix} 0 & 0 & \frac{1}{\beta} \\ 0 & v & u\\ 1 & 0 & 2v \end{bmatrix}\end{split}\]

  The resolved Eigen values of \(\left [ A \right ]\) and \(\left [ B \right ]\) matrices are \(u, u+a, u-a\) and \(v, v+a, v-a\), respectively. Taking \(\left [ A \right ]\) for example, the maximum convection velocity that transmit the numerical information can then be \(\left | u \right | + a\). Therefore, the Courant number for this case can also be determined by:

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.

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.

Problem1 - e, f¶

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

Re = 100

• Residual of u-velocity change in numerical iteration

  

• RMS error of u-velocity (reference: Ghia’s u-velocity data)

  

• Computational time with different grid spacing

  

Re = 500

• Residual of u-velocity change in numerical iteration

  

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?