# 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 - model, with steady method of manufactured solutions (MMS) (cf. [229]).

## 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)

or

(3.761)

Here, characteristic scaling, , 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)

Of particular interest to this research is the role of the stabilization term, , on formal time accuracy when (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 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)

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 and the interior element operator ). Above we note the independent interpolation of the density and velocity rather than , as is done in Stanford’s ASC Alliance code CDP. It does seem that the full interpolation of 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 [219], that was derived from the monolithic scheme of [218]. In this algorithm, the projection time scales are defined as

(3.764)

(3.765)

(3.766)

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

(3.767)

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

(3.768)

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)

This is equivalent to neglecting the term in (3.741), or by defining .

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)

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

Certainly, the pressure smoothing can be removed, i.e., , that leads to the following set of errors,

(3.771)

where is either the characteristic scale, , or the simulation time step, (with ). 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)

The unsmoothed algorithm is very similar to the staggered formulation of SIMPLE, (cf. [131]), with (the inverse of the diagonal matrix from operator {bf A}). However, by design, the staggered mesh arrangement holds the property that . 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 [230].

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

(3.773)

where 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, is unity while is negative unity. For a given variable time step, the BDF2 factors are,

(3.774)

(3.775)

(3.776)

(3.777)

For a fixed time step, the -factors reduce to the canonical set.

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

(3.778)

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

(3.779)

A linearization is given by

(3.780)

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)

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

(3.782)

where it is noted that the full time derivative at state is saved. The linearization is given by

(3.783)

In practice, the usage of the above formula for a second-order density derivative has proven unstable. As such, 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 as a provisional density for the step

  - Solve all scalar transport equations using

  - Solve momentum, with

  - Compute from the equation of state using the result of the scalar transport equations

  - Compute the density-time derivative using

  - 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 state consistent with the pressure Poisson equation.

For a simplified case of a single scalar determining density through an equation of state , predict a solution as , then iterate

(3.784)

(3.785)

(3.786)

(3.787)

(3.788)

(3.789)

(3.790)

where 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 depending on a state parameterization , we have

(3.791)

where 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 -passes of a nodal filter,

(3.792)

where and 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 () error with each pass.