Basic Principles

Fluid dynamics

To construct a model for fluid motion one looks at properties of matter in a control volume. These properties of a fluid are often expressed in terms of mass, momentum and energy. In most cases the fluid is considered to be a continuum, whereas for rarefied gases one needs to take into account the behaviour of molecules in a statistical way [15, p. 34].

To arrive at a mathematical model of fluid motion one looks at the properties of mass, momentum and energy in a control volume (Figure 2.1) and considers the following physical laws:

Figure 2.1: Fluid control volume

1. Conservation of mass,
2. Newton's second law of motion,
3. Conservation of energy, and
4. Second law of thermodynamics.

Conservation of mass

An expression for the conservation of mass is arrived at by considering the rate of change of mass in the control volume; that is, the flow of mass in and the flow of mass out of our control volume [15]. In differential form this is

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{V})=0$$ (1)

where is the density and is the velocity vector. Equation 2.1 is called the continuity equation.

If the flow is incompressible this becomes:

$$\nabla \cdot \mathbf{V}=0$$ (2)

Navier-Stokes equation

Equation 2.3, generally known as the Navier-Stokes equation, is obtained by considering the effect of Newton's second law of motion () on the control volume and by taking into account all the forces acting upon it [15, p. 54]. It is expressed here in cartesian tensor notation.

$$\rho \left( \frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial u_{i}}{\partial x_{j}}\right) = -\frac{\partial p}{\partial x_{i}}+B_{i}+\frac{\partial }{\partia... ...i}}-\frac{2}{3}\delta _{ij}\frac{\partial u_{k}}{\partial x_{k}}\right) \right]$$

$$+\frac{\partial }{\partial x_{i}}\left( \zeta \frac{\partial u_{k}}{\partial x_{k}}\right)$$ (3)

represents body forces, is the coefficient of viscosity and is the second coefficient of viscosity. is zero for a monatomic gas.

The Navier-Stokes equation is often simplified for various fluid dynamics problems by using the incompressible condition (Equation 2.2). This reduces the above equation to:

$$\rho \frac{D\mathbf{V}}{Dt}=-\nabla p+\mathbf{B}+\mu \nabla ^{2}\mathbf{V}$$

where is called the substantial or material derivative, and is given by:

$$\frac{D\mathbf{V}}{Dt}=\frac{\partial \mathbf{V}}{\partial t}+\left( \mathbf{V}\cdot \nabla \right) \mathbf{V}$$ (4)

Another way to simplify solutions to the Navier-Stokes equation is to keep constant.

Conservation of energy

The energy equation is derived by looking at the heat flux into and out of the control volume and by looking at the work done on the control volume by pressure and shear stresses [15, p. 56].

$$\rho C_{v}\frac{DT}{Dt}=-p\nabla \cdot \mathbf{V}+\kappa \nabla ^{2}T-\nabla \cdot \mathbf{q}_{r}+\Phi +q'''$$

Here, is the pressure, is the specific heat at constant volume, is the temperature, is the thermal conductivity, is the radiation heat flux vector, is the dissipation function and is the internal heat production rate per unit volume.

Equation of state

To fully track the state of a fluid, it is necessary to have an equation of state that relates fluid quantities to each other. For low and moderate density gases, the ideal gas equations (Equations 2.5, 2.6 and 2.7) can be used. In practice, these equations are also used in the solution of problems at higher densities, at the expense of complete accuracy.

The ideal gas equation of state is:

$$p=\rho RT$$ (5)

where is the gas constant. and are the specific heats at constant pressure and temperature respectively.

For an ideal gas the energy is related to the temperature by:

$$E=\rho C_{v}T$$ (6)

and pressure is related to density by:

$$\frac{p}{\rho ^{\gamma }}=constant$$ (7)

Current numerical methods

To exactly simulate fluid flow on a computer it would be necessary to solve the above equations with infinite accuracy. In reality, numerical researchers must choose a method to discretise the problem. Some of the general numerical methods used in computational fluid dynamics are described here.

The setting up of a numerical simulation begins with creating a computational grid. The flow variables are calculated at the node points of this grid and, in some methods, at some intermediate points as well. The spacing between grid points has to be fine enough to attain a high enough degree of accuracy. There are advantages, however, to keeping the number of grid points small since more grid points means more computer memory is required and a greater time is needed to perform each iteration of the calculation. Most computer architectures deal with smaller amounts of memory much more efficiently and cheaply than with large amounts of memory. The simplest computational grid is a rectangular lattice with fixed spacing between node points in each dimension. There are a range of methods that use unstructured grids where the density of the node points is not constant and is higher in the regions where more accuracy is required. Unstructured meshes often end up being connected in a triangular or tetrahedral fashion because these shapes fill space well and they require a minimum number of vertices. Some methods even use adaptive meshes where node points are created and destroyed as flow features move through the computational domain. This keeps the total number of nodes to a minimum, whilst still providing the required resolution for certain flow features.

A numerical simulation of a fluid flow problem requires that all flow variables that are tracked by the algorithm are given an initial value. In one scenario, just the free stream flow values are specified. Other scenarios need the specification of an `initial' condition which is specified at most points in the computational domain and an `inlet' condition which is specified only along an inlet boundary. Together with the initial conditions, boundary conditions define the nature of the flow problem being simulated. They include inlet conditions, various kinds of walls and outflow conditions. There are also special kinds of conditions that arise from the geometry of the problem like an axi-symmetric boundary.

Discretising the model equations

Finite differences are discrete approximations to the model equations. They are formed by performing a Taylor series expansion of the equations to be solved. All derivatives incorporated in the expansion need to be continuous. The Taylor series can be expanded to attain almost any order of accuracy, however practical methods are only first or second order accurate.

Table 2.1 shows three different approaches to discretising a partial derivative in a space variable [30, p. 18].

Table 2.1: Three methods for discretising

Method	expression	accuracy
Forward difference
Backward difference
Centred difference

As an example the one dimensional model transport equation (Equation 2.8) can be discretised to give Equation 2.9. Equation 2.8 is an expression for the transport of vorticity that is applicable in one dimensional incompressible flow.

$$\frac{\partial \zeta }{\partial t}=-u\frac{\partial \zeta }{\partial x}+\alpha \frac{\partial ^{2}\zeta }{\partial x^{2}}$$ (8)

$$\frac{\zeta _{i}^{n+1}-\zeta _{i}^{n}}{\Delta t}=-\frac{u\ze... ...eta _{i+1}^{n}+\zeta _{i-1}^{n}-2\zeta _{i}^{n}}{\Delta x^{2}}$$ (9)

This is known as the forward time, centred space (FTCS) form of this equation. It is second order accurate ( ).

Control volume approach to discretisation

Another approach to formulating the finite difference equations is to use a control volume. This process is similar to the way the continuous fluid dynamics equations themselves are formulated (Section 2.1). The advantage of control volume approaches are that they are based on macroscopic physical laws, rather than on continuum mathematics. When flow discontinuities occur (eg. shock waves) a discretised form of a continuous equation often does not represent the feature well. On the other hand, by directly representing each physical law in a finite difference equation, insights may be obtained on how to deal with difficult flow features like shock waves. This approach highlights the `numerical simulation' process, exemplified by the PIC (Particle-in-Cell) and FLIC (Fluid-in-Cell) codes developed at Los Alamos ([14],[12]). Roache [30, p. 28] gives an example where he arrives at the same equation as Equation 2.9 for vorticity transport by using a control volume approach.

The control volume approach, or direct numerical simulation, is essentially the approach Milthorpe used in formulating his algorithm (Chapter 4).

Method of Characteristics

Roache [30, p. 226] outlined another approach to solving the fluid equations, namely the `Method of Characteristics' [5]. He states that for supersonic inviscid shock-free flow, where the equations are purely hyperbolic, it is the natural choice for numerical computation. In this well-known method, the computational mesh is not rectangular, nor is its shape known beforehand, but it develops as the solution is explicitly stepped out during the computation. This method gives the most accurate results possible, since the calculation points lie along `characteristics,' across which derivatives may be discontinuous. The greatest limitation of this method is that viscous effects cannot be included, except by boundary-layer methods. Another disadvantage of this method is the complexity of the process of evaluating the location of solution points from the intersection of the characteristic lines.

Godunov type methods

Methods developed for dealing with incompressible flow can often be used for compressible flow problems. However, problems arise when there are strong shocks or high Reynolds' numbers. Errors, in the form of oscillations, occur close to discontinuities due to limitations in the approximation.

Godunov methods calculate the flow in two different reference frames. These methods are complicated but they do treat flow discontinuities better than some other methods. A solution is first found in a Lagrangian frame and then re-mapped to an Eulerian grid. This method can be thought of as operator splitting based on convective and sound waves.

The following is an outline of the method:

1. Initialisation Step - average the initial distribution over the computational cells:

   
   

   $$u_{j}^{0}=\frac{1}{\Delta x_{j}}\int _{x_{j}-\Delta x_{j}/2}^{x_{j}+\Delta x_{j}/2}u(x)dx.$$ (10)

   
   

2. Reconstruction Step - reconstruct the initial distribution as piecewise polynomials over the computational cells.

   
   

   $$u_{j}(x)=P_{j}(x),$$ (11)

   
   

   where is a polynomial in cell .

3. Solution in the Small Step - solve the initial-value problem at each cell interface where discontinuities can exist:

   
   

   $$u(x,t)=E(x,t-t^{n})\cdot u(x,t^{n}),$$ (12)

   
   

   where symbolically represents the evolution operator given by the solution to the Riemann shock-tube problem.

4. Averaging Step - re-average the solution over the grid cells. Given the solution operator in the previous step.
5. Go back to the reconstruction step.

Rider [29] reviewed a number of approximate Riemann solvers that can be used in Step 3. He showed that the choice of a good solver can be used to trade off accuracy with computational efficiency. Care must be taken with some of these solvers because of their inherent smearing of contact discontinuities.

Particle-in-Cell, PIC

Harlow [14] outlined the PIC method, developed at Los Alamos Laboratory in 1955. It is another combined Eulerian-Lagrangian scheme which has been applied successfully to solve a wide variety of multi-fluid problems where the fluid distortions are large. It utilises Lagrangian fluid particles to transport mass, momentum and energy through an Eulerian mesh of cells. While the use of these particles facilitates the calculation of multi-fluid problems, it also results in non-physical fluctuations of the fluid quantities. The dual-mesh system of the PIC method makes great demands on computer memory capacity and on calculation time.

Fluid-in-Cell, FLIC

Gentry, et al. [12] described an Eulerian finite difference method called `Fluid-in-Cell' that can be used to solve time-dependent equations of motion for the compressible flow of a fluid. It is based on the Particle-in-Cell method (Section 2.2.5), though the transport calculation is different and it does not require the use of particles. Elimination of particles allows solutions that are free of the fluctuations characteristic of PIC. The FLIC method will be described in some detail because it has similarities with Milthorpe's method.

The difference equations

Step 1

Pressure is calculated for each cell and intermediate values for , , and are calculated based on pressure gradient effects alone. An artificial viscosity in the form of a pressure term is included. This viscosity term is of the form:

$$\begin{array}{rll} q^{n}_{i+\frac{1}{2},j} & =Bc^{n}_{i+\fra... ...}>v^{n}_{i,j+1}\\ & =0\, \, \mathrm{otherwise} & \end{array}$$ (13)

The quantity determines the maximum value of the Mach number at a cell interface for which the artificial viscosity will be applied. The quantity determines the magnitude of the viscous pressure term. It should be large enough to ensure stability but small enough to avoid obscuring important details of the solution. Generally, need not exceed a value of 0.5.

Step 2

Transport effects are then calculated. The mass, which flows from cell to cell, is directly proportional to the density of the donor cell (the cell from which the fluid is flowing) and the velocity. Transport of momentum components and energy are then accomplished, assuming that the mass which crosses the cell boundaries carries the velocity components and the specific internal energy of the donor cell.

Boundary conditions

Three types of mesh boundaries have been used in the FLIC method: input boundaries, continuative output boundaries, and reflective boundaries. In all three cases the boundary criteria reduce to a determination of the values of the flow variables in fictitious cells outside the calculating mesh.

At any given time, input boundary cells are constant along the boundary, but they may vary with time to represent, for example, a decaying shock wave. At a continuative output boundary, the flow quantities in the fictitious boundary cells are defined in such a way that the normal space derivatives of the variables vanish at the boundary. If the output boundary coincides with cell boundaries, the flow variables in each fictitious exterior cell are set to equal the value in the adjacent interior cell.

The difference equations for those cells adjacent to the bounding surface of a rigid body must be modified to ensure that there is no flow of mass or energy across the boundary. This requires the normal velocity component to be zero at the body surface and is also achieved using fictitious cells for those cases where the reflective surface coincides with a cell boundary. A more complicated approach is used for a general body shape.

Stability and accuracy

The dominant error terms that are introduced by the numerical approximations can be obtained by performing a Taylor series expansion of the finite difference equations, neglecting all terms that are of order or . The form of these error terms is shown below for the case where the solution is a function of one space variable and time [12]:

$$\begin{array}{\vert r\vert l\vert l\vert} \hline & \mathrm{A... ...\partial x}\frac{\partial I}{\partial x} \\ \hline \end{array}$$ (14)

where,

$$\epsilon =\frac{1}{2}\vert u\vert\delta x$$ (15)

There are two types of error terms: artificial viscosity terms that contain and ; and truncation terms, which contain and [12]. The effect of the error terms is to cause an artificial diffusion of mass, momentum and energy to occur in regions where there are large variations.