AbstractA novel exact analytical solution is derived for the equation y(4)+xny=0 in the region x≥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.IntroductionThe 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)+xny=0, where x (≥  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≠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 F30 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, Es (or Gs). 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 ProblemThe 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 Ep, length L, diameter D, and second moment of area I (=πD4/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).In this paper, the following power-law Winkler stiffness profile is considered, shown in Fig. 1(b)(1a) where kD = Winkler spring stiffness at a depth of one diameter; n = power-law exponent; and z0 = characteristic length corresponding to the elevation above the pile head that k(z) (if extrapolated above the ground surface) reaches zero. z0 is shown graphically in Fig. 1(b) and is given by (1b) where k0 = 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 z0, 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 y0 and rotation (anticlockwise positive) θ0, which can be related to the applied loading using either a (2×2) stiffness matrix K or a corresponding flexibility matrix F as shown in Eq. (2) (2a) P=[PM]=Ku=[K11K12K21K22][y0θ0](2b) u=[y0θ0]=FP=[F11F12F21F22][P0M0]where K11, K22, and K12 (= K21) = swaying, rocking, and cross-swaying rocking stiffness coefficients, respectively; and F11, F22, and F12 (= F21) = corresponding flexibilities. These matrices are related by the simple expression (2c) Because both the pile and soil are modeled as linearly elastic, both matrices are symmetric by reciprocity (i.e., K12=K21 and F12=F21). 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 K12 and K21 are positive, whereas the off-diagonal flexibility terms F12 and F21 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) (3) EpIy⁗(z)=q(z)=−k(z)y(z)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(4a) (4b) (4c) Inputting the Winkler stiffness profile under consideration [Eq. (1a)] into Eq. (3) yields the following governing equation(5a) EpIy⁗(z)+kD(z+z0D+z0)ny(z)=0This can be simplified by introducing a Winkler parameter, λ (carrying units of length−1) such that (5b) y⁗(z)+(n+4)λn+4[z0+z]ny(z)=0(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 kD 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 k0 and z0 directly, neither of which depend on the choice of reference depth (5d) λn+4=(kD1/n−k01/n)n(n+4)EpIDn=k0(n+4)EpIz0nWith the introduction of λ, only three dimensionless parameters are required to describe any problem modeled as in Fig. 1: λL, λz0, and n. This is evident from dimensional analysis, considering the five relevant problem parameters (k0, z0, n, EI, 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 λz0 incorporates a pile-soil stiffness ratio dependent on k0, z0, and n and given by (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, λz0 can be interpreted as a dimensionless surface stiffness parameter.Homogeneous Stiffness ProfileA 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) (6a) (6b) (6c) y(z)=e−λz[A1cos(λz)+A2sin(λz)]+e+λz[A3cos(λz)+A4sin(λz)]where Am 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).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 La, 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 A3 and A4 must be set to zero for long piles, significantly simplifying the solution.Various definitions for La are available such as λLa=22, 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: •Y1(z), the displaced shape due to unit displacement at the pile head under zero rotation,•Y2(z), the displaced shape due to unit rotation at the pile head under zero displacement,•Y3(z), the displaced shape due to unit load at the pile head under zero moment, and•Y4(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, g1(z), and g2(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 A1 and A2 in Eq. (6) [i.e., g1(z)=e−λzcos(λz), g2(z)=e−λzsin(λz)], or any linear combination of them, and result in the expressions in Eq. (7) (Hetenyi 1946; Pender 1993). With the exception of Y1(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) (7a) Y1(z)=g1(z)g2′(0)−g2(z)g1′(0)g1(0)g2′(0)−g2(0)g1′(0)=e−λz[cos(λz)+sin(λz)](7b) Y2(z)=g1(z)g2(0)−g2(z)g1(0)g1′(0)g2(0)−g2′(0)g1(0)=1λe−λzsin(λz)(7c) Y3(z)=1EIg1(z)g2″(0)−g2(z)g1″(0)g1‴(0)g2″(0)−g2‴(0)g1″(0)=12EpIλ3e−λzcos(λz)(7d) Y4(z)=1EIg1(z)g2‴(0)−g2(z)g1‴(0)g1″(0)g2‴(0)−g2″(0)g1‴(0)=−12EpIλ2e−λz[cos(λz)−sin(λz)]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(8a) K=EpI[Y1‴(0)Y2‴(0)−Y1″(0)−Y2″(0)]=EpI[4λ32λ22λ22λ](8b) F=[Y3(0)Y4(0)Y3′(0)Y4′(0)]=12EpI[λ−3−λ−2−λ−22λ−1]Once these have been used to establish the pile head deflection y0 and rotation θ0, the overall displaced shape of the pile can be readily obtained from Eqs. (9a) or (9b)(9a) (9b) which can also be input into Eq. (4) to obtain the bending moment and shear force variation with depth.Linearly Varying Stiffness ProfileAnalysis 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 z0=0. In this case, Eq. (5) can be rewritten (10a) (10b) (10c) Hayashi (1921) and later Rifaat (1935) and Hetenyi (1946) proposed the following solution to Eq. (10): (11a) y(z)=B1h1(z)+B2h2(z)+B3h3(z)+B4h4(z)where the Bm terms = coefficients to be determined from the boundary conditions; and the hm(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)]: (11b) h1(z)=(λz)00!F03(;45,35,25;−(λz)553)(11c) h2(z)=(λz)11!F03(;45,35,65;−(λz)553)(11d) h3(z)=(λz)22!F03(;45,75,65;−(λz)553)(11e) h4(z)=(λz)33!F03(;85,75,65;−(λz)553)where F03(;b1,b2,b3;x) = generalized hypergeometric function of argument x and parameters b1, b2, and b3 (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 cm, 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.(12) fm(z)=∫cmexp(515tλz+t55)dtRepeated differentiation of this equation yields the expression for the jth derivative (13) fm(j)(z)=5j5λj∫cmtjexp(515tλz+t55)dtPerforming integration by parts on the fourth derivative as follows yields (14) fm⁗(z)=545λ4[exp(515tλz+t55)]cm−5λ5z∫cmexp(515tλz+t55)dtThis shows that if cm is chosen such that the integrated part (first term on the right-hand side) vanishes, y(z)=fm(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.This results in five functions, corresponding to computing fm(z) from Eq. (12) with each contour cm. 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[fm(z)], and imaginary, I[fm(z)], components: •R[f0(z)]=−R[f4(z)],•I[f0(z)]=+I[f4(z)],•R[f1(z)]=−R[f3(z)], and•I[f1(z)]=+I[f3(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 (15) y(z)=C1R[f0(z)]+C2I[f0(z)]+C3R[f1(z)]+C4I[f1(z)]where the Cm terms = coefficients to be determined from the boundary conditions.Taking the contours to follow exactly the limits of their defined region, the jth derivative of f0(z) and f1(z) are given by [modified from Franklin and Scott (1979)] (16a) f0(j)(z)=−eπi(j+1)5j5λj∫0∞tjexp(515eπitλz−t55)dt+e3πi5(j+1)5j5λj∫0∞tjexp(515e3πi5tλz−t55)dt(16b) f1(j)(z)=−e3πi5(j+1)5j5λj∫0∞tjexp(515e3πi5tλz−t55)dt+eπi5(j+1)5j5λj∫0∞tjexp(515eπi5tλz−t55)dtScott (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=51/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 f0(z) approach zero as z increases, and f1(z) is unbounded. In addition, f0(z) and f1(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 Bm terms in Eq. (11a) such that y(z)=f0(z) or f1(z)]. Both f0(z) and f1(z) can be evaluated analytically in terms of the gamma function Γ() at z=0 (Franklin and Scott 1979) (17a) f0(j)(0)=52j−45λjΓ(j+15)[−eπi(j+1)+e3πi5(j+1)](17b) f1(j)(0)=52j−45λjΓ(j+15)[−e3πi5(j+1)+eπi5(j+1)]The power-series solution in Eq. (11) has the following convenient property at z=0: (18) where y(j)(z) = jth 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 Bm terms and, consequently, analytical expressions for Eq. (16) can be obtained (19a) f0(z)=5−45[1+(−1)35]Γ(15)h1(z)−5−25[1+(−1)15]Γ(25)h2(z)+[1+(−1)−15]Γ(35)h3(z)−525[1+(−1)−35]Γ(45)h4(z)(19b) f1(z)=5−45[(−1)35−(−1)15]Γ(15)h1(z)+5−25[(−1)65−(−1)25]Γ(25)h2(z)+[(−1)95−(−1)35]Γ(35)h3(z)+525[(−1)25−(−1)45]Γ(45)h4(z)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 g1(z)=R[f0(z)] and g2(z)=I[f0(z)], yields the four basic response functions for n=1 as follows(20a) Y1(z)=h1(z)−5451+52Γ(35)Γ(15)h3(z)+5651+52Γ(45)Γ(15)h4(z)(20b) λY2(z)=h2(z)−5251+52Γ(35)Γ(25)h3(z)+545Γ(45)Γ(25)h4(z)(20c) EpIλ3Y3(z)=5−651+52Γ(15)Γ(45)h1(z)−5−451+52Γ(25)Γ(45)h2(z)+h4(z)(20d) EpIλ2Y4(z)=−5−45Γ(15)Γ(35)h1(z)+5−251+52Γ(25)Γ(35)h2(z)−h3(z)where h1(z), h2(z), h3(z), and h4(z) = hypergeometric functions in Eq. (11). It is worth stressing that although all four functions hm(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 (21a) K=EpI1+52[565Γ(45)Γ(15)λ3545Γ(35)Γ(15)λ2545Γ(35)Γ(15)λ2525Γ(35)Γ(25)λ]≅EpI[2.831λ31.902λ21.902λ22.068λ](21b) F=1EpI1+52[5−65λ3Γ(15)Γ(45)−5−45λ2Γ(25)Γ(45)−5−45λ2Γ(25)Γ(45)5−25λΓ(25)Γ(35)]≅1EpI[0.925λ−3−0.851λ−2−0.851λ−21.266λ−1]In which λ5=kD/(5 EpID) [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/Ep)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 Ep and a lower dependence on k (and therefore Es) 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 k0, Eq. (10) must be modified as follows(22a) (22b) (22c) where z0, 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→z+z0. 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 g1(z) and g2(z) taken as R[f0(z+z0)] and I[f0(z+z0)], 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.Power-Law Stiffness ProfileFor 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 (z0=0). Firstly, a simple monomial trial solution y(z)=zs can be substituted into Eq. (5b) to provide the indicial equation (23) s(s−1)(s−2)(s−3)zs−4+(n+4)λn+4zs+n=0From 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(24) y(z)=zs∑m=0∞αmzm(n+4)where αm = coefficient of the mth term. Inputting this solution into Eq. (5b) yields the relationship between sequential αm coefficients as follows (25) αm+1=−(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]αmwhich can be used to derive the following four solutions to Eq. (5b) that can be used with Eq. (11a) (26a) h1(z)=(λz)00!F03(;1−ν,1−2ν,1−3ν;−ν3(λz)1ν)(26b) h2(z)=(λz)11!F03(;1−ν,1−2ν,1+ν;−ν3(λz)1ν)(26c) h3(z)=(λz)22!F03(;1−ν,1+2ν,1+ν;−ν3(λz)1ν)(26d) h4(z)=(λz)33!F03(;1+3ν,1+2ν,1+ν;−ν3(λz)1ν)where (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 K11 and K22 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 K11 and K22 than a pile fixed or hinged at the tip. This, however, is possible for the cross term K12.•The opposite naturally holds for the flexibility terms F11, F22, and F12.•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 Bm 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 (27a) g1(z)=h1(z)Γ(ν)Γ(2ν)Γ(3ν)−cos(πν)ν1−3νh2(z)Γ(ν)Γ(2ν)Γ(1−ν)+cos(2πν)ν2−6νh3(z)Γ(ν)Γ(1−2ν)Γ(1−ν)−cos(3πν)ν3−9νh4(z)Γ(1−3ν)Γ(1−2ν)Γ(1−ν)(27b) g2(z)=sin(πν)ν1−3νh2(z)Γ(ν)Γ(2ν)Γ(1−ν)−sin(2πν)ν2−6νh3(z)Γ(ν)Γ(1−2ν)Γ(1−ν)+sin(3πν)ν3−9νh4(z)Γ(1−3ν)Γ(1−2ν)Γ(1−ν)The basic response functions are obtained, as before, by substituting Eq. (27) into the relationships in Eq. (7) (28a) Y1(z)=h1(z)−ν6ν−2Γ(1−2ν)Γ(1−ν)Γ(2ν)Γ(3ν)h3(z)+ν9ν−3Γ(1−3ν)Γ 2(1−ν)Γ 2(2ν)Γ(3ν)h4(z)(28b) λY2(z)=h2(z)−ν3ν−1Γ(1−ν)Γ(ν)Γ 2(2ν)h3(z)+ν6ν−2Γ(1−2ν)Γ(1−ν)Γ(2ν)Γ(3ν)h4(z)(28c) EpIλ3Y3(z)=ν3−9νΓ 2(ν)Γ(3ν)Γ(1−3ν)Γ 2(1−2ν)h1(z)−ν2−6νΓ(ν)Γ(2ν)Γ(1−3ν)Γ(1−2ν)h2(z)+h4(z)(28d) EpIλ2Y4(z)=−ν2−6νΓ(ν)Γ(2ν)Γ(1−3ν)Γ(1−2ν)h1(z)+ν1−3νΓ(1−ν)Γ(ν)Γ 2(2ν)h2(z)−h3(z)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(29a) K=EpI[ν9ν−3Γ(1−3ν)Γ 2(1−ν)Γ 2(2ν)Γ(3ν)λ3ν6ν−2Γ(1−2ν)Γ(1−ν)Γ(2ν)Γ(3ν)λ2ν6ν−2Γ(1−2ν)Γ(1−ν)Γ(2ν)Γ(3ν)λ2ν3ν−1Γ(ν)Γ(1−ν)Γ 2(2ν)λ](29b) F=1EpI[ν3−9νλ3Γ 2(ν)Γ(3ν)Γ(1−3ν)Γ 2(1−2ν)−ν2−6νλ2Γ(ν)Γ(2ν)Γ(1−3ν)Γ(1−2ν)−ν2−6νλ2Γ(ν)Γ(2ν)Γ(1−3ν)Γ(1−2ν)ν1−3νλΓ(ν)Γ(1−ν)Γ 2(1−2ν)]where λn+4=kD[(n+4)EpIDn] [given in Eqs. (5c) or (5d) with k0=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)=kD(z/D)n with zero surface stiffnessTable 1. Pile head stiffness terms for a power-law Winkler spring stiffness profile k(z)=kD(z/D)n with zero surface stiffnessnK11EpIλ3K12EpIλ2K22EpIλ1F11EpIλ−1F12EpIλ1F12EpIλ04.0002.0002.0002.0002.0001.0001/43.4911.9532.0151.5981.6490.9221/23.1751.9242.0321.3531.4280.8663/42.9691.9082.0491.1921.2800.82312.8311.9022.0681.0811.1760.7903/22.6741.9092.1060.9451.0420.74422.6091.9312.1450.8700.9660.715Nonzero surface stiffness can be incorporated, as previously, with the simple substitution z→z+z0. The response functions must be recalculated from the relationships in Eq. (7) with g1(z) and g2(z) taken as g1(z+z0) and g2(z+z0), 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.As the dimensionless surface stiffness k0 and the associated dimensionless parameter λz0 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 z0, 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 (30) λ04=k04EpI=λ4n+44(λz0)nAssessment of Winkler Spring StiffnessApplication 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)(31a) ηs2Gs(z)∂∂r[1r(∂(ru)∂r+∂v∂θ)]−Gs(z)1r2∂∂θ[∂(rv)∂r−∂u∂θ]+∂∂z[Gs(z)∂u∂z]=0(31b) Gs(z)∂∂r[1r(∂(ru)∂r−∂v∂θ)]+ηs2Gs(z)1r2∂∂θ[∂(rv)∂r+∂u∂θ]+∂∂z[Gs(z)∂u∂z]=0where 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; Gs(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 ηs2=(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(32a) u(r,z,θ)=u(r,θ)Y(z)/Y(0)(32b) v(r,z,θ)=v(r,θ)Y(z)/Y(0)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(33a) ηs2∂∂r[1r(∂(ru)∂r+∂v∂θ)]−1r2∂∂θ[∂(rv)∂r−∂u∂θ]−b2u=0(33b) ∂∂r[1r(∂(ru)∂r−∂v∂θ)]+ηs21r2∂∂θ[∂(rv)∂r+∂u∂θ]−b2v=0where 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) (33c) b2=λ2∫0∞Gs(λz)[dY(λz)d(λz)]2d(λz)/∫0∞Gs(λz)[Y(λz)]2d(λz)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) (34) k/Es=πs22(1+νs)4K1(q)K1(s)+sK1(q)K0(s)+qK0(q)K1(s)qsK0(q)K0(s)+sK1(q)K0(s)+qK0(q)K1(s)where s = bD/2 (instead of the dimensionless frequency ao 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≈Es. 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(35) k/Es=2πηs2/(1+νs)ln(ηs)−(1+ηs2)ln(χbD)where χ=eγ/4≈0.445; ηs2=(2−νs)/(1−νs); and γ (≈ 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)=Y1(z) from Eq. (7a) results in b=λ2/3; similarly, for a free-head pile (under load only), using Y(z)=Y3(z) from Eq. (7a) results in b=λ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)(36) (b/λ)2=6P2−16PλM+12(λM)23P2−4PλM+2(λM)2The 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, kF. 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)/y0 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), kF, an applied load only (free-head and zero moment), kP, and an applied moment only (free-head and zero lateral force), kM, the corresponding stiffness terms K11, F11, and F12 can be obtained, and the remaining terms calculated as follows(37a) (37b) K12=K21=(F11K11−1)/(−F12)(37c) Here, kP, kM, and kF are determined using their corresponding b terms, namely bF, bP, and bM, which are, in turn, calculated from Eq. (33c) using functions Y1(z), Y3(z), and Y4(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: bF (fixed-head), bP (free-head, force only), and bM (free-head, moment only)Table 2. Values of parameter b [Eq. (33)] for different soil stiffness profiles and boundary conditions at the pile head: bF (fixed-head), bP (free-head, force only), and bM (free-head, moment only)Stiffness profilenbF/λbP/λbM/λEs(z)=Es00.816 (2/3)1.414 (2)2.449 (6)Es(z)=EsDz/D1/21.0271.5252.276Es(z)=EsDz/D11.2261.6592.259The 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.2Es; and Syngros (2004), who provided separate expressions for kF and kP. 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 Ep/EsD axis by a factor (b/bF)(n+4). Therefore, by applying this factor to Ep/EsD in existing formulas for fixed-head piles, approximate solutions for other load-moment combinations can be obtained.For nonzero surface stiffness (z0>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.In summary, the approach detailed in this section to calculate k(z) consists of the following steps: 1.Determine Y1(z), Y3(z), and Y4(z) from Eq. (28) using the approximation k=Es. Use the λ definition in Eq. (5c) as well as the substitution z→z+z0 for nonzero surface stiffness.2.Input each of these into Eq. (33c) as Y(z) to get bF, bP, and bM, respectively.3.Calculate kF, kP, and kM, using Eq. (35) and the corresponding b term.4.Calculate K11, F11, and F12 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=Es) 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)/Es(z) values are independent of depth, and this is a soil–structure interaction property depending on the full problem configuration.Comparison with Continuum ResultsThe 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 kF, kP, and kM provided predictions of K11, F11, and F12 in very good agreement with the continuum results as well as the recent Winkler spring stiffness calibrations for K11 and F11 by Syngros (2004).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 K11 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 F11 [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 [Es(z)=EsDz/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.The predicted maximum moment, Mmax, 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, Mmax occurs at a depth zM=π/4λ≈0.79/λ, and for the linearly increasing soil stiffness [Es(z)=EsDz/D], Mmax occurs at a depth zM≈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, MF. This is obtained directly from the stiffness matrix using the expression MF=PK12/K11.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 Ep/Es<104 and diverged significantly above this value.ConclusionsThe 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/Es 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).Appendix. Illustrative ExampleScott (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 EpI of 69  MNm2. 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: Es(z)=Es1z, where Es1=35  MN/m3 is the rate of increase in Es 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 F11 is required to predict the pile head displacement. This is obtained directly using kp, corresponding to inputting Y3(z) from Eq. (28) into Eq. (33c). From Table 2, this results in bP=1.659λ, where λ is first estimated from Eq. (5c) [simplified for this case in Eq. (10c)] using k(z)≈Es(z) to give λ=0.63  m−1. Inputting these into Eq. (35) yields kP(z)=1.9Es(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 F11=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=Es, 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 StatementAll data, models, and code generated or used during the study appear in the published article.AcknowledgmentsThe 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.NotationThe following symbols are used in this paper:Am, Bm, Cm=integration constants determined from boundary conditions;a0 (=ωD/Vs)=dimensionless frequency;b, bP, bM, bF=Winkler stiffness calibration constants;cm=complex contour;D=pile diameter;Ep, Es, EsD=Young’s modulus of the pile, (homogeneous) soil, and soil at one diameter depth;Es1=rate of increase of soil modulus with depth;F=pile head flexibility matrix;Fnm()=generalized hypergeometric function;F11, F12, F21, F22=pile head flexibility terms;fm(z), gm(z), hm(z)=solutions to governing equation;Gs(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 μ;K11, K12, K21, K22=pile head stiffness terms;k(z), k0, kD=Winkler stiffness variation and values the surface and one diameter depth;kP, kM, kF=Winkler stiffness piles under load-only, moment-only, and fixed-head conditions;L, La=pile length and active length;M(z), M=bending moment distribution and applied head moment;Mmax, MF=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), Ym(z)=displacement shape functions;y(z), y0=pile displacement variation and pile head displacement;z=depth below ground level;zM=depth of maximum bending moment;z0=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.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. 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. 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. 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. Bender, C. M., and S. A. Orszag. 1978. Advanced mathematical methods for scientists and engineers. New York: McGraw-Hill. Crispin, J. J. 2017. Flexural elastic response of single piles in inhomogeneous soil. Undergraduate Research Rep. No. 1617RP005M. Bristol, UK: Univ. of Bristol. Crispin, J. J. 2022. “Static and dynamic analysis of piles in inhomogeneous soil.” Ph.D. thesis, Dept. of Civil Engineering, Univ. of Bristol. 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. 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. 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. 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. 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. 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. Gazetas, G. 1991. “Foundation vibrations.” In Foundation engineering handbook, 2nd ed. New York: Van Nostrand Reinhold. 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. Guo, W. D. 2012. Theory and practice of pile foundations. New York: Taylor & Francis. 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. Hayashi, K. 1921. Theorie des trägers auf elastischer unterlage: Und ihre anwendung auf den tiefbau, nebst einer tafel der kreisund hyperbelfunktionen. Berlin: Springer. Hetenyi, M. 1946. Beams on elastic foundation. Ann Arbor, MI: Univ. of Michigan Press. 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. Lam, I. P., and G. Martin. 1986. Seismic design for highway bridge foundations. Washington, DC: USDOT. 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. Miche, R. J. 1930. “Investigation of piles subject to horizontal forces, application to quay walls.” J. School Eng. Giza 4. 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. Mylonakis, G. 2004. Analytical investigation of flexural pile response in inhomogeneous soil. Unpublished Research Report. Patras, Greece: Univ. of Patras. 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. 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. 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. 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. Parashakis, H. 2020. “Stress wave propagation in inhomogeneous soil media.” Ph.D. thesis, Dept. of Civil Engineering, Univ. of Patras. Reese, L. C., and W. F. Van Impe. 2011. Single piles and pile groups under lateral loading. 2nd ed. London: CRC Press. 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. Roesset, J. M. 1980. “Stiffness and damping coefficients of foundations.” In Dynamic response of pile foundations, 1–29. New York: ASCE. Scott, R. F. 1981. Foundation analysis. Englewood Cliffs, NJ: Prentice Hall. 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. Trachanas, S. 1989. Differential equations. [In Greek.] Crete, Greece: Crete University Press. 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. Vlasov, V. Z., and N. N. Leontiev. 1966. Beams, plates and shells on elastic foundations. Washington, DC: Israel Program for Scientific Translations.

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *