Elements of Finite-Difference Methods

The method of characteristics described in the previous section legitimately can be considered a part of computational fluid dynamics because it uses discrete algebraic forms of the governing equations (such as Equations (13.17) and (13.18)) which are solved at discrete points in the flow (the characteristic mesh illustrated in Figure 13.5). However, most authors consider that CFD is represented by mainly finite difference and finite volume techniques, such as are discussed in Reference 64, and the method of characteristics is usually not included in the study of CFD. The purpose of this section is to give you the flavor of finite-difference techniques by describing one particular

method that is readily applicable to a number of compressible flow problems. The method discussed here is representative of mainstream CFD, but it is just the tip of the iceberg of CFD. The intensive research in CFD since 1960 has produced a multitude of different algorithms and philosophies, and it is far beyond the scope of this book to go into the details of such work. See Reference 64 for an in-depth presentation of CFD at the introductory level. In addition, you are strongly encouraged to read the current literature in this regard, in particular the AIAA Journal, Computers and Fluids, and the Journal of Computational Physics. Finally, in this chapter we are dealing with numerical solutions of inviscid supersonic flows. See Reference 21 for an expanded discussion of finite difference methods applied to supersonic flows.

How do we use the finite differences obtained here? Imagine that a flow in xy space is covered by the mesh shown in Figure 13.2b. Assume there are N grid points. At each one of these grid points, evaluate the continuity, momentum, and energy equations with their partial derivatives replaced by the finite-difference expressions derived above. For example, replacing the derivatives in Equations (7.40), (7.42a and b), and (7.44) with finite differences, along with Equations (7.1) and (7.6a), we obtain a system (over all N grid points) of 6N simultaneous nonlinear algebraic equations in terms of the 6N unknowns, namely, p, u, v, p, T, and e, at each of the N grid points. In principle, we could solve this system for the unknown flow variables at all the grid points. In practice, this is easier said than done. There are severe problems in solving such a large number of simultaneous nonlinear equations. Moreover, we have to deal with problems associated with numerical instabilities that sometimes cause such attempted solutions to “blow up’’ on the computer. Finally, and most importantly, we must properly account for the boundary conditions. These considerations make

all finite-difference solutions a nontrivial endeavor. As a result, a number of spe­cialized finite-difference techniques have evolved, directed at solving different types of flow problems and attempting to increase computational efficiency and accuracy. It is beyond the scope of this book to describe these difference techniques in detail. However, one technique in particular was widely used during the 1970s and 1980s. This is an approach developed in 1969 by Robert MacCormack at the NASA Ames Research Center. Because of its widespread use and acceptance at the time, as well as its relative simplicity, we will describe MacCormack’s technique in enough detail to give you a reasonable understanding of the method. This description will be carried out in the context of the following example.

Consider the two-dimensional supersonic how through the divergent duct shown in Figure 13.8a. Assume the how is supersonic at the inlet, and that all properties are known at the inlet. That is, the how-field variables at grid points (1, 1), (1, 2), (1, 3),

(1,4) , and (1,5) are known. The duct is formed by a fiat surface at the bottom and a specified contour, ys = f(x), at the top. In addition, assume that the how is inviscid, adiabatic, and steady, and with no body forces. It can be rotational or irrotational—the method of solution is the same. The governing equations are obtained from Equations (7.40), (7.42a and b), (7.44), (7.1), and (7.6a), which yield

[13.35]

Let us express these equations in slightly different form, as follows. Multiplying Equation (13.30) by u, and adding the results to Equation (13.31), we have

Similarly, multiplying Equation (13.30) by v, and adding the result to Equation

(13.32) , we obtain

d(puv) d (pv2 + p)

dx dy

Then, Equations (13.30) and (13.36) to (13.38) become

Equations (13.40) to (13.43) are the continuity, * and у momentum, and energy equations, respectively—but in a slightly different form from those we are used to seeing. The above form of these equations is frequently called the conservation form. Let us now treat F, G, H, and К as our primary dependent variables; these quantities are called flux variables, in contrast to the usual p, p, T, u, v, e, etc., which are called primitive variables. It is important to note that once the values of F, G, H, and К are known at a given grid point, the primitive variables at that point can be found from Equations (13.39a to d) and

p = pRT [13.44]

e = cvT [13.45]

V2 = u2 + v2 [13.46]

That is, Equations (13.39a to d) and (13.44) to (13.46) constitute seven algebraic

equations for the seven primitive variables, p, u, v, p, e, T, and V.

Let us return to the physical problem given in Figure 13.8a. Because the duct diverges, it is difficult to deal with an orthogonal, rectangular mesh; rather, a mesh which conforms to the boundary of the system will be curved, as shown in Figure 13.8a. On the other hand, to use our finite-difference quotients as given in Equation (13.27), (13.28), or (13.29), we desire a rectangular computational mesh. Therefore, we must transform the curved mesh shown in Figure 13.8a, known as the physical plane, to a rectangular mesh shown in Figure 13.8b, known as the computational plane. This transformation can be carried out as follows. Define

|=jc [13.47a]

У

h = —

ys

where ys = f(x) [13.47b]

In the above transformation, t] ranges from 0 at the bottom wall to 1.0 at the top wall. In the computational plane (Figure 13.8b), t] = constant is a straight horizontal line, whereas in the physical plane, t] = constant corresponds to the curved line shown in Figure 13.8. Because we wish to apply our finite differences in the computational plane, we need the governing equations in terms of f and r rather than x and у. To ac­complish this transformation, apply the chain rule of differentiation, using Equations

[13.48]

[13.49]

.43) become

Note in the above equations that the £ derivatives are on the left and the i] derivatives are all grouped on the right.

Let us now concentrate on obtaining a numerical, finite-difference solution of the problem shown in Figure 13.8. We will deal exclusively with the computational plane, Figure 13.8b, where the governing continuity, x and у momentum, and energy equations are given by Equations (13.50) to (13.53), respectively. Grid points (1, 1), (2, 1), (1, 2), (2, 2), etc., in the computational plane are the same as grid points (1, 1), (2, 1), (1, 2), (2, 2), etc., in the physical plane. All the flow variables are known at the inlet, including F, G, H, K. The solution for the flow variables downstream of the inlet can be found by using MacCormack’s method, which is based on Taylor’s series expansions for F, G, Я, and К as follows:

gm-j=gij+&Lm

/дН

Xw-j=Xij+(wL&t

In Equations (13.54a to d), F, G, Я, and К at point (i, j) are considered known, and these equations are used to find F, G, Я, and К at point (/ + 1, j) assuming that we can calculate the values of (9F/9£)ave, (9G/9£)ave, etc. The main thrust of MacCormack’s method is the calculation of these average derivatives. Examining

Equations (13.54a to d), we find that this finite-difference method is clearly a “down­stream marching” method; given the flow at point (i, j) we use Equations (13.54a to d) to find the flow at point (i + l, j). Then the process is repeated to find the flow at point (i + 2, j), etc. This downstream marching is similar to that performed with the method of characteristics.

The average derivatives in Equations (13.54a to d) are found by means of a straightforward “predictor-corrector” approach, outlined below. In carrying out this approach, we assume that the flow properties are known at grid point (i, j), as well as at all points directly above and below (г, j), namely, at (г, j + 1), (г, j +2), (г, j — 1), (ij -2), etc.