Coupled Flow Solver

Governing Equations in Vector Form

The system of governing equations in Cartesian integral form for an arbitrary control volume $V$ with differential surface area $d\mathbf{\text{a}}$ can be written as:

The vectors are defined as:

where:

• $\rho$ is the density of the fluid.
• $\mathbf{\text{v}}$ is the velocity of the fluid.
• $p$ is the pressure.
• $E$ is the total energy per unit mass.
• $\mathbf{\text{T}}$ is the viscous stress tensor.
• $\dot{q}$ is the heat flux vector.
• $\mathbf{\text{H}}$ is the vector of body forces.

Total energy is related to the total enthalpy H by:

where:

and $h=C_{p}T$

where $C_p$ is the specific heat capacity and $T$ is the temperature.

For low Mach number and generally incompressible flows, the system of equations can become numerically stiff, which in turn can lead to convergence problems. To remedy this problem, preconditioning is introduced.

Preconditioned Equations

To provide efficient solution of both compressible and incompressible flows at all speeds, a preconditioning matrix $\Gamma$ is incorporated into Eqn. (946) by multiplication with the transient term [294]:

where:

and ${\rho }_T=\frac{\partial \rho }{\partial T}{|}_p$ is the derivative of density with respect to temperature at constant pressure.

For $\delta$ = 1, this matrix becomes a member of Turkel’s family of preconditioners [290].

The parameter $\theta$ is defined as:

The reference velocity $U_r$ appearing in Eqn. (953) is chosen such that the eigenvalues of the system remain well-conditioned with respect to the convective and diffusive time scales [294]. This is accomplished by limiting $U_r$ such that it does not go below either the local convection velocity or the diffusion velocity. An additional limitation on $U_r$ considers local pressure differences and is designed to increase numerical stability in stagnation regions by prohibiting amplification of pressure perturbations there. Thus, the restriction on $U_r$ becomes:

where $\Delta x$ is the inter-cell length scale over which the diffusion occurs and $\delta p$ is the pressure difference between adjacent cells. For compressible flows, the maximum reference velocity is limited to the local speed of sound, $c$ . In this case, the parameter $U_{r,\text{max}}$ in Eqn. (954) is replaced by $c$ , and the value prescribed for $U_{r,\text{max}}$ is not considered. The scaling parameter $\epsilon$ appearing in Eqn. (954) is set to 2.

The non-preconditioned Navier-Stokes equations are recovered exactly from Eqn. (946) by setting:

In this case, $\Gamma$ reduces to the Jacobian $\partial W/\partial Q$ .

Although Eqn. (950) is conservative for steady-state flows, it is not strictly conservative for time-dependent flows [290]. A dual time-stepping method is used to recover a time-accurate solution of the preconditioned equations for unsteady flows with implicit time stepping.

Applying Eqn. (950) to a cell-centered control volume for cell-0, the following discretized system is obtained:

Inviscid Fluxes

Based on upwind-differencing, there are two schemes available for discretizing the inviscid fluxes in Eqn. (955)

Weiss-Smith Preconditioned Roe’s Scheme

By default, the inviscid fluxes appearing in Eqn. (955) are evaluated by using the Weiss-Smith preconditioned Roe’s flux-difference splitting scheme as described in [293] and [294]. This scheme acknowledges that ${\mathbf{\text{f}}}_f$ contains characteristic information propagating through the domain with speed and direction according to the eigenvalues of the system. By splitting ${\mathbf{\text{f}}}_f$ into parts, where each part contains information traveling in a particular direction (that is, characteristic information), and upwind differencing the split fluxes in a manner consistent with their corresponding eigenvalues, the following expression for the value of the flux at each face is obtained:

$${\mathbf{\text{f}}}_f=\frac{1}{2}({\mathbf{\text{f}}}_0+{\mathbf{\text{f}}}_1)-\frac{1}{2}\Gamma \left|\mathbf{\text{A}}\right|\Delta \mathbf{\text{Q}}$$

(956)

where the subscripts “0” and “1” refer to the cells on either side of face-f, $\Gamma$ is preconditioning matrix, and $\Delta \mathbf{\text{Q}}$ is:

$$\mathrm{\Delta }\mathrm{}Q=(Q_1^r-Q_0^r)$$

(957)

where ${\mathbf{\text{Q}}}_0^r$ and ${\mathbf{\text{Q}}}_1^r$ are the solution vectors from cell-0 and cell-1 interpolated to the face using reconstruction gradients.

The matrix $\left|\mathbf{\text{A}}\right|$ is defined as:

$$\left|\mathbf{\text{A}}\right|=\mathbf{\text{M }}\left|\Lambda \right|{\mathbf{\text{ M}}}^{-1}$$

(958)

where $\Lambda$ is the diagonal matrix of eigenvalues and $\mathbf{\text{M}}$ is the modal matrix that diagonalizes ${\Gamma }^{-1}(\partial \mathbf{\text{f}}/\partial \mathbf{\text{Q}})$ .

By using these reconstructed solution vectors, the discretization scheme becomes formally second-order accurate. In this form, Eqn. (956) can be viewed as a second-order central-difference plus an added matrix dissipation, combining the pressure and velocity differential dissipations. This added matrix dissipation term is not only responsible for producing an upwinding of the convected variables, and of pressure and flux velocity in supersonic flow. It also provides the pressure-velocity coupling that is required for stability and efficient convergence of low-speed and incompressible flows.

The mass-flux pressure dissipation term is modified to provide improved numerical dissipation for the entire range of Strouhal numbers, from Strouhal number zero (steady state) to high Strouhal numbers.

This modification affects the pressure differences term of the mass-flux dissipation. It is inspired by the method proposed by Hosangadi et al. [248], but with the following differences:

• It improves the accuracy of simulation for unsteady flows that require both sound and advected quantities (like total enthalpy and entropy).
• It allows you to specify a limit to the allowed amount of dissipation correction, for more robust simulation. The Unsteady Preconditioning Max Factor sets this limit. This factor has a range from 0 (no correction) to 1 (full correction).
• It bases the mass-flux dissipation corrections on local mesh and flow conditions.
• It does not require you to provide a reference length scale.

Liou’s AUSM+ Scheme

Inviscid fluxes can be evaluated using Liou’s AUSM+ flux-vector splitting scheme, as described in [230], [259] and other references. This scheme is also based on the upwind concept, which is applied consistently on the two physically distinct parts of the inviscid flux: convective and pressure. Thus, the numerical flux is written as:

$${\mathbf{\text{f}}}_{\mathbf{\text{f}}}={\mathbf{\text{f}}}^{\text{(c)}}+\mathbf{\text{P}}={\mathbf{\text{m}}}_{\mathbf{\text{i}}}^{\text{+}}{(1,u,v,\dots ,H)}_0^T+{\mathbf{\text{m}}}_{\mathbf{\text{i}}}^{\text{-}}{(1,u,v,\dots ,H)}_1^T+{\mathbf{\text{P}}}_i$$

(959)

Here ${\mathbf{\text{m}}}_{\mathbf{\text{i}}}$ is the mass flux, ${\mathbf{\text{m}}}_{\mathbf{\text{i}}}^{\text{+}}$ is defined as $0.5({\mathbf{\text{m}}}_{\mathbf{\text{i}}}+\left|{\mathbf{\text{m}}}_{\mathbf{\text{i}}}\right|)$ , ${\mathbf{\text{m}}}_{\mathbf{\text{i}}}^{\text{-}}$ is defined as $0.5({\mathbf{\text{m}}}_{\mathbf{\text{i}}}-\left|{\mathbf{\text{m}}}_{\mathbf{\text{i}}}\right|)$ , and ${\mathbf{\text{P}}}_i$ is the pressure flux. The mass flux and pressure flux are computed based on local flow characteristics, emphasizing on correct information propagation inside the fluid for convective and acoustic processes.

The mass-flux pressure dissipation term is modified to provide improved numerical dissipation for the entire range of Strouhal numbers, from Strouhal number zero (steady state) to high Strouhal numbers.

This modification follows the method proposed by Hosangadi et al. [248], but with the following differences:

• It bases the mass-flux dissipation corrections on local mesh and flow conditions.
• It does not require you to provide a reference length scale.

Steady Solution Algorithms

For steady-state simulations, the coupled system of equations Eqn. (955) is discretized in time and time-stepping is performed until a quasi steady-state solution obtained. Simcenter STAR-CCM+ provides an explicit and an implicit procedure to perform the time-integration.

Explicit Time-Stepping

An explicit multi-stage time-stepping scheme [255] is used to discretize the time-derivative in Eqn. (950). The solution is advanced from time $t$ to time $t+\Delta t$ with an m-stage Runge-Kutta scheme:

$$\begin{array}{c}{\mathbf{\text{Q}}}^{\left(0\right)}={\mathbf{\text{Q}}}_t\\ {\mathbf{\text{Q}}}^{\left(i\right)}={\mathbf{\text{Q}}}^{\left(0\right)}-{\alpha }_i\Delta t{\mathrm{\Gamma }}^{-1}{\mathbf{\text{R}}}^{\left(i-1\right)}\\ {\mathbf{\text{Q}}}_{t+\Delta t}={\mathbf{\text{Q}}}^{\left(m\right)}\end{array}$$

(960)

where $i=1,2,\dots ,m$ is the stage counter for the m-stage scheme and ${\alpha }_i$ is the multi-stage coefficient for the ith stage.

For explicit time-stepping, the value of the time-step $\mathrm{\Delta }t$ is not specified manually. The user inputs CFL (Courant number) and $\mathrm{\Delta }t$ is determined by Eqn. (968).

The residual ${\mathbf{\text{R}}}^{\left(i\right)}$ is computed from the intermediate solution ${\mathbf{\text{Q}}}^{\left(i\right)}$ and, for Eqn. (955), is:

$${\mathbf{\text{R}}}^{\left(i\right)}=\frac{1}{V}\begin{array}{c}\sum \end{array}_{\begin{array}{c}f\end{array}}\left\{\mathbf{\text{f}}\left({\mathbf{\text{Q}}}^{\left(i\right)}\right)-\mathbf{\text{g}}\left({\mathbf{\text{Q}}}^{\left(i\right)}\right)\right\}$$

(961)

Residual Smoothing

Residual smoothing is a mechanism for increasing the explicit time-step size by removing the high-wave number oscillations from the residuals. The residuals for cell $i$ are filtered through a Laplacian smoothing operator:

$${\stackrel{\_}{\mathbf{\text{r}}}}_i={\mathbf{\text{r}}}_i+\epsilon \sum _{\mathrm{n}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{b}\mathrm{o}\mathrm{r}\mathrm{s}}\left({\stackrel{\_}{\mathbf{\text{r}}}}_j-{\stackrel{\_}{\mathbf{\text{r}}}}_i\right)$$

(962)

where the subscript $j$ refers to the neighboring cells.

This equation is implemented using a Jacobi iteration as follows:

$${\stackrel{\_}{\mathbf{\text{r}}}}_i^m=\frac{{\text{r}}_i^0+\epsilon {\begin{array}{c}\\ \begin{array}{c}\begin{array}{c}\sum \end{array}\\ \text{neighbors}\end{array}\end{array}}_{}^{}{\stackrel{\_}{\mathbf{\text{r}}}}_j^{m-1}}{1+\epsilon \begin{array}{c}\\ \begin{array}{c}\begin{array}{c}\sum \end{array}\\ \text{neighbors}\end{array}\end{array}1}$$

(963)

where ${\mathbf{\text{r}}}_i^0$ represents the original (unsmoothed) residuals and ${\stackrel{\_}{\mathbf{\text{r}}}}_i^m$ represents the smoothed residuals after iteration $m$. The default number of smoothing operations is 2 with the under-relaxation factor $\epsilon$ = 0.5. This is typically sufficient to allow the Courant number to be doubled.

Implicit Time-Stepping

The coupled system of equations Eqn. (946) is discretized in time by use of an Euler implicit scheme in conjunction with a Newton-type linearization of the fluxes .

$$[{\mathbf{\text{D}}}_i+{\displaystyle \sum _j^{{\text{N}}_{\text{faces}}}}{\text{S}}_{j,k}]\Delta \mathbf{\text{Q}}=-{\mathbf{\text{R}}}_i^n$$

(964)

where $\Delta \mathbf{\text{Q}}\equiv {\mathbf{\text{Q}}}^{n+1}-{\mathbf{\text{Q}}}^n$. The center and off-diagonal coefficient matrices, ${\text{D}}_i$ and ${\text{S}}_{j,k}$, and the residual vector ${\mathbf{\text{r}}}_i^n$ are:

$${\text{D}}_i=\frac{V}{\Delta t}\Gamma +{\displaystyle \sum _j^{{\text{N}}_{\text{faces}}}}{\text{S}}_{j,k}$$

(965)

$${\text{S}}_{j,k}=\frac{\partial {\mathbf{\text{F}}}_j}{\partial {\mathbf{\text{Q}}}_k}-\frac{\partial {\mathbf{\text{G}}}_j}{\partial {\mathbf{\text{Q}}}_k}$$

(966)

$${\mathbf{\text{R}}}_i^n={\displaystyle \sum _j^{\text{Nfaces}}}{\left({\mathbf{\text{F}}}_j-{\mathbf{\text{G}}}_j\right)}^n$$

(967)

Transient Solution Algorithms

For transient simulations, Simcenter STAR-CCM+ provides explicit and implicit time-stepping schemes. The implicit time-stepping procedure is also known as dual-time formulation.

Explicit Time-Stepping

In Simcenter STAR-CCM+, explicit time-stepping for transient simulations corresponds to the explicit time-stepping scheme that is described in Explicit Time-Stepping for steady-state simulations.

The physical time-step $\Delta t$ for unsteady simulations is applied uniformly to all cells in the domain.

The time-step $\Delta t$ , appearing in Eqn. (960) for explicit, unsteady calculations, is computed by consideration of the Courant number and Von Neumann stability conditions:

$$\Delta t=\begin{array}{c}\underbrace{\text{min }}\\ \forall x\end{array}\left(\frac{\text{CFL }V\left(\mathbf{\text{x}}\right)}{{\lambda }_{\text{max}}\left(\mathbf{\text{x}}\right)},\frac{VNN\Delta x^2\left(\mathbf{\text{x}}\right)}{\nu \left(\mathbf{\text{x}}\right)}\right)$$

(968)

The chosen time-step is the minimum local time-step over all the fluid cells in the model, and is equal to the minimum value of the pseudo-time step $\Delta \tau$ (Eqn. (971)) in the domain.

Implicit Time-Stepping (Dual Time-Stepping)

Since time-derivative preconditioning destroys the time accuracy of the governing equations, solving for unsteady flows requires a dual time-stepping, with inner iterations in pseudo-time. A preconditioned pseudo-time derivative term is introduced into conservation Eqn. (946) as follows:

$$\frac{\partial }{\partial t}\underset{V}{\overset{}{{\int }_{}^{}}}\mathbf{\text{W}}dV+\Gamma \frac{\partial }{\partial \tau }\underset{V}{\overset{}{{\int }_{}^{}}}\mathbf{\text{Q}}dV+\oint \left[\mathbf{\text{F}}-\mathbf{\text{G}}\right]\cdot d\mathbf{\text{a}}=\underset{V}{\overset{}{{\int }_{}^{}}}\mathbf{\text{H}}dV$$

(969)

where $t$ denotes the physical time-step which is applied uniformly to all cells in the domain. $\tau$ is a locally varying pseudo-time used in the time-marching procedure. As inner iterations progress to convergence, the second term of Eqn. (969) vanishes and Eqn. (946) is recovered.

The pseudo-time derivative is driven to zero with the following iterative algorithm:

$$\begin{array}{c}{\mathbf{\text{Q}}}^{\left(0\right)}={\mathbf{\text{Q}}}_n\\ [\Gamma +\frac{3}{2}\frac{\Delta \tau }{\Delta t}\frac{\partial \mathbf{\text{W}}}{\partial \mathbf{\text{Q}}}]\Delta \mathbf{\text{Q}}=\Delta \tau \{{\mathbf{\text{R}}}^{\left(i-1\right)}+\frac{1}{2\Delta t}(3{\mathbf{\text{W}}}^{\left(i-1\right)}-4{\mathbf{\text{W}}}^{\left(n\right)}+{\mathbf{\text{W}}}^{\left(n-1\right)})\}\\ {\mathbf{\text{Q}}}_{n+1}={\mathbf{\text{Q}}}^{\left(m\right)}\end{array}$$

(970)

where $i$ is the inner iteration counter, $n$ represents any given physical time level, $\mathrm{\Delta }\mathbf{\text{Q}}\equiv {\mathbf{\text{Q}}}^{\left(i\right)}-{\mathbf{\text{Q}}}^{\left(0\right)}$ , and $\mathrm{\Delta }t$ is the time-step.

Throughout the pseudo-time iterations, ${\mathbf{\text{W}}}^{\left(n\right)}$ and ${\mathbf{\text{W}}}^{\left(n-1\right)}$ are held constant and ${\mathbf{\text{W}}}^{\left(i-1\right)}$ is computed from ${\mathbf{\text{Q}}}^{\left(i-1\right)}$ . When inner iterations converge, the solution at the next physical time level ${\mathbf{\text{W}}}^{(i+1)}$ is given by $\mathbf{\text{W}}\left({\mathbf{\text{Q}}}_{\tau +\Delta \tau }\right)$ . For each inner iteration, the linear system (block coupled for the explicit scheme or fully coupled for the implicit scheme) is solved to compute $\Delta \mathbf{\text{Q}}$ .

The local pseudo time-step $\Delta \tau$ associated with the iterations of all steady simulations, or within a time-step of a transient simulation, is computed considering the Courant number and Von Neumann stability conditions.

$$\Delta \tau =\text{min}\left(\frac{\text{CFL }V\left(\mathbf{\text{x}}\right)}{{\lambda }_{\text{max}}\left(\mathbf{\text{x}}\right)},\frac{VNN\Delta x^2\left(\mathbf{\text{x}}\right)}{\nu \left(\mathbf{\text{x}}\right)}\right)$$

(971)

where $\text{CFL}$ is the user-specified Courant number (dimensionless), $V$ is the cell volume ( ${\text{m}}^3$ ), $VNN$ is the Von Neumann number ( $VNN\approx 1$ ), $\Delta x$ is a characteristic cell length scale (m) and $\nu$ is kinematic viscosity ( ${\text{m}}^2\text{/s}$ ). ${\lambda }_{\text{max}}$ is the maximum eigenvalue of the system as follows:

$${\lambda }_{\text{max}}=\left|\mathbf{\text{u}}\mathbf{\cdot }\mathbf{\text{a}}\right|+c\left|\mathbf{\text{a}}\right|$$

(972)

where $\mathbf{\text{u}}$ is velocity and $c$ the speed of sound.

The maximum eigenvalue is determined from the pre-conditioned system, Eqn. (950). The field function Convective Courant Number provides the local convective Courant number for a given time-step.

For implicit unsteady computations, Eqn. (965), you specify the time-step.

Continuity Convergence Accelerator Algorithm

The Continuity Convergence Accelerator pressure-correction discretization uses the Coupled Flow numerical dissipation, not the Rhie-Chow dissipation of the Segregated solver. Choosing the Coupled Flow numerical dissipation formulation over the Rhie-Chow formulation prevents some of the known issues of the Segregated solver, especially in high-speed regimes, such as wiggles at flow discontinuities due to improper numerical dissipation.

When activated, the Continuity Convergence Accelerator algorithm adds a post-update conditioning step to help satisfy better the continuity equation.

The Coupled Implicit solver algorithm with Continuity Convergence Accelerator is:

1. As part of the normal Coupled Implicit iteration update, solve the discretized coupled flow equations, and perform the normal update to obtain an intermediate flow field state $p^{*}$, $v^{*}$...
2. Compute the new cell mass fluxes and mass imbalances for this intermediate state, $m_f^{*}$.
3. Discretize and solve the Continuity Convergence Accelerator pressure-correction equation with the new mass imbalances in the right-hand side:

   $$a_p{{\mathrm{p}}^{\prime }}_p+{\displaystyle \sum _n^{}}{a}_n{{\mathrm{p}}^{\prime }}_n=-{\displaystyle \sum _n^{}}{m}_f^{*}$$

   (973)

   to produce cell values of the pressure correction $\mathrm{p\prime }$. The subscript $n$ denotes all neighbors of cell $p$ with faces $f$.

4. Update cell pressures using the under-relaxation factor $\omega$:

   $$p^{n+1}=p^{*}+\omega p^{\prime }$$

   (974)
5. Update the boundary pressure corrections ${\mathrm{p\prime }}_b$, and then pressure $p_b$.
6. Correct the face mass fluxes:

   $${m_f}^{n+1}=m_f^{*}+{\text{m}}^{\prime }$$

   (975)

   where ${\text{m}}^{\prime }$ is computed using Eqn. (936).

   If the Enhanced Mass-Imbalance Calculations expert property is activated, the mass imbalances are recomputed using either Eqn. (956) or Eqn. (959).

7. Correct the cell velocities

   $$v^{n+1}=v^{*}-f\left(\omega \right)\frac{V\Delta {\text{p}}^{\prime }}{A_{p}v}$$

   (976)

   where:

   • $f\left(\omega \right)$ is a non-linear function of $\omega$
   • $\Delta {\text{p}}^{\prime }$ is the gradient of the pressure corrections
   • $A_p$ is the vector of central coefficients for the original coupled equations discretized system representing the velocity equations
   • $V$ is the cell volume
8. Update density, total enthalpy, and other variables with the pressure and velocity changes.