Manual sketching of flownets is an acceptable way of dealing with two-dimensional groundwater water flow and seepage in soils. However, as the geometry gets more complicated and the flow becomes anisotropic, manual sketching of flownets become more tedious and less accurate. In this case, finite difference solution to flow of water through soils can be adopted.

Many flow problems can be considered as two-dimensional, and in cases where permeability has the same value for all directions, the Laplace’s equation for flow becomes;

(∂²H/∂x²) + (∂²H/∂z²) = 0 ——— (1)

For problems with orthogonal geometry and Laplace’s equation, such as the flow around a pervious barrier shown in Figure 1, the finite difference method provides a quick and accurate solution. This entails replacing the continuous soil cross-section delimited by ABCD with a pattern of discrete points on an orthogonal grid within the cross-section. At each grid point, the governing differential equation can be expressed approximatively in terms of H values at that point and adjacent points.

Previously, we solved the elastic analysis of simply supported thin plates using the finite difference method. We will do the same to solve Laplace’s equation to determine two-dimensional confined flow through soils. Let us consider a grid of a flow domain, as shown in Figure 2, where (i, j) is a nodal point.

Using Taylor’s theorem, we have;

k_x(∂²H/∂x²) + k_z(∂²H/∂z²) = k_x/∆x² (h_i+1,j + h_i-1,j – 2h_i,j) + k_z/∆z² (h_i,j+1 + h_i,j-1 – 2h_i,j) = 0 ——— (2)

Let α = k_x/k_z and ∆x = ∆z (i.e., we subdivide the flow domain into a square grid). Then, solving for h_i,j from Equation (2) gives;

h_i,j = 1/2(1 + α) × (αh_i+1,j + αh_i-1,j + h_i,j+1 + h_i,j-1) ——— (3)

For isotropic conditions, α = 1 (k_x = k_z) and Equation (3) becomes;

h_i,j = 1/4 × (h_i+1,j + h_i-1,j + h_i,j+1 + h_i,j-1) ——— (4)

Since we are considering confined flow, one or more of the boundaries would be impermeable. Flow cannot cross impermeable boundaries and, therefore, for a horizontal impermeable surface;

∂h/∂x = 0 ——— (5)

The finite difference form of Equation (5) is;

∂h/∂x =1/2∆x(h_i,j+1 – h_i,j-1) = 0 ——— (6)

Therefore, h_i,j+1 = h_i,j-1 and, by substitution in Equation (4), we get;

h_i,j = 1/4 × (h_i+1,j + h_i-1,j + 2h_i,j-1) ——— (7)

Various types of geometry of impermeable boundaries are encountered in practice, three of which are shown in Figure 3. For Figure 3a, b, the finite difference equation is;

h_i,j = 1/2 × (h_i+1,j + h_i,j-1) ——— (8)

and, for Figure 3c,

h_i,j = 1/3 × (h_i,_j-1 + h_i+1,j + h_i,j+1 + 0.5h_i-1,j + 0.5h_i,j+1) ——— (9)

The pore water pressure at any node (u_i,j) is;

u_i,j = γ_w(h_i,j – z_i,j) ——— (10)

where z_i,j is the elevation head.

Contours of potential heads can be drawn from discrete values of h_i,j. The finite difference equations for flow lines are analogous to the potential lines; that is, ψ_s replaces h in the above equations and the boundary conditions are specified for ψ_s rather than for h.

The horizontal velocity of flow at any node (v_i,j) is given by Darcy’s law:

v_i,j = k_xi_i,j ——— (11)

Where i_i,j is the hydraulic gradient expressed as;

i_i,j = (h_i+1,j – h_i-1,j)/2∆x ——— (12)

Therefore;

v_i,j = k_x/2∆x × (h_i+1,j – h_i-1,j) ——— (13)

The flow rate, q, is obtained by considering a vertical plane across the flow domain. Let L be the top row and K be the bottom row of a vertical plane defined by column i (Figure 2). Then the expression for q is;

q = k_x/4[h_i+1,L – h_i-1,L + 2∑(h_i+1,j – h_i-1,j) + h_i+1,K – h_i-1,K] ——— (14)

Procedure of Finite Difference Solution to Flow of Water Through Soils

The procedure to determine the distribution of potential head, flow, and porewater pressure using the finite difference method is as follows (Budhu, 2011):

1. Divide the flow domain into a square grid. Remember that finer grids give more accurate solutions than coarser grids, but are more tedious to construct and require more computational time. If the problem is symmetrical, you only need to consider one-half of the flow domain.
   For example, the sheet pile wall shown in Figure 4 is symmetrical about the wall and only the left half may be considered. The total flow domain should have a width of at least four times the thickness of the soil layer. For example, if D is the thickness of the soil layer (Figure 4), then the minimum width of the left half of the flow domain is 2D.
2. Identify boundary conditions, for example, impermeable boundaries (flow lines) and permeable boundaries (equipotential lines).
3. Determine the heads at the permeable or equipotential boundaries. For example, the head along the equipotential boundary AB (Figure 4) is DH. Therefore, all the nodes along this boundary will have a constant head of DH. Because of symmetry, the head along nodes directly under the sheet pile wall (EF) is DH/2.
4. Apply the known heads to corresponding nodes and assume reasonable initial values for the interior nodes. You can use linear interpolation for the potential heads of the interior nodes.
5. Apply Equation (4), if the soil is isotropic, to each node except (a) at impermeable boundaries, where you should use Equation (7), (b) at corners, where you should use Equations (8) and (9) for the corners shown in Figure 3a–c, and (c) at nodes where the heads are known.
6. Repeat item 5 until the new value at a node differs from the old value by a small numerical tolerance, for example, 0.001 m.
7. Arbitrarily select a sequential set of nodes along a column of nodes and calculate the flow, q, using Equation (14). It is best to calculate q’ = q for a unit permeability value to avoid too many decimal points in the calculations.
8. Repeat items 1 to 6 to find the flow distribution by replacing heads by flow q’. For example, the flow rate calculated in item 7 is applied to all nodes along AC and CF (Figure 4). The flow rate at nodes along BE is zero.
9. Calculate the porewater pressure distribution by using Equation (10).

References

[1] Budhu M. (2011): Soil Mechanics and Foundations (3rd Edition). Pearson Education
[2] Tomlinson M. J. (2001): Foundation Design and Construction (7th Edition). John Wiley and Sons, Inc

Laplace’s equation is used to describe the movement of water through soils. By comparison, the flow of water through soils is analogous to the steady-state heat flow and steady-state current flow in homogeneous conductors. Flownets can be used to calculate the flow of water through soils based on the Laplace’s equation. The common form of Laplace’s equation for the flow of water through two-dimensional soils is:

k_x(∂²H/∂x²) + k_z(∂²H/∂z²) = 0 ——– (1)

where H is the total head and k_x and k_z are the hydraulic conductivities in the X and Z directions. The condition that the changes in hydraulic gradient in one direction are balanced by changes in the other directions is expressed by Laplace’s equation.

The assumptions in Laplace’s equation are:

• Darcy’s law is valid.
• Irrotational flow (vorticity) is negligible. This assumption leads to the following two-dimensional relationship in velocity gradients.

∂v_z / ∂Z = ∂v_x / ∂X

where v_z and v_x are the velocities in the Z and X directions, respectively. This relationship is satisfied for a uniform flow field and not a general flow field. Therefore, we will assume all flows in this chapter are uniform, i.e., v_z = v_x = constant.
• There is inviscid flow. This assumption means that the shear stresses are neglected.
• The soil is homogeneous and saturated.
• The soil and water are incompressible (no volume change occurs).

Laplace’s equation is also called the potential flow equation because the velocity head is neglected. If the soil is an isotropic material, then k_x = k_z and Laplace’s equation becomes;

(∂²H/∂x²) + (∂²H/∂z²) = 0 ——– (2)

Any differential equation requires knowledge of the boundary conditions in order to be solved. Since the boundary conditions of the majority of “real” structures are complex, an analytical or closed-form solution cannot be obtained for these structures. Using numerical techniques such as finite difference, finite element, and boundary element, it is possible to obtain approximate solutions.

We can also attempt to replicate the flow through the actual structure using physical models. There are two major techniques for solving Laplace’s equation. The first is an approximation known as flownet sketching, and the second is the finite difference method. In this article, we are going to focus on flownet sketching.

Flownet Sketching

The flownet sketching technique is straightforward and adaptable, and it represents the flow regime. It is the preferred method of analysing flow through soils for geotechnical engineers. Before delving into these solution techniques, however, we will establish a few key conditions necessary to comprehend two-dimensional flow.

The solution of Equation (1) is solely dependent on the total head values within the flow field in the XZ plane. Let us introduce a velocity potential (ξ) that describes the variation of total head in a soil mass as follows:

ξ = kH ——– (3)

where k is a generic hydraulic conductivity. The velocities of flow in the X and Z directions are;

v_x = k_x(∂H/∂x) = ∂ξ/∂x ——– (4a)
v_z = k_z(∂H/∂z) = ∂ξ/∂z ——– (4b)

The inference from Equations (4a) and (4b) is that the velocity of flow (v) is normal to lines of constant total head, as illustrated in Figure 1 The direction of v is in the direction of decreasing total head. The head difference between two equipotential lines is called a potential drop or head loss.

If we draw lines that are tangent to the flow velocity at each point in the flow field in the XZ plane, we will obtain a series of lines that are normal to the equipotential lines. These lines are known as streamlines or flow lines (Figure 1). A flow line represents the expected path of a particle of water in a steady-state flow. ψ_s is a stream function that represents a streamline family (x, z). According to the stream function, the components of velocity in the X and Z directions are as follows:

v_x = ∂ψ_s / ∂z ——– (5a)
v_z = ∂ψ_s / ∂x ——– (5b)

Since flow lines are normal to equipotential lines, there can be no flow across flow lines. The rate of flow between any two flow lines is constant. The area between two flow lines is called a flow channel (Figure 1). Therefore, the rate of flow is constant in a flow channel.

Criteria for Sketching Flownets

A flownet is a graphical representation of a flow field that satisfies Laplace’s equation and comprises a family of flow lines and equipotential lines. A flownet must meet the following criteria (Budhu, 2011):

1. The boundary conditions must be satisfied.
2. Flow lines must intersect equipotential lines at right angles.
3. The area between flow lines and equipotential lines must be curvilinear squares. A curvilinear square has the property that an inscribed circle can be drawn to touch each side of the square and continuous bisection results, in the limit, in a point.
4. The quantity of flow through each flow channel is constant.
5. The head loss between each consecutive equipotential line is constant.
6. A flow line cannot intersect another flow line.
7. An equipotential line cannot intersect another equipotential line.

An infinite number of flow lines and equipotential lines can be drawn to satisfy Laplace’s equation. However, only a few are required to obtain an accurate solution. The procedure for constructing a flownet is described next.

Flownet for Isotropic Soils

According to Budhu (2011), the procedure for constructing the flownet of isotropic soils are as follows;

1. Draw the structure and soil mass to a suitable scale.
2. Identify impermeable and permeable boundaries. The soil–impermeable boundary interfaces are flow lines because water can flow along these interfaces. The soil–permeable boundary interfaces are equipotential lines because the total head is constant along these interfaces.
3. Sketch a series of flow lines (four or five) and then sketch an appropriate number of equipotential lines such that the area between a pair of flow lines and a pair of equipotential lines (cell) is approximately a curvilinear square. You would have to adjust the flow lines and equipotential lines to make curvilinear squares. You should check that the average width and the average length of a cell are approximately equal by drawing an inscribed circle. You should also sketch the entire flownet before making adjustments.

The flownet in confined areas between parallel boundaries typically consists of elliptical and symmetrical flow lines and equipotential lines (Figure 2). Avoid abrupt changes between straight and curved flow and equipotential lines. Transitions should be smooth and gradual. For certain problems, portions of the flownet are enlarged, are not curvilinear squares, and do not satisfy Laplace’s equation.

For instance, the portion of the flownet beneath the base of the sheet pile in Figure 2 is not composed of curvilinear squares. Check these sections to ensure that repeated bisection results in a point for a precise flownet.

Another example of flownet are shown in Figures 3. Figure 2 shows a flownet for a sheet pile wall, and Figure 3 shows a flownet beneath a dam. In the case of the retaining wall, the vertical drainage blanket of coarse-grained soil is used to transport excess porewater pressure from the backfill to prevent the imposition of a hydrostatic force on the wall. The interface boundary, is neither an equipotential line or a flow line. The total head along the boundary is equal to the elevation head.

Flow Rate

Let the total head loss across the flow domain be ΔH, that is, the difference between upstream and downstream water level elevation. Then the head loss (Δh) between each consecutive pair of equipotential lines is;

Δh = ΔH/N_d ——– (6)

where N_d is the number of equipotential drops, that is, the number of equipotential lines minus one. In Figure 1, ΔH = H = 8 m and N_d = 18. From Darcy’s law, the flow through each flow channel for an isotropic soil is;

q = Aki = (b × 1)k(Δh/L) = kΔh(b/L) = k(ΔH/N_d)(b/L) ——– (7a)

where b and L are defined as shown in Figure 14.3. By construction, b/L = 1, and therefore the total flow is;

q = kΔH(N_f/N_d) ——– (7b)

where N_f is the number of flow channels (number of flow lines minus one). In Figure 1, N_f = 9. The ratio N_f /N_d is called the shape factor. Finer discretization of the flownet by drawing more flow lines and equipotential lines does not significantly change the shape factor. Both N_f and N_d can be fractional. In the case of anisotropic soils, the quantity of flow is;

q = ΔH(N_f/N_d)√(k_xk_z) ——– (8)

Summarily, Flow nets are typically designed for homogeneous, isotropic porous media undergoing saturated flow to known boundaries. There are extensions to the basic method that make it possible to solve the following cases:

• inhomogeneous aquifer: matching conditions at property boundaries
• anisotropic aquifer: drawing the flownet in a transformed domain and scaling the results differently in the principal hydraulic conductivity directions before returning the solution
• one boundary is a seepage face: iteratively solving for both the boundary condition and the solution throughout the domain

The method is typically applied to these types of groundwater flow problems, but it can be applied to any problem described by the Laplace equation, such as the flow of electric current through the earth.

References
Budhu M. (2011): Soil Mechanics and Foundations (3rd Edition). John Wiley & Sons, Inc.