Navier-Stokes solver (RANSOL module) ==================================== Governing Equations ^^^^^^^^^^^^^^^^^^^ The incompressible Navier-Stokes-Equations in a given three-dimensional domain :math:`\Omega` and time interval :math:`(0, t)` can be written as: .. math:: \begin{aligned} & \rho \left( \frac{\partial \boldsymbol{u}}{∂t} + \left(\boldsymbol{u}\cdot\nabla \right) \boldsymbol{u} \right) +\nabla p -\nabla \cdot \left(\mu\nabla \boldsymbol{u}\right) = \rho f & \text{in } \Omega \times (0,t) \\ &\nabla \boldsymbol{u} = 0 & \text{in } \Omega \times (0,t) \end{aligned} :label: where :math:`\boldsymbol{u} = \boldsymbol{u} (\boldsymbol{x}, t)` denotes the velocity vector, :math:`p = p (\boldsymbol{x}, t)` the pressure field, :math:`\rho ` the (constant) density, :math:`\mu` the dynamic viscosity of the fluid and :math:`f` the volumetric acceleration. The above equations need to be combined with the following boundary conditions: .. math:: \begin{aligned} & \boldsymbol{u} = \boldsymbol{u}_c & \text{in } \Gamma_D \times (0,t) \\ & p = p_c & \text{in } \Gamma_P \times (0,t) \\ & \boldsymbol{n} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{g}_1 = 0 & \text{in } \Gamma_M \times (0,T) \\ & \boldsymbol{n} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{g}_2 = 0 & \text{in } \Gamma_M \times (0,T) \\ & \boldsymbol{n} \cdot \boldsymbol{u} = \boldsymbol{u}_M & \text{in } \Gamma_M \times (0,T) \\ & \boldsymbol{u}(\boldsymbol{x}, 0) = \boldsymbol{u}_o(\boldsymbol{x}) & \text{in } \Omega_D \times (0) \\ & p(\boldsymbol{x}, 0) = p_o(\boldsymbol{x}) & \text{in } \Omega_D \times (0) \end{aligned} :label: In the above equations, :math:`\Gamma ≔∂Ω`, denotes the boundary of the domain :math:`\Omega`, with :math:`\boldsymbol{n}` the normal unit vector, and :math:`\boldsymbol{g}_1`, :math:`\boldsymbol{g}_2` the tangent vectors of the boundary surface :math:`\partial \Omega`. :math:`\boldsymbol{u}_c` is the velocity field on :math:`\Gamma_D` (the part of the boundary of Dirichlet type, or prescribed velocity type), :math:`p_c` the prescribed pressure on :math:`\Gamma_P` (prescribed pressure boundary). :math:`\sigma` is the stress field, :math:`\boldsymbol{u}_M` the value of the normal velocity and :math:`\boldsymbol{u}_0`, :math:`p_o` the initial velocity and pressure fields. The union of :math:`\Gamma_D`, :math:`\Gamma_P` and :math:`\Gamma_M` must be :math:`\Gamma`; their intersection must be empty, as a point of the boundary can only be part of one of the boundary types, unless it is part of the border between two of them. The spatial discretization of the Navier-Stokes equations has been done by means of the finite element method, while for the time discretization an iterative algorithm that can be considered as an implicit two steps "Fractional Step Method" has been used. Problems with dominating convection are stabilized by the so-called "Finite Increment Calculus" method, presented below. CFD algorithms and stability ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Using the standard Galerkin method to discretize the incompressible Navier-Stokes equations leads to numerical instabilities coming from two sources. Firstly, owing to the advective-diffusive character of the governing equations, oscillations appear in the solution at high Reynolds numbers, when the convection terms become dominant. Secondly, the mixed character of the equations limits the choice of velocity and pressure interpolations such that equal order interpolations cannot be used. In recent years, a lot of effort has been put into looking for ways to stabilize the governing equations, many of which involve artificially adding terms to the equations to balance the convection, for example, the artificial diffusion method [Oñate_2000]_. A new stabilization method, known as finite increment calculus, has recently been developed [Oñate_1999]_, [García_1998]_. By considering the balance of flux over a finite sized domain, higher order terms naturally appear in the governing equations, which supply the necessary stability for a classical Galerkin finite element discretization to be used with equal order velocity and pressure interpolations. This method has been Christianized “Finite Calculus” (FIC). Finite calculus (FIC) formulation --------------------------------- To demonstrate this method, consider the flux problem associated with the incompressible conservation of mass in a 2D source less finite sized domain defined by four nodes, as shown in Figure 1. .. figure:: ../figures/tdyn/fig1.png :align: center Mass balance in 2 dimensions Consider the flux problem (mass balance) through the domain, taking the average of the nodal velocities for each surface. .. math:: \begin{aligned} & \frac{h_y}{2} \left( u ( x - h_x, y ) + u ( x - h_x, y - h_y ) - u ( x, y ) + u ( x, y - h_y ) \right) + \\ & + \frac{h_x}{2} \left( u ( x, y - h_y ) + u( x - h_x, y - h_y ) - u ( x - h_x, y ) + u( x, y) \right) = 0 \end{aligned} :label: Now expand these velocities using a Taylor series expansion, retaining up to second order terms. Write :math:`u (x,y) = u`: .. math:: \begin{aligned} & u(x-h_x,y) = u - h_x \frac{∂u}{∂x} + \frac{h_x^2}{2} \frac{\partial ^2 u}{∂x^2} - O(h_x^3) \end{aligned} :label: and similarly for the other two components of the velocity. Substituting this back into the original flux balance equation gives, after simplification, the following stabilized form of the 2D mass balance equation: .. math:: \begin{aligned} & \frac{∂u}{∂x} + \frac{∂v}{∂y} - \underline { \frac{h_x}{2} \left( \frac{\partial ^2 u}{∂x^2} - \frac{\partial ^2 v}{∂x∂y} \right) } - \underline { \frac{h_y}{2} \left( \frac{\partial ^2 u}{∂x∂y} - \frac{\partial ^2 v}{∂y^2} \right) } = 0 \end{aligned} :label: Note that as the domain size tends to zero, i.e. :math:`h_x , h_y \rightarrow 0`, the standard form of the incompressible mass conservation equation in 2D is recovered. The underlined terms in the equation provide the necessary stabilization to allow a standard Galerkin finite element discretization to be applied. They come from admitting that the standard form of the equations is an unreachable limit in the finite element discretization, i.e. by admitting that the element size cannot tend to zero, which is the basis for the finite element method. It also allows equal order interpolations of velocity and pressure to be used. Stabilized Navier-Stokes equations ---------------------------------- The Finite increment Calculus methodology presented above is used to formulate stabilized forms of the momentum balance and mass balance equations and the Neumann boundary conditions. The velocity and pressure fields of an incompressible fluid moving in a domain :math:`\Omega \in R^d (d = 2,3)` can be described by the incompressible Navier Stokes equations: .. math:: \begin{aligned} & \rho \frac{∂u_i}{∂t} + \rho \frac{\partial u_i u_j}{∂x_j} + \frac{∂p}{∂x_i} - \frac{∂s_{ij}}{∂x_j} = \rho f_i & \\ & \frac{∂u_i}{∂x_i} = 0 & i, j = 1, N \\ \end{aligned} :label: Where :math:`1 ≤ i`, :math:`j ≤ d`, :math:`\rho` is the fluid density field, :math:`u_i` is the ith component of the velocity field :math:`u` in the global reference system :math:`x_i`, :math:`p` is the pressure field and :math:`s_{ij}` is the viscous stress tensor defined by: .. math:: \begin{aligned} & s_{ij} = 2ν \left( ϵ_{ij} - \frac{1}{3} \frac{∂u_k}{∂x_k} \Delta_{ij} \right) & \epsilon_{ij} = \frac{1}{2} \left( \frac{∂u_i}{∂x_j} + \frac{∂u_j}{∂x_i} \right) \end{aligned} :label: The stabilized FIC form of the governing differential equation can be written as .. math:: \begin{aligned} & r_{m_i} - \frac{1}{2} h_{ij}^2 \frac{∂r_{m_i}}{∂x_j} = 0 & \text{in } \Omega \end{aligned} :label: .. math:: \begin{aligned} & r_d - \frac{1}{2} h_j^d \frac{∂r_d}{∂x_j} = 0 & \text{in } \Omega \end{aligned} :label: Summation convention for repeated indices in products and derivatives is used unless otherwise specified. In above equations, terms :math:`r_{m_i}` and :math:`r_d` denote the residual of Eq. (2-6) , and :math:`h_{ij}^m`, :math:`h_i^d` are the characteristic length distances, representing the dimensions of the finite domain where balance of mass and momentum is enforced. Details on obtaining the FIC stabilized equations and recommendation for the calculation of the stabilization terms can be found in [Oñate_1998]_. Let :math:`n` be the unit outward normal to the boundary :math:`\partial \Omega`, split into two sets of disjoint components :math:`\Gamma_t`, :math:`\Gamma_u` where the Neumann and Dirichlet boundary conditions for the velocity are prescribed respectively. The boundary conditions for the stabilized problem to be considered are [Oñate_2004]_: .. math:: \begin{aligned} & n_j \sigma_{ij} - t_i + \frac{1}{2} h_{ij}^m n_j r_{m_i} = 0 & \text{on } \Gamma_t \\ & u_j = u_j^p & \text{on } \Gamma_u \end{aligned} :label: Where :math:`t_i`, :math:`u_j^p` are the prescribed surface tractions and :math:`\sigma_{ij}` is the total stress tensor, defined as .. math:: \begin{aligned} & \sigma_{ij} = s_{ij} - p \Delta _{ij} \end{aligned} :label: Equations (2-8) - (2-10) constitute the starting point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations. An interesting feature of the FIC formulation is that it allows using equal order interpolation for the velocity and pressure variables [García_2003]_ and [Oñate_2004]_. Stabilized integral forms ------------------------- From the momentum equations, the following relations can be obtained ([Oñate_2004]_) .. math:: \begin{aligned} & \frac{∂r_d}{∂x_i} = \frac{h_{ii}^m}{2a_i} \frac{∂r_{m_i}}{∂x_j} \\ & a_i = \frac{2\mu}{3} + \frac{u_i h_i^d}{2} & \text{no sum in } i \end{aligned} :label: Substituting the equation above into Eq. (2-9) and retaining only the terms involving the derivatives of with respect to , leads to the following alternative expression for the stabilized mass balance equation .. math:: \begin{aligned} & r_d - \sum{i=1}^{n_d} τ_i \frac{∂r_{m_i}}{∂x_i} = 0 \\ & τ_i = \left( \frac{8\mu}{3h_{jj}^m h_j^d} + \frac{2\rho u_i}{h_{ii}^m} \right)^{-1} & \text{no sum in } i \end{aligned} :label: The :math:`\tau_i` ’s in Eq. (2-13) when multiplied by the density are equivalent to the intrinsic time parameters, seen extensively in the stabilization literature. The interest of Eq. (2-13) is that it introduces the first space derivatives of the momentum equations into the mass balance equation. These terms have intrinsic good stability properties as explained next. The weighted residual form of the momentum and mass balance equations is written as .. math:: \begin{aligned} & \int_{\Omega} ν_i \left[ r_{m_i} - \frac{1}{2} h_j \frac{∂r_{m_i}}{∂x_j} \right] d\Omega + \int_{\Gamma _i} ν_i \left[ n_j \sigma_{ij} - t_i + \frac{1}{2} h_j n_j r_{m_i} \right] d\Gamma = 0 \\ & \int_{\Omega} q \left[ r_d - \sum_{i=1}^N τ_i \frac{∂r_{m_i}}{∂x_i} \right] d\Omega = 0 \end{aligned} :label: where :math:`\nu_i`, :math:`q` are generic weighting functions. Integrating by parts the residual terms in the above equations leads to the following weighted residual form of the momentum and mass balance equations: .. math:: \begin{aligned} & \int_{\Omega} ν_i r_{m_i} d\Omega + \int_{\Gamma _i} ν_i (n_j \sigma_{ij} - t_i) d\Gamma + \int_{\Omega} \frac{1}{2} h_j \frac{∂ν_i}{∂x_j} r_{m_i} d\Omega = 0 \\ & \int_{\Omega} q r_d d\Omega + \int_{\Omega} \left[ \sum_{i=1}^N τ_i \frac{∂q}{∂x_i} r_{m_i} \right] d\Omega - \int_{\Gamma} \left[ \sum{i=1}^N q r_i n_i r_{m_i} \right] d\Gamma = 0 \end{aligned} :label: We will neglect here onwards the third integral in Eq. (2-15) by assuming that :math:`r_{m_i}` is negligible on the boundaries. The deviatoric stresses and the pressure terms in the second integral of Eq. (2-15) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are: .. math:: \begin{aligned} & \int_{\Omega} ν_i \rho \left( \frac{∂u_i}{∂t} + u_j \frac{∂u_i}{∂x_j} \right) + \frac{∂ν_i}{∂x_j} \left( \mu \frac{∂u_i}{∂x_j} - \Delta _{ij} p \right)d\Omega - \\ & - \int_{\Omega} ν_i \rho f_i d\Omega - \int_{\Gamma} ν_i t_i d\Gamma + \int_{\Omega} \frac{h_j}{2} \frac{∂ν_i}{∂x_j} r_{m_i} d\Omega = 0 \end{aligned} :label: and .. math:: \begin{aligned} & \int_{\Omega} q \frac{∂u_i}{∂x_i} d\Omega + \int_{\Omega} \left[ \sum_{i=1}^N τ_i \frac{∂q}{∂x_i} r_{m_i} \right] d\Omega = 0 \end{aligned} :label: In the derivation of the viscous term in Eq. (2-16) we have used the following identity holding for incompressible fluids (prior to the integration by parts) .. math:: \begin{aligned} & \frac{∂s_{ij}}{∂x_j} = 2\mu \frac{∂ϵ_{ij}}{∂x_j} = \mu \frac{\partial ^2 u_i}{∂x_j ∂x_j} \end{aligned} :label: Convective and pressure gradient projections -------------------------------------------- The computation of the residual terms are simplified if we introduce the convective and pressure gradient projections :math:`c_i` and :math:`\pi_i`, respectively defined as .. math:: \begin{aligned} & c_i = r_{m_i} - u_j \frac{∂u_i}{∂x_j} \\ & π_i = r_{m_i} - \frac{∂p}{∂x_i} \end{aligned} :label: We can express the residual terms in Eq. (2-16) and Eq. (2-17) in terms of :math:`c_i` and :math:`\pi_i`, respectively which then become additional variables. The system of integral equations is now augmented in the necessary number by imposing that the residuals vanish (in average sense). This gives the final system of governing equations as: .. math:: \begin{aligned} & \int_{\Omega} ν_i \rho \left( \frac{∂u_i}{∂t} + u_j \frac{∂u_i}{∂x_j} \right) + \frac{∂ν_i}{∂x_j} \left( \mu \frac{∂u_i}{∂x_j} - \Delta _{ij} p \right) d\Omega - \int_{\Omega} ν_i \rho f_i d\Omega - \\ & - \int_{\Gamma _t} ν_i t_i d\Gamma + \int_{\Omega} \frac{h_k}{2} \rho \frac{∂ν_i}{∂x_k} \left( u_j \frac{∂u_i}{∂x_j} + c_i \right) d\Omega = 0 \end{aligned} :label: .. math:: \begin{aligned} & \int_{\Omega} q \frac{∂u_i}{∂x_i} d\Omega + \int_{\Omega} \left[\sum_{i=1} ^N τ_i \frac{∂q}{∂x_i} \left( \frac{∂p}{∂x_i} + \pi_i \right) \right]d\Omega = 0 \end{aligned} :label: .. math:: \begin{aligned} & \int_{\Omega} b_i \left( u_j \frac{∂u_i}{∂x_j} + c_i \right) d\Omega = 0 & \text{no sum in } i \end{aligned} :label: .. math:: \begin{aligned} & \int_{\Omega} w_i \left( \frac{∂p_i}{∂x_j} + \pi_i \right)d\Omega = 0 & \text{no sum in } i \end{aligned} :label: Being :math:`i`, :math:`j`, :math:`k = 1,N` and :math:`b_i`, :math:`w_i` the appropriate weighting functions. Monolithic time integration scheme ---------------------------------- In this section an implicit monolithic time integration scheme, based on a predictor corrector scheme for the integration of Eq. (2-20) is presented. Let us first discretize in time the stabilized momentum Eq. (2-8), using the trapezoidal rule (or \theta method) as (see [Zienkiewicz_1995]_): .. math:: \begin{aligned} & \rho \left( \frac{u_i^{n+1} - u_i^n}{\Delta t} + \frac{\partial}{∂x_i} \left( u_i u_j \right)^{n + \theta} \right) + \frac{∂p^{n+1}}{∂x_i} - \frac{∂s_{ij}^{n+\theta}}{∂x_j} - \rho f_i^{n+\theta} - \\ & - \frac{1}{2} \rho h_{m_j} \frac{\partial}{∂x_j} \left( u_j^{n+\theta} \frac{∂u_i^{n+\theta}}{∂x_j} + c_i^{n+\theta} \right) = 0 \end{aligned} :label: where superscripts :math:`n` and :math:`\theta` refer to the time step and to the trapezoidal rule discretization parameter, respectively. For :math:`\theta = 1` the standard backward Euler scheme is obtained, which has a temporal error of :math:`O(t)`. The value :math:`\theta = 0.5` gives a standard Crank Nicholson scheme, which is second order accurate in time :math:`O(t^2)`. An implicit fractional step method can be simply derived by splitting Eq. (2-21) as follows: .. math:: \begin{aligned} & \rho \left( \frac{u_i^{*,n+1} - u_i^n}{\Delta t} + \frac{\partial}{∂x_i} \left( u_i u_j \right)^{n+\theta} \right) + \frac{∂p^{n+1}}{∂x_i} - \frac{∂s_{ij}^{n+\theta}}{∂x_j} - \rho f_i^{n+\theta} - \\ & - \frac{1}{2} \rho h_{mj} \frac{\partial}{∂x_j} \left( u_j^{n+\theta} \frac{∂u_i^{n+\theta}}{∂x_j} + c_i^{n+\theta} \right) = 0 \end{aligned} :label: .. math:: \begin{aligned} & u_i^{n+1} = u_i^{*,n+1} - \frac{\Delta t}{\rho } \left( \frac{∂p^{n+1}}{∂x_i} -\frac{∂p^n}{∂x_i} \right) \end{aligned} :label: On the other hand, substituting the last equation above into Eq. (2-9) and after some algebra leads to the following alternative mass balance equation .. math:: \begin{aligned} & \rho \frac{∂u_i^*}{∂x_i} - \Delta t \frac{\partial^2}{∂x_i ∂x_j} \left( p^{n+1} - p^n \right) + \sum_{i=1}^N τ_i \frac{\partial}{∂x_i} \left( \frac{∂p^{n+1}}{∂x_i} + \pi_i^{n+1} \right) = 0 \end{aligned} :label: The weighted residual form of the above equations can be written as follows: .. math:: \begin{aligned} & \int_{\Omega} ν_i \rho \left( \frac{u_i^{*,n+1} - u_i^n }{\Delta t} + u_j^{*,n+\theta} \frac{\partial}{∂x_j} u_i^{*,n+\theta} \right) d\Omega + \int_{\Omega} \rho \frac{∂ν_i}{∂x_j} \left( \frac{∂u_i^{*,n+\theta}}{∂x_j} - \Delta _{ij} p^n \right) d\Omega - \\ & - \int_{\Omega} ν_i \rho f_i^{n+\theta} d\Omega - \int_{\Gamma _t} ν_i t_i^{n+\theta} d\Gamma + \int_{\Omega} \frac{h_k}{2} \rho \frac{∂ν_i}{∂x_k} \left( u_j^{*,n+\theta} \frac{∂u_i^{*,n+\theta}}{∂x_j} + c_i^n \right) d\Omega = 0 \end{aligned} :label: .. math:: \begin{aligned} & \int_{\Omega} q \left( \rho \frac{∂u_i^{*}}{∂x_i} - \Delta t \frac{\partial ^2}{∂x_i ∂x_j} \left( p^{n+1} - p^n \right) \right) d\Omega + \int_{\Omega} τ_i \frac{∂q}{∂x_i} \left( \frac{∂p^{n+1}}{∂x_i} + \pi_i^{n+1} \right) d\Omega = 0 \end{aligned} :label: .. math:: \begin{aligned} & \int_{\Omega} ν_i \left( u_i^{n+1} - u_i^* + \frac{\Delta t}{\rho } \left( \frac{∂p^{n+1}}{∂x_i} - \frac{∂p^n}{∂x_i} \right) \right) d\Omega = 0 & \text{no sum in } i \end{aligned} :label: At this point, it is important to introduce the associated matrix structure corresponding to the variational discrete FEM form of Eq. (2-24): .. math:: \begin{aligned} & M \frac{1}{\Delta t} \overline{U}^{n+1} - M \frac{1}{\Delta t} U^n + K \left( \overline{U}^{n+\theta} \right) \overline{U}^{n+\theta} + S_1 \overline{U}^{n+\theta} + S_2 C - G P^n + M_{\Gamma _t} T = F \end{aligned} :label: .. math:: \begin{aligned} & \left( \Delta t L + L^{τ_2} \right) P^{n+1} + D^{τ_2} \Pi + G^T U^{n+1} = \Delta t L P^n \end{aligned} :label: .. math:: \begin{aligned} & M \frac{1}{\Delta t} U^{n+1} - M \frac{1}{\Delta t} \overline{U}^{n+1} - G P^{n+1} - G P^n = 0 \end{aligned} :label: .. math:: \begin{aligned} & N U^{n+1} + M C = 0 \end{aligned} :label: .. math:: \begin{aligned} & G P^{n+1} + M Π = 0 \end{aligned} :label: Where :math:`U`, :math:`P` are the vectors of the nodal velocity and pressure fields, :math:`T` is the vector of prescribed tractions and :math:`\Pi`, :math:`C` the vectors of convective and pressure gradient projections. Terms denoted by over-bar identify the intermediate velocity obtained from the fractional momentum equation. The system of equations above includes an error due to the splitting of the momentum equation. This error can be eliminated by considering the analogous system of equations .. math:: \begin{aligned} & M \frac{1}{\Delta t} U_{i+1}^{n+1} - M \frac{1}{\Delta t} U^n + K \left( U_i^{n+\theta} \right) U_{i+1}^{n+\theta} + S_1 U_{i+1}^{n+\theta} + S_2 C_i - G P_i^{n+1} + M_{\Gamma _t}T = F \end{aligned} :label: .. math:: \begin{aligned} & \Delta t L \left( P_{i+1}^{n+1} - P_i^{n+1} \right) + L^{τ_2} P_{i+1}^{n+1} + D^{τ_2} \Pi + G^T U_{i+1}^{n+1} = 0 \end{aligned} :label: .. math:: \begin{aligned} & N U_{i+1}^{n+1} + M C_{i+1} = 0 \end{aligned} :label: .. math:: \begin{aligned} & G P_{i+1}^{n+1} + M \Pi_{i+1} = 0 \end{aligned} :label: Where :math:`i` is the iteration counter of the monolithic scheme. Basically, in this final formulation, the convergence of the resulting monolithic uncoupled scheme is enforced by the first term of the second equation of Eq. (2-26) (see [Soto_2001]_). Compressible flows ------------------ In the previous sections the basics of the incompressible flow solver implemented in Tdyn has been introduced. The extension of that algorithm to compressible flows is straightforward. Eq. (2-23) can be easily modified to take into account the compressibility of the flow. The resulting equation is (see M. [Vázquez_1999]_) .. math:: \begin{aligned} & \rho \frac{∂u_i^*}{∂x_i} + \frac{α}{\Delta t} \left( p^{n+1} - p^n \right) - \Delta t \frac{\partial^2}{∂x_i ∂x_i} \left( p^{n+1} - p^n \right) + \sum_{i=1}^N τ_i \frac{\partial}{∂x_i} \left( \frac{∂p^{n+1}}{∂x_i} + \pi_i^{n+1} \right) = f \end{aligned} :label: In the above equation the terms :math:`α` and :math:`f` depend on the compressibility law used. Turbulence solvers ^^^^^^^^^^^^^^^^^^ At high Reynolds numbers the flow will certainly be turbulent, and the resulting fluctuations in velocities need to be taken into account in the calculations. A process known as Reynolds averaging [3] is applied to the governing equations whereby the velocities, :math:`u_i` are split into mean and a fluctuating component, where the fluctuating component, :math:`u_i’`, is defined by: .. math:: \begin{aligned} & u_i = \overline{u}_i + u_i^{'} \end{aligned} :label: In the expression above: .. math:: \begin{aligned} & \overline{u}_i \left( x,t \right) = \frac{1}{T} \int_{–T/2}^{T/2} u_i \left( x,t+τ \right) dτ \end{aligned} :label: This leads to extra terms in the governing equations that can be written as a function of ‘Reynolds stresses’ tensor, defined in Cartesian coordinates as: .. math:: \begin{aligned} & τ_{ij}^R = - \overline{\rho u_i^{'} u_j^{'}} \end{aligned} :label: Since the process of Reynolds averaging has introduced extra terms into the governing equations, we need extra information to solve the system of equations. To relate these terms to the other flow variables, a model of the turbulence is required. A large number of turbulence models exist of varying complexity, from the simple algebraic models to those based on two partial differential equations. The more complex models have increased accuracy at the expense of longer computational time. It is therefore important to use the simplest model that gives satisfactory results. Several turbulence models as Smagorinsky, :math:`k`, :math:`k-\epsilon`, :math:`k-\omega`, :math:`k-kt` and Spalart Allmaras have been implemented in Tdyn. The final form of the so-called Reynolds Averaged Navier Stokes Equations (RANSE) using these models is: .. math:: \begin{aligned} & r_{mi} - \frac{h_{mj}}{2} \frac{∂r_{mi}}{∂x_j} = 0 & \text{on } \Omega \text{ ,i,j = 1,2,3. No sum in} i \end{aligned} :label: .. math:: \begin{aligned} & r_d - \frac{h_{dj}}{2} \frac{∂r_d}{∂x_j} = 0 & \text{on } \Omega \text{ ,j = 1,2,3.} \end{aligned} :label: Where: .. math:: \begin{aligned} & r_{mi} = \rho \left[ \frac{∂u_i}{∂t} + \frac{\partial (u_i u_j )}{∂x_j} \right] + \frac{∂p}{∂x_i} - \frac{\partial}{∂x_j} \left[ (\mu + \mu_T) \frac{∂u_i}{∂x_j} \right] - \rho f_i \end{aligned} :label: .. math:: \begin{aligned} & r_d = \rho \frac{∂u_i}{∂x_i} \end{aligned} :label: being :math:`\mu T` the so-called eddy viscosity. The theory behind the above mentioned models may be found in the references, but a basic outline of the three most representative models is given below. Zero equation (algebraic) model: the Smagorinsky model ------------------------------------------------------ Zero equation models are the most basic turbulence models. They assume that the turbulence is generated and dissipated in the same place, and so neglect the diffusion and convection of the turbulence. They are based on the idea of turbulent viscosity, and prescribe a turbulence viscosity, :math:`\nu_t`, either by means of empirical equations or experiment. The Smagorinsky turbulence model is based on the idea of large eddy simulations (LES) in which the coherent large-scale structures are modelled directly within the computational mesh, whilst the small scales are modelled with the concept of eddy viscosity [Kleinstreuer_1997]_. The Reynolds stress is written in terms of the turbulence kinetic energy, :math:`k`, and eddy viscosity :math:`\nu_t` as: .. math:: \begin{aligned} & \overline{u_i^{'} u_j^{'}} = \frac{2}{3} k \Delta _{ij} - 2 \nu_t ϵ_{ij} \end{aligned} :label: where .. math:: \begin{aligned} & \epsilon = \frac{1}{2} \left( \underline{\nabla} \cdot \underline{u^{'}}+ \underline{\nabla} \cdot \underline{u^{'}} \right)^T \end{aligned} :label: Smagorinsky proposed that the eddy viscosity depends on the mesh density and velocity gradients, and suggested the following expression: .. math:: \begin{aligned} & \nu_t = \frac{\mu_t}{\rho } = C {h^e}^2 \sqrt{\epsilon_{ij} \epsilon_{ij}} \end{aligned} :label: where :math:`\nu_T` is the kinematic eddy viscosity, :math:`C` is a constant of the order of :math:`0.001`, and :math:`h^e` is the element size. One-equation models: the kinetic energy model (:math:`k`-model) --------------------------------------------------------------- One-equation models attempt to model turbulent transport, by developing a differential equation for the transport of a turbulent quantity. In the kinetic energy model, the velocity scale of the turbulence is taken as the square root of the turbulence kinetic energy, :math:`k`. Kolmogorov and Prandtl independently suggested the following relationship between the eddy viscosity, this velocity scale, :math:`\sqrt{k}`, and the length scale, :math:`L`: .. math:: \begin{aligned} & \nu_t = C_D \sqrt{k L} \end{aligned} :label: where .. math:: \begin{aligned} & k = \frac{1}{2} \left( \overline{\underline{{u^{'}}^2}} = \frac{1}{2} \overline{\underline{{u^{'}}^2} + \underline{{v^{'}}^2} + \underline{{w^{'}}^2}} \right) \end{aligned} :label: and :math:`C_D` is an empirical constant, usually taken as :math:`1.0`. Manipulation of the Reynolds and momentum equations leads to a differential equation for :math:`k` [1]. The only gap in this analysis is that left by the length scale, :math:`L`. This has to be specified either from experiment, or empirical equation. In the analysis used here, it has been specified as :math:`1%` of the smallest mesh size. Clearly this is the largest downfall of the analysis, and many consider it insufficiently accurate [Hirsch_1991]_. However, it is more accurate than the zero-equation models and less computationally expensive than the more accurate 2-equation models, and so it can be a good compromise between computational expense and accuracy, depending on the problem. Two-equation models: the :math:`k-\epsilon` model ------------------------------------------------- In order to increase the accuracy of our turbulence modelling, it is therefore necessary to develop a second differential transport equation to provide a complete system of closed equations without the need for empirical relationships. The :math:`k-\epsilon` turbulence model develops two of such differential transport equations: one for the turbulence kinetic energy, :math:`k`, and a second one for the turbulent dissipation, :math:`\epsilon`. As for the :math:`k`-model, the :math:`k-\epsilon` model relies on the Prandtl-Kolmogorov expression for the eddy viscosity provided above. From dimensional arguments, the turbulent dissipation, can be written in terms of the turbulence kinetic energy and the turbulence length scale, and the eddy viscosity written in terms of k and \epsilon: .. math:: \begin{aligned} & \epsilon = \frac{C_ϵ k^{3/2}}{L} \end{aligned} :label: .. math:: \begin{aligned} & \nu_t = C_{\nu} \frac{k^2}{\epsilon} \end{aligned} :label: Where :math:`C_\epsilon` and :math:`C_{\nu}` are empirical constants, and are both taken to have a value of :math:`0.09`. In a similar way to that of the :math:`k`-model, differential transport equations can be formulated for :math:`k` and :math:`\epsilon` [Hirsch_1991]_, thus closing the system of equations, without having to empirically define any of the turbulent quantities. Since this method directly models the transport of all of the turbulent quantities, it is the most accurate. However, it involves solving 2 differential equations and is consequently the most computationally expensive. It is therefore necessary to evaluate the importance of turbulence in the problem before a model is selected. In many cases the easiest way to do this may be by trying several methods and using the simplest which gives the most accurate results. If turbulence is a relatively small feature of the flow then it is not necessary to waste computational expense by modelling it with a 2-equation model. At the same time, a zero-equation model would inadequately model a problem with large areas of turbulence, and in many cases this will lead to instabilities in the simulation. Implicit LES turbulence model ----------------------------- In recent years a significant progress has been carried out in the development of new turbulence models based on the fact that not the entire range of scales of the flow is interesting for the majority of engineering applications. In this type of applications information contained in "the large scales" of the flow is enough to analyze magnitudes of interest such as the velocity, temperature ... Therefore, the idea that the global flow behaviour can be correctly approximated without the necessity to approximate the smaller scales correctly, is seen by many authors as a possible great advance in the modelling of turbulence. This fact has originated the design of turbulence models that describe the interaction of small scales with large scales. These models are commonly known as Large Eddy Simulation models (LES). Numerous applications have shown that an extension of the FIC method allows to model low and high Reynolds number flows. In the standard large eddy simulation, a filtering process is applied to the Navier Stokes Eq. (2-6). After the filtering, a new set of equations are obtained, variables of this new set of equations are the filtered velocities or also called large scale velocities. As a consequence of the filtering process a new term called the subgrid scale tensor appears into the momentum equation. The subgrid scale tensor is defined in function of the large scale velocities and the subscale velocities too. Then it is necessary to model this term in function only of the large scale velocities. Several models exist for the subgrid scale tensor, but basically all of them propose an explicit description of the subgrid scale, an analytic expression in term of large scale velocities. In conclusion, all of these models define the subgrid scale tensor as a new nonlinear viscous term. It is our understanding that the stabilized FIC formulation, with an adequate evaluation of the characteristic lengths introduces the same effects of a LES model without giving an explicit expression of the subgrid scale tensor. All the effects related to turbulence modelling are included into the non-linear stabilization parameters. Computation of the characteristic lengths ----------------------------------------- The computation of the stabilization parameters is a crucial issue as they affect both the stability and accuracy of the numerical solution. The different procedures to compute the stabilization parameters are typically based on the study of simplified forms of the stabilized equations. Contributions to this topic are reported in [Oñate_1998]_ and [García_2005]_. Despite the relevance of the problem there still lacks a general method to compute the stabilization parameters for all the range of flow situations. .. figure:: ../figures/tdyn/fig2.png :align: center Decomposition of the velocity for characteristic length evaluation The application of the FIC/FEM formulation to convection-diffusion problems with sharp arbitrary gradients has shown that the stabilizing FIC terms can accurately capture the high gradient zones in the vicinity of the domain edges (boundary layers) as well as the sharp gradients appearing randomly in the interior of the domain [Oñate_1998]_. The FIC/FEM thus reproduces the best features of both, the so called transverse (cross-wind) dissipation or shock capturing methods. The approach proposed in this work is based on a standard decomposition of the characteristic lengths for every component of the momentum equation as .. math:: \begin{aligned} & h_i = \alpha \frac{\boldsymbol{u}}{|\boldsymbol{u}|} + β \frac{\nabla u_i}{|\nabla u_i|} \end{aligned} :label: Being :math:`\alpha` the projection of the characteristic length component in the direction of the velocity and :math:`\beta` the projection in the direction of the gradient of every velocity component. The decomposition defined by Eq. 40 can be demonstrated to be unique if the characteristic lengths are understood as tensors written in a system of coordinates aligned with the principal curvature directions of the solution. Obviously, Eq. (2-40) is a simplification of that approach. However, it can be easily demonstrated that the variational discrete FEM form of the resulting momentum equations are equivalent using one or the other approach. Then, the characteristic lengths are calculated as follows: .. math:: \begin{aligned} & \alpha = \left( coth⁡ \Gamma ^u - \frac{1}{\Gamma ^u} \right) u_j l_j \end{aligned} :label: .. math:: \begin{aligned} & \Gamma ^u = \frac{u_j l_j}{2\mu} \end{aligned} :label: Where :math:`l_j` are the maximum projections of the element on every coordinate axis as seen in the figure below: .. figure:: ../figures/tdyn/fig3.png :align: center Calculation of :math:`l_i` in a triangular element and .. math:: \begin{aligned} & \beta = \left( coth \gamma_{ij} - \frac{1}{\gamma _{ij}} \right) l^{\nabla u_i} \\ & \Gamma _{ij} = u_i \frac{\nabla u_i}{|\nabla u_i|} \frac{l^{\nabla u_i}}{2 \mu} \end{aligned} :label: .. math:: \begin{aligned} & l^{\nabla u_i} = \frac{\nabla u_i}{|\nabla u_i|} l_i & \text{, } j = 1,2 \end{aligned} :label: As for the length parameters :math:`h_i^d` in the mass conservation equation, the simplest assumption :math:`h_i^d = h^d` has been taken. The overall stabilization terms introduced by the FIC formulation presented above have the intrinsic capacity to ensure physically sound numerical solutions for a wide spectrum of Reynolds numbers without the need of introducing additional turbulence modelling terms. Flow through porous media ^^^^^^^^^^^^^^^^^^^^^^^^^ The incompressible Navier-Stokes quations, in a given domain :math:`\Omega` and time interval :math:`(0,t)` is given by Eq.(2-1). In the case of solids, small velocities can be considered and the term :math:`(u \cdot \nabla) u` can be neglected. Therefore, assuming and incompressible flow (constant density) in a certain domain :math:`\Omega` and considering the mass conservation equation, the Navier-Stokes equations can be written as follows: .. math:: \begin{aligned} & \rho \frac{∂u}{∂t} - \mu \nabla^2 u + \nabla p = f & \text{in } \Omega \times (0,t) \end{aligned} :label: .. math:: \begin{aligned} & \nabla \cdot u = 0 & \text{in } \Omega \times (0,t) \end{aligned} :label: The general form of the Navier-Stokes equation is valid for the flow inside pores of a porous media, but its solution cannot be generalized to describe the flow in porous region. Therefore, the general form of the Navier-Stokes equation must be modified to describe the flow through porous media. To this aim, Darcy's law is used to describe the linear relation between the velocity :math:`u` and the gradient of pressure :math:`p` within the porous media. It defines the permeability resistance of the flow in a porous media as: .. math:: \begin{aligned} & \nabla p = - \mu D u & \text{in } \Omega_p \times (0,t) \end{aligned} :label: where :math:`\Omega_p` is the porous domain, :math:`D` is the Darcy's law resistance matrix and :math:`u` the velocity vector. In the case of considering an homogeneous porous media, :math:`D` is a diagonal matrix with coefficients :math:`1/\alpha`, being :math:`\alpha` the permeability coefficient. As usual, the Reynolds number is defined as: .. math:: \begin{aligned} & Re = \frac{\rho U L}{\mu} \end{aligned} :label: being :math:`U` and :math:`L` a characteristic velocity and a characteristic length scale, respectively. In order to characterize the inertial effects, it is possible to define the Reynolds number associated to the pores as: .. math:: \begin{aligned} & Re_p = \frac{\rho U \delta}{\mu} \end{aligned} :label: where :math:`\delta` is the characteristic pore size. Whereas Darcy's law is reliable for values of :math:`Re_p < 1`, for larger values of :math:`Re_p` it is necessary to consider a more general model, which accounts also for the inertial effects, in the form: .. math:: \begin{aligned} & \nabla p = - \left( \mu D u + \frac{1}{2} \rho C u |u| \right) & \text{in } \Omega_p \times (0,t) \end{aligned} :label: where :math:`C` is the inertial resistance matrix. The modified Navier-Stokes equation in the whole domain results from considering the two source terms associated to the resistance induced by the porous medium (linear Darcy and inertial loss term). Hence, these source terms are added to the standard fluid flow momentum equation as follows: .. math:: \begin{aligned} & \rho \frac{∂u}{∂t} - \mu \nabla^2 u + \nabla p - \mu D u - \frac{1}{2} \rho C u |u| = 0 & \text{in } \Omega \times (0,t) \end{aligned} :label: Again, considering an homogeneous porous media, :math:`D` is a diagonal matrix with coefficients :math:`\frac{1}{\alpha}` and :math:`C` is also a diagonal matrix. Therefore, a modified Darcy's resistance matrix should be used in Tdyn as follows: .. math:: \begin{aligned} & \rho \frac{∂u}{∂t} - \mu \nabla^2 u + \nabla p = \mu D^* u & \text{in } \Omega \times (0,t) \end{aligned} :label: .. math:: \begin{aligned} & D^* = D + \frac{1}{2\mu} \rho C |u| I \end{aligned} :label: being :math:`I` the identity matrix. It should be noted that in laminar flows through porous media, the pressure :math:`p` is proportional to the velocity :math:`u`, and :math:`C` can be considered zero :math:`(D^* = D)`. Therefore, the Navier-Stokes momentum equation can be rewritten as: .. math:: \begin{aligned} & \rho \frac{∂u}{∂t} - \mu \nabla^2 u + \nabla p = - \mu D u & \text{in } \Omega \times (0,t) \end{aligned} :label: Note that this momentum equation is solved by Tdyn in the case of solids, where the :math:`(u \cdot \nabla) u` term vanishes due to small velocities. In the case that :math:`(u \cdot \nabla) u` cannot be neglected in the modelization (i.e. high velocities), then Tdyn should solve the followimg momentum equation within a fluid instead of a solid: .. math:: \begin{aligned} & \rho \left( \frac{∂u}{∂t} + \left(u \cdot \nabla \right) u \right) - \mu \nabla^2 u + \nabla p = - \mu D u & \text{in } \Omega \times (0,t) \end{aligned} :label: