Milthorpe's Algorithm

Milthorpe [23] noted that many existing methods at the time for computing supersonic flows were based on continuum mechanics. These methods predicted a flow by solving a set of partial differential equations subject to certain boundary conditions. When these methods were applied to hyperbolic problems, as those encountered in hypervelocity research, instabilities often occurred near shocks or other discontinuities.

Milthorpe set out to produce a method that would capture flow discontinuities accurately. His method does not attempt to solve the partial differential equations, but rather examines the behavior of the fundamental quantities: mass, momentum and energy. He accomplished direct numerical simulation where the physical laws are applied to the fluid variables on small volume elements or cells.

The approach is similar to a finite difference method derived from the same physical laws, except that conditions have been placed on many of the calculations. For example, diffusion is limited to a certain proportion of a fluid quantity in a cell. Another example is that pressure difference effects are applied with an upwinding coefficient, which depends on some locally averaged quantities. It could be said that the progress of the flow calculation is managed so that the calculation proceeds in a realistic fashion. This is in contrast to a brute force finite difference, or other equation solving method where complete trust is placed on the governing equations and only the solution of these equations is deemed important.

The algorithm used to calculate the flow values is based on conservation of mass, momentum and energy as the simulation proceeds over discrete time steps. The flow domain is divided into rectangular cells. Each cell contains, at any given time, a certain block of material that possesses a certain amount of momentum and energy. In the course of a time step this block is modified by pressure gradient effects and diffusion. This block, with its associated mass, momentum, and energy, is then convected a small distance depending on the local velocity vector. Greatest in magnitude, at the velocities modelled, is convective transport of mass, momentum and energy due to the velocity of the fluid. Following this are the effects of pressure gradients with the effect of diffusive transport being least significant in magnitude.

Conservation algorithm

Phase 1 - Diffusion

The fluid quantities of mass, momentum and energy are modified by the process of diffusion according to Equation 4.1:

$\begin{displaymath} dQ^{k}_{i,j}=\delta t\alpha ^{k}[dy(Q^{k}_{i-1,j}-2Q^{k}_{i,j}+Q^{k}_{i+1,j})+dx(Q^{k}_{i,j-1}-2Q^{k}_{i,j}+Q^{k}_{i,j+1})] \end{displaymath}$

(17)

**Table 4.1:** Diffusion coefficients
viscosity	$\mu$
molecular diffusivity	$D$
thermal conductivity	$k$

Sutherland's law (Equation 4.2) [16] is used to calculate the temperature dependence of the viscosity.

$\begin{displaymath} \mu =\mu _{0}\left( \frac{T}{T_{0}}\right) ^{3/2}\frac{T_{0}+S}{T+S} \end{displaymath}$

(18)

$\begin{displaymath} D=\frac{\mu }{\rho S_{c}}\end{displaymath}$

$\begin{displaymath} k=\frac{\mu c_{p}}{P_{r}}\end{displaymath}$

In addition, momentum and energy are modified by diffusion of those quantities along velocity gradients and energy is modified by diffusion along temperature gradients. These are calculated using an expression similar to Equation 4.1.

In all cases, the amount of a fluid quantity that is transported by diffusion is limited to half the amount of that quantity in the donor cell.

Phase 2 - Pressure gradient effects

The momentum and energy in the cell are modified by the action of the pressure gradients. The pressure difference between adjacent cells affects the momentum in the cells according to the following relations:

$\begin{displaymath} \delta (\rho u_{i,j})=\delta tdy(p_{i-1,j}-2p_{i,j}+p_{i+1,j}) \end{displaymath}$

(19)

$\begin{displaymath} \delta (\rho v_{i,j})=\delta tdx(p_{i,j-1}-2p_{i,j}+p_{i,j+1}) \end{displaymath}$

(20)

At the same time as the pressure gradients are transporting momentum between adjacent cells, if the local velocity is non-zero they are also doing work, and hence transporting energy between the cells. This energy transport is obtained from the momentum transport calculated in Equations 4.3 and 4.4 by multiplying by the local velocity.

Upwinding factor

Milthorpe found that by including some upwinding terms in the calculation of pressure gradient effects the accuracy of the solution was enhanced. Adding an upwinding factor to Equation 4.3 gives:

$\begin{displaymath} \delta (\rho u)=\frac{\delta tdy}{2}\left[ (1+F)p_{i-1,j}-2Fp_{i,j}-(1-F)p_{i+1,j}\right] \end{displaymath}$

(21)

The expression used to calculate the energy transport due to pressure gradients also included a second upwinding parameter, $e_{f}$ , that is used to combine different amounts of upstream and downstream velocity in the calculation. The rather complicated equation for this is:

$\begin{eqnarray*} \delta E & = & \frac{\delta tdy}{2}\left\{ -\left[ (1+F)p_{i,j... ...F)u_{i-1,j}+(1+F)u_{i,j}\right] +e_{f}u_{i-1,j}\right] \right\} \end{eqnarray*}$

Phase 3 - Convection

Convection of the cells is illustrated in Figure 4.1 for convection with positive $u$ - and $v$ -velocities. A cell is regarded as moving with unchanged dimensions relative to the mesh of cell boundaries, so that it overlaps adjacent cells. The letter $A$ indicates the area of the cell that remains in the original cell at the end of the time step $u\delta t$ , while $B$ , $C$ and $D$ represent the areas of a cell that are located in adjacent cells. The quantities of mass, momentum and energy in the cell are divided in proportion to the areas $A$ , $B$ , $C$ and $D$ . Those proportions represented by $B$ , $C$ and $D$ are subtracted from the quantities in the original cell and added to the appropriate adjacent cells.

**Figure 4.1:** Illustration of cell convection
$\includegraphics{newconvection_html.eps}$

Calculation of dependent flow variables

When the adjustments to mass, momentum and energy have been completed for all cells, the dependent flow variables density ( $\rho$ ), velocity ( $u$ and $v$ ), pressure ( $p$ ) and temperature ( $T$ ) can then be calculated. These values are then used in the next iteration. The variables $\rho$ , $u$ and $v$ are calculated in the usual way. However, $T$ is calculated from the cell's energy and velocity using:

$\begin{displaymath} T=\frac{E-\frac{1}{2}\left( u^{2}+v^{2}\right) }{mc_{p}}\end{displaymath}$

Stability limits

It is computationally more economical to use the largest possible time step in the computations, but there are a number of limits imposed by accuracy and stability considerations. The first limit is that since convection occurs only into the adjacent cells, the distance travelled in a time step must not exceed the cell dimension. The second limit is that due to over prediction of diffusion, as described by Roache [30]. This restricts the time step to:

$\begin{displaymath} \Delta t<\frac{\rho \Delta x^{2}}{2\alpha_k}, \end{displaymath}$

(22)

A third limit arises from pressure gradients causing a sufficiently large change in momentum to alter the sign of the momentum in a cell.

At each iteration, these three limits, and several others, are checked to see if the time step is too large. If this is so, the time step is reduced by a certain amount, otherwise it is increased by a smaller amount. This process stabilises the time step to a particular value while flow conditions are effectively steady and adjusts to changes, such as when the flow reaches an object or a wall. This self adjusting time step allows the calculation to progress as fast as possible whilst keeping the system stable.

Boundary conditions

Milthorpe implemented boundary conditions by testing for certain array indices as the algorithm looped through the arrays of fluid values.

Upstream boundary

The simplest boundary implemented was the `upstream' boundary. Most of the problems that Milthorpe's method has been applied to involve a transient solution of a flow where the fluid enters from one side of the computational domain. The side where the fluid enters is called the upstream boundary. Upstream boundary cells are reset to their initial values at the start of each iteration. They take part in the process of diffusion and convection into their neighbouring cells, though their own fluid values are not updated.

Symmetric boundary

A feature that often appears in the kind of problems Milthorpe's algorithm has been applied to is a symmetry surface or slip plane. There is no flow across this surface, no shear stress and no diffusion of energy or material. These conditions occur at a plane of symmetry, on the centre line of an axisymmetric flow and at a solid surface at which it is not required to enforce a no-slip condition. To implement this, values of all variables in a `dummy' layer immediately outside the computational domain are set equal to those inside, ensuring no diffusion across the boundary. The only restriction on the values in the surface layer is that the normal velocity component cannot be directed into the boundary. The values in any cell are, of course, only a representative average of the point values, so non-zero velocities into the computation domain do not violate the condition of zero flow outward across the boundary. The normal velocity in the dummy layer is set to zero and hence there is no flow inwards across the boundary.

Freestream boundary

This boundary condition is by far the most difficult to implement successfully. In contrast to the predetermined downstream condition, such as a set value of pressure in the freestream, no information is known about the conditions at a free downstream boundary, except that no influence from the boundary can be allowed to propagate into the flow upstream. Many experimental results are concerned with external flows and these require free downstream boundaries if the size of the computation domain is to be at all reasonable.

Roache [30, p. 165] discussed what he called the `downstream paradox'. Normally in supersonic flow, downstream disturbances cannot effect the flow upstream. The paradox in CFD is that the outflow condition, which is downstream, is more important in supersonic flow than in subsonic flow. Roache described an example scenario demonstrating how this can be the case. He outlined [30, p. 279] various approaches used to deal with a downstream boundary in compressible flow - most relying on some kind of extrapolation. Roache described an `upper' boundary separately and outlined a method that follows simple wave characteristics to determine a direction of extrapolation.

Milthorpe implemented a freestream boundary by projecting the contours of the flow variables across the boundary. The direction of the contours at the boundary was determined for each cell by calculating the direction that gave the minimum average deviation from the contour value. The value outside the boundary was taken as the mean value along that direction for a small distance into the region. This use of extrapolation in the direction of a contour sets Milthorpe's method apart from those described by Roache. Milthorpe's freestream boundary implementation produced quite a straight projection of the contours over the boundary [24].

Comparison with FLIC

Milthorpe's method is similar to the Fluid-in-Cell (FLIC) [12] method described in Section 2.2.6, both being derived from a control volume approach (Section 2.2.2). They are both intended to be used for compressible flow problems and the implementation of convection is very similar in the two methods. The main differences, however, are:

Pressure gradient effects

The difference equation used in the FLIC method to calculate pressure gradient effects in velocity is:

$\begin{displaymath} \tilde{u}^{n}_{i,j}=u^{n}_{i,j}-\frac{\delta t}{\rho ^{n}_{i,j}\delta x}(P_{i+\frac{1}{2},j}-P_{i-\frac{1}{2},j}) \end{displaymath}$

(23)

An analysis of Equations 4.3 and 4.7, shows that they differ by a factor of $\delta x\delta y$ , which is the cell volume. The FLIC version has an artificial pressure term and Milthorpe's has an upwinding factor, $F$ . The discretisation scheme is also a little different. The Milthorpe method calculates values at cell-centred node points, whereas FLIC calculates values at cell boundaries.

Transport effects

The FLIC method calculates a mass flux at each cell boundary proportional to the mass of the donor cell and the intermediate value of the velocity calculated ( $\tilde{u}^{n}_{i,j}$ ). Transport of velocity and internal energy is calculated based on the quantity of mass transported.

This differs from Milthorpe's convection algorithm in that Milthorpe transports a proportion of mass, momentum and energy independently of each other - depending only on the factor $u\Delta t$ . Here $u$ is equivalent to FLIC's $u^{n}$ not the intermediate value $\tilde{u}^{n}$ .