Seismic behavior of soil–pile–structure interaction with a modified Desai thin-layer interface element

Abstract

The Desai thin-layer interface element is widely utilized in the simulation of interaction between piles and soil under seismic load. Conventional seismic analysis using the interface element cannot simulate the process of energy dissipation because tangential damping is disregarded. In this study, Rayleigh damping is added to the interface element to simulate energy dissipation in a strong nonlinear contact behavior. A user-defined element program based on a modified Desai interface element is developed. A hyperbolic model is adopted to simulate normal and tangential interaction behaviors. Certain behavior pattern rules of the modified Desai element, such as bonding, slipping, gapping, and reclosing, are defined. A three-dimensional pile–soil–structure model with a modified Desai interface element is established to investigate the effect of contact patterns on the inner force responses of a superstructure and pile foundation to an earthquake. Numerical results show that the contact patterns significantly influence the shear force, bending moment, and torque of the superstructure, while axial force is unaffected. With regard to the pile foundation, shear force and bending moment are also significantly influenced.

Introduction

The scale of superstructures and group-pile foundations is increasing because of rapid urbanization. These large-scale structures can directly change the dynamic property of the soil surrounding foundation, and the changed dynamic property may influence the dynamic properties of a superstructure, such as its frequency and model. In addition, the dynamic characteristics of a superstructure and pile foundation may be significantly influenced by each other during an earthquake. ¹

Several studies on soil–pile–structure interaction have been conducted with different methods, such as finite element method (FEM) and boundary element method (BEM),^2–4 direct shearing test, ⁵ numerical analysis,^6–8 and interface element method. ⁹ During an earthquake, the deformation of pile is different from that of soil because of the differences between the material property of pile and that of soil. On the interface between the soil and pile, slipping or gapping may occur. The contact between soil and pile involves several complex mechanical problems, such as deformation discontinuity, deformation nonlinearity, and large deformation. ¹⁰ Strong nonlinear contact can affect the stress-state solution of soil near the interface and the seismic response level of superstructures. Therefore, the interface model is a key tool in studying the seismic responses of certain pile–soil–structure systems. ¹¹

The contact zone between a pile and soil can be simulated with interface elements. ¹² The zero-thickness Goodman element ¹³ and the Desai thin-layer element ¹⁴ are the most commonly used elements. The interface between a pile and soil is not a simple two-dimensional (2D) curved face with zero thickness but a contact zone with nonzero thickness. Because of this zone, the soil constrained by the rough surface of concrete structures presents different mechanical properties from the soil out of the contact zone. Thus, the Desai thin-layer element is the focus of this study. Desai and Drumm ¹⁵ simulated the contact behaviors on the surface between soil and structures, such as bonding, slipping, gapping, and reclosing by employing Desai thin-layer element. The numerical procedure is used to predict the behavior of a model structure tested in the field, and the results show that the Desai element can precisely describe the dynamic performance of the interface. The Desai thin-layer element can transmit force and can be employed to simulate contact behaviors, such as bonding, slipping, gapping, and reclosing. The specific process can be expressed as follows: ¹⁶ in the initial stage of analysis, the normal stress of the thin-layer element is smaller than tensile strength, and tangential stress is smaller than shear strength, so the interface element is bonded to the pile; when tangential stress is greater than shear strength and normal stress is smaller than tensile strength, interface elements start to slip, and tangential stress equals to shear strength, while normal stress is unchanged; when the normal stress of an element is greater than tensile strength, the interface element separates and the normal stiffness of the separated element is zero. Once the distance between the interface element node and the pile node reaches a very small value, normal and tangential stiffness recover, ¹⁷ which means that the interface element returns to the bond state. In previous studies on the Desai element, damping was not considered, while damping is crucial in simulating energy dissipation in strong nonlinear contact behavior.

The Desai thin-layer element has been widely applied. Muqtadir and Desai ¹⁸ utilized the thin-layer interface element in a pile-group foundation model to include interaction effects involving relative slipping and gapping. Fishman et al. ¹⁹ employed thin-layer elements to model the response of rock joints in finite element analyses of rock mechanics problems. This study provided valuable data that contribute to the verification of thin-layer finite elements for modeling the response of rock joints. Shao et al. ²⁰ proposed a new interface model, combining the advantage of Goodman element and Desai thin-layer element. In this model, the method of controlling penetration used in thin-layer element by Desai is introduced into the Goodman method. In this model, nonlinear characteristics of interface in normal and shear directions are considered, and a mixed iteration method of residual stress redistribution and penetration controlling is utilized.

In this study, the original Desai thin-layer element is modified by adding Rayleigh damping ²¹ to the constitutive relationship to make it more realistic and rational for seismic analysis. The modified Desai interface element is employed to simulate the contact zone between soil and pile. The three-directional ground motions are generated artificially. Seismic response analyses of the pile–soil–structure system with different contact patterns are conducted in ABAQUS, and the results are compared.

Modified thin-layer interface element

The three-dimensional (3D) Desai elements in true space and parameter space are shown in Figure 1, respectively; herein $\zeta$ denotes the normal direction of the elements. This type of element is anisotropic. The element displacement interpolation functions are

$$\begin{array}{ccc}u={\displaystyle \sum _{i=1}^8{N}_i{u}_i,}& v={\displaystyle \sum _{i=1}^8{N}_i{v}_i,}& w={\displaystyle \sum _{i=1}^8{N}_i{w}_i}\end{array}$$

(1)

where $u_i$ , $v_i$ , and $w_i$ denote the displacements of node i in 3D space and $N_i$ denotes the shape function which can be expressed as

$$N_i=\frac{1}{8}\left(1+{\xi }_i\xi \right)\left(1+{\eta }_i\eta \right)\left(1+{\zeta }_i\zeta \right)$$

(2)

Figure 1.

3D Desai interface element.

The coordinate transformation can be realized by

$$\begin{array}{c}x=\sum _{i=1}^8{N}_i{x}_i,y=\sum _{i=1}^8{N}_i{y}_i,z=\sum _{i=1}^8{N}_i{z}_i\hfill \end{array}$$

(3)

The Jacobian matrix of the 3D solid element can be expressed as

$$J=\left[\begin{array}{ccc}\frac{\partial x}{\partial \xi }\hfill & \frac{\partial y}{\partial \xi }\hfill & \frac{\partial z}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \eta }}& \frac{\partial y}{\partial \eta }\hfill & \frac{\partial z}{\partial \eta }\hfill \\ \multicolumn{1}{c}{\frac{\partial x}{\partial \zeta }}& \frac{\partial y}{\partial \zeta }\hfill & \frac{\partial z}{\partial \zeta }\hfill \end{array}\right]$$

(4)

The element stiffness matrix is

$${\left[K\right]}^e=\underset{V_e}{\int }{\left[B\right]}^T\left[D\right]\left[B\right]\mathrm{dV}=\int \int \int {\left[B\right]}^T\left[D\right]\left[B\right]\left|J\right|d\xi d\eta d\zeta$$

(5)

The incremental constitutive condition of the 3D Desai element can be expressed as

$$\left\{d\sigma \right\}={\left[D\right]}_i\left\{d\epsilon \right\}$$

(6)

which can be expanded as

$$\left\{\begin{array}{c}\mathrm{d}{\sigma }_x\\ \mathrm{d}{\sigma }_y\\ \mathrm{d}{\sigma }_z\\ d{\tau }_{xy}\\ d{\tau }_{yz}\\ d{\tau }_{zx}\end{array}\right\}=\left[\begin{array}{cccccc}{D}_1& D_2& D_2& & & \\ D_2& D_1& D_2& & & \\ D_2& D_2& k_{zz}& & & \\ & & & D_3& & \\ & & & & k_{yz}& \\ & & & & & k_{zx}\end{array}\right]\left\{\begin{array}{c}d{\epsilon }_x\\ d{\epsilon }_y\\ d{\epsilon }_z\\ d{\gamma }_{xy}\\ d{\gamma }_{yz}\\ d{\gamma }_{zx}\end{array}\right\}$$

(7)

where $\{d\sigma \}$ and $\{d\epsilon \}$ denote the matrix of stress and strain increment, respectively; $[D{]}_i$ is the elasticity matrix; x and y denote the tangential direction; and z denotes the normal direction of the interface; and $D_1$ , $D_2$ , and $D_3$ are $D_1=(E(1-\mu ))/((1+\mu )(1-2\mu ))$ , $D_2=E\mu /((1+\mu )(1-2\mu ))$ , and $D_3=E/(2(1+\mu ))$ , respectively.

Desai indicated that the normal stiffness of a thin-layer element is related to the material property of nearby soil elements and structural elements. Thus, the normal stiffness matrix of thin-layer element can be expressed as

$${\left[{\overline{D}}_n\right]}_i={\lambda }_1{\left[D_n^g\right]}_i+{\lambda }_2\left[D_n^g\right]+{\lambda }_3\left[D_n^{\mathrm{st}}\right]$$

(8)

where $[{\overline{D}}_n{]}_i$ , $[D_n^g]$ , and $[D_n^{\mathrm{st}}]$ denote the normal stiffness matrices of soil in contact zone, soil, and structural material. ${\lambda }_1$ , ${\lambda }_2$ , and ${\lambda }_3$ denote the weight factors. If ${\lambda }_1=1$ and ${\lambda }_2={\lambda }_3=0$ , then the material property of the elements on the sides of the Desai element is disregarded. According to the results of large number of cyclic loading tests, ${\lambda }_1$ , ${\lambda }_2$ , and ${\lambda }_3$ can be 0.75, 0.25, and 0, respectively.

According to the hyperbolic model, ²² normal stiffness can be expressed as

$$k_{\mathrm{zz}}=\frac{k_{\mathrm{zz}0}}{{(1+v/t)}^2}$$

(9)

where $k_{\mathrm{zz}0}$ denotes the initial normal stiffness on the interface, v denotes relative normal deformation, and t denotes the thickness of the interface. In this study, $k_{\mathrm{zz}0}$ is set to be $D_1$ .

The strain-hardening hyperbola model is utilized to simulate the tangential stiffness of the interface element, which can be expressed as

$$\{\begin{array}{c}{k}_{\mathrm{zx}}=k_{\mathrm{zx}0}{(1-\frac{1}{\cos \alpha }\frac{R{\tau }_{\mathrm{zx}}}{(c-f{\sigma }_{\mathrm{zz}})})}^2\hfill \\ \multicolumn{1}{c}{k_{\mathrm{yz}}=k_{\mathrm{yz}0}{(1-\frac{1}{\cos \alpha }\frac{R{\tau }_{\mathrm{yz}}}{(c-f{\sigma }_{\mathrm{zz}})})}^2}\end{array}$$

(10)

where $k_{\mathrm{zz}0}$ and $k_{\mathrm{yz}0}$ are initial tangential stiffnesses which are set to be $D_3$ and c, f, and $\alpha$ are the coherence force, the friction factor, and the friction shear slip angle, respectively. The asymptotic value of the hyperbola is set equal to $(c-f{\sigma }_{\mathrm{zz}})/R$ .

To simulate the process of energy dissipation caused by the contact between pile and soil during an earthquake, damping item is added to the stress–strain relationship. To truly embody the nonlinear variance of the soil near the pile and the stiffness change in the interface element, Rayleigh damping is applied. Then, the damping matrix can be expressed as

$${\left[C\right]}_{\mathrm{i}}=\alpha {\left[M\right]}_{\mathrm{i}}+\beta {\left[K\right]}_{\mathrm{i}}={\lambda }_{\mathrm{i}}{\omega }_{\mathrm{i}}{\left[M\right]}_{\mathrm{i}}+\frac{{\lambda }_i}{{\omega }_i}{\left[K\right]}_i$$

(11)

where ${\lambda }_i$ and ${\omega }_{\mathrm{i}}$ are the damping ratio and natural frequency of the soil near the pile, respectively. During pile–soil interaction, a negative correlation exists between the damping ratio of the soil near the pile and shear strain of soil. The nonlinearity of the soil damping ratio is addressed using equivalent linearization method.

A program of a user-defined element (UEL) is developed with Fortran language. The element node and Gauss integral point are arranged, and node force and node displacement are calculated. The modified Desai thin-layer element is 3D with eight nodes, which is compatible with the C3D8 element in ABAQUS. The developed program can be imported into ABAQUS directly through the subroutine interface.

Model verification

A simple pile–soil–structure interaction system was established to investigate the effect of damping in the modified Desai interface element on structural responses (Figure 2). The dimension of the soil model was 40 m × 40 m × 20 m. A concrete circular column was utilized as superstructure and pile. The superstructure part was 30 m long, and the pile part was 10 m. The elastic modulus of concrete was $3.3\times {10}^4\mathrm{MPa}$ , and the density was 2600 kg/m³. The cohesion of soil was 30 kPa, and internal friction angle was 20°. The soil damping ratio is assumed to be 0.09. Two horizontal cosine acceleration time histories $A_{x,y}=\cos (3.14t)$ and one vertical cosine acceleration time history $A_z=0.6\times \cos (3.14t)$ were input at the bottom of the model. The artificial viscoelastic boundary ²³ is utilized for the model edges and bottom. The modified and original Desai elements were employed to simulate the pile–soil interaction, respectively. Structural responses with different Desai interface elements were investigated.

Figure 2.

Simple pile–soil–structure system.

Figure 3 shows the acceleration time histories at the column cap. It can be seen that the trend of acceleration time histories of the model with damped Desai element is approximately the same as that of the model with undamped Desai element. The model with damped Desai element exhibits lower acceleration responses, and the maximum reduction can be 9.93%.

Figure 3.

Acceleration time histories at the column cap in the x-direction.

The envelope curves of acceleration time histories in the x-direction of the whole structure (column and pile) under damping and no damping occasions are shown in Figure 4. The acceleration response of the model with damped Desai element is lower than that of the model with undamped Desai element.

Figure 4.

The envelope curves of acceleration time histories in the x-direction.

Figure 5 shows the surface diagram of shear force and bending moment of column with damped and undamped Desai elements in the x-direction. It is shown that all inner force responses of the model with damped Desai element are lower than that of the model with undamped Desai element. The reduction in peak shear force and bending moment are 7.6% and 8.7%, respectively.

Figure 5.

The surface diagram of shear force and bending moment: (a) the shear force on no damping occasion, (b) the shear force on damping occasion, (c) the bending moment on no damping occasion, and (d) the bending moment on damping occasion.

Finite element model of the pile–soil–structure system

A pile–soil–structure system is established. The system includes soil and a frame structure with 8 floors, 4 columns and 4 piles for every column, as shown in Figure 6. The dimension of the soil is 50 m × 50 m × 30 m. For the piles, friction pile is utilized, and the diameter, length, and spacing distance are 0.4, 5, and 0.7 m, respectively. The pile cap dimension is 1.5 m × 1.5 m × 0.2 m. The dimensions of the frame are shown in Table 1.

Figure 6.

Pile–soil–structure system: (a) entire model and (b) frame structure and pile foundations.

Table 1.

Dimension of the frame.

Instance	Height of first floor (mm)	Height of standard floor (mm)	Beam (mm2)	Column (mm2)	Floor slab (mm)
Dimension	3800	3300	300 × 700	300 × 300	200

Material properties

Material properties for concrete

A concrete damage plastic (CDP) model is used for the superstructure and pile foundations to consider the effect of concrete failure on seismic responses. The specific factors of the CDP model are shown in Table 2. The density of concrete is $2600\mathrm{kg}/{\mathrm{m}}^3$ , and the elasticity modulus is $3.3\times {10}^4\mathrm{MPa}$ to consider the reinforcement effect of rebar.

Table 2.

Factors of the CDP model.

Dilation angle (°)	Eccentricity	Ratio of biaxial compressive strength to uniaxial compressive strength	Invariant stress ratio	Viscosity parameter
30	0.1	1.16	0.667	0.0005

Material properties for soil

The equivalent elastic model is utilized for the soil. The program of the equivalent elastic model is developed in Fortran language and imported into ABAQUS through a user-defined material (UMAT) interface. The specific factors of the equivalent elastic model are shown in Table 3.

Table 3.

Factors of the equivalent elastic model for soil.

Initial stiffness (N/m)	n	Damping ratio	Poisson’s ratio	Fundamental frequency (rad/s)	Maximum shear strain	Density (kg/m3)
1573.3	1	0.09	0.25	0.18	0.001	2000

Clarification of the equivalent elastic model

A simple free field model established in ABAQUS is used to clarify the equivalent elastic model as shown in Figure 7. The dimension of the model is 40 m × 40 m × 20 m. Three cosine acceleration time histories, as three-direction load, are input at the bottom of the model, which are $A_{x,y}=\cos (2.5t)$ for the two horizontal directions and $A_z=0.66\cos (2.5t)$ for the vertical direction. Besides, C3D8 element is employed for the whole model, and the element size is 2 m. The hypothetical confining pressure in the soil is shown in Table 4.

Figure 7.

Free field model.

Table 4.

The hypothetical confining pressure in the soil.

Depth (m)	5	10	15	20
Confining pressure (kPa)	100	200	300	400

Figure 8(a) and (b) shows the hysteretic stress–strain curve of top and bottom central elements, respectively. Residual plastic strain is not considered, which is consistent with the characteristic of the equivalent elastic model. It can be seen that the slope of hysteretic curve of the bottom element is larger than that of the top element, because the confining pressure at the bottom is larger.

Figure 8.

Hysteretic stress–strain curves of (a) top and (b) bottom central elements.

Figure 9 shows the acceleration time histories on the surface and bedrock, which are consistent basically. The consistence between these curves clarifies that the program of the equivalent elastic model can judge the inflection point accurately.

Figure 9.

Surface and bedrock acceleration response.

The artificial viscoelastic boundary is utilized for the model edges and bottom again. The velocities of S wave and P wave are $3500$ and $6000\mathrm{m}/\mathrm{s}$ , respectively. According to the National Code of China, ²⁴ three acceleration time histories, namely, one vertical and two horizontal ground motions, are generated artificially. The two horizontal ground motions include north–south (N-S) and a west–east (W-E) components. The peak values of the horizontal and vertical ground motions are $0.18$ and $0.12\mathrm{m}/{\mathrm{s}}^2$ , respectively. The duration time is 30 s, including a rise time of 2 s and an attenuation time of 10 s. All ground motions are inputted at the bottom of the model. For conciseness, only the N-S earthquake ground motion is shown in Figure 10. Figure 11 shows the response spectrum of the N-S ground motion.

Figure 10.

N-S acceleration time history.

Figure 11.

Response spectrum of N-S acceleration time history.

Discussion of the results

To study the effect of the modified Desai thin-layer element on the inner force responses of the 3D soil–pile–structure system, five different types of contact between pile and soil are applied in the model. A numerical analysis is performed for each case, and the results are compared. The study cases are shown in Table 5.

Table 5.

Matrix of the analyses performed.

Cases	Contact pattern
Case 1	Desai thin-layer element with damping
Case 2	Desai thin-layer element without damping
Case 3	Hard contact and coulomb friction
Case 4	Hard contact and no friction
Case 5	Nonmatching grid coupling

Analysis of the inner force responses of the superstructure

Figure 12 shows the normalized peak inner force at the top of the frame column for each floor, including peak shear force, peak bending moment, peak axial force, and peak torque. Normalization is realized by dividing the results of cases 1–4 by that of case 5. In Figure 12, the value of case 2 is greater than that of case 1, except for the ratio of peak axial force. The trend of cases 1 and 2 is definitely different from that of cases 3 and 4 except for the ratio of peak axial force. This is because that the interaction methods of cases 1 and 2 and cases 3 and 4 are significantly different. The method based on interface element takes the variation of material property of soil and pile and a lot of contact behaviors into consideration, including bonding, slipping, gaping, and reclosing, which can accurately simulate the interaction between pile and soil. However, the method based on hard contact, friction, or frictionless behavior regards the pile and soil as rigid body without considering the material property variation.

Figure 12.

Normalized peak inner force at the top of the frame columns for each floor: (a) ratio of N-S peak shear force of the frame column, (b) ratio of N-S peak bending moment of the frame column, (c) ratio of peak axial force of the frame column, and (d) ratio of peak torque of the frame column.

Figure 12(a) shows the comparison of the peak shear force of the frame columns, in which the trends of cases 1 and 2 and cases 3 and 4 are different. The curves of cases 1 and 2 fluctuate greatly, while the curves of cases 3 and 4 are stable. A comparison of the peak bending moments of the frame columns is shown in Figure 12(b). The curves of case 3 and case 4 descend monotonically from the first floor to the seventh floor, whereas the curves of cases 1 and 2 ascend, reach the maximum at the seventh floor, and then descend.

In Figure 12(c), the trend of each curve is identical, and all ratios are very close to 1; thus, the contact pattern has minimal influence on the axial force of the frame column. For the peak torque of the frame column, the results of cases 1–4 are all smaller than that of case 5. The results of cases 1 and 2 are obviously smaller than that of cases 3 and 4 from the first floor to the sixth floor. Numerical comparison of the peak inner force of the frame column between cases 1 and 2 is shown in Table 6.

Table 6.

Comparison of peak inner force of frame column between cases 1 and 2.

Difference	Shear force		Bending moment		Axial force	Torque
	N-S	W-E	N-S	W-E
Maximum decrease	7.9%	8.5%	7.4%	8.1%	0	16.1%
Minimum decrease	2.1%	2%	1.1%	0.8%	0	9.3%

Analysis of the displacement responses of the superstructure

The absolute displacement responses of superstructure of cases 1–5 are analyzed and compared as shown in Figure 13, in which the displacement responses of cases 1 and 2 are more significant than that of cases 3–5. The vibration components of cases 1 and 2 are more abundant. This is because that in the contact methods used in cases 3–5, two contact objects are treated as rigid approximately, and the effect of constitutive relation variance on the interaction is not considered. In addition, in case 1, the Rayleigh damping is considered which contributes to the energy dissipation in the soil–pile interaction and to some extent constraints the displacement of pile so that the peak displacement of case 1 is smaller than case 2. Although there is no friction in case 4, hard contact is used, which provides more constraint for the interaction between soil and pile. This is because that Desai thin-layer element-based method considers the gaping and closing thoroughly, which lowers the constraint of soil to pile. However, hard contact–based method does not. The envelope curve of N-S displacement time history of the frame column is shown in Figure 14, in which the displacement responses of cases 1 and 2 are larger than those of cases 3–5 in different levels. The peak displacement of case 1 is smaller than that of case 2. The difference decreases along with the distance from column cap increases.

Figure 13.

The comparison of N-S absolute displacement time history on the top corner point of the superstructure in the first 12 s.

Figure 14.

The envelope curve of N-S absolute displacement time history of the frame column.

Analysis of the peak inner force of the group-pile foundation

The inner force of group-pile foundations is analyzed to determine the effect of the improved Desai thin-layer element. Figures 15 and 16 show the normalized N-S peak shear force and peak bending moment of piles, respectively. Note that two diagonal piles of one of the four group-pile foundations are chosen to be investigated which are named piles 1 and 2. Piles 1 and 2 are at the outside corner and inside corner, respectively. Normalization is realized by dividing the results of cases 1 and 2 by that of case 5.

Figure 15.

Normalized N-S peak shear force of (a) pile 1 and (b) pile 2.

Figure 16.

Normalized N-S peak bending moment of (a) pile 1 and (b) pile 2.

The peak shear force of case 1 is smaller than that of case 2 because of the damping in the interface element. It can be concluded that the Desai interface element without damping overestimates the shear force response of the pile. The peak shear force of case 1 is very close to that of case 2 at the top and bottom of the pile. For both cases 1 and 2, the maximum shear force appears in the position 1 m from the top of the pile, and the difference between the peak shear forces of cases 1 and 2 is more obvious in this position. For piles 1 and 2, the peak shear force of the pile of case 1 is 34.9% and 39.6% lower than that of case 2 in the position 1 m from the top of the pile, respectively.

Figure 16 reveals that the peak bending moment of pile 1 and 2 of case 1 is smaller than that of case 2. In Figure 16(b), maximum disparity reaches 80%. This finding indicates that using the Desai element without damping can overestimate the bending moment response of the pile. The ratio of the peak bending moment at the bottom of the pile is meaningless because all the peak bending moments of cases 1, 2, and 5 are very small.

Conclusion

Seismic response analyses for a 3D finite element model were performed to simulate a soil–pile–structure system with an improved Desai interface element during an earthquake. The numerical results showed that the peak inner force of the frame column of the model with improved Desai element can be approximately10% lower than that of the model with the original Desai element. The pile peak shear force of the model with the improved Desai element is at least 19.2% lower than that of the model with the original Desai element. With regard to the pile bending moment, maximum disparity reached 80%. Using the Desai element without damping would overestimate the bending moment response of the pile.