Numerical Approach
In computational fluid dynamics, solving coupled systems of equations requires a robust numerical approach to handle the complexities of fluid flow and related physical phenomena. The following discussion outlines the semi-discretized formulation, preconditioning techniques, and temporal discretization used to solve such systems effectively.
Semi-discretized Form of the Coupled System
The governing equations for a coupled system are expressed as:
\[ \Gamma \frac{\partial Q}{\partial t} + \frac{\partial \left( F_c - F_d \right)_j}{\partial x_j} = G \]
Here, \( Q \) is the vector of primitive variables: \[ Q^T = [u, v, p, T, k, \varepsilon, Y_1, \ldots, Y_{N-1}], \] where \( u \) and \( v \) are velocity components, \( p \) is pressure, \( T \) is temperature, and \( Y_i \) are the species mass fractions.
The flux vector \( F_c \) represents convective and pressure terms, calculated using the AUSM+-up scheme in CMPS by defalt, while \( F_d \) accounts for diffusive terms. The source term vector \( G \) includes turbulent production and dissipation. The preconditioning matrix \( \Gamma \) modifies the system to improve numerical stability and convergence in low-Mach-number flows. The matrix is given as:
\[ \Gamma = \left[ \begin{matrix} \rho & 0 & \Theta u & \rho_T u & 0 & \cdots & \rho_{Y_1} u & \cdots & \rho_{Y_{N-1}} u \\ 0 & \rho & \Theta v & \rho_T v & 0 & \cdots & \rho_{Y_1} v & \cdots & \rho_{Y_{N-1}} v \\ 0 & 0 & \Theta & \rho_T & 0 & \cdots & \rho_{Y_1} & \cdots & \rho_{Y_{N-1}} \\ \rho u & \rho v & \Theta H - 1 & \rho_T H + \rho C_p & \cdots & \rho_{Y_1} H + \rho h_1 & \cdots & \rho_{Y_{N-1}} H + \rho h_{N-1} \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \end{matrix} \right] \]
In this matrix, \( \rho_T \) and \( \rho_{Y_k} \) are the partial derivatives of density with respect to temperature and species mass fractions, respectively. The term \( \Theta \) is defined as:
\[ \Theta = \left( \frac{1}{V_r^2} - \frac{\rho_T}{\rho C_p} \right), \] where \( V_r^2 \) is the local reference velocity. This velocity is determined by the following conditions:
\[ V_r = \begin{cases} 10^{-5} a & \text{if } |u| \leq 10^{-5} a, \\ |u| & \text{if } 10^{-5} a < |u| < a, \\ a & \text{if } |u| \geq a. \end{cases} \]
This formulation ensures numerical stability across a wide range of flow regimes. The AUSM+-up flux formulation, when combined with the preconditioning matrix, is particularly effective for steady subsonic problems.
Temporal Discretization
To advance the solution in time, the implicit temporal discretization of the coupled system is written as:
\[ \left[ \frac{\Gamma}{\Delta t} + \frac{\partial}{\partial Q} \left( \frac{\partial (F_c - F_d)_j}{\partial x_j} - G \right) \right]^n (Q^{n+1} - Q^n) = \left( \frac{\partial (F_c - F_d)_j}{\partial x_j} \right)^n + G^n. \]
The explicit version of this equation is given as:
\[ \frac{\Gamma^n}{\Delta t} (Q^{n+1} - Q^n) = \left( \frac{\partial (F_c - F_d)_j}{\partial x_j} \right)^n + G^n. \]
Integrating over a control volume \( \Omega_P \) and applying the divergence theorem, the equation becomes:
\[ \left[ \frac{1}{\Delta t} \int_{\Omega_P} \Gamma \, d\Omega + \frac{\partial}{\partial Q} \left( \oint_{\partial \Omega_P} (F_c - F_d) \cdot dS - \int_{\Omega_P} G \, d\Omega \right) \right]^n Q^* = \left( \oint_{\partial \Omega_P} (F_c - F_d) \cdot dS - \int_{\Omega_P} G \, d\Omega \right)^n. \]
The spatially discretized form is obtained as:
\[ \left[ \frac{\Omega_P \Gamma_P}{\Delta t} + \frac{\partial}{\partial Q} \left( \sum_{f \in nb(P)} \left( F_c - F_d \right)_f \cdot S_f - \Omega_P G_P \right) \right]^n Q^* = \left( \sum_{f \in nb(P)} \left( F_c - F_d \right)_f \cdot S_f - \Omega_P G_P \right)^n. \]
Here, \( f \) refers to the faces of the control volume \( \Omega_P \), \( nb(P) \) denotes neighboring cells, and \( S_f \) is the outward normal vector associated with each face. This formulation is essential for achieving accurate and stable solutions in finite-volume methods.