7. Nonlinear Stabilization Operator (NSO)

An alternative to classic Peclet number blending is the usage of a discontinuity capturing operator (DCO), or in the low Mach context a nonlinear stabilization operator (NSO). In this method, an artifical viscosity is defined that is a function of the local residual and scaled computational gradients. Viable usages for the NSO can be advection/diffusion problems in addition to the aforementioned RTE VMS approach.

The formal finite element kernel for a NSO approach is as follows,

(7.1)$$\sum _e{\int }_{\mathrm{\Omega }}\nu (\mathbf{R})\frac{\partial w}{\partial x_i}{g}^{ij}\frac{\partial \varphi }{\partial x_j}d\mathrm{\Omega },$$

where $\nu (\mathbf{R})$ is the artifical viscosity which is a function of the pde fine-scale residual and $g^{ij}$ is the covariant metric tensor).

For completeness, the covariant and contravarient metric tensor are given by,

(7.2)$$g^{ij}=\frac{\partial x_i}{\partial {\xi }_k}\frac{\partial x_j}{\partial {\xi }_k},$$

and

(7.3)$$g_{ij}=\frac{\partial {\xi }_k}{\partial x_i}\frac{\partial {\xi }_k}{\partial x_j},$$

where $\xi =({\xi }_1,{\xi }_2,{\xi }_3{)}^T$. The form of $\nu (\mathbf{R})$ currently used is as follows,

(7.4)$$\nu =\sqrt{\frac{{\mathbf{R}}_{\mathbf{k}}{\mathbf{R}}_{\mathbf{k}}}{\frac{\partial \varphi }{\partial x_i}{g}^{ij}\frac{\partial \varphi }{\partial x_j}}}.$$

The classic paper by Shakib ( [SHZ91]) represents the genesis of this method which was done in the accoustically compressible context.

A residual for a classic advection/diffusion/source pde is simply the fine scale residual computed at the gauss point,

(7.5)$$\hat{\mathbf{R}}=\frac{\partial \rho \varphi }{\partial t}+\frac{\partial }{\partial x_j}(\rho u_j\varphi -{\mu }^{eff}\frac{\partial \varphi }{\partial x_j})-S$$

Note that the above equation requires a second derivative whose source is the diffusion term. The first derivative is generally determined by using projected nodal gradients. As will be noted in the pressure stabilization section, the advection term carries the pressure stabilization terms. However, in the above equation, these terms are not explicity noted. Therefore, an option is to subtract the fine scale continuity equation to obtain the final general form of the source term,

(7.6)$$\mathbf{R}=\hat{\mathbf{R}}-\varphi (\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_j}{\partial x_j}).$$

An alternative to the fine-scale PDE is a form that is found by differencing the linearized form of the residual from the nonlinear residual,

(7.7)$$\mathbf{R}=\frac{\partial \rho u_j\varphi }{\partial x_j}-(\varphi \frac{\partial \rho u_j}{\partial x_j}+\rho u_j\frac{\partial \varphi }{\partial x_j}).$$

The above resembles a commutation error in the nonlinear advection term.

In general, the NSO-$\nu$ is prone to percision issues when the gradients are very close to zero. As such, the value of $\nu$ is limited to a first-order like value. This parameter is proposed as follows: if an operator were written as a Galerkin (un-stabilized) plus a diffusion operator, what is the value of the diffusion coefficient such that first-order upwind is obtained for the combined operator? This upwind limited value of $\nu$ provides the highest value that this coefficient can (or should) be. The current form of the limited upwind $\nu$ is as follows,

(7.8)$${\nu }^{upw}=C_{upw}(\rho u_i{g}_{ij}\rho u_j{)}^{\frac{1}{2}}$$

where $C_{upw}$ is generally taked to be  0.1.

Using a piecewise-constant test function suitable for CVFEM and EBVC schemes (the reader is refered to the VMS RTE section), Eq. (7.1) can be written as,

(7.9)$$-\sum _e{\int }_{\mathrm{\Gamma }}\nu (\mathbf{R})g^{ij}\frac{\partial \varphi }{\partial x_j}{n}_{i}dS.$$

A fourth order form, which writes the stabilization as the difference between the Gauss-point gradient and the projected nodal gradient interpolated to the Gauss-point, is also supported,

(7.10)$$-\sum _e{\int }_{\mathrm{\Gamma }}\nu (\mathbf{R})g^{ij}(\frac{\partial \varphi }{\partial x_j}-G_j\varphi )n_{i}dS.$$

7.1. NSO Based on Kinetic Energy Residual

An alternative formulation explored is to share the general kernal form shown in Equation (7.10), however, compute $\nu$ based on a fine-scale kinetic energy residual. In this formulation, the fine-scale kinetic energy residual is obtained from the fine-scale momentum residual dotted with velocity. As with the continuity stabilization approach, the fine-scale momentum residual is provided by Equation (7.11). Therefore, the fine-scale kinetic energy is written as,

(7.11)$${\mathbf{R}}_{ke}=\frac{u_j(\frac{\partial p}{\partial x_j}-G_{j}p)}{2},$$

while the denominator for $\nu$ now includes the gradient in ke,

(7.12)$$\nu =\sqrt{\frac{{\mathbf{R}}_{ke}{\mathbf{R}}_{ke}}{\frac{\partial ke}{\partial x_i}{g}^{ij}\frac{\partial ke}{\partial x_j}}}.$$

The kinetic energy is simply given by,

(7.13)$$ke=\frac{u_k{u}_k}{2}$$

The kinetic energy form of $\nu$ is used for all equation sets with transformation by usage of a turbulent Schmidt/Prandtl number.

7.2. Local or Projected NSO Diffusive Flux Coefficient

While the NSO kernel is certainly evaluated at the subcontrol surfaces, the evaluation of $\nu$ can be computed by a multitude of approaches. For example, the artificial diffusive flux coefficient can be computed locally (with local residuals and local metric tensors) or can be projected to the nodes (via an $L_{oo}$ or $L_2$ projection) and re-interpolated to the gauss points. The former results in a sharper local value while the later results in a more filtered-like value. The code base only supports a local NSO $\nu$ calculation.

7.3. General Findings

In general, the NSO approach seems to work best when running the fourth-order option as the second-order implementation still looks more diffuse. When compared to the standard MUSCL-limited scheme, the NSO is the preferred choice. More work is underway to improve stabilization methods. Note that a limited set of NSOs are activated in the code base with specific interest on scalar transport, e.g, momentum, mixture fraction and static enthalpy transport. When using the $4^{th}$ order method, the consistent mass matrix approach for the projected nodal gradients is supported for higher order.

7.4. NSO as a Turbulence Model

The kinetic energy residual form has been suggested to be used as a turbulence model (Guermond and Larios, 2015). However, inspection of the above NSO kernel form suggests that the model form is not symmetric. Rather, in the context of turbulence modeling, is closer to the metric tensor acting on the difference between the rate of strain and antisymmetric tensor. As such, the theory developed, e.g., for eigenvalue perturbations of the stress tensor (see Jofre and Domino, 2017) can not be applied. In this section, a new form of the NSO is provided in an effort to be used for an LES closure.

In this proposed NSO formulation, the subgrid stress tensor, ${\tau }_{ij}^{sgs}=\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j$, is given by,

(7.14)$${\tau }_{ij}^{sgs}=-2\rho \nu g^{ij}(S_{ij}-\frac{1}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij})=-2\rho \nu g^{ij}{S}_{ij}^{\ast }.$$

Interestingly, the units of $\nu$ are of an inverse time scale while the product $2\rho \nu g^{ij}$ can be viewed as an non-isotropic eddy viscosity, ${\mu }_{ij}^t$.

The first order clipping may be relaxed by defining $\nu$ as,

(7.15)$$\nu =\frac{|{\mathbf{R}}_{ke}|}{||ke|{|}_{\mathrm{\infty }}}.$$

The above form would be closer to what Guermond uses and would avoid the divide-by-zero noted in regions of uniform flow.