\begin{center}\underline{\Huge \bf Introduction}\end{center} \vspace{1cm} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Aims of the document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This chapter constitutes an introduction to the theory and developer's guide associated with the kernel of \CS\footnote{With some slight modifications, the algorithms are the same in versions~1.0, 1.1 and 1.2.}. The system of equations considered consists of the Navier-Stokes equations, with turbulence and passive scalars. To begin with, the continuous equations for mass, momentum, turbulence and passive scalars are presented. Then, information related to the time scheme is supplied. Finally, the spatial discretization is detailed: it is based on a co-located\footnote{All the variables are located at the centres of the cells.} finite volume scheme for unstructured meshes. To make the documentation immediately suitable to the developers' needs, it has been organized into sub-sections corresponding to the major steps of the algorithm and to some important subroutines of the code. In each sub-section (for each subroutine), the {\bf function} of the piece of code considered is indicated. Then, a description of the {\bf discretization} follows. To finish with, and more oriented towards the developers, details of the {\bf implementation} are provided and a list of open problems is given (improvements, limitations...). Several accessibilty levels have been defined for the documentation, so that it is possible to choose the information level best suited to the intended use. At the present time, on UNIX and Linux plat-forms, free access is granted to all information (level called "\texttt{Complet}"). The most restrictive level provides only the function of the subroutines ("\texttt{Fonction}" level). The intermediate level provides the function of the subroutines and the associated discretization ("\texttt{Discret}" level). % fa modification (no more 1.1.0.m, no more "only basic") %The current version of this documentation details the basic algorithms %of the code and it was therefore important to be well explained by an EDF report. During the development process of the code, the documentation is naturally updated as and when required by the evolution of the source code itself. Suggestions for improvement are welcome. In particular, it will be necessary to deal with some transverse subjects (such as parallelism, periodicity or memory management) which were voluntarily left out of the first versions, to concentrate efforts on the algorithms and their implementation. To make it easier for the developers to keep the documentation up to date during the development process, the files have been associated "physically" with the release of the code (each release of the code includes a directory containing the whole documentation). In practice, % fa modification : discard (from the development version 1.1.0.n on), the users of \CS\ will find the documentation at the following location % fa modification (MFTT-MFEE) (UNIX and Linux server at MFEE): % fa modification (SC -> CS) \begin{center} \texttt{\$CS\_HOME/doc/NOYAU/Postscripts/Base}, \end{center} The general command \texttt{info\_cs [noyau]} also provides this information. % fa modification (useless now) %(also in \texttt{Netscape},alwzys in the Unix MFTT server, % the page which resends the additional documents). The current report will show a %number of very limited examples and it is more favourable to consult the online % documentation. \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Continuous equations} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This section presents the continuous equations. It does not substitutes for the specific sub-sections of this documentation: the purpose here is mainly to provide an overview before more detailed reading. In the following, $\rho$ stands for the density, $\vect{u}$ for the velocity field and $\phi^i$ for the value of the variable $\phi$ injected when an additional mass source term is present ($\Gamma\, \phi^i$) on the right-hand side of the equations. %\clearpage {\bf Mass}\\ \begin{equation} \dive(\rho \vect{u}) = 0 \end{equation} In fact, to compute a given unknown $\phi$ (and in particular for the velocity prediction), the equation $\displaystyle \frac{\partial \rho}{\partial t} + \dive(\rho \vect{u}) = \Gamma$ is used to re-write the term $\displaystyle \frac{\partial \rho\,\phi}{\partial t}$ as follows: \mbox{$\displaystyle \rho\frac{\partial \phi}{\partial t} - \phi\,\dive(\rho \phi) + \Gamma\,\phi$}. The possible variations in time of the density are however not taken into account in the pressure correction step. \clearpage {\bf Momentum} \begin{equation} \left\{\begin{array}{l} \displaystyle\frac{\partial}{\partial t}(\rho \vect{u}) + \dive(\rho\, \vect{u} \otimes \vect{u}) =\dive(\tens{\sigma}) + \vect{TS} - \tens{K}\,\vect{u}\\ \dive(\rho \vect{u}) = 0 \end{array}\right. \end{equation} $[\,\vect{TS}-\tens{K}\,\vect{u}\,]$ stands for the additional momentum source terms which may be prescribed by the user (head loss, $\Gamma \vect{u}^i$ contribution associated with a user-prescribed mass source term...). $\tens{K}$ is a positive diagonal tensor, by definition (derived, for example, from the diagonal source term of the head loss tensor). $\bullet$ For laminar flow, $\tens{\sigma}$ is the stress tensor: \begin{equation} \tens{\sigma} = \tens{\tau} - P\tens{Id} \end{equation} where $\tens{\tau}$ is the viscous stress tensor, defined from $\mu=\mu_l$ (dynamic molecular viscosity) and $\tens{D}$ (strain rate tensor) as: \begin{equation} \tens{\tau} = 2\,\mu\ \tens{D} - \frac{2}{3}\mu\ tr({\tens{D}})\ \tens{Id} \text{ \hspace*{1cm}with \hspace*{1cm} } \tens{D} = \frac{1}{2}(\ggrad\vect{u}+\,^{t}\ggrad\vect{u}) \label{tau_eq}\end{equation} $\bullet$ For turbulent flow, $\displaystyle \tens{\sigma}$ also accounts for the turbulent Reynolds stress tensor (correlations of the velocity fluctuations arising from the non linear convective term). The modelling of the latter depends upon the turbulence model adopted: \hspace*{1cm} - with a $k-\varepsilon$ model, the closure requires a turbulent viscosity $\mu_t$. With the same definition for $\tens{\tau}$ as in equation~(\ref{tau_eq}), but with $\mu=\mu_l+\mu_t$, $\tens{\sigma}$ reads: \begin{equation} \tens{\sigma} = \tens{\tau} - (P+\frac{2}{3}\rho k)\tens{Id} \end{equation} \hspace*{1cm} - with a Reynolds stress model, the components of the Reynolds stress tensor $\tens{R}$ are variables of the system of equations considered, and so one obtains, with $\mu=\mu_l$ in the definition of $\tens{\tau}$ (equation~(\ref{tau_eq}))~: \begin{equation} \tens{\sigma} = \tens{\tau} - P\tens{Id} -\rho\tens{R} \end{equation} \hspace*{1cm} - with a $LES$ approach, the mixing length hypothesis is applied and leads to a modelling of the sub-grid effects through a turbulent viscosity $\mu_t$: the definition for $\tens{\sigma}$ remains the same as that used for laminar flow (equation~(\ref{tau_eq})), but with $\mu=\mu_l+\mu_t$. \clearpage {\bf Equations for the variables $k$ and $\varepsilon$ (standard $k-\varepsilon$ model)} \begin{equation} \left\{\begin{array}{l} \displaystyle %\rho\frac{\partial k}{\partial t} + \frac{\partial }{\partial t} (\rho k) + \dive\left[\rho \vect{u}\,k-(\mu+\frac{\mu_t}{\sigma_k}\grad{k})\right] = \mathcal{P}+\mathcal{G}-\rho\varepsilon %+k\dive(\rho\vect{u}) %+\Gamma(k_i-k)\\ +\Gamma k^i %\multicolumn{1}{c}{+\alpha_k k +\beta_k}\\ +TS_k\\ \displaystyle %\rho\frac{\partial \varepsilon}{\partial t} + \frac{\partial }{\partial t} (\rho\varepsilon) + \dive\left[\rho \vect{u}\,\varepsilon- (\mu+\frac{\mu_t}{\sigma_\varepsilon}\grad{\varepsilon})\right] = C_{\varepsilon_1}\frac{\varepsilon}{k}\left[\mathcal{P} +(1-C_{\varepsilon_3})\mathcal{G}\right] -\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k} %\varepsilon\dive(\rho\vect{u})+\\ %\multicolumn{1}{c}{+\Gamma(\varepsilon_i-\varepsilon)+\alpha_\varepsilon \varepsilon +\beta_\varepsilon} %\multicolumn{1}{c}{+\Gamma \varepsilon_i+TS_\varepsilon} +\Gamma \varepsilon^i+TS_\varepsilon \end{array}\right. \end{equation} $\mathcal{P}$ is the production term created by mean shear: \begin{displaymath} \begin{array}{rcl} \mathcal{P} & = & \displaystyle -\rho R_{ij}\frac{\partial u_i}{\partial x_j} = -\left[-\mu_t \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) +\frac{2}{3}\mu_t\frac{\partial u_k}{\partial x_k}\delta_{ij} +\frac{2}{3}\rho k\delta_{ij}\right] \frac{\partial u_i}{\partial x_j}\\ & = & \displaystyle \mu_t \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j} -\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})\\ & = & \displaystyle \mu_t \left[ 2\left(\frac{\partial u}{\partial x}\right)^2+ 2\left(\frac{\partial v}{\partial y}\right)^2+ 2\left(\frac{\partial w}{\partial z}\right)^2+ \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2+ \left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)^2+ \left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2 \right]\\ \multicolumn{3}{r}% {\displaystyle -\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})} \end{array} \end{displaymath} $\mathcal{G}$ is the production term created by gravity effects: $\displaystyle \mathcal{G}=-\frac{1}{\rho}\frac{\mu_t}{\sigma_t} \frac{\partial\rho}{\partial x_i}g_i$ The turbulent viscosity is: $\displaystyle \mu_t=\rho C_\mu\frac{k^2}{\varepsilon}$. $TS_\varphi$ ($\varphi=k$ or $\varepsilon$) stands for the additional source terms prescribed by the user (in rare cases only). The constants of the model are given below:\\ \begin{table}[htp] \begin{center} \begin{tabular}{|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|} \hline $C_\mu$ & $C_{\varepsilon_1}$ & $C_{\varepsilon_2}$ & $\sigma_k$ & $\sigma_\varepsilon$ \\ \hline $0,09$ & $1,44$ & $1,92$ & $1$ & $1,3$\\ \hline \end{tabular} \end{center} \end{table} $C_{\varepsilon_3}=0$ if $\mathcal{G}\geqslant0$ (unstable stratification) and $C_{\varepsilon_3}=1$ if $\mathcal{G}\leqslant0$ (stable stratification). \clearpage {\bf Equations for the Reynolds stress tensor components $R_{ij}$ and $\varepsilon$ (LRR model)} \begin{equation} \left\{\begin{array}{lll} \displaystyle \frac{\partial }{\partial t} (\rho R_{ij}) + \dive(\rho \vect{u}\,R_{ij} - \mu \grad{R_{ij}}) & = \mathcal{P}_{ij} + \mathcal{G}_{ij}+\Phi_{ij} + \it{d}_{ij} - \varepsilon_{ij} & \displaystyle + \Gamma R^{i}_{ij} + TS_{R_{ij}} \\ \displaystyle \frac{\partial }{\partial t}(\rho\varepsilon) + \dive\left[\rho \vect{u}\,\varepsilon- (\mu \grad{\varepsilon})\right] & = \displaystyle \it{d} + C_{\varepsilon_1}\frac{\varepsilon}{k}\left[\mathcal{P} +\mathcal{G}_{\varepsilon}\right] -\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k} & \displaystyle +\Gamma \varepsilon^{i} +TS_\varepsilon \end{array}\right. \end{equation} $\tens{\mathcal{P}}$ stands for the turbulence production tensor associated with mean shear and $\tens{\mathcal{G}}$ stands for the turbulence production tensor associated with gravity effects: \begin{equation} \displaystyle P_{ij} = \displaystyle -\rho \left[ R_{ik} \frac{\partial u_j}{\partial x_k} + R_{jk} \frac{\partial u_i}{\partial x_k} \right]\text{\hspace*{1cm}and\hspace*{1cm}} \mathcal{G}_{ij}= \left[ G_{ij} - C_3 (G_{ij}-\frac{1}{3} \delta_{ij} G_{kk}) \right] \end{equation} \begin{equation} \text{where\hspace*{0,3cm}} \displaystyle G_{ij} = - \frac{3}{2} \frac{C_{\mu}}{\sigma_{t}} \frac{k}{\varepsilon} (r_i g_j + r_j g_i) \text{\hspace*{0,3cm}with\hspace*{0,3cm}} \displaystyle k = \frac{1}{2} R_{ll} \text{\hspace*{0,3cm}and\hspace*{0,3cm}} \displaystyle r_i = R_{ik} \frac{\partial \rho}{\partial x_k} \end{equation} moreover, the following notations are introduced: \begin{equation} \mathcal{P} =\frac{1}{2} \mathcal{P}_{kk}\text{\hspace*{1cm}and\hspace*{1cm}} \mathcal{G}_{\varepsilon} = max(0,\frac{1}{2}G_{kk}) \end{equation} $\tens{\Phi}$ is the term associated with pressure-velocity correlations: \begin{equation} \displaystyle \Phi_{ij} - \varepsilon_{ij} = \phi_{ij,1} + \phi_{ij,2} + \phi_{ij,w} -\frac{2}{3} \rho \ \delta_{ij} \varepsilon \end{equation} \begin{equation} \text{with\hspace*{0,3cm}} \phi_{ij,1} = -\rho\,C_1 \frac{\varepsilon}{k} (R_{ij} - \frac{2}{3} k \delta_{ij})\text{\hspace*{1cm}and\hspace*{1cm}} \phi_{ij,2} = -\rho\,C_2 (\mathcal{P}_{ij} - \frac{2}{3} \mathcal{P}) \delta_{ij} \end{equation} The term $\phi_{ij,w}$ is called ``wall echo term'' (by default, it is not accounted for~: see \fort{turrij}). The turbulent diffusion terms are~: \begin{equation} \it{d}_{ij} = C_{S} \frac{\partial}{\partial x_k} (\rho \frac{k}{\varepsilon} R_{km} \frac{\partial R_{ij}}{\partial x_m})\text{\hspace*{1cm}and\hspace*{1cm}} \it{d} = C_{\varepsilon} \displaystyle \frac{\partial}{\partial x_k} \left( \rho \frac{k}{\varepsilon} R_{km} \frac{\partial \varepsilon}{\partial x_m} \right) \end{equation} $TS_\varphi$ stands for the additional source terms prescribed by the user (in rare cases only) for the variable $\phi$ ($R_{ij}$ or $\varepsilon$). The constants of the model are given below:\\ %\begin{table}[htp] \begin{center} \begin{tabular}{|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|p{0,8cm}|} \hline $C_\mu$ & $C_{\varepsilon}$ & $C_{\varepsilon_1}$ & $C_{\varepsilon_2}$ & $C_1$ & $C_2$ & $C_3$& $C_S$ & $C'_1$ & $C'_2$\\ \hline $0,09$ & $0,18$ & $1,44$ & $1,92$ & $1,8$ & $0,6$ & $0,55$& $0.22$ & $0,5$ & $0,3$\\ \hline \end{tabular} \end{center} %\end{table} \clearpage {\bf Definition of $\mu_t$ for $LES$} {\tiny$\bigstar$} \underline{Smagorinsky model} \begin{equation} \mu_{t} = \rho\,(C\,\overline{\Delta})^{2} \sqrt{2\overline{D_{ij}}\,\overline{D_{ij}}} \end{equation} With $\overline{D_{ij}}$ the filtered strain rate tensor components: \begin{equation} \overline{D_{ij}} = \frac{1}{2}\left( \frac{\partial \,\overline{u_{i}}}{\partial \, x_{j}}+\frac{\partial \,\overline{u_{j}}}{\partial \, x_{i}}\right) \end{equation} Here, $\overline{u_i}$ stands for the $i^{th}$ resolved velocity component\footnote{In the case of implicit filtering, the discretization in space introduces a spectral low pass filter: only the structures larger that twice the size of the cells are accounted for. Those structures are called "the resolved scales", whereas the rest, $u_i-\overline{u_i}$, is referred to as "unresolved scales" or "sub-grid scales". The influence of the unresolved scales on the resolved scales have to be modelled.}. \\\\ $C$ is the Smagorinsky constant. Its theoretical value is $0.18$ for homogenous isotropic turbulence, but the value $0.065$ is classic for channel flow. \\\\ $\overline{\Delta}$ is the filter width associated with the finite volume formulation (implicit filtering which corresponds to the integration over a cell). The value recommended for hexahedral cells is: \begin{equation} \overline{\Delta} = 2\,|\Omega|^{\frac{1}{3}}\\ \end{equation} where $|\Omega|$ is the volume of the cell. {\tiny$\bigstar$} \underline{Classic dynamic model} \\\\ A second filter is introduced: it is an explicit filter with a characteristic width $\widetilde{\Delta}$ superior to that of the implicit filter ($\overline{\Delta}$). If $a$ is a discrete variable defined over the computational domain, the variable obtained after applying the explicit filter to $a$ is noted $\tilde{a}$. Moreover, with: $$L_{ij} = \widetilde{\overline{u_i} \ \overline{u_j}} - \widetilde{\overline{u_i}}\ \widetilde{\overline{u_j}}$$ $$\tau_{ij} = \overline{u_iu_j} - \overline{u_i} \ \overline{u_j}$$ $$T_{ij} = \widetilde{\overline{u_iu_j}} - \widetilde{\overline{u_i}} \ \widetilde{\overline{u_j}}$$ Germano identity reads: $$L_{ij} = T_{ij} - \widetilde{\tau_{ij}}$$ Both dynamic models described herafter rely on the estimation of the tensors $\tens{T}$ and $\tens{\tau}$ as functions of the filter widths and of the strain rate tensor (Smagorinsky model). The following modelling is adopted\footnote{$\delta_{ij}$ stands for the Kroeneker symbol.}: $$T_{ij}-\frac{1}{3}\delta_{ij} T_{kk}= -2 C \widetilde{\Delta}^2 |\widetilde{\overline{D_{ij}}}| \widetilde{\overline{D_{ij}}} $$ $$\tau_{ij}-\frac{1}{3}\delta_{ij} \tau_{kk}= -2 C^* \overline{\Delta}^2 |\overline{D_{ij}}| \overline{D_{ij}} $$ $\overline{u}$ stands for the "implicit-filtered" value of a variable $u$ defined at the centres of the cells and $\overline{u}$ represents the "explicit-filtered" value associated with the variable $u$. It follows that the numerical computation of $L_{ij}$ is possible, since it requires the explicit filtering to be applied to implicitly filtered variables only ({\it i.e.} to the variables explicitly computed). On the following assumption: $$ C = C^* $$ and assuming that $C^*$ is only slightly non-uniform, so that it can be taken out of the explicit filtering operator, the following equation is obtained: $$L_{ij} -\frac{1}{3} \delta_{ij} L_{kk}= C (\alpha_{ij}-\widetilde{\beta}_{ij}) $$ with: $$\alpha_{ij} = -2 \widetilde{\Delta}^2 |\widetilde{\overline{D_{ij}}}| \widetilde{\overline{D_{ij}}} $$ $$ \beta_{ij} = -2 \overline{\Delta}^2 |\overline{D_{ij}}| \overline{D_{ij}}$$ Since we are left with six equations to determine one single variable, the least square method is used\footnote{$L_{kk}$ has no effect for incompressible flows.}. With: $$E_{ij} = L_{ij}-C(\alpha_{ij} - \widetilde{\beta}_{ij}) $$ the value for $C$ is obtained by solving the following equation: $$ \frac {\partial E_{ij}E_{ij}}{\partial C} = 0$$ Finally: $$C = \frac{M_{ij}L_{ij}} {M_{kl}M_{kl}} $$ with $$M_{ij} = \alpha_{ij} - \widetilde{\beta}_{ij}$$ This method allows to calculate the Smagorinsky "constant" dynamically at each time step and at each cell. However, the value obtained for $C$ can be subjected to strong variations. Hence, this approach is often restricted to flows presenting one or more homogeneous directions (Homogeneous Isotropic Turbulence, 2D flows presenting an homogeneous spanwise direction...): indeed, in such cases, the model can be (and is) stabilized by replacing $C$ by an average value of $C$ computed over the homogeneous direction(s). For a general case (without any homogeneous direction), a specific averaging is introduced to stabilize the model: for any given cell of the mesh, the averaged Smagorinsky constant is obtained as an average of $C$ over the "extended neighbourhood" of the cell (the set of cells that share at least one vertex with the cell considered). More precisely, the average value (also denoted $C$ hereafter) is calculated as indicated below: $$C = \frac{\widetilde{M_{ij}L_{ij}}} {\widetilde{M_{kl}M_{kl}}} $$ \clearpage {\bf Equations for the scalars} Two types of equations are considered: \\ {\tiny$\bigstar$} advection of a scalar with additional source terms: \begin{equation} \frac {\partial (\rho a)}{\partial t} + \underbrace{\,\dive\,((\rho \underline{u})\,a)}_{\text{convection}} - \underbrace{\,\dive\,(K \grad a)}_{\text{diffusion}} = TS_a +\Gamma\,a^i \end{equation} {\tiny$\bigstar$} advection of the variance $\widetilde{{a"}^2}$ with additional source terms ~: \begin{equation} \begin{array}{lcl} &\displaystyle \frac {\partial (\rho \widetilde{{a"}^2})}{\partial t} + \underbrace{\dive\,((\rho\,\underline{u})\ \widetilde{{a"}^2})}_{\text{convection}} - \underbrace{\dive\,(K\ \grad \widetilde{{a"}^2})}_{\text{diffusion}} = TS_{\widetilde{{a"}^2}} +\ \Gamma\,\widetilde{{a"}^2}^i\\ &\displaystyle \underbrace {+ 2\,\frac{\mu_t}{\sigma_t}(\grad \widetilde{a})^2 - \frac{\rho\,\varepsilon}{R_f k}\ \widetilde{{a"}^2}}_{\text{production and dissipation}} \end{array} \end{equation} The two previous equations can be unified formally as: \begin{equation}\label{depart} \frac {\partial (\rho f)}{\partial t} + \dive\,((\rho\,\underline{u}) f) - \dive\,(K \grad f) = TS_f + \Gamma\,f^i + T_s^{\,pd} \end{equation} with: \begin{equation} \begin{array}{lll} &\displaystyle T_s^{\,pd}= \begin{cases} 0 & \text{for $f=a$}, \\ 2\ \displaystyle \frac{\mu_t}{\sigma_t}(\grad \widetilde{a})^2 - \displaystyle \frac{\rho\,\varepsilon}{R_f k}\ \widetilde{{a"}^2} & \text{for $f=\widetilde{{a"}^2}$ } \end{cases} \end{array} \end{equation} $TS_f$ represents the additional source terms that may be prescribed by the user. \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Discretization in time} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% At first, the physical properties of the flow are computed (density, viscosity, specific heat...): indeed, they may depend upon the variables (such as the temperature for example). The time scheme is a $\theta$-scheme: \begin{equation} \left\{\begin{array}{ll} \theta = 1 & \text{for an implicit first order Euler scheme}\\ \theta = 1/2 & \text{for second order Crank-Nicolson scheme}\\ \end{array} \right. \end{equation} For the second order scheme, the time step is assumed to be constant. A fractional step scheme is used to solve the mass and momentum equations (Chorin). The first step (predictor step) provides predicted velocity components: they are determined sequentially and without coupling between each other. The mass equation is taken into account during the second step (corrector step): a pressure Poisson equation is solved and the mass fluxes at the cell faces are updated. If required, the equations for the turbulent variables are solved (turbulent kinetic energy and dissipation or Reynolds stresses and dissipation), using a $\theta$-scheme again. For the $k-\varepsilon$ model, an additional step is carried out to couple the source terms. For the Reynolds stress model, the variables (turbulent stresses and dissipation) are solved sequentially, without coupling. Next, the equations for the ``scalars'' (enthalpy, temperature, tracers, concentrations, mass fractions...) are solved, also with a $\theta$-scheme. To finish with, all the variables are updated and another time step may start. The general equation for advection (valid for the velocity components, the turbulent variables and the scalars) is re-written as follows in a condensed form; the mass equation ($\frac{\partial \rho } {\partial t}+ \dive(\rho \underline{u}) = \Gamma$) has been used to split the time derivative: \begin{equation}\label{simple} \rho \frac {\partial f}{\partial t} + \dive\,((\rho\,\underline{u}) f) - \dive\,(K \grad f) = S_{i}(\Phi,\varphi)\,f + S_{e}(\Phi,\varphi) + \dive\,(\rho\,\underline{u})\,f \end{equation} In this equation:\\ % \begin{tabular}{ll} $\Phi$ & : represents the physical properties $(\rho,K,\mu_{t},...)$\\ $\varphi$ & : represents the variables of the problem $(\vect{u},k,\epsilon,...)$\\ $S_{i}(\Phi,\varphi)\,f$ & : represents the linear part of the source terms\\ $S_{e}(\Phi,\varphi)$ & : includes all other source terms \\ $\dive\,(\rho\,\underline{u})\,f$ & : is the term associated with ``mass accumulation''\\ \end{tabular} \\\\ The time at which the different quantities are evaluated is indicated below:\\ $\bullet$ $\Phi$: the time considered is defined by the time scheme applied to the physical properties.\\ $\bullet$ $(\rho\,\underline{u})$: the time considered is defined by the time scheme applied to the mass flux.\\ $\bullet$ $S_{e}(\Phi,\varphi)$: the time considered is defined by the time scheme applied to the explicit source terms. If $\theta=1/2$, or if an extrapolation is used, the time step $\Delta t$ is constant in time and uniform in space. \subsection{Physical properties} The physical properties of the flow (density, viscosity, specific heat...) are: \begin{itemize} \item[-] explicit, defined at the time step $n$. \item[-] extrapolated at $n+\theta_{\Phi}$ using the Adam-Bashforth time scheme (in this case, the time step is assumed to be constant). \end{itemize} Under a more general form, this reads: \begin{equation} \Phi \equiv \Phi^{n+\theta_{\Phi}}=(1+\theta_{\Phi})\,\Phi^{n}- \theta_{\Phi}\,\Phi^{n-1} \end{equation} \begin{equation} \left\{\begin{array}{ll} \theta_{\Phi} = 0 & \text{standard explicit formulation}\\ \theta_{\Phi} = 1/2 & \text{second order extrapolation at } n+1/2\\ \theta_{\Phi} = 1 & \text{first order extrapolation at } n+1 \\ \end{array} \right. \end{equation} \subsection{Mass flux} For the mass flux, three time schemes are available. The mass flux may be: \begin{itemize} \item [-] explicit, taken at time step $n$ for the momentum equations and updated with its value at time step $n+1$ for the equations for turbulence and scalars (standard scheme).\\ \item [-] explicit, taken at time step $n$ for the momentum equations and also for the equations for turbulence and scalars. \item [-] taken at $n+\theta_{F}$ (second order if $\theta_{F}=1/2$). To solve the momentum equations, $(\rho\,\underline{u})^{n-2+\theta_{F}}$ and $(\rho\,\underline{u})^{n-1+\theta_{F}}$ are known. Hence, the value at $n+\theta_{F}$ is obtained as a result of the following extrapolation: \begin{equation} (\rho\,\underline{u})^{n+\theta_{F}}= 2\,\,(\rho\,\underline{u})^{n-1+\theta_{F}} -\,\,(\rho\,\underline{u})^{n-2+\theta_{F}} \end{equation} At the end of this phase (after the pressure correction step), $(\rho\,\underline{u})^{n+1}$ is known and the following interpolation is used to determine the mass flux at $n+\theta_{F}$ that will be adopted for the equations for turbulence and scalars: \begin{equation} (\rho\,\underline{u})^{n+\theta_{F}}= \frac{1}{2-\theta_{F}}\,(\rho\,\underline{u})^{n+1} +\frac{1-\theta_{F}}{2-\theta_{F}}\,(\rho\,\underline{u})^{n-1+\theta_{F}} \end{equation} \end{itemize} \subsection{Source terms} As for the physical properties, the explicit source terms are: \begin{itemize} \item[-] explicit: \begin{equation} [S_{e}(\Phi,\varphi)]^{n} = S_{e}(\Phi^{n+\theta_\Phi},\varphi^n)\\ \end{equation} \item[-] extrapolated at $n+\theta_{S}$ using the Adams-Bashforth scheme: \begin{equation} [S_{e}(\Phi,\varphi)]^{n+\theta_{S}} = (1+\theta_{S})\,S_{e}(\Phi^{n},\varphi^n)- \theta_{S}\,S_{e}(\Phi^{n-1},\varphi^{n-1})\\ \end{equation} \end{itemize} By default, to be consistent and preserve the order of convergence in time, the implicit source terms are discretized with the same scheme as that used for the unknown considered, {\it i.e.} taken at $n+\theta$: \begin{equation} [S_{i}(\Phi,\varphi)\,f]^{n+\theta} = S_{i}(\Phi^{n+\theta_\Phi},\varphi^n)\,[\theta\,f^{n+1}+(1-\theta)\,f^{n}]\\ \end{equation} \underline{Note:}\\ in the code, for $\theta_{S}=0$, the implicit source terms are taken at $n+1$ by convention (implicit Euler scheme, even if $\theta \ne 0$). \subsection{General discrete form} For the sake of clarity, it is assumed hereafter that, unless otherwise explicitly stated, the mass flux is taken at $n+\theta_F$ and the physical properties are taken at $n+\theta_\Phi$, with $\theta_F$ and $\theta_\Phi$ dependant upon the specific schemes selected for the mass flux and the physical properties respectively. Under a general form, the discrete counterpart of equation~(\ref{simple}) at $n+\theta$ reads: \begin{equation} \begin{array}{c} \displaystyle \frac{\rho}{\Delta t}(f^{n+1}-f^{n})+ \dive\,((\rho\,\underline{u}) f^{n+\theta})-\dive\,(K \grad f^{n+\theta}) = \\ \left[S_{i}(\Phi,\varphi)\,f\right]^{n+\theta} + \left[S_{e}(\Phi,\varphi)\right]^{n+\theta_{S}} + \dive\,(\rho\,\underline{u})\,f^{n+\theta}\\ \end{array} \end{equation} Using the standard $\theta$-scheme $f^{n+\theta}=\theta\,f^{n+1}+(1-\theta)\,f^{n}$, the equation reads: \begin{equation} \begin{array}{c} \displaystyle \frac{\rho}{\Delta t}(f^{n+1}-f^{n})+ \theta\,\dive\,((\rho\,\underline{u}) f^{n+1})-\theta\, \dive\,(K \grad f^{n+1}) = \\ - (1-\theta)\,\dive\,((\rho\,\underline{u}) f^{n})+(1-\theta)\,\dive\,(K \grad \,f^{n}) +\\ S_{i}(\Phi,\varphi^n)\left[\theta\,f^{n+1}+(1-\theta)\,f^{n}\right] + \left[S_{e}(\Phi,\varphi)\right]^{n+\theta_{S}} + \dive\,(\rho\,\underline{u})\,(\theta\,f^{n+1}+(1-\theta)\,f^{n}) \end{array} \end{equation} For numerical reasons, the system is solved in an iterative and incremental manner, with the help of the series $\delta f^{n+1,k+1}=f^{n+1,k+1}-f^{n+1,k}$ (with, by definition, $f^{n+1,0}=f^{n}$). In particular, this method allows to treat implicitly a portion of the advection terms associated with non orthogonal correction terms. Hence, the system actually solved reads: \begin{equation} \begin{array}{c} \displaystyle \underbrace{\left[\frac{\rho}{\Delta t}-\theta\,S_{i}(\Phi,\varphi^n)-\theta\,\dive\,(\rho\,\underline{u})\,\right]}_{ \displaystyle \var{ROVSDT}}\,(f^{n+1,k+1}-f^{n+1,k})\\ +\theta\,\dive\,((\rho\,\underline{u})\,(f^{n+1,k+1}-f^{n+1,k}))-\theta\, \dive\,(K \grad\,(f^{n+1,k+1}-f^{n+1,k})) =\\ \var{SMBRS} \left\{ \begin{array}{c} -\theta\,\dive\,((\rho\,\underline{u})\,f^{n+1,k})+\theta\, \dive\,(K \grad\,f^{n+1,k})\\ -(1-\theta)\,\dive\,((\rho\,\underline{u})\,f^{n})+(1-\theta)\, \dive\,(K \grad\,f^{n})\\ \displaystyle -\frac{\rho}{\Delta t}(f^{n+1,k}-f^{n})+S_{i}(\Phi,\varphi^n)\left[\theta\,f^{n+1,k}+(1-\theta)\,f^{n}\right]\\ +\dive\,(\rho\,\underline{u})\,(\theta\,f^{n+1,k}+(1-\theta)\,f^{n}) + \left[S_{e}(\Phi,\varphi)\right]^{n+\theta_{S}} \end{array} \right. \end{array} \end{equation} Whatever the equation considered (momentum equation, equations for turbulence or scalars,...) the system representation is always the same: a right-hand side (stored in the vector-array \var{SMBRS}) and a vector-array \var{ROVSDT} for the part linear with respect to $\delta f^{n+1,k+1}$.\\ The vector-array \var{ROVSDT} is computed by the subroutine \fort{covofi} for the scalars, by \fort{preduv} for the velocity and by \fort{turbke} or \fort{turrij} for the turbulence. The vector-array \var{SMBRS} is not computed at one go, but in two successive steps.\\ The first part is calculated by the subroutines \fort{covofi}, \fort{preduv}, \fort{turbke} and \fort{turrij} (respectively for the scalars, the velocity and the turbulence). At that point, the vector \var{SMBRS} is defined as follows: \begin{equation} \begin{array}{c} \var{SMBRS} = S_{i}(\Phi,\varphi^n)\,f^{n}+\dive\,(\rho\,\underline{u})\,f^{n}+\left[S_{e}(\Phi,\varphi)\right]^{n+\theta_{S}}\\ \end{array} \end{equation} then, the calculation of \var{SMBRS} is complemented at each sub-iteration during the resolution of the equation by the subroutine \fort{codits} as follows: \begin{equation} \begin{array}{c} \var{SMBRS} = \var{SMBRS} - \var{ROVSDT}\,(f^{n+1,k}-\,f^{n}) \\ -\theta\,\dive\,((\rho\,\underline{u})\,f^{n+1,k})+\theta\, \dive\,(K \grad\,f^{n+1,k})\\ -(1-\theta)\,\dive\,((\rho\,\underline{u})\,f^{n})+(1-\theta)\, \dive\,(K \grad\,f^{n})\\ \end{array} \end{equation} \newpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Discretization in space} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Within the framework of the finite volume approach, the equations are integrated over each cell of the mesh (or "control volume" $\Omega_i$). This section is limited to a brief description of the way the terms appearing in the equations are integrated. Specific attention is devoted to the calculation of gradients, since it is a major characteristic of the co-located finite volume method (all the variables are associated with the same point, namely the cell centre\footnote{The centre of a cell is a geometric point associated with the cell and located preferably inside the cell. Nevertheless, the word "centre" shall not be taken literally, %The use of the word "centre" is not absolutely correct, especially in the case of polyhedral cells that do not have a regular shape.}). The terms of {\bf order 0} ({\it i.e.} the terms that are not a space derivative) are integrated to introduce their average over the cell. For example, $\rho g$ becomes $|\Omega_i| \overline\rho_i\,g$. In this expression, $|\Omega_i|$ is the volume of $\Omega_i$ and $\overline\rho_i$ denotes the average of $\rho$ over the control volume (the cell) $\Omega_i$. For less clumsy notations, $|\Omega_i| \overline\rho_i\,g$ is simply written as $|\Omega_i| \rho_i\,g$. When ``reconstructions'' (in fact Taylor series) are required to reach a higher order in space, the average value $\rho_i$ is assumed to be associated with the centre I of $\Omega_i$. The "divergence" terms (or ``flux'' terms, or again ``conservative'' terms) are integrated using the Green relation to introduce cell faces values (and so, "fluxes" appear naturally). For example, a term such as $div(p)$ becomes\footnote{If the cell $i$ is at the domain boundary, the sum becomes $\sum\limits_{j\in Vois(i)}\overline p_{ij}S_{ij}+ \sum\limits_{k\in \gamma(i)}\overline p_{b_{ik}}S_{b_{ik}}$, with $b_{ik}$ referring to the faces of the cell $i$ which are at the domain boundary.} $\sum\limits_{j\in Vois(i)}\overline p_{ij}S_{ij}$. In this expression, $\overline p_{ij}$ is the average of $p$ on the interface between the neighbouring cells $i$ and $j$. As for the cell-average above, less clumsy notations are used: the sum is simply written $\sum\limits_{j\in Vois(i)} p_{ij}S_{ij}$. The value $p_{ij}$ is assumed to be associated with the centre F of the face $ij$ when ``reconstructions'' are needed to reach a higher order in space. The precision of the value $p_{ij}$ determines the precision of the calculation of $div(p)$. For $p_{ij}$, it is possible to take a non-centred and non-interpolated value (upwind scheme for convection) or a linear interpolation between the values at the centres I and J of the neighbouring cells. Both methods are relatively straighforward but may lack consistence and precision for arbitrary meshes (and in particular on non-orthogonal meshes). A higher order in space may be reached if reconstruction techniques are used. The idea is to compute the value for $p_{ij}$ more precisely: to do so, $p$ is interpolated at $F_{ij}$ (the centre of the face) using the values for $p$ at I and J and... the gradients of $p$ calculated at I and J. The reader will notice that it is precisely the calculation of the space derivatives\footnote{The first derivatives in space are required to obtain $div(p)$ in each cell.} of $p$ that motivated the need for a good interpolation of $p$ at $F_{ij}$. Hence, the computation of the space derivatives of $p$ requires to solve a system: this is done in an iterative manner. Doing so allows to calculate the {\bf "cell gradient"} of the variables. It is important to keep in mind that the "cell gradient" of a given variable represents the gradient of the variable in the cell and that it is obtained from the values of the variable interpolated at the cell faces. Similarly, the terms written as {\bf ``divergence of gradient''} are integrated to introduce face values. One has to calculate the gradient at each face (or "face gradient") in the direction normal to the face. This concept of "face gradient" is extremely important. The ``face gradient'' normal to a face can be easily calculated with just the values at the centres of the two cells that share the face (points I and J on figure~\ref{fig_geom}), but this method is limited to orthogonal meshes. % fa modification : reconstruction is compulsory for consistence For consistence and to reach a higher order in space, the values of the variables at points I' and J' have to be used. These are calculated by a Taylor series from the values at I and J and from the "cell gradient" previously determined. \begin{figure}[t] \parbox{8cm}{% \centerline{\includegraphics[height=4.5cm]{Images/facette.pdf}}} \parbox{8cm}{% \centerline{\includegraphics[height=4.5cm]{Images/facebord.pdf}}} \caption{\label{fig_geom}Identification of the geometric entities for internal (left) and boundary faces (right).} \end{figure} \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Calling tree} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Each sub-section of this document is associated with an important subroutine. The full list of the subroutines described here is the following: \fort{bilsc2} \fort{clptur} \fort{clsyvt} \fort{codits} \fort{condli} \fort{covofi} \fort{gradmc} \fort{gradrc} \fort{inimas} \fort{itrmas} \fort{matrix} \fort{navsto} \fort{preduv} \fort{recvmc} \fort{resolp} \fort{turbke} \fort{turrij} \fort{viscfa} \fort{visort} \fort{vissec}. The table~\ref{simple_calling_tree} presents their sequence within a time step. This calling tree is only partial. In particular, it does not account for the number of calls to each subroutine. Also, for the sake of clarity, no reference has been made to the subroutines dedicated to the gradient calculation (\fort{gradmc}, \fort{gradrc}), which are called very often. For the same reason, the calls to \fort{bilsc2} (advection fluxes) and \fort{matrix} (matrix calculation) which are made from within \fort{codits} (treatment of an advection equation with source terms) have not been reported. The sub-sections where important information can be found are indicated below:\\ \\ \nl \\ {\bf Calculation of gradients}\\ \hspace*{1cm}\fort{gradrc}\\ \hspace*{1cm}\fort{gradmc}\\ {\bf Least square method}\\ \hspace*{1cm}\fort{recvmc}\\ \hspace*{1cm}\fort{gradmc}\\ {\bf Convective schemes}\\ \hspace*{1cm}\fort{bilsc2}\\ {\bf Wall-laws} (for velocity and temperature)\\ \hspace*{1cm}\fort{clptur}\\ \hspace*{1cm}\fort{condli}\\ {\bf System solve} (incremental method)\\ \hspace*{1cm}\fort{codits}\\ {\bf Calculation of the values at the faces} (not exhaustive)\\ \hspace*{1cm}\fort{viscfa}\\ \hspace*{1cm}\fort{visort}\\ Finally, for the reader wishing to become more familiar with the methods implemented in \CS\, it is recommended to begin with the study of the advection equation for a scalar (\fort{covofi}) which is solved iteratively using an incremental method (\fort{codits}). It will then be useful to look at \fort{navsto} which briefly presents the solution of the system made up of the mass and momentum equations. \begin{table}[htp] \underline{\bf Calculation of the physical properties}\\ \underline{\bf Boundary Conditions}\\ \hspace*{1cm}\fort{condli} \\ \hspace*{1,5cm}\fort{clptur} \hspace*{1cm} ``turbulent'' conditions at the wall\\ \hspace*{1,5cm}\fort{clsyvt} \hspace*{1cm} symmetry conditions for the vectors and the tensors\\ \underline{\bf Navier-Stokes solution}\\ \hspace*{1cm}\fort{navsto}\\ \hspace*{1,5cm}{\bf Velocity prediction}\\ \hspace*{2,0cm}\fort{preduv} \\ \hspace*{2,5cm}\fort{vissec} \hspace*{1cm} momentum source terms related to the \\ \hspace*{4,5cm} \hspace*{1cm} transposed gradient of the velocity\\ \hspace*{2,5cm}\fort{viscfa} \hspace*{1cm} calculation of the viscosity at the faces\\ \hspace*{2,5cm}\fort{codits} \hspace*{1cm} iterative solution of the system using an incremental method\\ \hspace*{1,5cm}{\bf Pressure correction}\\ \hspace*{2,0cm}\fort{resolp} \\ \hspace*{2,5cm}\fort{viscfa} \hspace*{1cm}calculation of the time step at the faces...\\ \hspace*{2,5cm}\fort{visort} \hspace*{1,5cm} ...according to the selected options\\ \hspace*{2,5cm}\fort{matrix} \hspace*{1cm}calculation of the Poisson equation matrix\\ \hspace*{2,5cm}\fort{inimas} \hspace*{1cm}initialization of the mass flow rate\\ \hspace*{2,5cm}\fort{itrmas} \hspace*{1cm}update of the mass flow rate\\ \hspace*{1,5cm} {\bf Velocity correction}\\ \hspace*{3,2cm} \hspace*{1cm}standard method or...\\ \hspace*{2,0cm}\fort{recvmc} \hspace*{1cm}least square method\\ \underline{\bf $k-\varepsilon$ model}\\ \hspace*{1cm}\fort{turbke}\\ \hspace*{1,5cm}\fort{viscfa} \hspace*{1cm}preliminary steps before...\\ \hspace*{1,5cm}\fort{bilsc2} \hspace*{1,5cm}...source terms coupling\\ \hspace*{1,5cm}\fort{viscfa} \hspace*{1cm}calculation of the viscosity at the faces\\ \hspace*{1,5cm}\fort{codits} \hspace*{1cm}iterative solution of the systems using an incremental method\\ \underline{\bf Reynolds stress model}\\ \hspace*{1cm}\fort{turrij}\\ \hspace*{1,5cm}\fort{visort} \hspace*{1cm}calculation of the viscosity at the faces\\ \hspace*{1,5cm}\fort{codits} \hspace*{1cm}iterative solution of the systems using an incremental method\\ \underline{\bf Equations for the scalars}\\ \hspace*{1cm}\fort{covofi}\\ \hspace*{1,5cm}\fort{viscfa} \hspace*{1cm}calculation of the viscosity at the faces\\ \hspace*{1,5cm}\fort{codits} \hspace*{1cm}iterative solution of the systems using an incremental method\\ \caption{Partial and simplified calling tree associated with the successive stages within a time step.} \label{simple_calling_tree} \end{table}