Introduction
For strong viscous-inviscid interactions and large separated flow regimes (e.g., capturing shock-boundary layer interactions, wake flows, leading edge flows, etc.), the full Navier-Stokes equations must be solved. Two strategies are employed depending on whether compressibility effects exist. For compressible flows, the system is solved in a coupled manner with \rho as the primary variable. For incompressible flows, an iterative approach is used, called the segregated approach, with p as the primary variable.
Solving the compressible Navier-Stokes equations (both viscous and inviscid) is very similar to that of the Euler equations, as the unsteady compressible form of the Navier-Stokes equations are hyperbolic-parabolic, whereas the unsteady compressible Euler equations are hyperbolic. Therefore, both the viscous and inviscid form of the Navier-Stokes equations can be solved similarly with a marching procedure. The same is not true for incompressible: treating \rho as a constant causes the reduced Euler equations to become elliptic. When viscosity is included, the incompressible Navier-Stokes equations become elliptic-parabolic, which requires a marching problem. Therefore, the inviscid and viscous incompressible equations must be solved with different strategies.
A key difference between solving viscous and inviscid flows is the need to refine the mesh near the surface. Since viscous flows must satisfy no-slip at the walls, there is a significant velocity gradient near the solid boundaries that require small cell heights. Further away in the freestream, the cells may become larger as the velocity gradients are not as strong. Large aspect ratio cells must be used in the viscous regions of the flow to capture the fluid physics, which oftentimes decreases the rate of convergence.
Compressible Navier-Stokes Equations
Incompressible Navier-Stokes Equations
When heat transfer is negligible and the fluid’s properties remain relatively constant, the incompressible form of the Navier-Stokes equations may be used. The conservations are as follows.
\begin{gathered}
\nabla\cdot\textbf V=0 \\
\rho\frac {\mathrm D\textbf V}{\mathrm Dt}=-\nabla p+\mu\nabla^2\textbf V \\
\rho c_v\frac {\mathrm DT}{\mathrm Dt}=k\nabla^2T+\Phi
\end{gathered}
\tag{1}
The unknowns are velocity vector \textbf V, pressure p, and temperature T. Since fluid properties are constant and independent of T, the energy equation can be decoupled from continuity and momentum.
Incompressible, sometimes called constant density, flows arise when changes in fluid properties due to temperature are negligible. For these cases, the properties (i.e., \mu, \rho, c_p, etc.) can be considered independent of temperature, so the mass and momentum equations are independent of temperature. This permits either equation to be solved without knowing the thermal state of the fluid.
The 2D constant density conservation equations in primitive form (neglecting energy equation) are as follows.
\begin{gathered}
\frac {\partial u}{\partial x}+\frac {\partial v}{\partial y}=0 \\
\frac {\partial u}{\partial t}+u\frac {\partial u}{\partial x}+v\frac {\partial u}{\partial y}=-\frac 1{\rho}\frac {\partial p}{\partial x}+\nu\left(\frac {\partial^2u}{\partial x^2}+\frac {\partial^2u}{\partial y^2}\right) \\
\frac {\partial v}{\partial t}+v\frac {\partial u}{\partial x}+v\frac {\partial v}{\partial y}=-\frac 1{\rho}\frac {\partial p}{\partial y}+\nu\left(\frac {\partial^2v}{\partial x^2}+\frac {\partial^2v}{\partial y^2}\right)
\end{gathered}
\tag{2}
These equations may be solved using derived or primitive variables. Derived quantities, such as vorticity and stream functions, are most useful for 2D case as they introduce an extra layer of derivatives and the equations become very complicated in three dimensions.
Derived Variable Approach
Primitive Variable Approach
The derived variable approach (using \zeta and \psi), is rarely used when solving 3D incompressible flow fields. The primitive approach is more widely used for both 2D and 3D flows.
General
The nondimensional primitive form of Eqn. 1, ignoring energy, is
\begin{gathered}
\frac {\partial u^*}{\partial x^*}+\frac {\partial v^*}{\partial y^*}+\frac {\partial w^*}{\partial z^*}=0 \\
\frac {\partial u^*}{\partial t^*}+u^*\frac {\partial u^*}{\partial x^*}+v^*\frac {\partial u^*}{\partial y^*}+w^*\frac {\partial u^*}{\partial z^*}=-\frac {\partial p^*}{\partial x^*}+\frac 1{\mathrm{Re}_L}\left(\frac {\partial^2u^*}{\partial x^{*2}}+\frac {\partial^2u^*}{\partial y^{*2}}+\frac {\partial^2u^*}{\partial z^{*2}}\right) \\
\frac {\partial v^*}{\partial t^*}+u^*\frac {\partial v^*}{\partial x^*}+v^*\frac {\partial v^*}{\partial y^*}+w^*\frac {\partial v^*}{\partial z^*}=-\frac {\partial p^*}{\partial y^*}+\frac 1{\mathrm{Re}_L}\left(\frac {\partial^2v^*}{\partial x^{*2}}+\frac {\partial^2v^*}{\partial y^{*2}}+\frac {\partial^2v^*}{\partial z^{*2}}\right) \\
\frac {\partial w^*}{\partial t^*}+u^*\frac {\partial w^*}{\partial x^*}+v^*\frac {\partial w^*}{\partial y^*}+w^*\frac {\partial w^*}{\partial z^*}=-\frac {\partial p^*}{\partial z^*}+\frac 1{\mathrm{Re}_L}\left(\frac {\partial^2w^*}{\partial x^{*2}}+\frac {\partial^2w^*}{\partial y^{*2}}+\frac {\partial^2w^*}{\partial z^{*2}}\right)
\end{gathered}
\tag{3}
The nondimensional parameters are
\begin{array}{ccc}
u^*=\dfrac u{V_{\infty}} & x^*=\dfrac xL & p^*=\dfrac p{\rho_{\infty}V_{\infty}^2} \\ \\
v^*=\dfrac v{V_{\infty}} & y^*=\dfrac yL & t^*=\dfrac {V_{\infty}t}L \\ \\
w^*=\dfrac w{V_{\infty}} & z^*=\dfrac zL & \mathrm{Re}_L=\dfrac {V_{\infty}L}{\nu}
\end{array}
\tag{4}
The primitive variable approach can be divided into two categories. The first is the coupled approach, where all the unknowns are solved simultaneously. For the unsteady Navier-Stokes equations, this is implemented through the artificial compressibility method. An artificial time derivative of pressure term is added to the continuity equation so the solution to the coupled hyperbolic system may be advanced with time. Without it, the algebraic system is singular. When the solution converges, the artificial pressure-time term disappears.
The second approach is the pressure-correction approach, or the segregated method. Here, a separate equation is used to determine p and the momentum equations are solved for \textbf V.
Coupled Approach: The Method of Artificial Compressibility
Pressure-Correction Approach: SIMPLE
The semi-implicit method for pressure linked equations (SIMPLE) approach is based on an iterative process where the velocity is first solved from the momentum equations based off a guessed pressure field, and then corrected to satisfy continuity. This process repeats until all conservation equations are satisfied. The actual property A is written as the sum of an intermediate value, A_0, and corrected term, A'. For pressure and velocity, they become
\begin{aligned}
p & =p_0+p' \\
u & =u_0+u' \\
v & =v_0+v'
\end{aligned}
\tag{5}
Through an iterative approach, starting with a guess for the pressure field, the final p and velocity components, u and v, may be solved using the procedure below.
Start with an initial guess for the pressure, p_0, and solve the momentum equations to determine u_0 and v_0.
Solve the pressure correctione quation below to determine p' in each cell. Note that \textbf V_0 is the estimated velocity vector and consists of u_0 and v_0 velocity components.
\nabla^2 p'=\frac 1A(\nabla\textbf V_0)=\frac 1A\left(\frac {\partial u_0}{\partial x}+\frac {\partial v_0}{\partial y}\right)
\tag{6}
- Correct the p_0, u_0, and v_0 values obtained from Step 1 using the formulas below.
\begin{aligned}
p & =p_0+p' \\
u & =u_0-\frac A{2\Delta x}\left(p_{i+1,j}'-p_{i-1,j}'\right) \\
v & =v_0-\frac A{2\Delta y}\left(p_{i,j+1}'-p_{i,j-1}'\right)
\end{aligned}
\tag{7}
- Replace the p_0, u_0, and v_0 values calculated from Step 1 with the p, u, and v calculated from Eqn. 7 and return to Step 2. Repeat this process until the solution converges.
The pressure correction is related to the velocity correction through a linearized form of the momentum equation.
\begin{aligned}
\rho\frac {\partial u'}{\partial t} & =-\frac {\partial p'}{\partial x} \\
\rho\frac {\partial v'}{\partial t} & =-\frac {\partial p'}{\partial y}
\end{aligned}
\tag{8}
Since for any given iteration, u_k=u_{0,k}+u_k', then the previous iteration is u_{k-1}=u_{0,k-1} and the correction term for the previous step is zero. Therefore, Eqn. 8 can be written as
\begin{aligned}
u' & =-A\frac {\partial p'}{\partial x} \\
v' & =-A\frac {\partial p'}{\partial y}
\end{aligned}
\tag{9}
where A=\frac {\Delta t}{\rho}. The \Delta t in the A equation is a fictitious time for steady flows relevant only to the iterative sequence. Substituting Eqn. 9 and Eqn. 5 into mass conservation gives
\cancel{\frac {\partial u}{\partial x}+\frac {\partial v}{\partial y}}-\left(\frac {\partial u_0}{\partial x}+\frac {\partial v_0}{\partial y}\right)+\frac {\Delta t}{\rho}\left(\frac {\partial^2p'}{\partial x^2}+\frac {\partial^2p'}{\partial y^2}\right)=0
The first term disappears since \nabla\cdot\textbf V=0 and rearranging the remaining terms gives Eqn. 6.
In practice, Eqn. 6 tends to overestimate the needed p' to fully correct u_0 and v_0. To accelerate the rate of convergence, an underrelaxation constant, \alpha_p, is added to the pressure correction equation.
p=p_0+\alpha_p p'
Two notable extensions of the SIMPLE algorithm are the revised (SIMPLER) scheme (Patankar, 1981) and the SIMPLEC scheme (van Doormal and Raithby, 1984). The SIMPLEC scheme uses a more complete approximation of the momentum equation to compute p' to improve convergence. Recall that with SIMPLE, the convective and viscous (diffusive) terms of the velocity correction were neglected. Physically, SIMPLE assumes that all changes to the velocity are due to the error in the guessed pressure field. By incorporating some of these convective and diffusive terms back, the scheme more completely captures the changes in the flow field. This translates to a modification to the A coefficient in Eqn. 9. With this modification, van Doormal and Raithby reported that it was no longer needed to underrelax p'. SIMPLEC generally performs more efficiently than both SIMPLE and SIMPLER and is usually one of the default pressure correction schemes in most CFD softwares.