An Explicit Finite Difference Method for Dynamic Interaction of Damped Saturated Soil Site-Pile Foundation-Superstructure System and Its Shaking Table Analysis

: The saturated soil site-pile foundation-superstructure system, with large degrees of freedom or strong nonlinear problems, often involves a large amount of calculation and low computational efﬁciency. In this paper, a fully explicit ﬁnite difference method is proposed for a saturated soil site-pile foundation-superstructure system. Since the proposed method has the advantages of decoupling in both time and space, it does not need to solve the equations simultaneously, which grants it high computational efﬁciency. At the same time, the method is implemented on the open-source software OpenSees and is used to compare and analyze the dynamic responses of the shaking table of a liqueﬁable soil site-pile foundation-superstructure system. After the calculation and analysis, the numerical solutions were found to be in good agreement with the experimental solutions, which veriﬁes the proposed method and illustrates that the proposed method can reasonably simulate the seismic responses of the whole system. In addition, the proposed calculation platform in OpenSees can be used for the analysis of the liquefaction process and the possible large deformation of soil after liquefaction, as well as for analyzing the failure mode of the complex and nonlinear saturated soil sites and structures under the effects of an earthquake.


Introduction
Soft soils or liquefiable soils are often regarded as saturated two-phase porous media composed of a solid granular soil skeleton and pore fluids. Practical engineering in a saturated soil site involves the saturated soil-structure interaction problem under dynamic loads. Especially under cyclic loading, saturated soil may gradually change from a solid state to liquid state, and the pore water pressure increases while the effective stress decreases [1]. Liquefaction can easily lead to the failure of the site soil, resulting in the instability or destruction of the foundation and the coupled superstructure. There are many cases of structural failures caused by seismic liquefaction at home and abroad. In the 1906 San Francisco earthquake, the site soil liquefied, and the piles of the Salinas railway bridge moved along the river, leading to the bridge's collapse [2]; in the 1964 Niigata earthquake, the liquefaction of the loose sand layer 10 m underground, and the subsequent collapses of the Showa Bridge as well as the beams that fell off of the five-span bridge into the river were caused by the excessive deformation of the steel pipe piles [3]; in the 1993 Chile earthquake, the relative displacement between the bridge abutment and the site soil caused by soil liquefaction resulted in damage to the piles of the Llacohen bridge [4]; in the 1995 Kobe earthquake in Japan, a lateral flow of soil caused by the site liquefactionseriously damaged the piles of the Nishinomiya Bridge, leading to the collapse of the piers and the fall of a beam [5]; and in the 2008 Wenchuan earthquake, a large area of method. The estimation results were in good agreement with the centrifuge test and field test results. Su et al. [40,41] carried out a numerical analysis on the saturated soil-pile-quay wall system and discussed the failure mechanism of the quay wall and evaluated three different types of quay wall seismic design schemes. Ghasemi-Fare et al. [42] simulated the lateral expansion of the inclined liquefaction site based on the u-p equation. Considering the effect of the pile group, Sun et al. [43] carried out a numerical simulation of the pile-soilpier interaction based on the OpenSees software. Yu [44] studied the dynamic responses of a pile-nuclear island structure and analyzed the effects of the input seismic wave, the layered soil, and the mass and stiffness of the superstructure on the seismic response of the piles. The results show that when the input pulse period decreases, the mass of the superstructure increases and the stiffness decreases, and thus the bending point of the pile moves to the top of the pile. The mutation of a pile's internal force easily occurs in the interface between soft and hard soil.
It is difficult for experiments to reflect the real site soil properties, the self-weight stress of the soil and the structure, and the interaction between the piles and soil, and these experiments are limited by the test site, model conditions, costs, and labor factors; consequently, it is not practical to carry out a great amount of test research. The finite element method is usually used to analyze the seismic responses of the saturated soil site-pile-superstructure system, which often encounters many degrees of freedom and strong nonlinear problems. The problems of a large amount of calculation, a highly timeconsuming process, and a low computational efficiency will occur in the computational process when the implicit algorithm is adopted. In this paper, the saturated soil site-pilessuperstructure system is regarded as a whole system, and a fully explicit finite difference method for the whole system is proposed. In the case of the site liquefaction triggered by earthquakes, the proposed method is used to calculate and analyze the seismic responses of the shaking table of a liquefiable site-pile foundation-superstructure system. The proposed method is verified by comparing the numerical results with the experimental results. By analyzing the responses of the system, the liquefaction process and the pore pressure development mechanism, the lateral displacement and seismic subsidence mechanism of the site soil, and the influence of site liquefaction on the mechanical behavior and failure mode of the structure are in turn analyzed.

Explicit Difference Method of the System
In the system of saturated soil-pile foundation-superstructure, saturated soil can be approximately regarded as a saturated two-phase media composed of a soil skeleton and pore fluid. Due to the different control equations and elements, the dynamic equations of the soil and structure need to be calculated by different spatial discretization and time domain algorithms, respectively. The nonlinear dynamic interaction problem of the saturated soil-pile foundation-superstructure system often involves a large number of degrees of freedom, and computational efficiency is one of the main bottlenecks during the calculation. An explicit difference method is proposed to solve this problem.

Nonlinear Finite Element Equation of Damped System
After the finite element discretization, the u-p-formed nonlinear dynamic equations with damping for the saturated soil subsystem can be expressed as: where subscript 1 represents the site soil system; u, . u, and ..
u are the displacement, velocity, and acceleration of the soil, respectively; p is the pore pressure of soil; M and C are the mass and damping matrices, respectively; R represents the nonlinear internal force; Q is the coupling matrix; S and J are the compressibility and permeability matrix of the pore fluid, respectively; and fu and fq are the external mechanical loadings, respectively. For the pile foundation-superstructure subsystem, it is necessary to use appropriate elements to discretize the finite element model according to the different structural characteristics. The superstructure and pile's foundation are discretized by beam elements, while the pile cap is discretized by two-dimensional four node elements. Since the node number of each element and the degrees of freedom of each node are different, the nonlinear dynamic equation with damping of the subsystem is expressed as follows: where subscript 2 represents the superstructure, pile foundation, and pile cap subsystem, and f ui is the external loading of structure.
Combining Equations (1)-(3), the dynamic equation of the saturated soil-pile foundationpile cap-superstructure system can be expressed as: where the matrices and unknown vectors are expressed as follows: where p = pdt represents a meaningless scalar, which is only used for aiding the calculation.
In general, the finite element Equation (4) of the whole system is often solved by an implicit time-domain algorithm, such as the Newmark method or the Wilson-θ method. For problems with a large number of degrees of freedom or strong nonlinear problems, the implicit method needs to solve the equation system simultaneously, which will occupy larger computational costs and greatly affect the computational efficiency.

Explicit Difference Formulation
Divide the total load time into time intervals with time steps ∆t. Any time can be expressed as t k = t k−1 + ∆t, (k = 1, 2, 3 . . . ). The acceleration .. u is expanded by the central difference method .. u(t k ) = [u(t k+1 ) − 2u(t k ) + u(t k−1 )]/∆t 2 , and the pore pressure . p k and velocity . u k are expanded by the forward difference method . p(t k ) = [p(t k+1 ) − p(t k )]/∆t and the backward difference method . u(t k ) = [u(t k ) − u(t k−1 )]/∆t, respectively. The selfstarting explicit difference formulation for solving the finite element, Equation (4), of the saturated soil-pile foundation-superstructural system can then be obtained: Equation (7) has the following equivalent global equation form: where the matrices and vectors are shown as follows: Buildings 2022, 12, 1186 The real displacement and pore pressure of each node in the whole system are obtained by solving Equation (11): In Equation (8), the mass matrix M and the fluid compression matrix S of each element are diagonalized in the following form, before the system matrix is integrated into the overall matrix according to the conventional method. The mass and compression deformation are concentrated on each node, and the coupling between the adjacent nodes is ignored. The diagonal matrix formed in this way is only related to the value of the diagonal, which avoids the problem of simultaneously solving equations when solving the inverse of the M and S matrices.
where M e l and M e are the lumped and consistent mass matrices in an element, respectively, and S e l and S e are the lumped and consistent fluid compressibility matrices in an element, respectively.
Through Equations (8)- (11), the displacement and pore pressure of the current step are calculated recursively by the displacement and pore pressure solutions of the previous two steps, respectively. Therefore, compared with the implicit method, the proposed explicit difference formulation has the diagonalization of the matrix and the explicit stepby-step recursion formulation, which has the characteristics of decoupling both in time and space. This method does not need to solve the equations simultaneously, which saves computational time and demonstrates its high computational efficiency.

OpenSees Implementation of the Method
As a program platform with open-source code, OpenSees is mainly used for seismic response analysis and simulations in structural and geotechnical engineering.
To efficiently calculate and analyze the dynamic interaction of the nonlinear saturated soil site-pile foundation-superstructure system, the explicit difference formulation of the whole system is implemented in the open-source finite element analysis software OpenSees. The calculation flow chart is shown in Figure 1: The element quadUP program with the diagonal mass matrix and fluid compression matrix is embedded in Element module, and the UPExplicitdifferenceinteraction program is embedded in the Integrator module, which mainly includes the following calculation steps: (1) Distinguish the node numbers of the soil, pile foundation (superstructure), and pile cap subsystem expressed by different element models; (2) Based on the explicit difference formulation of different subsystems, namely, Equation (8), the solutions of each node at t k+1 are deduced according to the results at time t k−1 and t k ; (3) The real solutions are updated according to Equation (11).

Numerical Verification
The simplified model of the saturated soil site-pile foundation-superstructure system is shown in Figure 2. The bottom of the model is set as an impermeable rigid horizontal boundary, and the surface of the site is a free permeable boundary. A simplified shear Buildings 2022, 12, 1186 6 of 24 boundary is used for the left and right lateral boundaries. The horizontal acceleration history record of the El Centro wave, shown in Figure 2, was selected as the input ground motion. and space. This method does not need to solve the equations simultaneously, which saves computational time and demonstrates its high computational efficiency.

OpenSees Implementation of the Method
As a program platform with open-source code, OpenSees is mainly used for seismic response analysis and simulations in structural and geotechnical engineering.
To efficiently calculate and analyze the dynamic interaction of the nonlinear saturated soil site-pile foundation-superstructure system, the explicit difference formulation of the whole system is implemented in the open-source finite element analysis software OpenSees. The calculation flow chart is shown in Figure 1: The element quadUP program with the diagonal mass matrix and fluid compression matrix is embedded in Element module, and the UPExplicitdifferenceinteraction program is embedded in the Integrator module, which mainly includes the following calculation steps: (1) Distinguish the node numbers of the soil, pile foundation (superstructure), and pile cap subsystem expressed by different element models; (2) Based on the explicit difference formulation of different subsystems, namely, Equation (8), the solutions of each node at 1 k t + are deduced according to the results at time 1 k t − and k t ; (3) The real solutions are updated according to Equation (11).

Numerical Verification
The simplified model of the saturated soil site-pile foundation-superstructure system is shown in Figure 2. The bottom of the model is set as an impermeable rigid horizontal boundary, and the surface of the site is a free permeable boundary. A simplified shear boundary is used for the left and right lateral boundaries. The horizontal acceleration history record of the El Centro wave, shown in Figure 2, was selected as the input ground motion. In Figure 2, the model site is composed of four layers of saturated soil with different permeability, which are a loose sand layer 1, dense sand layers 2-3, and a clay layer. Based on the OpenSees software, the quad elementUP is used for the saturated soil element, and the PDMY and PIMY constitutive models [45] are used for the materials of saturated sand and clay. The element and material parameters are shown in Table 1. The pile foundation and superstructure are simulated by dispBeamColumn element, and the same reinforced concrete material is used. The material parameters are shown in Table 2.  In Figure 2, the model site is composed of four layers of saturated soil with different permeability, which are a loose sand layer 1, dense sand layers 2-3, and a clay layer. Based on the OpenSees software, the quad elementUP is used for the saturated soil element, and the PDMY and PIMY constitutive models [45] are used for the materials of saturated sand and clay. The element and material parameters are shown in Table 1. The pile foundation and superstructure are simulated by dispBeamColumn element, and the same reinforced concrete material is used. The material parameters are shown in Table 2.  The seismic responses of the model are calculated by the Newmark method and the proposed method, respectively. As shown in Figure 2, soil nodes with different depths are selected. The horizontal displacement and pore pressure time-history curves of each node are shown in Figure 3. The horizontal displacement time-history curves of the pile foundation and superstructure nodes are shown in Figure 4. By comparing the calculation results of the Newmark method and the proposed method, the displacement and pore pressure results of each node obtained by the two methods are in good agreement, which verifies the accuracy of the proposed method.
proposed method, respectively. As shown in Figure 2, soil nodes with different depths are selected. The horizontal displacement and pore pressure time-history curves of each node are shown in Figure 3. The horizontal displacement time-history curves of the pile foundation and superstructure nodes are shown in Figure 4. By comparing the calculation results of the Newmark method and the proposed method, the displacement and pore pressure results of each node obtained by the two methods are in good agreement, which verifies the accuracy of the proposed method.   proposed method, respectively. As shown in Figure 2, soil nodes with different depths are selected. The horizontal displacement and pore pressure time-history curves of each node are shown in Figure 3. The horizontal displacement time-history curves of the pile foundation and superstructure nodes are shown in Figure 4. By comparing the calculation results of the Newmark method and the proposed method, the displacement and pore pressure results of each node obtained by the two methods are in good agreement, which verifies the accuracy of the proposed method.

Experiment Description
The shaking table test model, shown in Figure 5a,b, shows the site's layer distribution, which is composed of a dense sand layer, liquefiable sand layer, and a clay layer from bottom to top. Figure 6 shows the arrangement of the sensors used in the model test, mainly measuring the acceleration and pore pressure of soil, and the acceleration of the structure.

Experiment Description
The shaking table test model, shown in Figure 5a,b, shows the site's layer distribution, which is composed of a dense sand layer, liquefiable sand layer, and a clay layer from bottom to top. Figure 6 shows the arrangement of the sensors used in the model test, mainly measuring the acceleration and pore pressure of soil, and the acceleration of the structure.

Experiment Description
The shaking table test model, shown in Figure 5a,b, shows the site's layer di tion, which is composed of a dense sand layer, liquefiable sand layer, and a clay from bottom to top. Figure 6 shows the arrangement of the sensors used in the mod mainly measuring the acceleration and pore pressure of soil, and the acceleration structure.

Experiment Description
The shaking table test model, shown in Figure 5a,b, shows the site's layer dis tion, which is composed of a dense sand layer, liquefiable sand layer, and a clay from bottom to top. Figure 6 shows the arrangement of the sensors used in the mode mainly measuring the acceleration and pore pressure of soil, and the acceleration o structure.
(a) (b)  The acceleration time-history record of Wolong station from the Wenchuan earthquake was selected as the input ground motion, and the peak acceleration was 0.3 g. The acceleration time-history and its Fourier spectrum are shown in Figure 7. The acceleration time-history record of Wolong station from the Wenchuan earthquake was selected as the input ground motion, and the peak acceleration was 0.3 g. The acceleration time-history and its Fourier spectrum are shown in Figure 7.

Element and Material Definition
The quadrilateral plane strain element quadUP was used for loose sand, dense sand, and clay. Each node has three degrees of freedom corresponding to horizontal displacement, vertical displacement, and pore pressure, respectively. The pile cap was simulated by the quadrilateral plane strain quad element. Each node includes two degrees of free-

Element and Material Definition
The quadrilateral plane strain element quadUP was used for loose sand, dense sand, and clay. Each node has three degrees of freedom corresponding to horizontal displacement, vertical displacement, and pore pressure, respectively. The pile cap was simulated by the quadrilateral plane strain quad element. Each node includes two degrees of freedom of the horizontal and vertical displacements. The element parameters are listed in Table 3. The pile foundation and the superstructure were simulated by the dispBeamColumn element, and the section fiber and section fiberSec were used, respectively. Table 3. Element parameters of the foundation soil and the pile cap.
To truly describe the liquefaction process of soil, the PDMY02 material model [45] was used for the liquefiable sand and the PIMY material model [45] was used for clay. The material parameters are shown in Table 4. The same concrete material was used for the pile foundation and pile cap, and the H-shaped steel was used for the superstructure. The material parameters are shown in Table 5.

Boundaries, Damp and Mass of the Superstructure
The bottom of the model was set as an impermeable boundary, and the surface of the saturated soil site was a free permeable boundary. The simplified shear boundary condi-tions were used for the two lateral boundaries. The mass of each layer of the superstructure was 450 kg.
In the case of small strain, soil materials usually show nonlinearity and have an energydissipation capacity. To simulate the dynamic nonlinear characteristics of soil, Rayleigh damping was introduced to ensure the stability of the calculation, and the damping ratio was 3%.

Numerical Results and Analysis
Seismic liquefaction easily leads to the large lateral displacement and vertical settlement of the site soil. The deformation of the site soil is one of the main reasons for the stress change and even damage of the structure. The shaking table test of the liquefaction site-pile foundation-superstructure system was carried out, and the proposed method was used to calculate and analyze the seismic responses of the test model. By comparing the experimental and numerical solutions, the proposed method was further verified. Through the model test and numerical results, the pore pressure development mechanism of the site soil, the horizontal displacement and seismic deposition mechanism of the soil, and the forces of the structure were analyzed.
(1) Responses and development mechanism of the pore pressure The pore-pressure responses of site soil are an important result to obtain for understanding whether sand liquefaction has occurred and the stress state of soil. Pore pressure sensors (W1-W3 and W7-W11) were set at far-field and near-pile locations, as shown in Figure 6. W1 and W7, W2 and W9, and W3 and W11 were located in the same depth of the sand layer, respectively. The W1 and W7 nodes were near the surface, W2 and W9 were located in the middle of the sand layer, and W3 and W11 were the bottom nodes of the sand layer. The time-history curves of the excess pore pressure of the far-field nodes W1-W3 and the near-pile nodes W7-W11 are shown in Figure 8. With the increase in the soil depth, the excess pore pressures of the nodes along the depth direction increased gradually. In Figure 8, it can be seen that the excess pore pressures of the soil nodes are divided into three stages: (1) from 0-5 s, the excess pore pressure is maintained at a low level and grows slowly; (2) from 5-40 s, with the increase in the input ground motion amplitude, the excess pore pressure increases rapidly. The amplitude of the seismic acceleration reaches a maximum, such as 9.5 s and 19.9 s, corresponding to the presence of two peaks of pore pressure at 9.5 s and 19.9 s; and (3) from 40-60 s, the amplitude decreases gradually and is maintained at a low amplitude, and the pore pressure of the sand layer shows a slow dissipation process due to the low permeability of the overlying clay layer.
As shown in Figure 7, the test results of the excess pore pressures of W7 and W8 are significantly greater in magnitude than the numerical results, because during the test, the shallow soil is always in a liquefied or nearly liquefied state, which causes the pore pressures sensors at these locations to sink, leading to the test recorders becoming greater; alternatively, since the measuring points W7 and W8 are adjacent to the piles and the pile cap, the structural vibration causes the densification of the surrounding soil, leading to a greater actual vertical effective-overburden stress. Additionally, for measuring points W3 and W11, the numerical results are significantly higher than the test results, and this difference was persistent despite our efforts to adjust the permeability of the saturated soil in the FE simulation. A possible explanation could be that during the experimental design and instrument arrangement, the pore pressure sensors were fixed in their position by strings (the upper ends were fixed) that were cut off after filling the soil and were thus left in the soil. During the test, the pore water moved upward and was drained through the weakest drainage channels in the soil, such as the interfaces between the pile cap and the soil as well as the positions of the strings, thereby changing the drainage conditions of the soil to some extent and accordingly leading to a decrease in the excess pore pressure, as shown in the test results. Additionally, in the test we directly observed some macroscopic phenomena, such as a large amount of fluid seeping out at the interface between the pile cap and the clay, as well as water spraying and sand emission at the weak part of the soil where the strings of the sensors were placed; therefore, the difference between the numerical and test results has been preliminarily explained. The placement precision of the measuring points in the test may be another reason for the difference between the test and numerical results.
phenomena, such as a large amount of fluid seeping out at the interface between the pile cap and the clay, as well as water spraying and sand emission at the weak part of the soil where the strings of the sensors were placed; therefore, the difference between the numerical and test results has been preliminarily explained. The placement precision of the measuring points in the test may be another reason for the difference between the test and numerical results. Figure 8 compares the excess pore pressures at the far-field (W1-W3) and near-pile (W7-W11) measuring points. The pore pressures at the near-pile measuring points exhibited relatively significant small-amplitude oscillations under both the main shock and aftershock (stages III and IV), indicating that the shear-induced dilation tendency of the soil occurred. This phenomenon has been reported in a previous experimental study [46].  Figure 8 compares the excess pore pressures at the far-field (W1-W3) and near-pile (W7-W11) measuring points. The pore pressures at the near-pile measuring points exhibited relatively significant small-amplitude oscillations under both the main shock and aftershock (stages III and IV), indicating that the shear-induced dilation tendency of the soil occurred. This phenomenon has been reported in a previous experimental study [46].
The numerical solutions and experimental solutions of the excess pore pressure are compared in Figure 8. The variation trends of the two curves agree well, and the time of the peak value is also consistent. With the accumulation of the pore pressure, the excess pore pressure of the saturated soil site reaches the maximum value of 13.189 kPa at 19.9 s, as shown in Figure 9. In the pore pressure dissipation stage, the results of the two methods tend to be consistent, indicating that the numerical solution of the proposed method is in good agreement with the experimental solution.
In the test, there was a large amount of water on the surface of the site soil, accompanied by water spraying and sand emission, indicating that the sand layer liquefied under the earthquake conditions. Figure 10 shows the liquefaction phenomenon of the test, and it can be seen that the sand liquefaction near the pile foundation is more obvious. There are multiple water-emitting points, accompanied by sand gushing.
To reflect the liquefaction of the sand layer, Figure 11 shows the pore pressure ratio of nodes W1-W3. Comparing the pore pressure ratio time-history curves of each node in the figure, the pore pressure ratios of W1 and W7 are close to 1.0, and the pore pressure ratios of other nodes are less than 0.7, indicating that the sand close to the clay layer liquefied. With the decrease in the seismic acceleration amplitude, the pore water is blocked by the impermeable upper clay layer, and the pore pressure ratios of the W1 and W7 nodes decrease slowly. Other nodes do not liquefy or liquefy incompletely. Their pore pressure ratio gradually decreases, and the soil returns to the stable state with the decrease in the seismic acceleration amplitude.
The numerical solutions and experimental solutions of the excess pore pressure are compared in Figure 8. The variation trends of the two curves agree well, and the time of the peak value is also consistent. With the accumulation of the pore pressure, the excess pore pressure of the saturated soil site reaches the maximum value of 13.189 kPa at 19.9 s, as shown in Figure 9. In the pore pressure dissipation stage, the results of the two methods tend to be consistent, indicating that the numerical solution of the proposed method is in good agreement with the experimental solution. In the test, there was a large amount of water on the surface of the site soil, accompanied by water spraying and sand emission, indicating that the sand layer liquefied under the earthquake conditions. Figure 10 shows the liquefaction phenomenon of the test, and it can be seen that the sand liquefaction near the pile foundation is more obvious. There are multiple water-emitting points, accompanied by sand gushing. To reflect the liquefaction of the sand layer, Figure 11 shows the pore pressure ratio of nodes W1-W3. Comparing the pore pressure ratio time-history curves of each node in the figure, the pore pressure ratios of W1 and W7 are close to 1.0, and the pore pressure ratios of other nodes are less than 0.7, indicating that the sand close to the clay layer liquefied. With the decrease in the seismic acceleration amplitude, the pore water is blocked by the impermeable upper clay layer, and the pore pressure ratios of the W1 and W7 nodes decrease slowly. Other nodes do not liquefy or liquefy incompletely. Their pore pressure  In the test, there was a large amount of water on the surface of the site soil, accompanied by water spraying and sand emission, indicating that the sand layer liquefied under the earthquake conditions. Figure 10 shows the liquefaction phenomenon of the test, and it can be seen that the sand liquefaction near the pile foundation is more obvious. There are multiple water-emitting points, accompanied by sand gushing. To reflect the liquefaction of the sand layer, Figure 11 shows the pore pressure ratio of nodes W1-W3. Comparing the pore pressure ratio time-history curves of each node in the figure, the pore pressure ratios of W1 and W7 are close to 1.0, and the pore pressure ratios of other nodes are less than 0.7, indicating that the sand close to the clay layer liquefied. With the decrease in the seismic acceleration amplitude, the pore water is blocked by the impermeable upper clay layer, and the pore pressure ratios of the W1 and W7 nodes decrease slowly. Other nodes do not liquefy or liquefy incompletely. Their pore pressure The pore pressure ratio of the test points was selected along the depth and horizontal direction of the soil to determine the main liquefaction area of the sand layer. The maximum pore pressure ratios of the sand nodes are compared along the soil depths in Figure  12. The pore pressure ratio of the sand nodes close to the upper clay layer is close to or equal to 1.0, and the other pore pressure ratio of the sand nodes is less than 1.0. The figure shows that the pore water accumulates at the interface between the clay layer and the loose sand layer. Until the end of the seismic input, the pore pressure ratio at this interface is always equal to or approximately equal to 1.0, which indicates that the sand near the interface is always in a liquefied state.  The pore pressure ratio of the test points was selected along the depth and horizontal direction of the soil to determine the main liquefaction area of the sand layer. The maximum pore pressure ratios of the sand nodes are compared along the soil depths in Figure 12. The pore pressure ratio of the sand nodes close to the upper clay layer is close to or equal to 1.0, and the other pore pressure ratio of the sand nodes is less than 1.0. The figure shows that the pore water accumulates at the interface between the clay layer and the loose sand layer. Until the end of the seismic input, the pore pressure ratio at this interface is always equal to or approximately equal to 1.0, which indicates that the sand near the interface is always in a liquefied state. mum pore pressure ratios of the sand nodes are compared along the soil depths in Figure  12. The pore pressure ratio of the sand nodes close to the upper clay layer is close to or equal to 1.0, and the other pore pressure ratio of the sand nodes is less than 1.0. The figure shows that the pore water accumulates at the interface between the clay layer and the loose sand layer. Until the end of the seismic input, the pore pressure ratio at this interface is always equal to or approximately equal to 1.0, which indicates that the sand near the interface is always in a liquefied state. At the same time, Figure 13 compares the pore pressure ratios of the near-pile nodes and far-field nodes at the same depth. The W1 and W7 nodes were located in the upper part of the sand layer, and the W4 and W12 nodes were located in the middle of the sand layer. The pore pressure ratios of the near-pile nodes W7 and W12 are significantly greater than those of the far-field nodes W1 and W4. The analysis shows that the saturated sand layer near the surface clay liquefy. Affected by the piles, the near-pile sand nodes were more prone to liquefaction than the far-field sand nodes, which is consistent with the phenomenon of the test shown in Figure 10. At the same time, the pore pressure ratio oscillations of the far-field soil nodes is small and the curve is relatively smooth, while the pore pressure oscillations of the near-pile soil nodes are obvious, indicating that the pile vibration affects the pore pressure development of the soil around the pile. At the same time, Figure 13 compares the pore pressure ratios of the near-pile nodes and far-field nodes at the same depth. The W1 and W7 nodes were located in the upper part of the sand layer, and the W4 and W12 nodes were located in the middle of the sand layer. The pore pressure ratios of the near-pile nodes W7 and W12 are significantly greater than those of the far-field nodes W1 and W4. The analysis shows that the saturated sand layer near the surface clay liquefy. Affected by the piles, the near-pile sand nodes were more prone to liquefaction than the far-field sand nodes, which is consistent with the phenomenon of the test shown in Figure 10. At the same time, the pore pressure ratio oscillations of the far-field soil nodes is small and the curve is relatively smooth, while the pore pressure oscillations of the near-pile soil nodes are obvious, indicating that the pile vibration affects the pore pressure development of the soil around the pile. (2) Lateral deformation and seismic subsidence mechanism of soil Figure 14 shows the numerical time-history curves of the relative horizontal displacements of the soil at different depths. It can be seen from the figure that the change trend of the relative horizontal displacement of each measuring point is consistent, and the maximum horizontal displacement of each point occurs at about 9.5 s of the main shock. With the decrease in the soil depth, the lateral displacement of each point gradually increased, and the relative displacement of the ground surface reached 8.6 mm.  (2) Lateral deformation and seismic subsidence mechanism of soil Figure 14 shows the numerical time-history curves of the relative horizontal displacements of the soil at different depths. It can be seen from the figure that the change trend of the relative horizontal displacement of each measuring point is consistent, and the maximum horizontal displacement of each point occurs at about 9.5 s of the main shock. With the decrease in the soil depth, the lateral displacement of each point gradually increased, and the relative displacement of the ground surface reached 8.6 mm. Figure 15a shows the contour map of the site horizontal displacement at 9.5 s, and the deformation effect is magnified 20 times. The soil displacement shows an increasing trend from bottom to top, and the shear capacity of the liquefied loose sand is poor, resulting in a relatively large lateral deformation. The largest horizontal displacement of the site soil occurs at 9.5 s when the input seismic acceleration peak reaches 0.3 g. The blue area in the figure shows that the maximum horizontal displacement of the interface between the sand layer and the clay layer reaches 8.88. mm, and the time-history curve is shown in Figure 15b. Moreover, the horizontal displacement of the clay layer presents the tendency of the overall lateral displacement in the figure, indicating that after the sand's liquefaction, the soil particles settle and rearrange, and the pore fluid is discharged upward. However, due to the poor permeability of the overlying clay, the drainage is blocked, and the pore water is continuously gathered at the interface between the clay layer and the sand layer, which is equivalent to a thin film with an extremely low shear strength formed at the interface. This leads to the large horizontal lateral displacement of the clay layer. (2) Lateral deformation and seismic subsidence mechanism of soil Figure 14 shows the numerical time-history curves of the relative horizontal displacements of the soil at different depths. It can be seen from the figure that the change trend of the relative horizontal displacement of each measuring point is consistent, and the maximum horizontal displacement of each point occurs at about 9.5 s of the main shock. With the decrease in the soil depth, the lateral displacement of each point gradually increased, and the relative displacement of the ground surface reached 8.6 mm.  Figure 15a shows the contour map of the site horizontal displacement at 9.5 s, and the deformation effect is magnified 20 times. The soil displacement shows an increasing trend from bottom to top, and the shear capacity of the liquefied loose sand is poor, resulting in a relatively large lateral deformation. The largest horizontal displacement of the site soil occurs at 9.5 s when the input seismic acceleration peak reaches 0.3 g. The blue area in the figure shows that the maximum horizontal displacement of the interface between the sand layer and the clay layer reaches 8.88. mm, and the time-history curve is shown in Figure 15b. Moreover, the horizontal displacement of the clay layer presents the tendency of the overall lateral displacement in the figure, indicating that after the sand's liquefaction, the soil particles settle and rearrange, and the pore fluid is discharged As shown in Figure 16a, the site soil produces a small compaction under gravity. Figure 16b shows the final vertical displacement of the site soil under the earthquake test conditions. It can be seen that the settlement of the upper loose sand layer is larger, and the vertical displacements of the pile cap and the soil below the pile cap are less than zero, indicating that the soil below the pile cap has a small subsidence after the earthquake; the maximum seismic subsidence value is 10.09 mm. Due to the support of the pile foundation, the vertical displacement of the pile cap is 2.4 mm, which is significantly less than the seismic subsidence of the soil below the pile cap. Moreover, the vertical displacement of the loose sand adjacent to the sand layer and clay layer interface is smaller than that of the lower sand (blue rectangle block), which further indicates that pore water gathered at the interface, and that a water film formed at the interface, which reduced the settlement of the upper soil.
To further understand the seismic subsidence and uplift of the soil surface, the vertical displacements of the J1-J3 soil nodes and the nodes of the pile cap surface are shown in Figure 17. The maximum vertical displacement of the pile cap reaches 1.8 mm, and node J1 near the pile cap has a small seismic subsidence of 1.1mm. In contrast, the vertical displacements of the J2 and J3 nodes are greater than zero, which illustrates that the surrounding soil was uplifted due to the dilatancy effect, and the maximum of uplift soil reached 1.9 mm.
(3) Acceleration responses and structural force mechanism To monitor the influence of the liquefaction on wave propagation, the acceleration sensors were arranged relatively far from the piles (SAA2 in Figure 6). Figure 18 shows the numerical and experimental acceleration time-history curves of SAA2- (1)(2)(3)(4)(5)(6). The trends of the two curves are consistent, and the peak acceleration times of the numerical responses are consistent with the experimental responses. Compared with the test results, the numerical solutions overestimate the acceleration amplitudes of SAA2. In the test phenomenon, it was found that during the test process, there was pore-water spraying and sand emission at the weak part of the soil by the strings of the sensors. The local soil particles were consolidated to a certain extent compared with the undrained situation. However, the drainage conditions are not defined in the numerical model. Therefore, the calculation results show that the soil acceleration amplitude of SAA2-6 is significantly smaller than the experimental solution, which also shows that the drainage conditions may affect the acceleration amplitude of the soil softening (after liquefaction). upward. However, due to the poor permeability of the overlying clay, the drainage is blocked, and the pore water is continuously gathered at the interface between the clay layer and the sand layer, which is equivalent to a thin film with an extremely low shear strength formed at the interface. This leads to the large horizontal lateral displacement of the clay layer. As shown in Figure 16a, the site soil produces a small compaction under gravity. Figure 16b shows the final vertical displacement of the site soil under the earthquake test conditions. It can be seen that the settlement of the upper loose sand layer is larger, and the vertical displacements of the pile cap and the soil below the pile cap are less than zero, indicating that the soil below the pile cap has a small subsidence after the earthquake; the maximum seismic subsidence value is 10.09 mm. Due to the support of the pile foundation, the vertical displacement of the pile cap is 2.4 mm, which is significantly less than the seismic subsidence of the soil below the pile cap. Moreover, the vertical displacement of the loose sand adjacent to the sand layer and clay layer interface is smaller than that of the lower sand (blue rectangle block), which further indicates that pore water gathered at the interface, and that a water film formed at the interface, which reduced the settlement of the upper soil.  upward. However, due to the poor permeability of the overlying clay, the drainage is blocked, and the pore water is continuously gathered at the interface between the clay layer and the sand layer, which is equivalent to a thin film with an extremely low shear strength formed at the interface. This leads to the large horizontal lateral displacement of the clay layer. As shown in Figure 16a, the site soil produces a small compaction under gravity. Figure 16b shows the final vertical displacement of the site soil under the earthquake test conditions. It can be seen that the settlement of the upper loose sand layer is larger, and the vertical displacements of the pile cap and the soil below the pile cap are less than zero, indicating that the soil below the pile cap has a small subsidence after the earthquake; the maximum seismic subsidence value is 10.09 mm. Due to the support of the pile foundation, the vertical displacement of the pile cap is 2.4 mm, which is significantly less than the seismic subsidence of the soil below the pile cap. Moreover, the vertical displacement of the loose sand adjacent to the sand layer and clay layer interface is smaller than that of the lower sand (blue rectangle block), which further indicates that pore water gathered at the interface, and that a water film formed at the interface, which reduced the settlement of the upper soil. To further understand the seismic subsidence and uplift of the soil surface, the vertical displacements of the J1-J3 soil nodes and the nodes of the pile cap surface are shown in Figure 17. The maximum vertical displacement of the pile cap reaches 1.8 mm, and node J1 near the pile cap has a small seismic subsidence of 1.1mm. In contrast, the vertical displacements of the J2 and J3 nodes are greater than zero, which illustrates that the To reflect the amplification or damping effects of site soil on the input earthquake, the comparisons between the experimental and the numerical acceleration-amplification factor of SAA2-(1-6) (the ratio of the peak acceleration of each node to the input peak acceleration of the bottom node) along the soil depths are shown in Figure 19. The earthquake was inputted from the bottom of the model. With the decrease in the soil depth, the acceleration-amplification factors of each measuring point show a decreasing trend. The peak accelerations of the measuring points in the dense sand layer are the largest, and the amplification factors of acceleration are close to 1.0. The acceleration-amplification factor of the loose sand layer was attenuated to 50% of the input seismic acceleration and reached to 0.65 in the clay layer. Since the loose sand layer was in a liquefied state due to the development of excess pore pressure, the seismic energy clearly dissipated in the liquefied soil, indicating that soil softening (after liquefaction) decreases the acceleration amplitude, which is consistent with the conclusion obtained in the references [47,48].

(3) Acceleration responses and structural force mechanism
To monitor the influence of the liquefaction on wave propagation, the acceleration sensors were arranged relatively far from the piles (SAA2 in Figure 6). Figure 18 shows the numerical and experimental acceleration time-history curves of SAA2-(1-6). The trends of the two curves are consistent, and the peak acceleration times of the numerical responses are consistent with the experimental responses. Compared with the test results, the numerical solutions overestimate the acceleration amplitudes of SAA2. In the test phenomenon, it was found that during the test process, there was pore-water spraying and sand emission at the weak part of the soil by the strings of the sensors. The local soil particles were consolidated to a certain extent compared with the undrained situation. However, the drainage conditions are not defined in the numerical model. Therefore, the calculation results show that the soil acceleration amplitude of SAA2-6 is significantly smaller than the experimental solution, which also shows that the drainage conditions may affect the acceleration amplitude of the soil softening (after liquefaction).

(3) Acceleration responses and structural force mechanism
To monitor the influence of the liquefaction on wave propagation, the acceleration sensors were arranged relatively far from the piles (SAA2 in Figure 6). Figure 18 shows the numerical and experimental acceleration time-history curves of SAA2-(1-6). The trends of the two curves are consistent, and the peak acceleration times of the numerical responses are consistent with the experimental responses. Compared with the test results, the numerical solutions overestimate the acceleration amplitudes of SAA2. In the test phenomenon, it was found that during the test process, there was pore-water spraying and sand emission at the weak part of the soil by the strings of the sensors. The local soil particles were consolidated to a certain extent compared with the undrained situation. However, the drainage conditions are not defined in the numerical model. Therefore, the calculation results show that the soil acceleration amplitude of SAA2-6 is significantly smaller than the experimental solution, which also shows that the drainage conditions may affect the acceleration amplitude of the soil softening (after liquefaction). To reflect the amplification or damping effects of site soil on the input earthquake, the comparisons between the experimental and the numerical acceleration-amplification factor of SAA2-(1-6) (the ratio of the peak acceleration of each node to the input peak acceleration of the bottom node) along the soil depths are shown in Figure 19. The earthquake was inputted from the bottom of the model. With the decrease in the soil depth, the acceleration-amplification factors of each measuring point show a decreasing trend. The peak accelerations of the measuring points in the dense sand layer are the largest, and the amplification factors of acceleration are close to 1.0. The acceleration-amplification factor of the loose sand layer was attenuated to 50% of the input seismic acceleration and reached to 0.65 in the clay layer. Since the loose sand layer was in a liquefied state due to the development of excess pore pressure, the seismic energy clearly dissipated in the liquefied soil, indicating that soil softening (after liquefaction) decreases the acceleration amplitude, which is consistent with the conclusion obtained in the references [47,48]. The comparisons between the experimental and numerical solutions of the pile acceleration are shown in Figure 20. The trends of the two curves agree well, and the peak acceleration times of the numerical responses are consistent with the experimental responses. Apart from SAA1-1, the numerical solution overestimates the acceleration response of the pile foundation compared with the experimental solution. Similar results were produced after several comparative computations, which might have been caused by the inaccurate placement of the experimental acceleration sensor or a lack of fully fixed connections between the piles and the pile cap in the experiment, and this is explained in the following analysis.
peak accelerations of the measuring points in the dense sand layer are the largest, and the amplification factors of acceleration are close to 1.0. The acceleration-amplification factor of the loose sand layer was attenuated to 50% of the input seismic acceleration and reached to 0.65 in the clay layer. Since the loose sand layer was in a liquefied state due to the development of excess pore pressure, the seismic energy clearly dissipated in the liquefied soil, indicating that soil softening (after liquefaction) decreases the acceleration amplitude, which is consistent with the conclusion obtained in the references [47,48]. The comparisons between the experimental and numerical solutions of the pile acceleration are shown in Figure 20. The trends of the two curves agree well, and the peak acceleration times of the numerical responses are consistent with the experimental responses. Apart from SAA1-1, the numerical solution overestimates the acceleration response of the pile foundation compared with the experimental solution. Similar results were produced after several comparative computations, which might have been caused by the inaccurate placement of the experimental acceleration sensor or a lack of fully fixed connections between the piles and the pile cap in the experiment, and this is explained in the following analysis. By comparing the experimental and numerical solutions, it was found that the acceleration-amplification factors of the pile decrease first and then increase with the decrease in the soil depth. At the interface of the loose sand layer and the clay layer, the acceleration-amplification factor of the pile head is larger than the other nodes. The pile acceleration-amplification factors in the loose sand layer attenuate to 60-70% of the input seismic acceleration. After the liquefaction of the loose sand layer, the seismic energy is clearly dissipated in the liquefied soil. Due to the soil-structure dynamic interaction, the seismic energy of the soil is fed back to the structure, and the structural acceleration responses decreases. Comparing Figures 18-21, it is clear that the waveforms and amplitudes of the acceleration time-history curves of the pile and the soil at the same depth are different, indicating that the vibrations of the pile and the soil are not consistent in the test and that the soil-structure dynamic interaction affects the seismic responses of the computational model. At the same time, the acceleration-amplification factors of the pile are shown in Figure 21. By comparing the experimental and numerical solutions, it was found that the acceleration-amplification factors of the pile decrease first and then increase with the decrease in the soil depth. At the interface of the loose sand layer and the clay layer, the acceleration-amplification factor of the pile head is larger than the other nodes. The pile acceleration-amplification factors in the loose sand layer attenuate to 60-70% of the input seismic acceleration. After the liquefaction of the loose sand layer, the seismic energy is clearly dissipated in the liquefied soil. Due to the soil-structure dynamic interaction, the seismic energy of the soil is fed back to the structure, and the structural acceleration responses decreases. Comparing Figures 18-21, it is clear that the waveforms and amplitudes of the acceleration time-history curves of the pile and the soil at the same depth are different, indicating that the vibrations of the pile and the soil are not consistent in the test and that the soil-structure dynamic interaction affects the seismic responses of the computational model. According to the existing tests and the seismic damage observation, the connections of the piles and the pile cap are usually the locations with the most prominent bending moments and are the most vulnerable parts. During the preparation and performance of this experiment, the embedment method with reserved holes was adopted to assemble the piles and the pile cap. The piles were not fully fixed in the pile cap in this experiment. Additionally, since the superstructure was designed with the natural vibration period as the controlling factor, and the design of its mass and inertia force was simplified, the feedback of the bending moment of the superstructure to the piles was relatively small, resulting in a slight bending moment at the connection between the pile heads and the pile cap [49]. Figure 22 shows the maximum bending moment envelope of the pile along the soil depth. The bending moment of the pile tip is small. With the decrease in the soil depth, the bending moment of the pile node close to the interface between the loose sand layer and the dense sand layer reaches the maximum; then, the bending moment in the middle of the loose sand layer decreases gradually from the bottom to the top of the soil. The experimental bending moment envelope exhibits a trend similar to the numerical results. Compared with the numerical and experimental bending moments, the numerical solutions are larger than the experimental solutions. The pile bending moment reaches the maximum in the middle to lower part of the loose sand layer (near the interface between the loose and dense sand layers) due to the sand liquefaction mainly occurring in the middle to upper part of the loose sand layer. At the same time, it can be seen from the figure that the shear force at the interface between the loose sand layer and the clay layer suddenly increases. Due to the liquefaction of the shallow sand, the liquefied soil generates a discontinuous horizontal displacement, resulting in a sudden increase in the shear value of the piles at the interface. According to the existing tests and the seismic damage observation, the connections of the piles and the pile cap are usually the locations with the most prominent bending moments and are the most vulnerable parts. During the preparation and performance of this experiment, the embedment method with reserved holes was adopted to assemble the piles and the pile cap. The piles were not fully fixed in the pile cap in this experiment. Additionally, since the superstructure was designed with the natural vibration period as the controlling factor, and the design of its mass and inertia force was simplified, the feedback of the bending moment of the superstructure to the piles was relatively small, resulting in a slight bending moment at the connection between the pile heads and the pile cap [49]. Figure 22 shows the maximum bending moment envelope of the pile along the soil depth. The bending moment of the pile tip is small. With the decrease in the soil depth, the bending moment of the pile node close to the interface between the loose sand layer and the dense sand layer reaches the maximum; then, the bending moment in the middle of the loose sand layer decreases gradually from the bottom to the top of the soil. The experimental bending moment envelope exhibits a trend similar to the numerical results. Compared with the numerical and experimental bending moments, the numerical solutions are larger than the experimental solutions. The pile bending moment reaches the maximum in the middle to lower part of the loose sand layer (near the interface between the loose and dense sand layers) due to the sand liquefaction mainly occurring in the middle to upper part of the loose sand layer. At the same time, it can be seen from the figure that the shear force at the interface between the loose sand layer and the clay layer suddenly increases. Due to the liquefaction of the shallow sand, the liquefied soil generates a discontinuous horizontal displacement, resulting in a sudden increase in the shear value of the piles at the interface.
The superstructure was approximately replaced by the lumped mass, and three acceleration sensors were set at the corresponding mass nodes. The numerical accelerations of the superstructure were compared with the experimental solutions, as shown in Figure 23. The accelerations of the superstructure do not show an obvious amplification effect, and the trend of the numerical solutions is consistent with the experimental solutions. The numerical acceleration amplitudes are slightly larger than those of the experimental solutions. Although the numerical solutions do not adequately estimate the accelerations of the experiment in terms of amplitude, it was found by comparing the acceleration waveforms of the pile foundation that the waveforms of the superstructure acceleration are relatively close to the experimental results, and the difference between the two results may be affected by the parameters and the model of the superstructure. The superstructure was approximately replaced by the lumped mass, and three acceleration sensors were set at the corresponding mass nodes. The numerical accelerations of the superstructure were compared with the experimental solutions, as shown in Figure  23. The accelerations of the superstructure do not show an obvious amplification effect, and the trend of the numerical solutions is consistent with the experimental solutions. The numerical acceleration amplitudes are slightly larger than those of the experimental solutions. Although the numerical solutions do not adequately estimate the accelerations of the experiment in terms of amplitude, it was found by comparing the acceleration waveforms of the pile foundation that the waveforms of the superstructure acceleration are relatively close to the experimental results, and the difference between the two results may be affected by the parameters and the model of the superstructure.    The superstructure was approximately replaced by the lumped mass, and three acceleration sensors were set at the corresponding mass nodes. The numerical accelerations of the superstructure were compared with the experimental solutions, as shown in Figure  23. The accelerations of the superstructure do not show an obvious amplification effect, and the trend of the numerical solutions is consistent with the experimental solutions. The numerical acceleration amplitudes are slightly larger than those of the experimental solutions. Although the numerical solutions do not adequately estimate the accelerations of the experiment in terms of amplitude, it was found by comparing the acceleration waveforms of the pile foundation that the waveforms of the superstructure acceleration are relatively close to the experimental results, and the difference between the two results may be affected by the parameters and the model of the superstructure.    There is significant nonlinearity in each element of the soil. The shallow sand exhibits a larger shear strain, but the shear stress of the deep sand is twice that of the shallow sand. It shows that the complete liquefaction of the shallow sand has a great influence on the deformation and mechanical properties of the soil. Note that during 0-5 s, the shear stress-shear strain curve is nearly a straight line, indicating that the soil has not shown an obvious nonlinear response and is in an elastic stage. During 8-15 s, both the shear stress and shear strain increase rapidly; specifically, the shear strain is approximately six times that during 0-5 s, and the shear stress multiplies. During 15-25 s, under the aftershock, the shear stress and shear strain decreases compared with that during 8-15 s. During 55-60 s, the seismic loading weakens, and the soil gradually approaches a stable state, restoring the approximate linear shear stress-shear strain relation of the soil. shallow sand exhibits a larger shear strain, but the shear stress of the deep sand is twice that of the shallow sand. It shows that the complete liquefaction of the shallow sand has a great influence on the deformation and mechanical properties of the soil.  . Note that during 0-5 s, the shear stress-shear strain curve is nearly a straight line, indicating that the soil has not shown an obvious nonlinear response and is in an elastic stage. During 8-15 s, both the shear stress and shear strain increase rapidly; specifically, the shear strain is approximately six times that during 0-5 s, and the shear stress multiplies. During 15-25 s, under the aftershock, the shear stress and shear strain decreases compared with that during 8-15 s. During 55-60 s, the seismic loading weakens, and the soil gradually approaches a stable state, restoring the approximate linear shear stressshear strain relation of the soil.

Conclusions
In this paper, a fully explicit finite difference method was proposed to analyze the dynamic responses of the saturated soil site-pile foundation-superstructure system. By numerical and experimental comparations, the proposed method was verified, and the following conclusions were obtained:

Conclusions
In this paper, a fully explicit finite difference method was proposed to analyze the dynamic responses of the saturated soil site-pile foundation-superstructure system. By numerical and experimental comparations, the proposed method was verified, and the following conclusions were obtained: (1) Compared with the implicit method, the proposed explicit difference formulation has the diagonalization of the matrix and the explicit step-by-step recursion formulation, which have the characteristics of decoupling both in time and space; (2) This proposed method does not need to solve the equations simultaneously, which saves computational time and has high computational efficiency; (3) Liquefaction significantly affects the acceleration amplitude of the soil, indicating that the liquefied soil dissipates the seismic energy; (4) For the liquefiable soil layer with an overlying clay layer with an extremely low permeability, a water film layer with a very low shear strength is easily formed, which directly affects the horizontal and vertical discontinuous displacements at the interface, indicating that this material discontinuity results in the discontinuous displacement of the soils on the two sides.
The proposed method is embedded in open-source software OpenSees and can be used for the analysis of the liquefaction process and the possible large deformation of soil after liquefaction, as well as for analyzing the failure mode of a complex and nonlinear saturated soil site and structure under earthquake.
In this paper, the computational model size was relatively small, and the advantage of the efficient computation of the explicit method is not obvious due to the small number of degrees of freedom. However, the purpose of this paper is to propose a method and illustrate its capability for analyzing the dynamic interaction problem of the soil-structure system. In addition, compared with the unconditionally stable implicit method, the explicit difference method proposed in this paper is conditionally stable, which means that the calculation time step must be within a certain range to ensure the calculations' convergence. However, with the feature of a small-time step, the explicit method is more suitable for calculating problems under high-frequency loads compared to the implicit method. Regarding the computational efficiency, accuracy, and stability conditions of the proposed method, the authors will conduct a further analysis in subsequent application studies.