Appendix VI

Appendix VI
Development of Finite-Difference Equation for Steady-State Flow in a Homogeneous, Isotropic Medium

The partial differential equation that describes steady-state flow in a two-dimensional, homogeneous, isotropic region of flow (Section 2.11) is Laplace’s equation:

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} (A6.1)

To find the finite-difference equation for an interior node in the nodal grid used to discretize the region of flow, we must replace the second-order partial derivatives in Eq. (A6.1) by differences. Let us consider the first term of the equation first. Recall that the definition of the partial derivative with respect to x of a function of two variables h(x ,z) is

\frac{\partial h}{\partial x} = \substack {lim \\ \Delta x \rightarrow 0} \frac{h(x + \Delta x, z) - h(x, z)}{\Delta x} (A6.2)

On a digital computer it is impossible to take the limit as \Delta x \rightarrow 0, but it is possible to approximate the limit by assigning to \Delta x some arbitrarily small value; in fact, we can do so by designing a nodal network with a mesh spacing of \Delta x.

For any value of z, say z0, we can expand h(x, z0) in a Taylor’s expansion about the point (x0, z0) as follows:

h(x, z_0) = h(x_0, z_0) + (x - x_0)\frac{\partial h}{\partial x}(x_0, z_0) + \frac{(x - x_0)^2}{2}\frac{\partial^2h}{\partial x^2}(x_0, z_0) + ... (A6.3)

If we let x = x_0 + \Delta x (this is known as a forward difference) and abandon all the terms of order greater than unity, we can approximate \partial h/\partial x by

\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - h(x_0, z_0)}{\Delta x} (A6.4)

The abandoned terms of the Taylor expansion represent the truncation error in the finite-difference approximation.

We can obtain a similar expression to Eq. (A6.4) by substituting the backward difference, x = x_0 - \Delta x, into Eq. (A6.3). This yields

\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0, z_0) - h(x_0 - \Delta x, z_0)}{\Delta x} (A6.5)

To obtain the approximation for \partial^2h/\partial x^2, we write the difference equation in terms of \partial h/\partial x using a forward-difference expression:

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{\frac{\partial h}{\partial x}(x_0 + \Delta x, z_0) - \frac{\partial h}{\partial x}(x_0, z_0)}{\Delta x} (A6.6)

and substitute the backward-difference expression Eq. (A6.5) in Eq (A6.6) to get

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - 2h(x_0, z_0) + h(x_0 + \Delta x, z_0)}{(\Delta x)^2} (A6.7)

In a similar manner, we can develop a difference expression for \partial^2h/\partial x^2 as

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0, z_0 + \Delta z) - 2h(x_0, z_0) + h(x_0, z_0 - \Delta z)}{(\Delta z)^2} (A6.8)

For a square grid, \Delta x = \Delta z; adding Eqs. (A6.7) and (A6.8) to for Laplace’s equation yields

\frac{1}{(\Delta z)^2}[h(x_0 + \Delta x, z_0) + h(x_0 - \Delta x, z_0) + h(x_0, z_0 + \Delta z) + h(x_0, z_0 - \Delta z) + 4h(x_0, z_0) = 0] (A6.9)

If we let (x0, z0) be the nodal point (i, j), Eq. (A6.9) can be rearranged to yield

h_{ij} = \frac{1}{4}(h_{i+1,j} + h_{i-1,j} + h_{i,j+1} + h_{i,j-1}) (A6.10)

which is identical to Eq. (5.24).