Buckling Analysis of Piles in Multi-Layered Soils

: Pile buckling is infrequent, but sometimes it can occur in slender piles (i.e., micropiles) driven into soils with soft layers and/or voids. Buckling analysis of piles becomes more complex if the pile is surrounded by multi-layered soil. In this case, the well-known Timoshenko’s solution for pile buckling is of no use because it refers to single-layered soils. A variational approach for buckling analysis of piles in multi-layered soils is herein proposed. The proposed method allows for the estimation of the critical buckling load of piles in any multi-layered soil and for any boundary condition, provided that the distribution of the soil coefﬁcient of the subgrade reaction is available. An eigenvalue-eigenvector problem is deﬁned, where each eigenvector is the set of coefﬁcients of a Fourier series describing the second-order displaced shape of the pile, and the related buckling load is the eigenvalue, thus obtaining the effective buckling load as the minimum eigenvalue. Besides the pile deformed shape, the stiffness distribution in the multi-layered soil is also described through a Fourier series. The Rayleigh–Ritz direct method is used to identify the Fourier development coefﬁcients describing the pile deformation. For validation, buckling analysis results were compared with those obtained from an experimental test and a ﬁnite element analysis available in the literature, which conﬁrmed this method’s reliability.


Introduction
The design problem of pile buckling is not usually faced in civil engineering practice because its occurrence is infrequent. However, this problem is becoming more important with the increasing use of micropiles, especially for retrofit interventions, and indeed elastic buckling of micropiles surrounded by soft soils is possible and is therefore worth investigating.
In the literature, different approaches have been adopted for studying pile buckling. Following Madhav and Davis' [1] study on buckling of piles in an elastic continuum and supported at the bottom by a rigid base, some other authors studied rod buckling through modelling the soil as an elastic continuum. Poulos and Davis [2] demonstrated that the soil models using Winkler's medium or elastic continuum show identical behaviour. Timoshenko and Gere [3] proposed a variational approach to study the elastic buckling of a bar on an elastic foundation, hinged at its ends and axially loaded on top. Timoshenko's study's conclusion is similar to the Engesser formula for calculating the buckling load of an elastically embedded beam.
Gabr et al. [4] developed a pile buckling model assuming a general power distribution of the soil's horizontal subgrade reaction to represent various soil conditions. The minimum potential energy method was adopted to develop the model using the Rayleigh-Ritz method to select a suitable deflection function. The effects of nine different boundary conditions on the buckling capacity and the equivalent buckling length as well as the effects of the distribution of the horizontal subgrade reaction were examined; the effect of the subgrade reaction was considered as fully embedded or partially embedded, so could only be used to evaluate two layers of soil [4]. Shields [5] compared pile buckling loads obtained from a semiempirical relationship developed by Bergfelt [6] with the allowable loads permitted by building codes and design guidelines [7]. A new formula (P cr = β (c u EI) 0.5 , where c u is undrained shear strength of the soil, EI is flexural rigidity, and β is a factor for regulation) derived from Engesser was derived to calculate the critical global buckling load. While the formula factor β was improved to calculate the pile imperfection, the method was for homogeneous soils [5]. Using the linear small-angle bending theory, Hegazy [8] provided a theoretical relationship between critical buckling load ratio and lateral deflection of a micropile subjected to vertical axial load and embedded in a weak homogeneous soil. Although the study covered different boundary conditions, it merely studied fully embedded piles.
Ofner et al. [9] extended buckling analysis of micropiles in homogeneous soils to micropiles driven into different soil layers. They based their assumptions on experimental buckling tests carried out on 4 m long micropiles embedded in soft clay. The research studied three types of pile-buckling calculations: 1. Code design calculations [10,11], 2. Experimental tests, and 3. FEM calculations using ABAQUS software. Despite the fact that the study considered multi-layered soil, the calculations were done for each layer separately, and the softest soil layer was the only candidate for the buckling position. Vogt et al. [12] highlighted how buckling failure can occur at higher undrained shear strengths in slender piles (i.e., micropiles) despite the codes' assumptions. In this case, the eccentricity of the pile was included in the momentum equilibrium of the forces under vertical load. Thus, this method considered pile imperfection in the derivations. Nevertheless, this method is only valid for single-layered soil.
Piles instability in liquefied soils is due to the absence of lateral support [13]. Buckling failure of piles, considering vertical and seismic loads, require a more critical regard of the undefined effects of forces [14]. In the recent decade, many research groups have studied liquefaction effects on buckling of piles. A method to calculate buckling critical load based on the beam-on nonlinear Winkler foundation (BNWF) model was proposed, considering multi-boundary conditions, the degree of soil layers liquefaction, and piles' nonlinearities. The critical buckling load of a pile in liquefied soils increases with the increase of the soil's relative density and the flexural rigidity of the pile [15,16]. Bhattacharya et al. [17] investigated the dynamic instability of piles in liquefiable soils. A correction factor was suggested to compensate for the natural frequency reduction of piles in a liquefied condition. Utilising the differential transformation method, Vega-Posada et al. [18] showed the effects of high and low soil stiffness and intermediate boundary conditions on pile buckling behaviour.
To reinforce the foundation of existing structures, micropiles require less effort to apply than do other driven piles, which are longer, and if they are precast, splicing is inevitable according to the site limit space [19]. Furthermore, micropiles enhance the load-bearing capacity of the foundations, even when soil is prone to liquefaction [20,21]. In this paper, the buckling of piles driven into multi-layer (non-homogeneous) soils for all boundary conditions is investigated. The uneven distribution of the coefficient of the subgrade reaction of a multi-layered soil is described by using a Fourier series. The boundary condition at the pile top can be given by assigning a suitable coefficient of the subgrade reaction to the part of the foundation crossed by the upper part of the pile. This means that the pile can be clamped or partially clamped at the top. The boundary condition at the bottom end is automatically assigned as it depends on the stiffness of the soil around the lower region of the pile. Using the principle of minimum potential energy, a variational approach to find buckling load through the Rayleigh-Ritz method has been proposed (Section 2.1). Two Fourier series were used, one describing the distribution of the coefficient of subgrade reaction and the other describing the displaced shape of the buckled pile. The algorithm was shown to be of general validity for any kind of multi-layer soil, and has been applied to different multi-layered soils, thus allowing for the investigation of the influence of position and thickness of layers of soft soils on buckling load (Section 3.1). The influence of the rotation of the pile sections in the closeness of the interface between layers with different stiffness, as well as of the rotation of the pile's lower end depending on the pile embedding in a stiffer soil layer close to its base were also considered (Section 3.1). This study shows how pile buckling can be an issue when slender piles pass through cavities or very soft soils, especially peaty soils, even if some silty sands and normally consolidated clays can be soft enough to allow for buckling of the slender piles driven in them (Section 3.1).

Methodology
Buckling analysis is often neglected in the design of foundation piles, but in the case of highly slender piles in soft soils, the pile load-bearing capacity can be reduced by buckling. Therefore, in this case, an elastic buckling analysis is equipped to calculate the critical buckling load of such slender compression members. Typically, buckling analysis involves solving an eigenvalue-eigenvector problem; that process was used in this article to perform buckling analysis of piles in multi-layered soils. In general, such a problem is faced using finite differences or finite elements, where one could even consider different soil stiffness values for each pile element. In contrast, by using a Fourier development of soil stiffness along the pile, this method considers continuous functions only, notwithstanding the stiffness discontinuity between different soil layers. Moreover, the variational approach was pursued using a direct method (the Rayleigh-Ritz method) to implement the analytical model. In the following, considering the applied load P as a function of the coefficients ak (k = 1, 2, . . . , m) of another Fourier series describing the pile modal shapes, the critical buckling load can be identified through minimising function P = P(a1, a2, . . . , ak, . . . , am) of m unknowns ak. Figure 1 shows the flowchart of the methodology adopted for buckling analysis of piles in multi-layered soils.

Formulation of the Problem
Consider a circular pile with a constant cross-section A, moment of inertia I of the cross-section, length L, and elastic modulus E. Modelling the soil as a Winkler's medium, the subgrade reaction K h can be assumed as either constant or linearly variable, depending on the soil type. The same coefficients of the subgrade reaction as those used in laterally loaded piles are considered [22,23]. A list of coefficients of subgrade reaction for different soils is provided by [2,24,25].
For a pile with a constant cross-section and, therefore, constant diameter d, considering an abscissa x with origin at the pile top, the soil reaction is The pile is axially loaded, hinged at its ends, and embedded in a Winkler's soil with varying rigidity K = K(x), ( Figure 2). Function K(x) is assumed to be a step function. Therefore, depending on step height and spacing, K(x) can vary almost linearly with K h (x) (depending on soil type) or be constant.
Through using the energetic criterium, the potential energy U 1 in the pile added to the potential energy in the Winkler's soil U 2 is equal to the decrease of potential energy U 3 of the external force P applied at the pile top, that is: Say y is the second-order displacement of the pile subject to the axial force P. In the following, the order of the derivatives of y are indicated through a Roman numeral superscript. The potential energy of the pile and that in the Winkler's soil are, respectively, expressed as: where E and I are the elastic modulus and moment of inertia of the pile cross-section, respectively, and y is the second-order displacement.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 18 Figure 1. Flowchart of the proposed method for buckling analysis of piles in multi-layered soils.

Formulation of the Problem
Consider a circular pile with a constant cross-section A, moment of inertia I of the cross-section, length L, and elastic modulus E. Modelling the soil as a Winkler's medium the subgrade reaction Kh can be assumed as either constant or linearly variable, depending on the soil type. The same coefficients of the subgrade reaction as those used in laterally loaded piles are considered [22,23]. A list of coefficients of subgrade reaction for different soils is provided by [2,24,25].
For a pile with a constant cross-section and, therefore, constant diameter d, consid-  Function K(x) is assumed to be a step function. Therefore, depending on step height and spacing, K(x) can vary almost linearly with Kh(x) (depending on soil type) or be constant.
Through using the energetic criterium, the potential energy U1 in the pile added to the potential energy in the Winkler's soil U2 is equal to the decrease of potential energy U3 of the external force P applied at the pile top, that is: Say y is the second-order displacement of the pile subject to the axial force P. In the following, the order of the derivatives of y are indicated through a Roman numeral superscript. The potential energy of the pile and that in the Winkler's soil are, respectively, expressed as: where E and I are the elastic modulus and moment of inertia of the pile cross-section, respectively, and y is the second-order displacement. The variation of the potential energy of the force P is: By substituting Equations (3)-(5) into Equation (2), we obtain: where K(x) is the soil horizontal stiffness along the pile. Therefore, the buckling load P is obtained from Equation (6) by finding the function y = y(x) for which the functional P of Equation (6) is a minimum, that is Pcr = min(P). For this aim, as y and K are variables dependent on the pile's depth, the Fourier series of y The variation of the potential energy of the force P is: By substituting Equations (3)-(5) into Equation (2), we obtain: where K(x) is the soil horizontal stiffness along the pile. Therefore, the buckling load P is obtained from Equation (6) by finding the function y = y(x) for which the functional P of Equation (6) is a minimum, that is P cr = min(P). For this aim, as y and K are variables dependent on the pile's depth, the Fourier series of y and K are developed, considering y and K as odd functions with period 2L in the interval (−L, L), thus obtaining: Coefficients a k , k = 1,2, . . . m, are unknown and will be obtained through solving the variational problem of calculating the buckling load. Coefficients b i , i = 1, 2, . . . , n, in Equation (8) depend on function K(x) distribution along with the pile. A multi-layered soil with s layers, each one with constant rigidity, is herein considered. This leads to: The generic coefficient b i is obtained within this equation, which is transformed to: where x r is the initial abscissa of the r-th layer, r = 1, 2, . . . , s. After substituting Equations (7) and (8) into Equations (3)-(5), the potential energy of the pile, soil, and external force are expressed as: After integration, U 1 and U 3 become, respectively: Taking into account the Fourier series of y(x) and K(x), it can be noted that the generic term of the integrating function of the expression of U 2 in Equation (12) where i = 1, 2, . . . , n, indicates the terms of the Fourier series expansion of soil rigidity K(x), and j, k = 1, 2, . . . , m, indicates the term of the Fourier series expansion of the pile deformed shape y(x). From the well-known trigonometric formulae, we have: Hence, for S 1 = i + j − k, S 2 = j + k − i, S 3 = k + i − j, S 4 = i + j + k, the expression of U 2 becomes: By integrating: and substituting: where Φ ijk = Φ ikj , the soil potential energy is thus obtained as: Finally, after substituting Equations (14), (15), and (21) into Equation (2), and extracting P, we obtain:

Conditions for Minimum P
Equation (22) synthetically shows that, for a given geometry and various mechanical properties of the pile driven into the multi-layered soil, the load P can vary through varying the coefficients a 1, a 2 , . . . , a k , . . . , a m that define the Fourier series of the deformed shape of the pile. To find the buckling load P cr and the values of the coefficients a k defining the related pile deformed shape, the Rayleigh-Ritz method is used. It is a so-called direct variational method that allows for finding of a solution for the minimum of the functional P[y(x), y"(x), y (x), x] defined in Equation (6) and re-arranged in Equation (22) after developing through a Fourier series the functions y(x) and K(x). Since for Equation (7), the function y(x) is a linear combination of m sinusoidal functions through the m unknown constants a k , k = 1, . . . , m, then, according to the Rayleigh-Ritz method, the functional P[y(x), y" (x), y (x), x] can be considered as a function P = P(a 1, a 2 , . . . , a k , . . . , a m ) of m unknowns a k . The stationary condition for this function is then imposed through the stationary conditions for P = P(a 1, a 2 , . . . , a k , . . . , a m ), that is: thus allowing us to find the m coefficients a k and therefore the related buckling deformation shape y(x). Of course, the higher the number of equations and therefore of coefficients a k used to describe the function y = y(x), the better is the accuracy of the determination of y(x) and therefore of the buckling load P cr . Since in Equation (6) After deriving Equation (24), we can obtain: By varying k, this is an equation system that defines an eigenvalue-eigenvector problem. It can be, in fact, noted that Equation (25) defines a characteristic system of equations, that is: where I is the unit matrix, and the square matrix A is the sum of a matrix A 1 is: and of a matrix A 2 with all terms zero except the m terms on the diagonal is: Therefore, from the characteristic Equation (26), m eigenvalues P cr,k , k = 1, . . . , m, can be obtained, and the minimum one is the effective buckling load P cr . For each eigenvalue P cr,k , an eigenvector a k = [a 1, a 2 , . . . , a k , . . . , a m ] T can also be obtained. Saying the eigenvector is related to the effective buckling load P cr , the related eigenvector can be obtained through solving the system: [A − P cr I] a = 0 (29) a = [a 1, a 2 , . . . , a k , . . . , a m ] T For P = P cr the deformed shape y(x) of the pile, as defined through its Fourier series in Equation (7), is thus obtained.

Discussion and Results
After having defined the methodology, the accuracy and efficiency of the calculations was to be adjusted. Since P cr depends on the number m of elements of the related eigenvector (where it is on the number m of coefficients a k ) and converges to an asymptotic value for increasing values of m, it is necessary to control the solution accuracy by checking if two successive values of P cr are sufficiently close for two respective values of m. Figure 3 shows an example of how the buckling load P cr of a pile with flexural stiffness EI = 100 kNm 2 converges to a sufficiently approximated value for an increasing number of equations. It is well evident that the larger is the number of equations, the closer is the convergence to the actual critical buckling load. For example, a sufficiently approximated buckling load with optimal deformed shape is reached with 25 coefficients a k , but 20 equations already provide essentially the same numerical result.
A similar consideration can be made regarding the number of coefficients of the Fourier series of the step function describing the distribution of the soil stiffness K(x). Figure 4 shows the distribution with the depth of function K(x) for soil with four layers in which K(x) is described through a Fourier series with 100 terms. The results have shown that even by halving the number of terms describing K(x), the value of P cr varies only slightly, although the given step function results are clearly less described by a Fourier series with a lower number of terms.

Applications and Validations
In general, elastic buckling may occur in foundation piles with lower flexural rigidity, as in micropiles for foundation underpinning, especially if surrounded by soft soils or crossing caves and voids.
Next, some applications are described. First, a comparison of the critical buckling load of micropiles in a single-layered (homogeneous) soil calculated both through the variational approach herein presented and through the Engesser [26] approach was made ( Figure 5). The outcomes showed identical results for any soil stiffnesses.
However, since this method can also be suitably applied to piles driven in more than one soil layer, next, some examples of pile buckling in multi-layered soils were analysed ( Figure 6). In addition, the outcome of the experimental test of a pile driven into multilayered soil [9] was compared to validate the proposed method. A similar consideration can be made regarding the number of coefficients of the Fourier series of the step function describing the distribution of the soil stiffness K(x) Figure 4 shows the distribution with the depth of function K(x) for soil with four layers in which K(x) is described through a Fourier series with 100 terms. The results have shown that even by halving the number of terms describing K(x), the value of Pcr varies only slightly, although the given step function results are clearly less described by a Fourie series with a lower number of terms.

Applications and Validations
In general, elastic buckling may occur in foundation piles with lower flexural rigid ity, as in micropiles for foundation underpinning, especially if surrounded by soft soils o crossing caves and voids.
Next, some applications are described. First, a comparison of the critical buckling load of micropiles in a single-layered (homogeneous) soil calculated both through the variational approach herein presented and through the Engesser [26] approach was made ( Figure 5). The outcomes showed identical results for any soil stiffnesses.  A similar consideration can be made regarding the number of coefficients of the Fourier series of the step function describing the distribution of the soil stiffness K(x). Figure 4 shows the distribution with the depth of function K(x) for soil with four layers in which K(x) is described through a Fourier series with 100 terms. The results have shown that even by halving the number of terms describing K(x), the value of Pcr varies only slightly, although the given step function results are clearly less described by a Fourier series with a lower number of terms.

Applications and Validations
In general, elastic buckling may occur in foundation piles with lower flexural rigidity, as in micropiles for foundation underpinning, especially if surrounded by soft soils or crossing caves and voids.
Next, some applications are described. First, a comparison of the critical buckling load of micropiles in a single-layered (homogeneous) soil calculated both through the variational approach herein presented and through the Engesser [26] approach was made ( Figure 5). The outcomes showed identical results for any soil stiffnesses. The critical buckling load of six micropiles made of steel hollow bars, length 12 m, embedded in different multi-layered soils were analysed. The geometrical and mechanical characteristics of the pile hollow bars are briefly described in Table 1 [27]. The micropile tube is usually surrounded by mortar. However, in this study, the flexural stiffness of the mortar was neglected due to the uncertain contribution in flexural resistance because of cracks forming in the mortar. Figure 6 illustrates the soil stratigraphy and the pile deformed shapes.
The load-bearing capacity of the micropiles under consideration would be N pl = A s · f y = 676 kN, unless buckling occurs before N pl is attained, depending on the mechanical characteristics of the different soil layers and on soil stratigraphy, in particular, if soft layers and caves are present. This is what happens in the following examples. The soil stiffness was calculated with two different formulas based on the soil type. For cohesive soils, K h = 60 ·c u /d, where K h is the subgrade reaction, and c u is the undrained shear strength of the soil; for cohesionless soils, K h = n h ·z/d, where n h is the coefficient of the subgrade reaction, z is the soil depth, and d is the diameter of the pile.
Note that when the stratigraphy includes soft soils and rock layers, high precision in the evaluation of rock rigidity is not required because the buckling load is mainly influenced by the presence of soft layers rather than by the variation of rock rigidity, that in any case is very high compared to the rigidity of soft soil. Moreover, the stiffness of the reinforced concrete foundation at the pile top end is given to be of the same magnitude as that assigned to the rock layers.  However, since this method can also be suitably applied to piles driven in more than one soil layer, next, some examples of pile buckling in multi-layered soils were analysed ( Figure 6). In addition, the outcome of the experimental test of a pile driven into multi-layered soil [9] was compared to validate the proposed method. The buckling load of the pile in Figure 6a is mainly affected by the not-embedded part of the pile between the foundation and the first soil layer. The pile also deviates from the straight configuration when the lower cave is passed through, but this deviation is so slight that it can be identified only by calculations, while it is not distinguishable in the plot. Figure 6b shows that when a pile, instead of being embedded in the limestone rock, is only embedded in the soft clay layer at the cave base, its bottom end is free to rotate, thus causing the buckling load to decrease dramatically and become much less than N pl . Figure 6c,d show that if a pile passes through a limestone cave with two layers of soft soil with low but different stiffnesses (peat and very soft clay) sedimented over the bottom of the cave, buckling load is highly affected by the inversion of the sedimentation order of the two soft soils. Figure 6e,f show that a micropile may buckle in soft soil layers with c u even higher than 15 kPa, which is more than the limit proposed by some design codes for checking pile buckling [11,12]. Even if there are no cavities or peaty soils where buckling can more easily occur, the low stiffness of these types of soils can make a pile prone to buckle. In fact, Figure 6f shows that buckling can also occur in soft soils without peat, such as very soft clays and silty sands. Buckling load results increased with the absence of peat, but were still lower than N pl . The buckling behaviour of one more pile was analysed to validate the proposed method through comparison with the experimental results obtained by testing a real micropile in situ [9]. The micropile characteristics are presented in Table 2 [27], and the soil layers specifications are described in Figure 7. The pile load-bearing capacity is N pl = A s ·f y = 2750 kN. The flexural stiffness of the grout was neglected, as in the previous buckling analyses. Loads were applied in successive steps of ∆F = 325 kN over a period of about five hours. The experimental result shows that the micropile experienced buckling between the fifth and sixth load steps, which is between 1625 to 1923 kN. In addition, the authors implemented in ABAQUS a finite element model (FEM) for which buckling occurred with P cr = 2000 kN, which is decreased to 1740 kN accounting for the eventual effect of imperfections. Through the variational approach proposed in this paper, P cr = 1923 kN was obtained. Moreover, the same pile (details in Table 2) was modelled in SAP2000 [28] to confirm the approach's accuracy. The buckling factor was obtained from the modal analysis, using horizontal springs on the pile representing Winkler's medium. For the sake of precision, spring spacing was 25 cm, and spring stiffness was selected according to the coefficient of the subgrade reaction in this interval. The output of the analysis showed a buckling factor equal to 1895 kN with the identical deformed shape. It must be noticed that while the variational method was a continuum-based analytical method, the FEA was performed on a discrete-spring model. Even if the experimental evaluation of the buckling load cannot be precise, the results obtained through the proposed method were shown to fit very well with those obtained from finite element analyses and quite well with the experimental results. Yield strength (f y ) MPa 500 Cross-section moment of inertia (I) mm 4 4,200,000 Therefore, pile buckling is significantly affected by the low rigidity of the softest layer, where second-order displacements can more easily occur between the two stiffer layers over and below the soft layer. In particular, pile buckling is considered critical in caves, both when the cave is between two stiff layers (i.e., made of rock) and when a soft soil layer fills up the bottom of the cave. The latter case is even more critical when the rotation of the bottom end of the pile is free, as it happens, for instance, if it is not embedded sufficiently at the bottom end of the cave, even if the cave is partially filled up by a soft soil layer.

Conclusions
As expected, the proposed approach was demonstrated to be of general use, because developing in Fourier series the stiffness distributions of the different layers of a nonhomogeneous soil allows for the handling of many different stiffnesses of the different layers as a single non-uniform stiffness distribution function. Additionally, this feature can be used to determine boundary conditions. In this way, pinned, clamped, or partially clamped ends of piles can be assigned by varying the stiffness and thickness of the soil at that level.
Furthermore, the Rayleigh-Ritz method successfully solved complex buckling analysis problems, where a pile usually crosses several soil layers. Indeed, adopting the Rayleigh-Ritz method effectively simplifies the variational problem of minimum to a minimisation problem of the function of the coefficients of the Fourier development describing the second-order displaced shape y(x) of the pile subject to buckling.
The robustness and efficiency of the proposed method have allowed for the solving of many different cases of buckling analyses of piles in multi-layered soils. In order to understand how different soil layers with varied depths and stiffnesses can affect the critical buckling load, some examples were reported to highlight some crucial cases. It has been confirmed that pile buckling is unlikely to occur unless some voids or very soft layers are interposed between stiffer layers. Unfavourable boundary conditions at the pile top (not enough of it clamped to the foundation) and at the pile bottom (insufficient embedding of the pile's lower end in a stiff layer) have the effect of decreasing buckling load.
Since the proposed method is valid for vertical piles with constant cross-sections, in future research the cases of tapered piles, inclined piles, as well as the effect of imperfections will be investigated.

Conflicts of Interest:
The authors declare no conflict of interest.