General solutions for the displacements, strains and changes in stress in a transversely isotropic linear elastic (TILE) half-space subject to a point load were published in 1998. They represent a massive advance on the nineteenth-century solutions by Boussinesq, Cerruti and others, which were for fully isotropic linear elastic half-spaces. The 1998 solutions were for loads at a general depth below the surface and were criticised for being difficult. The present paper deduces more accessible expressions for the special case of loads that are on the surface. This allows one clarification and one potential error in the original equations to be identified and avoided. With support from existing data of geological materials, agreement is found with a previous claim that a limit proposed in 1970 on one of the shear moduli for a TILE material is invalid. Corrected equations derived herein can be the basis of calculating ground stresses and movements due to linear responses to distributed surface loads, behaviours of flexible foundations and stiffnesses for rigid surface footings.
Notation
Principal stress and principal strain are taken positive in compression. ‘Vertical’ means parallel to the axis of material symmetry (the z-direction), and horizontal means at 90° to the axis. Displacements are assumed to be small enough that second-order effects due to the associated small changes of geometry can be neglected. Symbols for stress refer to changes of effective stress but do not carry the traditional prime marker. Total stresses and constant-volume parameters are identified by superscript u. The ± and ∓ signs are such that the upper sign is for subscript i = 1 and the lower is for i = 2.
- Aij
constrained modulus, Liao and Wang (1998) notation
- a, b
miscellaneous parameters defined where used
- Cab
(with letter subscripts) dimensionless compliance factor characterising the surface displacement in direction a due to a point load in direction b
- Cij
(with integer subscripts) constrained modulus, some sources in the literature
- d
Infinitesimal
- Eh, Ev
Young’s moduli in the horizontal and vertical directions, respectively
- G
shear modulus (un-subscripted for a fully isotropic linear elastic (FILE) material, subscripted for a transversely isotropic linear elastic (TILE) material)
- k
Liao and Wang (1998) factor with dimensions of inverse modulus
- K′
dimensionless factor
- mi
two (i = 1, 2) of Liao and Wang (1998) dimensionless factors
- P
point load
- q
dimensionless ratio A 11/A 33
- R
actual distance from the origin to a point
- Ri,
modified distances
- r, θ, z
cylindrical coordinates
- S
dimensionless parameter involved in solving the characteristic equation
- T
one of Liao and Wang (1998) factors with dimensions of inverse modulus
- U
displacement at a general point in the half-space
- u
change in pore pressure
- ua
(with letter subscript) surface displacement in direction a
- ui
(i = 1, 2, 3) three of Liao and Wang (1998) dimensionless factors that multiply the z-coordinate
- W
strain energy, computed as work done starting from an initial stress of zero
- x, y, z
Cartesian coordinates, z in the direction of the axis of material symmetry
- γij
small engineering shear strain in a plane containing the directions of the i and j axes
- Δ, Δ1, Δ2
determinant and two of the factors in the determinant of the TILE compliance matrix
- ϵ
dimensionless infinitesimal
- ϵr, ϵθ, ϵz
small radial, circumferential and vertical normal strains, respectively
- λi
dimensionless factor
- μ
Poisson’s ratio (un-subscripted for an FILE material, subscripted for a TILE material)
- σ
change in normal stress
- τ
change in shear stress
Subscripts
- ax
axial
- h
relating to a lateral direction
- hh
for shearing in a horizontal (z = constant) plane
- i
integer index
- ij
integer indices indicating directions; (123) = (rθz) for a cylindrical coordinate system, (123) = (xyz) for a Cartesian coordinate system
- q
deviatoric
- r, θ, z
associated with radial, circumferential and vertical directions, respectively
- rad
radial in the triaxial cell
- v
relating to the direction of material symmetry (taken to be vertical)
- vh
for shearing in a vertical plane (one containing a line parallel to the z-axis)
- vol
volumetric
- x, y, z
associated with Cartesian directions x, y and z
Superscripts
Introduction
Geomaterials, including soils and rocks, are known to be anisotropic in their small-strain elastic responses. Casagrande and Carillo (1944) made the distinction between induced anisotropy, which is generated by strain history and so is a feature of plasticity, and inherent anisotropy, associated with soil formation processes. Transverse isotropy, also called cross-anisotropy, is when there is an axis of rotational material symmetry in the material behaviour (e.g. Atkinson, 1975; Barden, 1963; Gibson, 1974; Graham and Houlsby, 1983; Simpson, 2017). Gazetas (1981) argued that transverse isotropy can be important in calculating the elastic settlements of footings.
Anisotropy can affect all aspects of soil behaviour (Carter, 2023), but the present paper focuses solely on small-strain responses that are conventionally termed ‘elastic’, for serviceability states at loads much lower than those expected to induce ultimate limits such as bearing capacity failure (e.g. Das, 2009).
Calculations for a transversely isotropic linear elastic (TILE) material are significantly more complicated than those for a fully isotropic linear elastic (FILE) half-space. Liao and Wang (1998) appear to have been the first to propose a general solution for stresses and displacements in a TILE half-space subject to three orthogonal point loads. Their solution method used a Hankel transform. They identified previous solutions for vertical loading, noting that some previous solutions were either partial or contained errors. They used notation that is easy to relate to subsequent work in geotechnical engineering by Lings et al. (2000) and others.
The solutions by Liao and Wang (1998) for three orthogonal loads, together with advances in soil testing for TILE properties, have the potential to improve significantly the accuracy of geotechnical calculations for settlements and more general ground motions. They may be integrated, analytically or numerically, to find displacements and stresses due to distributed loads or to loads on rigid footings, just as solutions by Boussinesq (1878), Cerruti (1884–1885) and others are used in this way for FILE materials (Davis and Selvadurai, 1996; Holl, 1941; Newmark, 1935). The methodology by Liao and Wang (1998) has also been applied to buried loads (Wang and Liao, 2002), piles (Wang, 2003; Wang and Pan, 2004; Wang et al., 2008) and strip loads (Wang, 2005).
Anyaegbunam (2014) reviewed previous work including solutions for vertical loading by Gerrard and Wardle (1973) and Liao and Wang (1998), but used similar notation to Mitchell (1900), which seems less easy to relate to modern work. An upper limit on shear modulus was proposed by Raymond (1970) on the basis of the use by Barden (1963) of the solution by Mitchell (1900) for point vertical loading. Anyaegbunam (2014) argued that this limit is theoretically invalid. This finding is supported by data examined in the present paper.
Marmo et al. (2020) presented an updated review, including Liao and Wang (1998), and proposed a completely different analysis based on a theorem by Almansi (1899). Their results were presented in terms of a differential operator weighted by a quasi-potentials or their derivatives. They require further processing to obtain relationships to notation in geomechanics. They presented fascinating graphics of the stresses due to lateral load on the surface of an E-glass/epoxy TILE half-space.
Dean (2019) explored some of published data on the transverse isotropy of soils. Measurement procedures and data for geomaterials tested in modified triaxial, hollow cylinder and other devices were described by Kuwano (1999), Sadek et al. (2007), Fioravante et al. (2013), Ratananikom et al. (2013), Wang et al. (2019), Liu et al. (2020), Shi et al. (2021), Gu et al. (2022), Zuo et al. (2022) and others. The present paper is illustrated using data by Wang (2002a, 2002b) for geological materials found at depths of up to several kilometres. These depths can be relevant to geophysical interpretation (e.g. Kearey et al., 2002, 2009), deep mining and deep waste and other facilities. Tectonic processes may move materials to or from the surface over geological timescales.
Constitutive properties: a brief review
Notation
Figure 1 shows Cartesian (x, y, z) and cylindrical coordinates (r, θ, z) and the sign conventions for directions of changes in stress that are positive in this paper. The origin is on the surface of the half-space, at the point at which loads are applied. The z-axis points downwards into the half-space. The surface is at z = 0. The z-axis is also the axis of material symmetry. R = (r 2 + z 2)1/2 is the shortest distance from the origin to a point at (r, θ, z).
Cartesian and cylindrical coordinate frames and directions for positive stresses. The origin is the point of loading. The flat horizontal surface of the half-space is at z = 0. The z-axis is the axis of material symmetry
Cartesian and cylindrical coordinate frames and directions for positive stresses. The origin is the point of loading. The flat horizontal surface of the half-space is at z = 0. The z-axis is the axis of material symmetry
Strain is assumed to be sufficiently small that second-order effects can be neglected. Except in the discussion of strain energy, stress symbols such as σ zz refer to changes in stress from an initial state at which strain is taken to be zero. Principal stresses and principal strains are taken positive in compression. Following Lings (2001), symbols superscripted ‘u’ (such as ) are used only for total stress and for constant-volume, undrained conditions. Symbols without such a superscript are for effective, drained conditions but can be regarded as general where convenient.
Constitutive equations for TILE materials
The present paper uses the notation by Lings et al. (2000), but with occasional exceptions. The relation between normal strains ϵ and changes in effective normal stresses σ in cylindrical coordinates with the axis of symmetry in the vertical, z-direction is expressed as
E h and E v are non-negative normal moduli and are equal for the special case of an FILE material. Poisson’s ratios μ hh and μ vh are equal for an FILE material. Lings et al. (2000) and others use μ hv/E h instead of the equivalent μ vh/E v in the third row. However, the reciprocity relations by Onsager (1931) imply that the matrix is symmetric, and the present paper follows Gibson (1974) and others by incorporating the symmetry automatically in the notation. For shear strains and changes of shear stresses
The two non-negative shear moduli G hh and G vh are equal for FILE materials. A combination of tensor rules and material symmetry about the z-axis gives
This reduces to a familiar relation for FILE materials (e.g. Davis and Selvadurai, 1996). For general TILE materials, there is no analogous relation for G vh.
Some authors express the compliance equations in a 6 × 6 matrix in which the 3 × 3 matrix of Equation 1 is a part and the shear equations in Equations 2–4 are also represented. The determinant of the 6 × 6 matrix is , where Δ = Δ1Δ2/E v is the determinant of the 3 × 3 matrix and
The inverse of the 3 × 3 matrix is
Liao and Wang (1998) presented the point load solutions that are the basis of proposals in the present paper. They used A ij for the 6 × 6 matrix of constrained moduli. The above is represented as
with A 12 = A 11 − 2A 66. The shear moduli are A 44 = A 55 = G vh and A 66 = G hh. All other A ij are zero. Some authors use C ij instead of A ij.
Strain energy per unit volume
Gibson (1974) proposed that the strain energy of a material must be non-negative and zero only when the stress is zero. This is consistent with conventional continuum mechanics (e.g. Spencer, 1980). For infinitesimal strains and if a TILE material modelling a soil under static conditions was initially under zero stress, the strain energy W per unit volume is the work done on a unit volume to reach the current state of stress, so that
where in this equation the stress symbols refer to stresses that are zero when the strain is zero. The sum is for i = zz, rr and θθ. Using Equation 1 gives
with
giving
Using Equation 8 to substitute for the changes of stress in terms of the strains gives
Recent work in lattice mechanics by Adhikari et al. (2020) and others use negative ‘equivalent’ dynamic Young’s moduli that incorporate effects other than stiffness. To be clear, therefore, it can be useful to derive carefully the constraints on ‘static’ TILE material properties of interest herein. The expressions in parentheses in the preceding equation can be considered to be independent variables whose squares cannot be negative, so the requirement by Gibson (1974) for non-negative strain energy implies that the factors outside the parentheses must be non-negative. The following can then be deduced:
to ensure that the first factor in the preceding equation is non-negative
Δ1 ≥ 0 to ensure that the second factor is non-negative
Δ2/E v ≥ 0 is then required to ensure that the third factor is non-negative
and Δ1 ≥ 0 taken together imply −1 ≤ μ hh ≤ 1 and E h ≥ 0
Δ2/E v ≥ 0 then implies , which implies E v ≥ 0
Δ2/E v ≥ 0 and E v ≥ 0 then implies Δ2 ≥ 0,
Δ1 ≥ 0 and Δ2 ≥ 0 are equivalent to limits on Young’s moduli and Poisson’s ratios attributed to Love (1927), Raymond (1970) and Pickering (1970) and explored further by Lings et al. (2000) and Lings (2001). The shear moduli must also be non-negative. Together these limits also ensure that the compliance matrix is positively definite, giving stable material behaviour.
Parameters for some geological materials
The present paper looks at data of the seismic anisotropy of sedimentary rocks published by Wang (2002a, 2002b). The data are summarised in Table 1, where the descriptions are geophysical. The ‘sands’ are classified as soft rocks and may what a geotechnical engineer might call lightly cemented sandstones. The data are presented as constrained moduli C ij (= A ij) and were analysed in various ways. The present paper looks at the data in new ways. Most of the published data are used, excepting where C 13 was not given or C 12 was found to be noticeably different from C 11 − 2C 66.
Data sources and key for data figures
| Materials | Symbol in figures | Number of samples | Sample depth: km |
|---|---|---|---|
| Brine-saturated African shales | | 5 | 2.3–2.5 |
| Brine-saturated North Sea shales | | 3 | 1.6–2.1 |
| Other brine-saturated shales | | 8 | 1.4–4.3 |
| Brine-saturated shaly coal | | 1 | 3.6 |
| Gas-saturated shaly coal | | 1 | 3.6 |
| Brine-saturated sands | | 8 | 2.4–4.1 |
| Gas-saturated African sands | | 2 | 2.4–2.5 |
| Gas-saturated tight sands | | 12 | 3.6–4.1 |
| Gas-saturated carbonatesa | | 35 | 3.3–4.3 |
| Brine-saturated carbonatesa | | 25 | 3.5–4.5 |
| Materials | Symbol in figures | Number of samples | Sample depth: km |
|---|---|---|---|
| Brine-saturated African shales | 5 | 2.3–2.5 | |
| Brine-saturated North Sea shales | 3 | 1.6–2.1 | |
| Other brine-saturated shales | 8 | 1.4–4.3 | |
| Brine-saturated shaly coal | 1 | 3.6 | |
| Gas-saturated shaly coal | 1 | 3.6 | |
| Brine-saturated sands | 8 | 2.4–4.1 | |
| Gas-saturated African sands | 2 | 2.4–2.5 | |
| Gas-saturated tight sands | 12 | 3.6–4.1 | |
| Gas-saturated carbonates | 35 | 3.3–4.3 | |
| Brine-saturated carbonates | 25 | 3.5–4.5 |
Limestone and dolomite
Source: Wang (2002b)
Figure 2 shows various aspects of the effective TILE parameters for the tabulated data. Figure 2(a) shows the relation between the effective normal modulus ratio E v/E h plotted horizontally against shear modulus ratio G vh/G hh. Both ratios are 1 for FILE materials, so the spread of values, most less than 1, is an indication of the degree of anisotropy of these soils. For clays, Pegah et al. (2021) suggested a unique relation of the following form, based on the paper by Mašín and Rott (2014):
with a = 1 and b = 1.25. There is some scatter about this curve for the data of Table 1.
Drained linear elastic properties of some cross-anisotropic geomaterials: (a) modulus ratios; (b–d) Poisson’s ratios; (e, f) anisotropy parameters. See Table 1 for data sources and key
Drained linear elastic properties of some cross-anisotropic geomaterials: (a) modulus ratios; (b–d) Poisson’s ratios; (e, f) anisotropy parameters. See Table 1 for data sources and key
Figures 2(b) and 2(c) show values of Poisson’s ratios μ hh and μ vh, plotted against the effective normal modulus ratio. Some authors report considerable difficulty in measuring the latter parameter for clays, partly due to strain dependence. Gibson (1974) showed that μ vh would be 1/2 in an undrained test, and μ hh would equal 1 − μ vh E h/E v. The latter are plotted against each other in Figure 2(d). This plot confirms that the data of Table 1 were not undrained data (see also Equations 50 and 51 later).
Measures of anisotropy of TILE materials
Figures 2(a)–2(d) indicate that E v/E h may be incomplete as an identifier of transverse isotropy. Other proposals include parameters by Thomsen (1986), Tsvankin and Thomsen (1994) and Alkhalifah and Tsvankin (1995) and others, which are focused on geophysical applications and are based on ratios of stress wave velocities.
The five independent material constants of a TILE material are reduced to two for an FILE material, so an unbiased measure of anisotropy needs to involve three components. One possibility might involve two further shear moduli from the geotechnical literature:
Both use drained material parameters and equal G hh for an FILE material. G iso is described by Lings et al. (2000) and G * by Graham and Houlsby (1983). They are one-third of the slopes of a graph of deviator stress against deviator strain in standard drained and undrained triaxial tests, respectively, on a vertically aligned TILE material (see Appendix). Three aspects of anisotropy can then be defined as follows:
Each would be between −1 and 1. A combined measure might be
If χ = 0, then the material is an FILE material. Otherwise, the material is anisotropic and may be a TILE material or have more complex anisotropy. Larger values of χ would be expected to imply greater inaccuracy if the material is modelled as fully isotropic. The largest possible value is 3.
Figure 2(e) shows χ plotted against E v/E h for the data of Table 1. The limit curve corresponds to χ 2 = χ 3 = 0. The graph indicates that χ correlates approximately but non-linearly with E v/E h and that this is 1 when χ is close to zero.
Absolute values of the three aspects might be interpreted as indicators of potential errors associated with assuming that their relevant determining moduli are equal. Figure 2(f) shows χ 2 plotted against χ sign (χ 2). Lines of slopes 0.5 and 1/21/2 correspond respectively to and . The graph indicates that the shear aspect is an important part of the overall anisotropy for the data of Table 1.
Liao and Wang (1998) results
Characteristic equation
The solutions by Liao and Wang (1998) for a TILE half-space subject to general point loads are the essential basis of the present paper. Their solutions involve dimensionless variables u i (i = 1 to 3) and m i (i = 1–2) that satisfy
Equation 22 is the ‘characteristic equation’ for TILE materials. It can be deduced by applying the physical constraints associated with internal equilibrium. It also serves as a definition of the dimensionless parameters m i. Using the equality between the two expressions for m i, Liao and Wang (1998) noted that u 1 and u 2 satisfy a quadratic with
The solutions for , for subscript i = 1 or 2, were deduced as follows:
Liao and Wang (1998) considered three cases, depending on the sign of s 2 − 4q. Sheet A4 in the online supplementary material verifies that u 1 and u 2 are either both real and positive for surface loads or else are complex conjugates of each other. The present paper takes the signs in the preceding expression such that .
Some new results
Using the correspondence between Equations 8 and 9 to substitute for the A ij in the expressions by Liao and Wang (1998) for s and q gives
Together with Equation 26, these give a way of calculating u 1 and u 2 without needing to invert the compliance matrix. Using the correspondence again, together with the second expression for m i in the characteristic equation, gives
This gives a way of calculating m i without needing to invert the compliance matrix. Using the first expression for m i in the characteristic equation to compute the product m 1 m 2 gives
Using Equations 24 and 25 to simplify the denominator gives
Liao and Wang (1998) defined a variable k with dimensions of inverse modulus. Its definition may be conveniently re-expressed using a new dimensionless variable k′ as follows:
with
where the second expression is deduced using Equations 6–9 and 26. It follows that k′ is either wholly real or wholly imaginary. Liao and Wang (1998) also defined four variables T, two of which are related as follows:
Using Equation 31 to simplify the definitions of the other two gives
Hence, T 2 = T 3. The following new dimensionless variables will be found useful:
with
λ 1 and λ 2 are real if u 1 and u 2 are real, and complex conjugates otherwise.
Effective properties of some geological materials
Figure 3 shows some of these results for the data of Table 1. Figure 3(a) shows that s is positive for all of these data and in most cases greater than 2. Figure 3(b) shows that s 2 − 4q can be positive or negative for these data. Raymond (1970) proposed that the shear modulus G vh would be limited by a value here denoted as G R (Lings et al., 2000):
Anyaegbunam (2014) asserted that this limit is invalid. Using Equations 27 and 28, the criterion G vh ≤ G R is found to be equivalent to insisting that s 2 − 4q ≥ 0. The data show that this is not satisfied by real materials. The limit by Raymond (1970) is further discussed later in this paper.
Drained point load properties: (a, b) s and s 2 − 4q; (c, d) u 1 and u 1; (e, f) m 1 and m 2. See Table 1 for data sources and key
Drained point load properties: (a, b) s and s 2 − 4q; (c, d) u 1 and u 1; (e, f) m 1 and m 2. See Table 1 for data sources and key
Figures 3(c) and 3(d) shows the solutions for u 1 and u 2 for the carbonates of Table 1. In Figure 3(c), the real parts of u 1 and u 2 are plotted on the vertical axis. When s 2 − 4q > 0, the real parts are greater for u 1 than for u 2. When s 2 − 4q < 0, the real parts are equal because then u 1 and u 2 are complex conjugates. Figure 3(d) shows the imaginary coefficients. The values are zero when s 2 − 4q > 0, and equal and opposite when s 2 − 4q < 0.
Figures 3(e) and 3(f) show corresponding results for m 1 and m 2. Real values of both m i and u i are approximately symmetric about the value of 1, which is the value for FILE materials. This is consistent with the values of s for the carbonates in Figure 3(a) being close to 2.
Surface loading: general cases
Motivation and notation
Surface loading is an application of the equations by Liao and Wang (1998) with practical relevance, since most soils are anisotropic and many calculations in practical engineering are for loaded areas of ground or for rigid footings on or close to the surface. However, those equations occupy a total of over three densely packed pages of their paper. Consequently, a first motivation for exploring them is to find simplifications. A second motivation is that the present author noticed what is believed to be a notational confusion and a typo in one of their equations, and these need to be clarified and avoided.
Regarding simplifications, the equations involve eight different adjusted z-coordinates, eight associated adjusted distance coordinates and eight other distances. Calculation sheet A2 in the online supplementary material shows that these 24 variables reduce to three groups of three for the special case of surface loading. In each group, the subscripts are 1, 2 and 3, with
The factors u i are 1 for FILE materials, and the parameters then reduce to z, R and r 2/(R + z). These are familiar from solutions by Boussinesq (1878), Cerruti (1884–1885) and others in the theory of fully isotropic elasticity (e.g. Davis and Selvadurai, 1996).
Regarding the notational confusion and typo, this is explained in the section headed ‘Lateral surface loading’ and is detailed in the online supplementary material calculation sheets.
Vertical surface loading
Calculations for vertical loading are common in geotechnical practice, and the simplest is for surface loading. The solution by Boussinesq (1878) is the basis for calculations for a surface footing on an FILE material, so a corresponding calculation for a TILE material is naturally of interest.
Liao and Wang (1998) present their results for displacements and stresses in two parts, which are combined to form the complete solution. Combining their Equations 47–49 with Equations 73–75, considering only a vertical point load P z at the origin of coordinates and applying the simplifications for surface loading gives
Using Equations 31–39, these simplify as in (a) in Table 2. Details are provided in sheets B2–B4 in the online supplementary material. Circumferential displacement U θ = 0 is zero by symmetry.
Results for vertical point load at the surface
| T2.1 | Note | and |
| (a) Displacements | ||
|---|---|---|
| T2.2 | Radial | |
| T2.3 | Circumferential | Uθ = 0 |
| T2.4 | Vertical | |
| (b) Strains | ||
| T2.5 | ||
| T2.6 | ||
| T2.7 | ||
| T2.8 | ||
| T2.9 | ||
| (c) Stresses | ||
| T2.10 | ||
| T2.11 | ||
| T2.12 | ||
| T2.13 | τrθ = τθr = 0 | |
| T2.14 | ||
| T2.1 | Note | |
| (a) Displacements | ||
|---|---|---|
| T2.2 | Radial | |
| T2.3 | Circumferential | Uθ = 0 |
| T2.4 | Vertical | |
| (b) Strains | ||
| T2.5 | ||
| T2.6 | ||
| T2.7 | ||
| T2.8 | ||
| T2.9 | ||
| (c) Stresses | ||
| T2.10 | ||
| T2.11 | ||
| T2.12 | ||
| T2.13 | τrθ = τθr = 0 | |
| T2.14 | ||
Based on the general solutions by Liao and Wang (1998)
Using standard equations for strains in cylindrical coordinates, which Liao and Wang (1998) quoted, gives the equations for strains in (b) in Table 2. Applying the constitutive equations for a TILE material then gives the equations for changes in stress in (b) in Table 2. These are found to be consistent with Equations 50–55 and 76–81 by Liao and Wang (1999). For completeness, sheets F1 and F2 in the online supplementary material confirms that the stresses in (c) in Table 2 also satisfy the standard equilibrium equations in cylindrical coordinates.
Raymond’s proposed limit
Raymond (1970) argued that the solution by Mitchell (1900) for vertical point loading, as quoted by Barden (1963), would produce unphysical complex-valued vertical stresses at depths below the load if G vh exceeded G R. In the equations by Liao and Wang (1998), r = 0 implies R i = u i z, so from the equation for vertical stress in Table 2
If λ i, m i and u i are real numbers, the summation is real. Reference to Equations 31–39 shows that if u 1 and u 2 are complex conjugates, each of the two terms in the preceding summation is the complex conjugate of the other, so their sum is real. It can readily be checked that this applies for all the stresses and at all positions in the FILE half-space. Hence, the solutions by Liao and Wang (1998) do not produce complex values of this stress. Raymond’s argument does not therefore apply.
As a check, Figure 4 shows the ratio G vh/G R plotted against the effective normal modulus ratio for the data of Table 1. Figure 4(a) uses the tabulated parameters directly, and Figure 4(b) shows undrained parameters inferred from the drained ones using algebra described later in this paper. The results indicate that Raymond’s proposed limit was exceeded in many of the tests. This seems unlikely to be explicable as experimental error. On these theoretical and experimental grounds, in agreement with Anyaegbunam (2014), and unless some other argument is found, it is proposed that the limit on G vh by Raymond (1970) be considered invalid.
Investigation of the proposed limit by Raymond (1970) on G vh for some sedimentary rocks: (a) drained parameters; (b) constant-volume parameters. See Table 1 for data sources and key
Investigation of the proposed limit by Raymond (1970) on G vh for some sedimentary rocks: (a) drained parameters; (b) constant-volume parameters. See Table 1 for data sources and key
Lateral surface loading
Lateral loading is relevant in many geotechnical applications, including when lateral loads from wind, wave and current forces are transmitted into foundations (e.g. Cassidy et al., 2004; Dean, 2009; Randolph and Gourvenec, 2011), or when transport or machine vibrations occur, or seismic events that often involve strong lateral shaking (Kramer, 1996; Srbulov and O’Brien, 2012).
The equations by Liao and Wang (1998) for displacements and stresses use symbols P r and P θ described as loads in the radial and circumferential directions, respectively. However, the present author noticed what appeared to be a notational confusion. For a lateral point load applied in some general azimuthal direction θ P, physical symmetry would indicate that displacements would have reflection symmetry in the plane θ = θ P. Vertical and radial displacements would have a symmetry like cos(θ − θ P), and circumferential displacements like sin(θ − θ P). However, the equations by Liao and Wang (1998) indicate that the responses to P r exhibit symmetries of cos θ and sin θ, respectively, indicating that θ P = 0. This is considered herein to imply that P r can be interpreted as a point load P x in the +x direction. Similarly, P θ is found to be a point load P y in the +y direction.
Additionally, a possible typo was noticed in one of the equations for stress. To resolve these issues, the equations by Liao and Wang (1998) for stress for the lateral loading case were ignored. Instead, their equations for displacements were assumed to be correct, and the equations for strains and changes in stress were then calculated explicitly. The results are in Table 3 and were further checked by checking equilibrium. These calculations are included in calculation sheets E1 to G4 in the online supplementary material and are briefly explained in the following.
Results for horizontal point loads at the surface
| (a) Displacements | ||
|---|---|---|
| T3.1 | Radial | |
| T3.2 | Circumferential | |
| T3.3 | Vertical | |
| (b) Strains | ||
| T3.4 | ||
| T3.5 | ||
| T3.6 | ||
| T3.7 | ||
| T3.8 | ||
| T3.9 | ||
| (c) Stresses | ||
| T3.10 | ||
| T3.11 | ||
| T3.12 | ||
| T3.13 | ||
| T3.14 | ||
| T3.15 | ||
| (a) Displacements | ||
|---|---|---|
| T3.1 | Radial | |
| T3.2 | Circumferential | |
| T3.3 | Vertical | |
| (b) Strains | ||
| T3.4 | ||
| T3.5 | ||
| T3.6 | ||
| T3.7 | ||
| T3.8 | ||
| T3.9 | ||
| (c) Stresses | ||
| T3.10 | ||
| T3.11 | ||
| T3.12 | ||
| T3.13 | ||
| T3.14 | ||
| T3.15 | ||
Deduced from the displacements in the general solutions by Liao and Wang (1998)
Liao and Wang (1998) solved the equations by considering two separate analyses whose combination would give the results sought. Combining their Equations 47–49 with Equations 73–75, setting P r = P x and P z = P θ = 0 and applying the simplifications associated with surface loading gives
Using Equations 31–39, these simplify as in (a) in Table 3. Details are given in sheets B5–B7 in the online supplementary material. Applying the standard equations for infinitesimal strains in cylindrical coordinates, as quoted by Liao and Wang (1998) and elsewhere, gives the strains in (b) in Table 3. Using material behaviour Equations 1–4 gives the stresses in (c) in Table 3. Details are in sheets C2–C7, D6 and D7 in the online supplementary material.
Except for the shear stress τ rz, the results are found to be consistent with Equations 50–55 and 76–80 by Liao and Wang (1998). The present results would be fully consistent if the first sign on the right of their Equation 79 was changed from plus to minus. Checks on equilibrium are presented in sheets F3–F5 in the online supplementary material. They use the same standard equilibrium equations in cylindrical coordinates as those quoted by Liao and Wang (1998) and elsewhere. They confirm that the equations of Table 3 do satisfy equilibrium. Hence, these equations are judged to be correct.
Special cases 1: constant-volume behaviour
Motivation and notation
The constant-volume condition is considered equivalent to the undrained condition if the soil particles are incompressible and Terzaghi’s principle of effective stress applies (Bowles, 1996; Dyvik et al., 1987; Knappett and Craig, 2012). The condition therefore reduces the degrees of freedom available in terms of total stress parameters, and this requires special treatment.
Undrained properties and behaviours of TILE materials were well explored by Gibson (1974), Lings (2001) and others, including empirical proposals for estimating the drained properties from the smaller number of independent undrained parameters (Pegah et al., 2021). For present purposes, following Lings (2001), an undrained parameter is signified by superscript ‘u’ and a change in the total stress is also denoted by a superscript. With this notation, two key relationships given by Gibson (1974) may be stated as follows:
Reference to Equation 7 reveals that these imply , which implies that the undrained version of the compliance matrix of Equation 1 is not invertible. The equations by Liao and Wang (1998) use the inverse, so a workaround is developed in the following to enable the equations to be used for undrained conditions.
Calculations for undrained material parameters
As noted earlier, these calculations are well established in the literature. The following can be a useful alternative approach. Changes in pore stress are herein denoted as ‘u’, and the absence of a subscript on this symbol distinguishes it from displacements and from parameters u i by Liao and Wang (1998). Terzaghi’s principle of effective stress implies that changes σ of effective stress in Equation 1 can be replaced by differences σ u − u (e.g. Knappett and Craig, 2012). When this is done, application of the constant-volume condition ϵ rr + ϵ θθ + ϵ zz = 0 gives
with
These imply 2α + β = 1. The condition Δ2 ≥ 0 implies that γ is non-negative. Then, changes in effective stress are given by
Using this to substitute for the changes in effective stress in Equation 1, a new compliance equation is obtained with the drained parameters replaced with undrained ones as follows:
These results are consistent with expressions given by Lings (2001) and others. Equations 57 and 60 ensure that the constant-volume moduli and are greater than the corresponding drained moduli. Shearing is unaffected by the constant-volume condition, so and , and an undrained version of Equation 5 applies.
Results for undrained material parameters
Figure 5 explores undrained parameters for the TILE data of Table 1. Figure 5(a) shows that the parameter α does not stray far from its value of 1/3 for FILE materials, suggesting that these materials are only weakly anisotropic. Figure 5(b) indicates that there is an approximate correlation between α and the effective anisotropy parameter χ.
Undrained linear elastic properties: (a, b) α; (c–e) ratios; (f) anisotropy parameters. See Table 1 for data sources and key
Undrained linear elastic properties: (a, b) α; (c–e) ratios; (f) anisotropy parameters. See Table 1 for data sources and key
Figures 5(c) shows that the undrained normal ratio clusters around 1, in contrast to the drained ratio E v/E h. Gibson (1974) quotes the result ‘E H ≤ 4E V’ by Ferrar (1941) for an incompressible orthotropic elastic medium. Using Equations 6 and 7 with undrained parameters gives
Hence, the expression on the left is non-negative. Putting and then gives the condition for undrained TILE materials by Ferrar (1941). Figure 5(c) confirms that the data of Table 1 satisfy this.
Figure 5(d) shows that is more than 1/2 for some of these materials. The data points here form a pattern that is rather similar to that in Figure 5(c). This is because of Gibson (1974) relations, Equations 50 and 51.
Liao and Wang (1998) equations for the undrained condition
Applying Equations 50 and 51 by Gibson (1974) to Equations 27 and 28 gives the undrained values of Liao and Wang (1998) parameters as
This implies that s u2 − 4q u is only zero when , which is represented by the line of equality in Figure 5(e). s u2 − 4q u is negative for points above the diagonal line of equality, and positive for points below it. This also confirms that undrained parameters u i and m i can be calculated even though the undrained compliance matrix is singular. Equation 29 reduces to
Using this with and with Gibson (1974) relations and Equations 33 and 37–39, and taking Equation 63 to imply , gives
The denominator in the last two equations is zero in the case of full isotropy, which is analysed later.
Equations 16 and 17 can be written using constant-volume parameters. In the latter, dominates , so that . Hence, the triaxial aspect of anisotropy is zero for the undrained condition, Equation 20, although the shear and normal aspects are not. Figure 5(f) shows the relation between these aspects. An FILE material would plot at the origin in this diagram, and a curve of constant χ u plots as a circle centred on the origin.
Special cases 2: full isotropy
Motivation and methodology
Barden (1963: p. 203) stated that ‘isotropy is simply a special case of anisotropy’. It can therefore be useful to check that the TILE equations of Tables 2 and 3 do reduce to known solutions for the FILE materials for appropriate material constants. However, the following problem is found. For an FILE material, Young’s moduli are equal, Poisson’s ratios are equal and the shear moduli are equal. The solutions by Liao and Wang (1998) give s = 2 and q = 1, so s 2 − 4q = 0 and u 1 = u 2 = m 1 = m 2 = 1. This makes k and k′ infinite (Equations 32 and 33) and renders λ 1 and λ 2 indeterminate (Equations 37 and 38). Special treatment is therefore needed for this material.
To find a workaround, the method by Barden (1963) of taking a limit is adapted herein. A TILE material is considered that is infinitesimally close to fully isotropic. There are several ways in which a TILE material, with five independent constants, can be close to being fully isotropic, with two. In the present paper, a material is considered with E v = E h = E and μ vh = μ hh = μ and with the following:
where ϵ is an infinitesimally small number. This will tend to become an FILE material when ϵ tends to zero. Since s 2 − 4q = 0 for the FILE material, it will be infinitesimal for the nearly FILE material, and parameters with (s 2 − 4q)1/2 in the denominator tend to infinity as ϵ → 0. Details are given in sheets H1–H3 in the online supplementary material. Key results are quoted in the following.
Vertical loading
Boussinesq (1878) developed the solution for a point vertical load on the surface of an FILE half-space. The solution is explored by Westergaard (1952), Davis and Selvadurai (1996), Podio-Guidugli and Favata (2014) and others. Material displacements are given by
where R = (r 2 + z 2)1/2 and G is the isotropic shear modulus. Displacements for a TILE material for this case, deduced from the paper by Liao and Wang (1998), are given by the equations in (a) in Table 2. Applying results from sheets H1–H3 in the online supplementary material gives
with η = z 2/R 2, where O(ϵ) represents terms of order ϵ and smaller. Here and below, the ± and ∓ signs are such that the upper sign is for i = 1 and the lower for i = 2. Consequently, the 1/ϵ terms cancel in the summations of (a) in Table 2. Taking the limit as ϵ → 0 then gives
and U θ = 0. Putting η = z 2/R 2 then gives Boussinesq (1878) equations, as required.
Horizontal loading
Cerruti (1884–1885) developed the solution for a point horizontal load on the surface of an FILE half-space. The solution is explored by Westergaard (1952), Davis and Selvadurai (1996), Podio-Guidugli and Favata (2014) and others. Material displacements are given in cylindrical coordinates by
Displacements for a TILE material for this case are given by the equations in (a) in Table 3. Applying results from sheets H1–H3 in the online supplementary material gives
The 1/η terms cancel in the summations of (a) in Table 3. Taking the limit as η → 0 then gives
Putting η = z 2/R 2 and setting u 3 = 1 for the FILE material (Equation 23) then gives the solution by Cerruti (1884–1885) for this case, as required.
Serendipitous result concerning accuracy
The preceding calculations can also be interpreted as an assessment of the consequences of a small inaccuracy in the assessment of the shear modulus G vh. Now ϵ 2 appears in Equation 69, but only terms of order ϵ resulted in the penultimate equations after the terms involving 1/ϵ cancel. Hence, small errors can have larger effects. For example, an error of order ϵ 2 = ±0.01 in G vh produces an error of order ϵ = ±0.1 in the final displacements. This suggests that, while Barden (1963) is correct that isotropy is a special case, an FILE material may also be singular in some sense.
Compliance factors for practical applications
Motivation
Davis and Selvadurai (1996) outlined how the solution by Boussinesq (1878) for point loads on an FILE half-space can be used to calculate displacements of loaded areas of ground or stiffnesses and interface stress distributions for rigid footings. Similar techniques can evidently be used with the point load solutions listed in Tables 2 and 3 for TILE materials. The concept of ‘compliance factors’ proposed in the following section is intended to facilitate these calculations.
Compliance factors: theory
Equations 72 and 71 and 79 and 78 gave the displacements respectively for the solution for a point load by Boussinesq (1878) and the solution for a horizontal load by Cerruti (1884–1885), both applied to the surface of the TILE half-space. The surface has z = 0, implying R = r (Equation 42). Hence, for combined loading, displacements (u r, u θ, u z) at a point (r, θ) on the surface can be expressed as
where G vh = G for the FILE material, and the compliance factors C ij are deduced directly from the displacement equations and are listed in Table 4.
Compliance factors for surface displacements
| Factor | Description | FILE material | TILE material | |
|---|---|---|---|---|
| General | Constant volume | |||
| Czz | Vertical displacement due to vertical force | 1 − μ | ||
| Czx | Vertical displacement due to horizontal force | 0 | ||
| Crz | Radial displacement due to vertical force | 0 | ||
| Crx | Radial displacement due to horizontal force | 1 | ||
| Cθz | Circumferential displacement due to vertical force | 0 | 0 | 0 |
| Cθx | Circumferential displacement due to horizontal force | μ − 1 | ||
| Factor | Description | FILE material | TILE material | |
|---|---|---|---|---|
| General | Constant volume | |||
| Czz | Vertical displacement due to vertical force | 1 − μ | ||
| Czx | Vertical displacement due to horizontal force | 0 | ||
| Crz | Radial displacement due to vertical force | 0 | ||
| Crx | Radial displacement due to horizontal force | 1 | ||
| Cθz | Circumferential displacement due to vertical force | 0 | 0 | 0 |
| Cθx | Circumferential displacement due to horizontal force | μ − 1 | ||
For any known loads on a surface, once the compliance factors are known, the preceding equations are all that are needed to determine surface displacements, by integration over the loaded area. Also, for any given displacements of a rigid footing, these equations and the interface friction conditions are all that are needed to determine interface stresses and overall footing stiffnesses. Details of what happens inside the half-space are not needed for such determinations, as long as the half-space is known to be either an FILE or a more general TILE material. Consequently, the factors may be of particular interest in calculations using the boundary-element method described by Brebbia (1978), Banerjee and Butterfield (1981), Katsikadelis (2016), Zhou et al. (2109) and others.
Undrained case
Using these results with the undrained results of Equation 68 gives . This implies that there is no interaction between vertical and lateral directions at the surface, which is also the case for an FILE material with μ = 1/2. Using Equations 65–68 gives
Multiplying top and bottom on the right by and using the undrained version of Equation 26 then gives the expression for the constant volume in Table 4. The result for follows from and .
Compliance factors: values for materials of Table 1
Figures 6(a)–6(c) show the values of three of the compliance factors for drained conditions for the data of Table 1. While some approximate trends are visible, there are no clear correlations. The values are mostly of the same order of magnitude as would be obtained for an FILE material with μ in the range 0–0.5. Figure 6(d) shows the directional coupling factor C zx plotted vertically against C zz. The line marked FILE is the relationship for FILE materials, and the data tend to cluster around this line. These results may suggest that an engineering approximation of TILE materials as FILE materials may not be too inaccurate as far as surface loading is concerned, provided an accurate estimate of G vh is used.
Compliance factors: (a–d) drained; (e) undrained; (f) compared. See Table 1 for data sources and key
Compliance factors: (a–d) drained; (e) undrained; (f) compared. See Table 1 for data sources and key
Figure 6(e) shows values of plotted vertically against the undrained normal modulus ratio. The data for the saturated carbonates and for the coal seem to have similar trends. Figure 6(f) shows plotted against the drained value. The plot shows that for all of the same of Table 1, indicating stiffer vertical responses in the undrained condition compared with those in the drained condition.
Discussion and concluding remarks
Achievements of this research
It has been well established that soils are anisotropic and that the form of anisotropy known as transverse isotropy is a common form. This paper has made use of the equations by Liao and Wang (1998) for TILE materials. These equations, together with advances in soil testing, have the potential to improve significantly the accuracy of geotechnical calculations for ground movements in small-strain linear elastic contexts.
Specifically, this paper has applied the equations by Liao and Wang (1998) to the special case of point loading on the surface of a TILE half-space. Many simplifications were achieved, making the results more accessible. Equations for the special case of undrained loading were also explored, showing that this is more complex than for fully isotropic materials. The present paper has not disagreed with the suggestion by Barden (1963) that isotropy is simply a special case of anisotropy, but it has also shown that fully isotropic material is in one sense a singular material.
Anyaegbunam (2014) tried to make the equations by Liao and Wang (1998) more accessible for vertical loading at any point below the surface of a TILE half-space. The present paper has hopefully made them accessible for vertical and lateral loadings at a point on the surface. The paper has clarified the physical meanings of symbols used by Liao and Wang (1998) for lateral loads and has identified (and avoided) a suspected error in one of their equations. Details are given in sheet E2 in the online supplementary material.
Practical applications of this research
For a linear elastic material, the effects of distributed loads over a surface can be calculated by suitable integration of effects of infinitesimal point loads. The results in Tables 2 and 3, and the associated compliance factors in Table 4, can be readily applied to such calculations for settlements and other ground movements for transversely isotropic soils. In particular, the compliance factors can be convenient for the boundary-element method. Results of such calculations would need to be interpreted within limitations outlined in the following.
It has been shown that small errors in shear modulus G vh for materials close to fully isotropic can produce larger errors in displacement predictions. In agreement with Anyaegbunam (2014), the present calculations indicate that the proposed limit on G vh by Raymond (1970) is invalid. This is also indicated by the data on soft rocks examined in the paper.
Limitations of this research
The primary practical limitations of the present work are that, firstly, the elastic properties of soils are not straightforward (Castellón and Ledesma, 2022; Jardine, 1992) and behaviour can be non-linear even at small strains (Atkinson, 2000; Houlsby et al., 2005). Consequently, elastic parameters must be selected carefully based on expected strain levels. The solutions here are limited to surface loading of a laterally and vertically uniform TILE half-space, with a flat surface and with the axis of material symmetry at right angles to the surface. More complex geometries exist, including layered soils and sloping strata. More complicated forms of anisotropy exist (Dean, 2019; Gibson, 1974).
Supplementary data availability
A supplementary file is available in PDF format as online supplementary material accompanying this article. The file presents calculations that support those in the main text and tables. All the data presented on figures in the paper are from tables in the paper by Wang (2002b). The bespoke Excel VBA software used to process the data for presentation here can be made available for research purposes by request to the author.
Appendix: shear moduli in standard triaxial tests
The purpose of this appendix is to determine the physical interpretations of the shear moduli G dtx and G utx of Equations 16 and 17.
The theory of the responses of TILE materials in the standard triaxial test are well documented – for example, in the papers by Barden (1963), Gibson (1974), Atkinson (1975), Graham and Houlsby (1983). If the axis of material symmetry of the sample aligned with the axial direction of the cell, then z corresponds to this direction and x and y are radial directions. Denoting changes in the axial and radial effective stresses as Δσ ax and Δσ rad and corresponding strains ϵ ax and ϵ rad, the compliance equation in Equation 1 for the TILE material then gives
(e.g. Nishimura and Magalona (2020), with μ hv/E h = μ vh/E v). A change Δq in deviatoric stress, a volume strain ϵ vol and a deviatoric strain ϵ q can be defined as (e.g. Schofield and Wroth, 1968).
Graham and Houlsby (1983) showed that the change of deviator stress is related to strains as follows:
G * can be expressed as in Equation 17. Hence, G * = G utx is one-third of the slope of the plot of deviatoric stress against deviatoric strain in an undrained, constant-volume test. In a drained test with constant radial stress, volume strain occurs but Δσ = 0, and Equation 90 implies
where G iso = G dtx is given by Equation 16. Hence, G dtx is one-third of the slope of the plot of deviatoric stress against deviatoric strain in this drained (constant radial stress) type of test.







