## David R. Steward

Print publication date: 2020

Print ISBN-13: 9780198856788

Published to Oxford Scholarship Online: November 2020

DOI: 10.1093/oso/9780198856788.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (oxford.universitypressscholarship.com). (c) Copyright Oxford University Press, 2021. All Rights Reserved. An individual user may print out a PDF of a single chapter of a monograph in OSO for personal use. Subscriber: null; date: 09 December 2021

# Analytic Element Method across Fields of Study

Chapter:
(p.1) 1 Analytic Element Method across Fields of Study
Source:
(p.iii) Analytic Element Method
Publisher:
Oxford University Press
DOI:10.1093/oso/9780198856788.003.0001

# Abstract and Keywords

This chapter introduces the philosophical perspective for solving problems with the Analytic Element Method, organized within three common types of problems: gradient driven flow and conduction, waves, and deformation by forces. These problems are illustrated by classic, well known solutions to problems with a single isolated element, along with their extension to complicated interactions occurring amongst collections of elements. Analytic elements are presented within fields of study to demonstrate their capacity to represent important processes and properties across a broad range of applications, and to provide a template for transcending solutions across the wide range of conditions occurring along boundaries and interfaces. While the mathematical and computational developments necessary to solve each problem are developed in later chapters, each figure documents where its solutions are presented.

# 1.1 Philosophical Perspective

The Analytic Element Method (AEM) conceptualization is quite intuitive. It is founded in the mathematical representations employed across fields of study in science and engineering, which describe a disciplinary perspective of important processes and properties. For many problems, the physical environment may be represented as a collection of elements with distinct properties. Each element interacts with its adjacent environment, and methods are formulated along an element to solve its conditions with nearly exact precision. Collectively, sets of elements interact to transform their environment, and these synergistic interactions are expanded upon for three common types of problems.

### Problem 1. Flow and Conduction

The first type of problem studies a vector field that is directed from high to low values of a function as illustrated in Fig. 1.1. The function is visualized using a gray color pallet that shades the domain from lighter tones where the value of the function is larger to darker tones where it is smaller, and solid black lines lie along equipotential contours that have constant function value. The vector field is computed at discrete locations and arrows are drawn from these points in the direction of the vector with lengths scaled by the vector’s magnitude. Many solutions may also be represented using streamlines, which are everywhere tangent to the vector field and drawn as white lines.

Figure 1.1 Example problems with flow and conduction in a background uniform vector field, with figure numbers indicating where solution methods are presented.

(p.2) The use of analytic elements to representation solutions is presented across broad fields of study. Groundwater flow is formulated in Section 1.2.1 with equipotential surfaces of uniform groundwater elevation (head), the vector field points in the direction of flow, and streamlines represent the trajectories of water particles. A similar representation is used to study flow through a vadose zone in Section 1.2.2 using equipotentials with uniform pressure, and vector fields are developed for incompressible fluid flow in Section 1.2.3. Analytic elements are formulated for thermal conduction in Section 1.2.4 using equipotentials with uniform temperature and the vector field represents conduction of heat. Electrostatics in Section 1.2.5 represent solutions for electric current using electric potentials related to voltage.

A common representation of solutions are formulated using relevant properties and processes from these disciplinary perspectives, but with cross-disciplinary interpretations of results (Bulatewicz et al., 2010, 2013). For example, the solution in Fig. 1.1a was used to study fluid flowing about an impermeable sphere by Lamb (1879, sec.102), and also represents the electric field near an object with dielectric constant equal to zero (Sommerfeld, 1952, p.62). The solution in Fig. 1.1b illustrates electric current through a dielectric inhomogeneity with the geometry of an ellipse or a prolate spheroid (Moon and Spencer, 1961a, sec.6), and has related applications in heat conduction (Carslaw and Jaeger, 1959, sec.16.4). The example in Fig. 1.1c illustrates a vector field for flow and conduction with many interacting elements that collectively shape the solution. A general framework is presented in Section 1.2 for such problems with flow and conduction to achieve a consistent representation with broadly applicable solutions.

### Problem 2. Periodic Waves

A second type of problem studies the interactions of analytic elements with waves. This is illustrated in Fig. 1.2 for a background field of plane waves traveling from left to right that encounter and interact with elements. Wave fields are visualized using two plots for each solution, where the upper panels illustrate the amplitude of waves, and the lower panels illustrate the phase of a wave field. The amplitude shown in Fig. 1.2 is scaled by the (p.3) amplitude of the background wave field and thicker lines are drawn where the amplitude would be equal to the background if the elements did not exist. The phase is contoured at equal position on a wave; for example, if the darker tones represent the highest elevation at the crest of a wave field at a specific time, then the lighter tones represent the lowest elevation at the trough of the wave field at the same time.

Figure 1.2 Example solutions for problems with periodic waves in a field of plane waves. (Reprinted from Steward, 2018, Wave resonance and dissipation in collections of partially reflecting vertical cylinders, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 144(4), Fig. 1, with permission from ASCE.) (Reprinted from Steward, 2020, Waves in collections of circular shoals and bathymetric depressions, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 146, Figs. 2,4, preproduction version DOI 10.1061/(ASCE)WW.1943-5460.0000570, with permission from ASCE.)

Analytic elements are developed for water waves in Section 1.3.1 to illustrate the amplitude and phase of the surface of a sea. The analytic elements for acoustics in Section 1.3.2 illustrate vibrations traveling through solids, liquids, or gases. Boundary conditions are developed for wave reflection where waves are scattered, for example, by the fully reflective element in Fig. 1.2a. This solution was developed by MacCamy and Fuchs (1954) to study the formation of partially standing waves on the windward side of an island, and wave diffraction occurs on the leeward side with calmer water. Analytic elements may (p.4) also be formulated to study wave refraction as illustrated in Fig. 1.2b to represent solutions when water waves travel across a shoal with shallower water, or when acoustic waves travel through an inhomogeneity with a wavelength shorter than that of the surrounding media. Wave refraction also occurs in Fig. 1.2c, where many elements interact to collectively shape the wave field. Analytic elements are formulated for waves in Section 1.3 for problems with reflection, diffraction, and refraction.

### Problem 3. Deformation by Forces

A third type of problem studies the interactions of analytic elements with stresses and displacements. The examples in Fig. 1.3 illustrate stress fields occurring as analytic elements interact with a background field generated by body forces. These solutions are visualized using two panels for each solution: the top panels show principle stresses and mean stress, and the lower panels show displacement and maximum shear. The principle stresses are computed at equally spaced positions and lines are drawn from these locations to indicate direction and magnitude, and isostress contours with uniform mean stress have darker tones, indicating a larger compression. Displacements are shown by vectors scaled (p.5) by its magnitude, and isoshear contours of uniform maximum shear show lighter tones where it has a larger value.

Figure 1.3 Example solutions for problems with deformation by forces in a field of biaxial compression.

The use of analytic elements to solve problems with deformation by forces is formulated in Section 1.4. The example in Fig. 1.3a illustrates the background stress field in geomechanics generated at a depth within a soil column by the body forces of the overlying material. The solution in Fig. 1.3b from Jaeger (1969) illustrates the stress field and displacement formed around a tunnel placed in this background field, where the boundary is traction free with no normal or shear stresses. The example in Fig. 1.3c shows a set of elements in a plate that satisfy a boundary condition of zero displacement where their fixed locations cause a distribution of the stress field. Solutions are developed in Section 1.4.1 for analytic elements with boundary conditions of specified stress or specified displacement. While the examples are shown for stress–strain studies of elasticity in structures and geomechanics, the methods are extensible to the studies of thermal stress and inviscid fluid flow.

# 1.2 Studies of Flow and Conduction

A common problem in engineering and science is to compute the value of a function and an associated vector field across a region of interest, where the vector field is directed from higher to lower values of the function (as in Fig. 1.1). For example, the function might represent groundwater elevation or temperature, with associated vector fields of groundwater flow or thermal conduction. Such problems are formulated in terms of a potential function Φ‎ and a vector field v, and vector calculus provides the basic operations to manipulate these functions.

The gradient, , operates on the scalar function and produces a vector with components equal to the partial derivative of the function. The gradient of a three-dimensional Φ‎(x, y, z) has components

$Display mathematics$
(1.1)

where $x^$, $y^$, and $z^$ are unit vectors pointing in the x, y, and z coordinate directions. The gradient is directed normal to surfaces of constant Φ‎ and directed towards increasing Φ‎, as illustrated in Fig. 1.4a. A second vector operator is the divergence, ·, which is used to quantify the net flux of a vector field through a control volume as illustrated in Fig. 1.4b. The divergence of v with Cartesian components (vx, vy, vz) is obtained by multiplying the normal components of the vector field directed out minus those directed into the control (p.6) volume times the surface area of each face of the parallelepiped with sides of length Δ‎x, Δ‎y, and Δ‎z, and dividing by the volume:

$Display mathematics$
(1.2)

Evaluating the limit of this central difference equation gives

$Display mathematics$
(1.3)

A third vector operator is the curl, ×, which may be written in Cartesian coordinates for the vector v:

$Display mathematics$
(1.4)

The curl is used to quantify infinitesimal rotation as illustrated in Fig. 1.4c, where the component of the curl in the z-direction, $∂vy∂x−∂vx∂y$, is zero when the angles $ϑ1=ϑ2$ in this figure are equal in the limit as Δ‎x and $Δy→0$, and the vector field is irrotational.

There are two functions commonly utilized to solve problems with flow and conduction, a scalar potential function, Φ‎, and a vector potential function, Ψ‎. The Helmholtz theorem states that any vector field, v, may be decomposed into a linear combination of the gradient of this scalar potential and the curl of this vector potential:

$Display mathematics$
(1.5)

Figure 1.4 Vector operators: gradient, divergence, and curl.

Figure 1.5 Two-dimensional vector fields and the stream function.

(p.7) The first corollary to the Helmholtz theorem states that a vector field may be represented using the gradient of the scalar potential if and only if the field is irrotational:

$Display mathematics$
(1.6)

and most problems will be formulated in terms of Φ‎. This sign convention aligns a vector field in the direction of decreasing potential. The second corollary to the Helmholtz theorem states that a vector field may be represented using the curl of the vector potential function if and only if the field is divergence-free:

$Display mathematics$
(1.7)

and use of vector potentials for three-dimensional problems is presented in Section 5.7.2. Vector fields that have both divergence and curl require use of both the scalar and vector potential, while those that are divergence-free and irrotational may be formulated using either function.

A simpler representation exists for divergence-free problems with a two-dimensional vector field that lies in a plane (e.g., $vz=0$) and is invariant in the direction normal to this plane (i.e., $∂vx∂z=∂vy∂z=0$). A streamline in such a flow is defined to be everywhere tangent to the vector field, and this property is satisfied in Fig. 1.5a when

$Display mathematics$
(1.8)

where the incremental step ds along the streamline has components dx and dy in the coordinate directions. The Lagrange (1781) stream function, Ψ‎, is defined to be constant along a streamline in the s-direction, and the chain rule for differentiation gives

$Display mathematics$
(1.9)

(p.8) Together, the last two equations give a relation between the partial derivatives of Ψ‎ and the components of the velocity vector:

$Display mathematics$
(1.10)

A vector potential for a two-dimensional vector field is given by multiplying the Lagrange stream function times $z^$, a unit vector pointing in the positive z-direction (Steward, 2002),

$Display mathematics$
(1.11)

and the curl of this vector in (1.7) reproduces the two-dimensional relationship (1.10). Consequently, any two-dimensional divergence-free vector field may be studied using a Lagrange stream function Ψ‎.

The stream function is important because contours of constant Ψ‎ may be visualized to illustrate streamlines. The stream function also provides the flux occurring between two points, as illustrated in Fig. 1.5b for the two streamlines that pass through the points A and B. This flux is obtained from the integral of the normal component of the vector field along the vertical and horizontal lines connecting A to B,

$Display mathematics$
(1.12)

Substituting the relationship between the vector field and stream function (1.10), integrating, and evaluating the stream function at the limits of integration gives

$Display mathematics$
(1.13)

where Ψ‎A and Ψ‎B are the value of the stream function at points A and B. Thus, the flux between two streamlines is equal to their difference in stream function.

The potential and stream functions both exist for two-dimensional vector fields that are irrotational and divergent-free. For this case, the relations between the vector field and the potential function (1.6) and the stream function (1.10) together give the Cauchy–Riemann equations

$Display mathematics$
(1.14)

Furthermore, such functions Φ‎ and Ψ‎ each satisfy the Laplace equation, which will be discussed further at (2.2a). A flow net exists for irrotational and divergent-free vector fields, where families of equipotentials with constant Φ‎ and families of streamlines with constant Ψ‎ are mutually orthogonal.

Problem 1.1 Develop expressions for the vector field (components vx and vy) associated with the following potential and stream functions, where $r=x2+y2$, and $θ=arctan(y/x)$. The equation numbers are listed for where these functions occur later.

1. A. $Φ=Q2πlnr$, (3.1)

2. B. $Φ=−Γ2πθ$, (3.3)

3. C. $Ψ=Γ2πlnr$, (3.3)

4. D. $Ψ=Q2πθ$, (3.1)

5. E. $Φ=−r24$, (3.114a)

6. F. $Ψ=r24$, (3.114b)

7. G. , (4.26)

8. H. , (4.26)

9. I. , (4.26)

10. J. , (4.26)

11. K. , (4.26)

12. L. , (4.26)

13. M. , (4.26)

14. N. , (4.26)

## 1.2.1 Groundwater Flow: Head and Discharge

Studies of groundwater flow examine the movement of water through geological media, and the important properties and processes are illustrated in Fig. 1.6. Groundwater enters the flow domain as the portion of precipitation that seeps downward through the vadose zone, where pores between soil particles are partially filled with water and partially with air, to the groundwater zone, where porous and fractured media become fully filled with water. Destinations for groundwater include water pumped by wells, baseflow from groundwater to surface water, and phreatophyte root uptake by plants that can directly tap groundwater. The geologic media are separated into layers through which water moves readily called aquifers (e.g., coarse sands and gravels), and aquifer layers may be separated by aquitards with low flow rates (e.g., clays and shales), as illustrated by the confining layer in Fig. 1.6. It is important to know whether an unsaturated zone with unconfined conditions exists at its top or whether an aquifer layer is completely saturated with groundwater and is confined, and these conditions are quantified by the

$Display mathematics$
(1.15)

where B is the base elevation and H is the vertical thickness of the aquifer layer. Different approaches are employed to formulate the equations for saturated groundwater flow in (p.10) this section (where horizontal flow predominates) and the unsaturated vadose zone flow in Section 1.2.2 (where vertical flow predominates).

Figure 1.6 Groundwater properties and processes shown in section view.

Solutions to groundwater flow problems are presented in terms of head, h, which is the level to which water will rise in a pipe and is a measure of potential energy. Groundwater moves from regions with high head to low head, which is represented in Darcy’s law (Darcy, 1856) as a constitutive relation between the specific discharge, (qx, qy, qz), and changes in head:

$Display mathematics$
(1.16)

This equation uses the hydraulic conductivity K, which is a property of the fluid (density, ρ‎, and absolute viscosity, μ‎) and the porous media (intrinsic permeability, κ‎), and g is the acceleration of gravity. The example solutions in Fig. 1.7 illustrate the distribution of groundwater level (head) and flow rates. The Dupuit assumption (Dupuit, 1863) approximates head and the horizontal components of flow as uniform in the vertical direction within an aquifer later, which enables the specific discharge to be vertically integrated across the saturated depth to give the two-dimensional discharge per unit width with components Qx and Qy:

$Display mathematics$
(1.17)

It is convenient to adopt a representation where equations take on a common form for both unconfined and confined flow. This is accomplished by substituting Darcy’s law, (1.16), into (p.11) the Dupuit form of discharge per width, (1.17), and formulating the problem using the discharge potential, Φ‎, from Strack (1989):

Figure 1.7 A regional model of groundwater, and approximating the regional flow for another problem as uniform in a local model. (Reprinted from Steward and Allen, 2013, The Analytic Element Method for rectangular gridded domains, benchmark comparisons and application to the High Plains Aquifer, Advances in Water Resources, Vol. 60, Fig. 5, with permission from Elsevier.)

$Display mathematics$
(1.18)

when K, B, and H are piecewise constant in an aquifer layer.

Solutions to groundwater problems must also satisfy conservation of mass, which, together with (1.18), leads to the Laplace equation, (2.2a), for divergence-free flows:

$Display mathematics$
(1.19)

The examples in Fig. 1.8 illustrate the use of discharge potentials that satisfy this equation to reproduce the drawdown associated with wells, and the interactions between groundwater and surface water bodies such as lakes. Boundary conditions are utilized in their construction by specifying the pumping rate of the wells, and by specifying the head h* at locations (xm, ym) along the boundary of a lake that is directly connected to the groundwater:

Figure 1.8 Groundwater interactions with wells and lakes.

(p.12)

$Display mathematics$
(1.20)

Solutions to the Laplace equation with a discharge vector that is both irrotational and divergence-free may be studied using a stream function, (1.14), that is constant along streamlines. These streamlines illustrate the zones where water is captured by a well, as well as those portions of the boundary of a lake fed by groundwater and those where surface water moves from a lake to groundwater. The straight white line to the left of each element is formed by branch cuts, which are discussed later in Section 3.1.

A primary source of groundwater is recharge from the terrestrial system. Conservation of mass requires the divergence of the discharge per width to equal the recharge rate R, which leads to the Poisson equation, (2.2b):

$Display mathematics$
(1.21)

The recharge illustrated in Fig. 1.9 is modeled using analytic elements that satisfy conservation of energy where the head (and the discharge potential) has the same values outside and inside the recharge area, and conservation of mass where the normal derivative in the $n$-direction is also continuous at locations (xm, ym) on the polygon’s boundary is

$Display mathematics$
(1.22)

The notation used to designate sides of interfaces outside an element (the + side) and inside an element (the − side) is further elaborated for the Riemann–Hilbert problem in Section 2.3.1. The uniform recharge within a polygon in Fig. 1.9a is extended to model the extraction of groundwater by a plant’s, roots over a circular area in Fig. 1.9b, and by many plants in Fig. 1.9c. This formulation is also extensible to studying the interactions occurring between aquifer layers interconnected by leaky aquitards (Strack, 1989, sec. 14).

The groundwater and surface water exchanges occurring through rivers are often approximated using line elements. A river may be directly connected to groundwater, with head hm specified at a set of M points located at (xm, ym) along the river,

Figure 1.9 Recharge to groundwater and root water uptake by plants (phreatophytes). (Reprinted from Steward and Ahring, 2009, An analytic solution for groundwater uptake by phreatophytes spanning spatial scales from plant to field to regional, Journal of Engineering Mathematics, Vol. 64(2), Figs. 2,5, used with permission from Springer Nature; distributed under Creative Commons license (Attribution-Noncommercial), https://creativecommons.org.)

(p.13)

$Display mathematics$
(1.23)

where the potential Φ‎m are obtained from hm using (1.18). This boundary condition is illustrated in Figs. 1.10a and 1.10b, where head varies linearly along each river segment. A river may also be separated from groundwater by riverbed sediments of hydraulic conductivity K* and vertical thickness H* that cause a difference in head between the surface water h* and the groundwater h. In this case, the discharge per length of the river flowing into the aquifer, $ΔQn$, is equal to minus the vertical component of flow through the riverbed times the width of the river w*,

$Display mathematics$
(1.24)

using Darcy’s law (1.16). The riverbed in Fig. 1.10c has sediments that satisfy boundary conditions at points m expressed in terms of the resistance of the river bed sediments, c*,

$Display mathematics$
(1.25)

where this head may be expressed in terms of the discharge potential, (1.18).

Groundwater flow is influenced by the geological medium through which water flows. Porous media may be grouped into heterogeneity zones where jumps in an aquifer layer may occur in the base elevation (B+ to B), thickness (H+ to H), and/or hydraulic conductivity (K+ to K), across the boundary of a zone. Heterogeneities may either be isolated (like the polygon elements in Figs. 1.11a and 1.11b) or form adjacent elements that tile a domain (like the rectangle elements in Fig. 1.11c). Each interface satisfies conservation of energy and conservation of mass, like the recharge zone in (1.22):

$Display mathematics$
(1.26a)

Figure 1.10 River networks with specified head or a riverbed with resistance. (Reprinted from Steward, 2015, Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resources Research, Vol. 51(11), Fig. 6, with permission from John Wiley and Sons. Copyright 2015 by the American Geophysical Union.)

Figure 1.11 Heterogeneity with the geometry of a polygon or interconnected rectangles. (Reprinted from Steward and Allen, 2013, The Analytic Element Method for rectangular gridded domains, benchmark comparisons and application to the High Plains Aquifer, Advances in Water Resources, Vol. 60, Fig. 3, with permission from Elsevier.)

(p.14) (p.15) Heterogeneities with a jump only in hydraulic conductivity lead to the following condition for the potential,

$Display mathematics$
(1.26b)

and conditions for discontinuous base and thickness follow from Steward (2007). The contours of head and the vector field in Fig. 1.11 illustrate the continuous head and normal component of flow across the interfaces of heterogeneities.

Thin geologic features with a small cross-sectional width w* and properties different than those in the surrounding media form heterogeneities that may be approximated using line elements. A thin heterogeneity with higher conductivity K* than K in the surrounding media provides a conduit through which groundwater travels more readily than in the surrounding media, with examples in Fig. 1.12. Such an element has a continuous potential across the element, and the discharge per width in the $s$-direction tangential to the heterogeneity may be expressed in terms of the stream function or the partial derivative of the potential by

$Display mathematics$
(1.27a)

In the limiting case of an infinitely conductive heterogeneity, the potential is uniform along the element:

$Display mathematics$
(1.27b)

Analytic elements may also be developed to satisfy the boundary conditions of thin heterogeneities with lower conductivity than the surrounding media. The flow through (p.16) such a leaky element is related to the jump in discharge potential across the element, and the stream function is continuous across the element:

Figure 1.12 Thin groundwater features with high conductivity. (Reprinted from Steward, 2015, Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resources Research, Vol. 51(11), Fig. 4, with permission from John Wiley and Sons. Copyright 2015 by the American Geophysical Union.)

$Display mathematics$
(1.28a)

In the limiting case of zero conductivity, there is no normal component of flow, and the stream function is uniform along the element

$Display mathematics$
(1.28b)

Such thin geologic features provide a barrier that impedes groundwater flow, as illustrated in Fig. 1.13.

Analytic elements may also be developed for aquifer features that generate three-dimensional flow. The specific discharge vector for such problems satisfy Darcy’s law, (1.16), which is represented here in terms of a three-dimensional specific discharge potential Φ‎(x, y, z), and this together with conservation of mass gives the Laplace equation:

$Display mathematics$
(1.29)

The average groundwater velocity v is related to the specific discharge by

$Display mathematics$
(1.30)

Figure 1.13 Thin groundwater features with low conductivity. (Reprinted from Steward, 2015, Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resources Research, Vol. 51(11), Fig. 5, with permission from John Wiley and Sons. Copyright 2015 by the American Geophysical Union.)

Figure 1.14 Three-dimensional streamlines along stream surfaces connecting points A, B, C, and D near fully penetrating and partially penetrating wells. (Reprinted from Steward, 1998, Stream surfaces in two-dimensional and three-dimensional divergence-free flows, Water Resources Research, Vol. 34(5), Figs. 3,5, with permission from John Wiley and Sons. Copyright 1998 by the American Geophysical Union.)

(p.17) where the porosity, n, is a ratio of the volume of voids between soil particles divided by the total volume. The examples in Fig. 1.14 illustrate three-dimensional pathlines existing around wells (Steward and Jin, 2001), which are obtained by tracking particles that originate along the straight lines connecting points A, B, C, and D and travel with the groundwater velocity v.

Three-dimensional pathlines occur as groundwater travels through and around heterogeneities in geologic media. This is illustrated in Fig. 1.15 for inclusions, where the conductivity inside the three-dimensional element is higher than that of the surrounding soil. Groundwater becomes transported faster in such heterogeneity, resulting in a lateral shift in the relative position of water particles as observed by deformations to the geometry of upgradient and downgradient stream surfaces. This phenomena is called advective mixing and contributes towards transverse dispersion in contaminant transport (Janković et al., 2009).

A vertical component of flow for quasi three-dimensional flow, qz, is also generated for the two-dimensional formulation of the discharge potential, (1.18), due to recharge R occurring into the top of an unconfined aquifer and due to the gradient of head (Anderson, 2005). This component is obtained by vertically integrating the three-dimensional conservation of mass equation (1.29) from the base of an aquifer layer to the point z where qz is to be evaluated (Strack, 1984)

$Display mathematics$
(1.31a)

where the Dupuit assumption (1.17) requires qx and qy to remain constant in the z-direction. These partial derivatives are expressed in terms of the discharge per width using (1.17)

$Display mathematics$
(1.31b)

Figure 1.15 Three-dimensional stream surfaces near a highly conductive heterogeneity. (Reprinted from Janković, Steward, Barnes, and Dagan, 2009, Is transverse macrodispersivity in three-dimensional transport through isotropic heterogeneous formation equal to zero? A counterexample, Water Resources Research, Vol. 45(8), Fig. 2, with permission from John Wiley and Sons. Copyright 2009 by the American Geophysical Union.)

(p.18) These terms may be simplified using conservation of mass, (1.19), and rearranging the relations between discharge potential, head and discharge per width (1.18):

$Display mathematics$
(1.31c)

Together, the equations in (1.31) give the vertical component of flow:

$Display mathematics$
(1.32)

The three-dimensional streamlines associated with groundwater uptake by phreatophytes are illustrated in Fig. 1.16.

Example 1.1 This example illustrates computation of head and flow in unconfined and confined aquifer layers for the flow in Fig. 1.17. with a thin aquitard separating zones 2 and 3. The discharge potential at the boundaries of the four zones in Fig. 1.17 may be obtained from the specified heads hA and hD and the computed heads:

$Display mathematics$

Within each zone, the general solution to the Poisson equation, (1.21), is

$Display mathematics$

(p.19) where the coefficients c0 and c1 are adjusted to give the correct value of discharge potential and head at the edges of each zone. This gives the following for zone 1 ($0≤x and $B1≤z≤h$)

Figure 1.16 Quasi three-dimensional capture zones for root water uptake from groundwater by a field of trees, tracing pathlines backwards from tree locations. (Reprinted from Steward and Ahring, 2009, An analytic solution for groundwater uptake by phreatophytes spanning spatial scales from plant to field to regional, Journal of Engineering Mathematics, Vol. 64(2), Figs. 6,7, used with permission from Springer Nature; distributed under Creative Commons license (Attribution-Noncommercial), https://creativecommons.org.)

$Display mathematics$

for zone 2 ($L1≤x and $B2≤z≤h$)

$Display mathematics$

for zone 3 ($L1≤x and $B3≤z≤H3$)

$Display mathematics$

and for zone 4 ($L1+L2≤x≤L1+L2+L4$ and $B4≤z≤h$)

$Display mathematics$

similar to Strack (1981).

Figure 1.17 Dupuit groundwater flow in a vertical plane with an impermeable aquitard, and both unconfined (zone 1, 2, and 4) and confined (zone 3) conditions.

(p.20) This groundwater head and velocity is illustrated in Fig. 1.17b for the variables , , , , , , $n=0.25$, , and . Note that recharge R enters the top of zones 1, 2, and 4, and no recharge seeps through the confining layer into zone 3. Conservation of energy is explicitly satisfied by requiring that the head be equal to hB and hC across respective zones. Conservation of mass may be verified at $x=L1$ by computing the discharge per width on the left side in zone 1 and on the right side in zones 2 and 3:

$Display mathematics$

Likewise, continuity of flow may be verified at $x=L1+L2$ on the left and right sides of the interface

$Display mathematics$

The values of Φ‎ and Qx may be computed at any point to provide solutions for head (1.18), specific discharge (1.17) and (1.32), and velocity (1.30).

Problem 1.2 Compute the head (h in m) and the horizontal and vertical components of the velocity (vx, vz in m/day) for the aquifer in example 1.1 at the following locations:

1. A. and $z=h$ (zone 1)

2. B. and $z=h/2$ (zone 1)

3. C. and $z=h$ (zone 1)

4. D. and $z=h$ (zone 1)

5. E. and $z=h$ (zone 2)

6. F. and $z=h$ (zone 2)

7. G. and $z=h$ (zone 2)

8. H. and $z=H3$ (zone 3)

9. I. and $z=H3$ (zone 3)

10. J. and $z=H3$ (zone 3)

11. K. and $z=h$ (zone 4)

12. L. and $z=h$ (zone 4)

13. M. and $z=h$ (zone 4)

14. N. and $z=h/2$ (zone 4)

## (p.21) 1.2.2 Vadose Zone Flow: Pressure and Seepage Velocity

Studies of vadose zone flow examine the movement of water from the land surface to groundwater as depicted earlier in Fig. 1.6. In the vadose zone, water seeps downward through pores between soil particles that are partially filled with water and partially filled with air. Solutions are presented in terms of the pressure head p from Bernoulli’s equation (Bernoulli, 1738)

$Display mathematics$
(1.33)

where head is partitioned into pressure, P, and elevation terms with the z-axis directed upward against gravity, and the kinetic velocity head $|v|2/(2g)$ is ignored since it is small. For example, the pressure head distribution is shown in Fig. 1.18 for three different soils experiencing the same rate of recharge above a saturated groundwater table. While $p=0$ at $z=0$, it decreases asymptotically at higher elevations towards a background suction pressure head (Raats and Gardner, 1974).

The governing equations for vadose zone processes contain non-linear relations since water moves more readily through moist soil (higher p) than dry soil (lower p). This is observed in Darcy’s law, (1.16),

$Display mathematics$
(1.34)

with hydraulic conductivity expressed as a function of pressure head. The Richards equation (Richards, 1931) is obtained by substituting this form of Darcy’s law into conservation of mass

$Display mathematics$
(1.35)

Figure 1.18 Pressure head distribution above saturated groundwater.

Table 1.1 Representative saturated hydraulic conductivity Ks and sorptive number a for the Gardner equation, (1.38), matric flux potential at saturation Fs, (1.40), and modified Helmholtz coefficient k, (1.39) from Steward (2016, table 1).

Soil type

Ks [m/day]

a [1/m]

Fs [m2/day]

k [1/m]

Silt

0.02

2

0.01

1

Fine sand

1

5

0.2

2.5

Coarse sand

35

25

1.4

12.5

(p.22) This equation may be rearranged using the matric flux potential, denoted F, defined by the Kirchhoff transformation (Kirchhoff, 1894) as

$Display mathematics$
(1.36)

to provide a linear form of the Richards equation

$Display mathematics$
(1.37)

where the last term was rearranged using the chain rule for differentiation $∂K∂z=dKdp∂p∂z$, with $∂p∂z$ from (1.36).

A solution to flow in the vadose zone requires specification of the soil properties relating the hydraulic conductivity K(p) to the pressure head. While many functional forms exist, it is convenient to use the Gardner equation (Gardner, 1958), which provides the “quasi-linear” formulation of the Richards equation (Philip, 1968; Pullan, 1990),

$Display mathematics$
(1.38)

where a is the sorptive number, and the saturated conductivity Ks occurs at $p=0$. The first derivative on the right-hand side of the last equation may be rearranged by changing variables to a potential function Φ‎ and constant k (Wooding, 1968; Knight et al., 1989), giving the modified Helmholtz equation, (2.2d):

$Display mathematics$
(1.39)

The Gardner equation (1.38) also provides a relation between the matric flux potential and pressure head by integrating the Kirchhoff transformation (1.36):

$Display mathematics$
(1.40)

where Fs is the matric flux potential at $p=0$. The representative values for soil properties in Table 1.1 were used to construct Fig. 1.18.

Figure 1.19 Typical vadose zone problems of seepage through horizontal soil layers and through embedded heterogeneities.

(p.23) While solutions to vadose zone problems are formulated in terms of the potential Φ‎, solutions are represented in terms of pressure head h and flow. Expressions relating pressure head and potential are obtained by equating the expressions for F in the last two equations, and rearranging terms as,

$Display mathematics$
(1.41)

Likewise, expressions that provide the specific discharge from Φ‎ are obtained from Darcy’s law (1.34) with the Kirchhoff transformation (1.36),

$Display mathematics$
(1.42)

using derivative of the matric flux, (1.39).

Soils may be placed in horizons with layered soils of different types as depicted in Fig. 1.19a. Such problems require the pressure head to be continuous across the interface of adjacent soil layers, and solutions are shown in Fig. 1.20 for horizontal layers of alternating soil types with the same recharge rate as the single soil in Fig. 1.18. The continuous pressure head and vertical component of flow between soil layers is explored shortly in Example 1.2.

Soils may also contain heterogeneities with one soil type with properties $Ks−$ and a embedded within a uniform soil with properties $Ks+$, a+ as illustrated in Fig. 1.19b. Such interfaces satisfy conservation of energy, with the same value of pressure head outside and inside the heterogeneity, which leads to a non-linear boundary condition

$Display mathematics$
(1.43)

(p.24) in terms of the potential function outside and inside the heterogeneity, Φ‎+ and Φ‎, using (1.41). Conservation of mass is satisfied if the normal component of the specific discharge is continuous across the boundary

Figure 1.20 Pressure head distribution in soil layers with interface conditions.

$Display mathematics$
(1.44)

which is expressed in terms of the potential using (1.42). Solutions to these two interface conditions are illustrated for heterogeneities formed by silt embedded within a fine sand in Fig. 1.21, and for a coarse sand embedded within a fine sand in Fig. 1.22. Note that (p.25) the recharge flowing through these layers of embedded heterogeneities is the same as that flowing through the continuous horizontal layers in Fig. 1.20.

Figure 1.21 Vadose zone seepage through silt heterogeneities embedded in a fine sand. (Reprinted from Steward, 2016, Analysis of vadose zone inhomogeneity toward distinguishing recharge rates: Solving the nonlinear interface problem with Newton method, Water Resources Research, Vol. 52(11), Figs. 3,4,5, with permission from John Wiley and Sons. Copyright 2016 by the American Geophysical Union.)

Figure 1.22 Vadose zone seepage through coarse sand heterogeneities embedded in a fine sand. (Reprinted from Steward, 2016, Analysis of vadose zone inhomogeneity toward distinguishing recharge rates:Solving the nonlinear interface problem with Newton method, Water Resources Research, Vol. 52(11), Figs. 3,4,5, with permission from John Wiley and Sons. Copyright 2016 by the American Geophysical Union.)

An example is presented next to illustrate how the equations for the vadose zone are used to compute pressure and velocity, as well as application of the conditions of conservation of energy and mass across the interfaces between soils with different properties. The velocity is computed from the specific discharge in the vadose zone by

$Display mathematics$
(1.45)

where the volumetric moisture content θ‎ is a ratio of the volume of pores filled with water to the total volume. When the soil moisture content θ‎(p) is approximated using the same form as K(p) in the Gardner equation, this leads to an expression for the specific moisture capacity

$Display mathematics$
(1.46)

where θ‎r is the residual (irreducible) volumetric water content and θ‎s is the value at saturation (Warrick, 1974; Basha, 1999).

Example 1.2 Vertical seepage of water through a vadose zone with horizontal soil layers is shown in Fig. 1.20. This solution is obtained from separation of variables, (4.19),

$Display mathematics$

with coefficients, (4.20), obtained by matching specified values of pressure head and specific discharge at elevation z0:

(p.26)

$Display mathematics$

In the limit as z increases in each layer the pressure approaches

$Display mathematics$
(1.47)

The pressure head distribution in a stacked set of layers may be computed successively following Raats and Gardner (1974), beginning by computing c1 and c2 at the the bottom of the lowest soil layer, where $p=0$ and $z=0$. The pressure head from this layer is computed at the interface with the next higher layer, which gives the coefficients c1 and c2 in that layer, and so on. The solution in Fig. 1.20a was constructed for recharge . The pressure head at the interfaces are:

$Display mathematics$

Problem 1.3 Compute the pressure head at elevations zA, zB, zC, and zD for Example 1.2, with the following recharge rates:

1. A.

2. B.

3. C.

4. D.

5. E.

6. F.

7. G.

8. H.

9. I.

10. J.

11. K.

12. L.

Figure 1.23 Incompressible fluid flow properties and processes: Wing shown in section view.

Figure 1.24 Fluid flow through a duct comprised of interconnected rectangles and around curvilinear elements. (Reprinted from Steward, Le Grand, Janković, and Strack, 2008, Analytic formulation of Cauchy integrals for boundaries with curvilinear geometry, Proceedings of The Royal Society of London, Series A, Vol. 464, Fig. 3, used with permission of The Royal Society; permission conveyed through Copyright Clearance Center, Inc.)

## (p.27) 1.2.3 Incompressible Fluid Flow: Pressure and Velocity

The study of fluid flowing around an element with velocity v is illustrated in Fig. 1.23. A wing imposes boundary conditions along the element, where the the normal component $vn$ is zero and fluid flows parallel to its surface. Wings are designed so the fluid moves at with a higher magnitude of velocity $|v1|$ at point 1 along the top surface than $|v2|2$ at point 2 along the bottom. This creates a lower pressure P along the top that may be quantified using the Bernoulli (1738) equation

$Display mathematics$
(1.48)

where z is the elevation and the density ρ‎ is uniform for incompressible fluids. The difference in pressure generates a lift force that is easily observed by placing your hand in a wind and varying this force by cupping ones hand and changing the attack angle relative to the background air velocity. Lift forces exist in slender bodies such as airfoils, hydrofoils, and helicopters; and are used in the motion of birds and fish.

Incompressible fluid flow with constant density may be studied using a velocity potential Φ‎ (Lamb, 1879, sec.23)

$Display mathematics$
(1.49)

For example, fluid flowing through an impermeable duct is shown in Fig. 1.24a where the fluid particles travel from higher to lower values of this velocity potential. Such flows satisfy conservation of mass and may be formulated using the Laplace equation, (2.2):

$Display mathematics$
(1.50)

Figure 1.25 Streamlines in half-plane and three-dimensional breakthrough planes near an impermeable element.

(p.28) Two-dimensional problems with incompressible fluid flow may also be formulated using the Lagrange stream function (Lamb, 1879, sec.69), presented earlier in (1.10)

$Display mathematics$
(1.51)

Such solutions are illustrated by streamlines with constant Ψ‎ in Figs. 1.24b and 1.24c for fluid flow around obstructions on a flat surface and potential flow in cavities.

Analytic elements may be formulated to study the three-dimensional flow of incompressible fluid around impermeable elements as illustrated in Fig. 1.25. The upper panels show the velocity potential Φ‎ and vector field around a sphere, a prolate spheroid, and an oblate spheroid; where each three-dimensional element is formed by rotation of a circle or ellipse about the z-axis. These axisymmetric flows may be represented by streamlines with uniform Stokes (1842) stream function, obtained for components of the vector field in cylindrical (r, θ‎, z) coordinates using

$Display mathematics$
(1.52)

The three-dimensional visualization in the lower panels of this figure illustrate stream surfaces formed by tracing the trajectories of fluid particles released along the lines at the left-most edge. These pathlines are traced to the right-most edge and stream surfaces are illustrated by connecting the locations where neighboring streamlines intersect breakthrough planes at the middle and edges of the displayed region.

Figure 1.26 Point vortexes and circle and line elements with circulations.

(p.29) The velocity potential and stream function may be used to study fluid flow with circulation as shown in Fig. 1.26. These examples illustrate circulation about isolated points, about a infinitely long cylinder perpendicular to this plane, and about a two-dimensional line element. For elements that generate circulation, a discontinuity in the velocity potential occurs along a branch cut, which has been aligned to occur on the left side of these elements. A branch cut is related to the net circulation Γ‎ generated by an element and will be analyzed in detail beginning in Section 3.1. The use of circulation to study flight is described next.

The geometry of a wing that is commonly studied was developed by Joukowsky (1912), and is illustrated in Fig. 1.27. The vector field in Fig. 1.27a satisfies the far-field condition of uniform vector field approaching the wing, and the impermeable condition that the vector field must be aligned tangent to the wing. Another vector field is illustrated in Fig. 1.27b where flow circulates about an impermeable wing, but with no flow at large distances from the element. However, both solutions contain a singularity at the sharp trailing edge where the vector field becomes infinite. The Kutta condition requires flow to remain finite at the trailing edge, and this condition is satisfied when the wing generates exactly the correct circulation so that the two singularities cancel. This gives the fluid flow in Fig. 1.27c where the far-field velocity passes around a wing that satisfies this Kutta condition.

Examples of flow about wings in Fig. 1.28 illustrate the capacity to study fluid flow around flat airfoils, curved wings, and the Joukowsky airfoil. While some flows may be described mathematically using a single analytic element, more complicated geometries required the wing to be discretized into a set of adjacent elements. Development of the AEM for solving these problems is founded on the panel method for flow about airplane wings developed by Hess and Smith (1967) and Hess (1990).

Example 1.3 The vector field for flow about the rotating cylinder in Fig. 1.26b is given by

$Display mathematics$
(1.53)

Figure 1.27 Flight boundary conditions: $vn=0$ and Kutta condition, from Section 3.5.

(p.30) where vx0 is the background vector field in the x-direction and Γ‎ is the circulation, adapted from (3.3) and (3.7). Note that the vertical z-axis is uniform for this problem and pressure is related to the fluid velocity by Bernoulli’s equation (1.48).

Problem 1.4 Compute the difference in pressure in units of N/m2 between points at the top and bottom of a rotating cylinder in water for the specified background velocity vx0 and circulation Γ‎:

1. A. , $Γ=0.25π/s$

2. B. , $Γ=0.25π/s$

3. C. , $Γ=0.25π/s$

4. D. , $Γ=0.5π/s$

5. E. , $Γ=0.5π/s$

6. F. , $Γ=0.5π/s$

7. G. , $Γ=π/s$

8. H. , $Γ=π/s$

9. I. , $Γ=π/s$

Figure 1.28 Wings with Kutta condition at left ends (shown in x-y plane).

## (p.31) 1.2.4 Thermal Conduction: Temperature and Heat Flux

Studies of thermal conduction examine the propagation of heat through media as illustrated in Fig. 1.29. Temperature gradients are generated as heat propagates from sources to sinks, and these processes satisfy Fourier’s law (Fourier, 1878):

$Display mathematics$
(1.54)

where the thermal conduction q, which is rate of heat flow per unit time per unit area, varies linearly with the gradient of temperature T, and K is the thermal conductivity. Problems with piecewise uniform thermal conductivity may be formulated in terms of the thermal potential Φ‎. Such problems are governed by the heat equation, which results from the conservation law for heat conduction together with Fourier’s law (Carslaw and Jaeger, 1959):

$Display mathematics$
(1.55)

where ρ‎ is the density, c is the specific heat, and $D$ is the thermal diffusivity. For steady heat conduction, the transient term on the right-hand side of the heat equation is zero, leading to the Laplace equation, (2.2a):

$Display mathematics$
(1.56)

Use of the AEM to formulate solutions for thermal conduction builds upon foundational studies by Carslaw and Jaeger (1959), Moon and Spencer (1961a), and MacRobert (1967).

Figure 1.29 Thermal conduction properties and processes shown in section view.

Figure 1.30 Thermal conduction: temperature and flux near a cooling rod.

(p.32) The temporal process of heat condition for a heat sink is illustrated in Fig. 1.30. For this problem, an infinitely long cooling rod oriented normal to this picture begins to withdraw heat from the domain from initial conditions of uniform temperature as shown in Fig. 1.30a. The temperature becomes lower in the vicinity of the sink in Fig. 1.30b, and at larger times the thermal condition q approaches the final steady-state conditions in Fig. 1.30c. This is the same problem as groundwater flow toward a pumping well studied by Theis (1935), which approaches the final steady-state conditions from Thiem (1906). The transition from initial to final steady states is explored later in Example 1.4 and Problem 1.5.

Extensions to problems for heat generation along elements with more complex geometries are illustrated next. For example, heat is conducted from four rods in a cylinder to its boundary of uniform temperature in Fig. 1.31a. Such solutions are extensible to problems where heat is generated across the domain using functions that satisfy the Poisson equation (Moon and Spencer, 1961a, p.427):

$Display mathematics$
(1.57)

where R is the rate of heat production. Such solutions may be used to represent, for example, propagation of the heat of hydration within a concrete cylinder inserted into a water bath. Heat production is also illustrated for elements with the geometry of a Bezier curve in Fig. 1.31b or a B-spline in Fig. 1.31c.

Analytic elements may be used to study three-dimensional problems with heat production, as illustrated in Fig. 1.32, where each plate is bounded on the top and bottom surfaces by insulators that hold heat within the domain. A heat source is embedded within the domain in Fig. 1.32a, and the three-dimensional trajectories of thermal conduction illustrate pathways as heat moves away from the source to further regions within the plate. Trajectories of thermal conduction are also illustrated for elements that apply heat to the boundary of the plate along a line in Fig. 1.32b, and across an area on the top boundary in Fig. 1.32c.

Figure 1.31 Thermal conductors with uniform temperature. (Reprinted from Steward, Le Grand, Janković, and Strack, 2008, Analytic formulation of Cauchy integrals for boundaries with curvilinear geometry, Proceedings of the Royal Society of London, Series A, Vol. 464, Fig. 4, used with permission of The Royal Society; permission conveyed through Copyright Clearance Center, Inc.)

(p.33) Thermal conduction around thin insulators is illustrated in Fig. 1.33 for both isolated elements and strings of connected elements. A perfect insulator through which heat cannot penetrate satisfies a boundary condition where the normal component of conduction is zero:

$Display mathematics$
(1.58)

Figure 1.32 Three-dimensional thermal conduction near a source or sink in a plate. (Reprinted from Steward, 1999, Three-dimensional analysis of the capture of contaminated leachate by fully penetrating, partially penetrating, and horizontal wells, Water Resources Research, Vol. 35(2), Fig. 2, with permission from John Wiley and Sons. Copyright 1999 by the American Geophysical Union.)

Figure 1.33 Thermal conduction around impermeable lines and interconnected line segments. (Reprinted from Steward, 2015, Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resources Research, Vol. 51(11), Fig. 5, with permission from John Wiley and Sons. Copyright 2015 by the American Geophysical Union.)

(p.34) Such a boundary also has a continuous stream function along the element, which is evident in the Fig. 1.33 solutions. Analytic elements will also be formulated for thin features with lower thermal conductivity, K*, than tthe surrounding media, K. Such elements satisfy a boundary condition where the normal component of conduction is related to the partial derivative of the thermal potential,

$Display mathematics$
(1.59)

which is approximated as the jump in thermal potential across the thin feature of width w*. This gives a Robin boundary condition,

$Display mathematics$
(1.60)

and in the limit as $δ→∞$ the normal component become zero, (1.58).

Thermal conduction is influenced by changes in the media through which heat propagates. For example, the circular element in Fig. 1.34a is impermeable to heat, and the object satisfies boundary conditions of a perfect insulator, (1.58). Interface conditions may be formulated for inhomogeneities with a different thermal conductivity, K, than the background, K+ into which it is inserted. Such elements satisfy continuity of the normal component of the temperature gradient at each point m on the interface

$Display mathematics$
(1.61a)

Figure 1.34 Thermal conduction around elements with low thermal conduction in a uniform vector field (showing contour intervals with uniform temperature).

(p.35) a condition that is also satisfied if the stream function is continuous at adjacent points across the interface. Additionally, the temperature must be continuous, which provides conditions for the potential function

$Display mathematics$
(1.61b)

Thermal conduction through and around a circular or elliptical element with lower thermal conductivity than the background is illustrated in Fig. 1.34.

These methods for solving problems with an inhomogeneity embedded into a uniform background may be extended to problems involving many elements. For example, the thermal conduction near circular elements is shown in Fig. 1.35a, and elliptical elements in Fig. 1.35b. These examples illustrate inhomogeneities with randomly distributed thermal conductivity, as evident by thermal trajectories diverging around elements with lower conductivity and converging within elements with higher conductivity. The interface conditions for inhomogeneities, (1.61), may also be applied to interconnected regions with different thermal conductivity as illustrated in Fig. 1.35c. These formulations facilitate studies of heterogeneous media with two different conceptualizations: embedded inhomogeneities within a uniform background, or a domain filled with interconnected regions each with unique thermal properties.

Analytic elements may also be formulated to study three-dimensional conduction near inhomogeneities. Such solutions are depicted in Fig. 1.36 by the trajectories of thermal conduction near spheres. Each element satisfies the interface conditions for an inhomogeneity, (1.61), and results are shown as the thermal conductivity varies from conditions of impermeable to infinitely conductive.

Example 1.4 The thermal potential and conduction for a point-sink in Fig. 1.30 that removes heat flux Q may be computed at distance r from the sink at time t using

Figure 1.35 Thermal conduction through heterogeneity with elements of randomly distributed conductivity embedded within a uniform background or with interconnected edges.

(p.36)

$Display mathematics$
(1.62)

where $Φ0=KT0$ represents the thermal potential at the initial condition, and the exponential integral, E1(u), and may be evaluated using Abramowitz and Stegun (1972, eqns 5.1.53, 5.1.56). Final steady-state conditions are given by

$Display mathematics$
(1.63)

Figure 1.36 Thermal conduction near impermeable, highly conductive, or infinitely conductive spheres.

(p.37) and extensions to transient problems with more complicated geometries are developed by Steward (1991), Zaadnoordijk and Strack (1993), Furman and Neuman (2003), Kuhlman and Neuman (2009), and Bakker (2013).

Problem 1.5 Compute thermal conduction, qr, toward a heat sink with unit flux $Q=1watt/m$ and thermal diffusivity at distance r from a point-sink in the specified material at initial conditions $t=0$, at , at , and at the final steady-state conditions. Thermal diffusivity values from (Moon and Spencer, 1961a, p.213).

1. A. copper, ,

2. B. copper, ,

3. C. cast iron, ,

4. D. cast iron, ,

5. E. concrete, ,

6. F. concrete, ,

7. G. dry soil, ,

8. H. dry soil, ,

## 1.2.5 Electrostatics: Voltage and Electric Fields

Studies of electric fields examine the propagation of electric currents through a conducting medium as illustrated in Fig. 1.37. Electrostatics may be formulated in terms of the electric field strength E, which is obtained from minus the gradient of the electric potential Φ‎. This, together with conservation of charge, provides solutions in terms of the Poisson equation (Sommerfeld, 1952, p.38):

$Display mathematics$
(1.64)

Figure 1.37 Electric field properties and processes.

(p.38) where ρ‎ is charge density, and ε‎ is the dielectric constant. Problems where steady charges propagate through the domain from sinks and sources along boundaries reduce to solutions of the Laplace equation, (2.2a): (Moon and Spencer, 1960)

$Display mathematics$
(1.65)

This formulation is grounded in Ohm’s law (Ohm, 1827), where the current I is equal to voltage V divided by resistance R:

$Display mathematics$
(1.66)

and these quantities are related to E and Φ‎ as follows. The current through a surface A may be obtained from the surface integral of electric current density $J=κE$, where κ‎ is the electric conductivity (Moon and Spencer, 1960, p.87). The voltage between points 1 and 2 is equal to the differences in electric potential at these points (Sommerfeld, 1952, p.39). The following solutions for problems in electrostatics build upon examples presented in the foundational works of Sommerfeld (1952), Moon and Spencer (1960, 1961a,b), and MacRobert (1967).

Electric currents generated by charged elements are illustrated in Fig. 1.38. The solution in Fig. 1.38a illustrates an electric field E directed towards a long charged wire, where the electric potential Φ‎ decreases in the direction of this sink. This solution may be extended to the problem of two charged rods with opposite charge in Fig. 1.38b, and the field is oriented from the positively to negatively charged rods. Each rod satisfies a boundary condition along the interface between the rod and the problem domain where the potential is uniform along its boundary

$Display mathematics$
(1.67)

In the limit as two charged rods with opposite strength become very thin and approach one another, the field takes the form of a point-dipole in Fig. 1.38c.

Figure 1.38 Sources and sinks of electric current, and a dipole.

Figure 1.39 Electric conduction through infinitely conductive and highly conductive circle and ellipse.

(p.39) The point-dipole is used in developing solutions for electric fields near inhomogeneities, where elements with a dielectric constant ε‎ are inserted into a uniform background with a different dielectric constant ε‎+. For example, the solution in Fig. 1.39a illustrates conduction through an infinitely conductive element in the limit as $ε−→∞$. Solutions are also shown for highly conductive elements for a circle in Fig. 1.39b (Sommerfeld, 1952, p.62) and for an ellipse in Fig. 1.39c (Moon and Spencer, 1961a, p.259). These elements have a larger dielectric constant inside the element than outside, $ε−>ε+$, but not large enough to approach the limiting value in Fig. 1.39a. These inhomogeneities satisfy conditions where voltage is continuous across its interface, so the electric potential is also continuous at point m on the interface:

$Display mathematics$
(1.68)

They also satisfy a condition where the normal component of the displacement, D, is continuous across the interface:

$Display mathematics$
(1.69)

This, together with the relation between the electric field and the gradient of the potential, (1.64), requires the normal component of the gradient of Φ‎ to jump across the interface (Sommerfeld, 1952, p.60). Note that the solutions in Fig. 1.39 visualize surfaces of constant Φ‎ and the displacement vector.

The electric fields near thin conducting elements with the geometry of a two-dimensional line are illustrated in Fig. 1.40. Each element satisfies a boundary condition of uniform potential, (1.68), similar to the infinitely conductive circle in Fig. 1.39a. This condition also requires that the electric field be oriented normal to each element and its tangential component is zero along the boundary

$Display mathematics$
(1.70)

Figure 1.40 Electric conduction through thin conductors and interconnected line segments. (Reprinted from Steward, 2015, Analysis of discontinuities across thin inhomogeneities, groundwater/surface water interactions in river networks, and circulation about slender bodies using slit elements in the Analytic Element Method, Water Resources Research, Vol. 51(11), Fig. 4, with permission from John Wiley and Sons. Copyright 2015 by the American Geophysical Union.)

(p.40) Examples are shown for an isolated element, for connected elements that transmit current along a string of elements, and for a group of non-connecting elements that collectively shape the electric field. It is evident that the potential is uniform along the thin conductors in each solution, and the vector field is directed normal to their boundaries.

The electrostatic formulation for studying steady electric fields and potential is extended to three-dimensional problems in Fig. 1.41. These examples illustrate inhomogeneities where the dielectric constant is larger inside the element than in the surrounding media, $ε−>ε+$, similar to the two-dimensional solutions in Fig. 1.39. The electric field near the sphere in Fig. 1.41a follows Moon and Spencer (1961a, pp.224–8), and solutions illustrate the pathways for current in half-planes. The intersection of these pathways with breakthrough planes is also shown for three planes, one at the left side of the review area, one at the right, and one midway between these. The electric field near a prolate spheroid in Fig. 1.41b follows (Moon and Spencer, 1961a, p.244–59), and has been used to study conduction through a grounding rod. The electric field is also illustrated in Fig. 1.41c for an object with the geometry of an oblate spheroid.

Example 1.5 The electric potential for the solution in Fig. 1.39a may be obtained from (3.6)

$Display mathematics$
(1.71)

where Ex0 is the background electric field in the x-direction, and the element is of radius and centered at the origin.

Figure 1.41 Electric fields near three-dimensional elements illustrating current pathways in half-planes and their intersections with breakthrough planes.

(p.41) Problem 1.6 Compute the electric field E for the field in Example 1.5 at the specified locations for unit $Ex0=1$ with .

1. A. $(x,y)=(−4,−1)cm,(−4,0)cm,and(−4,1)cm$

2. B. $(x,y)=(−3,−1)cm,(−3,0)cm,and(−3,1)cm$

3. C. $(x,y)=(−2,−1)cm,(−2,0)cm,and(−2,1)cm$

4. D. $(x,y)=(−1,−1)cm,(−1,0)cm,and(−1,1)cm$

5. E. $(x,y)=(1,−1)cm,(1,0)cm,and(1,1)cm$

6. F. $(x,y)=(2,−1)cm,(2,0)cm,and(2,1)cm$

7. G. $(x,y)=(3,−1)cm,(3,0)cm,and(3,1)cm$

8. H. $(x,y)=(4,−1)cm,(4,0)cm,and(4,1)cm$

# 1.3 Studies of Periodic Waves

The methods used to formulate solutions for wave fields and their interactions with analytic elements are presented next. The propagation of waves may be studied using the wave equation, (2.2g):

$Display mathematics$
(1.72)

where the constant c is related to the wave velocity (Airy, 1845; Lamb, 1879). The separation of variables subdivides the spatial and temporal variables in Φ‎ for periodic waves into functions that repeat with angular frequency ω‎ and a complex function φ‎ that satisfies the Helmholtz equation, (2.2c):

Figure 1.42 Background wave fields of plane waves and waves emanating from a point.

(p.42)

$Display mathematics$
(1.73)

where k is the wave number (Courant and Hilbert, 1962). This formulation utilizes complex variables with the imaginary number $i=−1$, and an introduction to complex functions with real and imaginary parts, $φ=ℜ(φ)+iℑ(φ)$, is provided in Section 2.3.2. This formulation of the Helmholtz equation is used to study periodic waves in acoustics (Morse and Ingard, 1968, p.728), electromagnetism (Sommerfeld, 1972, p.159), optics (Goodman, 1996, p.55), and surface water (Lamb, 1879, p.195).

The use of the complex function φ‎ for studying periodic wave fields is illustrated in Fig. 1.42 for the following examples:

$Display mathematics$
(1.74)

(p.43) and the exponential and Bessel functions are introduced later with the specified equations. These solutions are illustrated using lines of constant amplitude of the wave field, |φ‎|, and lines of constant cosine of the phase, cos[arg(φ‎)]. The plane waves in Fig. 1.42a have uniform amplitude |φ‎0|, and travel in the $θ0=30°$ direction. The waves originating from points in Figs. 1.42b and 1.42c have decreasing amplitude as they travel away from these sources. The wavelength, L, is related to the wave number and the wave period, T, is related to angular frequency by

$Display mathematics$
(1.75)

which provides an expression for the celerity, C, or wave velocity. The wavelength is evident in Fig. 1.42 by identifying locations with the same value of the cosine of the phase on neighboring waves.

Waves are transformed by their interactions at the boundaries and interfaces of elements. Wave reflection results as incoming waves approach a boundary in Fig. 1.43a and part of the energy is reflected back into the domain at a reduced amplitude factored by the reflection coefficient R. This leads to a Robin boundary condition

$Display mathematics$
(1.76)

where the coefficient α‎ is developed later as a function of the wave number k and the reflection coefficient R. Wave fields also also transformed by diffraction as waves bend around elements in Fig. 1.43b This requires the scattered waves to satisfy a radiation condition at large distances from the circular boundary condition (Sommerfeld, 1949, Section 28)

Figure 1.43 Wave interactions at boundaries and interfaces.

(p.44)

$Display mathematics$
(1.77)

so the scattered waves are outgoing and move away from the element. The wave field in Fig. 1.44 illustrates the reflected waves in front of the circle element and the calmer zone behind it. As the refection coefficient decreases, the wave absorption increases since less energy is reflected from the boundary.

A third mechanism for transformation of waves is refraction, where waves travel between regions with different wave numbers. This is illustrated in Fig. 1.43c for incoming waves that move across an interface normal to the $n$-direction. Two conditions must be satisfied at this interface: continuity of the complex wave function, and a condition relating the normal derivatives across the interface:

$Display mathematics$
(1.78)

Figure 1.44 Wave diffraction of plane waves around a circle element (Fig. 4.17). (Reprinted from Steward, 2018, Wave resonance and dissipation in collections of partially reflecting vertical cylinders, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 144(4), Fig. 1, with permission from ASCE.)

(p.45) The coefficients for this interface are developed later by satisfying Snell’s law

$Display mathematics$
(1.79)

relating the change in wave direction to the wave number on each side of the interface.

Problem 1.7 Compute the variable φ‎, its magnitude, |φ‎|, and the cosine of its phase, cos[arg(φ‎)] for the plane waves in Fig. 1.42a and Eq. (1.74) at the specified locations: , , and . This wave field is specified for , and the following direction θ‎0:

1. A. $θ0=0°$

2. B. $θ0=10°$

3. C. $θ0=20°$

4. D. $θ0=30°$

5. E. $θ0=40°$

6. F. $θ0=50°$

7. G. $θ0=60°$

8. H. $θ0=70°$

9. I. $θ0=80°$

10. J. $θ0=90°$

## 1.3.1 Water Waves: Amplitude and Phase

Waves are generated by winds blowing over an open sea that cause water particles in deep water to move in circular paths with the size of the orbits decreasing with depth, as illustrated in Fig. 1.45. These motions generate wave trains moving with a common period T in the same direction, and wave envelopes may travel hundreds to thousands of kilometers from their source. Waves are impacted by bathymetry as they move into coastal regions with shallower depth h, and the motion of water particles takes on the form of elliptical orbits. The wavelength L decreases and the amplitude A increases near the coastal boundary until the maximum slope of the water surface reaches $17$, at which point surface tension can no longer hold the wave together and it breaks and dissipates energy (Bascom, 1980).

The circular orbits of water particles are illustrated in Fig. 1.46 using plots of the water surface profile and water particle velocity shown at specific time intervals. Solutions are formulated by separating the spatial coordinates of the scalar potential in (1.73) into a real function Z of the z-coordinate and a complex function φ‎ of x and y:

Figure 1.45 Water waves properties and processes: The periodic motion of water particles generates a wave train.

Figure 1.46 Monochromatic waves in water of constant depth.

(p.46)

$Display mathematics$
(1.80)

and v is the water particle velocity. This gives solutions to the Helmholtz equation for φ‎, and an ordinary differential equation of Z

$Display mathematics$
(1.81)

where the solution for Z is given for waves of small amplitude in water of uniform depth (Dean and Dalrymple, 1984). Water waves satisfy a dispersion relation from Airy (1845, p. 290)

$Display mathematics$
(1.82)

which fixes the wave number k for a specific angular frequency and water depth. This gives expressions for the wave velocity (1.75), and the group velocity Cg, which is velocity at which energy gets propagated by waves (Dean and Dalrymple, 1984, p. 98)

$Display mathematics$
(1.83)

The hyperbolic functions in these relations approach asymptotic limits for shallow water (where $tanhkh→kh$ and $sinh2kh→2kh$) and deep water (where $tanhkh→1$ and $sinh2kh→e2kh/2$) that are summarized in Table 1.2.

Waves are transformed by reflection of wave energy by the coast back to the sea. The incoming plane waves in Fig. 1.43a are directed towards a coastline at angle γ‎, and it is natural to express their potential, (1.74), in terms of the tangential and normal coordinates (p.47) $s$-$n$ along the coast with incoming wave amplitude A0 proportional to φ‎0. These incoming waves are superimposed with reflected waves directed at angle $−γ$ away from the coast with amplitude A0R,

Table 1.2 Water wave properties for water depth that is shallow, intermediate, or deep

Variable

Shallow

Intermediate

Deep

Range of wavelength

$20h

$2h

$L<2h$

Range of wave number

$kh<π/10$

$π/10

$π

Dispersion relation, (1.82)

$ω2=k2gh$

$ω2=kgtanhkh$

$ω2=kg$

Wavelength, (1.75)

$L=2πghω$

$L=2πk$

$L=2πgω2$

Wave velocity, (1.83)

$C=gh$

$C=gktanhkh$

$C=gk$

Group velocity, (1.83)

$Cg=C$

$Cg=C2(1+2khsinh2kh)$

$Cg=C2$

$Display mathematics$
(1.84)

where the reflection coefficient R takes on values from $R=0$ for a coast that fully absorbs wave energy to $R=1$ for full reflection of wave energy, and β‎ is a phase shift between incoming and reflected waves at the boundary $n=0$ (Berkhoff, 1976). A coastal boundary condition may be developed by evaluating the partial derivative with respect to the normal direction on the boundary

$Display mathematics$
(1.85)

It is common to approximate the coastal boundary condition as

$Display mathematics$
(1.86)

with incoming waves normal to the coast and no phase shift (Steward and Panchang, 2000). Wave reflection is illustrated for a fully absorbing coast with $γ=30°$ in Fig. 1.47a, for partial wave reflection in Fig. 1.47b, and for a fully reflective coast in Fig. 1.47c. It is observed that standing waves form with increased coastal reflection and generate a sea with amplitude varying between 0 and twice that of the incoming waves, and phase shifts occur with full reflection across locations where $|φ|=0$.

Diffraction of water waves occur when the coastal boundary condition (1.86) is satisfied along elements, and the scattered waves travel away from the element, (1.77). For example, the envelope of incoming water waves that travel around a circular island in Fig. 1.44a is extended to multiple islands in Fig. 1.48. These elements interact to form a complicated sea with some regions experiencing wave amplification and others with reduction in amplitude (Steward, 2018). The boundary condition (1.86) may be separated into its amplitude and phase parts (Berkhoff, 1976):

$Display mathematics$
(1.87)

Figure 1.47 Water wave reflection with incoming waves at an oblique angle $γ=30°$ (Fig. 4.3).

(p.48) and the close-up figures illustrate that lines of constant amplitude and phase correctly intersect the elements orthogonally along their fully reflected boundaries.

Wave refraction occurs as waves propagate across changes in depth as illustrated in Fig. 1.43c for a region of uniform depth h1 on the + side of an interface to a depth of h2 on the − side. Continuity of the free surface elevation requires that the potential φ‎ takes on the same values across this interface,

$Display mathematics$
(1.88)

A second condition is developed from the component of the vector for water particle movement (1.80) normal to the interface with (1.81)

$Display mathematics$
(1.89)

which is vertically integrating over the depth of the water column to give

$Display mathematics$
(1.90)

Conservation of mass is satisfied when

$Display mathematics$
(1.91a)

with coefficients

$Display mathematics$
(1.91b)

Figure 1.48 Water wave diffraction around regularly spaced circle elements with fully reflective boundaries (Fig. 4.18). (Reprinted from Steward, 2018, Wave resonance and dissipation in collections of partially reflecting vertical cylinders, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 144(4), Figs. 5,10, with permission from ASCE.)

(p.49) that take on simpler forms for shallow water waves where $sinhkh→kh$ and $coshkh→1$. Examples of wave refraction are illustrated in Fig. 1.49 with the same incident wave angle $γ1=30°$ as shown for reflection in Fig. 1.47.

Example 1.6 Surface water waves with period and amplitude in water of depth are shown in Fig. 1.46. The frequency (1.75) is $ω=1.5708/sec$, and the wave number k may be determined iteratively by applying Newton’s method, (2.56) to the dispersion relation (1.82),

$Display mathematics$
(1.92)

Figure 1.49 Water wave refraction by a change in depth for shallow water waves (Fig. 4.4).

(p.50) where k|l is the current estimate for k and k|l+1 is the next estimate for k, and k|0 is chosen to be the wave number for deep water given in Table 1.2. The wavelength (1.75) is , and the wave velocity and group velocity (1.83) are and . The potential and velocity of water particles are

$Display mathematics$
(1.93)

giving maximum velocities at the water surface and .

Problem 1.8 Compute the frequency (ω‎ in 1/sec), wave number (k in 1/m), wavelength (L in m), phase velocity (C in m/sec), group velocity (Cg in m/sec), and maximum horizontal and vertical components of the velocity of water particles ($maxvx$ and $maxvz$ in m/sec) for water waves with the following period, amplitude, and water depth:

1. A. , , and

2. B. , , and

3. C. , , and

4. D. , , and

5. E. , , and

6. F. , , and

7. G. , , and

8. H. , , and

9. I. , , and

10. J. , , and

11. K. , , and

12. L. , , and

## (p.51) 1.3.2 Acoustics: Intensity and Vibration

The field of acoustics studies the mechanical propagation of waves through gases, liquids, and solids. Acoustics may be formulated using the wave equation for pressure waves (Airy, 1871) and separation of variables leads to solutions in terms of the Helmholtz equation (Moon and Spencer, 1961a, Chap.16)

$Display mathematics$
(1.94)

Waves may originate at a point-source as shown in Fig. 1.50, or they may be generated a large distance from the study region and be treated as plane waves. The reflection, diffraction, and refraction processes formulated for water waves are extended in this presentation of acoustics to illustrate solutions to wave propagation through and around inhomogeneities with different wave properties than those of the surrounding media.

The reflection of acoustic waves by a fully reflective boundary at $x=0$ is illustrated in Fig. 1.51. This extends the solution for plane waves and point-sources in free space from Fig. 1.42 and Eq. (1.74) using the following expressions for the complex potential

$Display mathematics$
(1.95)

Figure 1.50 Acoustics properties and processes: waves translated by the periodic motion of particles.

Figure 1.51 Wave reflection by a fully reflective boundary on the right side.

(p.52) which are developed later. Each solution satisfies a boundary condition for fully reflective surfaces with $R=1$ in

$Display mathematics$
(1.96)

where the normal component of the derivative is zero at the boundary where $x=0$ in this figure.

The diffraction produced by collection of partially reflecting elements is illustrated in Fig. 1.52. Each element reproduces the partially reflective boundary condition (1.96) with the specified reflection coefficient R, and the scattered waves generated by each element travels in the outward direction, (1.77). The wavelength for an acoustic field may be computed from (1.75)

$Display mathematics$
(1.97)

Figure 1.52 Randomly distributed circle elements, from Fig. 4.18. (Reprinted from Steward, 2018, Wave resonance and dissipation in collections of partially reflecting vertical cylinders, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 144(4), Fig. 9, with permission from ASCE.)

(p.53) for waves with specified frequency, ω‎, where the speed of wave propagation, C, is specified for the media through which waves propagate (Morse, 1936; Morse and Ingard, 1968). These results illustrate the resonance occurring within a collection of elements, where regions exist with amplitudes higher than the incoming waves, and dissipation occurs in other areas with lower amplitude. They also show the impact of reflection coefficient, where resonance is larger for more reflective objects.

The transmission of waves through elements is illustrated in Figs. 1.53 and 1.54. These elements require that waves have the same amplitude and phase across their interface

$Display mathematics$
(1.98)

similar to (1.88) for water wave transmission across changes in water depth. The additional constraints on the normal component of the partial derivative across an interface is similar to (1.91). The coefficients β‎+ and β‎ may be chosen to satisfy continuity of the gradient (p.54) of the potential in the wave equation (1.72) following Moon and Spencer (1961a). The transmission of waves through inhomogeneities is illustrated in Fig. 1.53 for elements with a shorter wavelength than that on the outside, and for larger wavelengths inside the element in Fig. 1.54. Individually, the interactions of an element with a field of plane waves generate locations within and outside the inhomogeneity with higher and lower amplification of the waves. Collectively, groups of inhomogeneities have the capacity to create zones of resonance and dissipation of wave energy, which are influenced by the acoustic properties of the medium.

Figure 1.53 Waves through elements with wavelengths approximately half the exterior wavelength. (Reprinted from Steward, 2020, Waves in collections of circular shoals and bathymetric depressions, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 146, Figs. 2,3,4, preproduction version DOI 10.1061/(ASCE)WW.1943-5460.0000570, with permission from ASCE.)

Problem 1.9 Compute the intensity, |φ‎|, for the reflected plane waves of Fig. 1.51a and Eq. (1.95) at the specified locations: , , and . These waves are specified for unit $φ0=1$ with the following wavelength and orientation.

1. A. , $θ0=0°$

2. B. , $θ0=30°$

3. C. , $θ0=60°$

4. D. , $θ0=90°$

5. E. , $θ0=0°$

6. F. , $θ0=30°$

7. G. , $θ0=60°$

8. H. , $θ0=90°$

Figure 1.54 Waves through elements with wavelengths twice the exterior wavelength. (Reprinted from Steward, 2020, Waves in collections of circular shoals and bathymetric depressions, Journal of Waterways, Ports, Coastal, and Ocean Engineering, Vol. 146, Figs. 2,3,4, preproduction version DOI 10.1061/(ASCE)WW.1943-5460.0000570, with permission from ASCE.)

# (p.55) 1.4 Studies of Deformation by Forces

The deformation of material by forces is formulated in terms of the stress (force per area) and displacement, as illustrated by the three-dimensional control block in Fig. 1.55. Surface forces at the boundary with stress components σ‎x, σ‎y, and σ‎z act normal to planes of constant x, y, and z, and a sign convention is adopted where stress is positive for tension and negative for compression. Shear stresses act tangentially to these planes (e.g., τ‎xy is the y component of shear stress acting on the x-face). Body forces with components Fx, Fy, and Fz in the directions of the coordinate axes may also act internally to deform the material. The displacements ux, uy, and uz generate strain (displacement per length) with normal strain ε‎x, ε‎y, and ε‎z, and shear strain with components that act in the same direction as shear stress (e.g., γ‎xy is the shear strain generated by shear stress τ‎xy).

A set of equations that enables solutions are given by Newton’s first law, where the sum of forces in each of the coordinate directions is zero, together with stress–strain relations for elastic deformations from a generalized Hooke’s law (Hooke, 1678),

Figure 1.55 Stresses generated by surface and body forces result in displacement and strain.

(p.56)

$Display mathematics$
(1.99a)

where E is the modulus of elasticity, ν‎ is the Poisson ratio, and strain.

$Display mathematics$
(1.99b)

is the modulus of rigidity (or shear modulus). An additional set of moment equilibrium conditions about the center of the control block in Fig. 1.55 requires shear stress and strain components in coordinate planes to be equal:

$Display mathematics$
(1.100)

The remaining components of strain are related to the displacement by

$Display mathematics$
(1.101)

Together, Eqs (1.99) and (1.101) provide 15 conditions that must be satisfied along with appropriate boundary conditions prescribed in terms of stress and/or displacement to solve for the 6 stress, 6 strain, and 3 displacement components.

(p.57) Two-dimensional problems in the xy plane may be studied using two different assumptions that reduce the equilibrium conditions (1.99) to

$Display mathematics$
(1.102)

One assumption is that of plane strain, where displacements occur in a plane and there is no strain in the z-direction. This leads to the following stress–strain relations in Hooke’s law (1.99) that, together with the compatibility conditions (1.101), directly relate displacements to stress:

$Display mathematics$
(1.103a)

Another assumption is that of plane stress, where the stresses act in a plane and there is no stress in the z-direction,

$Display mathematics$
(1.103b)

Together, Eq. (1.102) and either form of (1.103a) or (1.103b) provide five equations to solve for the five unknown components of stress (σ‎x, σ‎y, τ‎xy), and displacement (ux, uy).

Example 1.7 Body forces are important for many problems, and their solutions provide a background into which other solutions may be superimposed to study stress and displacement around elements. For example, gravity force directed in the negative y-direction in a semi-infinite dry soil below ys with bulk density ρ‎s (Jaeger, 1969, p.120) is illustrated in Fig. 1.56 and contains the following body force, stress, and displacements

$Display mathematics$
(1.104)

where these stress components satisfy the equilibrium equations for two-dimensional plane strain in (1.102) and the relations between displacement and stress in (1.103a) are satisfied. (p.58) This solution may be adapted to elasticity problems with hydrostatic pressure using $ν=12$ (Jaeger and Cook, 1976, p.112), where the density of water ρ‎w occurs below yw, giving

Figure 1.56 Body force due to gravity for dry and wet soil conditions, showing the principal stress components.

$Display mathematics$
(1.105)

Individually, these solutions for body force may be applied to the elasticity problems in saturated soils where either dry soil lies above a saturated zone in $ysat or water lies above the saturated zone in $ysat. This gives a plane strain solution in this saturated domain where water fills the voids between soil particles with porosity n beneath ysat, with Poisson ratio of ν‎s for the body forces due to the soil and $νw=12$ for the hydrostatic forces, giving

$Display mathematics$
(1.106)

where σ‎xsat, σ‎ysat, and uysat are the values occurring at $y=ysat$ due to overlying dry soil or water.

Problem 1.10 Compute the principal stresses and maximum shear stress at the specified locations for a soil below $ys=0$ with porosity $n=0.30$, density $2650(1−n)kg/m3$, Poisson ration $ν=0.3$, coefficient of elasticity $E=80×106kg/(ms2)$ and water below $yw=y0=−4m$ with density 1000kg/m3.

1. A.

2. B.

3. C.

4. D.

5. E.

6. F.

7. G.

8. H.

9. I.

10. J.

(p.59) Two-dimensional problems may be further simplified to a single partial differential equation using the Airy (1863) stress function, $F$, defined as

$Display mathematics$
(1.107)

Displacement is related to the Airy stress function by substituting these terms into (1.103a) for plane strain,

$Display mathematics$
(1.108a)

or into (1.103b) for plane stress,

$Display mathematics$
(1.108b)

Integration of $F$ in these equations provides the displacement. A governing equation for $F$ is obtained by substituting either of these sets of displacement expressions into the the equilibrium conditions (1.102),

$Display mathematics$
(1.109)

to give the biharmonic equation, (2.2e):

$Display mathematics$
(1.110)

This formulation of the Airy stress function is used in the study of elasticity (Muskhelishvili, 1953b), thermal stress (Westergaard, 1952), and viscous fluid flow (Jaeger, 1969).

The biharmonic equation may be formulated using complex functions following Muskhelishvili (1953b). Mathematics are expressed in terms of the complex variable $z$ and its conjugate $z¯$, and readers unfamiliar with complex functions are referred to Section 2.3.2, and specifically Section 3.6 where complex functions of both $z$ and $z¯$ are presented. The biharmonic equation is expressed in a complex form using (3.111),

$Display mathematics$
(1.111)

where this arrangement of complex functions $ϕ(z)$ and $χ(z)$ with their complex conjugates $ϕ¯$ and $χ¯$ provides a real function $F$ and satisfies the biharmonic equation. The Kolosov–Muskhelishvili formulas are obtained through rearrangement of these functions with $ψ=χ′$ to provide stress and displacement related to these functions by

(p.60)

$Display mathematics$
(1.112)

The constant κ‎ is related to the Poisson ratio as follows (Muskhelishvili, 1953b, p.112):

$Display mathematics$
(1.113)

Note that derivation of the Kolosov–Muskhelishvili formulas follows from insertion of the complex form of the Airy stress function (1.111) into the relations (1.107) and (1.108) and rearranging terms, and a detailed formulation is presented by Muskhelishvili (1953b). Solutions for stress and displacement are developed either in terms of $F$, or of ϕ‎ and ψ‎.

Visualization of these solutions requires the use of eigenvalue σ‎ and eigenvector v for the the stress tensor, A, which satisfies the fundamental equation

$Display mathematics$
(1.114)

where I is the identity matrix, defined to have ones on the diagonal and zeros elsewhere. Thus, matrix multiplication of A times the eigenvector v produces the same result as scaling v by the factor σ‎. The determinant of $A−σI$ must be zero since the columns in this matrix are linearly dependent, and this gives the characteristic equation

$Display mathematics$
(1.115)

where the determinant for the two-dimensional stress matrix is given by Cramer’s rule. The solution to this second-order polynomial gives the two principal stresses

$Display mathematics$
(1.116)

The principal axes along which the principal stresses operate lie in the direction of the eigenvectors of the stress matrix. These eigenvalues may be substituted into (1.114) to give a set of equations for the two eigenvectors, $v1$ and $v2$:

$Display mathematics$
(1.117)

(p.61) Each set of equations is linearly dependent, and thus provide the directions of the eigenvectors but not their length. It is common to scale eigenvectors to have length of 1, and equations for the unit eigenvectors are given by

$Display mathematics$
(1.118)

Example 1.8 Determine eigenvalues and eigenvectors for the two-dimensional matrix with coefficients

$Display mathematics$

Eigenvalues are given by (1.116)

$Display mathematics$

Unit eigenvectors are given by (1.118)

$Display mathematics$

and

$Display mathematics$

Problem 1.11 Compute the eigenvalues and unit eigenvectors for the following matrices.

1. A. $[1223]$

2. B. $[2334]$

3. C. $[3445]$

4. D. $[4556]$

5. E. $[5667]$

6. F. $[6778]$

7. G. $[7889]$

8. H. $[89910]$

9. I. $[9101911]$

10. J. $[10111112]$

## 1.4.1 Elasticity: Stress and Displacement

Solutions to boundary value problems in elasticity may be represented using the principal stresses and the maximum shear stress. The principal stresses follow directly from the eigenvectors of the Cauchy stress tensor, in (1.116):

(p.62)

$Display mathematics$
(1.119)

These are illustrated near a point where a force is applied in Figs. 1.57a and 1.57b, where regions in tension are indicated by light eigenvectors while compression uses dark eigenvectors. This figure also illustrates a set of isostress contour lines with uniform mean stress, $(σ1+σ2)/2$. The maximum shear stress is given by

$Display mathematics$
(1.120)

This is illustrated using isoshear contour lines with uniform maximum shear stress τ‎max, which are shown along with the displacement vector with components ux and uy.

Figure 1.57 Stress and displacement near a point force and a point moment (Fig. 3.35).

Figure 1.58 Background field of uniform stress (Fig. 3.34).

(p.63) A background stress field into which analytic elements are placed is illustrated in Fig. 1.58 for uniaxial tension, biaxial compression, and shear stress. Each solution is developed later at the locations indicated in the figure using the complex Kolosov functions (1.112), with the visualized stress and displacement components obtained directly from ϕ‎ and ψ‎ using (1.112). These background solutions may be superimposed with the Kolosov formulas for an analytic element to study the influence of an element in the distribution of stress and displacement; for example, the circular element in Fig. 1.60 is inserted into the settings shown in Fig. 1.58. Each element will satisfy a boundary condition of specified stress, (1.121a), where the normal and shear stress components are both zero on the boundary and the only non-zero principle stress component will act in the tangent direction.

Solutions to elastic deformation require specification of boundary conditions, as illustrated by the cantilever beam in Fig. 1.59 (with stress and displacement computed later in Example 1.9). A boundary condition may specify the normal and shear stress components

Figure 1.59 Elasticity: Cantilever beam with a force applied to the free end.

(p.64)

$Display mathematics$
(1.121a)

It may specify displacement:

$Display mathematics$
(1.121b)

where $s$ and $n$ are the tangential and normal directions to the boundary. Or, it may specify elastic support where a linear relation exists between the the normal and shear stresses and the displacement in these directions

$Display mathematics$
(1.121c)

Figure 1.60 Circle with zero traction in uniform stress (Fig. 3.37).

Figure 1.61 Circles with zero traction in uniform stress (Fig. 3.38).

(p.65) and k is the stiffness. For example, the left end of the beam in Fig. 1.59 has specified displacement and stresses are specified along the other boundaries.

Boundary conditions along circle elements are formulated by adapting the Kolosov–Muskhelishvili functions (1.112) to solve problems with boundary conditions specified along boundaries at an angle θ‎. This is accomplished using conservation of momentum in the rotated two-dimensional coordinates which gives the Mohr’s representation of stress in the r and θ‎ directions: (Muskhelishvili, 1953b, Eq. 8.7)

$Display mathematics$
(1.122)

The components of the displacement vector in the rotated coordinates are related to those in the xy plane by

$Display mathematics$
(1.123)

A form of the Kolosov–Muskhelishvili formulas in rotated coordinates is given by substituting (1.112) into these equations to give (Jaeger, 1969, p.197) (Muskhelishvili, 1953b, p.138)

$Display mathematics$
(1.124)

(p.66) Subtracting the first two equations gives an expression

$Display mathematics$
(1.125)

that enables the solution to boundary value problems with specified radial and tangential components along a circular boundary.

Boundary conditions for analytic elements may be specified in terms of the value of the complex functions ϕ‎ and ψ‎ along the boundary. A boundary may specify the normal and shear stress components σ‎0 and τ‎0 using (1.124)

$Display mathematics$
(1.126)

where θ‎ is oriented so the radial and tangential components are aligned in the normal and tangential directions. The boundary may specify displacement with components ux0 and uy0 by evaluating the functions in (1.112) at locations along the boundary:

$Display mathematics$
(1.127)

This solution is shown in Fig. 1.62 for one circle, and in Fig. 1.63 for many circles. A Robin boundary condition may specify elastic support as a linear relation between the the normal and shear stresses and the displacement in these directions

$Display mathematics$
(1.128)

and k is the stiffness. The analytic elements in this section illustrate the application of these boundary conditions, and their functions ϕ‎ and ψ‎ are formulated later in Section 3.7.

Example 1.9 The Airy stress function is applied to the plane stress problem of a cantilever beam in Fig. 1.59 with a fixed left side and a net force per unit width F applied on the right. (p.67) A solution from separation of variables (Chapter 4) follows with stress components obtained by evaluating the derivatives in (1.107)

Figure 1.62 Circle with zero displacement in uniform stress (Fig. 3.39).

Figure 1.63 Circles with zero displacement in uniform stress (Fig. 3.40).

$Display mathematics$

with $I=2/3Ly3$ being the y moment of inertia per width in the z-direction. This satisfies the traction boundary condition along the top and bottom, where the normal σ‎y and shear τ‎xy components are zero at $y=±Ly$, and along the right-hand side, where the normal stress is $σx=0$ and the integral of shear stress is equal to the prescribed force: $∫−LyLyτxydy=−F$. The displacement is obtained by substituting these partial derivatives of $F$ into (1.108) and integrating

$Display mathematics$

where use has been made of $2μ(1+ν)=E$ from (1.99). The functions f1 and f2 may be obtained fro this problem by substituting ux, uy, and $F$ into the third equation in (1.108)

$Display mathematics$

then separating terms contain x and y, and integrating the separated terms

$Display mathematics$

(p.68) The coefficients in f1 and f2 require three additional conditions that may be chosen by approximating the fixed conditions at $x=y=0$ with displacement $ux=uy=0$ (so $b1=b2=0$). The additional constant may be chosen as $a2=0$ so $∂uy/∂x=0$ and there is no rotation of the body about the origin. Substituting these coefficients and f1 and f2 into ux and uy gives the displacement vector

$Display mathematics$

The solution in Fig. 1.59 is illustrated for a lifting hook with dimensions and , coefficients for steel and $ν=0.3$, and force per width in the direction normal to the plane . The top half of the bar is in tension, the bottom half is in compression, and the maximum shear stress occurs along the left side.

Problem 1.12 For the problem in Example 1.9 and Fig. 1.59, compute the stresses σ‎x, σ‎y, τ‎xy, σ‎1, σ‎2, and τ‎max in units of N/m2, and the displacement ux and uy in units of m at the following points:

1. A. $x=0$, $y=−Ly$

2. B. $x=0$, $y=−Ly/2$

3. C. $x=0$, $y=0$

4. D. $x=0$, $y=Ly/2$

5. E. $x=0$, $y=Ly$

6. F. $x=Lx/2$, $y=−Ly$

7. G. $x=Lx/2$, $y=−Ly/2$

8. H. $x=Lx/2$, $y=0$

9. I. $x=Lx/2$, $y=Ly/2$

10. J. $x=Lx/2$, $y=Ly$

11. K. $x=Lx$, $y=−Ly$

12. L. $x=Lx$, $y=−Ly/2$

13. M. $x=Lx$, $y=0$

14. N. $x=Lx$, $y=Ly/2$

15. O. $x=Lx$, $y=Ly$

Section 1.2.1, Groundwater Flow

• Anderson (2005)

• Aravin and Numerov (1965)

• Bear (1972)

• Bulatewicz, Yang, Peterson, Staggenborg, Welch, and Steward (2010)

• Bulatewicz, Allen, Peterson, Staggenborg, Welch, and Steward (2013)

• Dagan (1981)

• Darcy (1856)

• Dupuit (1863)

• Freeze and Cherry (1979)

• Haitjema (1995)

• Janković, Steward, Barnes, and Dagan (2009)

• Muskat (1937)

• (p.69) Polubarinova-Kochina (1962)

• Steward (1998)

• Steward and Jin (2001)

• Steward (2002)

• Steward (2007)

• Steward and Ahring (2009)

• Steward and Allen (2013)

• Steward (2015)

• Strack (1981)

• Strack (1984)

• Strack (1989)

• Theis (1935)

• Thiem (1906)

• Verruijt (1982)

Section 1.2.2, Vadose Zone Flow

• Bakker and Nieber (2004a)

• Basha (1999)

• Bernoulli (1738)

• Gardner (1958)

• Kirchhoff (1894)

• Knight, Philip, and Waechter (1989)

• Mualem (1976)

• Philip (1968)

• Pullan (1990)

• Raats (1970)

• Raats (1971)

• Raats and Gardner (1974)

• Richards (1931)

• Steward (2016)

• van Genuchten (1980)

• Warrick (1974)

• Wooding (1968)

Section 1.2.3, Incompressible Fluid Flow

• Abbott and von Doenhoff (1949)

• Anderson, Jr. (1978)

• Batchelor (1967)

• Bernoulli (1738)

• Chou and Pagano (1967)

• Corke (2003)

• Hess and Smith (1967)

• Hess (1990)

• Jaeger (1969)

• Joukowsky (1912)

• Kraus (1978)

• Lagrange (1781)

• Lamb (1879)

• Lamb (1916)

• Maxwell (1870)

• Steward (2002)

• Steward, Le Grand, Janković, and Strack (2008)

• Stokes (1842)

• Stokes (1880)

• von Mises (1945)

Section 1.2.4, Thermal Conduction

• Abramowitz and Stegun (1972)

• Bakker (2013)

• Carslaw and Jaeger (1959)

• Davies (1978)

• Fourier (1878)

• Furman and Neuman (2003)

• Goodier and Hodge, Jr. (1958)

• Kuhlman and Neuman (2009)

• MacRobert (1967)

• Maxwell (1904a)

• Moon and Spencer (1961a)

• Muskhelishvili (1953b)

• Planck (1903)

• Steward (1991)

• Steward (1999)

• Steward, Le Grand, Janković, and Strack (2008)

• Steward (2015)

• Theis (1935)

• Thiem (1906)

• Westergaard (1952)

• (p.70) Zaadnoordijk (1998)

Section 1.2.5, Electrostatics

Section 1.3.1, Water Waves

• Airy (1845)

• Bascom (1980)

• Berkhoff (1972)

• Berkhoff (1976)

• Bernoulli (1738)

• Boussinesq (1871)

• Dean and Dalrymple (1984)

• Ehrenmark (1998)

• Homma (1950)

• Jonsson, Skovgaard, and Brink-Kjaer (1976)

• Lamb (1879)

• Liu, Lin and Shankar (2004)

• Liu and Li (2007)

• Longuet-Higgins (1967)

• MacCamy and Fuchs (1954)

• Mei (1989)

• Sarpkaya and Isaacson (1981)

• Smith and Sprinks (1975)

• Sommerfeld (1949)

• Steward and Panchang (2000)

• Steward (2018)

• Steward (2020)

• Stokes (1847)

• Xu, Panchang, and Demirbilek (1996)

Section 1.3.2, Acoustics

• Airy (1871)

• Brekhovskikh and Godin (1990)

• Byerly (1893)

• Courant and Hilbert (1962)

• Goodman (1996)

• Kuttruff (2007)

• Lamb (1879)

• Morse (1936)

• Morse and Ingard (1968)

• Moon and Spencer (1961a)

• Morse and Ingard (1968)

• Sommerfeld (1949)

• Sommerfeld (1972)

• Stephens and Bate (1966)

• Stewart and Lindsay (1930)

• Steward (2018)

• Steward (2020)

Section 1.4.1, Elasticity

• Airy (1863)

• Boresi (1965)

• Goodier and Hodge, Jr. (1958)

• Goursat (1898)

• Green and Zerna (1968)

• Hooke (1678)

• Jaeger (1969)

• Jaeger and Cook (1976)

• Love (1927)

• Mogilevskaya and Crouch (2001)

• Muskhelishvili (1953b)

• Sanford (2003)