Three-dimensional Navier-Stokes Calculations next up previous contents
Next: Convergence criteria Up: Types of analysis Previous: Turbulent Flow in Open   Contents

Three-dimensional Navier-Stokes Calculations

The solution of the three-dimensional Navier-Stokes equations has been implemented following the Characteristic Based Split (CBS) Method of Zienkiewicz and co-workers [69],[66]. The present implementation does give reasonable answers for a lot of problems (e.g. the lid-driven cavity), however, the accuracy is not sufficient for industrial use of solutions for compressible media. The problems are:

Further research is going on to locate the reason for these problems.

Fluid problems are of a quite different nature than structural problems. What we particularly noticed in fluid problems is that

The CBS Method transforms a transport equation of the form

$\displaystyle \frac{\partial x}{\partial t} = -(v_k C)_{,k} + D_{k,k}+F,$ (149)

where C stands for the convective term, D for the diffusion term and F for the external forces, into

$\displaystyle \sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta x_{\beta}=$ $\displaystyle - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k \varphi_{\beta})_{,k} C_{\beta} \right] dV$    
  $\displaystyle - \Delta t \int_V \varphi_{\alpha,k} D_k dV$    
  $\displaystyle + \Delta t \int_V \varphi_{\alpha} F dV$    
  $\displaystyle - \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} \left [ \sum_{\beta} (v_k \varphi_{\beta})_{,k} C_{\beta} \right] dV$    
  $\displaystyle + \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} F dV$    
  $\displaystyle + \Delta t \int_A \varphi_{\alpha} D_k n_k dA.$ (150)

Here $ \varphi_{\alpha}$ and $ \varphi_{\beta}$ are the shape functions. The first, second and third term on the right hand side correspond to convection, diffusion and external forces, respectively. The fourth and fifth terms are the stabilization terms for convection and external forces, while the last term is the area term corresponding to diffusion. It is the result of partial integration. The stabilization terms were obtained through partial integration too. In agreement with the CBS Method the corresponding area terms are neglected. Furthermore, third-order and higher order terms are neglected as well (particularly the stabilization terms corresponding to diffusion).

This method is now applied to the transport equations for mass, momentum and energy. Furthermore, the resulting momentum equation is split into two parts (Split scheme A in [69]), one part of which is calculated at the beginning of the iteration scheme. Subsequently, the conservation of mass equation is solved, followed by the second part of the momentum equation. To this end the correction to the momentum $ \Delta V_k=\rho \Delta v_k$ in direction $ k$ is written as a sum of two corrections:

$\displaystyle \Delta V_k = \Delta V_k^* + \Delta V_k^{**}.$ (151)

This results in the following steps:

Step 1: Conservation of Momentum (first part)

The partial differential equation reads:

$\displaystyle \frac{\partial V_i}{\partial t} = -\frac{\partial}{\partial x_k} ...{\partial t_{ik}}{\partial x_k} - \frac{\partial p}{\partial x_i} + \rho g_i.$ (152)

Applying the CBS method to all terms except the pressure term leads to:

$\displaystyle \sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta V_{\beta i}^*=$ $\displaystyle - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k \varphi_{\beta})_{,k} V_{\beta i} \right] dV$    
  $\displaystyle - \Delta t \int_V \varphi_{\alpha,k} (t_{ik}+t_{ik}^R) dV$    
  $\displaystyle + \Delta t \int_V \varphi_{\alpha} \rho g_i dV$    
  $\displaystyle - \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} \left [ \sum_{\beta} (v_k \varphi_{\beta})_{,k} V_{\beta i} \right] dV$    
  $\displaystyle + \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} \rho g_i dV$    
  $\displaystyle + \Delta t \int_A \varphi_{\alpha} (t_{ik}+t_{ik}^R) n_k dA.$ (153)

$ V_i$ is the momentum, $ t_{ik}$ is the diffusive stress and $ t_{ik}^R$ is the Reynolds stress multiplied by $ \rho$ (only for turbulent flow), all evaluated at time t. $ g_i$ is the gravity acceleration at time $ t+ \Delta t$. The diffusive stress satisfies

$\displaystyle t_{ik} = \mu (v_{i,k}+v_{k,i}-\frac{2}{3} v_{m,m} \delta_{ik})$ (154)

whereas $ t_{ik}^R$ is defined by

$\displaystyle t_{ik}^R = \mu_t (v_{i,k}+v_{k,i}-\frac{2}{3} v_{m,m} \delta_{ik})-\frac{2}{3} \rho k \delta_{ik}.$ (155)

Here, $ \mu_t$ is the turbulent viscosity and $ k$ is the turbulent kinetic energy. What is lacking in equation (153) to be equivalent to the momentum transport equation is the pressure term.

Step 2: Conservation of mass

The partial differential equation reads:

$\displaystyle \frac{\partial \rho}{\partial t} = - \frac{\partial V_i}{\partial x_i}$ (156)

Applying Galerkin, this leads to:

$\displaystyle \sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta \rho_{\beta}$ $\displaystyle + \theta_1 \theta_2 (\Delta t)^2 \sum_{\beta} \left[ \int_V \varphi_{\alpha,i} \varphi_{\beta,i} dV \right] \Delta p_{\beta}$    
  $\displaystyle = \Delta t \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta} V_{\beta i} \right] dV$    
  $\displaystyle + \theta_1 \Delta t \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta} \Delta V_{\beta i}^* \right] dV$    
  $\displaystyle -\theta_1 (\Delta t)^2 \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta,i} p_{\beta} \right] dV$    
  $\displaystyle -\Delta t \int_A \varphi_{\alpha} V_i n_i dA.$ (157)

In agreement with [66] the following approximation was made:

$\displaystyle \Delta t \int_A \varphi_{\alpha} [V_i+\theta_1(\Delta V_i^* + \Delta V_i^{**})] n_i dA \approx \Delta t \int_A \varphi_{\alpha} V_i n_i dA,$ (158)

leading to the last term in equation (157). The velocity in the mass conservation equation is calculated at time $ t+\theta_1 \Delta t$, whereas the pressure in the momentum transport equation is expressed at time $ t+\theta_2
\Delta t$ ( $ 0 \le \theta_1, \theta_2 \le 1$). If $ \theta_2 = 0$ the scheme is called explicit, else it is semi-implicit (in the latter case it is not fully implicit, since the diffusion term in the momentum equation is still expressed at time t). For compressible fluids (gas) an explicit scheme is taken. This means that the second term on the left hand side of equation (157) disappears and the only unknowns are $ \Delta \rho_{\beta}$. For incompressible fluids the density is constant and consequently the first term is zero: the unknowns are now the pressure terms $ \Delta p_{\beta}$.

An additional difference between compressible and incompressible fluids is that the left hand side of equation (157) for incompressible fluids (liquids) is usually not lumped: a regular sparse linear equation solver is used. For compressible fluids it is lumped, leading to a diagonal matrix. Lumping is also applied to all other equations (momentum,energy..), irrespective whether the fluid is a liquid or not.

Step 3: Conservation of Momentum (second part)

This equation takes care of the pressure term in the momentum equation, which was not covered by step 1:

$\displaystyle \sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta V_{\beta i}^{**}=$ $\displaystyle -\Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} \varphi_{\beta,i} p_{\beta} \right] dV$    
  $\displaystyle -\theta_2 \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} \varphi_{\beta,i} \Delta p_{\beta} \right] dV$    
  $\displaystyle -(1-\theta_2)\frac{\Delta t^2}{2} \int_V (\varphi_{\alpha} v_k)_{,k} \left[ \sum_{\beta} \varphi_{\beta,i} p_{\beta} \right] dV.$ (159)

Notice that for compressible fluids the second term on the right hand side disappears ( $ \theta_2 = 0$). Consequently, $ \Delta p$ is not needed for gases. This is good news, since only $ \Delta\rho$ is known at this point (conservation of mass).

Step 4: Conservation of Energy

The governing differential equations runs:

$\displaystyle \frac{\partial \rho \epsilon_t}{\partial t} = - [ v_k (\rho \epsi...
...p)]_{,k} + [t_{km}v_m + \kappa T_{,k}]_{,k} + [\rho f_k v_k + \rho h^{\theta}],$ (160)

where $ \epsilon_t$ is the total internal energy per unit of volume, $ \kappa$ is the conduction coefficient, $ f_k$ are the external forces and $ h^{\theta}$ represents volumetric heat sources. $ \epsilon_t$ satisfies

$\displaystyle \epsilon_t = \epsilon + c_v (v_i v_i)/2.$ (161)

Straightforward application of the CBS method yields

$\displaystyle \sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] (\Delta \rho \epsilon_t)_{\beta}=$ $\displaystyle - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k \varphi_{\beta})_{,k} (\rho \epsilon_t + p)_{\beta} \right] dV$    
  $\displaystyle - \Delta t \int_V \varphi_{\alpha,k} (t_{km}v_m + \kappa T_{,k}) dV$    
  $\displaystyle + \Delta t \int_V \varphi_{\alpha} [\rho f_k v_k + \rho h^{\theta}] dV$    
  $\displaystyle - \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} \le...
...sum_{\beta} (v_k \varphi_{\beta})_{,k} (\rho \epsilon_t + p)_{\beta} \right] dV$    
  $\displaystyle + \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l )_{,l} [\rho f_k v_k + \rho h^{\theta}] dV$    
  $\displaystyle + \Delta t \int_A \varphi_{\alpha} (t_{km}v_m + \kappa T_{,k}) n_k dA.$ (162)

For turbulent flows $ t_{km}$ has to be complemented by $ t_{km}^R$. For liquids the energy equation is uncoupled from the other equations, unless the temperature leads to motion due to differences in the density (buoyancy). For gases, however, there is a strong coupling with the other equations through the equation of state:

$\displaystyle p=\rho r T,$ (163)

where $ r$ is the specific gas constant.

Step 5: Turbulence

The turbulence implementation closely follows the equations in [42]. There are basically two extra variables: the turbulent kinetic energy $ k$ and the turbulence frequency $ \omega$. The governing differential equations read

$\displaystyle \frac{\partial \rho k}{\partial t} = - [ v_k (\rho k)]_{,k} + [(\mu+\sigma_k \rho \nu_t)k_{,k}]_{,k} + (t_{ij}^R u_{i,j} - \beta^* \rho \omega k)$ (164)


$\displaystyle \frac{\partial \rho \omega}{\partial t} = - [ v_k (\rho \omega)]_...
... \omega^2 + \frac{2}{\omega}(1-F_1) \rho \sigma_{\omega 2} k_{,j} \omega_{,j}).$ (165)

For the meaning of the constants the reader is referred to [42]. The turbulence equations are in a standard form clearly showing the convective, diffusive and external force terms. Consequently, application of the CBS scheme is straightforward.

Because of the problems occurring for laminar flow, the turbulence step has not been activated yet.

Notice that the unknowns in the systems of equations in all steps are the conservative variables $ V_i$, $ \rho$ (or $ p$ for liquids) and $ \rho \epsilon_t$. The physical variables the user usually knows and for which boundary conditions exist are $ v_i$, $ p$ and $ T$. So at the start of the calculation the initial physical values are converted into conservative variables, and within each iteration the newly calculated conservative variables are converted into physical ones, in order to be able to apply the boundary conditions.

The conversion of conservative variables into physical ones can be obtained using the following equations for gases:

$\displaystyle T=\frac{1}{\rho (c_p(T)-r)} \left[ (\rho \epsilon_t)-\frac{V_i V_i}{2 \rho} \right ] ,$ (166)

$\displaystyle v_i=V_i/\rho,$ (167)

and $ p=\rho r T$. For liquids $ \rho$ is a function of the temperature T and the first equation has to be replaced by

$\displaystyle T=\frac{1}{\rho(T) c_p(T)} \left[ (\rho(T) \epsilon_t)-\frac{V_i V_i}{2 \rho(T)} \right ] ,$ (168)

since $ c_v=c_p$. T in all equations above is the static temperature. For gases the total temperature and Mach number can be calculated by:

$\displaystyle T_t = T + V_i V_i /(2 c_p)$ (169)


$\displaystyle M=\sqrt{\frac{v_i v_i}{\gamma r T}}$ (170)

where $ \gamma = c_p/c_v$. Notice that the equations for the static temperature are nonlinear equations which have to be solved in an iterative way, e.g. by the Newton-Raphson procedure.

The semi-implicit procedure for fluids and the explicit procedure for liquids are conditionally stable. For each node $ i$ a maximum time increment $ \Delta
t_i$ can be determined. For the semi-implicit procedure it obeys:

$\displaystyle \Delta t_i = \min \left \lbrace \frac{h_i}{\sqrt{v_i v_i}} , \fra...
...i^2}{2 \mu(T_i)} , \frac {\rho_i h_i^2 Pr_i(T_i)}{2 \mu_i(T_i)} \right \rbrace,$ (171)


$\displaystyle Pr_i = \frac{\mu(T_i) c_p(T_i)}{\kappa(T_i)}$ (172)

is the Prandl number, and for the explicit procedure it reads

$\displaystyle \Delta t_i = \frac {h_i}{c_i+ \sqrt{v_i v_i}} ,$ (173)


$\displaystyle c_i = \sqrt \frac{c_p(T_i) r T_i}{c_p(T_i)-r}$ (174)

is the speed of sound. In the above equations $ h_i$ is the smallest distance from node i to all neighboring nodes. The overall value of $ \Delta t$ is the minimum of all nodal $ \Delta

Feasible elements are all linear volumetric elements (F3D4, F3D6 and F3D8). The Navier-Stokes procedure as not been tested yet for quadratic elements.

For gases a shock capturing technique has been implemented following [69]. This is essentially a smoothing procedure. To this end a field $ Sa_i$ is determined for each node i as follows:

$\displaystyle Sa_i=\frac{\vert \sum_i (p_i -p_j)\vert}{\sum_i \vert p_i - p_j \vert},$ (175)

where the sum is over all neighboring nodes and p is the static pressure. It can be verified that $ Sa_i=1$ for a local maximum and $ Sa_i=0$ if the pressure varies linearly. So $ Sa_i$ is a measure for discontinuous pressure changes. The smoothing procedure is such that the smoothed field $ \bar{x}$ is derived from the field $ x$ by

$\displaystyle {\bar{x}}_i = x_i + \frac{\Delta t C_e Sa_i}{\Delta t_i} [M_L]_{ii}^{-1} ([M]_{ij} - [M_L]_{ij}) x_j.$ (176)

$ [M]$ is the left hand side matrix for the variable involved, $ [M_L]$ is the lumped matrix (i.e. the matrix [M] where all values in each row are summed and put on the diagonal, all off-diagonal terms are zero) and $ C_e$ is a parameter between 0 and 2. The bigger $ C_e$, the stronger the smoothing. This procedure was elaborated in on [69]. After the regular calculation of $ \rho v_i$, $ \rho$ and $ \rho \epsilon_t$, the temperature $ T$ and the pressure $ p$ are calculated, the field $ Sa$ is determined and all conservative variables are smoothed. This leads to new values after which the boundary conditions for the velocity, the static pressure and static temperature are enforced again. If no convergence is reached, a new iteration is started.

It is important to note that for CFD calculations adiabatic boundary conditions have to be specified explicitly by using a *DFLUX card with zero heat flux. This is different from solid mechanics applications, where the absence of a *DFLUX or *DLOAD card automatically implies zero distributed heat flux and zero pressure, respectively.

Finally, it is worth noting that the construction of the right hand side of the systems of equations to solve has been parallelized (multithreading). Therefore you need the lpthread library at linking time. By setting the OMP_NUM_THREADS environment variable you can specify how many CPUs you would like to use (see Section 2).

next up previous contents
Next: Convergence criteria Up: Types of analysis Previous: Turbulent Flow in Open   Contents
guido dhondt 2012-10-06