Energy-Based Approach: Analysis of a Vertically Loaded Pile in Multi-Layered Non-Linear Soil Strata

Abstract

Numerous studies have been reported in the published literature on analytical solutions for a vertically loaded pile installed in a homogeneous single soil layer. However, piles are rarely installed in an ideal homogeneous single soil layer. This study presents an analytical model based on the energy-based approach to obtain displacements in an axially loaded pile embedded in multi-layered soil considering soil non-linearity. The developed analytical model incorporating Euler-Bernoulli beam theory proved to be an effective way in estimating the load-displacement responses of piles embedded in multi-layered non-linear elastic soil strata. The differential equations are solved analytically and numerically using the variational principle of mechanics. A parametric study investigated the effect of explicit incorporation of soil properties and layering in order to understand the importance of predicting appropriate pile displacement responses in linear elastic soil system. It is clear from the results that the analyses which consider the soil as a single homogeneous layer will not be able to produce an accurate estimation of the pile stiffnesses. Therefore, it is highly important to account for the effect of soil layering and the non-linear response. The pile displacement response is obtained using the software MATLAB R2019a and the results from the energy-based method are compared with those obtained from the field test data as well as the Finite Element Analysis (FEA) based on the software ANSYS 2019R3. The non-linear elastic constitutive relationship which described the variation of secant shear modulus with strain through a power law has shown reasonably accurate predictions when compared to the published field test data and the FEA. The developed mathematical framework is also more computationally efficient than the three-dimensional (3D) FEA.

1. Introduction

Many studies are available in the literature that explains the pile deformation response under different loading conditions with the pile embedded in multi-layered soil, which treat the soil as linear elastic. However, for practical applications in the real world, it is important to study the response of the non-linear behavior of the soil and the variation of soil stiffness based on its strain levels. A discretized continuum approach (numerical method) or the energy-based method is preferred over the beam-on-elastic foundation approach since the approach considers the surrounding soil behavior in three dimensions: vertical, radial and circumferential directions. Several researchers estimated the pile deformations based on the energy-based method [1,2,3,4,5,6,7]. There are many numerical methods available in the published literature with various pile geometry and constitutive models using the Finite Element Method (FEM) [8,9,10,11,12,13], the Finite Difference Method (FDM) [8,14,15,16,17,18] and the Boundary Element Method (BEM) [3,4,19,20,21,22]. The present study is based on the FDM since it yields results faster when compared to the FEM.

1.1. Linear Elastic Soil Model

An analytical solution based on the energy-based method was used by several researchers [8,16,23,24,25,26,27,28] to estimate laterally loaded and axially loaded pile deformation for linear elastic soil. Independent functions describing the soil displacement have been used; these functions vary in vertical, radial and circumferential directions. The linear elastic analysis has been developed by employing variational principles and minimization of energy, called Hamilton’s principle, to derive the governing equation and boundary conditions. Hamilton’s equation can be expressed as

$${{\displaystyle \int }}_{t_1}^{t_2}\left(\delta T-\delta U\right)dt+{{\displaystyle \int }}_{t_1}^{t_2}\delta W dt$$

(1)

where T and U are the kinetic and potential energies of the pile and soil; W is the work done by the applied load; and t₁ and t₂ are the initial and final times of loading [29].

The governing equations for pile deflection are obtained by minimizing kinetic and potential energies. These governing equations can be solved either numerically or analytically for a given set of boundary conditions. Each of the displacement components are expressed as a multiplication of one-dimensional functions when minimizing the energy to obtain a set of one-dimensional equation. These equations are solved numerically using the finite difference technique. In these studies, the soil is treated as linear elastic where the shear modulus G and the Lame’s constant λ are treated as constants for all the layers [30].

1.2. Constitutive Model for Non-Linear Soil Behavior

In order to consider the non-linear effect of soil, the soil moduli G and λ are assumed to vary in radial, circumferential and vertical directions according to the strain levels [7]. The strain decay with increasing distances in the radial direction is accompanied by an increase in soil stiffness G. In other words, soil stiffness degrades with increasing strain and hence it varies in both radial and vertical directions. The analysis of the soil behavior using a single constitutive model is very idealistic since the undrained strength of soil depends on a number of factors such as the soil anisotropy, failure mode, strain rate, stress paths and the mode of loading effects of stress–strain non-linearity which make the undrained strength dependent on the test type [31,32]. Consideration of the non-linear soil behavior will be more realistic in the analysis of pile displacement response. For piles subjected to external loads, the decay of the soil stiffness varies with strain which in turn depends on the type of the soil. The stiffness of soil was high at a very small strain level and decreases with the increase in the strain [33]. Many researchers [34,35,36,37,38,39] conducted triaxial tests and reported high values of soil stiffness when the shear strains are less than 10⁻⁵. Several factors including the mean effective stress, void ratio, stress history, rate of loading, soil plasticity for silts and clays, stress anisotropy for sands and the effective confining stress affect the small strain stiffness G_max [40,41,42,43]. The decay of the soil stiffness with the increase of strain levels can be defined using the power law [44,45,46].

1.3. Research Significance

The non-linear responses of piles subjected to external loading conditions are analyzed using the 3D Finite Element (FE) software packages with inbuilt constitutive soil models based on the theory of elasticity and plasticity which are beneficial when determining the ultimate load capacity. However, such models might not be necessary for analyzing the pile head displacements as non-linear elastic models may be sufficient to predict considerably accurate and faster pile responses. Therefore, an analytical model has been developed based on the energy-based approach to predict the load-displacement responses of the piles embedded in multi-layered soil strata subjected to static axial loads. The non-linear elastic constitutive relationship which described the variation of secant shear modulus with strain through a power law (proposed by Gunn [46]) has shown reasonably accurate predictions when compared to the published field test data and the FEA.

2. Problem Definition

The present study is based on an energy-based approach (simplified continuum model) to study the displacement profile of an axially loaded pile in multi-layered soil strata considering the soil’s non-linear elastic behavior. A mathematical framework is developed that considers the soil as a 3D elastic continuum and the pile is modeled following the elastic Euler-Bernoulli beam theory. The differential equations governing the soil and pile displacements have been obtained utilizing the principle of minimum potential energy assuming rational soil displacement functions. An iterative algorithm is adopted that solves the differential equations analytically and numerically. The study investigates the effect of explicit incorporation of soil properties and layering in order to understand the importance of predicting appropriate pile displacement responses in linear elastic soil system. The responses indicated that the effect of soil layers and their thicknesses, pile properties and the variation in soil moduli have a direct impact on the displacements of piles subjected to axial loading. Hence, proper emphasis has to be given to account for the soil non-linearity.

An axially loaded pile in an isotropic non-linear elastic multi-layered soil medium is shown in Figure 1. This study considers a pile of length L with circular cross section of radius r₀. The pile is embedded in n horizontal soil layers and is subjected to an axial (vertical) load P_t. The horizontal soil layers extend to infinity in the radial direction and the bottom n^th layer extends to infinity in the vertical direction.

The terms H₁, H₂, H₃…H_n−1 denote the vertical height from the ground surface to the bottom of any layer i. Therefore, the thickness of any layer i is H_i–H_i−1 with H₀ = 0. Due to the axisymmetric problem behavior, a system of cylindrical coordinates (r-θ-z) is chosen with the origin coinciding with the center of the pile cross section at the pile head and the z axis coinciding with the pile axis. The pile head is considered to be free and the tip of the pile is clamped. Another important assumption to be noted is that there is no slippage or separation between the pile and the surrounding soil and between soil layers. The stresses and the displacement within a soil continuum are shown below in Figure 2a,b.

At each point in the soil domain, G and λ have been calculated (Figure 3).

The goal of the analysis is to obtain pile head deflection caused by the action of the vertical load at the pile head.

2.1. Basic Assumptions

For an axially loaded pile, the horizontal and tangential displacements can be neglected as they are accompanied by very small strains [47]. For the case of a pile with circular cross-section, there are two functions to be considered: v(z), which will represent the vertical displacement at depth z and the dimensionless functions ϕ(r), describing the variation of soil displacements in the radial direction.

The vertical displacement at any point of the soil is represented as a function in (r, z):

$$v_{\mathrm{r}}=0$$

(2)

$$v_{\theta }=0$$

(3)

$$v_{\mathrm{z}}\left(\mathrm{r}, \mathrm{z}\right)=\varphi \left(\mathrm{r}\right) v\left(\mathrm{z}\right)$$

(4)

For a given uniform cross-sectional area of the pile along the length, ϕ(r) = 1 when r = 0 to r0, whereas ϕ(r) = 0 when r→∞. This explains the decay of the function ϕ(r) with an increase in the radial direction.

2.2. Governing Differential Equation

The pile and its surrounding elastic medium are subjected to a vertical displacement of the pile soil system when it is acted upon by a vertical load. The total potential energy of the pile and the soil is a summation of internal potential energy and external potential energy [8], which is given as:

$$\Pi =\frac{1}{2} E_p{A}_p{{\displaystyle \int }}_0^L{\left(\varphi \frac{dv}{dz}\right)}^2+\frac{1}{2} {{\displaystyle \int }}_0^L{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_{r_0}^{\infty }{\sigma }_{ij}{\epsilon }_{ij}rdrd\theta +\frac{1}{2} {{\displaystyle \int }}_L^{\infty }{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^{\infty }{\sigma }_{ij}{\epsilon }_{ij}rdr-v{P}_{t at z=0}$$

(5)

where E_p denotes the elastic Young modulus of the pile, A_p denotes the cross-section of the pile, v represents the vertical displacement, P_t is the vertical load and σ_ij, ε_ij are stress and strain components, respectively. The first term of the equation represents potential pile energy, the second and third terms are potential energy from the surrounding soil and the soil below the pile, respectively.

The soil non-linearity is considered by varying the soil elastic parameters (G and λ) at each discretized nodal point in the soil domain. The stress–strain and strain–displacement relationships at any given nodal point in the soil medium are idealized by the following relationships. The stress–strain relationship is expressed as:

$$\left[ \begin{array}{c}{\sigma }_r\\ {\sigma }_{\theta }\\ {\sigma }_z\\ {\tau }_{r\theta }\\ {\tau }_{rz}\\ {\tau }_{\theta z}\end{array}\right]=\left[ \begin{array}{cccccc}\lambda +2G& \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2G& \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2G& 0& 0& 0\\ 0& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G\end{array}\right] \left[ \begin{array}{c}{\epsilon }_r\\ {\epsilon }_{\theta }\\ {\epsilon }_z\\ {\gamma }_{r\theta }\\ {\gamma }_{rz}\\ {\gamma }_{\theta z}\end{array}\right].$$

(6)

where G and λ are the elastic constants of the soil. The strain–displacement relationship is given by:

$$\left[ \begin{array}{c}{\epsilon }_r\\ {\epsilon }_{\theta }\\ {\epsilon }_z\\ {\gamma }_{r\theta }\\ {\gamma }_{rz}\\ {\gamma }_{\theta z}\end{array}\right]= \left[ \begin{array}{c}-\frac{\partial u_r}{ \partial r}\\ -\frac{u_r}{r}-\frac{1}{r}\frac{\partial u_{\theta }}{ \partial \theta }\\ - \frac{\partial u_z}{ \partial z}\\ -\frac{1}{r}\frac{\partial u_r}{ \partial \theta }-\frac{\partial u_{\theta }}{ \partial r}+\frac{u_{\theta }}{r}\\ -\frac{\partial u_z}{ \partial r}-\frac{\partial u_r}{ \partial z}\\ -\frac{1}{r}\frac{\partial u_z}{ \partial \theta }-\frac{\partial u_{\theta }}{ \partial z}\end{array}\right]= [\begin{array}{c}0\\ 0\\ -\varphi \left(r\right)\frac{dv\left(z\right)}{dz}\\ 0\\ -v\left(z\right)\frac{d\varphi \left(r\right)}{dr}\\ 0\end{array}]$$

(7)

By substituting Equation (7) into Equation (6), the strain energy density function $W=\frac{{\sigma }_{ij}{\epsilon }_{ij}}{2}$ is obtained, where the summation implies the repetition of the indices i and j as required in indicial notation:

$$\frac{1}{2} {\sigma }_{ij} {\epsilon }_{ij}=\frac{1}{2} \left[\left(\lambda +2G\right){\left(\varphi \frac{dv}{dz}\right)}^2+G{\left(v\frac{d\varphi }{dr}\right)}^2\right]$$

(8)

Substituting Equation (8) into Equation (5) and integrating with respect to θ, the potential energy equation becomes:

$$\begin{array}{ll}\Pi =\frac{1}{2} E_p{A}_p{{\displaystyle \int }}_0^L& {\left(\varphi \frac{dv}{dz}\right)}^{2}dz\\ & +\pi {{\displaystyle \int }}_0^L{{\displaystyle \int }}_{r_0}^{\infty }\left(\left(\lambda +2G\right){\left(\varphi \frac{dv}{dz}\right)}^2+G{\left(v\frac{d\varphi }{dr}\right)}^2\right)rdrdz\\ & +\pi {{\displaystyle \int }}_L^{\infty }{{\displaystyle \int }}_0^{\infty }\left(\left(\lambda +2G\right){\left(\varphi \frac{dv}{dz}\right)}^2+G{\left(v\frac{d\varphi }{dr}\right)}^2\right)rdrdz-v{P}_{t at z=0}\end{array}$$

(9)

The variational principle has been used to calculate the potential energy δU and the external energy δW [48,49]. As a result, the governing equations of the pile-soil system are obtained by minimizing the potential energy of soil and pile. The expression of potential energy contains different functions, such as v(z), ϕ(r), dv(z)/dz and dϕ(r)/dr, and so minimizing the potential energy gives:

$$\delta \Pi =\left[A\left(v\right)\delta v+B\left(v\right)\delta \left(\frac{dv}{dz}\right)\right]+\left[C(\varphi )\delta \varphi \right]$$

(10)

where A, B and C are the terms associated with variations $\delta v$, $\delta$ (dv/dz) and $\delta \varphi$

The variation of Equation (9) becomes

$$\begin{array}{ll}\delta \Pi =\frac{1}{2}{E}_p{A}_p{{\displaystyle \int }}_0^L& \varphi \frac{dv}{dz}\delta \varphi \frac{dv}{dz}dz\\ & +\pi {{\displaystyle \int }}_0^L{{\displaystyle \int }}_{r_0}^{\infty }\left(\left(\lambda +2G\right)\left(\varphi \frac{dv}{dz}\delta \varphi \frac{dv}{dz}\right)+G\left(v\frac{d\varphi }{dr}\delta v\frac{d\varphi }{dr}\right)\right)rdrdz\\ & +\pi {{\displaystyle \int }}_L^{\infty }{{\displaystyle \int }}_0^{\infty }\left(\left(\lambda +2G\right)\left(\varphi \frac{dv}{dz}\delta \varphi \frac{dv}{dz}\right)+G\left(v\frac{d\varphi }{dr}\delta v\frac{d\varphi }{dr}\right)\right)rdrdz-v{P}_{t at z=0}\end{array}$$

(11)

2.3. Output Parameters

(i.)

Pile Displacement

The governing equation of the pile is obtained for 0 < z < L by collecting terms associated with $\delta$vdz, $\delta$(dv/dz) dz and its derivative, $\delta$v and $\delta$(dv/dz)dz ≠ 0. The governing equation is obtained as follows:

$$k\frac{d^{2}v}{d{z}^2}+C\frac{dv}{dz}+mv=0$$

(12)

where

$$C=2\pi {{\displaystyle \int }}_{r_0}^{\infty }r\left[\left(\lambda +2G\right){\varphi }^2\right]dr$$

(13)

$$k=E_p{I}_p+2\pi {{\displaystyle \int }}_{r_0}^{\infty }r\left[\left(\lambda +2G\right){\varphi }^2\right]dr$$

(14)

$$m=-2\pi {{\displaystyle \int }}_{r_0}^{\infty }rG{\left(\frac{d\varphi }{dr}\right)}^{2}dr$$

(15)

For this study, the tip of a pile is assumed to be clamped, which means that the displacement and the curvature are equal to zero at the base of the pile. The boundary conditions are obtained by collecting δv and $\delta$(dv/dz). At the head of the pile (z = 0):

$$-E_p{A}_p-2\pi {{\displaystyle \int }}_{r_0}^{\infty }\left[\lambda +G\left(2+{\varphi }^2\right)\right]dr\frac{dv}{dz}+P_t=0$$

(16)

The displacement at the tip of the pile (z = L):

$$v=0$$

(17)

The second order differential Equation (12) can be solved using a central finite difference scheme. Equation (12) becomes

$$k\left(\frac{v_{i-1}-2v_i+v_{i+1}}{\mathsf{\Delta }{z}^2}\right)+C\left(\frac{v_{i-1}+v_{i+1 }}{\mathsf{\Delta }z}\right)+m{v}_i=0$$

(18)

where i denotes the ith node in z direction, and Δz is the distance between two nodes. This discretized analysis is then solved using the software MATLAB R2019a.

(ii.)

Soil Displacement

The governing equation of the soil is obtained for r₀ ≤ r ≤ ∞ by collecting terms associated with $\delta$ϕdr:

$$r\frac{d^2\varphi }{d{r}^2}+{\gamma }_1\frac{d\varphi }{dr}+{\gamma }_2\varphi =0$$

(19)

where

$${\gamma }_1=\frac{{{\displaystyle \int }}_0^{L}G+rG}{{{\displaystyle \int }}_0^{L}G}$$

(20)

$${\gamma }_2=\frac{-{{\displaystyle \int }}_0^L\left(\lambda +rG\right){\left(\frac{dv}{dz}\right)}^2}{{{\displaystyle \int }}_0^{L}G{v}^2}$$

(21)

Note that the dimensionless parameter $\gamma$, defined in the above equation, describes the function $\varphi$. Similar to the solution of Equation (12) by the central finite difference scheme, the governing differential Equation (19) is also solved using the FDM (using the software MATLAB R2019a).

2.4. Soil Non-Linearity

The variation of the shear stress with strain can be described using two parameters A and n that have been obtained experimentally using a pressuremeter test as shown in the equation below:

$$q=A{\left({\epsilon }_q\right)}^n$$

(22)

where q represents the equivalent shear stress and ε_q is the deviator shear strain. Atkinson [33] shows that the decay curves of the soil stiffness with strain can be divided into three regions as shown in Figure 4. The first region in Figure 4 represents the very small strain where the stiffness, G₀ is constant, the second region comprising small strains starts from ε₀ until ε = 0.1% and the third region exceeds ε = 0.1%, which is indicative of large strains.

In the second region, the stiffness decays rapidly and in the third region with large strain levels, the stiffness is the smallest, which concludes that the soil stiffness is high at the small strain and decreases with the large strain [33]. Figure 5 shows the degradation of soil stiffness with increasing strains for different clay types.

The present study assumes decay of soil stiffness with strain using a power law to describe the stress–strain behavior of soil [44,46]:

$$G=G_0{\left(\frac{{\epsilon }_q}{{\epsilon }_{q0}}\right)}^n$$

(23)

$$G=a {\epsilon }_q^n$$

(24)

where $a=\frac{G_0}{{\epsilon }_{q0}^n}$ is a constant determined empirically; n describes soil nonlinearity which is equal to (−0.5) according to the experimental data analyzed by Osman et al. [39] (Figure 6); ${\epsilon }_q=\sqrt{\frac{2}{9}\left[{\left({\epsilon }_{rr}-{\epsilon }_{\theta \theta }\right)}^2+{\left({\epsilon }_{\theta \theta }-{\epsilon }_{zz}\right)}^2+{\left({\epsilon }_{zz}-{\epsilon }_{rr}\right)}^2\right]+\frac{4}{3}\left({\epsilon }_{r\theta }^2+{\epsilon }_{\theta z}^2+{\epsilon }_{zr}^2\right)}$ represents the deviatoric strain; and ε_q0 is the maximum deviatoric strain with linear elastic behavior which is equal to 10⁻⁵. The soil stiffness G is estimated by calculating the strain at each location followed by the power law.

3. Iterative Solution Methodology

The pile deflection equation can be solved when the soil and geometry related parameters k, C and m are known; however, these parameters depend on the unknown dimensionless soil function ϕ, which can be estimated by calculating ${\gamma }_1$ and ${\gamma }_2$. Soil displacement is obtained when the initial numbers of these values ${\gamma }_1$ and ${\gamma }_2$, are inserted into Equation (19), from which the parameters k, C and m are obtained as a result of the pile displacement. New values of ${\gamma }_1$ and ${\gamma }_2$ (Equations (20) and (21)) are determined and then inserted into Equation (19) to evaluate ϕ then v; therefore, an iteration technique is needed to obtain the condition $\frac{{\gamma }_{old}-{\gamma }_{new}}{{\gamma }_{old}}<0.001$. This iterative solution methodology (Figure 7) is used as input to the software MATLAB R2019a to obtain the pile head displacements.

4. Results

4.1. Effect of Explicit Incorporation of Soil Characteristics and Layering

A parametric study was performed in order to study the effect of explicit incorporation of soil characteristics and layering to predict appropriate pile displacement responses in linear elastic soil system using the proposed energy-based method. The study is performed in terms of normalized pile head stiffness, $\frac{P_t}{v_{z=0} E_{p}D}$ in which D (=2r₀) is the diameter of the pile. Evaluations were made for the cases of two and three-layered soil strata with different pile and soil parameters as shown in

Table 1. Pile and soil parameters for understanding the effect of explicit incorporation.

S. No.	Soil and Pile Parameters		Heights of Soil Layers, Hi (m)	Shear Moduli, G0 (MPa)
Case 1	(a)	L; D = 2r0; Ep; Pt Ep/G02 = 1000 L/D = 25	Two-layer soil system, H1 = 0 to L	G01/G02 = 0.2; G01/G02 = 0.5; G01/G02 = 1; G01/G02 = 2; G01/G02 = 5
	(b)	L; D = 2r0; Ep; Pt L/D = 25	Two-layer soil system, H1 = H2 = 0.5 L	G01/G02 = 0.2; G01/G02 = 0.5; G01/G02 = 1; G01/G02 = 2; G01/G02 = 5
Case 2	(a)	L; D = 2r0; Ep; Pt Ep/G = 1000 L/D = 25	Three-layer soil system, H1 = H2 = H3 = L/3	G01 = G; G02 = 2G01; G03 = 4G01
	(b)	L; D = 2r0; Ep; Pt Ep/G = 1000 L/D = 25	Three-layer soil system, H1 = H2 = H3 = L/3	G01 = 2G02; G02 = G; G03 = 4G02
	(c)	L; D = 2r0; Ep; Pt Ep/G = 1000 L/D = 25	Three-layer soil system, H1 = H2 = H3 = L/3	G01 = 4 G03; G02 = 2 G03; G03 = G

.

Two-layer soil system

The two-layered deposits included two different cases (Case (a) and Case (b)) with results obtained for five different soil modulus ratios: G₀₁/G₀₂ = 0.2, 0.5, 1, 2 and 5 and varying soil layer thicknesses (refer to Table 1). Figure 8a,b show the variation of the normalized pile head stiffness as a function of h/L (with and E_p/G₀₂ = 1000) and E_p/G₀₂ (with h = 0.5 L) when the pile slenderness ratio L/D = 25.

Three-layer soil system

The three-layer soil deposits included three different cases with the same thickness (L/3) but different equivalent shear modulus G, such that the average value G_avg = G₀₁ + G₀₂ + G₀₃/3 remains the same for all the cases. Case (a) has a variation of shear modulus increasing with depth with the smallest value observed for uppermost layers representing the case of end-bearing piles. Case (b) represents a profile with an intermediate weak layer. Case (c) has the shear modulus decreasing with depth with the largest value observed for the uppermost layer representing friction piles. Figure 9a shows the variation of the normalized pile head stiffness as a function of the pile slenderness ratio L/D with E_p/G = 1000 (for all the three cases) and Figure 9b shows the variation of the normalized pile head stiffness as a function of E_p/G with L/D = 25 (for all the three cases).

Results indicate that the normalized pile head stiffness decreases with the value of slenderness ratio L/D for end-bearing piles and increases for friction piles. The difference in the observed normalized pile stiffness becomes smaller as the soil becomes weaker (i.e., as E_p/G becomes larger) and when the soil becomes extremely weak, there is no difference observed in the normalized pile head stiffness for all the cases. It is clear from the results that the analyses which consider the soil as a single layer will not be able to produce an accurate estimation of the pile stiffnesses. Therefore, it is highly important to account for the effect of soil layering and the non-linear response.

4.2. Accuracy of the Model Considering the Soil Non-Linearity

(a).

Energy-based Method versus the field test data

“Full-scale load tests on instrumented micropiles”: Russo [51]

The present study makes a comparison of Russo’s field test data to validate the superiority of the energy-based method for non-linear elastic soil over the linear elastic solution of Salgado et al. [47]. Russo [51] outlined the case history of micropiles used for underpinning a historical building in Naples, Italy. A micropile is a type of drilled shaft with a diameter in the range of 75 to 300 mm [52]. A load test was carried out on a micropile, 0.2 m in diameter and 19 m in length with an elastic modulus of 27 GPa, as shown in Figure 10. The pile and soil characteristics used in the present study are adopted from the published literature by Seo and Prezzi [27].

The soil characteristics are given in

Table 2. Soil characteristics.

Soil Characteristics	Soil Thickness (m)	Young’s Modulus, E (MPa)	Shear Modulus, G0 (MPa)	Poisson’s Ratio
Layer 1	12	50	35	0.3
Layer 2	9	117	18	0.3

.

The power law equation (Equation (24)) is used to model the soil stiffness as a function of strain to account for soil non-linearity. In Equation (24), $a=\frac{G_0}{\epsilon q_0^n}$ is a constant value and the value of G₀ for each layer is taken as reported from Table 1, where $\epsilon$_qdonates deviatoric strain, $\epsilon$_q0 is the maximum deviatoric strain with linear elastic behavior which is equal to 10⁻⁵ and n represents a constant. The present study uses the published data for n = −0.5 proposed by Osman et al. [39].

A series of axial loads were applied to the pile head to obtain the corresponding vertical displacements. Figure 11 compares the pile displacements from the energy-based method with the field data and the linear elastic solution.

The predicted pile displacements in the present study are in good agreement with the field data. It can be seen that the pile displacements from the present study and the published field data are greater than the computed values based on the linear elastic soil solution when the axial load in the pile exceeds 300 kN.

“An instrumented driven pile in Dublin boulder clay”: Farrell et al. [53]

Farrell et al. [54] conducted field experiments with an instrumented tubular steel pile driven into Dublin boulder clay (DBC). A steel closed-end tubular driven pile, 0.273 m in diameter, 10 mm (0.01 m) in thickness and 7.5 m in length, was embedded in the boulder clay (Figure 12). A flat 20 mm (0.02 m) thick, 0.273 m diameter steel disc was used to close the pile base. The area’s geotechnical characteristics were investigated by various researchers [54,55,56,57]. A detailed description of the geology at the test site is given by Skipper et al. [54], whereas a brief summary is given by Long and Menkiti [55]. The geology of the area as investigated by (Long and Menkiti [55]) classifies the soil layers into four categories: Upper Brown Boulder Clay (UBrBC), Upper Black Dublin Boulder Clay (UBkBC), Lower Brown Dublin Boulder Clay (LBrBC) and Lower Black Dublin Boulder Clay. Long and Menkiti [55] reported G₀ = 98 MPa for the first layer and G₀ = 83 MPa for the second layer. A single ‘operational’ value of Eu = 100 MPa is used, as derived from field observations (Farrell et al. [58]).

The degradations of soil stiffness for UBrBC and UBkBC based on the present study are shown in Figure 13a,b, respectively. The initial soil stiffness used in the calculations is based on Long’s [55] study, and was also calculated using the power law relation in Equation (24). The present study uses the published data for n = −0.5 proposed by Osman et al. [39]. This value is in good agreement with published values for clay by Long and Menkiti [55].

A series of axial loads were applied to the pile in order to obtain the pile-soil deformation. The field test results are then compared with the analytical solutions from the energy-based approach. Figure 14 shows axial load versus observed and predicted pile head deformations.

Figure 14 shows that the energy-based method predicts the pile head displacements, which are in good agreement with the field measurements.

(b).

Energy-based Method versus FEA

A FEA was conducted using the software ANSYS Workbench 2019 R3 for the Soil-Pile problem as discussed in the above section [54]. A static structural analysis was chosen and the input data for the soil characteristics was taken from Long and Menkiti [55]. The pile characteristics were added according to the steel pile used by Farrell et al. [53]. The soil non-linearity was considered based on the actual soil stress–strain curves reported by Long and Menkiti [55] for each soil layer (as shown in Figure 15a,b).

The boundary conditions for the pile included a free head and fixed tip in order to replicate the assumption adopted in the proposed energy-based method. The soil block’s length was chosen as 15 times the pile diameter (15D) [30] and the depth was according to the depth of each soil layer. The soil block was fixed at the bottom and free on the surroundings. The interactions between the soil (contact element) and the pile (target element) have been defined manually to establish proper contact regions [59]. The discretized finite element model with a coarse mesh of 6978 nodes and 2039 elements is shown in Figure 16. A Directional Deformation along the Z-axis is used to calculate the displacements of pile head for the respective force applied.

Figure 17 compares the pile head displacements obtained by the energy-based method and the FEM Analysis. Both the methods underestimate the pile head displacements when compared to the field measured values. However, the displacements obtained by the energy-based method are in good agreement up to 20 kN. On the other hand, the FEA gives displacements close to field measured values only up to 10 kN and beyond, where the displacements appear to be smaller. Such an observation could be associated with poor quality meshes or improper surface interactions between the soil and pile model. More refined models may be necessary for future validations.

5. Discussions

Several researchers have studied the pile responses using numerical methods including BEM, FEM and FDM. Although the FEM using appropriate soil constitutive relationship, elements and domains for the soil and pile gives realistic results, the method requires enormous computation time and the resources required for such an analysis stands out as a major limitation. This gives rise to the need for a more efficient method which has both the strength of a rigorous three-dimensional non-linear approach for the pile-soil interaction and potentials for obtaining a faster computational effort. The present study presents an analytical solution based on an energy method (discretized continuum approach) to predict pile-soil displacements, where the soil is assumed to be non-linear elastic (soil parameters vary in radial, tangential and vertical directions). The analysis was conducted on an axially (vertical) loaded pile embedded in multi-layered soil. Governing equations for pile and soil have been obtained by applying the variational principle to the potential energy which are then solved using the software MATLAB R2019a. The soil nonlinearity is considered using non-linear constitutive relationships expressed in the form of a power law and a hyperbolic equation applicable to undrained clay deposits where the degradation of secant shear modulus is expressed as a function of the induced strain in soil. Comparisons have been made with the published field data and the FEM. It is observed that the energy-based method described in this study is in good agreement with the field data when compared to the linear elastic solution that does not consider the soil non-linearity. The FEA has been carried out using the software ANSYS Workbench 2019R3 in which the soil non-linearity is considered by inserting the stress–strain curves of each soil layer available from the published literature. It is observed that the results obtained from the energy-based method are in good agreement with the field measured values and those obtained from the FEA.

6. Conclusions

• An analytical model has been developed based on the energy-based approach to predict the load-displacement responses of the piles embedded in multi-layered soil strata subjected to static loading conditions.
• The analytical technique included the pile conventionally modeled as a Euler-Bernoulli beam and the soil as a 3D continuum considering the effect of non-linearity through constitutive relationships (describes the variation of secant modulus with strain).
• The differential equations are solved analytically and numerically using the variational principle of mechanics. A parametric study investigated the effect of explicit incorporation of soil properties and layering in order to understand the importance of predicting appropriate pile displacement responses in linear elastic soil system. Results indicate that the normalized pile head stiffness decreases with the value of slenderness ratio L/D for end-bearing piles and increases for friction piles. The difference in the observed normalized pile stiffness becomes smaller as the soil becomes weaker (i.e., as E_p/G becomes larger) and when the soil becomes extremely weak, there is no difference observed in the normalized pile head stiffness for all the cases. It is clear from the results that the analyses which consider the soil as a single layer will not be able to produce an accurate estimation of the pile stiffnesses. Therefore, it is highly important to account for the effect of soil layering and the non-linear response.
• The present energy-based method considering the non-linear response of the soil gives a good approximation of the field data when compared to the linear elastic solution and the FEA.
• The developed mathematical framework is also more computationally efficient than the 3D FEA.

7. Limitations and Recommendations for Future Work

• The main goal of this study is to extend the same approach to laterally loaded piles and to the combined action of lateral and axial loading on the pile.
• The non-linear analysis framework should be extended to include different constitutive soil models such as the elastic-plastic model especially for piles in sandy soil deposits.
• The proposed analytical models of the present study are subjected to static loads and can be extended to the effect of dynamic loading.
• The present analysis can be utilized to obtain the response of piles for other structures where the effect of pile-soil separation and slippage may have a significant effect on the non-linear pile response. Hence, the current nonlinear analysis framework can be extended to include the effect of pile-soil separation and slippage.
• In the present study, the load-displacement responses of single piles have been considered. However, the deformation responses of group piles is larger than the displacement of isolated single piles. Hence, the present work needs to be extended to understand the group action of the piles when subjected to several external loads.