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

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.


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. e soil layer bottom of the foundation is the bottleneck controlling the bearing capacity and also the main place where the creep occurs [1]. e 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]. e 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 spatialitynonlinear. erefore, 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]. e 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]. e 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]. e 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]. rough these researches, Gao considered the beam and foundation as two different element types, and each can use an own type of finite element [18]. is idea inspires us to study further.
e key to contact problem lies in finding the solution of unilateral problems. e structure and soft foundation are two different contacting bodies, and each can be solved independently.
erefore, when the time-affection relationship of soft foundation is established, the same of the contact problems can be solved. e 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]. ese equations have been well applied in complex foundation so far [21][22][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.
e 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].
is analytical method has few parameters and can be easily used in soft foundation. erefore, 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.
is paper introduces a methodology to analyze the time affection on structural soft foundation system. e 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. e analytical method is based on the Maxwell model and Lee's theory, extended to three dimensions and applied to the soft foundation. en, the deformation coordination equation and force-method asymmetric matrix equation for quasilinear problems of the contact surface can be established. e contact situations, contact forces, and deformation at each time increment step are determined finally.

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]. erefore, 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. e deformation coordination equations of contact surface can be put on solving the contact forces. en, these contact forces can act on the structure, and the stress state of the structure is obtained. e deformation of structure contains the statically determinate structural deformation and rigid body displacement. e 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. irdly, the statically determinate structural deformation is shown as where U st is the deformation obtained by all forces in the statically determinate structural system; U f0 is the deformation obtained by the external forces; δ 1 . . . δ n . are the deformation obtained by each unit force separately; and X 1 . . . X n are the unknown contact forces. e rigid body displacement of the structure can be shown by the deformation of six supports. rough the Lee's method, the deformation of contact surface can also be shown by n 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 2 Mathematical Problems in Engineering coordination equations and 6 equilibrium equations. It can be simplified as follows: where U rb is the rigid body displacement of structural contact surface; U st is the contact surface deformation of soft foundation; and F � 0 M � 0 are composited by the forces and bending moments in three directions. e deformation coordination equations of the contact surface are built in the first line in (2). e force and bending moment equilibrium equations in three directions are added. e contact forces and rigid body displacement are calculated. en, the stress state of the structure can be calculated.

Analytical Solution of Soft Foundation
e 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. en, 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 Figure 2, there is a rectangle with length a and width b on the boundary of the semi-infinite foundation. e normal uniform load (Figure 2(a)) and tangential uniform load (Figure 2(b)) are concentration of 1/ab. ey act on the rectangle. e elastic deformation acted normal concentrated force is shown: where u zr is the radial deformation; u zz 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. e normal differential force dP shown based on Lagrangian coordinate system is dηdξ/ab. e elastic deformation acted normal concentrated force is shown: 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:   Kelvin viscoelastic solution is derived by inverse Laplace transformation.
where K is the bulk modulus; G 1 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]: rough the Laplace transformation and inverse Laplace transformation, the viscoelastic deformation can also be solved and shown in the following equation: where G 2 , G 3 are the parameters of shear stiffness. e deformation at any point of the contact surface on soft foundation acted on a force can be shown as Mathematical Problems in Engineering where U sfx , U sfy , and U sfz are the deformation of soft foundation in three directions. F x , F y , and F z are the contact forces of soft foundation in three directions. (10) shows the deformation of a single point under a distribution force. e real contact surface deformation is more complicated.

Statically Determinate Structural Deformation.
e contact forces are treated as unknown variables and act on the structure and foundation separately. As Figure 3 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. e resultant point of distributed force is under the bottom of element. e deformation of the resultant point can be shown by deformation of element nodes. e equation is as follows: where N 1 , N 2 , N 3 , N 4 are the shape function of the four nodes which are at the bottom of the hexahedron element and q e is the deformation of element. e programming of these deformations can be finished with DLOAD subroutine and UTRACLOAD subroutine in Abaqus. e DLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in z-axis direction. e UTRACLOAD subroutine is used to get the solution of the deformation acted by each unit distribution force in x and y direction. e xleft, xright, yleft, and yright are the four boundaries of the unit distributed force. e xstart, xend, ystart, and yend are the first and final coordinates of the unit distributed force. e T_user(1), T_user(2), and T_user (3) are the direction of the unit distribution force. e v and time are the speed and time of the unit distributed force. e 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) e parameters of xleft, xright, yleft, and yright are updated. (3) e position of force is determined. (4) e time is updated and the above process is repeated.
A unit distributed force moves from the bottom of the first element to the bottom of the last element on the underside element. e deformation during the movement is recorded. According to equation (11), the deformation of the resultant point is calculated by numerical software. en, the statically determinate structural deformation can also be expressed in the same form as (10).

Rigid Body Displacement
e 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 where ε x , ε y , ε z are the normal strain in x, y, z directions and c yz , c zx , c xy are the shearing strain in x, y, z directions. Solving the above formulas, the rigid body displacement can be obtained as 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 u 0 , v 0 , w 0 are the translational deformations in x, y, z directions. ese 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: 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.

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. ese equations are as follows:  Mathematical Problems in Engineering where n k�1 X k , n k�1 Y k , n k�1 Z k are the sum of contact forces in x, y, z directions; a k , b k , c k are the coordinates of the contact forces' functional point; x F , y F , z F are the coordinates of the external forces' functional point; and F x , F y , F z are the external forces in x, y, z directions. When c � 0, (15) can be simplified as linear matrix equation: 1 . .. 1 0 ... 0 0 ... 0   0 ... 0 1 ... 1 0 ... 0   0 ... 0 0 ... 0 1 ... 1   0 ... 0 0 ... 0 b 1

Mixed Finite Element Methodology
rough the analysis of Sections 3, the time affection of soft foundation can be resolved. When the parameters of foundation are determined, the contact forces can be discussed.
is contact problem is solved based on the deformation coordination equations. e eight-node isoparametric elements are used to make the structure discrete. ree connecting rods at the bottom center point of the underside element are set along the x, y, and z directions. ey are used to connect the structure with foundation. e 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. e deformation and deflection are caused by the external loads and unit link force.
e normal incremental equation of structure on viscoelastic foundation is solved by combining finite element with the analytical method. It is shown as where δ (l) ki is the deformation of i position under unit distribution force acting on k position in z direction in the time step l; Δ (l) ip is the deformation of i position under external forces at z direction in the time step l; a i is the coordinate of y-axis at i position in the time step l; b i is the coordinate of xaxis at i position in the time step l; ΔZ (l) k is the element connecting rod force at z direction in the time step l; ΔF (l) z is the composition of forces at z direction in the time step l; ΔM (l) x , ΔM (l) y are the composition of bending moments at x and y directions in the time step l. e basic unknown quantities can be solved by the Gaussian elimination with partial pivoting method. en, the increment of each unknown quantity of the system in the time step is obtained. e equations of contact surface in three-dimensional are also shown as 8 Mathematical Problems in Engineering 3n formulas in total, i � 1, 2, . . . , n, n k�1  e augmented Lagrangian (AL) method is calculated by Abaqus. e time affection of soft foundation is described by Maxwell UMAT [31]. e length, width, and height of structure is separately 10 m, 6 m, and 3 m. e pressure acting on the upper surface in z direction is 1000 Pa. e pressure acting on the left surface in x direction is 500 Pa. e total time is 1200 d. e time increment is 30 d. e structure can be dispersed into multiple hexahedrons. e 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). ere is a total of 720 contact forces in three DOF. From (18), the contact forces equations are shown at time step l: where a, b is the coordinate of center point in x and y directions, respectively. e matrices in (19) can be portioned into seven blocks in (20) for programming, as shown in (20). e block A stands for the sum of the deformation of structure and foundation resulted by unit distribution force. e block B * G stands for the rigid body displacement. e block C * F stands for the external forces and bending moment. e parameters of structure and foundation are shown in Table 1. Figure 7 shows the deformation of soft foundation in three directions under unit distribution force at 30th day and 1200th day. e deformation results of δ sfxy and δ sfyx are  e deformation results of δ sfxx , δ sfxz , δ sfyy , δ sfyz , δ sfzx , δ sfzy , and δ sfzz are symmetrical about one axis. e influential sphere of deformation caused by the unit distribution force is restricted from 0 to 3.5 m. is 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 8 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.
rough 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. e contact forces in x direction are fan-shaped distribution and focused on the boundary. e contact forces in y direction are symmetrical about x-axis. e contact forces in z direction are basin-shaped distribution. Our method has larger values of corner points than the AL method. e 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. e maximum contact force in z 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. e difference between our method and the AL method is mainly influenced by the solution of foundation. e AL method is based on the finite element method and related to the size of the foundation model. e analytical method is based on the Laplace transformation and irrelevant to the size of the foundation model. erefore, the results are influenced by two reasons: (1) the principle of calculation; (2) the size of foundation model. e heavy loading combination and large size of the foundation model can make these differences not obvious. In this example, the external loads are in x and z directions and the length, width, and height of the foundation are three times of the structure. e 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.

Example 2.
Example 1 shows that our method has well adaptability. e local disengagement of contact surface is shown in example 2. e basic situation of this example is the same as example 1, but the external forces are different with example 1. e upper surface of the structure is affected by pressure and tension. e pressure acting on the left square (4.5 m * 6 m) and right square (4.5 m * 6 m) is 10 kPa. e tension acting on the middle square (1.0 m * 6m) is 8.9 kPa. Figure 11 shows the contact forces in three directions on the 30th and the 1200th day. e 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. e total time of the process is usually less than 10 seconds. e deformation of the center point in the contact surface is shown in Figure 12, where the center point is in the disengagement area. e 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. e disengagement area does not always exist.

Conclusion
is paper reports a methodology to analyze the time affection on a structural soft foundation system. e contact forces are treated as unknown variables and act on the structure and foundation separately. e deformation of structure contains the statically determinate structural deformation and rigid body displacement. It can be solved by the finite method. e creep deformation of soft foundation is solved by Lee's method. e final equations are formed by the deformation coordination equations and equilibrium equations. It is solved by the Gaussian elimination with the partial pivoting method. e augmented Lagrangian method is used to verify the accuracy of the method. e following important conclusions are summarized: (1) e 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) is 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. e redistribution phenomenon of contact forces is well shown with time changing.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.