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:

- The accuracy for compressible flow is not satisfactory. In particular, the variation of the total temperature in airfoil calculations is too strong.
- Without the smoothing effect of the shock capturing technique compressible flow problems do not converge.

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 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 CBS Method transforms a transport equation of the form

| (149)

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

|  | (150)

Here and 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 in direction is written as a sum of two corrections:

| (151)

This results in the following steps:

Step 1: Conservation of Momentum (first part)

The partial differential equation reads:

| (152)

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

|  | (153)

is the momentum, is the diffusive stress and is the Reynolds stress multiplied by (only for turbulent flow), all evaluated at time t. is the gravity acceleration at time . The diffusive stress satisfies

| (154)

whereas is defined by

| (155)

Here, is the turbulent viscosity and 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:

| (156)

Applying Galerkin, this leads to:

|  | (157)

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

| (158)

leading to the last term in equation (157). The velocity in the mass conservation equation is calculated at time , whereas the pressure in the momentum transport equation is expressed at time ( ). If 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 . For incompressible fluids the density is constant and consequently the first term is zero: the unknowns are now the pressure terms .

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:

|  | (159)

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

Step 4: Conservation of Energy

The governing differential equations runs:

| (160)

where is the total internal energy per unit of volume, is the conduction coefficient, are the external forces and represents volumetric heat sources. satisfies

| (161)

Straightforward application of the CBS method yields

|  | (162)

For turbulent flows has to be complemented by . 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:

| (163)

where 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 and the turbulence frequency . The governing differential equations read

| (164)

and

| (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 , (or for liquids) and . The physical variables the user usually knows and for which boundary conditions exist are , and . 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:

| (166)

| (167)

and . For liquids is a function of the temperature T and the first equation has to be replaced by

| (168)

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

| (169)

and

| (170)

where . 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 a maximum time increment can be determined. For the semi-implicit procedure it obeys:

| (171)

where

| (172)

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

| (173)

where

| (174)

is the speed of sound. In the above equations is the smallest distance from node i to all neighboring nodes. The overall value of is the minimum of all nodal 's.

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 is determined for each node i as follows:

| (175)

where the sum is over all neighboring nodes and p is the static pressure. It can be verified that for a local maximum and if the pressure varies linearly. So is a measure for discontinuous pressure changes. The smoothing procedure is such that the smoothed field is derived from the field by

| (176)

is the left hand side matrix for the variable involved, 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 is a parameter between 0 and 2. The bigger , the stronger the smoothing. This procedure was elaborated in on [69]. After the regular calculation of , and , the temperature and the pressure are calculated, the field 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: Convergence criteria Up: Types of analysis Previous: Turbulent Flow in Open Contents