# Combining Finite Element and Analytical methods to Contact Problems of 3D Structure on Soft Foundation

Authors: Chao Su; Heng Zhang; Shaopei Hu; Jiawei Bai; Jianjian Dai
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The computational efficiency and nonconvergence of the iteration are two main difficulties in contact problems, especially in the creep of the foundation. This paper proposes a method to analyze the structural soft foundation system affected by time. The methodology is an explicit method, combining the finite element method with the analytical method. The creep deformation of soft foundation is obtained based on Laplace transforms. The structural deformation contains the statically determinate structural deformation and rigid body displacement, solved by the finite method. The contact forces are calculated by the deformation coordination equations and equilibrium equations. The methodology is validated through augmented Lagrangian (AL) method. The results can clearly illustrate the local disengagement phenome, greatly overcome the nonconvergence of the iteration, and significantly improve computing efficiency.


---

## Body

## 1. Introduction

Hydraulic structures in plain areas are mostly built on soft soil which have high water content, low strength, and low bearing capacity. It often needs taking measures to strengthen the foundation during the construction period. The soil layer bottom of the foundation is the bottleneck controlling the bearing capacity and also the main place where the creep occurs [1]. The creep property results in the significant alteration of the foundation deformation, which leads to the uneven settlement in the upper structure [2]. In contrast to traditional problem of beam or plate on an elastic foundation, the influence of creep can be recognized as a necessary long-term stability method to analyze stress state of the structure [3]. Long-term monitoring data show that the creep deformation continues to increase during the operation period, and it can reach more than twice at the end of construction period [4]. The creep creates unique contact changing laws of the soft foundation and is cognized as the most serious factor that directly affects stress redistribution of the structure. In addition, the computational efficiency and calculation accuracy of the contact forces also face great challenge.In structural soft foundation system, the contact surfaces are deterministic and searching of contact surfaces is not necessary. Contact forces are time-sensitive and spatiality-nonlinear. Therefore, the contact problem can be simplified to finding the solution of contact forces, which mainly include analytical method, direct iterative method, mathematical programming method, penalty method, augmented Lagrangian method, and contact element method. In 1881, Hertz obtained an analytical solution in two contact bodies [5]. Signorini added the general formulation and defined it as a unilateral contact problem [6]. Winkler established the liner relationship between contact forces and deflection deformation, and got the analytical solution of contact problem in beam on elastic foundation [7]. Analytical method is still studied today and is widely used to solve the problem of dynamic load in rail transportation [8, 9]. The direct iterative method is numerical. For typical example, Francavilla obtained flexibility matrices in terms of contact forces at possible contact surfaces of two bodies and solved the quasilinear problem [10]. The direct iterative method has a clear concept, but the computational complexity is heavy owing to large quantities of the possible contact forces. In mathematical programming method, it is regarded as linear complementarity problem (LCP) or parametric linear complementarity problem (PLCP). Anders and Gunnar formulated the contact problems including varying contact surfaces and friction through the mathematical programming method [11, 12]. The advantage of the penalty method is imposing the contact conditions without increasing the number of variables. A penalty term is used to enforce the contact constraints. Luenberger discussed the ill-conditioned numerical problem caused by large values for penalty parameters [13]. Augmented Lagrangian method combines the penalty method and the Lagrangian method. It transforms the problem of original constraints into optimization so that it is well suited to the method of finite element and widely used in related software. Simo proposed the technique updating Lagrangian multipliers with penalty parameters to inherit the advantage of Lagrangian method [14]. Contact element method describes the contact behavior in interfaces. It expresses the contact stresses regarded as a function of the relative displacements in the mean planes of the microscopically rough surfaces. In another word, it assumes that the interface is a kind of element type and has constitutive laws. Goodman proposed the four-node planar and nonthickness contact surface element [15].A series of studies in contact problems of beam on elastic foundation make excellent contributions to long-span engineering programs [8]. It assumes that the reactive force of the foundation carrying a loaded beam at every point is proportional to the corresponding deflection of the beam [7]. Higher-order nonlinear partial differential equation can be used to reflect the contact phenomena [16, 17]. Through these researches, Gao considered the beam and foundation as two different element types, and each can use an own type of finite element [18]. This idea inspires us to study further. The key to contact problem lies in finding the solution of unilateral problems. The structure and soft foundation are two different contacting bodies, and each can be solved independently. Therefore, when the time-affection relationship of soft foundation is established, the same of the contact problems can be solved.The creep deformation of soft foundation has received extensive attention. Any foundation needs to quantitatively descript soil structure, rheological characteristic, spatial distribution, and mechanical property, evaluating accurately the bearing capacity, the effective stress, and settlement deformation. Terzaghi and Biot derived the three-dimensional consolidation equations based on the principle of effective stress [19, 20]. These equations have been well applied in complex foundation so far [21–23]. Singh adopted creep rate and creep function to construct a famous empirical creep model [24]. When the characterization of deformation is obtained, it is necessary to develop viscoelastic or elastic-viscoplastic constitutive theories based on micromechanics. The component models are applied to describe the rheological property and comprised of varying Hooke units and Kelvin units, such as the Kelvin model, the Maxwell model, and the Burgers model [25, 26]. Some elastic-viscoplastic constitutive models are built based on the Can-Clay creep model and modified Can-Clay creep model [27, 28]. Furthermore, Lee proposed an analytical method to solve the viscoelastic deformation by Laplace transform [29]. This analytical method has few parameters and can be easily used in soft foundation. Therefore, this method can be used to describe the contact forces and deformation of the soft foundation in the unilateral constraint problem.In this paper, the contact problem is considered as a unilateral constraint problem. This paper introduces a methodology to analyze the time affection on structural soft foundation system. The methodology is presented as a numerical implementation and combining finite element method with analytical method. It aims to solve the contact problem and focus on the contact forces varying with time. In this methodology, the finite element method is used to solve the structure deformation. The analytical method is based on the Maxwell model and Lee’s theory, extended to three dimensions and applied to the soft foundation. Then, the deformation coordination equation and force-method asymmetric matrix equation for quasilinear problems of the contact surface can be established. The contact situations, contact forces, and deformation at each time increment step are determined finally.


## 2. Analysis of the Problem Formulations

Similar to the direct method, contact forces are split into forces and reaction forces in the structural soft foundation system [7]. Therefore, the contact problems of the structural soft foundation system are divided into two unilateral constraint problems. As Figure 1 shows, contact forces are treated as unknown variables and act on the structure and foundation separately. The deformation coordination equations of contact surface can be put on solving the contact forces. Then, these contact forces can act on the structure, and the stress state of the structure is obtained.Figure 1
Structural soft foundation system. (a) Structure and soft foundation. (b) Load combination.
(a)(b)The deformation of structure contains the statically determinate structural deformation and rigid body displacement. The force method based on the Boltzmann superposition principle is used to solve the statically determinate structural deformation. Firstly, the basic deformation caused by each unit force of contact surface is calculated. Secondly, the deformation caused by external forces of contact surface is calculated. Thirdly, the statically determinate structural deformation is shown as(1)Ust=Uf0+δ1×X1+⋯+δn×Xn,where Ust is the deformation obtained by all forces in the statically determinate structural system; Uf0 is the deformation obtained by the external forces; δ1…δn. are the deformation obtained by each unit force separately; and X1…Xn are the unknown contact forces.The rigid body displacement of the structure can be shown by the deformation of six supports. Through the Lee’s method, the deformation of contact surface can also be shown byn contact force on the soft foundation multiplied by the deformation caused by the corresponding unit force. A total of n + 6 variables are formed. It still needs to add the force and bending moment equilibrium equations. Finally, the contact surface equations include n deformation coordination equations and 6 equilibrium equations. It can be simplified as follows:(2)Ust+Urb=Usf,∑Fx=0;∑Fy=0;∑Fz=0,∑Mx=0;∑My=0;∑Mz=0,where Urb is the rigid body displacement of structural contact surface; Ust is the contact surface deformation of soft foundation; and ∑F=0∑M=0 are composited by the forces and bending moments in three directions. The deformation coordination equations of the contact surface are built in the first line in (2). The force and bending moment equilibrium equations in three directions are added. The contact forces and rigid body displacement are calculated. Then, the stress state of the structure can be calculated.


## 3. Analytical Solution of Soft Foundation

The analytical solution of viscoelastic deformation is obtained by Lee’s method. Taking Kelvin model as an example, the elastic solution of a semi-infinite space body subjected to concentrated load is integrated to derive the solution of distributed force based on the principle of superposition. Then, the viscoelastic deformation in the Laplace space is obtained based on the elastic-viscoelastic correspondence principle. Finally, the solution of the viscoelastic deformation under distributed force can be obtained by inverse Laplace transformation.As shown in Figure2, there is a rectangle with length a and width b on the boundary of the semi-infinite foundation. The normal uniform load (Figure 2(a)) and tangential uniform load (Figure 2(b)) are concentration of 1/ab. They act on the rectangle. The elastic deformation acted normal concentrated force is shown:(3)uzz=1+μP2πER21−μ+z2R2,uzr=1+μP2πERrzR2−1−2μrR+z,where uzr is the radial deformation; uzz is the normal deformation; P is the normal concentrated force; R is the distance from a point to the origin of coordinates; r is the distance from a point to the normal line; E is the elasticity modulus; and μ is Poisson’s ratio. (3) has no solution when R = 0 so that unit distribution force is considered. The normal differential force dP shown based on Lagrangian coordinate system is dηdξ/ab. The elastic deformation acted normal concentrated force is shown:(4)δzz=1−μ2abπE∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdηξ2+η2,δzx=1+μ2μ−1abπE∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdηξ2+η2xx2+y2,δzy=1+μ2μ−1abπE∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdηξ2+η2yx2+y2,where δzx, δzy, δzz are the deformation in x, y, z directions under unit distribution force in z direction. According to elastic-viscoelastic correspondence principle, after Laplace transforms, the viscoelastic deformation formula is obtained:(5)δzzs=1−μs2abπEse−tssFk,δzxs=1+μs2μs−1abπEse−tssxx2+y2Fk,δzys=1+μs2μs−1abπEse−tssyx2+y2Fk,where Fk=∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdη/ξ2+η2 is obtained by the numerical integration method.Figure 2
Normal and tangential uniform loads act on space semi-infinite foundation. (a) Normal uniform load. (b) Tangential uniform load.
(a)(b)Kelvin viscoelastic solution is derived by inverse Laplace transformation.(6)δzz=36K+2G11−e−6K+2G1/2η1+12G11−e−G1/η1Fk2abπ,δzx=−36K+2G11−e−6K+2G1/2η1Fk2abπxx2+y2,δzy=−36K+2G11−e−6K+2G1/2η1Fk2abπyx2+y2,where K is the bulk modulus; G1 is the shear stiffness; and η1 is the coefficient of viscosity. Tangential load causes the deformation; the elastic deformation is shown in the following equation [30]:(7)uxx=1+μP2πER1+x2R2+1−2μRR+z−x2R+z2,uxy=1+μP2πERxyR2−1−2μxyR+z2,uxz=1+μP2πERxyR2+1−2μxR+z.Through the Laplace transformation and inverse Laplace transformation, the viscoelastic deformation can also be solved and shown in the following equation:(8)δxx=12πab12G1e−G1/η1t+36K+2G1e−6K+2G1/2η1tFk+12πab12G1e−G1/η1t−36K+2G1e−6K+2G1/2η1tx2Fg,δxy=xyFg2πab12G1e−G1/η1t−36K+2G1e−6K+2G1/2η1t,δxz=xFi2πab36K+2G1e−6K+2G1/2η1t.Fi=∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdη/ξ2+η2Fg=∫η=x−a/2η=x+a/2∫ξ=y−b/2ξ=y+b/2dξdη/ξ2+η23/2 is obtained by the numerical integration method, where δxx, δxy, δxz are the deformation in x, y, z directions under unit distribution force in x direction. The analytical solution of viscoelastic deformation based on the Maxwell model can be calculated by MATLAB. The process is given in the appendix. The final equation is shown as follows:(9)δxx=4G3Fi3KFi−2G3Fgx2+3KFgx24G3abπG3+3K−3η1etG2G3+3G2K/G2η1+G3η1+3Kη1G2Fi−G2Fgx24abπG3+3KG2η1+G3η1+3Kη1−η1eG2G3t/G2η1+G3η1G2Fgx2+G2Fi4G3abπG2η1+G3η1,δxy=3G2η1FgxyetG2G3+3G2K/G2η1+G3η1+3Kη14abπG3+3KG2η1+G3η1+3Kη1−Fgxy2G3−3K4G3abπG3+3K−G2η1FgxyeG2G3t/G2η1+G3η14G3abπG2η1+G3η1,δxz=3Fkx4abπG3+3K_3G2η1FgxyetG2G3+3G2K/G2η1+G3η1+3Kη14abπG3+3KG2η1+G3η1+3Kη1,δzx=3G2η1FgxyetG2G3+3G2K/G2η1+G3η1+3Kη14abπG3+3KG2η1+G3η1+3Kη1−3Fix4abπx2+y2G3+3K,δzy=3G2η1FgxyetG2G3+3G2K/G2η1+G3η1+3Kη14abπx2+y2G3+3KG2η1+G3η1+3Kη1−3Fiy4abπx2+y2G3+3K,δzz=Fi4G3+3K4G3abπG3+3K−3G2η1FgxyetG2G3+3G2K/G2η1+G3η1+3Kη14abπx2+y2G3+3KG2η1+G3η1+3Kη1−G2η1FkeG2G3t/G2η1+G3η14abπG2η1+G3η1,where G2, G3 are the parameters of shear stiffness. The deformation at any point of the contact surface on soft foundation acted on a force can be shown as(10)UsfxUsfyUsfz=δxxδyxδzxδxyδyyδzyδxzδyzδzzFxFyFz,where Usfx,Usfy,andUsfz are the deformation of soft foundation in three directions. Fx,Fy,andFz are the contact forces of soft foundation in three directions. (10) shows the deformation of a single point under a distribution force. The real contact surface deformation is more complicated.


## 4. Structure Deformation

### 4.1. Statically Determinate Structural Deformation

The contact forces are treated as unknown variables and act on the structure and foundation separately. As Figure3 shows, in order to ensure the structure statically determinate, the structure has six constraints on the contact surface. It is divided into multiple hexahedrons. Each unit distribution force acts on the bottom of underside element. Note that the unit distribution force acting on the contact surface separating two adjacent volume elements in the structure and the soft foundation must be equal and opposite. The resultant point of distributed force is under the bottom of element. The deformation of the resultant point can be shown by deformation of element nodes. The equation is as follows:(11)uvw=N100…N4000N10...0N4000N1…00N4qe,where N1,N2,N3,N4 are the shape function of the four nodes which are at the bottom of the hexahedron element and qe is the deformation of element. The programming of these deformations can be finished with DLOAD subroutine and UTRACLOAD subroutine in Abaqus. The DLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in z-axis direction. The UTRACLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in x and y direction. The xleft, xright, yleft, and yright are the four boundaries of the unit distributed force. The xstart, xend, ystart, and yend are the first and final coordinates of the unit distributed force. The T_user(1), T_user(2), and T_user(3) are the direction of the unit distribution force. The v and time are the speed and time of the unit distributed force. The coord (1) and coord (2) are the x and y coordinates of each position of the contact surface. It determines where the unit distributed force is applied and where the force is zero. Figure 4 shows the flow chart of subroutines in Abaqus. (1) All parameters are input in the subroutine. (2) The parameters of xleft, xright, yleft, and yright are updated. (3) The position of force is determined. (4) The time is updated and the above process is repeated.Figure 3
Statically determinate structure deformation.Figure 4
Flow chart of subroutines in Abaqus. (a) UTRACLOAD subroutine. (b) DLOAD subroutine.
(a)(b)A unit distributed force moves from the bottom of the first element to the bottom of the last element on the underside element. The deformation during the movement is recorded. According to equation (11), the deformation of the resultant point is calculated by numerical software. Then, the statically determinate structural deformation can also be expressed in the same form as (10).


## 4.1. Statically Determinate Structural Deformation

The contact forces are treated as unknown variables and act on the structure and foundation separately. As Figure3 shows, in order to ensure the structure statically determinate, the structure has six constraints on the contact surface. It is divided into multiple hexahedrons. Each unit distribution force acts on the bottom of underside element. Note that the unit distribution force acting on the contact surface separating two adjacent volume elements in the structure and the soft foundation must be equal and opposite. The resultant point of distributed force is under the bottom of element. The deformation of the resultant point can be shown by deformation of element nodes. The equation is as follows:(11)uvw=N100…N4000N10...0N4000N1…00N4qe,where N1,N2,N3,N4 are the shape function of the four nodes which are at the bottom of the hexahedron element and qe is the deformation of element. The programming of these deformations can be finished with DLOAD subroutine and UTRACLOAD subroutine in Abaqus. The DLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in z-axis direction. The UTRACLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in x and y direction. The xleft, xright, yleft, and yright are the four boundaries of the unit distributed force. The xstart, xend, ystart, and yend are the first and final coordinates of the unit distributed force. The T_user(1), T_user(2), and T_user(3) are the direction of the unit distribution force. The v and time are the speed and time of the unit distributed force. The coord (1) and coord (2) are the x and y coordinates of each position of the contact surface. It determines where the unit distributed force is applied and where the force is zero. Figure 4 shows the flow chart of subroutines in Abaqus. (1) All parameters are input in the subroutine. (2) The parameters of xleft, xright, yleft, and yright are updated. (3) The position of force is determined. (4) The time is updated and the above process is repeated.Figure 3
Statically determinate structure deformation.Figure 4
Flow chart of subroutines in Abaqus. (a) UTRACLOAD subroutine. (b) DLOAD subroutine.
(a)(b)A unit distributed force moves from the bottom of the first element to the bottom of the last element on the underside element. The deformation during the movement is recorded. According to equation (11), the deformation of the resultant point is calculated by numerical software. Then, the statically determinate structural deformation can also be expressed in the same form as (10).


## 5. Rigid Body Displacement

The rigid body displacement of statically determinate structure can be described by Cauchy equations. It means that the strain of the structure is zero. For infinitesimal motion, the relationship between strain and displacement is(12)εx=∂u∂x,εy=∂v∂y,εz=∂w∂z,γyz=∂w∂y+∂v∂z,γzx=∂u∂z+∂w∂x,γxy=∂v∂x+∂u∂y,εx=εy=εz=γyz=γzx=γxy=0,where εx,εy,εz are the normal strain in x, y, z directions and γyz,γzx,γxy are the shearing strain in x, y, z directions. Solving the above formulas, the rigid body displacement can be obtained as(13)u=u0+φyz−φzy,v=v0+φzx−φxz,w=w0+φxy−φyx,where u, v, w are the rigid body displacement in x, y, z directions; φx,φy,φz are the rigid body rotation angles in x, y, z directions; and u0, v0, w0 are the translational deformations in x, y, z directions. These equations can be comprehended through geometric transformation in Figure 5. It is important that the sign of rigid body rotation angles follows the right hand’s spiral rule. When z = 0, the equation shows the deformation of undersurface. (13) can be simplified as linear matrix equation:(14)uvw=10000−y01000x001yx0u0v0w0φxφyφz.Figure 5
Rigid body displacement of contact surface.After the statically determinate structural deformation and rigid body displacement are obtained, the contact surface deformation can be described as the sum of them in the structure.


## 6. Force Equilibrium Equations and Bending Moment Equilibrium Equations

By means of rational mechanics, the force equilibrium equations and bending moment equilibrium equations of 3D structure can be elucidated in three directions. These equations are as follows:(15)∑k=1nXk+Fx=0,∑k=1nYk+Fy=0,∑k=1nZk+Fz=0,∑k=1nbkZk−∑k=1nckYk+yFFz−zFFy=0,∑k=1nckXk−∑k=1nakZk+zFFx−xFFz=0,∑k=1nakYk−∑k=1nbkXk+xFFy−yFFx=0,where ∑k=1nXk,∑k=1nYk,∑k=1nZk are the sum of contact forces in x, y, z directions; ak,bk,ck are the coordinates of the contact forces’ functional point; xF,yF,zF are the coordinates of the external forces’ functional point; and Fx,Fy,Fz are the external forces in x, y, z directions. When c = 0, (15) can be simplified as linear matrix equation:(16)1...10...00...00...01...10...00...00...01...10...00...0b1...bn0...00...0−a1...−an−b1...−bna1...an0...0X1...XnY1...YnZ1...Zn=−Fx−Fy−Fz−yFFz+zFFy−zFFx+xFFz−xFFy+yFFx.


## 7. Mixed Finite Element Methodology

Through the analysis ofSections 3, the time affection of soft foundation can be resolved. When the parameters of foundation are determined, the contact forces can be discussed. This contact problem is solved based on the deformation coordination equations.The eight-node isoparametric elements are used to make the structure discrete. Three connecting rods at the bottom center point of the underside element are set along thex, y, and z directions. They are used to connect the structure with foundation. The increment of the normal deformation Δz and the increment of the rigid body rotation angles ΔφxΔφy are taken as unknown quantities. Taking time step l for example, the finite element method is used to obtain the deformation of nodes and the deflection of the bottom center point of the underside elements. The deformation and deflection are caused by the external loads and unit link force. The normal incremental equation of structure on viscoelastic foundation is solved by combining finite element with the analytical method. It is shown as(17)∑k=1nδkilΔZkl+Δz+biΔφx−aiΔφy+Δipl=0,nformulas in total,i=1,2,…,n,∑k=1nΔZkl+∑ΔFzl=0,∑k=1nbkΔZkl+∑ΔMxl=0,∑k=1nakΔZkl+∑ΔMyl=0,where δkil is the deformation of i position under unit distribution force acting on k position in z direction in the time step l; Δipl is the deformation of i position under external forces at z direction in the time step l; ai is the coordinate of y-axis at i position in the time step l; bi is the coordinate of x-axis at i position in the time step l; ΔZkl is the element connecting rod force at z direction in the time step l; ΔFzl is the composition of forces at z direction in the time step l; ΔMxl,ΔMyl are the composition of bending moments at x and y directions in the time step l. The basic unknown quantities can be solved by the Gaussian elimination with partial pivoting method. Then, the increment of each unknown quantity of the system in the time step is obtained. The equations of contact surface in three-dimensional are also shown as(18)∑k=1nδxxkilΔXkl+δyxkilΔYkl+δzxkilΔZkl+Δxl−biΔφzl+Δipxl=0,∑k=1nδxykilΔXkl+δyykilΔYkl+δzykilΔZkl+Δyl−aiΔφzl+Δipyl=0,∑k=1nδxzkilΔXkl+δyzkilΔYkl+δzzkilΔZkl+Δzl+aiΔφxl−biΔφyl+Δipzl=0,3nformulas in total,i=1,2,…,n,∑k=1nΔXkl+∑ΔFxl=0,∑k=1nΔYkl+∑ΔFyl=0,∑k=1nΔZkl+∑ΔFzl=0,∑k=1nakΔZkl+∑ΔMxl=0,∑k=1nbkΔZkl+∑ΔMyl=0,∑k=1nbkΔXkl+∑k=1nakΔYkl+∑ΔMzl=0,δmkil=δstmkil+δsfmkil,m=xx,yx,zx,xy,yy,zy,xz,yz,zz,where δmkil is the deformation of i position under unit distribution force acting on k position at m-axis direction in the time step l; δsfmkil,δstmkil are the deformation of foundation and structure under unit distribution force; ΔXkl,ΔYkl,ΔZkl are the contact forces at k position in x, y, and z directions in the time step l.The program flow chart is shown in Figure6. (1) The structure is meshed with hexahedral elements and the position of each connecting rod on the contact surface is set. (2) The analytical solution for soft foundation of each unit distribution force and FEM solution for the structure are obtained. (3) The coefficient matrix and typical equation are formed. (4) The increment of contact forces and foundation deformation are solved. (5) The contact state is judged and the redistribution of contact forces are calculated.Figure 6
Program flow chart.


## 8. Model Validation

### 8.1. Example 1

The model validation can be finished through comparative analysis of our method and the augmented Lagrangian (AL) method. The augmented Lagrangian (AL) method is calculated by Abaqus. The time affection of soft foundation is described by Maxwell UMAT [31]. The length, width, and height of structure is separately 10 m, 6 m, and 3 m. The pressure acting on the upper surface in z direction is 1000 Pa. The pressure acting on the left surface in x direction is 500 Pa. The total time is 1200 d. The time increment is 30 d. The structure can be dispersed into multiple hexahedrons. The labels of the bottom center point of the underside element on the contact surface are 1, 2,...,240. Each label has three degrees of freedom (DOF). There is a total of 720 contact forces in three DOF. From (18), the contact forces equations are shown at time step l:(19)δxx,1,1…δxx,240,1δyx,1,1…δyx,240,1δzx,1,1…δzx,240,110000−b|x=0.25y=0.25……………………………………δxx,1,240…δxx,240,240δyx,1,240…δyx,240,240δzx,1,240…δzx,240,24010000−b|x=9.75y=5.75δxy,1,1…δxy,240,1δyy,1,1…δyx,240,1δzy,1,1…δzy,240,101000a|x=0.25y=0.25………………………………………δxy,1,240…δxy,240,240δyy,1,240…δyy,240,240δzy,1,240…δzy,1,24001000a|x=9.75y=5.75δxz,1,1…δxz,240,1δyz,1,1…δyz,240,1δzz,1,1…δzz,240,1001b|x=9.75y=5.75a|x=9.75y=5.750………………………………………δxz,1,240…δxz,240,240δyz,1,240…δyz,240,240δyz,1,240…δzz,240,240001b|x=9.75y=2.75−a|x=9.75y=2.7501…10…00…00000000…01…10…00000000…00…01…00000000…00…0a|x=0.25y=0.25…00000000…00…0−a|x=9.75y=2.75…0000000−b|x=0.25y=0.25…−b|x=9.75y=5.75a|x=0.25y=0.25…a|x=9.75y=2.750…0000000X1…X240Y1…Y240Z1…Z240l1l2l3l4l5l6=−Δ1,x…−Δ240,x−Δ1,y…−Δ240,yZ1…−Δ240,z−Fx−Fy−Fz−yFFzxFFz−xFFy+yFFx,(20)ABC0FG=DE,where a, b is the coordinate of center point in x and y directions, respectively. The matrices in (19) can be portioned into seven blocks in (20) for programming, as shown in (20). The block A stands for the sum of the deformation of structure and foundation resulted by unit distribution force. The block B∗G stands for the rigid body displacement. The block C∗F stands for the external forces and bending moment. The parameters of structure and foundation are shown in Table 1.Table 1
The parameters of structure and foundation.
StructureSoft foundationE3.425∗1010 PaE11.7546∗107 Paμ0.167E21.552∗107 Pa——μ0.31——η18.139∗108 Pa∗d——K1.5391∗107 Pa——G26.6969∗106 Pa——G35.9237∗106 PaFigure7 shows the deformation of soft foundation in three directions under unit distribution force at 30th day and 1200th day. The deformation results of δsfxy and δsfyx are center symmetrical about the point of force. The deformation results of δsfxx,δsfxz,δsfyy,δsfyz,δsfzx,δsfzy, and δsfzz are symmetrical about one axis. The influential sphere of deformation caused by the unit distribution force is restricted from 0 to 3.5 m. This is mainly attributed to the elastic modulus and Poisson’s ratio of the soft foundation. Maximum deformation ranges from 7.93e-8m at the first step to 1.26e-7m at the final step. It shows a 59.52% increment within the time. Overall, the Lee’s method can be well used to describe the creep property of soft foundation.Figure 7
The deformation of soft foundation in three directions under unit distribution force. (a) On 30th day. (b) On 1200th day.
(a)(b)Figure8 shows the contact forces in three directions by two methods on the 30th day. Figure 9 shows the contact forces in three directions by two methods on the 1200th day. Figure 10 shows the deformation of the center point in the contact surface. Through the contrastive analysis of our method and AL method, some similarities and differences can be shown:(1)
From the distribution of contact forces perspective, both of the two methods almost have the same distribution of contact forces. The contact forces inx direction are fan-shaped distribution and focused on the boundary. The contact forces in y direction are symmetrical about x-axis. The contact forces in z direction are basin-shaped distribution. Our method has larger values of corner points than the AL method. The reason behind this scenario may be that the deformation of soft foundation in Lee’s method is different with finite method. Overall, contact forces can be well solved by our method.(2)
From the creep property perspective, our method shows the redistribution of the contact forces. The maximum contact force inz direction increases from 796N to 868N, with 9.01% increment. However, the maximum contact forces in x and y directions, respectively, decrease from 464N to 410N and 286N to 236N, with 13.17% and 21.19% reduction.(3)
From deformations perspective, the displacements of our method are more viscous than the AL method with time changing. Both of them almost have the same final settlement displacement and displacement inx direction.Figure 8
Distribution of contact forces in three directions by two methods on 30th day. (a) Our method. (b) Augmented Lagrangian method.
(a)(b)Figure 9
Distribution of contact forces in three directions by two methods on 1200th day. (a) Our method. (b) Augmented Lagrangian method.
(a)(b)Figure 10
Deformation of the center point in the contact surface.The difference between our method and the AL method is mainly influenced by the solution of foundation. The AL method is based on the finite element method and related to the size of the foundation model. The analytical method is based on the Laplace transformation and irrelevant to the size of the foundation model. Therefore, the results are influenced by two reasons: (1) the principle of calculation; (2) the size of foundation model. The heavy loading combination and large size of the foundation model can make these differences not obvious. In this example, the external loads are inx and z directions and the length, width, and height of the foundation are three times of the structure. The results show good agreement in the final settlement displacement and displacement in x and z directions with the AL method and proposed method, while there exists large difference in y direction.


### 8.2. Example 2

Example 1 shows that our method has well adaptability. The local disengagement of contact surface is shown in example 2. The basic situation of this example is the same as example 1, but the external forces are different with example 1. The upper surface of the structure is affected by pressure and tension. The pressure acting on the left square (4.5 m∗6 m) and right square (4.5 m∗6 m) is 10 kPa. The tension acting on the middle square (1.0 m∗6m) is 8.9 kPa.Figure11 shows the contact forces in three directions on the 30th and the 1200th day. The contact forces in z direction are directly affected by the external forces. It represents that there are negative contact forces in the middle of the contact surface. It means the local disengagement phenomenon appears. Actually, these negative contact forces in the disengagement area are zeros. By comparing the positive with negative values, the disengagement area can be intuitively found. After the contact forces redistribute, the real contact state and contact forces are explicitly calculated. The total time of the process is usually less than 10 seconds. The deformation of the center point in the contact surface is shown in Figure 12, where the center point is in the disengagement area. The displacement of the point in z direction has creep property. As time increases, the disengagement area will gradually close, and the contact forces will regenerate. The disengagement area does not always exist.Figure 11
Contact forces in three directions. (a) Contact forces at 30th d. (b) Contact forces at 1200th d.
(a)(b)Figure 12
Deformation of the center point in the contact surface.


## 8.1. Example 1

The model validation can be finished through comparative analysis of our method and the augmented Lagrangian (AL) method. The augmented Lagrangian (AL) method is calculated by Abaqus. The time affection of soft foundation is described by Maxwell UMAT [31]. The length, width, and height of structure is separately 10 m, 6 m, and 3 m. The pressure acting on the upper surface in z direction is 1000 Pa. The pressure acting on the left surface in x direction is 500 Pa. The total time is 1200 d. The time increment is 30 d. The structure can be dispersed into multiple hexahedrons. The labels of the bottom center point of the underside element on the contact surface are 1, 2,...,240. Each label has three degrees of freedom (DOF). There is a total of 720 contact forces in three DOF. From (18), the contact forces equations are shown at time step l:(19)δxx,1,1…δxx,240,1δyx,1,1…δyx,240,1δzx,1,1…δzx,240,110000−b|x=0.25y=0.25……………………………………δxx,1,240…δxx,240,240δyx,1,240…δyx,240,240δzx,1,240…δzx,240,24010000−b|x=9.75y=5.75δxy,1,1…δxy,240,1δyy,1,1…δyx,240,1δzy,1,1…δzy,240,101000a|x=0.25y=0.25………………………………………δxy,1,240…δxy,240,240δyy,1,240…δyy,240,240δzy,1,240…δzy,1,24001000a|x=9.75y=5.75δxz,1,1…δxz,240,1δyz,1,1…δyz,240,1δzz,1,1…δzz,240,1001b|x=9.75y=5.75a|x=9.75y=5.750………………………………………δxz,1,240…δxz,240,240δyz,1,240…δyz,240,240δyz,1,240…δzz,240,240001b|x=9.75y=2.75−a|x=9.75y=2.7501…10…00…00000000…01…10…00000000…00…01…00000000…00…0a|x=0.25y=0.25…00000000…00…0−a|x=9.75y=2.75…0000000−b|x=0.25y=0.25…−b|x=9.75y=5.75a|x=0.25y=0.25…a|x=9.75y=2.750…0000000X1…X240Y1…Y240Z1…Z240l1l2l3l4l5l6=−Δ1,x…−Δ240,x−Δ1,y…−Δ240,yZ1…−Δ240,z−Fx−Fy−Fz−yFFzxFFz−xFFy+yFFx,(20)ABC0FG=DE,where a, b is the coordinate of center point in x and y directions, respectively. The matrices in (19) can be portioned into seven blocks in (20) for programming, as shown in (20). The block A stands for the sum of the deformation of structure and foundation resulted by unit distribution force. The block B∗G stands for the rigid body displacement. The block C∗F stands for the external forces and bending moment. The parameters of structure and foundation are shown in Table 1.Table 1
The parameters of structure and foundation.
StructureSoft foundationE3.425∗1010 PaE11.7546∗107 Paμ0.167E21.552∗107 Pa——μ0.31——η18.139∗108 Pa∗d——K1.5391∗107 Pa——G26.6969∗106 Pa——G35.9237∗106 PaFigure7 shows the deformation of soft foundation in three directions under unit distribution force at 30th day and 1200th day. The deformation results of δsfxy and δsfyx are center symmetrical about the point of force. The deformation results of δsfxx,δsfxz,δsfyy,δsfyz,δsfzx,δsfzy, and δsfzz are symmetrical about one axis. The influential sphere of deformation caused by the unit distribution force is restricted from 0 to 3.5 m. This is mainly attributed to the elastic modulus and Poisson’s ratio of the soft foundation. Maximum deformation ranges from 7.93e-8m at the first step to 1.26e-7m at the final step. It shows a 59.52% increment within the time. Overall, the Lee’s method can be well used to describe the creep property of soft foundation.Figure 7
The deformation of soft foundation in three directions under unit distribution force. (a) On 30th day. (b) On 1200th day.
(a)(b)Figure8 shows the contact forces in three directions by two methods on the 30th day. Figure 9 shows the contact forces in three directions by two methods on the 1200th day. Figure 10 shows the deformation of the center point in the contact surface. Through the contrastive analysis of our method and AL method, some similarities and differences can be shown:(1)
From the distribution of contact forces perspective, both of the two methods almost have the same distribution of contact forces. The contact forces inx direction are fan-shaped distribution and focused on the boundary. The contact forces in y direction are symmetrical about x-axis. The contact forces in z direction are basin-shaped distribution. Our method has larger values of corner points than the AL method. The reason behind this scenario may be that the deformation of soft foundation in Lee’s method is different with finite method. Overall, contact forces can be well solved by our method.(2)
From the creep property perspective, our method shows the redistribution of the contact forces. The maximum contact force inz direction increases from 796N to 868N, with 9.01% increment. However, the maximum contact forces in x and y directions, respectively, decrease from 464N to 410N and 286N to 236N, with 13.17% and 21.19% reduction.(3)
From deformations perspective, the displacements of our method are more viscous than the AL method with time changing. Both of them almost have the same final settlement displacement and displacement inx direction.Figure 8
Distribution of contact forces in three directions by two methods on 30th day. (a) Our method. (b) Augmented Lagrangian method.
(a)(b)Figure 9
Distribution of contact forces in three directions by two methods on 1200th day. (a) Our method. (b) Augmented Lagrangian method.
(a)(b)Figure 10
Deformation of the center point in the contact surface.The difference between our method and the AL method is mainly influenced by the solution of foundation. The AL method is based on the finite element method and related to the size of the foundation model. The analytical method is based on the Laplace transformation and irrelevant to the size of the foundation model. Therefore, the results are influenced by two reasons: (1) the principle of calculation; (2) the size of foundation model. The heavy loading combination and large size of the foundation model can make these differences not obvious. In this example, the external loads are inx and z directions and the length, width, and height of the foundation are three times of the structure. The results show good agreement in the final settlement displacement and displacement in x and z directions with the AL method and proposed method, while there exists large difference in y direction.


## 8.2. Example 2

Example 1 shows that our method has well adaptability. The local disengagement of contact surface is shown in example 2. The basic situation of this example is the same as example 1, but the external forces are different with example 1. The upper surface of the structure is affected by pressure and tension. The pressure acting on the left square (4.5 m∗6 m) and right square (4.5 m∗6 m) is 10 kPa. The tension acting on the middle square (1.0 m∗6m) is 8.9 kPa.Figure11 shows the contact forces in three directions on the 30th and the 1200th day. The contact forces in z direction are directly affected by the external forces. It represents that there are negative contact forces in the middle of the contact surface. It means the local disengagement phenomenon appears. Actually, these negative contact forces in the disengagement area are zeros. By comparing the positive with negative values, the disengagement area can be intuitively found. After the contact forces redistribute, the real contact state and contact forces are explicitly calculated. The total time of the process is usually less than 10 seconds. The deformation of the center point in the contact surface is shown in Figure 12, where the center point is in the disengagement area. The displacement of the point in z direction has creep property. As time increases, the disengagement area will gradually close, and the contact forces will regenerate. The disengagement area does not always exist.Figure 11
Contact forces in three directions. (a) Contact forces at 30th d. (b) Contact forces at 1200th d.
(a)(b)Figure 12
Deformation of the center point in the contact surface.


## 9. Conclusion

This paper reports a methodology to analyze the time affection on a structural soft foundation system. The contact forces are treated as unknown variables and act on the structure and foundation separately. The deformation of structure contains the statically determinate structural deformation and rigid body displacement. It can be solved by the finite method. The creep deformation of soft foundation is solved by Lee’s method. The final equations are formed by the deformation coordination equations and equilibrium equations. It is solved by the Gaussian elimination with the partial pivoting method. The augmented Lagrangian method is used to verify the accuracy of the method. The following important conclusions are summarized:(1)
The simulations of contact bodies show that our method shows more viscous than the AL method with time changing. Our method has larger values of corner points than the AL method.(2)
This methodology is an explicit calculating method, avoiding the occurrence of iteration nonconvergence. It is simpler and more efficient in solving the contact forces compared with the AL method.(3)
It is shown that our method can directly describe local disengagement of two contact bodies. The redistribution phenomenon of contact forces is well shown with time changing.



---

*Source: 8827681-2020-10-30.xml*