3.4.3. Smoothing algorithms defined

Now that the smoothing and splitting errors have been formally defined, it is useful to consider three projection algorithms that have been implemented and verified within SIERRA/Fuego in the context of the classic two equation $k$ - $\epsilon$ model, with steady method of manufactured solutions (MMS) (cf. [228]).

3.4.3.1. Fourth-order smoothing with characteristic or time step scaling

In this algorithm, the projection time scales are defined by either

(3.760) $\tau = \tau_1 = \tau_2 = \tau_3 = \tau_{char},$

or

(3.761) $\tau = \tau_1 = \tau_2 = \tau_3 = {\bf I}\Delta t.$

Here, characteristic scaling, $\tau_{char}$ , is a diagonal matrix that represents a time scale based on convection and diffusion contributions, while for time step scaling, the time scale is based on the local time step. The characteristic scaling very closely follows the standard finite element method stabilization parameter.

The smoothing and splitting errors are now given by

(3.762) $\begin{bmatrix} {\bf A} & {\bf G} \\ {\bf D} & {\bf 0} \end{bmatrix} \begin{bmatrix} {\bf u}^{n+1} \\ p^{n+1/2} \end{bmatrix} = \begin{bmatrix} {\bf f} \\ {\bf b} \end{bmatrix} + \begin{bmatrix} ({\bf I}- {\bf A }\tilde \tau ){\bf G}(p^{n+{1\over2}}-p^{n-{1\over2}}) \\ ( {\bf L_{\tau}}-{\bf D}\tilde \tau{\bf G })p^{n+{1\over2}} \end{bmatrix} .$

Of particular interest to this research is the role of the stabilization term, $({\bf L_{\tau}}-{\bf D}\tilde \tau{\bf G })p^{n+{1\over2}}$ , on formal time accuracy when ${\bf \tau} = {\bf I}\Delta t$ (a scheme that has been shown to display more appealing stability characteristics). Clearly, a scheme that uses explicit pressure stabilization with time step scaling is first-order accurate. Expanding this stabilization term shows the fourth-order pressure derivative scaled by a length scale cubed. Therefore, by refining the time step $and$ mesh, one might be able to show a second-order accuracy for sufficiently resolved meshes.

In practice, the stabilization terms are carried within the mass flow rate that forms part of the right-hand side of the Pressure Poisson Equation solve and the convection term for the transport of any scalar field. The mass flow rate is defined as

(3.763) $\dot m^k = \left(\overline{\rho}\overline{ \hat {\bf u}} + \overline{\tau {{\bf G}p^{n-{1\over2}}} \over {V}} -\overline \tau \nabla^h p^{n+{1\over2}}\right)d{\bf A},$

where the introduction of the over bar is noted to represent interpolation of a nodal field to an integration point. Note that in the bulk of the collocated unstructured finite volume literature, the form of the mass flow rate defines the stabilization (the difference between the nodal gradient operator ${\bf G}$ and the interior element operator $\nabla^h$ ). Above we note the independent interpolation of the density and velocity rather than $\overline{\rho \hat{\bf u}}$ , as is done in Stanford’s ASC Alliance code CDP. It does seem that the full interpolation of $\rho\hat{\bf u}$ may be more consistent, although the effect of this algorithmic detail has not been explored.

3.4.3.2. Stabilized smoothing

The stabilized projection algorithm is based on the work of [218], that was derived from the monolithic scheme of [217]. In this algorithm, the projection time scales are defined as

(3.764) $\tau_1 = \Delta t {\bf I} + \tau_{char}.$

(3.765) $\tau_2 = \tau_{char} .$

(3.766) $\tau_3 = \tau_{char}.$

With the above definitions, the smoothing and splitting errors are now defined as

(3.767) $\begin{bmatrix} {\bf A} & {\bf G} \\ {\bf D} & {\bf 0} \end{bmatrix} \begin{bmatrix} {\bf u}^{n+1} \\ p^{n+{1\over2}} \end{bmatrix} = \begin{bmatrix} {\bf f} \\ {\bf b} \end{bmatrix} + \begin{bmatrix} ({\bf I}- {\bf A }\tilde \tau_{char} ){\bf G}(p^{n+{1\over2}}-p^{n-{1\over2}}) \\ ( {\bf L_{\tau_{char}}}-{\bf D}\tilde \tau_{char}{\bf G })p^{n+{1\over2}} + \Delta t {\bf L}(p^{n+{1\over2}}-p^{n-{1\over2}}) \end{bmatrix} .$

The mass flow rate now includes an additional stabilization factor and is now defined as

(3.768) $\dot m^k = \left(\overline{\rho}\overline{ \hat {\bf u}} + \overline{\tau {{\bf G}p^{n-{1\over2}}} \over {V}} -\overline \tau \nabla^h p^{n+{1\over2}}-\Delta t {\bf L}\Delta p^{n+{1\over2}}\right)d{\bf A} .$

Note that at full convergence, the stabilized scheme reduces to the fourth-order characteristic scaling algorithm.

3.4.3.3. Second-order smoothing with characteristic or time step scaling

In fact, the scaled nodal gradient need not be included in the mass flow rate equation, e.g.,

(3.769) $\dot m^k = \left(\overline{\rho}\overline{ \hat {\bf u}} -\overline \tau \nabla^h p^{n+{1\over2}}\right)d{\bf A}.$

This is equivalent to neglecting the $\tilde \tau_2{\bf G}p^{n-{1 \over2}}$ term in (3.741), or by defining $\tilde {\bf u} = \hat {\bf u}$ .

The smoothing for this algorithm is provided by the local Laplacian operator. The smoothing and splitting errors for this method are now given by

(3.770) $\begin{bmatrix} {\bf A} & {\bf G} \\ {\bf D} & {\bf 0} \end{bmatrix} \begin{bmatrix} {\bf u}^{n+1} \\ p^{n+1/2} \end{bmatrix} = \begin{bmatrix} {\bf f} \\ {\bf b} \end{bmatrix} + \begin{bmatrix} ({\bf I}- {\bf A }\tilde \tau ){\bf G}(p^{n+{1\over2}}-p^{n-{1\over2}}) \\ ({\bf L_{\tau}} - {\bf D}\tilde \tau{\bf G })\Delta p^{n+{1\over2}} + {\bf L_{\tau}} p^{n-{1\over2}} \end{bmatrix} .$

3.4.3.4. Zeroth-order smoothing with time step or characteristic scaling

Certainly, the pressure smoothing can be removed, i.e., $\tau_2 = 0$ , that leads to the following set of errors,

(3.771) $\begin{bmatrix} {\bf A} & {\bf G} \\ {\bf D} & {\bf 0} \end{bmatrix} \begin{bmatrix} {\bf u}^{n+1} \\ p^{n+{1\over2}} \end{bmatrix} = \begin{bmatrix} {\bf f} \\ {\bf b} \end{bmatrix} + \begin{bmatrix} ({\bf I}- {\bf A } \tilde \tau ){\bf G}(p^{n+{1\over2}}-p^{n-{1\over2}}) \\ ( {\bf L_{\tau}}-{\bf D}\tilde \tau{\bf G })(p^{n+{1\over2}} -p^{n-{1\over2}}) \end{bmatrix} .$

where $\tau$ is either the characteristic scale, $\tau_{char}$ , or the simulation time step, ${\bf I}\Delta t$ (with $\tau_1 = \tau_3$ ). Although the converged error is zero, this lack of smoothing can lead to a decoupled pressure field in certain flows.

Here, the mass flow rate reduces to a simple interpolation of nodal velocities within the element

(3.772) $\dot m^k = \left(\overline{\rho}\overline{ \hat {\bf u}} -\overline{\tau} \nabla^h \Delta p^{n+{1\over2}} \right)d{\bf A}.$

The unsmoothed algorithm is very similar to the staggered formulation of SIMPLE, (cf. [130]), with $\tau = A_p^{-1}$ (the inverse of the diagonal matrix from operator {bf A}). However, by design, the staggered mesh arrangement holds the property that $( {\bf L_{\tau}}-{\bf D}\tau{\bf G } ) = 0$ . In this method, no stabilization is added as none is required.

3.4.3.5. Time integration scheme

At present, three time integration schemes are supported in the code base that include: 1) first-order backward Euler, 2) second-order BDF2 and, 3), the Crank-Nicholson method described in [229].

The general two- and three-state scheme time derivatives are simply written as:

(3.773) $\int \frac{\partial \rho \phi }{\partial t} dV = \int \frac{ (\gamma_1 \rho^{n+1} \phi^{n+1} + \gamma_2 \rho^{n} \phi^{n} + \gamma_3 \rho^{n-1} \phi^{n-1})} {\Delta t} dV$

where $\gamma_i$ represent the appropriate factors for either Backward Euler or a three-point BDF2 scheme. The above time derivative is either nodally lumped or evaluated at the subcontrol volume quadrature points. For a two-state backward Euler scheme, $\gamma_1$ is unity while $\gamma_2$ is negative unity. For a given variable time step, the BDF2 factors are,

(3.774) $\tau = \Delta t^N/\Delta t^{N-1},$

(3.775) $\gamma_1 = \frac{(1+2\tau)}{(1+\tau)},$

(3.776) $\gamma_2 = (1+\tau),$

(3.777) $\gamma_3 = \frac{\tau^2}{(1+\tau)}.$

For a fixed time step, the $\gamma$ -factors reduce to the canonical $(\frac{3}{2}, -2, \frac{1}{2})$ set.

In the Crank-Nicholson implementation, the generalized method is written as

(3.778) ${{\partial \phi} \over {\partial t}}^{n+1} = \eta {{\phi^{n+1}-\phi^n} \over {\Delta t}} + (1-\eta) {{\partial \phi} \over { \partial t}}^{n},$

where $\eta$ is a blending coefficient between 1 and 2. Values of $\eta$ of unity result in first-order backward Euler, while values of 2 result in second order Crank-Nicholson, i.e.,

(3.779) ${{\partial \phi} \over {\partial t}}^{n+1} = 2{{(\phi^{n+1} - \phi^{n})} \over {\Delta t} } - {{\partial \phi} \over {\partial t}}^n.$

A linearization is given by

(3.780) ${{\partial \phi} \over {\partial t}}^{n+1} = 2{{( \phi^{k} - \phi^{n})} \over {\Delta t} } - {{\partial \phi} \over {\partial t}}^n,$

where the old time derivative is computed based on the old solution of the partial differential equation of interest. The above algorithm is especially useful in that it avoids the need to evaluate complex right-hand side source terms at the n+1 and n state, e.g., simulations that include the need to compute turbulence production, reaction rate terms, etc.

3.4.3.6. Variable density

In the case of variable density, the same set of options exists in the code. Specifically, the time derivative is written as (for backward Euler or BDF2),

(3.781) $\int \frac{\partial \rho }{\partial t} dV = \int \frac{ (\gamma_1 \rho^{n+1} + \gamma_2 \rho^{n} + \gamma_3 \rho^{n-1} )} {\Delta t} dV$

For the Crank/Nicolson scheme, the full time term is

(3.782) ${{\partial \rho \phi} \over {\partial t}}^{n+1} = \eta {{ \rho^{n+1} \phi^{n+1}- \rho^{n} \phi^n} \over {\Delta t}} + (1-\eta) {{\partial \rho \phi} \over {\partial t}}^{n},$

where it is noted that the full time derivative at $n^{th}$ state is saved. The linearization is given by

(3.783) ${{\partial \rho \phi} \over {\partial t}}^{n+1} = \eta {{ \rho^{k} \phi^{k}- \rho^{n} \phi^n} \over {\Delta t}} + (1-\eta) {{\partial \rho \phi} \over {\partial t}}^{n}.$

In practice, the usage of the above formula for a second-order density derivative has proven unstable. As such, $\eta$ is set to unity.

3.4.3.7. Shifted density iteration

A modified iteration originating from Shunn, Ham and Moin:cite:shunn:2012 is available. The iteration procedure proceeds in words as

Copy forward density state as $\rho_0 = \rho^n$ as a provisional density for the step
Solve all scalar transport equations using $\rho_0$
Solve momentum, with $\rho_0$
Compute $\rho_1$ from the equation of state using the result of the scalar transport equations
Compute the density-time derivative using $\rho_1$
Solve the pressure Poisson equation using the newly computed density-time derivative and the provisional density and velocity
Update the mass flux using the newly computed pressure consistent with the pressure Poisson equation
Project momentum to the new $\rho_1 u_1$ state consistent with the pressure Poisson equation.

For a simplified case of a single scalar $Z$ determining density through an equation of state $\rho_R$ , predict a solution as $Z_0 = Z^n, \, \rho_0 =\rho^n, \, u_0 = u^n, \, p_0 = p^n$ , then iterate

(3.784) $\int \frac{\gamma_1}{\Delta t} \rho_0 Z_1 \, {\rm d}V + \int \left( \dot{m}_{0} Z_1 - {\rm Sc}^{-1} \mu \nabla_h Z_1 \cdot A \right) \,{\rm d} S = -\int \frac{1}{\Delta t}\left( \gamma_2 \rho^n Z^n + \gamma_3 \rho^{n-1} Z^{n-1} \right) \, {\rm d}V$

(3.785) $\int \frac{\gamma_1}{\Delta t} \rho_0 u^\star \, {\rm d}V + \int \left( \dot{m}_{0} u^\star - 2\mu \nabla_{s, h} u^\star \cdot A \right) \,{\rm d} S = -\int \left[ \frac{1}{\Delta t}\left( \gamma_2 \rho^n u^n + \gamma_3 \rho^{n-1} u^{n-1} \right) + G p_0 \right] \, {\rm d}V$

(3.786) $\rho_1 = \rho_R(Z_1), \, \left(\partial_t \rho \right)_1 = \frac{1}{\Delta t} \left( \gamma_1 \rho_1 + \gamma_2 \rho^n + \gamma_3 \rho^{n-1} \right)$

(3.787) $\int \tau \nabla_h p_1 \, {\rm d}S = -\int \left(\partial_t \rho \right)_1 \, {\rm d}V + \int \left( \rho_0 u^\star + \tau G p_0 \right) \cdot A \, {\rm d} S$

(3.788) $\dot{m}_{1}|_{\rm ip} = \left[\left(\rho_0 u^\star + \tau G p_0\right)|_{\rm ip} - \tau|_{\rm ip} \nabla_h p_1|_{\rm ip} \right] \cdot A$

(3.789) $\int G p_1 \, {\rm d}V = \int p_1 A \, {\rm d S}$

(3.790) $u_1 = \frac{1}{\rho_1}\left[\rho_0 u^\star - \tau \left(G p_1 - G p_0\right) \right],$

where $\nabla_s$ is the symmetric gradient. The procedure helps the iteration procedure to converge more robustly with large density time-derivatives. Numerical evidence:cite:shunn:2012 indicates that recovering second-order accuracy with the iteration process depends strongly on the nature of the state equation, with an MMS at Courant number unity requiring 5 outer-iterations for a density ratio of 5 but 20 with a density ratio of 20.

3.4.3.8. Volume-averaged properties

Beta Capability

Volume-averaged properties is a beta capability, and will substantially increase the cost of property evaluation.

Different options for computing volume-averaged properties are available. For some property $\psi$ depending on a state parameterization $\phi$ , we have

(3.791) $\overline{\psi_i} = \frac{1}{\operatorname{vol}(\Omega_i)}\int_{\Omega_i} \psi \left(\phi_h \right) \rm{d V},$

where $\phi_h$ is the state parameterization represented in the nodal finite element basis. Using a nodal quadrature for the integral is the default and is equivalent to direct nodal evaluation of the properties. Alternatively, the quadrature of the consistent mass matrix can be used. A more accurate quadrature for the properties is useful to reduce aliasing errors arising due to the nonlinear nature of the state relationship. Additionally properties maybe smoothed directly with $n$ -passes of a nodal filter,

(3.792) $\hat{\psi}^n = \left(M^{-1}_L M_c\right)^n \psi.$

where $M_L$ and $M_c$ are the lumped and consistent mass matrices respectively, described in Sec. Discrete system of equations. This is equivalent to a box-filter on the dual mesh in CVFEM and helps reduce spurious oscillations in the property evaluation, but also introduces an additional ( $\mathcal{O} (h^2)$ ) error with each pass.