CoreFlows/Documentation/numericalPage.md

   1 The numerical methods
   2 =====================
   3
   4 CoreFlows proposes a variety of finite volume methods (see \ref leveque for an introduction). The method  can be explicit or implicit (enum \ref TimeScheme), and the convection fluxes can be approximated using the  \ref roe,  \ref vfroe or \ref vffc formulations (enum \ref NonLinearFormulation).
   5
   6 The basic method for fluid models is the \ref roe scheme with entropic correction (see \ref kieu) and source upwinding (see\ref wbsourceterms). In order to increase precision it is possible to use centered, staggered, pressureCorrection or lowMach schemes through the enum \ref SpaceScheme.
   7
   8 The finite volume discretisation allows an easy handling of general geometries and meshes generated by \ref salome .
   9
  10 Explicit schemes are used in general for fast dynamics solved with small time steps while implicit schemes allow the use of large time step to quickly reach the stationary regime. The implicit schemes result in nonlinear systems that are solved using a Newton type method.
  11
  12 The upwind scheme is the basic scheme but options are available to use a centered scheme (second order in space) or entropic corrections.
  13
  14 Our models can be written in generic form as a nonlinear system of balance laws:
  15
  16 $$
  17 \frac{\partial U}{\partial t} + \nabla \cdot \left(\mathcal{F}^{conv}(U)\right)+\nabla \cdot \left(\mathcal{F}^{diff}(U)\right) = 0,
  18 $$
  19
  20 where
  21 - $U$ is the vector of conservative unknowns,
  22 - $\mathcal{F}^{conv}$ is the convective flux
  23 - and $\mathcal{F}^{diff}$ the diffusive flux.
  24
  25 We decompose the computational domain into $N$ disjoint cells $C_i$ with volume $v_i$.
  26 + Two neighboring cells $C_i$ and $C_j$ have a common boundary $\partial C_{ij}$ with area $s_{ij}$.
  27 + We denote $N(i)$ the set of neighbors of a given cell $C_i$ and $\vec n_{ij}$ the exterior unit normal vector of $\partial C_{ij}$ .
  28
  29 Integrating the system (\ref NStokesEq) over $C_{i}$ and setting $U_i(t)=  \frac{1}{v_i} \int_{C_i} U(x,t) dx$, the semi-discrete equations can be written:
  30
  31 $$
  32  \frac{\mathrm{d} U_i}{\mathrm{d} t} + \sum_{j \in N(i)} \frac{s_{ij}}{v_i}\left(\overrightarrow \Phi^{conv}_{ij} +  \overrightarrow{\Phi}^{diff}_{ij}\right)= S_i(U,x).
  33 $$
  34 with:
  35 - $\overrightarrow{\Phi}^{con}_{ij}= \frac{1}{s_{ij}}\int_{\partial C_{ij}}\mathcal F^{conv}(U).\vec n_{ij}ds $,
  36 - $\overrightarrow \Phi^{diff}_{ij}= \frac{1}{s_{ij}}\int_{\partial C_{ij}}\mathcal {F}^{diff}(U).\vec n_{ij}ds $
  37
  38 To approximate the convection numerical flux $\overrightarrow \Phi^{conv}_{ij}$ we solve an  approximate Riemann problem at the interface $\partial C_{ij}$. There are three possible formulations for the convection fluxes.
  39 - Using the \ref roe local linearisation of the fluxes, we obtain the following formula:
  40 $$
  41 \begin{array}{lll}
  42 \overrightarrow{\Phi}^{conv, Roe}_{ij}&= &\frac{\mathcal{F}^{conv}(U_i) + \mathcal{F}^{conv}(U_j)}{2} \vec{n}_{ij}- \mathcal{D}(U_i,U_j) \frac{U_j-U_i}{2}\\
  43 &=&\mathcal{F}^{conv}(U_i) \vec{n}_{ij} + A^-(U_i,U_j) (U_j - U_i),
  44 \end{array}
  45 $$
  46 - Using the \ref vfroe local linearisation of the fluxes, we obtain the following formula:
  47 $$
  48 \overrightarrow{\Phi}^{conv, VFRoe}_{ij}=\mathcal{F}^{conv}\left(\frac{U_i + U_j}{2} - \mathcal{D}(U_i,U_j) \frac{U_j-U_i}{2}\right)\vec{n}_{ij},
  49 $$
  50 - Using the \ref vffc local linearisation of the fluxes, we obtain the following formula:
  51 $$
  52 \overrightarrow{\Phi}^{conv, VFFC}_{ij}= \frac{\mathcal{F}^{conv}(U_i) + \mathcal{F}^{conv}(U_j)}{2} \vec{n}_{ij}- \mathcal{D}(U_i,U_j) \frac{\mathcal{F}^{conv}(U_j)-\mathcal{F}^{conv}(U_i)}{2} \vec{n}_{ij},
  53 $$
  54 where
  55 - $\mathcal{D}$ is an upwinding matrix,
  56 - $A(U_i,U_j)$ the Roe matrix
  57 - and $A^-= \frac{A - \mathcal{D}}{2}$.
  58
  59  The choice $\mathcal{D}= 0$ gives the \ref centered upwinding, $\mathcal{D}= |A|$ for the \ref roe formulation and  $\mathcal{D}= sign(A)$ for the \ref vfroe and \ref vffc formulations give the full \ref upwind upwinding. The \ref lowMach, \ref pressureCorrection and \ref staggered upwindings allow more precision for  Mach number flows.
  60
  61 The diffusion numerical flux $\overrightarrow\Phi_{ij}^{diff}$ is approximated on structured meshes using the formula:
  62 $$
  63 \overrightarrow \Phi_{ij}^{diff}= D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i),
  64 $$
  65
  66 where
  67 - $D(U,\vec{n}_{ij})=\nabla\mathcal{F}^{diff}(U).\vec{n}_{ij}$ is the matrix of the diffusion tensor.
  68 $(\ref{eq:flux diff})$ is not accurate for highly non structured or non conforming meshes. However, since we are mainly interested in convection driven flows, we do not ask for a very precise scheme.
  69
  70 Finally, since $\sum_{j \in N(i)}\mathcal {F}^{conv}(U_i). \vec{n}_{ij}=0$, using $(\ref{eq:flux roe})$ and $(\ref{eq:flux diff})$ the equation $(\ref{eq:numer scheme})$ of the semi-discrete scheme becomes:
  71 $$
  72 \frac{\mathrm{d} U_{i}}{\mathrm{d} t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i,U_j)\}(U_j-U_i)} = S_i(U,x),
  73 $$
  74
  75 The source term in $(\ref{eq:reduced scheme})$ can be approximated using either a
  76 $$
  77  \textrm{ Centered source } S_i=S(U_i)\nonumber
  78 $$
  79 or an
  80 $$
  81 \begin{array}{lll}
  82  \textrm{ Upwind source } S_i&=&\frac{1}{2}(Id-signe(A^{Roe}_{i,i+1}))\frac{S(U_i)+S(U_{i+1})}{2}\\
  83                               &&+\frac{1}{2}(Id+signe(A^{Roe}_{i-1,i}))\frac{S(U_{i-1})+S(U_i)}{2}.\nonumber
  84 \end{array}
  85 $$
  86
  87
  88
  89 Explicit schemes
  90 ----------------
  91
  92 In explicit schemes, in order to compute the values $U_i^{n+1}$, the fluxes $\Phi^{conv}_{ij}$, $\Phi^{diff}_{ij}$ and the source term $S(U,x)$ in $(\ref{eq:numer scheme})$ are evaluated at time $n$ :
  93 $$
  94 \begin{array}{lll}
  95 \frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} &+& \sum_{j\in N(i)} \frac{s_{ij}}{v_i}\left(\frac{1}{2}(\mathcal{F}^{conv}(U_i^n) + \mathcal{F}^{conv}(U_j^n)). \vec{n}_{ij}- \mathcal{D}(U_i^n,U_j^n,\vec{n}_{ij}) \frac{U_j^n-U_i^n}{2}\right)\\
  96 &+&\frac{s_{ij}}{v_i} D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i)= S(U^n,x_i), \nonumber
  97 \end{array}
  98 $$
  99 or equivalently using $(\ref{eq:flux roe})$ and $(\ref{eq:flux diff})$
 100 $$
 101 \frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i^{n},U_j^{n},\vec{n}_{ij})\}(U^{n}_j-U^{n}_i)} =  S(U^n,x_i).
 102 $$
 103
 104 From the system $(\ref{explicitscheme})$ we can obtain $U_i^{n+1}$ easily using matrix-vector products and vector sum.
 105 However the time step of explicit methods is constrained by the CFL condition for stability reasons.
 106
 107 Implicit schemes
 108 ----------------
 109
 110 In implicit schemes, in order to compute the values $U_i^{n+1}$, the fluxes $\Phi^{conv}_{ij}$, $\Phi^{diff}_{ij}$ and the source term $S(U,x)$ in $(\ref{eq:numer scheme})$ are evaluated at time $n+1$ :
 111 $$
 112 \begin{array}{lll}
 113 \frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} &+& \sum_{j\in N(i)} \frac{s_{ij}}{v_i}\left(\frac{1}{2}(\mathcal{F}^{conv}(U_i^{n+1}) + \mathcal{F}^{conv}(U_j^{n+1})). \vec{n}_{ij}- \mathcal{D}(U_i^{n+1},U_j^{n+1},\vec{n}_{ij}) \frac{U_j^{n+1}-U_i^{n+1}}{2}\right)\\
 114 &+&\frac{s_{ij}}{v_i} D (\frac{U_i+U_j}{2},\vec{n}_{ij})(U_j-U_i) = S(U^{n+1},x_i),
 115 \end{array}
 116 $$
 117 or equivalently using $(\ref{eq:flux roe})$ and $(\ref{eq:flux diff})$
 118 $$
 119 \frac{U_{i}^{n+1} - U_{i}^{n}}{\Delta t} + \sum_{j\in N(i)} {\frac{s_{ij}}{v_i}\{(A^-+ D)(U_i^{n+1},U_j^{n+1},\vec{n}_{ij})\}(U^{n+1}_j-U^{n+1}_i)} =  S(U^{n+1},x_i).
 120 $$
 121
 122 The system $(\ref{implicitscheme})$ is nonlinear. The computation of $U_i^n$ is more expensive but we can expect to use larger time steps than would be allowed with the explicit scheme.
 123
 124 We use the following Newton iterative method to obtain the required solutions:
 125 $$
 126 \begin{array}{lll}
 127 \frac{\delta U_i^{k+1}}{\Delta t}& + &\sum_{j \in N(i)} \frac{s_{ij}}{v_i} \left[( A^-+ D)(U_i^k,U_j^k) \right] \left(\delta U_j^{k+1}- \delta U_i^{k+1} \right)\\
 128 & = &- \frac{U^k_i-U^n_i}{\Delta t} - \sum_{j \in N(i)} \frac{s_{ij}}{v_i} \left[( A^-+ D)(U_i^k,U_j^k)\right] (U^k_j-U^k_i),
 129 \end{array}
 130 $$
 131 where :
 132 - $\delta U_i^{k+1} = U_i^{k+1} - U_i^{k}$ is the variation of the $k$-th iterate that approximate the solution at time $n+1$.
 133
 134 Defining the unknown vector $\mathcal U = (U_1,\dots,U_N)^t$, each Newton iteration for the computation of $\mathcal U$ at time step $n+1$ requires the numerical solution of the following linear system:
 135 $$
 136  \mathcal A(\mathcal U^k)\delta \mathcal U^{k+1} =  b(\mathcal U^n, \mathcal U^{k}).
 137 $$
 138
 139
 140
 141 Numerical scheme for the Navier-Stokes equations
 142 ------------------------------------------------
 143
 144 For the Navier-Stokes equations $U=(\rho, \vec q, \rho E)^t$ and the fluxes in $(\ref{Navierstokes})$ write
 145 $$
 146 \mathcal{F}^{conv}(U)= \left(\begin{array}{c}
 147 \vec{q} \\
 148 \vec{q} \otimes \frac{\vec{q}}{\rho} + p {I}_d \\
 149 \left( \rho E + p \right) \frac{\vec{q}}{\rho}
 150 \end{array}
 151  \right) ,\;
 152 \mathcal{F}^{diff}(U)=\left(\begin{array}{c}
 153 0\\
 154 -\nu \vec \nabla (\frac{\vec{q}}{\rho}) \\
 155 -\lambda \vec \nabla T
 156 \end{array}
 157  \right).
 158 $$
 159
 160 For the Euler equations, we can build the \ref roe matrix $A(U_i,U_j)$ explicitly using the Roe averaged state $U_{Roe}(U_i,U_j)=(\tilde\rho, \tilde \rho\tilde u, \tilde\rho \tilde E=\tilde \rho\tilde H-\tilde p)^t$ defined by
 161  $$
 162 \begin{array}{lll}
 163 \tilde{\rho}&=&\sqrt{\rho_{i}\rho_{j}}\\
 164 \tilde{u}&=&\frac{\sqrt{\rho_{i}}u_{i}+\sqrt{\rho_{j}}u_{j}}{\sqrt{\rho_{i}}+\sqrt{\rho_{j}}}\\
 165 \tilde{H}&=&\frac{\sqrt{\rho_{i}}H_{i}+\sqrt{\rho_{j}}H_{j}}{\sqrt{\rho_{i}}+\sqrt{\rho_{j}}}.
 166 \end{array}
 167 $$
 168 The Roe matrix writes (see \ref leveque)
 169 $$
 170 A_{Roe}(U_i,U_j)=\nabla\mathcal{F}^{conv}(U_{Roe}(U_i,U_j))\vec{n}_{ij}=
 171 \left(\begin{array}{ccc}
 172        0  & 1 & 0\\
 173        \tilde{\chi}+(\frac{1}{2}\tilde{\kappa} -1)\tilde{u}^2 & (2-\tilde{\kappa})\tilde{u} & \tilde{\kappa}\\
 174        (\tilde{\chi}+\frac{1}{2}\tilde{\kappa} \tilde{u}^2- \tilde{H})\tilde{u} & \tilde{H}-\tilde{\kappa} \tilde{u}^2 & (\tilde{\kappa}+1)\tilde{u}
 175       \end{array}\right)
 176 $$
 177
 178 The diffusion numerical flux $\overrightarrow\Phi_{ij}^{diff}$ is approximated with the formula:
 179 $$
 180 \overrightarrow \Phi_{ij}^{diff}= D (\frac{U_i+U_j}{2})(U_j-U_i)
 181 $$
 182  with the matrix
 183 $$
 184 D(U)=
 185 \left(
 186 \begin{array}{ccc}
 187  0&\vec 0& 0\\
 188 \frac{\nu \vec q}{\rho^2}& \frac{-\nu}{\rho} I_d&0\\
 189 \frac{\lambda}{c_v}\left(\frac{c_vT}{\rho}-\frac{||\vec q||^2}{2\rho^3}\right)&\quad \frac{{\vec q}^{\;t} \lambda}{\rho^2 c_v}&\quad-\frac{\lambda}{c_v  \rho}
 190 \end{array}
 191 \right) , $$
 192 where $c_v$ is the heat capacity at constant volume.
 193