#### Abstract

In order to solve the time step setting problem in seismic liquefaction numerical analysis, the time step optimization method was introduced into seismic liquefaction numerical analysis platform, and the seismic response of earth dam was analyzed. The results show that the time step optimization method can effectively solve the time step setting problem. In addition, in order to evaluate the computational efficiency improvement of the optimized seismic response, a new method was proposed, and an example with analytical solution was selected to verify it. And the proposed method was used to evaluate the calculation efficiency of the optimized seismic response. The results show that the computing efficiency is greatly improved after the seismic liquefaction numerical analysis program is optimized by the time step optimization method.

#### 1. Introduction

In the process of dynamic time history analysis, the time step is too small in order to improve the calculation accuracy, resulting in high calculation cost. The time step is too large in order to save the calculation time, resulting in the divergence of the calculation results and even sudden interruption of the calculation. In order to improve the computational efficiency, it is difficult to find a suitable time step to save the computational cost and improve the computational accuracy. How to effectively improve dynamic numerical computational efficiency by adjusting the time step has always been a problem for researchers.

However, for time-discrete methods with a fixed time step, a completely opposite relationship exists between the time step and the computational accuracy. Hence, when applying time-discrete methods with a fixed time step, the problem of how to improve the computational efficiency is less likely to be effectively solved.

Time adaptive method can effectively solve the contradiction between time step and computational accuracy. It can save calculation time as much as possible on the basis of ensuring the calculation accuracy, so as to effectively improve the calculation efficiency. Many scholars have made efforts in the field of time adaptation. Zienkiewicz and Xie [1] proposed a simple time adaptive method that, at the time, was used for a dynamic analysis based on truncation error estimation. Later, Zeng et al. [2] and Wiberg and Li [3] modified the method individually. Sloan and Abbo [4] applied a time adaptive algorithm to a consolidation analysis in the field of geotechnical engineering. In addition, Kavetski et al. [5] adopted a heuristic algorithm and an error estimation method to carry out a time adaptive analysis. Tang et al. [6] employed an automatic time step to conduct a groundwater flow analysis by controlling global errors. Kuo Shyhrong et al. [7] developed a robust time-integration algorithm for solving nonlinear dynamic problems with large rotations and displacements. Wang Kaijia et al. [8] applied a time adaptive method to the numerical simulation of an offshore floating wind power generation platform. An adaptive time step method was also developed by Mao Weinan et al. [9] and applied to a water-heat coupling model of soil freezing. Ahmed and John [10] applied an adaptive time step control to convection–diffusion–reaction equations to achieve higher-order variational time discretization. In 2015, Marc Laforest et al. [11] created a preface to the special issue on error estimation and adaptivity for nonlinear and time-dependent problems. As mentioned in this preface, the current time adaptive methods are mostly posterior error evaluation methods [1–10, 12–14]. However, due to the defects of the posteriori error estimation time adaptive method, there are fewer and fewer research studies in recent years.

Different from the posteriori error estimation time adaptive method, Li Bin et al. [15, 16] proposed a time step optimization method, which was a useful attempt at an Apriori method. The time step optimization method not only improves the computational efficiency but also enriches the time adaptive theory. However, as analytical solutions are usually difficult to be obtained, it is impossible to evaluate the improvement degree of computational efficiency when the time step optimization method is used to improve computational efficiency. In order to solve this problem, a new method for evaluating the improvement degree was proposed, and a Bernoulli–Euler beam example with an analytical solution was used to verify the method. Then on this basis, the time step optimization method was introduced into the seismic liquefaction numerical analysis platform, and the seismic response of earth dam was analyzed. Finally, the seismic response computational efficiency after time step optimization was evaluated with the proposed method. The results show that the time step optimization method can effectively solve the time step setting problem in the seismic liquefaction numerical analysis and improve the calculation efficiency effectively while realizing the time step optimization.

#### 2. Determination of Optimization Time Step

For dynamic structural numerical analysis methods, a differential equation for dynamic control is usually written as follows:where ,, and are the displacement, speed, and acceleration vectors of a system, respectively; is a column vector of the external force; and ,, and are the mass, damping, and stiffness matrices.

Based on the Newmark- method, suppose the acceleration between *t*_{n} and *t*_{n+1} is the constant , and , where is the acceleration at *t*_{n} and is the acceleration at *t*_{n+1}. To obtain a stable and high-precision algorithm, can be expressed aswhere is the controlled variable, and . By integrating from time *t*_{n} to *t*_{n+1}, you can get the displacement at time *t*_{n+1} as follows:where and are the displacement and speed at *t*_{n}.

After equation (2) is substituted into equation (3), the following formulae can be obtained:

The Taylor series of the acceleration at *t*_{n+1} can be expanded into the following formula:where refers to the rate of acceleration change at *t*_{n}.

Equation (5) can be substituted into equation (4) as

The local absolute error of the time increment from *t*_{n} to *t*_{n+1} is defined as a difference in value between an approximate solution based on the Newmark- method and an exact solution and can be determined as

Considering that absolute error fails to reflect the degree of deviation, and because the related evaluation criteria cannot be unified, relative errors are used to carry out the estimations in this study. The relevant relative error is the specific value between the absolute error and the exact solution . The local relative error is given as follows:

An exact solution is needed to determine the relative error. However, it is quite difficult to obtain such an exact solution in general cases. Here, the Taylor series expansion method is employed to determine , as a solution to the displacement at *t*_{n+1}; this solution is then further utilized to substitute the exact solution for error estimation. Finally, the corresponding expression is as shown:

After equations (6), (7), (10), and (11) are substituted into equation (9), the following formulae can be obtained:

The upper bound on relative errors denoted as is given to obtain a time step that satisfies such a limiting value. If the time step is , items of infinitesimally higher order, and , in equation (12) can be omitted. Equation (12) can then be expressed as

#### 3. Evaluation Method for Computational Efficiency Improvement

In general, analytical solutions are difficult to obtain, so accurate relative errors cannot be obtained, and the calculation efficiency cannot be effectively compared. Thus, the degree of improvement of the calculation efficiency cannot be evaluated. The theoretical basis of the time step optimization method can solve this problem.

When a single node is selected as the research object, all time steps correspond to the same relative error limit in the time step optimization process. The minimum time step after optimization is taken as the fixed time step before optimization, so the relative error corresponding to each fixed time step is not higher than . Since the calculation accuracy is determined by the worst case of the calculation result, the calculation accuracy before and after time step optimization is consistent. On this basis, the calculation time before and after optimization are compared, and an effective comparison of the calculation efficiency is achieved.

##### 3.1. Bernoulli–Euler Beam Example

To compare the improvements in computational efficiency, an example with an analytic solution was selected. This is simply a supported beam with homogeneous uniform sections. Its initial displacement and original speed are both zero. In addition, a concentrated force load traverses the length of the beam at a uniform speed. Figure 1 shows the corresponding numerical model.

Under a moving load, the differential equation of movement for the Bernoulli–Euler beam with homogeneous uniform sections can be written as follows [17]:where and are the flexural rigidity and mass per unit length of the beam, respectively; is the damping coefficient; , , and are the lateral displacement response of the beam at *t* and the position and size of its moving load, respectively; and is the Dirac function. Thus, the relevant analytical solution is [18]

The beam is divided into 20 equal units. All calculation parameters are assumed to be dimensionless (among which and are the beam length and unit length, respectively).

In this paper, the steps to verify the evaluation method are as follows:(1)According to equation (14), the optimization time steps are determined, and the relative errors after optimization are determined.(2)The minimum optimization time step is selected as the fixed time step before optimization and the relative errors before optimization are determined.(3)Compare the maximum relative error in step 1 with that in step 2; if they are consistent, the calculation accuracy before optimization is consistent with that after optimization.(4)Comparing calculation time on the basis of consistent calculation accuracy, the comparison of computational efficiency is achieved.

##### 3.2. Time Steps and Relative Errors after Optimization

In this paper, three representative nodes (No. 5, No. 11, and No. 17 of the Bernoulli–Euler beam) are selected as the investigation objects. Select the relative error limit and substitute it into equation (14); the optimized time steps are shown in Figure 2 and the corresponding relative errors are shown in Figure 3.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

When Node 5 is chosen as the research object, as shown in Figure 2(a), the number of points is 1843, and each point corresponds to one time step, which has a total of 1843 time steps. The minimum time step of Node 5 after time step optimization is 0.0211 and the maximum time step is 0.6569. The maximum time step is 31.13 times as many as the minimum time step.

In Figure 2(b), the number of points is 1015, the minimum time step of Node 11 after time step optimization is 0.0282, and the maximum time step is 0.4665. The maximum time step is 16.54 times as many as the minimum time step.

In Figure 2(c), the number of points is 1605, the minimum time step of Node 17 after time step optimization is 0.0223, and the maximum time step is 0.3729. The maximum time step is 16.72 times as many as the minimum time step.

It can be seen from Figure 3 that in the process of time step optimization, the true relative error obtained by most of the time steps is not equal to the relative error limit but is very close. There are a few points with different degrees of deviation because of the change in the external load, the change in the node deflection, and other factors that affect the time step to varying degrees, thus affecting the relative error.

##### 3.3. Relative Errors before Optimization

According to the evaluation method proposed in this paper, the minimum optimization time step is selected as the fixed time step before optimization. That is, the minimum time Step 0.0211 is the fixed time step before optimization of Node 5，the minimum time Step 0.0282 is the fixed time step of Node 11, and the minimum time Step 0.0223 is the fixed time step of Node 17. And the corresponding relative errors before optimization are shown in Figure 4. The maximum relative error in Figure 4(a) is 0.0118, 0.0124 in Figure 4(b), and 0.0135 in Figure 4(c).

**(a)**

**(b)**

**(c)**

##### 3.4. Comparison of Calculation Efficiency

It can be seen from the comparison between Figures 3 and 4 that the maximum relative error before and after optimization has a small deviation. In this paper, it is considered that the deviation of the maximum relative error is acceptable, and the computational accuracy is consistent before and after optimization. Therefore, on this basis, the comparison result of the calculation efficiency can be obtained by comparing the calculation time.

When the calculation process of the program is timed, the computation time of different nodes before and after optimization are listed in Table 1.

The average computation time of three times for Node 5 is 14.376 s, which is calculated at a fixed time step of 0.0211 before optimization. And the average computation time of three times is 6.587 s after optimization, which reduces the computation time by 54.18% and improves the calculation efficiency by 54.18%. Similarly, the calculation efficiency of Node 11 is increased by 72.35%, and the calculation efficiency of Node 17 is increased by 59.31%.

To sum up, for a numerical calculation from which the analytical solution cannot be obtained, the following method can be used to evaluate the improvement degree of calculation efficiency: the minimum time step in the time step optimization process is used as the fixed time step before optimization in order to ensure that the calculation efficiency before and after optimization is consistent. On this basis, both calculation time can be compared to obtain the improvement degree of the calculation efficiency after adopting the time step optimization method.

#### 4. Application of Time Step Optimization Method in Seismic Liquefaction Analysis

##### 4.1. Numerical Method for Dynamic Liquefaction of Saturated Soil

Based on the theories of Biot porous media elastic wave propagation and two-phase consolidation [19–21], Tamura and Akai proposed a FEM–FDM (finite element method–finite difference method) coupling scheme in 1978 [22]. Subsequently, the scheme was extended by Oka et al. to the analysis of dynamic liquefaction of saturated soil based on the assumption of infinitesimal strain [23]. Finally, the system is established by Zienkiewicz et al. [24]. This numerical scheme is discretized by the finite element method in the spatial domain and by the Newmark method in the temporal domain.

###### 4.1.1. Solid–Fluid Governing Equations

The mixture should satisfy the mass balance and momentum balance equations. The equilibrium equation for the solid–fluid mixture can be derived aswhere is the mixture’s density; is the total stress tensor; is the body force vector; and is the displacement vector of the soil skeleton.

The continuity equation for pore water is derived aswhere is the porosity; is the fluid bulk modulus; is the pore water pressure; is the coefficient of permeability; is the volumetric strain of solid; and is the fluid intrinsic density.

The cyclic elastoplastic constitutive model [25] is used to establish the relationship between stress increment and strain increment as shown in the following formula:where is the strain increment; is the elastoplastic stiffness matrix of the material; and is the effective stress increment.

By adopting the FEM–FDM coupling scheme, the equilibrium equation and continuity equation for one element can be written aswhere , and are the mass, damping, and stiffness matrix; is the internal force vector; is the external force vector; is the modulus vector for the pore water pressure; and and are coefficients for the items related to the pore water. Other variable parameters in equation (21) are shown in the literature [13].

###### 4.1.2. Temporal Discretization

For the soil skeletons, the displacement, velocity, and acceleration at time are represented according to the variables at time by the Newmark- method as follows:where and are controlling parameters for the Newmark- method.

The pore water pressure at time is represented in terms of the variables at time by the backward difference method as follows:

Upon substituting equations (22)–(25) into governing equation (21) of the form, the final governing equation is obtained as

In equation (26), variable parameters are shown in literature [13].

##### 4.2. Earth Dam Example

In this paper, an example of earth dam under earthquake load is selected as the research object. The calculation model and meshing of earth dam are shown in Figures 5 and 6 [14]. A 4-node quadrilateral element is used, and there are 1186 nodes and 1100 elements in total. 50 m super-long elements are added at each horizontal side to reduce the boundary effect caused by the horizontal vibration. The boundary conditions are set as follows: the vertical and horizontal directions of the bottom nodes are fixed; the left and right nodes are fixed in horizontal direction and free in vertical direction; and the top surface is the drainage boundary. In the model, the dam crest node N1 and the dam toe node N2 are taken as the calculated output nodes.

The cyclic elastoplastic constitutive model is adopted for soil mass. Soil parameters are shown in Table 2. The input horizontal seismic wave is shown in Figure 7.

##### 4.3. Calculation Results after Time Step Optimization

The finite element program for calculating seismic liquefaction of earth dam was optimized by the time step optimization method, and the earth dam model was calculated by the optimized program. The calculation process is as follows: select the initial fixed time step and calculate the first two steps to obtain the initial variable value (, and ); the rate of acceleration change is determined by the following equation:set the relative error limit and calculate the optimization time step by equation (14); then substitute into the optimized program to determine , and ; then calculate the optimization time step by equation (14); and so on.

When determining the optimization time step, it is necessary to calculate the time steps of the two components of the coordinate axis, respectively, because the earth dam example is a two-dimensional model, and take the minimum time step as the final optimization time step. The calculation results of the optimization time steps are shown in Figure 8: there are a total of 3026 points, the minimum time step is , the maximum time step is 0.144s, and the maximum time step is 960 times the minimum time step.

The vertical displacement time history curve of the dam crest is shown in Figure 9, the horizontal displacement time history curve of the dam toe is shown in Figure 10, the distribution of excess pore water pressure is shown in Figure 11, and the distribution of vertical displacement is shown in Figure 12.

**(a)**

**(b)**

##### 4.4. Calculation Results before Time Step Optimization

According to the evaluation method for computational efficiency improvement proposed in this paper, the minimum optimization time step is selected as the fixed time step before optimization, so as to ensure the consistency of the calculation accuracy before and after optimization. The minimum time step in Figure 8 is selected as the fixed time step before optimization. The vertical displacement time history curve of the dam crest is shown in Figure 9, the horizontal displacement time history curve of the dam toe is shown in Figure 10, the distribution of the excess pore water pressure is shown in Figure 13, and the distribution of vertical displacement is shown in Figure 14.

As shown in Figures 9 and 10, the vertical displacement of the dam crest and the horizontal displacement of the dam toe are basically consistent before and after the time step optimization. The distribution of excess pore water pressure (Figures 11 and 13) and the vertical displacement distribution (Figures 12 and 14) of earth dam before and after the time step optimization are also basically consistent. This also shows that the calculation accuracy is basically consistent before and after the time step optimization, which further validates the evaluation method proposed in this paper.

##### 4.5. Evaluation of Computational Efficiency after Time Step Optimization

According to the evaluation method proposed in this paper, the minimum optimization time step is used as the fixed time step before optimization, which can guarantee the consistency of the calculation accuracy. On this basis, the evaluation of computing efficiency can be achieved by comparing the computing time.

On the IBM Intel(R) Xeon(R) 3.00 GHz server, the calculation time before optimization is 2716s, and the calculation time after optimization is 675 s. The calculation efficiency is improved by

Due to the drastic change in the external load, the time steps vary in a very large range: the minimum time step is , the maximum time step is 0.144 s, and the maximum time step is 960 times the minimum time step. Therefore, the improvement of computational efficiency is obviously higher than that of Bernoulli–Euler beam, which indicates that the time step optimization method is very suitable for the dynamic numerical calculation of impact load, vibration load, explosion load, or earthquake load.

#### 5. Future Research

From an academic point of view, it is very important to know how much computational efficiency has been improved by using the time step optimization method. However, from an application point of view, when processes with significant computational burden may require several days, weeks, or even longer to complete the calculations, using the proposed method to evaluate the improvement of computational efficiency will inevitably waste enormous computation time, which is not worth it. Therefore, when applying the time step optimization method, the improvement of computational efficiency can be obtained without additional calculation at the end of optimization, which is the next step to be studied.

By comparing Figures 7 and 8, it can be seen that the time steps are obviously affected by the earthquake load: the load changes relatively gently at the early stage (0-1 s) and the late stage (20–22 s), and the time steps are larger at this time; when the load changes dramatically (2–13 s), the time steps are very small. Therefore, it is necessary to study the correlation between optimization time steps and external load.

#### 6. Conclusion

In the dynamic numerical calculation with a fixed time step, especially when dealing with the impact load, vibration load, explosion load, or seismic load, the calculation result would be divergent and even the calculation would be interrupted due to the improper setting of time step. In this paper, the time step optimization method was introduced into seismic liquefaction numerical analysis platform to solve the time step setting problem successfully while realizing the time step optimization. Therefore, this paper provides a useful reference to solve the time step setting problem of dynamic numerical calculation.

In addition, for the case in which the analytical solution could not be obtained, a new method was proposed to evaluate the calculation efficiency improvement. And the proposed method was verified by an example with an analytical solution. Finally, the seismic response computational efficiency after time step optimization was evaluated with the proposed method. The results show that the computing efficiency is greatly improved after the seismic liquefaction numerical analysis program is optimized by the time step optimization method.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The work presented in this paper was part of the research sponsored by the Key Program of National Natural Science Foundation of China (Grant no. 41672258), the Natural Foundation of Shandong Province (Grant no. ZR2018LE017), and the Science and Technology Development Plan of Taian City (Grant no. 2019NS075).