Jump to ContentJump to Main Navigation
Analytic Element MethodComplex Interactions of Boundaries and Interfaces$

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

Analytic Element Method across Fields of Study

Chapter:
(p.1) 1 Analytic Element Method across Fields of Study
Source:
(p.iii) Analytic Element Method
Author(s):

David R. Steward

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.

Keywords:   groundwater, vadose zone, hydrology, thermal, conduction, electrostatics, water waves, acoustics, elasticity, geomechanics

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.

Analytic Element Method across Fields of Study

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.

Analytic Element Method across Fields of Study

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.

Analytic Element Method across Fields of Study

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

Φ=Φxx^+Φyy^+Φzz^
(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:

v=limΔx0vx(x+Δx2,y,z)ΔyΔzvx(xΔx2,y,z)ΔyΔzΔxΔyΔz+limΔy0vy(x,y+Δy2,z)ΔxΔzvy(x,yΔy2,z)ΔxΔzΔxΔyΔz+limΔz0vz(x,y,z+Δz2)ΔxΔyvz(x,y,zΔz2)ΔxΔyΔxΔyΔz
(1.2)

Evaluating the limit of this central difference equation gives

v=vxx+vyy+vzz
(1.3)

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

×v=(vzyvyz)x^+(vxzvzx)y^+(vyxvxy)z^
(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, vyxvxy, is zero when the angles ϑ1=ϑ2 in this figure are equal in the limit as Δ‎x and Δy0, 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:

v=Φ×Ψ
(1.5)

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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:

v=Φ×v=0
(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:

v=×Ψv=0
(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., vxz=vyz=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

tanϑ=vyvx=dydx=dydsdxdsvydxdsvxdyds=0
(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

dΨds=Ψxdxds+Ψydyds=0
(1.9)

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

vx=Ψy,vy=+Ψx
(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),

Ψ=Ψz^=(0,0,Ψ)
(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,

V=ACvxdy+BCvydx
(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

V=ACΨydy+BCΨxdx=ACdΨ+BCdΨ=(ΨCΨA)+(ΨCΨB)V=ΨAΨB
(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

vx=Φx=Ψyvy=Φy=+Ψx2Φx2+2Φy2=02Ψx2+2Ψy2=0
(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. Φ=sinh 2πnxmaxxymaxyminsinh 2πnxmaxxminymaxymincos 2πnyyminymaxymin, (4.26)

  8. H. Φ=sinh 2πnxmaxxymaxyminsinh 2πnxmaxxminymaxyminsin 2πnyyminymaxymin, (4.26)

  9. I. Φ=cos 2πnxxminxmaxxminsinh 2πnymaxyxmaxxminsinh 2πnymaxyminxmaxxmin, (4.26)

  10. J. Φ=sin 2πnxxminxmaxxminsinh 2πnymaxyxmaxxminsinh 2πnymaxyminxmaxxmin, (4.26)

  11. K. Φ=sinh 2πnxxminymaxyminsinh 2πnxmaxxminymaxymincos 2πnyyminymaxymin, (4.26)

  12. L. Φ=sinh 2πnxxminymaxyminsinh 2πnxmaxxminymaxyminsin 2πnyyminymaxymin, (4.26)

  13. M. Φ=cos 2πnxxminxmaxxminsinh 2πnyyminxmaxxminsinh 2πnymaxyminxmaxxmin, (4.26)

  14. N. Φ=sin 2πnxxminxmaxxminsinh 2πnyyminxmaxxminsinh 2πnymaxyminxmaxxmin, (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

saturated depth={hBunconfined (0 <hBH)Hconfined (H<hB),
(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).

Analytic Element Method across Fields of Study

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:

qx=Khx , qy=Khy , qz=KhzK=κρgμ
(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:

Qx={(hB)qxHqx,qx={Qx(hB)(0<hBH)QxH(H<hB)Qy={(hB)qyHqy,qy={Qy(hB)(0<hBH)QyH(H<hB)
(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):

Analytic Element Method across Fields of Study

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.)

Q=Φ,Φ={K(hB)22(0<hBH)KH(hB)KH22(H<hB)h=B+{2ΦK(0<ΦKH22)ΦKH+H2(KH22<Φ)
(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:

Q=Qxx+Qyy=02Φ=2Φx2+2Φy2=0
(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:

Analytic Element Method across Fields of Study

Figure 1.8 Groundwater interactions with wells and lakes.

(p.12)

h(xm,ym)=h*Φ(xm,ym)=Φ*
(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):

Q=Qxx+Qyy=R2Φ=2Φx2+2Φy2=R
(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

Φ+(xm,ym)=Φ(xm,ym),Qn+(xm,ym)=Qn(xm,ym)
(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,

Analytic Element Method across Fields of Study

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)

h(xm,ym)=hm Φ(xm,ym)=Φm(m=1,,M)
(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*,

Qn+Qn=ΔQn=w*qz=w*K*(hh*)H*
(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*,

(ΔQn)m=w*c*(hmh*),c*=K*H*
(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):

h+(xm,ym)=h(xm,ym),Qn+(xm,ym)=Qn(xm,ym)
(1.26a)

Analytic Element Method across Fields of Study

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.)

Analytic Element Method across Fields of Study

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,

Φ+(xm,ym)K+=Φ(xm,ym)K=0B+=B , H+=H
(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

Φ+=Φ,Qs=KK*w*(ΨΨ+)=w*K*KΦs
(1.27a)

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

K*Φs=0
(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:

Analytic Element Method across Fields of Study

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.)

Qn=K*Kw*(Φ+Φ),Ψ+=Ψ
(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

K*=0Qn=Ψn=0
(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:

q=Φ,Φ=Khq=qxx+qyy+qzz=02Φ=2Φx2+2Φy2+2Φz2=0
(1.29)

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

vx=qxn ,vy=qyn ,vz=qzn
(1.30)

Analytic Element Method across Fields of Study

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.)

Analytic Element Method across Fields of Study

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)

qz(x,y,z)=Bz(qxx+qyy)dz˜=(zB)(qxx+qyy)
(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)

qxx+qyy=1hB(Qxx+Qyy)1(hB)2[Qx(hB)xQy(hB)y]
(1.31b)

Analytic Element Method across Fields of Study

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):

Qxx+Qyy=R,Qx(hB)x=QxK(hB)Φx=Qx2K(hB)Qy(hB)y=QyK(hB)Φy=Qy2K(hB)
(1.31c)

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

qz=zBhB(R+Qx2+Qy2K(hB)2)
(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:

hB=RKL1(L2+L4)+L2+L4L1+L2+L4hA2+L1L1+L2+L4hD2hC=RKL4(L1+L2)+L4L1+L2+L4hA2+L1+L2L1+L2+L4hD2

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

d2Φdx2=RΦ=R2x2+c1x+c0Qx=Rxc1

(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 (0x<L1 and B1zh)

Analytic Element Method across Fields of Study

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.)

Φ=R2x(xL1)ΦA1ΦB1L1x+ΦA1Qx=R(xL12)+ΦA1ΦB1L1,ΦA1=K2hA2ΦB1=K2hB2

for zone 2 (L1x<L1+L2 and B2zh)

Φ=R2(xL1)(xL1L2)ΦB2ΦC2L2(xL1)+ΦB2Qx=R(xL1L22)+ΦB2ΦC2L2,ΦB2=K2(hBB2)2ΦC2=K2(hCB2)2

for zone 3 (L1x<L1+L2 and B3zH3)

Φ=ΦB3ΦC3L2(xL1)+ΦB3Qx=ΦB3ΦC3L2,ΦB3=KH3hBKH322ΦC3=KH3hCKH322

and for zone 4 (L1+L2xL1+L2+L4 and B4zh)

Φ=R2(xL1L2)(xL1L2L4)ΦC4ΦD4L4(xL1L2)+ΦC4Qx=R(xL1L2L42)+ΦC4ΦD4L4,ΦC4=K2hC2ΦD4=K2hD2

similar to Strack (1981).

Analytic Element Method across Fields of Study

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 L1=1000 m, L2=2000 m, L4=1000 m, B2=H3=20 m, K=30 m/day, R=0.0003 m/day, n=0.25, hA=28 m, and hD=20 m. 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:

Qx=Qx|Zone 1=1.14 m2/dayQx+=Qx|Zone 2+Qx|Zone 3=1.14 m2/day(x=L1)

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

Qx=Qx|Zone 2+Qx|Zone 3=1.74 m2/dayQx+=Qx|Zone 4=1.74 m2/day(x=L1+L2)

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. x=0 m and z=h (zone 1)

  2. B. x=0 m and z=h/2 (zone 1)

  3. C. x=500 m and z=h (zone 1)

  4. D. x=1000 m and z=h (zone 1)

  5. E. x=1000+ m and z=h (zone 2)

  6. F. x=2000 m and z=h (zone 2)

  7. G. x=3000 m and z=h (zone 2)

  8. H. x=1000+ m and z=H3 (zone 3)

  9. I. x=2000 m and z=H3 (zone 3)

  10. J. x=3000 m and z=H3 (zone 3)

  11. K. x=3000+ m and z=h (zone 4)

  12. L. x=3500 m and z=h (zone 4)

  13. M. x=4000 m and z=h (zone 4)

  14. N. x=4000 m 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)

h=p+z ,p=Pρg
(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),

qx=K(p)px,qy=K(p)py,qz=K(p)pzK(p)
(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

q=x[K(p)px]+y[K(p)py]+z[K(p)pz]+K(p)z=0
(1.35)

Analytic Element Method across Fields of Study

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

F=pK(p˜) dp˜Fx=Kpx,Fy=Kpy,Fz=Kpz
(1.36)

to provide a linear form of the Richards equation

2F=2Fx2+2Fy2+2Fz2=dKdp1KFz
(1.37)

where the last term was rearranged using the chain rule for differentiation Kz=dKdppz, with pz 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),

K(p)=KseapdKdp1K=a2F=aFz
(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):

F=Φekza=2k2Φ=2Φx2+2Φy2+2Φz2=k2Φ
(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):

F=pKseap˜ dp˜=Fseap,Fs=Ksa
(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.

Analytic Element Method across Fields of Study

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,

Φ=Fsek(z+2p)p=12klnΦekzFs=z2+12klnΦFs
(1.41)

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

qx=Fx=ekzΦxqy=Fy=ekzΦyqz=FzaF=ekz(Φz+kΦ)
(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

p+=p(ek+zFs+Φ+)12k+=(ekzFsΦ)12k
(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

Analytic Element Method across Fields of Study

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

qn+=qn ek+z(Φ+n+nzk+Φ+)=ekz(Φn+nzkΦ)
(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.

Analytic Element Method across Fields of Study

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.)

Analytic Element Method across Fields of Study

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

vx=qxθ ,vy=qyθ ,vz=qzθ
(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

θ=θr+(θsθr)eap
(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),

p=12klnc1e2kz+c2Fs,qz=2kc2

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

(p.26)

p(z0)=p0qz(z0)=Rc1=(Fse2kp0R2k)e2kz0,c2=R2k

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

p12klnR2kFs=1alnRKs
(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 R=0.1 m/year. The pressure head at the interfaces are:

zA=0 mzB=2 mzC=3 mzD=5 m,Ks1=0.02 m/dayKs2=1 m/dayKs3=0.02 m/day,a1=2/ma2=5/ma3=2/m,pA=0.0 mpB=1.72 mpC=1.64 mpD=2.13 m

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. R=0.01 m/yr

  2. B. R=0.02 m/yr

  3. C. R=0.03 m/yr

  4. D. R=0.04 m/yr

  5. E. R=0.05 m/yr

  6. F. R=0.06 m/yr

  7. G. R=0.07 m/yr

  8. H. R=0.08 m/yr

  9. I. R=0.09 m/yr

  10. J. R=0.011 m/yr

  11. K. R=0.012 m/yr

  12. L. R=0.013 m/yr

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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

P1ρg+|v1|22g+z1=P2ρg+|v2|22g+z2
(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)

v=Φ
(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):

v=02Φ=2Φx2+2Φy2+2Φz2=0
(1.50)

Analytic Element Method across Fields of Study

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)

vx=Ψy,vy=Ψx,vz=0
(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

vr=1rΨz,vθ=0,vz=1rΨr
(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.

Analytic Element Method across Fields of Study

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

vx=vx0vx0x2y2(x2+y2)2Γ2πyx2+y2,vy=vx02xy(x2+y2)2+Γ2πxx2+y2
(1.53)

Analytic Element Method across Fields of Study

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 (x=0,y=1 m) and bottom (x=0,y=-1 m) of a rotating cylinder in water ρ=1 g/cm3 for the specified background velocity vx0 and circulation Γ‎:

  1. A. vx0=1 m/s, Γ=0.25π/s

  2. B. vx0=2 m/s, Γ=0.25π/s

  3. C. vx0=3 m/s, Γ=0.25π/s

  4. D. vx0=1 m/s, Γ=0.5π/s

  5. E. vx0=2 m/s, Γ=0.5π/s

  6. F. vx0=3 m/s, Γ=0.5π/s

  7. G. vx0=1 m/s, Γ=π/s

  8. H. vx0=2 m/s, Γ=π/s

  9. I. vx0=3 m/s, Γ=π/s

Analytic Element Method across Fields of Study

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):

q=KTq=Φ,Φ=KT
(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):

q=ρcTt2Φ=1DΦt,D=Kρc
(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):

2Φ=2Φx2+2Φy2+2Φz2=0
(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).

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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):

2Φ=2Φx2+2Φy2+2Φz2=R
(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.

Analytic Element Method across Fields of Study

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:

qn=0,Ψ uniform
(1.58)

Analytic Element Method across Fields of Study

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.)

Analytic Element Method across Fields of Study

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,

qn=K*Tn=K*KΦn=K*K(Φ+Φw*)
(1.59)

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

qn=1δ(Φ+Φ),δ=Kw*K*
(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

qn+(xm,ym,zm)=qn(xm,ym,zm),Ψ+(xm,ym,zm)=Ψ(xm,ym,zm)
(1.61a)

Analytic Element Method across Fields of Study

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

T+(xm,ym,zm)=T(xm,ym,zm)1K+Φ+(xm,ym,zm)=1KΦ(xm,ym,zm)
(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

Analytic Element Method across Fields of Study

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

(p.36)

Φ=Q4πE1(u)+Φ0qr=Q2πeur,u=r24Dt,E1(u)=ueu˜u˜du˜
(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

Φ=Q2πlnr+constant,qr=Φr=Q2π1r
(1.63)

Analytic Element Method across Fields of Study

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 t=106 sec, at t=1 sec, and at the final steady-state conditions. Thermal diffusivity values from (Moon and Spencer, 1961a, p.213).

  1. A. copper, D=113 m2/sec, r=1 mm

  2. B. copper, D=113 m2/sec, r=1 cm

  3. C. cast iron, D=13.3 m2/sec, r=1 mm

  4. D. cast iron, D=13.3 m2/sec, r=1 cm

  5. E. concrete, D=0.58 m2/sec, r=1 mm

  6. F. concrete, D=0.58 m2/sec, r=1 cm

  7. G. dry soil, D=0.31 m2/sec, r=1 mm

  8. H. dry soil, D=0.31 m2/sec, r=1 cm

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):

E=ΦE=ρε2Φ=ρε
(1.64)

Analytic Element Method across Fields of Study

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)

2Φ=2Φx2+2Φy2+2Φz2=0
(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:

I=VR,I=AκEdA,V12=Φ1Φ2
(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

Φ uniform along boundary
(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.

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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:

Φ+(xm,ym,zm)=Φ(xm,ym,zm)
(1.68)

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

D=εE,Dn+=ε+En+Dn+=ε+En+ε+Φ+n=εΦn
(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

Es(xm,ym,zm)=0
(1.70)

Analytic Element Method across Fields of Study

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)

Φ=Ex0x+Ex0r02xx2+y2
(1.71)

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

Analytic Element Method across Fields of Study

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 r0=1 cm.

  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):

2Φx2+2Φy2+2Φz2=1c22Φt2
(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):

Analytic Element Method across Fields of Study

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

(p.42)

2φx2+2φy2+2φz2=k2φ,Φ(x,y,z,t)=[φ(x,y,z)eiωt]
(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:

φ=φ0eik(xcosθ0+ysinθ0)plane waves from (4.9)φ=S2πH0(1)(kr)cylindrical point source from (4.42)φ=S4πh0(1)(kr)spherical point source from (4.105)
(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

L=2πk,k=2πLT=2πω,ω=2πTC=LT=ωk
(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

φn+αφ=0
(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)

Analytic Element Method across Fields of Study

Figure 1.43 Wave interactions at boundaries and interfaces.

(p.44)

limrr(φrikφ)=0
(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:

φ+=φ,β+φ+n=βφn
(1.78)

Analytic Element Method across Fields of Study

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

k1sinγ1=k2sinγ2
(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: (x1,y1)=(0,0) m, (x2,y2)=(10,10) m, and (x3,y3)=(20,20) m. This wave field is specified for φ0=10 m2/s, k=0.25/ m 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:

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

Figure 1.46 Monochromatic waves in water of constant depth.

(p.46)

v=Φ,Φ(x,y,z,t)=Z(z)[φ(x,y)eiωt]=Z(z)[(φ)cosωt+(φ)sinωt]
(1.80)

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

2φx2+2φy2=k2φ,d2Zdz2=k2ZZ=coshk(z+h)coshkh
(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)

ω2=kgtanhkh
(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)

C=ωk=gktanhkh,Cg=ωk=C2(1+2khsinh2kh)
(1.83)

The hyperbolic functions in these relations approach asymptotic limits for shallow water (where tanhkhkh and sinh2kh2kh) and deep water (where tanhkh1 and sinh2khe2kh/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<L

2h<L<20h

L<2h

Range of wave number

kh<π/10

π/10<kh<π

π<kh

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

φ=φ0[eik(ssinγncosγ)+Reik(ssinγ+ncosγ+β)],φ0=gA0ω
(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

φn=ikcosγφ0[eik(ssinγncosγ)+Reik(ssinγ+ncosγ+β)]=ikcosγ1Reikβ1+Reikβφ(n=0)
(1.85)

It is common to approximate the coastal boundary condition as

φn+αφ=0,α=ik1R1+R(n=0)
(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):

|φ|n=0,argφn=k1R1+R
(1.87)

Analytic Element Method across Fields of Study

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,

φ+=φ
(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)

vn=Φn=coshk(z+h)coshkhn(φeiωt)
(1.89)

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

Vn=h0vndz=h0coshk(z+h)coshkhdzn(φeiωt)=sinhkhkcoshkhn(φeiωt)
(1.90)

Conservation of mass is satisfied when

Vn+=Vnβ+φ+n=βφn
(1.91a)

with coefficients

β+=sinhk1h1k1coshk1h1β=sinhk2h2k2coshk2h2,β+=h1(k1h1<π10)β=h2(k2h2<π10)
(1.91b)

Analytic Element Method across Fields of Study

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 sinhkhkh and coshkh1. 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 T=4 sec and amplitude A0=0.4 m in water of depth h=10 m 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),

k|l+1=k|l+ω2k|lgtanhk|lhk|lghcosh2k|lh+gtanhk|lh k|0=ω2g=(1.5708/sec)29.80665m/sec2=0.2516/mk|1=0.2547/mk|2=0.2547/m
(1.92)

Analytic Element Method across Fields of Study

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 L=24.67 m, and the wave velocity and group velocity (1.83) are C=6.167 m/sec and Cg=3.276 m/sec. The potential and velocity of water particles are

Φ=coshk(z+h)coshkhgωA0sin(kxωt)vx=Φx=coshk(z+h)coshkhgkωA0cos(kxωt)vz=Φz=sinhk(z+h)coshkhgkωA0sin(kxωt)
(1.93)

giving maximum velocities at the water surface maxvx=0.6361 m/sec and maxvz=0.6283 m/sec.

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. T=3 sec, A0=0.04 m, and h=0.1 m

  2. B. T=3 sec, A0=0.04 m, and h=1 m

  3. C. T=3 sec, A0=0.04 m, and h=10 m

  4. D. T=5 sec, A0=0.1 m, and h=0.5 m

  5. E. T=5 sec, A0=0.1 m, and h=5 m

  6. F. T=5 sec, A0=0.1 m, and h=50 m

  7. G. T=10 sec, A0=0.5 m, and h=2 m

  8. H. T=10 sec, A0=0.5 m, and h=20 m

  9. I. T=10 sec, A0=0.5 m, and h=200 m

  10. J. T=20 sec, A0=1 m, and h=5 m

  11. K. T=20 sec, A0=1 m, and h=50 m

  12. L. T=20 sec, A0=1 m, and h=500 m

(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)

2φ=2φx2+2φy2+2φz2=k2φ
(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

φ=φ0eik(xcosθ0+ysinθ0)+φ0eik(xcosθ0+ysinθ0),plane waves from (4.11)φ=S2πH0(1)(k(xd)2+y2)+S2πH0(1)(k(x+d)2+y2),cylindrical point source from (4.43)φ=S2πh0(1)(k(xd)2+y2)+S2πh0(1)(k(x+d)2+y2),spherical point source from (4.106)
(1.95)

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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

φn+αφ=0,α=ik1R1+R(n=0)
(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)

L=2πCω
(1.97)

Analytic Element Method across Fields of Study

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

φ+=φ,β+φ+n=βφn
(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.

Analytic Element Method across Fields of Study

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: (x1,y1)=(0,0) m, (x2,y2)=(25,0) m, and (x3,y3)=(50,0) m. These waves are specified for unit φ0=1 with the following wavelength and orientation.

  1. A. L=20 m, θ0=0°

  2. B. L=20 m, θ0=30°

  3. C. L=20 m, θ0=60°

  4. D. L=20 m, θ0=90°

  5. E. L=70 m, θ0=0°

  6. F. L=70 m, θ0=30°

  7. G. L=70 m, θ0=60°

  8. H. L=70 m, θ0=90°

Analytic Element Method across Fields of Study

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),

Analytic Element Method across Fields of Study

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

(p.56)

σxx+τxyy+τxzz+Fx=0,σyy+τxyx+τyzz+Fy=0,σzz+τxzx+τyzy+Fz=0,εx=1E[σxν(σy+σz)],εy=1E[σyν(σz+σx)],εz=1E[σzν(σx+σy)],γxy=τxyμ,γyz=τyzμγxz=τxzμ
(1.99a)

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

μ=E2(1+ν)
(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:

τxy=τyxτyz=τzyτzx=τxzγxy=γyzγyz=γzyγzx=γxz
(1.100)

The remaining components of strain are related to the displacement by

εx=uxx,γxy=uxy+uyxεy=uyy,γyz=uyz+uzyεz=uzz,γxz=uzx+uxz
(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

σxx+τxyy+Fx=0σyy+τxyx+Fy=0
(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:

εz=0σz=ν(σx+σy)γyz=0γyz=0γxz=0γxz=0uxx=εx=12μ[σxν(σx+σy)]uyy=εy=12μ[σyν(σx+σy)]uxy+uyx=γxy=1μτxy
(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,

σz=0εz=ν(σx+σy)2μ(1+ν)τyz=0γyz=0τxz=0γxz=0uxx=εx=12μ[σxν1+ν(σx+σy)]uyy=εy=12μ[σyν1+ν(σx+σy)]uxy+uyx=γxy=1μτxy
(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

Fx=Fz=0Fy=ρsg{σx=σz=ν1νρsg(yys)σy=ρsg(yys)τxy=τyz=τxz=0{ux=uz=0uy=12μ12ν1νρsg(yys)22(yys)
(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

Analytic Element Method across Fields of Study

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

Fx=Fz=0Fy=ρwg{σx=σy=σz=ρwg(yyw)τxy=τyz=τxz=0ux=uy=uz=0(yyw)
(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<yys or water lies above the saturated zone in ysat<yyw. 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

Fx=Fz=0Fy=(ρs+nρw)g(yysat){σx=σz=[νs(ρsρw+nρw)1νs+ρw]g(yyw)+σxsatσy=(ρs+nρw)g(yy0)+σysatτxy=ux=uz=0uy=12μ12νs1νs(ρsρw+nρw)g(yys)22+uysat
(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(1n)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. y=1 m

  2. B. y=2 m

  3. C. y=3 m

  4. D. y=4 m

  5. E. y=5 m

  6. F. y=6 m

  7. G. y=7 m

  8. H. y=8 m

  9. I. y=9 m

  10. J. y=10 m

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

σx=2Fy2,σy=2Fx2,τxy=2Fxy
(1.107)

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

uxx=12μ[2Fy2ν2F]uyy=12μ[2Fx2ν2F],uxy+uyx=1μ2Fxy
(1.108a)

or into (1.103b) for plane stress,

uxx=12μ(2Fy2ν1+ν2F)uyy=12μ(2Fx2ν1+ν2F),uxy+uyx=1μ2Fxy
(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),

2 y2(uxx)+2 x2(uyy)=2 xy(uxy+uyx)
(1.109)

to give the biharmonic equation, (2.2e):

4F=4Fx4+24Fx2y2+4Fy4=0
(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),

4F=164Fz2z¯2=0F=z¯ϕ+zϕ¯+χ+χ¯2=(z¯ϕ+χ)
(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)

σx+σy2=2ϕσyσx2+iτxy=z¯ϕ+ψux+iuy=12μ(κϕzϕ¯ψ¯)σx=(2ϕz¯ϕψ)σy=(2ϕ+z¯ϕ+ψ)τxy=(z¯ϕ+ψ)
(1.112)

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

κ={34νplane strain3ν1+νplane stress
(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

Av=σv,A=[σxτxyτxyσy](AσI)v=0,I=[1001]
(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

det(AσI)=|σxστxyτxyσyσ|=(σxσ)(σyσ)τxy2=0=0
(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

σ1=σx+σy2+(σxσy2)2+τxy2σ2=σx+σy2(σxσy2)2+τxy2
(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:

(Aσ1I)v1=[σxσ1τxyτxyσyσ1]v1=0(Aσ2I)v2=[σxσ2τxyτxyσyσ2]v2=0
(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

v^1=±[τxyτxy2+(σxσ1)2σxσ1τxy2+(σxσ1)2]=±[σyσ1τxy2+(σyσ1)2τxyτxy2+(σyσ1)2]v^2=±[τxyτxy2+(σxσ2)2σxσ2τxy2+(σxσ2)2]=±[σyσ2τxy2+(σyσ2)2τxyτxy2+(σyσ2)2]
(1.118)

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

A=[13335]

Eigenvalues are given by (1.116)

σ1=13+52+(1352)2+3×3=9+5=14σ2=13+52(1352)2+3×3=95=4

Unit eigenvectors are given by (1.118)

v^1=±[332+(1314)2131432+(1314)2]=±[51432+(514)2332+(514)2]=±[310110]

and

v^2=±[332+(134)213432+(134)2]=±[5432+(54)2332+(54)2]=±[110310]

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)

σ1=σx+σy2+(σxσy2)2+τxy2σ2=σx+σy2(σxσy2)2+τxy2
(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

τmax=σ1σ22=(σxσy2)2+τxy2
(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.

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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

(p.64)

shear specifiedσn(s)=σ0τns(s)=τ0sD
(1.121a)

It may specify displacement:

displacement specifiedux(s)=ux0uy(s)=uy0sD
(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

elastic support specifiedσn(s)=kun0τns(s)=kus0sD
(1.121c)

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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)

σr=σx+σy2+σxσy2cos2θ+τxysin2θσθ=σx+σy2σxσy2cos2θτxysin2θτrθ=σxσy2sin2θ+τxycos2θ
(1.122)

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

ur+iuθ=(ux+iuy)eiθ
(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)

σr+σθ2=2ϕσθσr2+iτrθ=e2iθ(z¯ϕ+ψ)ur+iuθ=eiθ2μ(κϕzϕ¯ψ¯)
(1.124)

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

σriτrθ=2ϕe2iθ(z¯ϕ+ψ)
(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)

2ϕ(s)[z¯ϕ(s)+ψ(s)]e2iθ(s)=σ0iτ0
(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:

12μ[κϕ(s)zϕ(s)¯ψ(s)¯]=ux0+iuy0
(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

elastic support specifiedσn(s)=kun0τns(s)=kus0sD
(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)

Analytic Element Method across Fields of Study

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

Analytic Element Method across Fields of Study

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

F=F2I(xLx)(y33Ly2y)σx=2Fy2=FI(xLx)yσy=2Fx2=0τxy=2Fxy=F2I(y2Ly2)

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

uxx=FEI(xLx)yuyy=FEIν(xLx)yux=FEI(x22Lxx)y+f1(y)uy=FEIν(xLx)y22+f2(x)

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)

uxy+uyx=1μ2FxyFEI(x22Lxx)+f1y+FEIνy22+f2x=FEI(1+ν)(y2Ly2)

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

f1yFEI(2+ν)y22=a1f2xFEI(x22Lxx)=a2f1=FEI(2+ν)y36+a1y+b1f2=FEI(x36Lxx22)+a2x+b2a1+a2=FEI(1+ν)Ly2

(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

ux=FEI(x22Lxx)y+FEI(2+ν)y36FEI(1+ν)Ly2yuy=FEIν(xLx)y22+FEI(x36Lxx22)

The solution in Fig. 1.59 is illustrated for a lifting hook with dimensions Lx=0.5 m and Ly=0.05 m, coefficients for steel E=200×109 N/m2 and ν=0.3, and force per width in the direction normal to the plane F=1×106 N/m. 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

Further Reading

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

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)

  • Selvadurai (2000)

  • Westergaard (1952)