The equations of motion

The large-scale model

To be verbose and specific, SOMAR solves the incompressible, Boussinesq Navier-stokes equations in the following form.

\begin{align} \frac{\partial u}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i u + 2 J \left( \nu + \nu_{sgs} \right) S^{i x} \right] - \frac{\partial \xi}{\partial x} \frac{\partial p}{\partial \xi} + fv + F_u \\[1em] \frac{\partial v}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i v + 2 J \left( \nu + \nu_{sgs} \right) S^{i y} \right] - \frac{\partial \eta}{\partial y} \frac{\partial p}{\partial \eta} -fu + F_v \\[1em] \frac{\partial w}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i w + 2 J \left( \nu + \nu_{sgs} \right) S^{i z} \right] - \frac{\partial \zeta}{\partial z} \frac{\partial p}{\partial \zeta} - b' + F_w \\[1em] \frac{\partial T}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i T + J \left( \kappa_T + \kappa_{sgs,T} \right) \frac{\partial T'}{\partial \xi^i} \right] + F_T \\[1em] \frac{\partial S}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i S + J \left( \kappa_S + \kappa_{sgs,S} \right) \frac{\partial S'}{\partial \xi^i} \right] + F_S \\[1em] \frac{\partial \lambda}{\partial t} & = \frac{1}{J} \sum_i \frac{\partial}{\partial \xi^i}\left[ -U^i \lambda + J \left( \kappa_{\lambda} + \kappa_{sgs,\lambda} \right) \frac{\partial \lambda}{\partial \xi^i} \right] + F_{\lambda} \\[1em] 0 & = \frac{\partial U}{\partial \xi} + \frac{\partial V}{\partial \eta} + \frac{\partial W}{\partial \zeta} \\[1em] b' & = \Gamma(T,S) - \Gamma(\overline{T}, \overline{S}) \end{align}

The \((\xi,\eta,\zeta)\) are logical coordinates and the \((x,y,z)\) are physical coordinates. The logical and physical coordinates are equal unless using a stretched coordinate system. See Coordinate systems and Stretching the grids for more details.

All primed quantities represent deviations from a static vertical stratification. That is, \(T'(x,y,z,t) = T(x,y,z,t) - \overline{T}(z)\), and so on. When a stratification is provided, splitting the active scalars avoids diffusion of the background state and avoids repeated computation of the dominant hydrostatic pressure component.

\(S^{i \alpha}\) is the rate-of-strain tensor, made curvilinear-based in the first index and Cartesian-based in the second. Also, the capital velocity components represent the divergence-free curvilinear-based advecting velocity components, \(\left(U,V,W\right) = \left(J \frac{\partial \xi}{\partial x} u, J \frac{\partial \eta}{\partial y} v, J \frac{\partial \zeta}{\partial z} w \right)\).

By default, the custom forcing functions, \(F_{\bullet}\) are zero, but these can be modified to fit your simulation’s needs. SOMAR also allows for custom passive scalars, \(\lambda\), which can be made active by modifying the custom forcing functions. Although we only show one such scalar, SOMAR allows for any number of them.

Warning

SOMAR uses a projection method to enforce incompressibility. This means the variable \(p\) is technically not the physical pressure. It is whatever it needs to be to make the advecting velocity divergence-free. In the domain’s interior, \(p\) is a reasonable, low-order estimate of the pressure. However, it is computed using unphysical homogeneous Neumann boundary conditions, causing it to deviate from reality in boundary layers!

The subgrid-scale model

Without an LES, the sub-grid scale data \(\nu_{sgs},\,\kappa_{sgs,\bullet}\) are identically zero. When using an LES, these are computed as decribed in Ducros, et. al.

\begin{align} \nu_{sgs} &= 0.0014\, C_K^{-3/2}\, \overline{\Delta}\, \sqrt{\overline{F_2}} \\[1em] \kappa_{sgs,\bullet} &= \nu_{sgs} / \text{Pr}_{sgs, \bullet} \\[1em] \overline{F_2} &= \frac{1}{6} \left[ \left( \tilde{u}_{i+1,j,k} - \tilde{u}_{i,j,k} \right)^2 + \left( \tilde{u}_{i-1,j,k} - \tilde{u}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{u}_{i,j+1,k} - \tilde{u}_{i,j,k} \right)^2 + \left( \tilde{u}_{i,j-1,k} - \tilde{u}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{u}_{i,j,k+1} - \tilde{u}_{i,j,k} \right)^2 + \left( \tilde{u}_{i,j,k-1} - \tilde{u}_{i,j,k} \right)^2 \right] \\ & + \frac{1}{6} \left[ \left( \tilde{v}_{i+1,j,k} - \tilde{v}_{i,j,k} \right)^2 + \left( \tilde{v}_{i-1,j,k} - \tilde{v}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{v}_{i,j+1,k} - \tilde{v}_{i,j,k} \right)^2 + \left( \tilde{v}_{i,j-1,k} - \tilde{v}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{v}_{i,j,k+1} - \tilde{v}_{i,j,k} \right)^2 + \left( \tilde{v}_{i,j,k-1} - \tilde{v}_{i,j,k} \right)^2 \right] \\ & + \frac{1}{6} \left[ \left( \tilde{w}_{i+1,j,k} - \tilde{w}_{i,j,k} \right)^2 + \left( \tilde{w}_{i-1,j,k} - \tilde{w}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{w}_{i,j+1,k} - \tilde{w}_{i,j,k} \right)^2 + \left( \tilde{w}_{i,j-1,k} - \tilde{w}_{i,j,k} \right)^2 \right. \\ & \hspace{2em} \left. + \left( \tilde{w}_{i,j,k+1} - \tilde{w}_{i,j,k} \right)^2 + \left( \tilde{w}_{i,j,k-1} - \tilde{w}_{i,j,k} \right)^2 \right] \\[1em] \tilde{\varphi} &= \mathscr{L}^3\left[ \varphi \right] \\[1em] \mathscr{L}\left[\varphi\right] &= \left( \varphi_{i+1,j,k} -2 \varphi_{i,j,k} + \varphi_{i-1,j,k} \right) \\ & + \left( \varphi_{i,j+1,k} -2 \varphi_{i,j,k} + \varphi_{i,j-1,k} \right) \\ & + \left( \varphi_{i,j,k+1} -2 \varphi_{i,j,k} + \varphi_{i,j,k-1} \right), \end{align}

where \(C_K=0.5\), \(\overline{\Delta}=(J \, \Delta\xi \, \Delta\eta \, \Delta \zeta)^{1/3}\), and each scalar’s subgrid Prandtl numbers are provided by the user through the input file.

When using multiple levels of adaptive refinement, the sub-grid parameterizations \(\nu_{sgs}\) and \(\kappa_{sgs,\bullet}\) can be computed at all levels, or only at the finest then interpolated down to coarser levels.

Coordinate systems

SOMAR bounces back and forth between two coordinate systems, logical and physical. The logical coordinates of each node are always defined by

\[\begin{split}\xi_i &= i \Delta\xi &= i L_x/N_x, \\ \eta_j &= j \Delta\eta &= j L_y/N_y, \\ \zeta_k &= k \Delta\zeta &= k L_z/N_z.\end{split}\]

Where \((i,j,k)\) are the nodal indices, \(L_{\bullet}\) is the domain’s length and \(N_{\bullet}\) is the number of cells in each direction. The logical coordinates at cell-centers are defined similarly, but with 1/2 added to the indices.

The physical coordinates are computed from the logical coordinates via three stretching functions,

\[\left( x, y, z \right) = \left( f(\xi), g(\eta), h(\zeta) \right).\]

The identity map, \(\left( x, y, z \right) = \left( \xi, \eta, \zeta \right)\), leads to a simple, unstretched Cartesian coordinate system. The user can customize these stretching functions to meet the needs of the simulation.

The Jacobian matrix is \(\rm{diag}\left( \frac{\partial x}{\partial \xi}, \frac{\partial y}{\partial \eta}, \frac{\partial z}{\partial \zeta} \right)\) and it’s determinant is \(J = \frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \eta} \frac{\partial z}{\partial \zeta}\). With this system, volume integrals of scalar fields are, for example,

\[M = \iiint \rho\,J\,d\xi\,d\eta\,d\zeta = \sum_{i,j,k} \left( \rho \, J \right)_{i,j,k} \, \Delta\xi \Delta \eta \Delta \zeta.\]