## 2 Subgrid Stress Tensor Equation Closure and Calibration
### Closure  
We are not aiming to construct the best possible closure model for unclosed terms in (3). Our goal is rather to optimize the coefficients of some typical form of closures based on DNS data. According to this, we use locally isotropic assumption, leading to  
\[2\
u\frac{\overline{\partial u_{i}^{\prime}}\overline{\partial u_{j}^{\prime}} }{\partial x_{k}}=\frac{2}{3}\varepsilon\delta_{ij},\quad\varepsilon=\
u \frac{\overline{\partial u_{i}^{\prime}}}{\partial x_{k}}\frac{\overline{ \partial u_{j}^{\prime}}}{\partial x_{k}},\]  
where \(\varepsilon\) is the subgrid kinetic energy dissipation rate. It is natural to express \(\varepsilon\) in terms of subgrid kinetic energy, \(k=\overline{u_{i}^{\prime}u_{j}^{\prime}}/2\), and characteristic filter size, \(\Delta\), as follows:  
\[\varepsilon=C_{E}\frac{k^{3/2}}{\Delta}. \tag{4}\]  
The simplest model assumes \(C_{E}=const\); however, due to the importance of accurate representation of the dissipation process, we will consider a dependence of \(C_{E}\) on invariants of subgrid stress anisotropy tensor \(\mathbf{a}=\big{(}a_{ij}\big{)}\):  
\[C_{E}=C_{E}^{(0)}+C_{E}^{(1)}A_{2}+C_{E}^{(2)}A_{3}. \tag{5}\]  
Here, \(a_{ij}\), \(A_{2}\) and \(A_{3}\) are defined as  
\[a_{ij}=\frac{\overline{u_{i}^{\prime}u_{j}^{\prime}}}{k}-\frac{2}{3}\delta_{ ij},\quad\ A_{2}=a_{i}\mu_{ji},\quad\ A_{3}=a_{ij}a_{jk}a_{ki}. \tag{6}\]  
Equation (5) can be considered as Taylor expansion of \(C_{E}\) in the corresponding variables.  
The pressure-strain correlation tensor, denoted hereinafter as \(\mathbf{\Phi}=\big{(}\Phi_{ij}\big{)}\), is usually expanded in powers of \(\mathbf{a}\). The most general form of linear model for \(\mathbf{\Phi}\) is the following [21]:  
\[\mathbf{\Phi}=-C_{\Phi 1}\mathbf{a}\varepsilon+k\bigg{[}C_{\Phi 2}\mathbf{S}+C_{\Phi 3}\bigg{(}\mathbf{a}\mathbf{S}+\mathbf{a}- \frac{2}{3}\mathbf{tr}(\mathbf{a}\mathbf{S})\mathbf{I}\bigg{)}+C_{\Phi 4}(\mathbf{\Omega a}-\mathbf{a}\mathbf{\Omega})\bigg{]}, \tag{7}\]  
with tensor multiplication obeying the matrix rule, \((\mathbf{A}\mathbf{B})_{ij}=A_{ik}B_{kj}\). Here, \(\mathbf{I}=\big{(}\delta_{ij}\big{)}\) is a unity tensor; \(\mathbf{S}\) and \(\mathbf{\Omega}\) are resolved strain rate and rotation tensors, respectively,  
\[S_{ij}=\frac{1}{2}\Bigg{(}\frac{\partial\overline{u}_{i}}{\partial x_{j}}+\frac {\partial\overline{u}_{j}}{\partial x_{i}}\bigg{)},\quad\ \ \Omega_{ij}=\frac{1}{2}\Bigg{(}\frac{\partial\overline{u}_{i}}{\partial x_{j}}- \frac{\partial\overline{u}_{j}}{\partial x_{i}}\Bigg{)}.\]  
Again, the simplest model assumes \(C_{\Phi 1},\dots,C_{\Phi 4}\) to be constants, but an inclusion of dependence on \(A_{2}\) and \(A_{3}\) may be considered as means to improve the model.  
For the triple-velocity correlation, we use a gradient model modified by means of resolved velocity gradients, the presence of which follows from exact the equation for \(\overline{u_{i}^{\prime}u_{j}^{\prime}u_{k}^{\prime}}\)[22], as follows:  
\[\overline{u_{i}^{\prime}u_{j}^{\prime}u_{k}^{\prime}}=\sum_{\{i\to j\to k \}}\bigg{[}-k^{1/2}\Delta(C_{T1}\delta_{kl}+C_{T2}a_{kl})+C_{T3}\Delta^{2} \frac{\partial\overline{u}_{k}}{\partial x_{l}}\bigg{]}\frac{\overline{u_{i} ^{\prime}u_{j}^{\prime}}}{\partial x_{l}}. \tag{8}\]Summation is performed over cyclic permutations of the \(i,j\) and \(k\) indices:  
\[\sum_{\{i\to j\to k\}}A_{ijk}=A_{ijk}+A_{jki}+A_{kij}.\]  
The coefficients \(C_{T1},C_{T2},C_{T3}\) can be either constants or functions of \(A_{2}\) and \(A_{3}\).  
It remains to consider the pressure-velocity correlation. According to [23], one can approximate it by contracting the triple-velocity correlation tensor:  
\[\frac{1}{\rho^{\prime}u_{i}^{\prime}}=-C_{PT}\overline{u_{i}^{\prime}u_{m}^{ \prime}u_{m}^{\prime}}. \tag{9}\]  
There is one empirical coefficient: \(C_{PT}\).