Nonlinear Equilibrium and Stability Analysis of Axially Loaded Piles Under Bilateral Contact Constraints

This paper presents a nonlinear stability analysis of piles under bilateral contact constraints imposed by a geological medium (soil or rock). To solve this contact problem, the paper proposes a general numerical methodology, based on the finite element method (FEM). In this context, a geometrically nonlinear beam-column element is used to model the pile while the geological medium can be idealized as discrete (spring) or continuum (Winkler and Pasternak) foundation elements. Foundation elements are supposed to react under tension and compression, so during the deformation process the structural elements are subjected to bilateral contact constraints. The errors along the equilibrium paths are minimized and the convoluted nonlinear equilibrium paths are made traceable through the use of an updated Lagrangian formulation and a Newton-Raphson scheme working with the generalized displacement technique. The study offers stability analyses of three problems involving piles under bilateral contact constraints. The analyses show that in the evaluation of critical loads a great influence is wielded by the instability modes. Also, the structural system stiffness can be highly influenced by the representative model of the soil.


INTRODUCTION
The structural elements used to transfer load from a super-structure to the geological medium include such elements as beams, piles, arches, plates, and shells.These elements, during the deformation process, can lose contact with the geological medium due to the medium's incapacity to react under tension.This kind of problem, known as a unilateral contact problem, leads to important differences in the reactions of foundation elements and in the internal stress of structural elements.Such differences result in the remaining contact area being subjected to high stress concentrations.Researchers have, since the 70's, actively studied the unilateral contact problem (Silveira et al., 2013;Celep et al., 2011;Sapountzakis and Kampitsis, 2010;Silveira et al., 2008;Holanda and Gonçalves, 2003;Silva et al., 2001;Silva, 1998;Silveira, 1995;Stein and Wriggers, 1984;Weitsman, 1970).A separate type of problem is the bilateral contact problem (Yu et al., 2013;Maciel, 2012;Shen, 2011 and2009;Kien, 2004;Naidu and Rao, 1995).Bilateral contact problems involve support systems found in civil and mine engineering where structural elements are completely attached to the geological medium.Examples include piles embedded in soils, underground pipelines, and circular or parabolic tunnel roofs.In the context of physical and geometric linearity, these engineering problems can be treated as simple structural minimization problems; numerical solutions, and even analytical solutions, can be easily achieved (Tsudik, 2013;Selvadurai, 1979;Hetényi, 1946).Some contexts, however, involve non-linear effects.In recent years, many studies have considered such effects in their analyses of bilateral contact problems under extreme static and dynamic loads (Fazelzadeh and Kazemi-Lari, 2013;Challamel, 2011;Sofiyev, 2011;Roy and Dutta, 2010;Dash et al., 2010;Manna and Baidya, 2010).In nonlinear structural and geotechnical problems with bilateral constraints, the property and characteristics of the geological medium can highly influence the support systems' equilibrium paths.Researchers trying to resolve these kinds of issues have made significant advances in several engineering fields leading to the development of new materials and techniques of analyses.Perhaps the greatest contribution in these advances has been the development of more precise computational tools, whose models allow a realistic simulation of problems (Silva, 2009).To come up with slenderer and lighter support systems and constructions, without compromising quality and security standards, engineers have had to evaluate a structure's behavior by considering its nonlinear effects as well as the structure's interaction with the soil-structure (Kausel, 2010).
Nonlinear bilateral structural contact problems are usually formulated by using numerical techniques (Maciel, 2012;Mullapudi and Ayoub, 2010;Matos Filho et al., 2005;Horibe and Asano, 2001), such as finite element method (FEM) or boundary element method (BEM).Regarding the contact constraints, one is able to transform the contact problem-by using the usual formulations of structural mechanics-into a minimization problem without constraint.As the problem's solution depends heavily on the constitutive model of geological medium, one can expect to find the best results by adopting a more rigorous mechanical model.Kausel (2010), Wang et al. (2005), and Dutta and Roy (2002) provide a state-of-the-art review for soil-structure interaction.To solve this kind of contact problem, these researchers suggest bars and plates on elastic foundation and soil modeling as well as analytical and numerical possibilities.
Latin A m erican Journal of Solids and Structures 12 (2015) 250-270 Classical references on structural stability (Brush and Almroth, 1975;Simitses and Hodges, 2006) involve analytical solutions of piles under contact constraints imposed by Winkler-type foundations.These solutions provide the pinned-pile critical load expression as a function of the number of half-wave imperfection mode and the dimensionless foundation parameter.Ai and Han (2009) and Ai and Yue (2009) studied the behavior of piles embedded in a multi-layered soil and subjected to axial load.Researchers have also examined such subjects as load-deflection response, stability analysis, and the buckling and initial post-buckling behavior of bilaterally constrained piles, as found in Fazelzadeh and Kazemi-Lari (2013), Challamel (2011), Chen and Baker (2003), Sironic et al. (1999), Chai (1998), Budkowska (1997a,b), Budkowska and Szymczak (1997), West et al. (1997), and Naidu and Rao (1995).Recently, Tzaros and Mistakidis (2011) presented the critical loads and buckling modes of columns under unilateral contact constraints.Limkatanyu andSpacone (2002, 2006) presented three formulations of frame elements with nonlinear lateral deformable supports.To analyze the problem of a pile partially buried in an elastic medium, Aljanabi et al. (1990) developed a contact finite element that took into account the normal stiffness of soil and the soil-structure friction.Badie and Salmon (1996) solved the same problem but by using a quadratic order elastic foundation finite element.Besides the normal stiffness of soil and the soil-structure friction, they considered the interaction between the individual springs, such as those used in the Pasternak-type foundation model.
The importance of considering the Pasternak's model to represent the ground was explored by Shirima and Giger (1992) using a Timoshenko-beam element that incorporated both stiffness parameters of the base.The model's importance was also investigated by Mullapudi and Ayoub (2010) who analyzed a bar of finite size in contact with sandy clay based on a mixed formulation (approximations of forces and displacements).The stability of piles on a Pasternak foundation has also been examined in Naidu and Rao (1995), Kien (2004), andShen (2011).Horibe and Asano (2001) used the boundary element method (BEM) to study bars in contact with a Pasternaktype base (or Filonenko-Borodich), considering large displacements.
The main objective of this paper is twofold.First, it proposes a numerical methodology for the geometrically nonlinear analysis of piles under bilateral contact constraints.Second, it studies the influence of the geological medium and its stiffness on the piles' nonlinear equilibrium path and buckling behavior.To model the piles, we use a nonlinear beam-column element (Alves, 1995;Silva, 2009); to approximate the geological medium's behavior, we use individual springs (discrete model) or a bed of springs (continuum models) that exhibit a non-sign-dependent forcedisplacement relationship (Silva, 1998;Maciel, 2012).
An updated Lagrangian formulation is adopted to follow the system's movement, ignoring the influence of friction in the contact area.The contact regions between the bodies are known a priori and the nonlinear problem involves as a variable only the pile-soil displacement field.To solve the resulting algebraic nonlinear equations with contact constraints, and to obtain the structural equilibrium configuration at each load step, the present work presents an iteration solution strategy using the Newton-Raphson scheme coupled with the generalized displacement method (Crisfield, 1991;Yang and Kuo, 1994).
To verify the proposed numerical solution, this paper presents three problems.The first one investigates the influence of the position and stiffness of a spring elastic support on the pile's crit-Latin A m erican Journal of Solids and Structures 12 (2015) 250-270 ical load (Brush and Almroth, 1975); the second analyzes the elastic stability behavior of slender pinned piles under contact constraint imposed by a Winkler-type foundation throughout its length (Brush and Almroth, 1975;Simitses and Hodges, 2006); the third studies the nonlinear behavior and stability of piles with different end conditions in contact with a Pasternak-type foundation (Naidu and Rao, 1995;Kien, 2004;Shen, 2011).These examples demonstrate the accuracy and versatility of the present numerical strategy in the nonlinear solution of piles under bilateral constraints.

THE BILATERAL CONTACT PROBLEM
Consider the support system shown in Figure 1a and its mathematical model in Figure 1b.It consists of a pile and an elastic foundation.Assume that both bodies may undergo large deflections and rotations but with only small strains that are within the material's elastic range.Assume also that the contact surface is bonded and frictionless, and the region S c corresponds to the region where contact occurs, which is known a priori.Consider now that the variables are known at 0 and t equilibrium configurations, and the solution at t + t configuration is required (Figure 1c).Particularly well suited for numerical analysis, the nonlinear bilateral contact problem can be solved through the following minimization problem (Maciel, 2012): with the functional written as: in which  is the Cauchy's stress vector; S is the 2 nd Piola-Kirchhoff's stress vector;  is the Green-Lagrange's strain vector; rb is the reaction vector of the elastic foundation; ub is the displacements vector of the elastic foundation; u is the displacements vector of the structure and F is the external forces vector specified on Sf and assumed to be independent of the bodies' deformations.In addition,  represents an incremental quantity while superscripts t and t + t define the reference configuration.
The equality (2) gives the contact condition  the gap in the potential contact area is always zero  after the increment of the displacements, that must be satisfied on Sc in the configuration t + t.According to Figure 2, Sc can be one contact region, some regions, or even some discrete points.Thus, for a given load increment, the unknown variables in the configuration t + t may be obtained by solving the above minimization problem.What makes the problem difficult to solve directly is the geometrically nonlinear nature of this analysis.The following sections demonstrate how this kind of nonlinearity may be treated.Solids and Structures 12 (2015) 250-270 3 THE GEOMETRICALLY NONLINEAR ANALYSIS This section presents the numerical solution process of the geometrically nonlinear problem.The first step is to discretize the structural problem by using the finite element method.In this context, one can assume that, for a generic structural element, the incremental displacement field u within the element is related to the incremental nodal displacements u  by: where H is the usual FE interpolation functions matrix.
To evaluate the corresponding incremental strains and stresses, one can write the Green-Lagrange increment tensor and the 2 nd Piola-Kirchhoff stress increment tensor using the following expressions (Bathe, 1996):

SC
(6) where B L is the same strain-displacement matrix as in the linear infinitesimal strain analysis and it is obtained by appropriately differentiating and combining the rows of H ; B NL depends on H and the incremental displacements; and C is the constitutive matrix.
The foundation soil's action, or simply the foundation's action, on the structural elements depends on the displacements and stresses on the foundation soil itself.However, as in many engineering applications the designer is interested in only the response of the foundation at the contact area and not in the stresses and displacements in the ground foundation, it is possible to develop a simple mathematical model to describe, with a reasonable degree of accuracy, the response of the foundation soil at the contact zone.This can be done by using the well-known Winkler's model (Kerr, 1964) or even the formulation for the elastic half-space (Cheung, 1977).Thus, the elastic foundation's behavior for a generic element can be written by the following discrete equilibrium equation: where rb and ub are the incremental elastic foundation reaction and displacement nodal values, respectively; and, C b is the constitutive matrix of the elastic foundation.Thus, for a generic elastic foundation element, the incremental displacement field ub may be related to its nodal where B b is the matrix containing the interpolation functions that describe the elastic base deformation.
Using the previous definitions and assuming that in the contact area the structure and elastic foundation nodal displacements are identical (i.e., ˆb uu ) and do not change, one arrives at the discretized functional of the contact problem, Eq. ( 3), at the element level: Latin A m erican Journal of Solids and Structures 12 (2015) 250-270 Taking the appropriate variations of  with respect to the incremental nodal displacements, and adding the contributions of each finite element, one can write: or, more concisely: where U contains the global nodal incremental displacements and t+  t Fi is the internal generalized forces vector of the support system in the load step t + t.Equation (10a) or (10b) is the equation that must be satisfied in an incremental process to obtain the system equilibrium.
On the left side of (10a), K L is the global stiffness matrix for small displacement; the matrix K  is the initial stress matrix or geometric matrix; the matrix K NL is the large displacement matrix, which contains only linear and quadratic terms in the incremental displacement; and K b is the stiffness matrix of the elastic foundation (Maciel, 2012).All these matrices contribute on the assemblage of the tangent stiffness matrix used in numerical solution procedure presented on the next section.Vectors t Fis represents the internal force vector of the structure and vector t Fib the elastic foundation in the equilibrium configuration t.They are typically computed by integrating the generalized stress resultants through the volume of each element and then summing the elemental contributions (Bathe, 1996).
On the right side of (10a), vector t+  t R is the nodal external load vector in the step t+t and is given by: which is assumed to be independent of the structure's deformation; R r is a fixed load vector, termed reference vector defining the load direction, and t+  t  is a scalar load multiplier which defines the intensity of the applied load.

NUMERICAL SOLUTION PROCEDURE
This section presents the main characteristics of the numerical solution strategy adopted for the minimization problem defined by Eqs. ( 1) and ( 2), considering the geometric nonlinearity and the bilateral contact constraints.
To obtain nonlinear equilibrium paths, this study considers an incremental and iterative solution strategy.It is assumed that perfect convergence was achieved in the previous load steps 0,…,t (Figure 1c), i.e., the solution of the previous steps satisfied the equilibrium equations and all contact constraints.Therefore, considering the updated Lagrangian formulation  which is more appropriate to obtain the tracking of interfaces between different materials or surfaces and also due the nonlinear finite element formulation adopted (Section 5) , the known displacements and stresses obtained at the conclusion of load step t are used as information to obtain the adjacent equilibrium configuration t + t.A cycle of the incremental and iterative strategy can be summarized in two steps: 1.As a starting point, the strategy employs an approximate solution, called a "tangential incremental solution."This approximate solution involves selecting an initial increment of the load.To calculate the initial load increment, the strategy then considers an additional constraint equation, which uses the "generalized stiffness parameter" equation (Yang and Kuo, 1994); Sc defines the contact regions between the bodies and during the incremental process remains constant.To calculate the initial increment of the displacements, the strategy employs this initial load increment approximation,  0 ; it also uses it to define what is called "tangential incremental solution," a solution that rarely satisfies the equilibrium equations.Hence, the strategy must use a correction.2. This correction deals with the geometric nonlinearity of the problem.It can be achieved by coupling the Newton-Raphson method and path-following techniques (continuation methods; Crisfield, 1991).The strategy uses the generalized displacement method developed by Yang and Kuo (1994) to allow limit points to be passed and, consequently, to identify snap buckling phenomena.If the convergence criterion is not satisfied, the correction procedure is repeated until it is.The numerical solution procedure summarized above is detailed in Figure 3.

NUMERICAL APPLICATIONS
In this section, the methodology presented for nonlinear static analysis is used to obtain the responses of three stability problems involving piles under bilateral contact constraints.Consistent units are used in all problems and the proposed numerical solution strategy considers a convergence tolerance factor equal to 10 -4 .To model the piles, it is adopted the nonlinear beam-column element developed by Alves (1995) and improved by Silva (2009).Basically, the nonlinear beamcolumn element formulation presents as main characteristics: the adoption of the Bernoulli-Euler simplifications (all sections remain normal to element axis and cross-sections remain plane after loading application); the shear deformation effect is disregarded; the geometrically nonlinear kinematics follow, as already mentioned, the Green-Lagrange tensor; the internal force vector increment is obtained using the natural displacement approach; and linear material behavior.
The proposed finite element contact problem formulation as well as the numerical solution methodology (Section 4) were implemented and adapted in FORTRAN language in a new computational system for advanced structural analysis, CS-ASA, which performs the nonlinear static and dynamic analyses of steel members and framed system structures (Silva, 2009).Brush and Almroth (1975) presented the analytical solution for the critical load (Pcr) of a pinned pile, represented by a bar with length L and bending stiffness EI, with an intermediate elastic spring, represented by its linear stiffness Kx in the horizontal direction X, located at a distance c from the pile top, as depicted in Figure 4.

Piles with Intermediate Elastic Spring Support
To investigate the influence of the intermediate elastic spring on the pile's critical load, several analyses were performed, maintaining as constant the length L and the bending stiffness EI of the pile and varying the position c and the linear stiffness Kx of the intermediate elastic spring.
Latin A m erican Journal of Solids and Structures 12 (2015) 250-270 A lgorithm : SUPPORT SYSTEM NONLINEAR STABILITY ANALYSIS 1: Read the input data: geometric, material and loading properties of the structural system 2: Define the contact regions Sc and obtain the reference nodal load vector, Fr (loading direction) 3: t = 0 4: t1 = t 5: Select the initial load increment ( 0 )1, Crisfield (1981) 6: for each load increment do

18:
Evaluate the initial incremental displacement vector: Evaluate the internal forces vector:

21:
Evaluate the residual force vector:
Presented in Figure 4 are the results obtained by Brush and Almroth (1975), using different dimensionless stiffness parameter  =  2 KxL/PE, in terms of the normalized position (c/L) and normalized critical load (Pcr/PE), where PE =  2 EI/L 2 is the Euler's critical load.
As can be observed in Figure 4, the numerical results obtained by this work are in good agreement with those presented by Brush and Almroth (1975).For a normalized position of 0.5 and  higher than 150, the elastic spring can be considered as a rigid support, with Pcr  4PE, which is the critical load of a pinned pile of length L/2.To behave as a rigid support, an intermediate elastic spring located at the normalized position less than 0.5 should present a parameter higher than 150.For an elastic spring with high stiffness and located near the pile top (c/L  0), the value of Pcr tends to 2.05 PE, which is equal to the fixed-pinned pile buckling load.
The same analyses performed with a pinned pile were conducted while adopting a fixedfree end condition.Results are presented in Figure 5.It should be noted that for  higher than 100, the discrete elastic spring already behaves as if it were a rigid support, regardless of the spring location.As expected, for the elastic spring located near of top pile (c/L = 0) and with high values of  the critical load obtained is about 2.05 PE, which is the critical load of a fixed-pinned pile.

Pinned Piles in Contact with a Winkler Type Foundation
Let us now look at the classical structural and geotechnical problem shown in Figure 6, i.e, the elastic stability analysis of slender pinned piles under a bilateral contact constraint imposed by a Winkler-type foundation throughout its length.Brush and Almroth (1975) and Simitses and Hodges (2006) showed that, as in the elastic instability of plates and shells, the buckling mode plays an important role in the stability of this type of problem.This means that the number of semi-waves' modes on the analytical solution has great influence on the pile's critical load value.According to these authors, the pile's critical load could be calculated by the following expression: where n is the number of semi-waves considered, w = kL 2 /( 2 PE) is the dimensionless stiffness parameter, and PE =  2 EI/L 2 is Euler's critical load.
For a pinned pile such as the one illustrated in Figure 6 with characteristics of L = 5, EI = 100 and k = 10 (all in compatible units) and considering n = 1 (one semi-wave), it is, according to Eq. ( 12), possible to obtain a pile critical load of 64.8.
To verify the influence of the discretization on the critical load, different analyses were conducted varying the number of finite elements and the elastic base models (discrete and continuum or Winkler-type).Table 1 presents the pile critical loads and the errors (related to the analytical solution of 64.8) observed when varying the number of finite elements from 4, 6, 8, 10, and 20.Note the good agreement for both models for the mesh with eight finite elements.Figure 7 shows the support system equilibrium paths for the two models by using a mesh with four finite elements.To investigate the influence of the elastic foundation stiffness of a Winkler base model on the pile's critical load, two analyses were performed maintaining as constant the length (L = 10) and the bending stiffness (EI = 100) of the pile and varying the dimensionless stiffness parameter βw. Figure 8 presents the numerical results in terms of the equilibrium paths for different numbers of semi-waves considering the dimensionless stiffness parameter βw of 16 and 48.The numerical results are in a good agreement with those presented by Brush and Almroth (1975).Table 2 presents the results in terms of pile critical load and mode.It can be observed that for w = 16, the critical mode obtained was n = 2 and the critical load was Pcr = 78.51.For w = 48, the critical mode obtained was n = 3 and the critical load Pcr = 141.41.These results show how important initial imperfections and their shapes are regarding stability analysis of piles under contact constraints.

Piles under Contact Constraint Imposed by a Pasternak Type Foundation
Finite element solutions of piles, with different end conditions, and under bilateral contact constraints of the Pasternak type (Figure 9) were initially presented by Naidu and Rao (1995).Kien (2004) and Shen (2011) also presented the buckling loads for the particular case of the pinned pile with bilateral constraints imposed by Winkler-and Pasternak-type foundations.These references are then used to validate the stability analyses carried out using the proposed numerical methodology, where for all piles the following were considered: 20 finite elements, The nonlinear equilibrium paths for these analyses are presented in Figs.8a-c, in terms of the dimensionless pile load  (= PL 2 /EI) considering four combinations of dimensionless parameters 1 and 2: (1 = 1; 2 = 0); (1 = 100; 2 = 0); (1 = 100; 2 = 0.5); and (1 = 100; 2 = 2.5).The figures reveal the good agreement between the results obtained in this article with those of literature.Tables 3-5 present the dimensionless pile critical load cr for different combination of the 1 and 2 stiffness parameter, different end conditions, and different authors.
For the combination of 1 and 2: (0; 0), which represents the classical column stability problems (without contact constraints), note (Table 3-5) that the values obtained for the dimensionless critical load cr through the present work, as well as in the literature, are quite close to those of the analytical solution for the three columns, namely: 2.4674, 9.8696, and 39.4784.
For the 2 parameter equaling zero, combinations (1; 0) and (100; 0), the Pasternak models are reduced to Winker models.In this situation, for pinned pile, Brush and Almoth (1975) obtained the following values of the dimensionless critical load, cr: 9.9681 and 20.0051.The values found through the present work for cr show good agreement with those from the literature.
For the elastic foundation represented by the Pasternak model and combination 1 and 2, (100; 2.5), the following values of the ratio cr (Pasternak)/ cr (without contact) were obtained for the three piles: 14.8, 4.5, and 1.8.Thus, the fixed-free pile was more sensitive to additional stiffness provided by the Pasternak elastic foundation.Good agreement can also be observed between the results of the present paper and those found in the literature.

CONCLUSIONS
This article evaluated the equilibrium and stability of piles under bilateral contact constraints imposed by a geological medium.The numerical strategy proposed for solving the geometric nonlinear problem was based on the finite element method and the Newton-Raphson method and the generalized displacement approach.The examples presented above validate the computational implementations.
Basically, the nonlinear analyses carried out had as objectives: to evaluate the critical load of piles with a discrete elastic intermediate support; to evaluate the imperfection (unstable modes) influence in the assessment of piles critical loads in contact with a Winkler-type foundation; and verify the influence on the stiffness system when considering the elastic base second parameter, i.e., by adopting the Pasternak model in representation of the soil.The nonlinear results presented have made it possible to establish some general conclusions that are summarized below: i. the higher the pile's critical buckling load-in which sideways movement is constrained by a single elastic spring support-the higher the value of elastic spring stiffness that acts as a rigid support and, consequently, the geometric condition, specific to each case, changes the buckling mode;   ii. to analyze the pile's instability with continuum bilateral contact constraint imposed by an elastic foundation, initial imperfections, together with the mathematical model used to represent the foundation, has great influence in determining the pile's critical load; iii. the numerical formulation presented in this work was able to define with precision the pile's critical load as well as the post-critical response; iv.when adopting the Pasternak model to represent the soil, the higher the pile's critical buckling load, the smaller the influence of the model's second parameter.

Figure 1 :
Figure 1: Pile under bilateral contact constraints imposed by an elastic foundation.
Figure 2: Different domain situations for the definition of the contact region Sc.

Figure 7 :
Figure 7: Nonlinear equilibrium paths of pile -discrete and continuum models.

Table 1 :
Pile critical load: influence of the discretization.

Table 2 :
Pile critical load and mode.

Table 3
Fixed-free pile: dimensionless critical load cr.

Table 5
Fixed pile: dimensionless critical load cr.