Exact Winkler Solution for Laterally Loaded Piles in Inhomogeneous Soil

Crispin, Jamie J.; Mylonakis, George

doi:10.1061/(ASCE)EM.1943-7889.0002125

Open access

Technical Papers

Sep 6, 2022

Exact Winkler Solution for Laterally Loaded Piles in Inhomogeneous Soil

Authors: Jamie J. Crispin, Ph.D. https://orcid.org/0000-0003-3074-8493 [email protected], and George Mylonakis, Ph.D., M.ASCEAuthor Affiliations

Publication: Journal of Engineering Mechanics

Volume 148, Issue 11

https://doi.org/10.1061/(ASCE)EM.1943-7889.0002125

PDF

Abstract

A novel exact analytical solution is derived for the equation

y^{(4)} + x^{n} y = 0

in the region

x \geq 0

, which is important for the analysis of piles in soil with stiffness varying with depth. To date, exact solutions for long piles are available only for the cases where

n = - 4

, 0, and 1. For other values of the exponent

n

, solutions have formerly been obtained numerically, mainly by the finite-difference method or approximate analytical solutions. An inherent difficulty in obtaining solutions for long beams (which are used to model flexible piles) lies in the inability to isolate the regular, converging part of the solution over the singular part that diverges with increasing

x

. In this paper, an exact solution is derived for

n > - 4

, focusing on the important case of semi-infinite beams. Key aspects of the problem such as the stiffness and flexibility matrices at the pile head, and the peak bending moments due to eccentrically acting lateral loads, are discussed. A novel approach for deriving Winkler spring moduli for combined force and moment loading is proposed and shown to provide good agreement with rigorous numerical continuum results.

Introduction

The problem of a loaded beam supported by an inhomogeneous Winkler spring bed has been studied for over 100 years. The earliest analytical solutions known to the authors were provided by Hayashi (1921), who considered a linear inhomogeneity function for the variable Winkler modulus. The problem is still relevant in the modern day, where solutions can be employed to determine the deflection of a pile in an inhomogeneous soil under a lateral head load (Chang 1937), with solutions for power-law soil stiffness profiles being of particular utility. The problem is expressed by the equation

y^{(4)} + x^{n} y = 0

, where

x

(

\geq 0

) and

y

are the independent and dependent variables, respectively, and

n

is a pertinent inhomogeneity parameter. Contrary to the corresponding problem of axial pile response in inhomogeneous soil which is straightforward to solve (Scott 1981; Guo 2012; Crispin et al. 2018), the problem of lateral pile response is harder to tackle. The underlying difficulties in obtaining analytical solutions lie in (1) the high order of the governing differential equation, (2) the variable coefficient (often associated with a noninteger power

n

), and (3) the divergent behavior of the solution at infinity.

Hayashi (1921) provided two solutions for

n = 1

. The first employed a Laplace transform to obtain an expression in an integral form. However, due to the complexities associated with evaluating (analytically or numerically) this integral, a second solution was provided as a formal power series. Similar power-series solutions have been reported by Rifaat (1935) and Hetenyi (1946); however, other than a similar solution attributed to Miche (1930) by Rifaat (1935), a workable integral solution was not provided prior to the seminal paper by Franklin and Scott (1979) using contour integration. Instead, numerical finite-difference schemes available in graphical form (Palmer and Thomson 1948; Gleser 1953; Matlock and Reese 1960) have been widely adopted in books and design manuals (Lam and Martin 1986; Reese and Van Impe 2011).

Despite the importance of the problem in engineering theory and practice and the intense research in the subject, until recently no significant progress in obtaining exact analytical solutions for

n \neq 1

has been achieved. The special case of

n = - 4

has recently been solved by Froio and Rizzi (2016), who provided a constructive discussion on the physics of the solution. However, this is an inhomogeneous bed of limited practical importance for the pile problem considered here. Approximate analytical solutions obtained for arbitrary values of

n

derived by means of asymptotic methods and Wentzel-Kramers-Brillouin (WKB) perturbation series (Franklin and Scott 1979; Liang et al. 2014; Fradelos 2016) provide insight as to the attenuation of the solution with depth, yet do not offer dependable results for the stiffness and deflection at the pile head. Such results are available using the aforementioned numerical approaches or approximate analytical techniques based on energy principles and/or multilayer approximations of inhomogeneous soil (Mylonakis 1995; Mylonakis and Roumbas 2001; Karatzia and Mylonakis 2012, 2017). A summary of relevant energy methods and some refinements pertaining to dynamic loads is available in a recent review article by Mylonakis and Crispin (2021).

In recent years, interest in analytical treatment of the problem has re-emerged and some general analytical Winkler solutions encompassing arbitrary values of

n

have been formulated in terms of four generalized hypergeometric functions of the

{}_{0}{F_{3}}

type (Mylonakis 2004; Fradelos 2016; Crispin 2017; Froio and Rizzi 2017; Parashakis 2020; Rosendo and Albuquerque 2021; Cucuzza et al. 2022). These results extended the aforementioned classical power-series solutions, yet they are limited by their inability to handle the important case of long beams, which is used to model flexible piles. The problem lies in the divergent (singular) behavior of all four hypergeometric functions and the difficulty in extracting nondivergent (regular) solutions that tend to zero at infinity in the form of a linear combination of the four functions. Consequently, only piles of finite length can be treated, and specific boundary conditions at the pile tip should be considered. This limitation significantly complicates the solution and often leads to numerical instabilities for long piles (Fradelos 2016), providing the motivation for the herein reported work.

In the following, the problem is solved in an exact manner, with a focus on the important case of long piles (treated as semi-infinite beams) and arbitrary values of the inhomogeneity exponent

n

. To this end, the following steps have been undertaken:

1.

The problem is formulated with the specific power-law inhomogeneity function defined and a novel definition of the Winkler parameter

λ

introduced. The desired form of the solution is detailed.

2.

The existing solution for

n = 0

(homogeneous soil) is rederived to introduce the solution procedure used subsequently for long piles.

3.

Two existing solutions for

n = 1

(linear soil profile) are compared and combined to develop a novel solution for long piles.

4.

The power-series solution for arbitrary

n

(power-law soil profile) and short piles are rederived and extended to long piles through careful consideration of the asymptotic properties of the series.

5.

Application of the Winkler solutions in Points 2, 3, and 4 to the pile problem is considered, including comparisons with numerical continuum solutions for the homogeneous and linear soil profiles.

Regarding this final point, the main challenge in employing Winkler models lies in the assessment of the moduli of the relevant springs,

k (z)

. Various approaches are available to select

k (z)

, which is usually taken to be proportional to soil stiffness,

E_{s}

(or

G_{s}

). Dependable expressions for

k (z)

from various sources have been summarized by Shadlou and Bhattacharya (2014), Anoyatis and Lemnitzer (2017), and Mylonakis and Crispin (2021). These are commonly obtained by calibrating the Winkler model against a more rigorous solution, usually a numerical continuum model implemented via finite or boundary elements. Calibrations of this type have been carried out by Biot (1937), Vesić (1961), Francis (1964), Yoshida and Yoshinaka (1972), Roesset (1980), Dobry et al. (1982), Gazetas and Dobry (1984), and Syngros (2004). Rigorous numerical continuum solutions for inhomogeneous soil are available (Banerjee and Davies 1978; Randolph 1981; Budhu and Davies 1988; Gazetas 1990; El-Marsafawi et al. 1992; Syngros 2004) that can be (and have been) used to calibrate the springs for inhomogeneous soils. However, the resulting Winkler modulus is dependent on both the soil stiffness profile as well as the pile head fixity conditions, thus requiring separate calibrations for different response parameters (i.e., head displacement, head rotation, and peak moment) and problem configurations (i.e., inhomogeneity parameters).

As an alternative, Mylonakis (2001) proposed a rational analytical approach based on a variant of the method by Vlasov and Leontiev (1966), originally developed for spread footings, that does not focus on calibrating a specific parameter. This solution was explored further by Karatzia and Mylonakis (2017), who approximated shape functions for inhomogeneous soil profiles using those determined from the Winkler solution when

n = 0

. However, this was limited to only two boundary conditions: a pile under a head force without head rotation (fixed head) and a pile under a head force without head moment (free head); moment loading was not considered. In addition, the solution does not take advantage of the exact Winkler solutions presented here for inhomogeneous soils. In this work, the Vlasov and Leontiev (1966) solution is revisited to derive the Winkler spring coefficient for combined force and moment loading, properly accounting for the boundary conditions at the pile head. This solution employs the exact shape functions that are developed in the first part of this paper for inhomogeneous Winkler spring beds.

In summary, in this work, a solution is provided for the response of long piles in inhomogeneous soils under both lateral and moment head loading. This consists of (1) an exact solution to the relevant Winkler problem for the important case of semi-infinite beams, and (2) a systematic approach for obtaining the Winkler spring stiffness in inhomogeneous soil profiles for combined pile head loading. The proposed novel solutions can be used in engineering applications as well as for assessing other related models in geotechnics and engineering mechanics. Comparisons with alternative Winkler calibrations and numerical continuum solutions are provided.

The Problem

The classical Winkler model for a laterally loaded pile is shown in Fig. 1(a). The pile is treated as a cylindrical Euler-Bernoulli elastic solid beam of Young’s modulus

E_{p}

, length

L

, diameter

D

, and second moment of area

I

(

= π D^{4} / 64

). The soil reaction is modeled using distributed Winkler springs of stiffness

k (z)

, which varies with depth

z

and is measured in units of force per length squared (modulus of subgrade reaction).

Fig. 1. (a) Winkler model for a laterally loaded pile; and (b) Winkler stiffness profile under consideration.

In this paper, the following power-law Winkler stiffness profile is considered, shown in Fig. 1(b)

k (z) = k_{D} {(\frac{z + z_{0}}{D + z_{0}})}^{n}

(1a)

where

k_{D}

= Winkler spring stiffness at a depth of one diameter;

n

= power-law exponent; and

z_{0}

= characteristic length corresponding to the elevation above the pile head that

k (z)

(if extrapolated above the ground surface) reaches zero.

z_{0}

is shown graphically in Fig. 1(b) and is given by

\frac{z_{0}}{D} = \frac{k_{0}^{1 / n}}{k_{D}^{1 / n} - k_{0}^{1 / n}}

(1b)

where

k_{0}

= Winkler spring stiffness at the ground surface.

This flexible formulation allows a variety of different soil profiles to be modeled. In particular, existing solutions are available for some special cases, including

n = 0

, which produces a constant (homogeneous) stiffness profile (independent of

z_{0}

, which is undefined for this case);

n = 1 / 2

, which produces a profile with stiffness varying with the square root of depth, representative of some granular deposits (Darendeli 2001); and

n = 1

, which produces a profile with stiffness varying linearly with depth often referred to as Gibson soil (Gibson 1974).

When designing a pile, of most interest to an engineer is usually the response at the pile head. Consider the combined lateral loading at the pile head shown in Fig. 1, consisting of a horizontal shear load

P

and a bending moment

M

. These actions result in a pile head deflection

y_{0}

and rotation (anticlockwise positive)

θ_{0}

, which can be related to the applied loading using either a (

2 \times 2

) stiffness matrix

K

or a corresponding flexibility matrix

F

as shown in Eq. (2)

P = [\begin{matrix} P \\ M \end{matrix}] = K u = [\begin{matrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{matrix}] [\begin{matrix} y_{0} \\ θ_{0} \end{matrix}]

(2a)

u = [\begin{matrix} y_{0} \\ θ_{0} \end{matrix}] = F P = [\begin{matrix} F_{11} & F_{12} \\ F_{21} & F_{22} \end{matrix}] [\begin{matrix} P_{0} \\ M_{0} \end{matrix}]

(2b)

where

K_{11}

,

K_{22}

, and

K_{12}

(

= K_{21}

) = swaying, rocking, and cross-swaying rocking stiffness coefficients, respectively; and

F_{11}

,

F_{22}

, and

F_{12}

(

= F_{21}

) = corresponding flexibilities. These matrices are related by the simple expression

F = K^{- 1}

(2c)

Because both the pile and soil are modeled as linearly elastic, both matrices are symmetric by reciprocity (i.e.,

K_{12} = K_{21}

and

F_{12} = F_{21}

). In addition, the sign of each term depends on the overall sign convention chosen. For the sign convention shown in Fig. 1 (used throughout this paper), the off-diagonal stiffness terms

K_{12}

and

K_{21}

are positive, whereas the off-diagonal flexibility terms

F_{12}

and

F_{21}

are negative.

These matrices are obtained by solving the governing equation of the problem and applying the pertinent boundary conditions for each term. With reference to the problem in Fig. 1, the governing equation for pile deflection

y (z)

can be obtained from the well-known strength-of-materials solution of an Euler-Bernoulli beam under the action of a distributed load,

q (z)

(Chang 1937)

E_{p} I y^{⁗} (z) = q (z) = - k (z) y (z)

(3)

where the prime = differentiation with respect to depth

z

.

For the sign convention in Fig. 1 (positive bending moments inducing negative curvatures) the relations among pile rotation

θ (z)

, bending moment

M (z)

, and shear force

P (z)

are given by

θ (z) = + y^{'} (z)

(4a)

M (z) = - E_{p} I y^{″} (z)

(4b)

P (z) = + E_{p} I y^{‴} (z)

(4c)

Inputting the Winkler stiffness profile under consideration [Eq. (1a)] into Eq. (3) yields the following governing equation

E_{p} I y^{⁗} (z) + k_{D} {(\frac{z + z_{0}}{D + z_{0}})}^{n} y (z) = 0

(5a)

This can be simplified by introducing a Winkler parameter,

λ

(carrying units of

{length}^{- 1}

) such that

y^{⁗} (z) + (n + 4) λ^{n + 4} {[z_{0} + z]}^{n} y (z) = 0

(5b)

λ^{n + 4} = \frac{k_{D}}{(n + 4) E_{p} I {(D + z_{0})}^{n}}

(5c)

Eq. (5c) is a more general form of the commonly employed Winkler parameter for homogeneous soil (Hetenyi 1946; Franklin and Scott 1979; Scott 1981) (also discussed in the following section). This parameter can be interpreted as a characteristic inverse length (“wave number”), which enables significant insights into the nature of the problem and properly normalizes results. Note that the selection of

D

and

k_{D}

as reference values in Eq. (5c) is an arbitrary choice; any other depth and corresponding stiffness would lead to the same

λ

. This is clear from the alternative form of Eq. (5c) provided in Eq. (5d), which uses

k_{0}

and

z_{0}

directly, neither of which depend on the choice of reference depth

λ^{n + 4} = \frac{{(k_{D}^{1 / n} - k_{0}^{1 / n})}^{n}}{(n + 4) E_{p} I D^{n}} = \frac{k_{0}}{(n + 4) E_{p} I z_{0}^{n}}

(5d)

With the introduction of

λ

, only three dimensionless parameters are required to describe any problem modeled as in Fig. 1:

λ L

,

λ z_{0}

, and

n

. This is evident from dimensional analysis, considering the five relevant problem parameters (

k_{0}

,

z_{0}

,

n

,

E I

, and

L

) and the two fundamental dimensions (force and length), yield three governing dimensionless groups (i.e., 5–2). The product

λ L

incorporates both pile-soil relative stiffness and geometric pile slenderness and can be interpreted as a “mechanical pile length” (Di Laora et al. 2012), and the product

λ z_{0}

incorporates a pile-soil stiffness ratio dependent on

k_{0}

,

z_{0}

, and

n

and given by

{(λ z_{0})}^{n + 4} = \frac{k_{0} z_{0}^{4}}{(n + 4) E_{p} I}

(5e)

For zero surface stiffness, this term is zero, and only

λ L

and

n

are required to describe the problem (as demonstrated subsequently in this work). Therefore,

λ z_{0}

can be interpreted as a dimensionless surface stiffness parameter.

Homogeneous Stiffness Profile

A reference application that deserves discussion is that of a homogeneous Winkler bed, i.e., when

k (z) = k

. In this case the governing Eq. (5b) can be simplified to Eq. (6a) with the substitution

n = 0

; the latter equation has the simple solution given by Eq. (6c)

y^{⁗} (z) + 4 λ^{4} y (z) = 0

(6a)

λ^{4} = \frac{k}{4 E_{p} I}

(6b)

y (z) = e^{- λ z} [A_{1} \cos (λ z) + A_{2} \sin (λ z)] + e^{+ λ z} [A_{3} \cos (λ z) + A_{4} \sin (λ z)]

(6c)

where

A_{m}

terms = integration constants to be determined from the boundary conditions at the pile head and base. Eq. (6b) is obtained from Eq. (5c) using the same substitution, which produces the well-known Winkler parameter employed by many authors (e.g., Hetenyi 1946; Scott 1981; Pender 1993; Mylonakis 1995; Mylonakis and Gazetas 1999). The factor of 4 is normally included by convention to simplify the solution in Eq. (6c).

Fig. 2 shows the pile head stiffness and flexibility as a function of dimensionless pile length

λ L

for different base conditions (the plotted results for inhomogeneous stiffness profiles will be discussed subsequently). These are obtained using the pile head boundary conditions from the definition of each stiffness/flexibility term and three different pile base conditions: floating (i.e., zero load and moment), fixed (i.e., zero displacement and rotation), and hinged (i.e., zero displacement and moment). A systematic approach for obtaining these results is included in the Supplemental Materials and has been discussed by Crispin (2022).

Fig. 2. Variation of pile head stiffness with pile length, base conditions, and Winkler stiffness profile: (a) $K_{11}$ ; (b) $K_{12}$ ; (c) $K_{22}$ ; (d) $F_{11}$ ; (e) $F_{12}$ ; and (f) $F_{22}$ .

From Fig. 2, it is evident that beyond a certain length, the boundary conditions at the pile base have negligible effect on the pile head response. It is therefore common to treat piles beyond this length, termed the active pile length

L_{a}

, as long piles (i.e., semi-infinite beams). The first two terms in Eq. (5), associated with the negative exponent

e^{- λ z}

, approach zero as

z

increases, whereas the second two terms, associated with the positive exponent

e^{+ λ z}

, are unbounded. Therefore, integration constants

A_{3}

and

A_{4}

must be set to zero for long piles, significantly simplifying the solution.

Various definitions for

L_{a}

are available such as

λ L_{a} = 2 \sqrt{2}

, suggested by Hetenyi (1946) for homogeneous soils (see also Di Laora and Rovithis 2015; Karatzia and Mylonakis 2017; Mylonakis and Crispin 2021 for a more recent discussion). However, regardless of definition, this category of long piles normally includes most piles installed in practice (with some notable exceptions such as large-diameter offshore monopiles) and is, therefore, the focus of this paper.

The pile head stiffness and/or flexibility terms can be determined from Eq. (5) using the appropriate pairs of boundary conditions at the pile head. For this purpose, it is convenient to introduce four basic response functions describing the pile displacement under certain conditions. In this paper, the following four functions are employed:

•

Y_{1} (z)

, the displaced shape due to unit displacement at the pile head under zero rotation,

•

Y_{2} (z)

, the displaced shape due to unit rotation at the pile head under zero displacement,

•

Y_{3} (z)

, the displaced shape due to unit load at the pile head under zero moment, and

•

Y_{4} (z)

, the displaced shape due to unit moment at the pile head under zero load.

These functions can be obtained using any pair of independent solutions to the governing equation [Eq. (5)] that approach zero at large

z

,

g_{1} (z)

, and

g_{2} (z)

. These are herein referred to as the regular part of the solution. For homogeneous soil, these can be the functions multiplying the integration constants

A_{1}

and

A_{2}

in Eq. (6) [i.e.,

g_{1} (z) = e^{- λ z} \cos (λ z)

,

g_{2} (z) = e^{- λ z} \sin (λ z)

], or any linear combination of them, and result in the expressions in Eq. (7) (Hetenyi 1946; Pender 1993). With the exception of

Y_{1} (z)

, these functions are not dimensionless or unitary. These are plotted in Fig. 3 (the plotted results for inhomogeneous stiffness profiles will be discussed subsequently)

Y_{1} (z) = \frac{g_{1} (z) g_{2}^{'} (0) - g_{2} (z) g_{1}^{'} (0)}{g_{1} (0) g_{2}^{'} (0) - g_{2} (0) g_{1}^{'} (0)} = e^{- λ z} [\cos (λ z) + \sin (λ z)]

(7a)

Y_{2} (z) = \frac{g_{1} (z) g_{2} (0) - g_{2} (z) g_{1} (0)}{g_{1}^{'} (0) g_{2} (0) - g_{2}^{'} (0) g_{1} (0)} = \frac{1}{λ} e^{- λ z} \sin (λ z)

(7b)

Y_{3} (z) = \frac{1}{E I} \frac{g_{1} (z) g_{2}^{″} (0) - g_{2} (z) g_{1}^{″} (0)}{g_{1}^{‴} (0) g_{2}^{″} (0) - g_{2}^{‴} (0) g_{1}^{″} (0)} = \frac{1}{2 E_{p} I λ^{3}} e^{- λ z} \cos (λ z)

(7c)

Y_{4} (z) = \frac{1}{E I} \frac{g_{1} (z) g_{2}^{‴} (0) - g_{2} (z) g_{1}^{‴} (0)}{g_{1}^{″} (0) g_{2}^{‴} (0) - g_{2}^{″} (0) g_{1}^{‴} (0)} = \frac{- 1}{2 E_{p} I λ^{2}} e^{- λ z} [\cos (λ z) - \sin (λ z)]

(7d)

Fig. 3. Basic response functions for different Winkler stiffness profiles: (a) unit head displacement under zero rotation, $Y_{1} (z)$ ; (b) unit head rotation under zero displacement, $Y_{2} (z)$ ; (c) unit head displacement under zero moment and finite load, $Y_{3} (z) / Y_{3} (0)$ ; and (d) unit head displacement under zero load and finite moment, $Y_{4} (z) / Y_{4} (0)$ .

Based on these definitions and the sign convention chosen (Fig. 1), the familiar stiffness and flexibility terms can be obtained directly (e.g., Pender 1993) as follows

K = E_{p} I [\begin{matrix} Y_{1}^{‴} (0) & Y_{2}^{‴} (0) \\ - Y_{1}^{″} (0) & - Y_{2}^{″} (0) \end{matrix}] = E_{p} I [\begin{matrix} 4 λ^{3} & 2 λ^{2} \\ 2 λ^{2} & 2 λ \end{matrix}]

(8a)

F = [\begin{matrix} Y_{3} (0) & Y_{4} (0) \\ Y_{3}^{'} (0) & Y_{4}^{'} (0) \end{matrix}] = \frac{1}{2 E_{p} I} [\begin{matrix} λ^{- 3} & - λ^{- 2} \\ - λ^{- 2} & 2 λ^{- 1} \end{matrix}]

(8b)

Once these have been used to establish the pile head deflection

y_{0}

and rotation

θ_{0}

, the overall displaced shape of the pile can be readily obtained from Eqs. (9a) or (9b)

y (z) = y_{0} Y_{1} (z) + θ_{0} Y_{2} (z)

(9a)

y (z) = P Y_{3} (z) + M Y_{4} (z)

(9b)

which can also be input into Eq. (4) to obtain the bending moment and shear force variation with depth.

Linearly Varying Stiffness Profile

Analysis of piles in other soil profiles requires repeating this process (notably solving the governing equation) for different

k (z)

profiles. To the authors’ knowledge, the homogeneous case is the only commonly encountered soil profile that can be solved with elementary functions. However, solutions are available for soil stiffness varying linearly with depth, one by Hayashi (1921) [repeated by Rifaat (1935) and Hetenyi (1946)], who used an infinite power-series solution, and another by Franklin and Scott (1979), who employed contour integration [also similar to an alternative solution by Hayashi (1921)]; both of these solutions are revisited here.

Consider first a Winkler modulus function increasing linearly from

k (0) = 0

, obtained from Eq. (1) by setting

n = 1

and

z_{0} = 0

. In this case, Eq. (5) can be rewritten

E_{p} {I y}^{⁗} (z) + k_{D} \frac{z}{D} y (z) = 0

(10a)

y^{⁗} (z) + 5 λ^{5} z y (z) = 0

(10b)

λ^{5} = \frac{k_{D}}{5 E_{p} I D}

(10c)

Hayashi (1921) and later Rifaat (1935) and Hetenyi (1946) proposed the following solution to Eq. (10):

y (z) = B_{1} h_{1} (z) + B_{2} h_{2} (z) + B_{3} h_{3} (z) + B_{4} h_{4} (z)

(11a)

where the

B_{m}

terms = coefficients to be determined from the boundary conditions; and the

h_{m} (z)

terms were given in the original work as a power series, which, following more recent solutions (Mylonakis 2004; Fradelos 2016; Crispin 2017; Froio and Rizzi 2017; Parashakis 2020; Rosendo and Albuquerque 2021), can be expressed as generalized hypergeometric functions as follows [modified based on the definition for

λ

in Eq. (10c)]:

h_{1} (z) = \frac{{(λ z)}^{0}}{0!} {{}_{0}F}_{3} (; \frac{4}{5}, \frac{3}{5}, \frac{2}{5}; \frac{- {(λ z)}^{5}}{5^{3}})

(11b)

h_{2} (z) = \frac{{(λ z)}^{1}}{1!} {{}_{0}F}_{3} (; \frac{4}{5}, \frac{3}{5}, \frac{6}{5}; \frac{- {(λ z)}^{5}}{5^{3}})

(11c)

h_{3} (z) = \frac{{(λ z)}^{2}}{2!} {{}_{0}F}_{3} (; \frac{4}{5}, \frac{7}{5}, \frac{6}{5}; \frac{- {(λ z)}^{5}}{5^{3}})

(11d)

h_{4} (z) = \frac{{(λ z)}^{3}}{3!} {{}_{0}F}_{3} (; \frac{8}{5}, \frac{7}{5}, \frac{6}{5}; \frac{- {(λ z)}^{5}}{5^{3}})

(11e)

where

{{}_{0}F}_{3} (; b_{1}, b_{2}, b_{3}; x)

= generalized hypergeometric function of argument

x

and parameters

b_{1}

,

b_{2}

, and

b_{3}

(Olver et al. 2010). The derivatives of these functions are provided in the Supplemental Materials for convenience.

These functions are now (relative to when the solution was first proposed) much easier to compute, with efficient built-in functions available in many modern programming languages. For piles of finite length, the head stiffness and flexibility terms for different pile lengths and base conditions are plotted in Fig. 2, obtained in the same way as the solutions for the homogeneous stiffness profile [see Supplemental Materials and Crispin (2022)]. However, all four functions in Eq. (11) are unbounded at large argument, making them problematic for application to long piles. In principle, an arbitrarily large length, greater than the active length, can be employed to get suitable results; however, this approach is cumbersome and requires enforcing specific boundary conditions at the pile tip.

Franklin and Scott (1979) introduced an alternative solution involving integrating a complex variable

t

along a contour

c_{m}

, as shown in Eq. (12) [modified to use the definition of

λ

in Eq. (10c)]. Hayashi (1921) also proposed a similar solution derived using the Laplace transform.

f_{m} (z) = \int_{c_{m}} \exp (5^{\frac{1}{5}} t λ z + \frac{t^{5}}{5}) d t

(12)

Repeated differentiation of this equation yields the expression for the

j

th derivative

f_{m}^{(j)} (z) = 5^{\frac{j}{5}} λ^{j} \int_{c_{m}} t^{j} \exp (5^{\frac{1}{5}} t λ z + \frac{t^{5}}{5}) d t

(13)

Performing integration by parts on the fourth derivative as follows yields

f_{m}^{⁗} (z) = 5^{\frac{4}{5}} λ^{4} {[\exp (5^{\frac{1}{5}} t λ z + \frac{t^{5}}{5})]}_{c_{m}} - 5 λ^{5} z \int_{c_{m}} \exp (5^{\frac{1}{5}} t λ z + \frac{t^{5}}{5}) d t

(14)

This shows that if

c_{m}

is chosen such that the integrated part (first term on the right-hand side) vanishes,

y (z) = f_{m} (z)

is a solution to Eq. (10). This naturally occurs when both limits tend to infinity with the same complex argument as a root of

{(- 1)}^{1 / 5}

. The five contours with this property chosen by Franklin and Scott (1979) are shown in Fig. 4.

Fig. 4. Contour integration paths for Eq. (12). (Reprinted from Franklin and Scott 1979, © ASCE.)

This results in five functions, corresponding to computing

f_{m} (z)

from Eq. (12) with each contour

c_{m}

. These are easiest to calculate by taking the contours to follow exactly the limits of their defined region (dashed lines in Fig. 4). Each of these can be separated into two arms, the first approaching zero from one limit, and the second diverging from zero toward the second limit. Because the contours are drawn to use the opposite signs along each shared arm, the functions sum to zero at any

z

and are therefore dependent. Furthermore, from the symmetry shown in Fig. 4, the following relations between the functions can be derived in terms of each function’s real,

R [f_{m} (z)]

, and imaginary,

I [f_{m} (z)]

, components:

•

R [f_{0} (z)] = - R [f_{4} (z)]

,

•

I [f_{0} (z)] = + I [f_{4} (z)]

,

•

R [f_{1} (z)] = - R [f_{3} (z)]

, and

•

I [f_{1} (z)] = + I [f_{3} (z)]

.

These properties have been discussed further by Franklin and Scott (1979).

Franklin and Scott (1979) proposed the following general solution to Eq. (10) and showed that each of the four solutions are independent

y (z) = C_{1} R [f_{0} (z)] + C_{2} I [f_{0} (z)] + C_{3} R [f_{1} (z)] + C_{4} I [f_{1} (z)]

(15)

where the

C_{m}

terms = coefficients to be determined from the boundary conditions.

Taking the contours to follow exactly the limits of their defined region, the

j

th derivative of

f_{0} (z)

and

f_{1} (z)

are given by [modified from Franklin and Scott (1979)]

f_{0}^{(j)} (z) = - e^{π i (j + 1)} 5^{\frac{j}{5}} λ^{j} \int_{0}^{\infty} t^{j} \exp (5^{\frac{1}{5}} e^{π i} t λ z - \frac{t^{5}}{5}) d t + e^{\frac{3 π i}{5} (j + 1)} 5^{\frac{j}{5}} λ^{j} \int_{0}^{\infty} t^{j} \exp (5^{\frac{1}{5}} e^{\frac{3 π i}{5}} t λ z - \frac{t^{5}}{5}) d t

(16a)

f_{1}^{(j)} (z) = - e^{\frac{3 π i}{5} (j + 1)} 5^{\frac{j}{5}} λ^{j} \int_{0}^{\infty} t^{j} \exp (5^{\frac{1}{5}} e^{\frac{3 π i}{5}} t λ z - \frac{t^{5}}{5}) d t + e^{\frac{π i}{5} (j + 1)} 5^{\frac{j}{5}} λ^{j} \int_{0}^{\infty} t^{j} \exp (5^{\frac{1}{5}} e^{\frac{π i}{5}} t λ z - \frac{t^{5}}{5}) d t

(16b)

Scott (1981) provided tabulated results for the real and imaginary components of these expressions obtained with a numerical integration technique (those results employed the dimensionless independent variable

x = 5^{1 / 5} λ z

). Liang et al. (2014) proposed an approximate analytical approach for evaluating these expressions using a power series and the WKB approximation (Bender and Orszag 1978).

From the sign of the terms multiplying

z

, it is evident that

f_{0} (z)

approach zero as

z

increases, and

f_{1} (z)

is unbounded. In addition,

f_{0} (z)

and

f_{1} (z)

are uniquely defined as solutions to Eq. (10) and by four criteria on the function and its derivatives. Therefore, it is possible to obtain these functions as a linear combination of the generalized hypergeometric functions in Eq. (11) [i.e., select

B_{m}

terms in Eq. (11a) such that

y (z) = f_{0} (z)

or

f_{1} (z)

]. Both

f_{0} (z)

and

f_{1} (z)

can be evaluated analytically in terms of the gamma function

Γ ()

at

z = 0

(Franklin and Scott 1979)

f_{0}^{(j)} (0) = 5^{\frac{2 j - 4}{5}} λ^{j} Γ (\frac{j + 1}{5}) [- e^{π i (j + 1)} + e^{\frac{3 π i}{5} (j + 1)}]

(17a)

f_{1}^{(j)} (0) = 5^{\frac{2 j - 4}{5}} λ^{j} Γ (\frac{j + 1}{5}) [- e^{\frac{3 π i}{5} (j + 1)} + e^{\frac{π i}{5} (j + 1)}]

(17b)

The power-series solution in Eq. (11) has the following convenient property at

z = 0

:

y^{(j)} (0) = λ^{j} B_{j + 1}

(18)

where

y^{(j)} (z)

=

j

th derivative of

y (z)

with respect to

z

, as previously.

Therefore, by setting Eq. (18) equal to Eq. (17) and letting

j

vary from 0 to 3, each of the

B_{m}

terms and, consequently, analytical expressions for Eq. (16) can be obtained

f_{0} (z) = 5^{- \frac{4}{5}} [1 + {(- 1)}^{\frac{3}{5}}] Γ (\frac{1}{5}) h_{1} (z) - 5^{- \frac{2}{5}} [1 + {(- 1)}^{\frac{1}{5}}] Γ (\frac{2}{5}) h_{2} (z) + [1 + {(- 1)}^{- \frac{1}{5}}] Γ (\frac{3}{5}) h_{3} (z) - 5^{\frac{2}{5}} [1 + {(- 1)}^{- \frac{3}{5}}] Γ (\frac{4}{5}) h_{4} (z)

(19a)

f_{1} (z) = 5^{- \frac{4}{5}} [{(- 1)}^{\frac{3}{5}} - {(- 1)}^{\frac{1}{5}}] Γ (\frac{1}{5}) h_{1} (z) + 5^{- \frac{2}{5}} [{(- 1)}^{\frac{6}{5}} - {(- 1)}^{\frac{2}{5}}] Γ (\frac{2}{5}) h_{2} (z) + [{(- 1)}^{\frac{9}{5}} - {(- 1)}^{\frac{3}{5}}] Γ (\frac{3}{5}) h_{3} (z) + 5^{\frac{2}{5}} [{(- 1)}^{\frac{2}{5}} - {(- 1)}^{\frac{4}{5}}] Γ (\frac{4}{5}) h_{4} (z)

(19b)

This combined solution can then be input into Eq. (15) and the response functions (as defined for the homogeneous case) obtained by enforcing the corresponding boundary conditions, i.e., employing the relationships shown in Eq. (7) with

g_{1} (z) = R [f_{0} (z)]

and

g_{2} (z) = I [f_{0} (z)]

, yields the four basic response functions for

n = 1

as follows

Y_{1} (z) = h_{1} (z) - 5^{\frac{4}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{3}{5})}{Γ (\frac{1}{5})} h_{3} (z) + 5^{\frac{6}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{4}{5})}{Γ (\frac{1}{5})} h_{4} (z)

(20a)

λ Y_{2} (z) = h_{2} (z) - 5^{\frac{2}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{3}{5})}{Γ (\frac{2}{5})} h_{3} (z) + 5^{\frac{4}{5}} \frac{Γ (\frac{4}{5})}{Γ (\frac{2}{5})} h_{4} (z)

(20b)

E_{p} I λ^{3} Y_{3} (z) = 5^{- \frac{6}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{1}{5})}{Γ (\frac{4}{5})} h_{1} (z) - 5^{\frac{- 4}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{2}{5})}{Γ (\frac{4}{5})} h_{2} (z) + h_{4} (z)

(20c)

E_{p} I λ^{2} Y_{4} (z) = - 5^{- \frac{4}{5}} \frac{Γ (\frac{1}{5})}{Γ (\frac{3}{5})} h_{1} (z) + 5^{- \frac{2}{5}} \frac{1 + \sqrt{5}}{2} \frac{Γ (\frac{2}{5})}{Γ (\frac{3}{5})} h_{2} (z) - h_{3} (z)

(20d)

where

h_{1} (z)

,

h_{2} (z)

,

h_{3} (z)

, and

h_{4} (z)

= hypergeometric functions in Eq. (11). It is worth stressing that although all four functions

h_{m} (z)

diverge with increasing

z

, the basic response functions in Eq. (20) tend to zero as

z

approaches infinity. These results complement and extend the solution of Franklin and Scott (1979).

The preceding solutions are plotted in Fig. 3, and tabulated results for the functions and their derivatives are included in the Supplemental Materials. The stiffness and flexibility matrices can be obtained directly from the derivatives of these functions with the relationships in Eq. (8), leading to an explicit solution for the stiffness and flexibility matrices

K = E_{p} I \frac{1 + \sqrt{5}}{2} [\begin{matrix} 5^{\frac{6}{5}} \frac{Γ (\frac{4}{5})}{Γ (\frac{1}{5})} λ^{3} & 5^{\frac{4}{5}} \frac{Γ (\frac{3}{5})}{Γ (\frac{1}{5})} λ^{2} \\ 5^{\frac{4}{5}} \frac{Γ (\frac{3}{5})}{Γ (\frac{1}{5})} λ^{2} & 5^{\frac{2}{5}} \frac{Γ (\frac{3}{5})}{Γ (\frac{2}{5})} λ \end{matrix}] ≅ E_{p} I [\begin{matrix} 2.831 λ^{3} & 1.902 λ^{2} \\ 1.902 λ^{2} & 2.068 λ \end{matrix}]

(21a)

F = \frac{1}{E_{p} I} \frac{1 + \sqrt{5}}{2} [\begin{matrix} \frac{5^{- \frac{6}{5}}}{λ^{3}} \frac{Γ (\frac{1}{5})}{Γ (\frac{4}{5})} & - \frac{5^{- \frac{4}{5}}}{λ^{2}} \frac{Γ (\frac{2}{5})}{Γ (\frac{4}{5})} \\ - \frac{5^{- \frac{4}{5}}}{λ^{2}} \frac{Γ (\frac{2}{5})}{Γ (\frac{4}{5})} & \frac{5^{- \frac{2}{5}}}{λ} \frac{Γ (\frac{2}{5})}{Γ (\frac{3}{5})} \end{matrix}] ≅ \frac{1}{E_{p} I} [\begin{matrix} 0.925 λ^{- 3} & - 0.851 λ^{- 2} \\ - 0.851 λ^{- 2} & 1.266 λ^{- 1} \end{matrix}]

(21b)

In which

λ^{5} = k_{D} / (5 E_{p} I D)

[Eq. (10c)]. Each of the terms in Eq. (21) has the same dependence on

λ

as in Eq. (8) (for homogeneous soil). Yet,

λ

in this case is proportional to

{(k / E_{p})}^{1 / 5}

instead of an exponent of 1/4 for homogeneous soil [Eq. (10c)]. Therefore, for linearly varying spring stiffness with depth, the pile response has a higher dependence on

E_{p}

and a lower dependence on

k

(and therefore

E_{s}

) than for a homogeneous soil. This trend naturally gets more pronounced with increasing

n

, as discussed subsequently in this paper.

For a nonzero surface stiffness

k_{0}

, Eq. (10) must be modified as follows

y^{⁗} (z) + 5 λ^{5} [z_{0} + z] y (z) = 0

(22a)

\frac{z_{0}}{D} = \frac{k_{0}}{k_{D} - k_{0}}

(22b)

λ^{5} = \frac{k_{D}}{5 E_{p} I (z_{0} + D)}

(22c)

where

z_{0}

, as introduced in Eq. (1b), is the elevation above ground level where

k (z)

(when extrapolated) would reach zero [Fig. 1(b)].

This governing equation can be transformed to Eq. (10) with the simple substitution

z \to z + z_{0}

. Therefore, the solutions in Eqs. (10) and (15) are still valid with this substitution employed. However, the response functions in Eq. (20) [and consequently the matrices in Eq. (21)] must be recalculated from the relationships in Eq. (6) with

g_{1} (z)

and

g_{2} (z)

taken as

R [f_{0} (z + z_{0})]

and

I [f_{0} (z + z_{0})]

, respectively [or any pair of response functions from Eq. (20) with the same substitution]. The resulting stiffness and flexibility terms are plotted in Figs. 5(b and c), respectively. Explicit closed-form expressions for these terms are possible; however, these are long and impractical and therefore not provided here. Naturally, for Winkler spring stiffness decreasing with depth, solutions for long piles are not available due to the stiffness becoming negative beyond a certain depth.

Fig. 5. Variation of (a) stiffness; and (b) flexibility terms with surface stiffness for a linear Winkler stiffness profile [Eq. (1) with $n = 1$ ; $λ$ is also provided in Eq. (22c)].

Power-Law Stiffness Profile

For the power-law Winkler profile given in Eq. (1), a solution to Eq. (5) can be obtained using Frobenius series (Trachanas 1989), i.e., the same approach employed by Hayashi (1921), Rifaat (1935), and Hetenyi (1946). This approach has previously been applied to this problem by Mylonakis (2004), Fradelos (2016), Crispin (2017), Parashakis (2020), and Rosendo and Albuquerque (2021) and is extended here.

Consider the simplest case with zero surface stiffness (

z_{0} = 0

). Firstly, a simple monomial trial solution

y (z) = z^{s}

can be substituted into Eq. (5b) to provide the indicial equation

s (s - 1) (s - 2) (s - 3) z^{s - 4} + (n + 4) λ^{n + 4} z^{s + n} = 0

(23)

From this result, it can be observed that for

n

greater than

- 4

, a solution to Eq. (5b) in terms of Frobenius series of step (

n + 4

) is possible when

s = 0

, 1, 2, or 3 as follows

y (z) = z^{s} \sum_{m = 0}^{\infty} α_{m} z^{m (n + 4)}

(24)

where

α_{m}

= coefficient of the

m

th term. Inputting this solution into Eq. (5b) yields the relationship between sequential

α_{m}

coefficients as follows

α_{m + 1} = \frac{- (n + 4) λ^{n + 4}}{[s + (m + 1) (n + 4)] [s + (m + 1) (n + 4) - 1] [s + (m + 1) (n + 4) - 2] [s + (m + 1) (n + 4) - 3]} α_{m}

(25)

which can be used to derive the following four solutions to Eq. (5b) that can be used with Eq. (11a)

h_{1} (z) = \frac{{(λ z)}^{0}}{0!} {{}_{0}F}_{3} (; 1 - ν, 1 - 2 ν, 1 - 3 ν; - ν^{3} {(λ z)}^{\frac{1}{ν}})

(26a)

h_{2} (z) = \frac{{(λ z)}^{1}}{1!} {{}_{0}F}_{3} (; 1 - ν, 1 - 2 ν, 1 + ν; - ν^{3} {(λ z)}^{\frac{1}{ν}})

(26b)

h_{3} (z) = \frac{{(λ z)}^{2}}{2!} {{}_{0}F}_{3} (; 1 - ν, 1 + 2 ν, 1 + ν; - ν^{3} {(λ z)}^{\frac{1}{ν}})

(26c)

h_{4} (z) = \frac{{(λ z)}^{3}}{3!} {{}_{0}F}_{3} (; 1 + 3 ν, 1 + 2 ν, 1 + ν; - ν^{3} {(λ z)}^{\frac{1}{ν}})

(26d)

where

ν = \frac{1}{n + 4}

(26e)

The derivatives of these functions, obtained using standard differentiation rules for hypergeometric functions, are provided in the Supplemental Materials for convenience. When

n = 1

, these expressions are the same as those in Eq. (11), and when

n = 0

, they are solutions to Eq. (6).

For piles of finite length, the head stiffness and flexibility terms for different pile lengths and base conditions are plotted in Fig. 2 for

n = 1 / 2

and

n = 1

. These were obtained in the same way as the solutions for the homogeneous stiffness profiles [Supplemental Materials and Crispin (2022)]. The following important trends are worthy of note:

•

The

K_{11}

and

K_{22}

terms are strictly increasing functions of pile length for floating piles, and decreasing functions of depth for piles fixed or hinged at the tip.

•

A floating pile can never have higher stiffness

K_{11}

and

K_{22}

than a pile fixed or hinged at the tip. This, however, is possible for the cross term

K_{12}

.

•

The opposite naturally holds for the flexibility terms

F_{11}

,

F_{22}

, and

F_{12}

.

•

Beyond a certain pile length, the changes on pile head stiffness become asymptotic, and solutions for different conditions at the tip converge.

Like in the case of linearly varying stiffness, the four hypergeometric functions in Eq. (26) diverge with increasing

z

. Accordingly, application of this general solution to long piles requires extracting the regular part of the solution, i.e., obtaining at least two combinations of

B_{m}

terms such that Eq. (11a) approaches zero for large argument. This can be done by inspection based on the asymptotic properties of the hypergeometric functions to get the regular part of the solution, resulting in

g_{1} (z) = \frac{h_{1} (z)}{Γ (ν) Γ (2 ν) Γ (3 ν)} - \frac{\cos (π ν) ν^{1 - 3 ν} h_{2} (z)}{Γ (ν) Γ (2 ν) Γ (1 - ν)} + \frac{\cos (2 π ν) ν^{2 - 6 ν} h_{3} (z)}{Γ (ν) Γ (1 - 2 ν) Γ (1 - ν)} - \frac{\cos (3 π ν) ν^{3 - 9 ν} h_{4} (z)}{Γ (1 - 3 ν) Γ (1 - 2 ν) Γ (1 - ν)}

(27a)

g_{2} (z) = \frac{\sin (π ν) ν^{1 - 3 ν} h_{2} (z)}{Γ (ν) Γ (2 ν) Γ (1 - ν)} - \frac{\sin (2 π ν) ν^{2 - 6 ν} h_{3} (z)}{Γ (ν) Γ (1 - 2 ν) Γ (1 - ν)} + \frac{\sin (3 π ν) ν^{3 - 9 ν} h_{4} (z)}{Γ (1 - 3 ν) Γ (1 - 2 ν) Γ (1 - ν)}

(27b)

The basic response functions are obtained, as before, by substituting Eq. (27) into the relationships in Eq. (7)

Y_{1} (z) = h_{1} (z) - ν^{6 ν - 2} \frac{Γ (1 - 2 ν) Γ (1 - ν)}{Γ (2 ν) Γ (3 ν)} h_{3} (z) + ν^{9 ν - 3} \frac{Γ (1 - 3 ν) Γ^{2} (1 - ν)}{Γ^{2} (2 ν) Γ (3 ν)} h_{4} (z)

(28a)

λ Y_{2} (z) = h_{2} (z) - ν^{3 ν - 1} \frac{Γ (1 - ν) Γ (ν)}{Γ^{2} (2 ν)} h_{3} (z) + ν^{6 ν - 2} \frac{Γ (1 - 2 ν) Γ (1 - ν)}{Γ (2 ν) Γ (3 ν)} h_{4} (z)

(28b)

E_{p} I λ^{3} Y_{3} (z) = ν^{3 - 9 ν} \frac{Γ^{2} (ν) Γ (3 ν)}{Γ (1 - 3 ν) Γ^{2} (1 - 2 ν)} h_{1} (z) - ν^{2 - 6 ν} \frac{Γ (ν) Γ (2 ν)}{Γ (1 - 3 ν) Γ (1 - 2 ν)} h_{2} (z) + h_{4} (z)

(28c)

E_{p} I λ^{2} Y_{4} (z) = - ν^{2 - 6 ν} \frac{Γ (ν) Γ (2 ν)}{Γ (1 - 3 ν) Γ (1 - 2 ν)} h_{1} (z) + ν^{1 - 3 ν} \frac{Γ (1 - ν) Γ (ν)}{Γ^{2} (2 ν)} h_{2} (z) - h_{3} (z)

(28d)

Tabulated results for these functions and their derivatives are included in the Supplemental Materials. The stiffness and flexibility matrices can be obtained in explicit form directly from these functions by using the relationships shown in Eq. (8) as follows

K = E_{p} I [\begin{matrix} ν^{9 ν - 3} \frac{Γ (1 - 3 ν) Γ^{2} (1 - ν)}{Γ^{2} (2 ν) Γ (3 ν)} λ^{3} & ν^{6 ν - 2} \frac{Γ (1 - 2 ν) Γ (1 - ν)}{Γ (2 ν) Γ (3 ν)} λ^{2} \\ ν^{6 ν - 2} \frac{Γ (1 - 2 ν) Γ (1 - ν)}{Γ (2 ν) Γ (3 ν)} λ^{2} & ν^{3 ν - 1} \frac{Γ (ν) Γ (1 - ν)}{Γ^{2} (2 ν)} λ \end{matrix}]

(29a)

F = \frac{1}{E_{p} I} [\begin{matrix} \frac{ν^{3 - 9 ν}}{λ^{3}} \frac{Γ^{2} (ν) Γ (3 ν)}{Γ (1 - 3 ν) Γ^{2} (1 - 2 ν)} & - \frac{ν^{2 - 6 ν}}{λ^{2}} \frac{Γ (ν) Γ (2 ν)}{Γ (1 - 3 ν) Γ (1 - 2 ν)} \\ - \frac{ν^{2 - 6 ν}}{λ^{2}} \frac{Γ (ν) Γ (2 ν)}{Γ (1 - 3 ν) Γ (1 - 2 ν)} & \frac{ν^{1 - 3 ν}}{λ} \frac{Γ (ν) Γ (1 - ν)}{Γ^{2} (1 - 2 ν)} \end{matrix}]

(29b)

where

λ^{n + 4} = \frac{k_{D}}{[(n + 4) E_{p} I D^{n}]}

[given in Eqs. (5c) or (5d) with

k_{0} = 0

];

ν = 1 / (n + 4)

[given in Eq. (26e)]. The resulting head stiffness terms for different

n

values between 0 and 2, the range most likely to be encountered in practice, are provided in Table 1 with a precision of four significant digits. The preceding explicit solutions can be used instead of the (limited) tabulated numerical results in the literature to compute the pile head stiffness and flexibility for any power-law spring bed.

Table 1. Pile head stiffness terms for a power-law Winkler spring stiffness profile

k (z) = k_{D} {(z / D)}^{n}

with zero surface stiffness

$n$	$\frac{K_{11}}{E_{p} I λ^{3}}$	$\frac{K_{12}}{E_{p} I λ^{2}}$	$\frac{K_{22}}{E_{p} I λ}$	$\frac{1}{F_{11} E_{p} I λ}$	$\frac{- 1}{F_{12} E_{p} I λ}$	$\frac{1}{F_{12} E_{p} I λ}$
0	4.000	2.000	2.000	2.000	2.000	1.000
1/4	3.491	1.953	2.015	1.598	1.649	0.922
1/2	3.175	1.924	2.032	1.353	1.428	0.866
3/4	2.969	1.908	2.049	1.192	1.280	0.823
1	2.831	1.902	2.068	1.081	1.176	0.790
3/2	2.674	1.909	2.106	0.945	1.042	0.744
2	2.609	1.931	2.145	0.870	0.966	0.715

Nonzero surface stiffness can be incorporated, as previously, with the simple substitution

z \to z + z_{0}

. The response functions must be recalculated from the relationships in Eq. (7) with

g_{1} (z)

and

g_{2} (z)

taken as

g_{1} (z + z_{0})

and

g_{2} (z + z_{0})

, respectively, from Eq. (27) [or any pair of response functions from Eq. (28) with the same substitution]. The resulting stiffness and flexibility terms are obtained, again, from the relationships in Eq. (8). For the important case of a square-root Winkler stiffness profile, where

n = 1 / 2

, possessing finite stiffness at the surface, these are plotted in Fig. 6 alongside the result from the numerical method employed by Franklin and Scott (1979) [tabulated by Scott (1981)]. Evidently, the results of the analytical and numerical solutions are practically indistinguishable.

Fig. 6. Variation of (a) stiffness; and (b) flexibility terms with surface stiffness for a square-root Winkler stiffness profile [Eq. (1) with $n = 1 / 2$ ].

As the dimensionless surface stiffness

k_{0}

and the associated dimensionless parameter

λ z_{0}

increase, it is expected that the pile head stiffness and flexibility terms approach those for homogeneous soil because the Winkler stiffness will vary less over the active pile length. However, this is not apparent in Figs. 5 and 6 because, despite the soil profile approaching being homogeneous for increasing

z_{0}

, the

λ

value from Eq. (5c) used to normalize results does not reduce to the homogeneous value in Eq. (6b). Instead, it remains dependent on the

n

value initially employed, which no longer affects the soil profile. An alternative normalization approach must therefore be employed to observe the natural trend. One approach, shown in Fig. 7, is to use

λ_{0}

to normalize the stiffness terms in the ordinates of the graphs. This stiffness is defined in Eq. (30) and is analogous to

λ

for the constant portion of the Winkler spring stiffness. Evidently, the pile stiffness and flexibility terms now approach those for homogenous soil with increasing surface stiffness

λ_{0}^{4} = \frac{k_{0}}{4 E_{p} I} = λ^{4} \frac{n + 4}{4} {(λ z_{0})}^{n}

(30)

Fig. 7. Comparison of (a) stiffness; and (b) flexibility terms for power-law Winkler stiffness profiles [Eq. (1)] with a homogeneous profile ( $n = 0$ ).

Assessment of Winkler Spring Stiffness

Application of the Winkler model to the pile problem for any soil inhomogeneity function requires calibration of the Winkler spring stiffness,

k (z)

. Therefore, the soil continuum must be considered in full, adding significant complexity to the problem. Existing solutions (discussed further in “Introduction” section) often employed boundary or finite-element analysis; however, in this work, an approximate analytical solution is employed. Considering the equilibrium of a soil element in an inhomogeneous half-space yields the following governing expression in cylindrical coordinates (Veletsos and Younan 1994)

η_{s}^{2} G_{s} (z) \frac{\partial}{\partial r} [\frac{1}{r} (\frac{\partial (r u)}{\partial r} + \frac{\partial v}{\partial θ})] - G_{s} (z) \frac{1}{r^{2}} \frac{\partial}{\partial θ} [\frac{\partial (r v)}{\partial r} - \frac{\partial u}{\partial θ}] + \frac{\partial}{\partial z} [G_{s} (z) \frac{\partial u}{\partial z}] = 0

(31a)

G_{s} (z) \frac{\partial}{\partial r} [\frac{1}{r} (\frac{\partial (r u)}{\partial r} - \frac{\partial v}{\partial θ})] + η_{s}^{2} G_{s} (z) \frac{1}{r^{2}} \frac{\partial}{\partial θ} [\frac{\partial (r v)}{\partial r} + \frac{\partial u}{\partial θ}] + \frac{\partial}{\partial z} [G_{s} (z) \frac{\partial u}{\partial z}] = 0

(31b)

where

r

= distance from the center of the pile;

θ

= aperture angle from the direction of loading;

u

and

v

= radial and tangential displacements of the soil, respectively;

G_{s} (z)

= soil shear modulus profile with depth; and

η_{s}

= compressibility constant, which is a function of the soil Poisson’s ratio,

ν_{s}

. The compressibility constant

η_{s}

depends on the choice of restraint employed to simplify the three-dimensional problem; in this case Mylonakis (2001) proposed assuming the additional vertical normal stress arising from pile movement,

σ_{z}

, to be small, to get

η_{s}^{2} = (2 - ν_{s}) / (1 - ν_{s})

.

Evidently, obtaining analytical solutions for the three-dimensional problem described by the partial differential equation in Eq. (31) is far more complex than solving the governing equation of the simplified Winkler model in Eq. (5). Nevertheless, a rational simplification to Eq. (31), using information obtained from the solution to the Winkler model, can be employed to reduce the problem to two dimensions. To this end, the horizontal soil displacement vector is first expressed in separable form using the planar radial and tangential soil displacement components

u (r, θ)

and

v (r, θ)

(defined at an arbitrary depth), times a dimensionless unitary shape function of depth,

Y (z) / Y (0)

, which is derived from the (nonunitary) basic response function

Y (z)

in Eq. (28), upon division by its value at the pile head,

Y (0)

as follows

u (r, z, θ) = u (r, θ) Y (z) / Y (0)

(32a)

v (r, z, θ) = v (r, θ) Y (z) / Y (0)

(32b)

Vlasov and Leontiev (1966) employed a similar approximation, but instead, they used a shape function in the direction perpendicular to the foundation-soil interface. Previous authors have used this approach for laterally loaded piles with a shape function for the attenuation in the radial coordinate,

r

(Vallabhan and Das 1988; Guo and Lee 2001; Basu and Salgado 2008; Basu et al. 2009). However, this approach results in a higher-order Winkler model of the Pasternak type with an extra term proportional to the second derivative of pile displacement. This can be interpreted either as an additional bed of distributed rotational springs, a membrane under tension connecting the Winkler springs at their base, or a shear beam connecting the base of the springs. This term adds significant complexity to solving the Winkler problem; to the best of the authors’ knowledge, analytical solutions to the Pasternak model are only available for homogeneous soil.

By inputting Eq. (32) into Eq. (31) and integrating over the vertical coordinate, a reduced two-dimensional pair of equilibrium equations is obtained (Mylonakis 2001) as follows

η_{s}^{2} \frac{\partial}{\partial r} [\frac{1}{r} (\frac{\partial (r u)}{\partial r} + \frac{\partial v}{\partial θ})] - \frac{1}{r^{2}} \frac{\partial}{\partial θ} [\frac{\partial (r v)}{\partial r} - \frac{\partial u}{\partial θ}] - b^{2} u = 0

(33a)

\frac{\partial}{\partial r} [\frac{1}{r} (\frac{\partial (r u)}{\partial r} - \frac{\partial v}{\partial θ})] + η_{s}^{2} \frac{1}{r^{2}} \frac{\partial}{\partial θ} [\frac{\partial (r v)}{\partial r} + \frac{\partial u}{\partial θ}] - b^{2} v = 0

(33b)

where

b

is a constant (dimensions of

{length}^{- 1}

) encompassing the properties and response of the soil and pile medium with depth. For the common case of zero tractions at the soil surface and zero displacements at an infinite depth,

b

is obtained from (modified from Mylonakis 2001)

b^{2} = λ^{2} \int_{0}^{\infty} G_{s} (λ z) {[\frac{d Y (λ z)}{d (λ z)}]}^{2} d (λ z) / \int_{0}^{\infty} G_{s} (λ z) {[Y (λ z)]}^{2} d (λ z)

(33c)

Remarkably, Eq. (33) perfectly matches the well-known plane-strain problem of Baranov and Novak (Novak 1974) governing the dynamic response of a horizontal soil slice under harmonic loading, with the

b

term retaining information from the suppressed third dimension being equivalent to a wave number normalizing the radial coordinate

r

or the pile radius

D / 2

. The Baranov-Novak solution for the Winkler spring stiffness is (Novak 1974)

k / E_{s} = \frac{π s^{2}}{2 (1 + ν_{s})} \frac{4 K_{1} (q) K_{1} (s) + s K_{1} (q) K_{0} (s) + q K_{0} (q) K_{1} (s)}{q s K_{0} (q) K_{0} (s) + s K_{1} (q) K_{0} (s) + q K_{0} (q) K_{1} (s)}

(34)

where

s

=

b D / 2

(instead of the dimensionless frequency

a_{o}

in the original model);

q = s / η_{s}

; and

K_{μ} ()

is the modified Bessel function of the second kind and order

μ

.

Because

b

is dependent on

Y (z)

which, in turn, depends on the Winkler parameter

λ

, the modulus

k

is present on both sides of Eq. (34) so, in principle, it must be determined iteratively. Karatzia and Mylonakis (2017) suggested that a good approximation can be obtained using only one iteration with an initial value of

k \approx E_{s}

. In addition, they provided a simplified version of Eq. (34) using the asymptotic forms of the Bessel functions for small argument pertaining to static conditions as follows

k / E_{s} = \frac{2 π η_{s}^{2} / (1 + ν_{s})}{\ln (η_{s}) - (1 + η_{s}^{2}) \ln (χ b D)}

(35)

where

χ = e^{γ} / 4 \approx 0.445

;

η_{s}^{2} = (2 - ν_{s}) / (1 - ν_{s})

; and

γ

(

\approx 0.557

) = Euler’s constant.

The shape function used to calculate

b

,

Y (z)

, depends on both the soil stiffness profile and the head loading conditions. For a fixed-head pile in homogeneous soil, using

Y (z) = Y_{1} (z)

from Eq. (7a) results in

b = λ \sqrt{2 / 3}

; similarly, for a free-head pile (under load only), using

Y (z) = Y_{3} (z)

from Eq. (7a) results in

b = λ \sqrt{2}

(Karatzia and Mylonakis 2017). For a general combination of head load

P

and moment

M, (z)

can be taken from the displaced shape given in Eq. (9b). For homogeneous soil, this can be resolved analytically resulting in Eq. (36), which is plotted in Fig. 8(a)

{(b / λ)}^{2} = \frac{6 P^{2} - 16 P λ M + 12 {(λ M)}^{2}}{3 P^{2} - 4 P λ M + 2 {(λ M)}^{2}}

(36)

Fig. 8. Variation of (a) $b$ [Eq. (36)] and normalized Winkler spring stiffness $k / k_{F}$ [Eq. (34)]; and (b) combined pile head stiffness with the applied load-moment ratio for a pile in homogeneous soil $n = 0$ . Note the nonlinear horizontal axis for $| 2 λ M / P | > 2$ ( $E_{p} / E_{s} = 1,000$ , and $ν_{s} = 0.4$ ).

The resulting Winkler spring stiffness

k

, obtained using Eq. (36), is also plotted in Fig. 8(a), normalized by the corresponding value for fixed-head piles,

k_{F}

. Evidently, despite it being common to only use the fixed-head value (corresponding to an applied load to moment ratio of

2 λ

) (Fig. 8), there is significant variation in Winkler spring stiffness depending on the load to moment ratio. Fig. 8(b) shows the difference in a combined stiffness term

(P - λ M) / y_{0}

when using the common assumption of the fixed-head Winkler spring stiffness compared with properly accounting for the load-moment ratio via Eq. (36). This combined stiffness was underestimated by around 25% for free-head piles with no head moment using the traditional approach and by as much as 60% for other load-moment ratios.

It is inconvenient to calculate

k

for each specific load-moment ratio to be considered, particularly for inhomogeneous soils, where analytical solutions to Eq. (33c) are hard to obtain. Instead, it is possible to take advantage of the property that only three terms are required to uniquely define both the stiffness and flexibility matrices in Eq. (2). Indeed, by establishing the spring moduli for a fixed-head pile (lateral force only, zero rotation),

k_{F}

, an applied load only (free-head and zero moment),

k_{P}

, and an applied moment only (free-head and zero lateral force),

k_{M}

, the corresponding stiffness terms

K_{11}

,

F_{11}

, and

F_{12}

can be obtained, and the remaining terms calculated as follows

F_{22} = K_{11} F_{12}^{2} / (F_{11} K_{11} - 1)

(37a)

K_{12} = K_{21} = (F_{11} K_{11} - 1) / (- F_{12})

(37b)

K_{22} = \frac{F_{11}}{F_{12}^{2}} (F_{11} K_{11} - 1)

(37c)

Here,

k_{P}

,

k_{M}

, and

k_{F}

are determined using their corresponding

b

terms, namely

b_{F}

,

b_{P}

, and

b_{M}

, which are, in turn, calculated from Eq. (33c) using functions

Y_{1} (z)

,

Y_{3} (z)

, and

Y_{4} (z)

, respectively, as the shape function

Y (z)

. The

b

values from numerical evaluation of these integrals are given in Table 2 for homogeneous, square-root, and linear soil profiles, the latter two possessing zero surface stiffness.

Table 2. Values of parameter

b

[Eq. (33)] for different soil stiffness profiles and boundary conditions at the pile head:

b_{F}

(fixed-head),

b_{P}

(free-head, force only), and

b_{M}

(free-head, moment only)

Stiffness profile	$n$	$b_{F} / λ$	$b_{P} / λ$	$b_{M} / λ$
$E_{s} (z) = E_{s}$	0	0.816 $(\sqrt{2 / 3})$	1.414 $(\sqrt{2})$	2.449 $(\sqrt{6})$
$E_{s} (z) = E_{s D} \sqrt{z / D}$	1/2	1.027	1.525	2.276
$E_{s} (z) = E_{s D} z / D$	1	1.226	1.659	2.259

The resulting Winkler spring stiffnesses are plotted in Fig. 9 and compared against some existing solutions, including those of Francis (1964), who doubled the result of Vesic (1961) for a beam on an elastic foundation; Roesset (1980), who suggested the simple expression

k = 1.2 E_{s}

; and Syngros (2004), who provided separate expressions for

k_{F}

and

k_{P}

. In general, the iterative solution to Eq. (34) showed good agreement with the existing solutions, and the simplified approach proposed by Karatzia and Mylonakis (2017) in Eq. (35) provided very similar results. Also, for Eqs. (34) or (35), the differences in Figs. 9(a–c) are simply a stretch of the

E_{p} / E_{s D}

axis by a factor

{(b / b_{F})}^{(n + 4)}

. Therefore, by applying this factor to

E_{p} / E_{s D}

in existing formulas for fixed-head piles, approximate solutions for other load-moment combinations can be obtained.

Fig. 9. Winkler spring stiffnesses (a) $k_{F}$ ; (b) $k_{P}$ ; and (c) $k_{M}$ ( $ν_{s} = 0.4$ ).

For nonzero surface stiffness (

z_{0} > 0

),

b

values from numerical evaluation of the integral in Eq. (33) are shown in Fig. 10. By normalizing these results with

λ_{0}

[Fig. 10(b)] and letting the dimensionless surface stiffness increase, the results approached the values obtained for a homogeneous spring bed.

Fig. 10. Values of $b$ for a nonzero surface stiffness: (a) normalized by $λ$ ; and (b) normalized by $λ_{0}$ [Eq. (30)].

In summary, the approach detailed in this section to calculate

k (z)

consists of the following steps:

1.

Determine

Y_{1} (z)

,

Y_{3} (z)

, and

Y_{4} (z)

from Eq. (28) using the approximation

k = E_{s}

. Use the

λ

definition in Eq. (5c) as well as the substitution

z \to z + z_{0}

for nonzero surface stiffness.

2.

Input each of these into Eq. (33c) as

Y (z)

to get

b_{F}

,

b_{P}

, and

b_{M}

, respectively.

3.

Calculate

k_{F}

,

k_{P}

, and

k_{M}

, using Eq. (35) and the corresponding

b

term.

4.

Calculate

K_{11}

,

F_{11}

, and

F_{12}

using the solution in Eq. (29) and the corresponding

k

term. Use the modified solution for nonzero surface stiffness.

5.

Calculate the remaining stiffness and flexibility terms using Eq. (37).

If high precision is required, the improved

k

values obtained at Step 3 can be input back into Step 1 (replacing

k = E_{s}

) and the process repeated until they converge within a chosen tolerance. Eq. (34) can also be used in place of Eq. (35), although, as shown subsequently, neither of these changes are necessary to get good results. An illustrative example of applying this method to a free-head pile is presented in the Appendix.

The solution, shown here for power-law soil stiffness profiles, can be modified for other soil profiles (including multiple layers) by using appropriate shape functions in Step 1 and an appropriate Winkler solution in Step 4. The resulting

k (z) / E_{s} (z)

values are independent of depth, and this is a soil–structure interaction property depending on the full problem configuration.

Comparison with Continuum Results

The performance of the proposed solution for predicting pile head stiffness for homogeneous soil is shown in Fig. 11, where results are presented against available Winkler solutions and numerical continuum solutions (including some fitted functions). Employing Eq. (34) to predict

k_{F}

,

k_{P}

, and

k_{M}

provided predictions of

K_{11}

,

F_{11}

, and

F_{12}

in very good agreement with the continuum results as well as the recent Winkler spring stiffness calibrations for

K_{11}

and

F_{11}

by Syngros (2004).

Fig. 11. Predicted stiffness terms (a) $K_{11}$ ; (b) $K_{12}$ ; (c) $K_{22}$ ; (d) $1 / F_{11}$ ; (e) $1 / F_{12}$ ; and (f) $1 / F_{22}$ in a homogeneous soil [ $E_{s} (z) = E_{s}$ ] for different Winkler stiffness calibrations compared with continuum results and corresponding fitted functions ( $ν_{s} = 0.4$ ).

The simplified one-iteration approach using Eq. (35) performed similarly, often indistinguishable from the results obtained from Eq. (34). The remaining head stiffness terms, obtained from Eq. (37) for both methods also showed very good agreement with the continuum results. Although the more rigorous solution from Francis (1964) provided better results for

K_{11}

than the simplified finite-element based expression by Roesset (1980), both solutions highlighted the problem with using the fixed-head Winkler calibration for all head stiffness terms. This was most pronounced when calculating

F_{11}

[Fig. 11(d)], which has the highest dependence on Winkler spring stiffness, but is also evident in Figs. 11(b, e, and f).

The performance of the proposed solution for a linear soil stiffness profile [

E_{s} (z) = E_{s D} z / D

] is shown in Fig. 12. The results, employing Eq. (21), were similar to those for homogeneous soil. However, because each term has a lower power dependence on Winkler spring stiffness than the homogeneous case, the impact of differences in the Winkler spring stiffness employed was lower, and all results were more dependent on pile stiffness, resulting in less variation in both the Winkler and continuum results.

Fig. 12. Predicted stiffness term (a) $K_{11}$ ; (b) $K_{12}$ ; (c) $K_{22}$ ; (d) $1 / F_{11}$ ; (e) $1 / F_{12}$ ; and (f) $1 / F_{22}$ for a linearly varying stiffness profile [ $E_{s} (z) = E_{s D} z / D$ ] for different Winkler stiffness calibrations compared with continuum results and corresponding fitted functions ( $ν_{s} = 0.4$ ).

The predicted maximum moment,

M_{\max}

, for a free-head pile (

M = 0

) and a fixed-head pile (

θ_{0} = 0

) are shown in Fig. 13. For the former, this was obtained from the Winkler model using Eqs. (4b) and (9b) with the corresponding shape functions for each soil profile and was negative in the sign convention shown in Fig. 1(a). For homogeneous soils,

M_{\max}

occurs at a depth

z_{M} = π / 4 λ \approx 0.79 / λ

, and for the linearly increasing soil stiffness [

E_{s} (z) = E_{s D} z / D

],

M_{\max}

occurs at a depth

z_{M} \approx 0.96 / λ

(compared with numerical continuum results in Fig. 14). For the latter, the maximum moment occurs directly at the pile head and is the fixing moment,

M_{F}

. This is obtained directly from the stiffness matrix using the expression

M_{F} = P K_{12} / K_{11}

.

Fig. 13. Predicted maximum moment for free-head piles ( $M = 0$ ) and fixed-head piles ( $θ_{0} = 0$ ) embedded in (a) a homogeneous; and (b) a linearly varying [ $E_{s} (z) = E_{s D} z / D$ ] stiffness profile ( $ν_{s} = 0.4$ ).

Fig. 14. Predicted depth of maximum moment for free-head piles ( $M = 0$ ) embedded in (a) a homogeneous; and (b) a linearly varying [ $E_{s} (z) = E_{s D} z / D$ ] stiffness profile ( $ν_{s} = 0.4$ ).

The results for different Winkler spring stiffness calibrations are compared with functions fitted to available numerical continuum results in Figs. 13 and 14. In general, the Winkler model predicted higher peak moments than the continuum formulas for free-head piles, possibly due to the discretization in the continuum results missing the true peak. Results from this paper and fitted formulas by Syngros (2004) for both free- and fixed-head piles matched well. However, the calibration by Roesset (1980) (for fixed-head piles) overpredicted the maximum moments for the free-head case, particularly in the linear stiffness profile. The peak moment depth predictions for both profiles showed good agreement with the continuum results well, except for the function provided by Davies and Budhu (1986) for homogeneous soil, which was only fitted for

E_{p} / E_{s} < 10^{4}

and diverged significantly above this value.

Conclusions

The Winkler model for the response of laterally loaded piled foundations embedded in an inhomogeneous soil with stiffness varying according to a power law with depth [Eq. (1)] has been investigated. To this end, the following conclusions can be drawn:

•

Existing solutions for spring stiffness profiles varying linearly with depth were extended to obtain solutions for semi-infinite piles. This included head stiffness and flexibility matrices [Eq. (21) and Fig. 5], deflections, bending moments, and shear forces with depth [Eqs. (4), (9), and (20)]. The solutions were expressed in terms of generalized hypergeometric functions, which are available in modern programming languages.

•

A closed-form solution [Eq. (26)] was obtained for the general case of a power-law stiffness profile [Eq. (1) and

n > - 4

], also in terms of four generalized hypergeometric functions. To this end, the regular part of these functions (which does not diverge at infinity) was identified and extracted [Eq. (27)], and then used to calculate head stiffness and flexibility matrices [Eq. (29), Table 1 and Fig. 6], deflections, bending moments, and shear forces with depth [Eqs. (4), (9), and (28)] for long piles. This was shown to perfectly match a numerical Winkler solution evaluated for

n = 0.5

(Fig. 6).

•

A generalized Winkler parameter

λ

[Eq. (5c)] was introduced to allow convenient normalization of results and comparison between solutions for different stiffness profiles. Alternative normalizations were achieved using a modified Winkler parameter

λ_{0}

[Eq. (30)], which allowed the results from the aforementioned solutions to approach those for homogeneous soil as the surface stiffness increased (Fig. 7).

•

Four basic response functions obtained from the Winkler solution [Eq. (28)] were employed with an approximate energy solution for a continuum model to determine the normalized Winkler modulus

k / E_{s}

for inhomogeneous soils [Eqs. (33c), (34), and (35), and Fig. 10], for which few solutions are available. Because the Winkler modulus varies with boundary conditions at the pile head, a novel approach was introduced where three separate

k

values were calculated to obtain the stiffness and flexibility matrices covering all cases [Table 2, Fig. 9, and Eq. (37)]. This was shown to be a significant improvement compared with selecting a single

k

value (Fig. 8) and in much better agreement with numerical continuum results for predicting head stiffness terms and the peak bending moment (Figs. 11–14).

Notation

The following symbols are used in this paper:

$A_{m}$ , $B_{m}$ , $C_{m}$: integration constants determined from boundary conditions;
$a_{0}$ ( $= ω D / V_{s}$ ): dimensionless frequency;
$b$ , $b_{P}$ , $b_{M}$ , $b_{F}$: Winkler stiffness calibration constants;
$c_{m}$: complex contour;
$D$: pile diameter;
$E_{p}$ , $E_{s}$ , $E_{s D}$: Young’s modulus of the pile, (homogeneous) soil, and soil at one diameter depth;
$E_{s 1}$: rate of increase of soil modulus with depth;
$F$: pile head flexibility matrix;
${{}_{n}F}_{m} ()$: generalized hypergeometric function;
$F_{11}$ , $F_{12}$ , $F_{21}$ , $F_{22}$: pile head flexibility terms;
$f_{m} (z)$ , $g_{m} (z)$ , $h_{m} (z)$: solutions to governing equation;
$G_{s} (z)$: soil shear modulus;
$I$: pile second moment of area;
$K$: pile head stiffness matrix;
$K_{μ} ()$: modified Bessel function of the second kind and of order $μ$ ;
$K_{11}$ , $K_{12}$ , $K_{21}$ , $K_{22}$: pile head stiffness terms;
$k (z)$ , $k_{0}$ , $k_{D}$: Winkler stiffness variation and values the surface and one diameter depth;
$k_{P}$ , $k_{M}$ , $k_{F}$: Winkler stiffness piles under load-only, moment-only, and fixed-head conditions;
$L$ , $L_{a}$: pile length and active length;
$M (z)$ , $M$: bending moment distribution and applied head moment;
$M_{\max}$ , $M_{F}$: maximum bending moment and fixed-head moment;
$n$: Winkler stiffness profile exponent;
$P (z)$ , $P$: pile shear force distribution and applied head load;
$t$: complex variable;
$u$ , $P$: pile head displacement and load vector;
$u$ , $v$: radial and tangential soil displacement;
$Y (z)$ , $Y_{m} (z)$: displacement shape functions;
$y (z)$ , $y_{0}$: pile displacement variation and pile head displacement;
$z$: depth below ground level;
$z_{M}$: depth of maximum bending moment;
$z_{0}$: elevation above ground surface such that $k (z)$ is zero;
$Γ ()$: gamma function;
$γ$: Euler’s constant;
$η_{s}$: compressibility constant;
$θ (z)$ , $θ_{0}$: pile rotation variation and pile head rotation, respectively;
$λ$ , $λ_{0}$: Winkler wavelength parameters;
$ν$: $1 / (n + 4)$ ; and
$ν_{s}$: soil Poisson’s ratio.

Supplemental Materials

File (supplemental_materials_em.1943-7889.0002125_crispin.pdf)

Download
148.03 KB

Appendix. Illustrative Example

Scott (1981) analyzed the results of a number of pile tests undertaken in the Arkansas River. For one of the tests, Pile 2, Scott (1981) proposed predicting the pile response using a linear soil stiffness profile. To demonstrate the method proposed in this work, this test will be used as an illustrative example for predicting the pile head displacement at a given load. The pile is a 16-m-long steel pipe with 0.41-m external diameter, 7.9-mm wall thickness, and a reported

E_{p} I

of

69 {MNm}^{2}

. Based on knowledge of the site and the pile test results, Scott (1981) estimated the soil stiffness varies according to the following linear function of depth:

E_{s} (z) = E_{s 1} z

, where

E_{s 1} = 35 MN / m^{3}

is the rate of increase in

E_{s}

per meter depth. In addition, strain gauges showed that the bending moments were negligible at a depth of 4.6 m, indicating that the pile can be treated as a semi-infinite beam and the solution developed in this work can be employed.

The load was applied at a height of 30 mm, so the applied head moment can be neglected, and only

F_{11}

is required to predict the pile head displacement. This is obtained directly using

k_{p}

, corresponding to inputting

Y_{3} (z)

from Eq. (28) into Eq. (33c). From Table 2, this results in

b_{P} = 1.659 λ

, where

λ

is first estimated from Eq. (5c) [simplified for this case in Eq. (10c)] using

k (z) \approx E_{s} (z)

to give

λ = 0.63 m^{- 1}

. Inputting these into Eq. (35) yields

k_{P} (z) = 1.9 E_{s} (z)

, which corresponds to

λ = 0.72 m^{- 1}

. A single iteration was deemed suitable here.

Finally, the result in Eq. (29b) [simplified for his case in Eq. (21b)] can be employed with these values to give

F_{11} = 36 mm / MN

. This means that at an applied load of 191 kN, the pile is predicted to deflect 6.9 mm. The prediction of Scott (1981) was based on a value of

k = E_{s}

, which is closer to values employed for fixed-head piles; therefore, the resulting predicted head displacement of 10 mm is not directly comparable with the value obtained here.

Data Availability Statement

All data, models, and code generated or used during the study appear in the published article.

Acknowledgments

The first author would like to thank the Engineering and Physical Sciences Research Council for their support (Grant No. EP/N509619/1). The authors are grateful to the anonymous reviewers who provided constructive feedback that improved the work.

References

Anoyatis, G., and A. Lemnitzer. 2017. “Dynamic pile impedances for laterally–loaded piles using improved Tajimi and Winkler formulations.” Soil Dyn. Earthquake Eng. 92 (3): 279–297. https://doi.org/10.1016/j.soildyn.2016.09.020.

Google Scholar

Banerjee, P. K., and T. G. Davies. 1978. “The behaviour of axially and laterally loaded single piles embedded in nonhomogeneous soils.” Géotechnique 28 (3): 309–326. https://doi.org/10.1680/geot.1978.28.3.309.

Google Scholar

Basu, D., and R. Salgado. 2008. “Analysis of laterally loaded piles with rectangular cross sections embedded in layered soil.” Int. J. Numer. Anal. Methods Geomech. 32 (7): 721–744. https://doi.org/10.1002/nag.639.

Google Scholar

Basu, D., R. Salgado, and M. Prezzi. 2009. “A continuum-based model for analysis of laterally loaded piles in layered soils.” Géotechnique 59 (2): 127–140. https://doi.org/10.1680/geot.2007.00011.

Google Scholar

Bender, C. M., and S. A. Orszag. 1978. Advanced mathematical methods for scientists and engineers. New York: McGraw-Hill.

Google Scholar

Biot, M. A. 1937. “Bending of an infinite beam on an elastic foundation.” J. Appl. Mech. 4 (1): 1–7. https://doi.org/10.1115/1.4008739.

Google Scholar

Budhu, M., and T. G. Davis. 1988. “Analysis of laterally loaded piles in soft clays.” J. Geotech. Eng. 114 (1): 21–39. https://doi.org/10.1061/(ASCE)0733-9410(1988)114:1(21).

Google Scholar

Chang, Y. L. 1937. “Discussion on lateral pile-loading tests by Lawrence B. Feagin.” Trans. Am. Soc. Civ. Eng. 102 (1): 272–278. https://doi.org/10.1061/TACEAT.0004822.

Google Scholar

Crispin, J. J. 2017. Flexural elastic response of single piles in inhomogeneous soil. Undergraduate Research Rep. No. 1617RP005M. Bristol, UK: Univ. of Bristol.

Google Scholar

Crispin, J. J. 2022. “Static and dynamic analysis of piles in inhomogeneous soil.” Ph.D. thesis, Dept. of Civil Engineering, Univ. of Bristol.

Google Scholar

Crispin, J. J., C. P. Leahy, and G. Mylonakis. 2018. “Winkler model for axially loaded piles in inhomogeneous soil.” Géotech. Lett. 8 (4): 290–297. https://doi.org/10.1680/jgele.18.00062.

Google Scholar

Cucuzza, R., G. Devillanova, A. Aloisio, M. M. Rosso, and G. C. Marano. 2022. “Analytical solutions for piles' lateral deformations: The nonlinear stiffness case.” Int. J. Mech. Sci. 229: 107505. https://doi.org/10.1016/j.ijmecsci.2022.107505.

Google Scholar

Darendeli, M. B. 2001. “Development of a new family of normalized modulus reduction and material damping curves.” Ph.D. thesis, Univ. of Texas at Austin.

Google Scholar

Davies, T. G., and M. Budhu. 1986. “Non-linear analysis of laterally loaded piles in heavily over-consolidated clays.” Géotechnique 36 (4): 527–538. https://doi.org/10.1680/geot.1986.36.4.527.

Google Scholar

Di Laora, R., A. Mandolini, and G. Mylonakis. 2012. “Insight on kinematic bending of flexible piles in layered soil.” Soil Dyn. Earthquake Eng. 43 (52): 309–322. https://doi.org/10.1016/j.soildyn.2012.06.020.

Google Scholar

Di Laora, R., and E. Rovithis. 2015. “Kinematic bending of fixed-head piles in nonhomogeneous soil.” J. Geotech. Geoenviron. Eng. 141 (4): 04014126. https://doi.org/10.1061/(ASCE)GT.1943-5606.0001270.

Google Scholar

Dobry, R., E. Vicente, M. J. O’Rourke, and J. M. Roesset. 1982. “Horizontal stiffness and damping of single piles.” J. Geotech. Eng. 108 (3): 439–459. https://doi.org/10.1061/AJGEB6.0001259.

Google Scholar

El-Marsafawi, H., A. M. Kaynia, and M. Novak. 1992. Interaction factors and the superposition method for pile group dynamic analysis. London: Univ. of Western Ontario.

Google Scholar

Fradelos, J. J. 2016. “Pile stiffness matrix in homogeneous and inhomogeneous soil under transverse loading.” [In Greek.] Diploma thesis, Dept. of Civil Engineering, Univ. of Patras.

Google Scholar

Francis, A. J. 1964. “Analysis of pile groups with flexural resistance.” J. Soil Mech. Found. Div. 90 (3): 1–32. https://doi.org/10.1061/JSFEAQ.0000615.

Google Scholar

Franklin, J. N., and R. F. Scott. 1979. “Beam equation with variable foundation coefficient.” J. Eng. Mech. Div. 105 (5): 811–827. https://doi.org/10.1061/JMCEA3.0002525.

Google Scholar

Froio, D., and E. Rizzi. 2016. “Analytical solution for the elastic bending of beams lying on a variable Winkler support.” Acta Mech. 227 (4): 1157–1179. https://doi.org/10.1007/s00707-015-1508-y.

Google Scholar

Froio, D., and E. Rizzi. 2017. “Analytical solution for the elastic bending of beams lying on a linearly variable Winkler support.” Int. J. Mech. Sci. 128–129 (5): 680–694. https://doi.org/10.1016/j.ijmecsci.2017.04.021.

Google Scholar

Gazetas, G. 1991. “Foundation vibrations.” In Foundation engineering handbook, 2nd ed. New York: Van Nostrand Reinhold.

Crossref

Google Scholar

Gazetas, G., and R. Dobry. 1984. “Horizontal response of piles in layered soil.” J. Geotech. Eng. 110 (1): 20–40. https://doi.org/10.1061/(ASCE)0733-9410(1984)110:1(20).

Google Scholar

Gibson, R. E. 1974. “The analytical method in soil mechanics.” Géotechnique 24 (2): 115–140. https://doi.org/10.1680/geot.1974.24.2.115.

Google Scholar

Gleser, S. M. 1953. “Lateral load tests on vertical fixed-head and free-head piles.” In Proc., Symp. on Lateral Load Tests on Piles: Presented at the Fifty-sixth Annual Meeting, 75–93. Reston, VA: ASTM.

Google Scholar

Guo, W. D. 2012. Theory and practice of pile foundations. New York: Taylor & Francis.

Crossref

Google Scholar

Guo, W. D., and F. H. Lee. 2001. “Load transfer approach for laterally loaded piles.” Int. J. Numer. Anal. Methods Geomech. 25 (11): 1101–1129. https://doi.org/10.1002/nag.169.

Google Scholar

Hayashi, K. 1921. Theorie des trägers auf elastischer unterlage: Und ihre anwendung auf den tiefbau, nebst einer tafel der kreisund hyperbelfunktionen. Berlin: Springer.

Google Scholar

Hetenyi, M. 1946. Beams on elastic foundation. Ann Arbor, MI: Univ. of Michigan Press.

Google Scholar

Karatzia, X., and G. Mylonakis. 2012. “Horizontal response of piles in inhomogeneous soil: Simple analysis.” In Proc., 2nd Int. Conf. on Performance Based Design in Earthquake Geotechnical Engineering. London: International Society for Soil Mechanics and Geotechnical Engineering.

Google Scholar

Karatzia, X., and G. Mylonakis. 2017. “Horizontal stiffness and damping of piles in inhomogeneous soil.” J. Geotech. Geoenviron. Eng. 143 (4): 04016113. https://doi.org/10.1061/(ASCE)GT.1943-5606.0001621.

Google Scholar

Lam, I. P., and G. Martin. 1986. Seismic design for highway bridge foundations. Washington, DC: USDOT.

Google Scholar

Liang, F., Y. Li, L. Li, and J. Wang. 2014. “Analytical solution for laterally loaded long piles based on Fourier–Laplace integral.” Appl. Math. Modell. 38 (21): 5198–5216. https://doi.org/10.1016/j.apm.2014.03.052.

Google Scholar

Matlock, H., and L. C. Reese. 1960. “Generalized solutions for laterally loaded piles.” J. Soil Mech. Found. Div. 86 (5): 63–92. https://doi.org/10.1061/JSFEAQ.0000303.

Google Scholar

Miche, R. J. 1930. “Investigation of piles subject to horizontal forces, application to quay walls.” J. School Eng. Giza 4.

Google Scholar

Mylonakis, G. 1995. “Contributions to static and seismic analysis of piles and pile-supported bridge piers.” Ph.D. thesis, Dept. of Civil Engineering, State Univ. of New York at Buffalo.

Google Scholar

Mylonakis, G. 2001. “Elastodynamic model for large-diameter end-bearing shafts.” Soils Found. 41 (3): 31–44. https://doi.org/10.3208/sandf.41.3_31.

Google Scholar

Mylonakis, G. 2004. Analytical investigation of flexural pile response in inhomogeneous soil. Unpublished Research Report. Patras, Greece: Univ. of Patras.

Google Scholar

Mylonakis, G., and G. Gazetas. 1999. “Lateral vibration and internal forces of grouped piles in layered soil.” J. Geotech. Geoenviron. Eng. 125 (1): 16–25. https://doi.org/10.1061/(ASCE)1090-0241(1999)125:1(16).

Google Scholar

Mylonakis, G., and D. Roumbas. 2001. “Lateral impedance of single piles in inhomogeneous soil.” In Proc., 4th Int. Conf. on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics. Paper 6.27. Rolla, MI: Univ. of Missouri-Rolla, CD-ROM.

Google Scholar

Mylonakis, G. E., and J. J. Crispin. 2021. “Simplified models for lateral static and dynamic analysis of pile foundations.” In Analysis of pile foundations subject to static and dynamic loading, 185–245. London: CRC Press. https://doi.org/10.1201/9780429354281-6.

Google Scholar

Novak, M. 1974. “Dynamic stiffness and damping of piles.” Can. Geotech. J. 11 (4): 574–598. https://doi.org/10.1139/t74-059.

Google Scholar

Olver, F. W. J., D. W. Lozier, R. F. Boisvert, and C. W. Clark. 2010. NIST handbook of mathematical functions. New York: Cambridge University Press.

Google Scholar

Palmer, L. A., and J. B. Thompson. 1948. “The earth pressure and deflection along the embedded lengths of piles subjected to lateral thrust.” In Proc., 2nd Int. Conf. on Soil Mechanics and Foundation Engineering, 156–161. London: International Society for Soil Mechanics and Geotechnical Engineering.

Google Scholar

Parashakis, H. 2020. “Stress wave propagation in inhomogeneous soil media.” Ph.D. thesis, Dept. of Civil Engineering, Univ. of Patras.

Google Scholar

Pender, M. 1993. “Aseismic pile foundation design analysis.” Bull. N. Z. Nat. Soc. Earthquake Eng. 26 (1): 49–160. https://doi.org/10.5459/bnzsee.26.1.49-160.

Google Scholar

Randolph, M. F. 1981. “The response of flexible piles to lateral loading.” Géotechnique 31 (2): 247–259. https://doi.org/10.1680/geot.1981.31.2.247.

Google Scholar

Reese, L. C., and W. F. Van Impe. 2011. Single piles and pile groups under lateral loading. 2nd ed. London: CRC Press.

Google Scholar

Rifaat, I. 1935. “Die Spundwand als Erddruckproblem. Das Spundwandproblem mit Berücksichtigung der Erddeformation und der Wandelastizität.” Ph.D. thesis, Swiss Federal Institute of Technology, Zurich.

Google Scholar

Roesset, J. M. 1980. “Stiffness and damping coefficients of foundations.” In Dynamic response of pile foundations, 1–29. New York: ASCE.

Google Scholar

Rosendo, D. C., and P. J. R. Albuquerque. 2021. “General analytical solution for laterally-loaded pile-based Miche model.” Geotech. Geol. Eng. 39 (2): 765–782. https://doi.org/10.1007/s10706-020-01520-1.

Google Scholar

Scott, R. F. 1981. Foundation analysis. Englewood Cliffs, NJ: Prentice Hall.

Google Scholar

Shadlou, M., and M. Bhattacharya. 2014. “Dynamic stiffness of pile in a layered elastic continuum.” Géotechnique 64 (4): 303–319. https://doi.org/10.1680/geot.13.P.107.

Google Scholar

Syngros, K. 2004. “Seismic response of piles and pile-supported bridge piers evaluated through case histories.” Ph.D. thesis, City College of the City Univ. of New York.

Google Scholar

Trachanas, S. 1989. Differential equations. [In Greek.] Crete, Greece: Crete University Press.

Google Scholar

Vallabhan, C. V. G., and Y. C. Das. 1988. “Parametric study of beams on elastic foundations.” J. Eng. Mech. 111 (5): 664–679. https://doi.org/10.1061/(ASCE)0733-9399(1988)114:12(2072).

Google Scholar

Veletsos, A. S., and A. H. Younan. 1994. “Dynamic soil pressures on rigid cylindrical vaults.” Earthquake Eng. Struct. Dyn. 23 (6): 645–669. https://doi.org/10.1002/eqe.4290230606.

Google Scholar

Vesić, A. B. 1961. “Beams on elastic subgrade and the Winkler’s hypothesis.” In Proc., 5th Int. Conf. on Soil Mechanics and Foundation Engineering. London: International Society for Soil Mechanics and Geotechnical Engineering.

Google Scholar

Vlasov, V. Z., and N. N. Leontiev. 1966. Beams, plates and shells on elastic foundations. Washington, DC: Israel Program for Scientific Translations.

Google Scholar

Yoshida, I., and R. Yoshinaka. 1972. “A method to estimate modulus of horizontal subgrade reaction for a pile.” Soils Found. 12 (3): 1–17. https://doi.org/10.3208/sandf1972.12.3_1.

Google Scholar

Information & Authors

Information

Published In

Journal of Engineering Mechanics

Volume 148 • Issue 11 • November 2022

Copyright

This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.

History

Received: Nov 6, 2021

Accepted: Mar 26, 2022

Published online: Sep 6, 2022

Published in print: Nov 1, 2022

Discussion open until: Feb 6, 2023

Authors

Affiliations

Jamie J. Crispin, Ph.D. https://orcid.org/0000-0003-3074-8493 [email protected]

Formerly, Ph.D. Student, Dept. of Civil Engineering, Univ. of Bristol, Bristol BS8 1TR, UK (corresponding author). ORCID: https://orcid.org/0000-0003-3074-8493. Email: [email protected]

View all articles by this author

George Mylonakis, Ph.D., M.ASCE

Professor of Geotechnics, Dept. of Civil Infrastructure and Environmental Engineering, Khalifa Univ., United Arab Emirates; Chair in Geotechnics and Soil-Structure Interaction, Dept. of Civil Engineering, Univ. of Bristol, Bristol BS8 1TR, UK; Adjunct Professor, Dept. of Civil and Environmental Engineering, Univ. of California, Los Angeles, CA 90024.

View all articles by this author

Metrics & Citations

Metrics

Citations

Download citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Exact Winkler Solution for Laterally Loaded Piles in Inhomogeneous Soil

Abstract

Introduction

The Problem

Homogeneous Stiffness Profile

Linearly Varying Stiffness Profile

Power-Law Stiffness Profile

Assessment of Winkler Spring Stiffness

Comparison with Continuum Results

Conclusions

Notation

Supplemental Materials

Appendix. Illustrative Example

Data Availability Statement

Acknowledgments

References

Information & Authors

Information

Published In

Copyright

History

Authors

Affiliations

Metrics & Citations

Metrics

Citations

Download citation

View Options

Media

Figures

Other

Tables

PREVIOUS ARTICLE

NEXT ARTICLE

Verify Phone

Congrats!

Abstract

Introduction

The Problem

Homogeneous Stiffness Profile

Linearly Varying Stiffness Profile

Power-Law Stiffness Profile

Assessment of Winkler Spring Stiffness

Comparison with Continuum Results

Conclusions

Notation

Supplemental Materials

Appendix. Illustrative Example

Data Availability Statement

Acknowledgments

References

Information

Published In

Copyright

History

Authors

Affiliations

Metrics

Citations

Download citation

Figures

Other

Share

Copy the content Link

Share with email

Share

Request Username

Create a new account

Change Password

Password Changed Successfully

Verify Phone

Congrats!