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 [99],[96].The present implementation does include laminar and turbulent calculations for compressible and incompressible fluids. The calculations are transient, however, they are pursued up to steady state or up to the number of iterations specified by the user.

The input deck format for CFD-calculations is very similar to structural calculations. Noticable differences are:

• boundary conditions are specified by the *BOUNDARY card. The velocity degrees of freedom are labeled 1 up to 3, the thermal degree is 11 and the pressure degree is 8.
• the maximum number of iterations is specified by the INCF parameter on the *STEP card. The writing frequency on e.g. the *NODE FILE card is specified by the FREQUENCYF parameter.

For incompressible flows the following additional comments are due:

• thermal calculations do not influence the velocity and the pressure field. A calculation is considered thermal if initial conditions have been specified for the temperature.
• the material properties are introduced by the *DENSITY and the *FLUID CONSTANTS card. In case temperatures are to be calculated the *CONDUCTIVITY card is needed.

For compressible flows the following additional information is needed:

• for compressible flow the temperature is strongly linked to the velocity and the pressure. Therefore, initial conditions for all these fields are needed.
• the *CONDUCTIVITY and *SPECIFIC GAS CONSTANT card underneath the *MATERIAL card are required, the *DENSITY card must not be used.
• the *PHYSICAL CONSTANTS card is required for the definition of absolute zero
• the *VALUES AT INFINITY card is needed for the calculation of $C_p$ and for the turbulence models.
• the *SHOCK SMOOTHING parameter on the *STEP card may be needed to obtain convergence (only for explicit calculations).
• the COMPRESSIBLE parameter on the *STATIC card is required.

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

• The solution can be mesh dependent, i.e. the fluid flow sometimes follows the element edges although this may be wrong. This is particularly true for coarse meshes.
• Coarse meshes may produce a solution which is completely wrong. Taking this solution as starting point for gradually finer meshes frequently leads to the right solution.
• If the boundaries of the mesh are too close to the area of interest the solution may not be unique. For instance, turbulent flow may lead to an undefined reentry at the exit of your mesh. Consequently, the boundaries of your mesh must be far enough away.

The basic idea of the CBS method is to formulate the governing equation in a coordinate system moving with the characteristics of the flow, leading to a disappearance of the convective first order terms. To illustrate this, we start from a one-dimensional equation in the non-conservative form (the velocity $v$ is brought outside the partial differentiation)

$$\frac{\partial \phi }{\partial t} + v \frac{\partial \phi }{\part... ...{\partial x} \left ( \kappa \frac{\partial \phi }{\partial x} \right ) - Q = 0,$$ (421)

exhibiting a transient, convective, diffusive and source term ($\phi$ is some dependent quantity such as temperature). Applying a change of variables from x to x':

$$d x' = d x - v d t,$$ (422)

where $x'$ moves with the fluid, this equation is transformed into:

$$\frac{\partial \phi }{\partial t} (x'(t),t) - \frac{\partial }{\partial x'} \left( \kappa \frac{\partial \phi }{\partial x'} \right ) - Q(x')=0,$$ (423)

i.e. the convective term disappears.

Figure 138: Formulating the equations along characteristics

Applying Finite Differences along the characteristic from time $t$ (superindex n) to time $t+ \Delta t$ (superindex n+1) leads to (Figure 138):

$$\frac{1}{\Delta t} \left ( \phi^{n+1} - \left . \phi^n \right \vert _{x-\delta} \right ) \approx \theta \left . \left [ \frac{\partial }{\partial x} \left ( \kappa \frac{\partial \phi }{\partial x} \right ) + Q \right ]^{n+1} \right \vert _x$$

$$+ (1-\theta) \left . \left [ \frac{\partial }{\partial x} \left ( \... ...\partial \phi }{\partial x} \right ) + Q \right ]^{n} \right \vert _{x-\delta},$$ (424)

where $\theta$ takes a value between 0 and 1. Now, by applying a Taylor series expansion the values at $x-\delta$ can be written as a function of values at $x$:

$$\left . \phi^n \right \vert _{x-\delta} = \left . \phi^n \right \... ...\left . \frac{\partial^2 \phi }{\partial x^2} ^n \right \vert _x + O(\delta^3),$$ (425)

$$\left . \frac{\partial }{\partial x } \left ( \kappa \frac{\parti... ...partial \phi }{\partial x} \right ) \right ] ^ n \right\vert _x + O(\delta^2 ),$$ (426)

and

$$\left . Q^n \right\vert _ {x-\delta} = \left . Q^n \right\vert _ ... ...\delta \left . \frac{\partial Q}{\partial x } ^n \right \vert _x + O(\delta^2).$$ (427)

Therefore, Equation (424) now yields (from now on the subindex $x$ is dropped to simplify the notation):

$$\frac{1}{\Delta t} ( \phi^{n+1} - \phi^n + \delta \frac{\partial \phi }{\partial x} ^n - \frac{\delta^2}{2} ... ...x} \left ( \kappa \frac{\partial \phi }{\partial x} \right ) + Q \right ]^{n+1}$$

$$+ (1-\theta) \left [ \frac{\partial }{\partial x } \left ( \kappa \... ...\left ( \kappa \frac{\partial \phi }{\partial x} \right ) \right ] ^ n \right ]$$

$$+ (1-\theta)\left [Q^n - \delta \frac{\partial Q}{\partial x } ^n \right ]$$ (428)

Now, $\delta = \overline{v} \Delta t$, where $\overline{v}$ can be approximated by:

$$\overline{v} \approx \frac{1}{2} \left ( v ^{n+1} + \left . v^n \right\vert _ {x-\delta} \right).$$ (429)

Since

$$\left . v^n \right\vert _{x-\delta} = v^n - \delta \frac{\partial v^n}{\partial x} + O(\delta^2)$$ (430)

one obtains:

$$\overline{v} = \frac{1}{2} (v^{n+1} + v^n) - \frac{\delta }{2} \frac{\partial v^n }{\partial x} + O(\Delta t^2)$$

$$= v^{n+1/2} - \left ( \frac{\overline{v} \Delta t }{2} \right ) \frac{\partial v^n}{\partial x} + O(\Delta t^2)$$

$$= v^{n+1/2} - \left ( \frac{v^{n+1/2} \Delta t }{2} \right ) \frac{\partial v^n}{\partial x} + O(\Delta t^2),$$ (431)

where

$$v^{n+1/2} := \frac{1}{2} (v^{n+1} + v^n)$$ (432)

was defined. Consequently:

$$\delta = v^{n+1/2} \Delta t - \left ( \frac{v^{n+1/2} \Delta t^2 }{2} \right ) \frac{\partial v^n}{\partial x} + O(\Delta t^3).$$ (433)

Substituting this in Equation (428) and setting $\theta=1/2$ leads to:

$$\frac{1}{\Delta t} ( \phi^{n+1} - \phi^n) = -\frac{1}{\Delta t} \left [ v^{n+1/2} \Delta t - \frac{v^{n+1/2}}... ...al v^n}{\partial x} + O(\Delta t^3) \right] \frac{\partial \phi }{\partial x}^n$$

$$+ \frac{1}{2 \Delta t} \left [ \left( v^{n+1/2} \right) ^2 \Delta t^2 + O(\Delta t ^3) \right] \frac{\partial^2 \phi }{\partial x^2 }^n$$

$$+ \frac{1}{2} \left [ \frac{\partial }{\partial x} \left ( \kappa \frac{\partial \phi }{\partial x} \right ) + Q \right ]^{n+1}$$

$$+ \frac{1}{2} \left [ \frac{\partial }{\partial x} \left ( \kappa \frac{\partial \phi }{\partial x} \right ) + Q \right ]^{n}$$

$$- \frac{\Delta t}{2} v^{n+1/2} \frac{\partial }{\partial x} \left [... ... \kappa \frac{\partial \phi }{\partial x} \right ) \right ] ^ n + O(\Delta t^2)$$

$$- \frac{\Delta t}{2} v^{n+1/2} \frac{\partial Q}{\partial x } ^n + O(\Delta t^2).$$ (434)

Since

$$v^{n+1/2}=v^n + \frac{\partial v}{\partial t}^n \Delta t/2 + O(\Delta t^2),$$ (435)

$v^n$ in the first line of Equation (434) can be replaced by $v^{n+1/2}$ without loss of accuracy. Therefore, the terms quadratic in $\Delta t$ in the first two lines can be merged into:

$$v^{n+1/2} \frac{\Delta t^2}{2} \Delta t^2 \frac{\partial v^{n+1/2... ...ial }{\partial x } \left( v^{n+1/2} \frac{\partial \phi }{\partial x} \right ),$$ (436)

and one now obtains:

$$( \phi^{n+1} - \phi^n) = - \Delta t \left \lbrace v^{n+1/2} \frac{\partial \phi }{\partial... ... \frac{\partial \phi }{\partial x} \right ) + Q \right ]^{n+1/2} \right \rbrace$$

$$+ \frac{\Delta t ^2}{2} v^{n+1/2} \frac{\partial }{\partial x} \lef... ...ppa \frac{\partial \phi }{\partial x} \right ) - Q \right ] ^n + O(\Delta t^3).$$ (437)

In the last equation $v^{n+1/2}$ can be replaced by an extrapolation of $v$ at time $t_n + \Delta t/2$ based on its values in iteration $n-1$ and $n$ without loss of accuracy. Indeed, combining

$$v^{n+1/2} = v^n + \frac{\partial v}{\partial t} ^n \frac{ \Delta t_{n}}{2}+ O(\Delta t_{n}^2),$$ (438)

( $\Delta t_n := t_{n+1}-t_n$) and

$$v^{n-1} = v^n - \frac{\partial v}{\partial t} ^n \Delta t_{n-1} + O(\Delta t_{n-1}^2),$$ (439)

( $\Delta t_{n-1} := t_{n}-t_{n-1}$) or, equivalently,

$$\frac{\partial v}{\partial t}^n = \frac{v_n - v_{n-1}}{\Delta t_{n-1}} + O(\Delta t_{n-1}),$$ (440)

one obtains

$$v^{n+1/2} = v^n + \frac{v^n - v^{n-1}}{\Delta t_{n-1}} \left( \frac{\Delta t_n}{2} \right ) + O(\Delta t_n \Delta t_{n-1}) + O (\Delta t_n^2),$$ (441)

or for $\Delta t_{n-1}=\Delta t_n= \Delta t$:

$$v^{n+1/2} = v^n + \frac{v^n - v^{n-1}}{2} + O (\Delta t^2).$$ (442)

In the same way the diffusive and source terms at time $t_{n+1/2}$ are evaluated based on a similar extrapolation of the velocity and temperature (for the momentum and energy equation, respectively).

Generalizing Equation (437) to three dimensions and writing the equation in conservative form (i.e. replacing $v^{n+1/2} \partial \phi / \partial x$ by $\partial v^{n+1/2} \phi / \partial x$) finally yields:

$$( \phi^{n+1} - \phi^n) \approx - \Delta t \left \lbrace \frac{\partial (v_j^{n+1/2} \phi^n) }{\p... ...frac{\partial \phi }{\partial x_j} \right ) + Q \right ]^{n+1/2} \right \rbrace$$

$$+ \frac{\Delta t ^2}{2} v_k^{n+1/2} \frac{\partial }{\partial x_k} ... ...} \left ( \kappa \frac{\partial \phi }{\partial x_j} \right ) - Q \right ] ^n .$$ (443)

The last three terms can be viewed as stabilization terms. Usually, only terms up to the second order derivative are taken into account. Therefore, the stabilization term for the diffusion is usually neglected.

The corresponding weak formulation is obtained by multiplying the above equation with the shape function $\varphi_\alpha$ for a concrete node and integrating over the volume. Therefore, the CBS Method transforms a transport equation of the form

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

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

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta C_{\beta}= - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} C_{\beta}^n \right] dV$$

$$- \Delta t \int_V \varphi_{\alpha,k} D_k^{n+1/2} dV$$

$$+ \Delta t \int_V \varphi_{\alpha} F^{n+1/2} dV$$

$$- \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_... ... \left [ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} C_{\beta}^n \right] dV$$

$$+ \Delta t \int_A \varphi_{\alpha} D_k^{n+1/2} n_k dA$$

$$+ \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_{,l} F^n dV.$$ (445)

Notice that the integral over the total volume in reality is a sum of the integrals over each element. For each element the local shape functions are used in expressions such as $C=\sum_{\beta } \varphi_{\beta} C_{\beta}$.

The first, second and third term on the right hand side correspond to convection, diffusion and external forces, respectively. The fourth and sixth terms are the stabilization terms for convection and external forces, while the fifth 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 [99]), 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:

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

This results in the following steps:

Step 1: Conservation of Momentum (first part)

The partial differential equation reads:

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

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

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta V_{\beta i}^*= - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} V_{\beta i}^n \right] dV$$

$$- \Delta t \int_V \varphi_{\alpha,k} (t_{ik}+t_{ik}^R)^{n+1/2} dV$$

$$+ \Delta t \int_V \varphi_{\alpha} \rho g_i^{n+1/2} dV$$

$$- \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_... ...left [ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} V_{\beta i}^n \right] dV$$

$$+ \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_{,l} \rho g_i^n dV$$

$$+ \Delta t \int_A \varphi_{\alpha} (t_{ik}+t_{ik}^R)^{n+1/2} n_k dA.$$ (448)

$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

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

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

$$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}.$$ (450)

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

Step 2: Conservation of mass

The partial differential equation reads:

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

This can be approximated by:

$$\frac{\Delta \rho }{\Delta t} \approx - \frac{\partial V_i^n}{\pa... ...V_i^{*}}{\partial x_i}- \theta_1 \frac{\partial \Delta V_i^{**}}{\partial x_i},$$ (452)

where $\theta_1$ is a parameter leading to an explicit scheme for $\theta_1=0$ and an implicit scheme for $\theta_1=1$. Now, for $\Delta V_i^{**}$ one can use the gradient of the pressure in the momentum equation (this term can be treated in a way similar to a source term):

$$\Delta V_i^{**} \approx - \Delta t \left( \frac{\partial p}{\part... ...frac{\partial }{\partial x_k} \left( \frac{\partial p}{\partial x_i} \right)^n.$$ (453)

Before substituting Equation (453) into Equation (452) the stabilization term is dropped (leads to a third order derivative) and the pressure gradient at $n+1/2$ is changed into a gradient in between $n$ and $n+1$ by use of a parameter $\theta_2$ ($\theta_2$ is equivalent to $\theta$ in Equation(428):

$$\Delta V_i^{**} \approx - \Delta t \left( \frac{\partial p}{\partial x_i} \right ) ^{n} - \theta_2 \Delta t \frac{\partial \Delta p}{\partial x_i}.$$ (454)

For $\theta_2=0$ one obtains an explicit scheme (used for compressible media), for $\theta_2=1$ an implicit scheme (used for incompressible media). Now one obtains for Equation (452):

$$\frac{\Delta \rho }{\Delta t} \approx - \frac{\partial V_i^n}{\pa... ..._i}^n + \theta_2 \frac{\partial^2 \Delta p}{\partial x_i \partial x_i}\right ].$$ (455)

Applying Galerkin and partial integration to all terms on the right, this leads to:

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta \rho_{\beta} + \theta_1 \theta_2 (\Delta t)^2 \sum_{\beta} \left[ \int_V \varphi_{\alpha,i} \varphi_{\beta,i} dV \right] \Delta p_{\beta}$$

$$= \Delta t \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta} V_{\beta i}^n \right] dV$$

$$+ \theta_1 \Delta t \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta} \Delta V_{\beta i}^* \right] dV$$

$$-\theta_1 (\Delta t)^2 \int_V \varphi_{\alpha,i} \left[ \sum_{\beta} \varphi_{\beta,i} p_{\beta}^n \right] dV$$

$$-\Delta t \int_A \varphi_{\alpha} V_i^n n_i dA.$$ (456)

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

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

leading to the last term in equation (456). 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 (456) 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 (456) 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. Now, the terms are evaluated at $n+\theta_2$:

$$\Delta V_i^{**} \approx - \Delta t \left( \frac{\partial p}{\part... ...frac{\partial }{\partial x_k} \left( \frac{\partial p}{\partial x_i} \right)^n.$$ (458)

In weak form this leads to (applying partial integration to the stabilization term):

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] \Delta V_{\beta i}^{**}= -\Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} \varphi_{\beta,i} p_{\beta} \right] dV$$

$$-\theta_2 \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} \varphi_{\beta,i} \Delta p_{\beta} \right] dV$$

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

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 equation runs:

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

where $\varepsilon_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. $\varepsilon_t$ satisfies

$$\varepsilon_t = \varepsilon + (v_i v_i)/2.$$ (461)

The energy equation in the above form can be directly obtained from Equation (28). Indeed, the right hand side in both equations is identical. The left hand side of Equation (28) can be written as:

$$\rho \frac{D\varepsilon_t }{Dt} = \frac{D(\rho \varepsilon_t)}{Dt... ...frac{\partial \rho \varepsilon_t }{\partial t } + (\rho \varepsilon_t v_k)_{,k}$$ (462)

in which the conservation of mass was used in the form

$$\frac{D \rho }{D t} = \frac{\partial \rho }{\partial t} + \rho_{,k} v_k = - \rho v_{k,k}.$$ (463)

Straightforward application of the CBS method yields

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] (\Delta \rho \varepsilon_t)_{\beta}= - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho \varepsilon_t + p)_{\beta}^n \right] dV$$

$$- \Delta t \int_V \varphi_{\alpha,k} (t_{km}v_m + \kappa T_{,k})^{n+1/2} dV$$

$$+ \Delta t \int_V \varphi_{\alpha} [\rho f_k v_k + \rho h^{\theta}]^{n+1/2} dV$$

$$- \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_... ...v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho \varepsilon_t + p)_{\beta}^n \right] dV$$

$$+ \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_{,l} [\rho f_k v_k + \rho h^{\theta}]^n dV$$

$$+ \Delta t \int_A \varphi_{\alpha} (t_{km}v_m + \kappa T_{,k})^{n+1/2} n_k dA.$$ (464)

If a heat flux boundary condition $\boldsymbol{q}$ is specified the term $\kappa \boldsymbol{\nabla}T$ is replaced by $\boldsymbol{q}$. Furthermore, for turbulent flows $t_{km}$ is replaced by $t_{km}+t_{km}^R$ and $\kappa$ by $\kappa + c_p \mu_t/$Pr$_t$, where Pr$_t$ is the turbulent Prandl number (for air Pr$_t$=0.9). 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:

$$p=\rho r T,$$ (465)

where $r$ is the specific gas constant.

Step 5: Turbulence

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

$$\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)$$ (466)

and

$$\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}).$$ (467)

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

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] (\Delta \rho k)_{\beta}= - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho k)_{\beta}^n \right] dV$$

$$- \Delta t \int_V \varphi_{\alpha,k} (\mu+\sigma_k \rho \nu_t)k_{,k}^{n+1/2} dV$$

$$+ \Delta t \int_V \varphi_{\alpha} [t_{ij}^R u_{i,j} - \beta^* \rho \omega k]^{n+1/2} dV$$

$$- \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_... ...[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho k)_{\beta}^n \right] dV$$

$$+ \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_{,l} [t_{ij}^R u_{i,j} - \beta^* \rho \omega k]^n dV$$

$$+ \Delta t \int_A \varphi_{\alpha} (\mu+\sigma_k \rho \nu_t)k_{,k}^{n+1/2} n_k dA.$$ (468)

$$\sum_{\beta} \left[ \int_V \varphi_{\alpha} \varphi_{\beta} dV \right] (\Delta \rho \omega )_{\beta}= - \Delta t \int_V \varphi_{\alpha} \left[ \sum_{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho \omega )_{\beta}^n \right] dV$$

$$- \Delta t \int_V \varphi_{\alpha,k} (\mu+\sigma_{\omega} \rho \nu_t)\omega_{,k}^{n+1/2} dV$$

$$+ \Delta t \int_V \varphi_{\alpha} [\frac{\gamma}{\nu_t} t_{ij}^R u_{i,j} - \beta \rho \omega^2$$

$$+ \frac{2}{\omega}(1-F_1) \rho \sigma_{\omega 2} k_{,j} \omega_{,j}]^{n+1/2} dV$$

$$- \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_... ..._{\beta} (v_k^{n+1/2} \varphi_{\beta})_{,k} (\rho \omega )_{\beta}^n \right] dV$$

$$+ \frac{{\Delta t}^2}{2} \int_V ( \varphi_{\alpha} v_l^{n+1/2} )_{,l} [\frac{\gamma}{\nu_t} t_{ij}^R u_{i,j} - \beta \rho \omega^2$$

$$+ \frac{2}{\omega}(1-F_1) \rho \sigma_{\omega 2} k_{,j} \omega_{,j}]^n dV$$

$$+ \Delta t \int_A \varphi_{\alpha} (\mu+\sigma_{\omega} \rho \nu_t)\omega_{,k}^{n+1/2} n_k dA.$$ (469)

The above equations were slightly modified according to [88] in order to avoid a non-physical decay of the turbulence variables at the freestream boundary conditions. To this end the terms $-\beta^* \rho \omega k$ and $- \beta \rho \omega^2$ were replaced by $-\beta^* \rho \omega k + \beta^* \rho \omega_{\text{free}} k_{\text{free}}$ and $- \beta \rho \omega^2+\beta \rho \omega_{\text{free}}^2$, respectively, where the subscript “free” denotes the freestream values.

Turbulent equations require the definition of the nodal set SOLIDSURFACE containing all nodes belonging to a solid surface and the nodal set FREESTREAMSURFACE containing all nodes belonging to free stream conditions such as inlet and outlet. For each of these surfaces CalculiX assigns specific boundary conditions to the turbulence parameters $k$ and $\omega$ according to the publication by Menter.

Note that the conservative variables $\rho k$ and $\rho \omega$ should always be positive. Therefore, when the calculated change of these variables in each increment is added to their previous values only if the sum of both is positive, else the change is set to zero for that increment.

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 \varepsilon_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:

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

$$v_i=V_i/\rho,$$ (471)

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

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

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

$$T_t = T + V_i V_i /(2 c_p)$$ (473)

and

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

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:

$$\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,$$ (475)

where

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

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

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

where

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

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 t_i$'s.

Feasible elements are all linear volumetric elements (F3D4, F3D6 and F3D8).

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

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

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

$${\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.$$ (480)

$[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 on in [99]. After the regular calculation of $\rho v_i$, $\rho$ and $\rho \varepsilon_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).