Analytic Element Method across Fields of Study
Analytic Element Method across Fields of Study
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.
(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.
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.
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
where , , and 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:
Evaluating the limit of this central difference equation gives
A third vector operator is the curl, ∇ ×, which may be written in Cartesian coordinates for the vector v:
The curl is used to quantify infinitesimal rotation as illustrated in Fig. 1.4c, where the component of the curl in the z-direction, , is zero when the angles in this figure are equal in the limit as Δx and , 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:
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:
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., ) and is invariant in the direction normal to this plane (i.e., ). 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
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
A vector potential for a two-dimensional vector field is given by multiplying the Lagrange stream function times , a unit vector pointing in the positive z-direction (Steward, 2002),
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,
Substituting the relationship between the vector field and stream function (1.10), integrating, and evaluating the stream function at the limits of integration gives
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
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 , and . The equation numbers are listed for where these functions occur later.
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
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).
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:
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:
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):
when K, B, and H are piecewise constant in an aquifer layer.
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:
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):
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 -direction is also continuous at locations (xm, ym) on the polygon’s boundary is
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,
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, , is equal to minus the vertical component of flow through the riverbed times the width of the river w*,
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):
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 -direction tangential to the heterogeneity may be expressed in terms of the stream function or the partial derivative of the potential by
In the limiting case of an infinitely conductive heterogeneity, the potential is uniform along the element:
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:
In the limiting case of zero conductivity, there is no normal component of flow, and the stream function is uniform along the element
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:
The average groundwater velocity v is related to the specific discharge by
(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)
Together, the equations in (1.31) give the vertical component of flow:
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:
Within each zone, the general solution to the Poisson equation, (1.21), is
for zone 2 ( and )
for zone 3 ( and )
and for zone 4 ( and )
similar to Strack (1981).
(p.20) This groundwater head and velocity is illustrated in Fig. 1.17b for the variables , , , , , , , , 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 by computing the discharge per width on the left side in zone 1 and on the right side in zones 2 and 3:
Likewise, continuity of flow may be verified at on the left and right sides of the interface
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:
A. and (zone 1)
B. and (zone 1)
C. and (zone 1)
D. and (zone 1)
E. and (zone 2)
F. and (zone 2)
G. and (zone 2)
H. and (zone 3)
I. and (zone 3)
J. and (zone 3)
K. and (zone 4)
L. and (zone 4)
M. and (zone 4)
N. and (zone 4)
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)
where head is partitioned into pressure, P, and elevation terms with the z-axis directed upward against gravity, and the kinetic velocity head 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 at , 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),
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
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).
(p.22) This equation may be rearranged using the matric flux potential, denoted F, defined by the Kirchhoff transformation (Kirchhoff, 1894) as
to provide a linear form of the Richards equation
where the last term was rearranged using the chain rule for differentiation , with 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),
where a is the sorptive number, and the saturated conductivity Ks occurs at . 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):
(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,
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 and a− embedded within a uniform soil with properties , 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
(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
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.
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
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
with coefficients, (4.20), obtained by matching specified values of pressure head and specific discharge at elevation z0:
In the limit as z increases in each layer the pressure approaches(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 and . 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:
Problem 1.3 Compute the pressure head at elevations zA, zB, zC, and zD for Example 1.2, with the following recharge rates:
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 is zero and fluid flows parallel to its surface. Wings are designed so the fluid moves at with a higher magnitude of velocity at point 1 along the top surface than at point 2 along the bottom. This creates a lower pressure P along the top that may be quantified using the Bernoulli (1738) equation
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)
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):
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
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.
(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(1.53)
(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 Γ:
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):
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):
where ρ is the density, c is the specific heat, and 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):
(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):
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.
(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:
(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,
which is approximated as the jump in thermal potential across the thin feature of width w*. This gives a Robin boundary condition,
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
(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
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(1.62)
where 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(1.63)
(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 and thermal diffusivity at distance r from a point-sink in the specified material at initial conditions , at , at , and at the final steady-state conditions. Thermal diffusivity values from (Moon and Spencer, 1961a, p.213).
A. copper, ,
B. copper, ,
C. cast iron, ,
D. cast iron, ,
E. concrete, ,
F. concrete, ,
G. dry soil, ,
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):
(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)
This formulation is grounded in Ohm’s law (Ohm, 1827), where the current I is equal to voltage V divided by resistance R:
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 , 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
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.
(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:
They also satisfy a condition where the normal component of the displacement, D, is continuous across the interface:
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
(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.
where Ex0 is the background electric field in the x-direction, and the element is of radius and centered at the origin.
(p.41) Problem 1.6 Compute the electric field E for the field in Example 1.5 at the specified locations for unit with .
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):
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):
where k is the wave number (Courant and Hilbert, 1962). This formulation utilizes complex variables with the imaginary number , and an introduction to complex functions with real and imaginary parts, , 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:
(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 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
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
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)
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 -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:
relating the change in wave direction to the wave number on each side of the interface.
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 , 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:
and v is the water particle velocity. This gives solutions to the Helmholtz equation for φ, and an ordinary differential equation of Z
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)
The hyperbolic functions in these relations approach asymptotic limits for shallow water (where and ) and deep water (where and ) 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) - 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
Range of wavelength
Range of wave number
Dispersion relation, (1.82)
Wave velocity, (1.83)
Group velocity, (1.83)
where the reflection coefficient R takes on values from for a coast that fully absorbs wave energy to for full reflection of wave energy, and β is a phase shift between incoming and reflected waves at the boundary (Berkhoff, 1976). A coastal boundary condition may be developed by evaluating the partial derivative with respect to the normal direction on the boundary
It is common to approximate the coastal boundary condition as
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 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 .
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):
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,
which is vertically integrating over the depth of the water column to give
Conservation of mass is satisfied when
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 , and the wave number k may be determined iteratively by applying Newton’s method, (2.56) to the dispersion relation (1.82),(1.92)
(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(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 ( and in m/sec) for water waves with the following period, amplitude, and water depth:
A. , , and
B. , , and
C. , , and
D. , , and
E. , , and
F. , , and
G. , , and
H. , , and
I. , , and
J. , , and
K. , , and
L. , , and
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)
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 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
where the normal component of the derivative is zero at the boundary where 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)
(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.
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.
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 with the following wavelength and orientation.
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),
where E is the modulus of elasticity, ν is the Poisson ratio, and strain.
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:
The remaining components of strain are related to the displacement by
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 x–y plane may be studied using two different assumptions that reduce the equilibrium conditions (1.99) to
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:
Another assumption is that of plane stress, where the stresses act in a plane and there is no stress in the z-direction,
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(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 (Jaeger and Cook, 1976, p.112), where the density of water ρw occurs below yw, giving(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 or water lies above the saturated zone in . 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 for the hydrostatic forces, giving(1.106)
where σxsat, σysat, and uysat are the values occurring at 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 with porosity , density , Poisson ration , coefficient of elasticity and water below with density 1000kg/m3.
(p.59) Two-dimensional problems may be further simplified to a single partial differential equation using the Airy (1863) stress function, , defined as
Displacement is related to the Airy stress function by substituting these terms into (1.103a) for plane strain,
or into (1.103b) for plane stress,
Integration of in these equations provides the displacement. A governing equation for is obtained by substituting either of these sets of displacement expressions into the the equilibrium conditions (1.102),
to give the biharmonic equation, (2.2e):
The biharmonic equation may be formulated using complex functions following Muskhelishvili (1953b). Mathematics are expressed in terms of the complex variable and its conjugate , and readers unfamiliar with complex functions are referred to Section 2.3.2, and specifically Section 3.6 where complex functions of both and are presented. The biharmonic equation is expressed in a complex form using (3.111),
where this arrangement of complex functions and with their complex conjugates and provides a real function 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
The constant κ is related to the Poisson ratio as follows (Muskhelishvili, 1953b, p.112):
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 , 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
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 must be zero since the columns in this matrix are linearly dependent, and this gives the characteristic equation
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
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, and :
(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
Example 1.8 Determine eigenvalues and eigenvectors for the two-dimensional matrix with coefficients
Eigenvalues are given by (1.116)
Unit eigenvectors are given by (1.118)
Problem 1.11 Compute the eigenvalues and unit eigenvectors for the following matrices.
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):
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, . The maximum shear stress is given by
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.
(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
It may specify displacement:
where and 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
(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)
The components of the displacement vector in the rotated coordinates are related to those in the x–y plane by
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)
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:
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
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)
with 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 , and along the right-hand side, where the normal stress is and the integral of shear stress is equal to the prescribed force: . The displacement is obtained by substituting these partial derivatives of into (1.108) and integrating
then separating terms contain x and y, and integrating the separated terms
(p.68) The coefficients in f1 and f2 require three additional conditions that may be chosen by approximating the fixed conditions at with displacement (so ). The additional constant may be chosen as so 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
The solution in Fig. 1.59 is illustrated for a lifting hook with dimensions and , coefficients for steel and , 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.
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)
• Selvadurai (2000)
• 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
• Choudhury (1989)
• Grant and Phillips (1990)
• Hammond (1986)
• MacRobert (1967)
• Maxwell (1890a)
• Maxwell (1890b)
• Maxwell (1904b)
• Muskat (1932)
• Moon and Spencer (1960)
• Moon and Spencer (1961a)
• Moon and Spencer (1961b)
• Ohm (1827)
• Slater and Frank (1947)
• Sommerfeld (1952)
• Steward (2015)
• Thomson (1904)
• Wait (1970)
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