Application of the Time Step Optimization Method in Seismic Liquefaction Analysis and Its Evaluation Method for Computational Efficiency Improvement

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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:$$\begin{array}{}\text{(1)}& \mathbf{M}\ddot{u}t+\mathbf{C}\dot{u}t+\mathbf{K}ut=\mathbf{F}t,\end{array}$$where $ut$,$\dot{u}t$, and $\ddot{u}t$ are the displacement, speed, and acceleration vectors of a system, respectively; $\mathbf{F}t$ is a column vector of the external force; and $\mathbf{M}$,$\mathbf{C}$, and $\mathbf{K}$ are the mass, damping, and stiffness matrices.

Based on the Newmark-$\beta$ method, suppose the acceleration between t_n and t_n+1 is the constant $a$, and ${\ddot{\overline{y}}}_n\le a\le {\ddot{\overline{y}}}_{n+1}$, where ${\ddot{\overline{y}}}_n$ is the acceleration at t_n and ${\ddot{\overline{y}}}_{n+1}$ is the acceleration at t_n+1. To obtain a stable and high-precision algorithm, $a$ can be expressed as$$\begin{array}{}\text{(2)}& a=1-2\beta {\ddot{\overline{y}}}_n+2\beta {\ddot{\overline{y}}}_{n+1},\end{array}$$where $\beta$ is the controlled variable, and $0\le \beta \le 0.5$. By integrating $a$ from time t_n to t_n+1, you can get the displacement ${\overline{y}}_{n+1}$ at time t_n+1 as follows:$$\begin{array}{}\text{(3)}& {\overline{y}}_{n+1}={\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5\Delta t^{2}a,\end{array}$$where${\overline{y}}_n$ and ${\dot{\overline{y}}}_n$ are the displacement and speed at t_n.

After equation (2) is substituted into equation (3), the following formulae can be obtained:$$\begin{array}{}\text{(4)}& {\overline{y}}_{n+1}={\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5-\beta {\ddot{\overline{y}}}_n\Delta t^2+\beta {\ddot{\overline{y}}}_{n+1}\Delta t^2.\end{array}$$

The Taylor series of the acceleration at t_n+1 can be expanded into the following formula:$$\begin{array}{}\text{(5)}& {\ddot{\overline{y}}}_{n+1}={\ddot{\overline{y}}}_n+{\ddot{\overline{y}}}_n\Delta t+{\rm O}\Delta t,\end{array}$$where ${\ddot{\overline{y}}}_n$ refers to the rate of acceleration change at t_n.

Equation (5) can be substituted into equation (4) as$$\begin{array}{}\text{(6)}& {\overline{y}}_{n+1}={\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5{\ddot{\overline{y}}}_n\Delta t^2+\beta {\stackrel{⃛}{\overline{y}}}_n\Delta t^3+{\rm O}\Delta t^3.\text{(7)}& {\rm O}\Delta t^3=\beta \frac{1}{2}{\overline{y}}_n^4\Delta t^4+\cdots +\frac{1}{m-2!}{\overline{y}}_n^m\Delta t^m+\cdots .\end{array}$$

The local absolute error $e$ of the time increment $\mathrm{\Delta }t$ from t_n to t_n+1 is defined as a difference in value between an approximate solution ${\overline{y}}_{n+1}$ based on the Newmark-$\beta$ method and an exact solution ${\overline{y}}_{n+1}^{\ast }$ and can be determined as$$\begin{array}{}\text{(8)}& e={\overline{y}}_{n+1}-{\overline{y}}_{n+1}^{\ast }.\end{array}$$

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 $e$ and the exact solution ${\overline{y}}_{n+1}^{\ast }$. The local relative error $\eta$ is given as follows:$$\begin{array}{}\text{(9)}& \eta =\frac{e}{{\overline{y}}_{n+1}^{\ast }}=\frac{{\overline{y}}_{n+1}-{\overline{y}}_{n+1}^{\ast }}{{\overline{y}}_{n+1}^{\ast }}.\end{array}$$

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 ${\overline{y}}_{n+1}^{\ast }$, 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:$$\begin{array}{}\text{(10)}& {\overline{y}}_{n+1}^{\ast }={\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5{\ddot{\overline{y}}}_n\Delta t^2+\frac{1}{6}{\ddot{\overline{y}}}_n\Delta t^3+{\rm O}{\Delta t^3}^{\prime }.\text{(11)}& {\rm O}{\Delta t^3}^{\prime }=\frac{1}{4!}{\overline{y}}_n^4\Delta t^4+\cdots +\frac{1}{m!}{\overline{y}}_n^m\Delta t^m+\cdots .\end{array}$$

After equations (6), (7), (10), and (11) are substituted into equation (9), the following formulae can be obtained:$$\begin{array}{}\text{(12)}& \eta =\frac{{\overline{y}}_{n+1}-{\overline{y}}_{n+1}^{\ast }}{{\overline{y}}_{n+1}^{\ast }}=\frac{\beta -1/6{\ddot{\overline{y}}}_n\Delta t^3+{\rm O}{\Delta t^3}^{″}}{{\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5{\ddot{\overline{y}}}_n\Delta t^2+1/6{\ddot{\overline{y}}}_n\Delta t^3+{\rm O}{\Delta t^3}^{\prime }}.\text{(13)}& {\rm O}{\Delta t^3}^{″}=\frac{1}{2}\beta -\frac{1}{4!}{\overline{y}}_n^4\Delta t^4+\cdots +\frac{1}{m-2!}\beta -\frac{1}{m!}{\overline{y}}_n^m\Delta t^m+\cdots .\end{array}$$

The upper bound on relative errors denoted as $\overline{\eta }$ is given to obtain a time step that satisfies such a limiting value. If the time step is $\mathrm{\Delta }t\ll 1$, items of infinitesimally higher order, $\mathrm{{\rm O}}{\mathrm{\Delta }{t}^3}^{\prime }$ and $\mathrm{{\rm O}}{\mathrm{\Delta }{t}^3}^{″}$, in equation (12) can be omitted. Equation (12) can then be expressed as$$\begin{array}{}\text{(14)}& \overline{\eta }=\frac{1/6-\beta {\ddot{\overline{y}}}_n\Delta t^3}{{\overline{y}}_n+{\dot{\overline{y}}}_n\Delta t+0.5{\ddot{\overline{y}}}_n\Delta t^2+1/6{\ddot{\overline{y}}}_n\Delta t^3}.\end{array}$$

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 $\overline{\eta }$ 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 $\overline{\eta }$. 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.

[figure omitted; refer to PDF]

Under a moving load, the differential equation of movement for the Bernoulli–Euler beam with homogeneous uniform sections can be written as follows [17]:$$\begin{array}{}\text{(15)}& \text{EI}x,t\frac{{\partial }^{4}yx,t}{{\partial }^{4}x}+\overline{m}x,t\frac{{\partial }^{2}yx,t}{{\partial }^{2}t}+cx,t\frac{\partial yx,t}{\partial t}=\delta x-x_{p}tPt,\end{array}$$where $\text{EI}x,t$ and $\overline{m}x,t$ are the flexural rigidity and mass per unit length of the beam, respectively; $cx,t$ is the damping coefficient; $yx,t$, $x_{p}t$, and $Pt$ are the lateral displacement response of the beam at t and the position and size of its moving load, respectively; and $\delta \cdot$ is the Dirac function. Thus, the relevant analytical solution is [18]$$\begin{array}{}\text{(16)}& yx,t=-{\displaystyle \sum _{i=1}^{\infty }\frac{2P}{\overline{m}{\omega }_i}\frac{\sin i\pi /Lx}{{i\pi v/L}^2-{\omega }_i^2}}{\omega }_i\sin \frac{i\pi v}{L}t-\frac{i\pi v}{L}\sin {\omega }_{i}t,\hspace{1em}0\le t\le \frac{L}{v},{\omega }_i={\frac{i\pi }{L}}^2\sqrt{\frac{\text{EI}}{\overline{m}}}.\end{array}$$

The beam is divided into 20 equal units. All calculation parameters are assumed to be dimensionless (among which $L$ and $l$ are the beam length and unit length, respectively).$$\begin{array}{c}\text{(17)}& EI=100,c=0.0,\\ L=10,\\ \overline{m}=1.0,\\ P=1.0,\\ v=\mathrm{0.1.}\end{array}$$

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 $\overline{\eta }=0.01$ 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.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

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 $\overline{\eta }=0.01$ 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).

[figures omitted; refer to PDF]

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.

Table 1

Comparison of computation time.

		First computation time (s)	Second computation time (s)	Third computation time (s)	Average computation time (s)
Node 5	Before optimization	14.375	14.388	14.366	14.376
	After optimization	6.592	6.595	6.579	6.587
Node 11	Before optimization	14.003	13.990	13.876	13.956
	After optimization	3.832	3.873	3.871	3.859
Node 17	Before optimization	14.169	14.149	14.184	14.167
	After optimization	5.769	5.758	5.766	5.764

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 $u-p$ 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 as$$\begin{array}{}\text{(18)}& \frac{\partial {\sigma }_{ij}}{\partial x_j}+\rho b_i=\rho {\ddot{u}}_i^s,\end{array}$$where $\rho$ is the mixture’s density; ${\sigma }_{ij}$ is the total stress tensor; $b_i$ is the body force vector; and ${\ddot{u}}_i^s$ is the displacement vector of the soil skeleton.

The continuity equation for pore water is derived as$$\begin{array}{}\text{(19)}& \frac{n}{K^w}\frac{\text{d}{p}^w}{\text{d}t}-\frac{\text{nk}}{{\gamma }_w}\frac{{\partial }^2{p}^w}{\partial x_i^2}-{\rho }^w{\ddot{\epsilon }}_{ii}^s-{\dot{\epsilon }}_{ii}^s=0,\end{array}$$where $n$ is the porosity; $K^w$ is the fluid bulk modulus; $p^w$ is the pore water pressure; $k$ is the coefficient of permeability; ${\epsilon }_{ii}^s$ is the volumetric strain of solid; and ${\rho }^w$ 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:$$\begin{array}{}\text{(20)}& d{\sigma }_{ij}^{\prime }=D_{\text{ijkl}}^{\text{ep}}d{\epsilon }_{\text{kl}}^s,\end{array}$$where $d{\epsilon }_{\text{kl}}^s$ is the strain increment; $D_{\text{ijkl}}^{\text{ep}}$ is the elastoplastic stiffness matrix of the material; and $d{\sigma }_{ij}^{\prime }$ is the effective stress increment.

By adopting the FEM–FDM coupling scheme, the equilibrium equation and continuity equation for one element can be written as$$\begin{array}{}\text{(21)}& \begin{array}{c}M\cdot {\ddot{\mathbf{u}}}_N+C\cdot {\dot{\mathbf{u}}}_N+K\cdot \Delta {\mathbf{u}}_N+{\mathbf{K}}_v\cdot p^w=\mathbf{F}-{\mathbf{R}}_t,\\ {\rho }^w{\mathbf{K}}_v^T\cdot {\ddot{\mathbf{u}}}_N-\frac{{\gamma }_w}{nk}{\mathbf{K}}_v^T\cdot {\dot{\mathbf{u}}}_N+A\cdot {\dot{p}}^w-{\displaystyle \sum _{i=1}^m{\alpha }_i}{p}_i^w+\alpha p^W=0,\end{array}\end{array}$$where $M,C$, and $K$ are the mass, damping, and stiffness matrix; ${\mathbf{R}}_t$ is the internal force vector; $\mathbf{F}$ is the external force vector; ${\mathbf{K}}_v$ is the modulus vector for the pore water pressure; and $\alpha$ and $A$ 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 $t+\mathrm{\Delta }t$ are represented according to the variables at time $t$ by the Newmark-$\beta$ method as follows:$$\begin{array}{}\text{(22)}& {\ddot{\mathbf{u}}}_{t+\Delta t}={\ddot{\mathbf{u}}}_t+\Delta \ddot{\mathbf{u}},\text{(23)}& {\dot{\mathbf{u}}}_{t+\Delta t}={\dot{\mathbf{u}}}_t+1-\gamma {\ddot{\mathbf{u}}}_t+\gamma {\ddot{\mathbf{u}}}_{t+\Delta t}\times \Delta t,\text{(24)}& {\mathbf{u}}_{t+\Delta t}={\mathbf{u}}_t+{\dot{\mathbf{u}}}_t\times \Delta t+\frac{1}{2}1-2\beta {\ddot{\mathbf{u}}}_t+2\beta {\ddot{\mathbf{u}}}_{t+\Delta t}\times \Delta t^2,\end{array}$$where $\gamma$ and $\beta$ are controlling parameters for the Newmark-$\beta$ method.

The pore water pressure at time $t+\mathrm{\Delta }t$ is represented in terms of the variables at time $t$ by the backward difference method as follows:$$\begin{array}{}\text{(25)}& p_{t+\Delta t}^w=p_t^w+{\dot{p}}_{t+\Delta t}^w\cdot \Delta t.\end{array}$$

Upon substituting equations (22)–(25) into governing equation (21) of the $u-p$ form, the final governing equation is obtained as$$\begin{array}{c}\text{(26)}& \begin{array}{cc}M+\Delta t\gamma C_{t+\Delta t}+\Delta t^2\beta K_{t+\Delta t}\hfill & \Delta t{K}_v\hfill \\ {\rho }^w-\frac{{\gamma }_w\Delta t\gamma }{nk}{K}_v^T\hfill & A+\alpha \Delta t\hfill \end{array}\begin{array}{c}{\ddot{\mathbf{u}}}_{t+\Delta t}\\ {\dot{p}}_{t+\Delta t}^w\end{array}+\begin{array}{c}0\\ \Delta t{\displaystyle \sum _{i=1}^m{\alpha }_i}{\dot{p}}_{i,t+\Delta t}^w\end{array}={\mathbf{F}}_{t+\Delta t}-{\mathbf{R}}_t-C_{t+\Delta t}{\dot{\mathbf{u}}}_t+\Delta t1-\gamma {\ddot{\mathbf{u}}}_t-K_{t+\Delta t}\Delta t{\dot{\mathbf{u}}}_t+\Delta t^2\frac{1}{2}-\beta {\ddot{\mathbf{u}}}_t-K_v{p}_t^w\frac{{\gamma }_w}{nk}{K}_v^T{\dot{\mathbf{u}}}_t+\Delta t1-\gamma {\ddot{\mathbf{u}}}_t+{\displaystyle \sum _{i=1}^m{\alpha }_i{p}_{i,t}^w-\alpha p_t^w}.\end{array}$$

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.

[figure omitted; refer to PDF]

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.

Table 2

Soil properties of 2D earth dam.

Soil parameters	Saturated soil	Dry soil
Initial void ratio, e0	0.80	0.75
Compression index, $\lambda$	0.015	0.018
Density, $\mathrm{\Delta }t$(t/m3)	1.951	1.85
Swelling index	0.0045	0.0055
Initial shear modulus ration, $G_0/G_{m0}^{\prime }$	930.0	930.0
Permeability coefficient, $km/s$	$2.6\times {10}^{-5}$
Phase transformation stress ratio, $M_m$	0.91	0.91
Failure stress ratio, $M_f$	1.12	1.12
Hardening parameter $B_0$	2500	3000
Hardening parameter $B_1$	30	60
Bulk modulus of pore water, $K_w$(N/m2)	$2.0\times {10}^6$

[figure omitted; refer to PDF]

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 $\mathrm{\Delta }{t}_0=1.0\times {10}^{-4}\text{s}$ and calculate the first two steps to obtain the initial variable value ($u_{t_2},{\dot{u}}_{t_2},{\ddot{u}}_{t_2}$, and ${\ddot{u}}_{t_1}$); the rate of acceleration change is determined by the following equation:$$\begin{array}{}\text{(27)}& {\stackrel{\dots }{u}}_{t_2}=\frac{{\ddot{u}}_{t_2}-{\ddot{u}}_{t_1}}{\Delta t_0},\end{array}$$set the relative error limit $\overline{\eta }=1.0\times {10}^{-4}$ and calculate the optimization time step $\mathrm{\Delta }{t}_1$ by equation (14); then substitute $\mathrm{\Delta }{t}_1$ into the optimized program to determine $u_{t_3},{\dot{u}}_{t_3},{\ddot{u}}_{t_3}$, and ${\stackrel{⃛}{u}}_{t_3}$; then calculate the optimization time step $\mathrm{\Delta }{t}_2$ 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 $1.5\times {10}^{-3}$, the maximum time step is 0.144s, and the maximum time step is 960 times the minimum time step.

[figure omitted; refer to PDF]

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.

[figure omitted; refer to PDF]

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$$\begin{array}{}\text{(28)}& \frac{2716-675}{2716}\times 100\%=75.15\%.\end{array}$$

Due to the drastic change in the external load, the time steps vary in a very large range: the minimum time step is $1.5\times {10}^{-3}\text{s}$, 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.

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).

References

[1] O. C. Zienkiewicz, Y. M. Xie, "A simple error estimator and adaptive time stepping procedure for dynamic analysis," Earthquake Engineering & Structural Dynamics, vol. 20 no. 9, pp. 871-887, DOI: 10.1002/eqe.4290200907, 1991.

[2] L. F. Zeng, N.-E. Wiberg, X. D. Li, Y. M. Xie, "A posteriori local error estimation and adaptive time-stepping for newmark integration in dynamic analysis," Earthquake Engineering & Structural Dynamics, vol. 21 no. 7, pp. 555-571, DOI: 10.1002/eqe.4290210701, 1992.

[3] N.-E. Wiberg, X. D. Li, "A post-processing technique and ana posteriori error estimate for the newmark method in dynamic analysis," Earthquake Engineering & Structural Dynamics, vol. 22 no. 6, pp. 465-489, DOI: 10.1002/eqe.4290220602, 1993.

[4] S. W. Sloan, A. J. Abbo, "Biot consolidation analysis with automatic time stepping and error control part 2: applications," International Journal for Numerical and Analytical Methods in Geomechanics, vol. 23 no. 6, pp. 493-529, DOI: 10.1002/(sici)1096-9853(199905)23:6<493::aid-nag950>3.0.co;2-p, 1999.

[5] D. Kavetski, P. Binning, S. W. Sloan, "Adaptive backward Euler time stepping with truncation error control for numerical modelling of unsaturated fluid flow fluid flow," International Journal for Numerical Methods in Engineering, vol. 53 no. 6, pp. 1301-1322, DOI: 10.1002/nme.329, 2002.

[6] G. Tang, A. N. Alshawabkeh, M. A. Mayes, "Automatic time stepping with global error control for groundwater flow models," Journal of Hydrologic Engineering, vol. 13 no. 9, pp. 803-810, DOI: 10.1061/(asce)1084-0699(2008)13:9(803), 2008.

[7] S.-R. Kuo, J. D. Yau, Y. B. Yang, "A robust time-integration algorithm for solving nonlinear dynamic problems with large rotations and displacements," International Journal of Structural Stability and Dynamics, vol. 12 no. 6,DOI: 10.1142/s0219455412500514, 2012.

[8] K. J. Wang, J. S. Cheng, J. Q. Wang, "Numerical simulation analysis of offshore floating wind power generation platform based on time step adaptive method," Naval Architecture and Ocean Engineering, vol. 2, pp. 52-57, 2013.

[9] W. N. Mao, J. K. Liu, "Application of adaptive time step method to water-heat coupling model of soil freezing," Journal of Disaster Prevention and Mitigation Engineering, vol. 34 no. 4, pp. 510-516, 2014.

[10] N. Ahmed, V. John, "Adaptive time step control for higher order variational time discretizations applied to convection-diffusion-reaction equations," Computer Methods in Applied Mechanics and Engineering, vol. 285, pp. 83-101, DOI: 10.1016/j.cma.2014.10.054, 2015.

[11] M. Laforest, S. Prudhomme, K. G. van der Zee, "Preface to the special issue on error estimation and adaptivity for nonlinear and time-dependent problems," Computer Methods in Applied Mechanics and Engineering, vol. 288,DOI: 10.1016/j.cma.2014.12.012, 2015.

[12] V. Ziaei-Rad, Y. Shen, "Massive parallelization of the phase field formulation for crack propagation with time adaptivity," Computer Methods in Applied Mechanics and Engineering, vol. 312, pp. 224-253, DOI: 10.1016/j.cma.2016.04.013, 2016.

[13] X. W. Zhang, Study on the Adaptive Time Stepping Method in Seismic Liquefaction Analysis of Saturated sand, 2015.

[14] X. W. Tang, X. W. Zhang, R. Uzuoka, "Novel adaptive time stepping method and its application to soil seismic liquefaction analysis," Soil Dynamics and Earthquake Engineering, vol. 71, pp. 100-113, DOI: 10.1016/j.soildyn.2015.01.016, 2015.

[15] B. Li, X. Tang, G. Yang, "Time-step optimization method to improve calculation efficiency and feasibility validation reliance," Engineering Optimization, vol. 51 no. 8, pp. 1324-1335, DOI: 10.1080/0305215x.2018.1519074, 2019.

[16] B. Li, X. W. Tang, "“A time step’s automatic adjustment method with load synchronous variation," Journal of Vibration and Shock, vol. 37 no. 11, 2018.

[17] Y.-H. Lin, M. W. Trethewey, "Finite element analysis of elastic beams subjected to moving dynamic loads," Journal of Sound and Vibration, vol. 136 no. 2, pp. 323-342, DOI: 10.1016/0022-460x(90)90860-3, 1990.

[18] L. Fryba, Vibration of Solids and Structures under Moving Loads, 1999.

[19] M. A. Biot, "Theory of propagation of elastic waves in a fluid‐saturated porous solid. II. Higher frequency range," The Journal of the Acoustical Society of America, vol. 28 no. 2, pp. 179-191, DOI: 10.1121/1.1908241, 1956.

[20] M. A. Biot, "Generalized theory of acoustic propagation in porous dissipative media," The Journal of the Acoustical Society of America, vol. 34 no. 9, pp. 1254-1264, DOI: 10.1121/1.1918315, 1962.

[21] M. A. Biot, "General theory of three‐dimensional consolidation," Journal of Applied Physics, vol. 12 no. 2, pp. 155-164, DOI: 10.1063/1.1712886, 1941.

[22] K. Akai, T. Tamura, "Numerical analysis of multi-dimensional consolidation accompanied with elasto-plastic constitutive equation," Proceedings of the Japan Society of Civil Engineers, vol. 1978 no. 269, pp. 95-104, DOI: 10.2208/jscej1969.1978.95, 1978.

[23] F. Oka, A. Yashima, T. Shibata, M. Kato, R. Uzuoka, "FEM-FDM coupled liquefaction analysis of a porous soil using an elasto-plastic model," Applied Scientific Research, vol. 52 no. 3, pp. 209-245, DOI: 10.1007/bf00853951, 1994.

[24] O. C. Zienkiewicz, T. Shiomi, "Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution," International Journal for Numerical and Analytical Methods in Geomechanics, vol. 8 no. 1, pp. 71-96, DOI: 10.1002/nag.1610080106, 1984.

[25] F. Oka, A. Yashima, A. Tateishi, Y. Taguchi, A. Yamashita, "A cyclic elasto-plastic constitutive model for sand considering a plastic-strain dependence of the shear modulus," Géotechnique, vol. 49 no. 5, pp. 661-680, DOI: 10.1680/geot.1999.49.5.661, 1999.