Numerical Method

Generalized Curvilinear Coordinates

TOSCA solves the governing equations in generalized curvilinear coordinates \((l_1, l_2, l_3) \equiv (\xi, \eta, \zeta)\), related to the Cartesian coordinates \((x_1, x_2, x_3)\) through a differentiable, invertible mapping \(x_i = x_i(\xi,\eta,\zeta)\). The Jacobian of the transformation is

\[J = \det\!\left(\frac{\partial x_i}{\partial l_q}\right).\]

The face area vectors are defined as

\[S^q_i = \frac{1}{J}\,\frac{\partial l_q}{\partial x_i}, \qquad q,i \in \{1,2,3\}.\]

Geometrically, \(\mathbf{S}^q \equiv (S^q_1, S^q_2, S^q_3)\) is the outward unit normal of the \(q\)-th cell face scaled by the ratio of face area to cell volume. The contravariant flux through the \(q\)-th face is

\[V^q = S^q_i\, u_i \equiv \,\mathbf{S}^q \cdot \mathbf{u}.\]

The contravariant metric tensor, multiplied by a factor \(1/J^2\) for numerical convenience, is the symmetric positive-definite tensor

\[S^{qr} = S^q_i S^r_i.\]

Differential Operators in Curvilinear Coordinates

The four fundamental operators used throughout the discretization are:

Gradient of a scalar \(\phi\):

\[\frac{\partial \phi}{\partial x_i} = J S^k_i\,\frac{\partial \phi}{\partial l_k}\]

Gradient of a vector field \(\mathbf{a}\):

\[\frac{\partial a_j}{\partial x_i} = J S^k_i\,\frac{\partial a_j}{\partial l_k}\]

Divergence of a vector field \(\mathbf{a}\):

\[\frac{\partial a_i}{\partial x_i} = J\,\frac{\partial}{\partial l_k}\!\left(J S^k_i a_i\right) = J\,\frac{\partial A^k}{\partial l_k},\]

where \(A^k = S^k_i a_i\) is the contravariant flux of \(\mathbf{a}\). For an incompressible velocity field this is equal to \(\partial V^q/\partial l_q = 0\).

Divergence of a rank-2 tensor \(\tau_{ij}\):

\[\frac{\partial}{\partial x_i}\left(\tau_{ij}\right) = J\,\frac{\partial}{\partial l_k} \left(S^k_i \tau_{ij}\right) = J\,\frac{\partial}{\partial l_k} T^k_j,\]

The quantity \(T^k_j\) is the \(k\)-th contravariant flux of the \(j\)-th row of the tensor. Using the above formula, the divergence of the advection term in the momentum equation can be written as:

\[\frac{\partial}{\partial x_i}\!\left(u_i u_j\right) = J\,\frac{\partial}{\partial l_k} \left(S^k_i u_i u_j\right) = J\,\frac{\partial}{\partial l_k} \left(V^k u_j\right).\]

Using the above rules, the full viscous term in the momentum equation can be rewritten as:

\[\frac{\partial}{\partial x_j}\left[\nu_\text{eff}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\right] = J\frac{\partial}{\partial l_k}\left[\nu_\text{eff}\left(JS^{kn}\frac{\partial u_i}{\partial l_n} + JS^k_jR_{ij}\right)\right].\]

where \(R_{ij}\) is a non-symmetric tensor defined as \(R_{ij}=S^m_i\partial u_j/\partial l_m\), and \(S^{kn}\) is the contravariant metric tensor evaluated with the face area vectors, namely \(S^{kn} = S^k_i S^n_i = 1/J^2 \partial l_k/\partial x_i \partial l_n / \partial x_i\).

Governing Equations in Curvilinear Form

Applying the operators above, the governing equations solved by TOSCA take the partially-transformed form (Ge & Sotiropoulos, 2007). The momentum equation is written for the contravariant flux \(V^q\), which eliminates Christoffel-symbol correction terms and maintains second-order accuracy on any smoothly varying mesh. The system reads:

\begin{align} \frac{\partial V^q}{\partial l_q} &= 0, \\[6pt] \frac{\partial V^q}{\partial t} + J S^q_i \frac{\partial}{\partial l_r}\!\left(V^r u_i\right) &= -\frac{J^2}{\rho_\text{ref}}\frac{\partial p}{\partial l_r}S^{rq} + J S^q_i \frac{\partial}{\partial l_r} \!\left[J\nu_\text{eff}\!\left(S^{rn}\frac{\partial u_i}{\partial l_n} + S^k_j R_{ij}\right)\right] \\ &\quad - \frac{J^2}{\rho_\text{ref}}\frac{\partial p_\infty}{\partial l_r}S^{rq} + S^q_i\frac{\rho_k}{\rho_\text{ref}}g_i - 2 S^q_i\epsilon_{ijk}\Omega_j u_k \\ &\quad + S^q_i\!\left(f_i + s^v_i + s^h_i\right), \\[6pt] \frac{\partial\theta}{\partial t} + \frac{\partial}{\partial l_r}\!\left(V^r\theta\right) &= J\frac{\partial}{\partial l_r} \!\left(J\kappa_\text{eff}S^{rk}\frac{\partial\theta}{\partial l_k}\right) + s_\theta, \end{align}

where \(R_{ij} = S^k_j\,\partial u_i/\partial l_k\) is the transpose-gradient term that ensures symmetry of the viscous flux, \(p_\infty\) is the large-scale background driving pressure, \(\rho_k/\rho_\text{ref}\) is the Boussinesq density ratio, and \(f_i\), \(s^v_i\), \(s^h_i\), \(s_\theta\) are body forces from turbines, vertical damping, horizontal damping, and the temperature controller, respectively (see Governing Equations).

Staggered Grid Arrangement

TOSCA employs a hybrid collocated–staggered discretization (Ge & Sotiropoulos, 2007). Contravariant fluxes \(V^q\) are stored at the center of their respective cell face (staggered), ensuring exact discrete mass conservation. Pressure \(p\), potential temperature \(\theta\), and the sub-grid viscosity \(\nu_t\) are stored at cell centers (collocated).

Cartesian velocity \(\mathbf{u}\) is required at cell faces for the advection and viscous operators. Rather than solving three Cartesian momentum equations at each face (which would triple the cost), TOSCA first assembles the momentum right-hand side at cell centers and then interpolates it to the faces using second-order linear interpolation along each grid direction. After each step, Cartesian velocity is reconstructed at cell centers from the updated contravariant fluxes by inverting \(V^q = S^q_i u_i\).

Fractional-Step Projection Method

Pressure–velocity coupling uses a second-order fractional-step projection method. Given the state \((\mathbf{V}^n, p^n, \theta^n)\) at time level \(n\), one time step proceeds through the four stages below.

Stage 1 — Momentum Predictor:

Advance the contravariant fluxes to a predicted state \(\mathbf{V}^{*}\) that satisfies the momentum equation but not yet the divergence-free constraint. Schematically, assuming to choose the Crank–Nicolson scheme (other time-stepping schemes, described in Time Integration Schemes, follow in a similar manner), the predictor step reads:

\[\frac{V^{q,*} - V^{q,n}}{\Delta t} = \frac{1}{2}\!\left[H^q_U(\mathbf{V}^{*},\mathbf{u}^{*}) + H^q_U(\mathbf{V}^n,\mathbf{u}^n)\right] + S^{q,n} + G^{q,n} + B^{q,*}.\]

Here \(H^q_U\) contains the advection and viscous fluxes (see Advection Schemes); \(S^q\) contains buoyancy, Coriolis, turbine, and damping source terms; and \(G^q\) is the discrete pressure-gradient operator. The buoyancy source \(B^{q}_U\) in the momentum equation depends on the face-interpolated potential temperature through the Boussinesq density ratio. Because \(\theta\) is solved only in Stage 4, it cannot be used implicitly in Stage 1. Hence, before the momentum predictor step, the buoyancy body-force vector \(B^{q,n}\) is computed from the current temperature field \(\theta^n\) and stored. The previous step’s buoyancy \(B^{q,n-1}\) is also retained, so that buoyancy can be extrapolated to the mid-step level \(n+\tfrac{1}{2}\) using an Adams–Bashforth 2 formula, which reads:

\[B^{q,*} = \tfrac{3}{2}\,B^{q}(\theta^n) - \tfrac{1}{2}\,B^{q}(\theta^{n-1}).\]

This approach reduces the temporal coupling error to \(\mathcal{O}(\Delta t^2)\) without the need to perform a temperature predictor step. On the first time step, a simple Forward Euler formula is used.

Stage 2 — Pressure Correction:

The scalar correction \(\phi^{n+1}\) is calculated as:

\[D\!\left(G^q(\phi^{n+1})\right) = -\frac{1}{\Delta t}\,D\!\left(V^{q,*}\right),\]

where the discrete divergence operator at cell \((i,j,k)\) is

\[D_{ijk}(\gamma) = \gamma_{i+\frac{1}{2},j,k} + \gamma_{i,j+\frac{1}{2},k} + \gamma_{i,j,k+\frac{1}{2}} - \gamma_{i-\frac{1}{2},j,k} - \gamma_{i,j-\frac{1}{2},k} - \gamma_{i,j,k-\frac{1}{2}}.\]

Substituting the pressure-gradient operator into the divergence gives a sparse 19-point linear system per cell (1 diagonal, 6 face-centred, 12 edge-centred cross-derivative contributions from \(S^{qr}\)). The system is solved with HYPRE’s BoomerAMG algebraic multigrid preconditioner (GMRES or PCG outer Krylov method), providing \(\mathcal{O}(N)\) scalable performance. The constant mode of \(\phi\) is removed after the solve by subtracting the domain-averaged value, ensuring that the pressure is defined up to a constant.

Stage 3 — Contravariant Flux Projection:

The predicted contravariant fluxes are corrected to enforce mass conservation:

\[V^{q,n+1} = V^{q,*} + \Delta t\, G^q(\phi^{n+1}).\]

The pressure is updated as \(p^{n+1} = p^n + \phi^{n+1}\).

Stage 4 — Temperature Solution:

With the divergence-free \(\mathbf{V}^{n+1}\) available, \(\theta^n\) is advanced to \(\theta^{n+1}\) using the selected temperature time-integration scheme (see Time Integration Schemes).