Momentum Advection

The momentum advection process represents the advective transport of $u_i$ (i.e. momentum normalized by density) in rotation form, i.e.

\[\frac{∂u_i}{∂t} = … − u_j \left(\frac{∂u_i}{∂x_j} − \frac{∂u_j}{∂x_i}\right) − \frac{∂}{∂x_i} \frac{u_j u_j}{2} \,.\]

Only the first term is computed here while the second term is assumed to become part of a modified pressure handled separately. See Consistent Discrete Energy Term for considerations on how to evaluate this term numerically e.g. when computing budgets.

The process is discretized by computing the vorticity in the frequency domain, then computing the velocity–vorticity product in the physical domain, and finally transforming the result back to a frequency-domain representation. If other processes (e.g. a StaticSmagorinskyModel) require physical-domain velocity gradients at the same resolution, the vorticity is computed from those to minimize the number of transforms.

If the level of dealiasing is set to at least :quadratic, the physical-domain products are equivalent to frequency-domain convolution and no additional aliasing errors are introduced.

The process relies on boundary conditions for $u₃$.

MomentumAdvection(; dealiasing = :quadratic)

Discretization

To compute the non-linear advection term, the velocity components $u_i$ as well as the vorticity components $ω_i$ are transformed to the physical domain.

The contributions of the advection term become

\[\begin{aligned} \frac{∂}{∂t} u₁(ζ_C) &= … + u₂(ζ_C) \, ω₃(ζ_C) - \frac{u₃(ζ_C⁻) \, ω₂(ζ_C⁻) + u₃(ζ_C⁺) \, ω₂(ζ_C⁺)}{2} \\ \frac{∂}{∂t} u₂(ζ_C) &= … + \frac{u₃(ζ_C⁻) \, ω₁(ζ_C⁻) + u₃(ζ_C⁺) \, ω₁(ζ_C⁺)}{2} - u₁(ζ_C) \, ω₃(ζ_C) \\ \frac{∂}{∂t} u₃(ζ_I) &= … + \frac{u₁(ζ_I⁻) + u₁(ζ_I⁺)}{2} \, ω₂(ζ_I) - \frac{u₂(ζ_I⁻) + u₂(ζ_I⁺)}{2} \, ω₁(ζ_I) \end{aligned}\]

where the values that are defined at the “wrong” set of nodes have been interpolated.

Contributions to Budget Equations

Contribution to the instantaneous momentum equation:

\[\begin{aligned} \frac{\partial}{\partial t} u_i &= … - u_j \frac{\partial u_i}{\partial x_j} \quad &\text{Convection Form} \\ &= … - \frac{\partial u_i u_j}{\partial x_j} \quad &\text{Divergence Form} \\ &= … - \frac{1}{2} \left( \frac{\partial u_i u_j}{\partial x_j} + u_j \frac{\partial u_i}{\partial x_j} \right) \quad &\text{Skew-Symmetric Form} \\ &= … - u_j \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) - \frac{\partial}{\partial x_i} \frac{u_j u_j}{2} \quad &\text{Rotation Form} \end{aligned}\]

Contribution to the mean momentum equation:

\[\begin{aligned} \frac{\partial}{\partial t} \overline{ u_i } = … - \overline{ u_j \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) } - \frac{\partial}{\partial x_i} \frac{\overline{u_j u_j}}{2} \end{aligned}\]

\[\begin{aligned} \frac{\partial}{\partial t} \overline{ u_i } = … + \frac{\partial}{\partial x_j} \left( − \overline{u_i}\,\overline{u_j} − \overline{u_i^\prime u_j^\prime} \right) \end{aligned}\]

\[\begin{aligned} \frac{\partial}{\partial t} \overline{ u_i } = … - \frac{\partial \overline{u_i}\,\overline{u_j}}{\partial x_j} - \frac{\partial \overline{u_i^\prime u_j^\prime}}{\partial x_j} \end{aligned}\]

Contribution to turbulent momentum equation:

\[\begin{aligned} \frac{\partial}{\partial t} u_i^\prime = … - \frac{\partial}{\partial x_j} \left( u_i^\prime \overline{u_j} + \overline{u_i} u_j^\prime + u_i^\prime u_j^\prime - \overline{u_i^\prime u_j^\prime} \right) \end{aligned}\]

Contribution to the mean kinetic energy equation:

\[\begin{aligned} \frac{\partial}{\partial t} \frac{\overline{u_i}^2}{2} &= … - \overline{u_i} \frac{\partial \overline{u_i}\,\overline{u_j}}{\partial x_j} - \overline{u_i} \frac{\partial \overline{u_i^\prime u_j^\prime}}{\partial x_j} \\ &= … \frac{\partial}{\partial x_j} \left( - \frac{\overline{u_i}^2}{2} \overline{u_j} - \overline{u_i} \overline{u_i^\prime u_j^\prime} \right) \underbrace{+ \overline{u_i^\prime u_j^\prime} \frac{\partial \overline{u_i}}{\partial x_j}}_{\text{TKE production}} \end{aligned}\]

Contribution to the turbulent kinetic energy equation:

\[\begin{aligned} \frac{\partial}{\partial t} \frac{\overline{u_i^\prime u_i^\prime}}{2} = … + \frac{\partial}{\partial x_j} \left( - \frac{\overline{u_i^\prime u_i^\prime}}{2} \overline{u_j} - \frac{1}{2} \overline{u_i^\prime u_i^\prime u_j^\prime} \right) \underbrace{- \overline{u_i^\prime u_j^\prime} \frac{\partial \overline{u_i}}{\partial x_j}}_{\text{TKE production}} \end{aligned}\]

Consistent Discrete Energy Term

Since the rotational form of the advection term groups a part of the advection term with the pressure gradient, we need to undo this “grouping” if we want to obtain the actual, physical advection and pressure terms from the quantities that are computed numerically. In particular, we need to decide how to evaluate $∂(u_j u_j)/∂x_i$ numerically, ideally in a manner that is consistent with the numerical methods of the advection and pressure term. It appears that this is indeed possible.

First, we can recognize that the vorticity–velocity product consists of the two parts $u_j\,∂u_i/∂x_j$ and $u_j\,∂u_j/∂x_i$ (ignoring the signs). The first part is exactly the convection form of the momentum advection term while the second part is what we are aiming to define numerically, as it is analytically equal to the gradient of $(u_j u_j)/2$.

Does this mean that the simulation already defines a numerical evaluation of $∂(u_j u_j)/∂x_i$? Not quite – there are two remaining issues: First, the “diagonal” terms with $i = j$ are never evaluated numerically, as they appear in both the “convection” and “energy” terms and therefore cancel. Second, the simulation also includes a definition of how the $∂/∂x_i$-gradient is evaluated in the pressure term, so we still have the potential for inconsistent numerics. The question therefore becomes the following: Can we define numerical expressions for the scalar $u_j^2$ on pressure-nodes as well as $u_i ∂u_i/∂x_i$ on the $u_i$-nodes such that $2 u_j\,∂u_j/∂xi = ∂u_j^2/∂x_i$ when evaluated with the numerical methods of the pressure term and the rotational advection term?

We can examine each of the nine $i,j$-pairs individually on the $u_i$-nodes, where both the advection and pressure terms are evaluated.

For the four combinations of the indices $1$ and $2$, all terms are defined on $\zeta_C$-nodes, so there are no interpolations involved and all the derivatives are spectral. If we define the numerical evaluation of $u_1^2$ and $u_2^2$ as the fully-dealiased product on the (“native”) $\zeta_C$-nodes, then

\[2 u_1\,∂u_1/∂x_2 = ∂u_1^2/∂x_2 \quad \text{and} \quad 2 u_2\,∂u_2/∂x_1 = ∂u_2^2/∂x_1\]

hold for the off-diagonal terms since all computations are exact up to the cut-off wavenumbers. Note that this is only true with full dealiasing, since otherwise the aliased wavenumbers are handled inconsistently. This “natural” definition of $u_1^2$ and $u_2^2$ also implies a “natural” definition of $u_1\,∂u_1/∂x_1$ and $u_2\,∂u_2/∂x_2$ computed with exact spectral derivatives, full dealiasing, and no interpolations.

Next up are the terms with $i=3$ for $j=1$ and $j=2$, which need to be compatible on the $\zeta_I$-nodes. Here we need

\[2\,\frac{u_j(\zeta_I^-)+u_j(\zeta_I^+)}{2} \, ∂_3(\zeta_I) \left(u_j(\zeta_I^+)-u_j(\zeta_I^-)\right) = ∂_3(\zeta_I) \left(u_j(\zeta_I^+)^2 - u_j(\zeta_I^-)^2\right)\]

which boils down to $(a+b)(a-b)=a²-b²$ and is therefore satisfied with the numerical expressions we are using.

This leaves the terms involving the vertical velocity component with $j=3$. Here, the numerical definition of $u_3^2$ on $\zeta_C$-nodes will inevitably involve interpolations, but it is not necessarily clear whether there is any such definition that is exactly equivalent to the expressions in the advection term for $i=1$ and $i=2$. These are evaluated on $\zeta_C$-nodes and interpolated after evaluating the product, so we can write

\[2 \, \frac{1}{2} \left( u_3(\zeta_C^-) \left. \frac{∂u_3}{∂x_i} \right|_{\zeta_C^-} + u_3(\zeta_C^+) \left. \frac{∂u_3}{∂x_i} \right|_{\zeta_C^+} \right) = \frac{1}{2} \left( \left. \frac{∂u_3^2}{∂x_i} \right|_{\zeta_C^-} + \left. \frac{∂u_3^2}{∂x_i} \right|_{\zeta_C^+} \right) = \frac{∂}{∂x_i} \frac{u_3(\zeta_C^-)^2 + u_3(\zeta_C^+)^2}{2}\]

i.e. the expressions are compatible (with spectral derivatives and full dealiasing) as long as we evaluate $u_3^2$ on $\zeta_C$-nodes by squaring first and interpolating second.

The only question that remains is what this implies for the numerical expression of $u_3\,∂u_3/∂x_3$ (on $\zeta_I$-nodes). We can work our way backwards from the $∂u_3^2/∂x_3$-term of the pressure gradient to obtain

\[\begin{aligned} \frac{∂u_3^2}{∂x_3} &= ∂_3(\zeta_I) \left( \frac{u_3(\zeta_I)^2 + u_3(\zeta_I^+)^2}{2} - \frac{u_3(\zeta_I^-)^2 + u_3(\zeta_I)^2}{2} \right) \\ &= \frac{∂_3(\zeta_I)}{2} \left( u_3(\zeta_I^+)^2 - u_3(\zeta_I^-)^2 \right) = \frac{∂_3(\zeta_I)}{2} \left( u_3(\zeta_I^+) + u_3(\zeta_I^-) \right) \left( u_3(\zeta_I^+) - u_3(\zeta_I^-) \right) \\ &= 2 \frac{ u_3(\zeta_I^-) + u_3(\zeta_I^+) }{2}\, ∂_3(\zeta_I) \left( \frac{u_3(\zeta_I) + u_3(\zeta_I^+)}{2} - \frac{u_3(\zeta_I^-) + u_3(\zeta_I)}{2} \right) = 2 u_3 \frac{∂u_3}{∂x_3} \end{aligned}\]

as the expression for the implied advection term. This can be understood as first interpolating $u_3$, then applying the $∂/∂x_3$-stencil, and finally multiplying with a “wide” interpolation of $u_3$. We can also think of the derivative as just the wider centered stencil with values a full (rather than a half) grid spacing away.

In summary, we can therefore say that for our numerical formulation on the staggered grid, the rotational advection term is exactly (to machine precision) equivalent to the convection form if we use full (quadratic) dealiasing everywhere, define the $u_j^2/2$-term on $\zeta_C$-nodes by first squaring each $u_j$ (with dealiasing) and then interpolating $u_3$, and define the “missing” $u_i ∂u_i/∂x_i$-terms of the convection form with the “obvious” formulation for $i=1$ and $i=2$ along with the above formulation for $i=3$. This last part is the only piece that is a bit less straightforward; otherwise everything is consistent in a rather clean and satisfying manner.

Discrete Conservation Properties

With the definition

\[2K\!E = \sum_{ζ_C} \left( u₁(ζ_C)² + u₂(ζ_C)² + \frac{u₃(ζ_C⁻)² + u₃(ζ_C⁺)²}{2} \right) = \sum_{ζ_C} u₁(ζ_C)² + \sum_{ζ_C} u₂(ζ_C)² + \sum_{ζ_I} u₃(ζ_I)²\]

we need

\[\frac{∂}{∂t} K\!E = \sum_{ζ_C} u₁ \frac{∂u₁}{∂t} + \sum_{ζ_C} u₂ \frac{∂u₂}{∂t} + \sum_{ζ_I} u₃ \frac{∂u₃}{∂t} = 0\]

for discrete energy conservation (to machine precision) in the absence of time-differencing errors. If grid stretching is used, we need to include the vertical grid spacing $Δζ \, ∂x₃/∂ζ$ in the sum. We get slightly different definitions depending on whether we evaluate the grid spacing at $ζ_C$ or $ζ_I$ for the $u₃²$-contribution.

Refer to Mansour, Moin, Reynolds & Ferziger 1979 for a general discussion of discrete conservation properties and how to derive them.

Analytically, we can easily show that the velocity–vorticity cross-product in the rotational form of the advection term conserves energy since

\[\frac{∂}{∂t} \frac{u_i^2}{2} = … − u_i u_j \left(\frac{∂u_i}{∂x_j} − \frac{∂u_j}{∂x_i}\right)\]

results in a net zero contribution as each term in the sum over $i$ and $j$ is either zero or appears with both positive and negative sign. The contribution is therefore zero at every point (i.e. locally) and therefore globally too.

For the discretized term, this is less obvious as the $ij$ and $ji$ contributions can have different discretization errors. With the numerical methods employed here, the global kinetic energy is conserved to machine precision as long as no grid stretching is used.

To demonstrate the net zero contribution on the staggered grid, we can rely on

\[\begin{aligned} &\sum_{ζ_I} u₃(ζ_I) \left( \frac{u₁(ζ_I⁻) + u₁(ζ_I⁺)}{2} \, ω₂(ζ_I) - \frac{u₂(ζ_I⁻) + u₂(ζ_I⁺)}{2} \, ω₁(ζ_I) \right) = \\ &\sum_{ζ_C} u₁(ζ_C) \frac{u₃(ζ_C⁻) \, ω₂(ζ_C⁻) + u₃(ζ_C⁺) \, ω₂(ζ_C⁺)}{2} - \sum_{ζ_C} u₂(ζ_C) \frac{u₃(ζ_C⁻) \, ω₁(ζ_C⁻) + u₃(ζ_C⁺) \, ω₁(ζ_C⁺)}{2} \end{aligned}\]

to show that the contributions still cancel in the global sum (assuming $u_3$ is zero at the boundaries).

Note that we also rely on results being unaffected by the order in which triple-products are computed, e.g. $u₁(ω₂u₃) = (u₁ω₂)u₃$. Due to truncated wavenumbers this is not true in general, but it does hold for the horizontal average, i.e. the $(0, 0)$ wavenumber pair, since the truncated wavenumbers of the first product are too large to contribute to the mean of the second product. The discrete energy conservation therefore holds independent of whether the advection term is dealiased or not.

In the presence of grid stretching, the $∂x₃/∂ζ$-factors do not cancel and the energy conservation is not exact.