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.
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.
(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 aswhere is the deformation obtained by all forces in the statically determinate structural system; is the deformation obtained by the external forces; . are the deformation obtained by each unit force separately; and 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 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 coordination equations and 6 equilibrium equations. It can be simplified as follows:where is the rigid body displacement of structural contact surface; is the contact surface deformation of soft foundation; and 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 Figure 2, 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:where is the radial deformation; 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:where , , 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:where is obtained by the numerical integration method.
(a)
(b)
Kelvin viscoelastic solution is derived by inverse Laplace transformation.where K is the bulk modulus; G1 is the shear stiffness; and is the coefficient of viscosity. Tangential load causes the deformation; the elastic deformation is shown in the following equation [30]:
Through the Laplace transformation and inverse Laplace transformation, the viscoelastic deformation can also be solved and shown in the following equation:
is obtained by the numerical integration method, where , , 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: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 aswhere are the deformation of soft foundation in three directions. 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 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. 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:where are the shape function of the four nodes which are at the bottom of the hexahedron element and 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 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.
(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 iswhere are the normal strain in x, y, z directions and are the shearing strain in x, y, z directions. Solving the above formulas, the rigid body displacement can be obtained aswhere u, , are the rigid body displacement in x, y, z directions; are the rigid body rotation angles in x, y, z directions; and u0, , 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:
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:where are the sum of contact forces in x, y, z directions; are the coordinates of the contact forces’ functional point; are the coordinates of the external forces’ functional point; and are the external forces in x, y, z directions. When c = 0, (15) can be simplified as linear matrix equation:
7. Mixed Finite Element Methodology
Through 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. 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 the x, y, and z directions. They are used to connect the structure with foundation. The increment of the normal deformation and the increment of the rigid body rotation angles 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 aswhere is the deformation of i position under unit distribution force acting on k position in z direction in the time step l; is the deformation of i position under external forces at z direction in the time step l; is the coordinate of y-axis at i position in the time step l; is the coordinate of x-axis at i position in the time step l; is the element connecting rod force at z direction in the time step l; is the composition of forces at z direction in the time step l; 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 aswhere is the deformation of i position under unit distribution force acting on k position at m-axis direction in the time step l; are the deformation of foundation and structure under unit distribution force; 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 Figure 6. (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.
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: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.
Figure 7 shows the deformation of soft foundation in three directions under unit distribution force at 30th day and 1200th day. The deformation results of and are center symmetrical about the point of force. The deformation results of , and 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.
(a)
(b)
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. 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 in x 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 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.(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 in x direction.
(a)
(b)
(a)
(b)
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 in x 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.
Figure 11 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.
(a)
(b)
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.
Appendix
The analytical solution of viscoelastic deformation based on the Maxwell model can be calculated by MATLAB. The codes are as follows:
Data Availability
The data used to support the findings of this study are included in the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 51579089).