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

Authors: Qi Li; Bin Li; Shaorui Sun; Nan Yang
Journal: Shock and Vibration
 (2020)

---

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


---

## Body

## 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:(1)Mu¨t+Cu˙t+Kut=Ft,where ut,u˙t, and u¨t are the displacement, speed, and acceleration vectors of a system, respectively; Ft is a column vector of the external force; and M,C, and K are the mass, damping, and stiffness matrices.Based on the Newmark-β method, suppose the acceleration between tn and tn+1 is the constant a, and y¯¨n≤a≤y¯¨n+1, where y¯¨n is the acceleration at tn and y¯¨n+1 is the acceleration at tn+1. To obtain a stable and high-precision algorithm, a can be expressed as(2)a=1−2βy¯¨n+2βy¯¨n+1,where β is the controlled variable, and 0≤β≤0.5. By integrating a from time tn to tn+1, you can get the displacement y¯n+1 at time tn+1 as follows:(3)y¯n+1=y¯n+y¯˙nΔt+0.5Δt2a,wherey¯n and y¯˙n are the displacement and speed at tn.After equation (2) is substituted into equation (3), the following formulae can be obtained:(4)y¯n+1=y¯n+y¯˙nΔt+0.5−βy¯¨nΔt2+βy¯¨n+1Δt2.The Taylor series of the acceleration attn+1 can be expanded into the following formula:(5)y¯¨n+1=y¯¨n+y¯¨nΔt+ΟΔt,where y¯¨n refers to the rate of acceleration change at tn.Equation (5) can be substituted into equation (4) as(6)y¯n+1=y¯n+y¯˙nΔt+0.5y¯¨nΔt2+βy¯⃛nΔt3+ΟΔt3.(7)ΟΔt3=β12y¯n4Δt4+⋯+1m−2!y¯nmΔtm+⋯.The local absolute errore of the time increment Δt from tn to tn+1 is defined as a difference in value between an approximate solution y¯n+1 based on the Newmark-β method and an exact solution y¯n+1∗ and can be determined as(8)e=y¯n+1−y¯n+1∗.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 errore and the exact solution y¯n+1∗. The local relative error η is given as follows:(9)η=ey¯n+1∗=y¯n+1−y¯n+1∗y¯n+1∗.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 determiney¯n+1∗, as a solution to the displacement at tn+1; this solution is then further utilized to substitute the exact solution for error estimation. Finally, the corresponding expression is as shown:(10)y¯n+1∗=y¯n+y¯˙nΔt+0.5y¯¨nΔt2+16y¯¨nΔt3+ΟΔt3′.(11)ΟΔt3′=14!y¯n4Δt4+⋯+1m!y¯nmΔtm+⋯.After equations (6), (7), (10), and (11) are substituted into equation (9), the following formulae can be obtained:(12)η=y¯n+1−y¯n+1∗y¯n+1∗=β−1/6y¯¨nΔt3+ΟΔt3″y¯n+y¯˙nΔt+0.5y¯¨nΔt2+1/6y¯¨nΔt3+ΟΔt3′.(13)ΟΔt3″=12β−14!y¯n4Δt4+⋯+1m−2!β−1m!y¯nmΔtm+⋯.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 Δt≪1, items of infinitesimally higher order, ΟΔt3′ and ΟΔt3″, in equation (12) can be omitted. Equation (12) can then be expressed as(14)η¯=1/6−βy¯¨nΔt3y¯n+y¯˙nΔt+0.5y¯¨nΔt2+1/6y¯¨nΔt3.


## 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. Figure1 shows the corresponding numerical model.Figure 1
Simply supported beam under moving load and node distribution.Under a moving load, the differential equation of movement for the Bernoulli–Euler beam with homogeneous uniform sections can be written as follows [17]:(15)EIx,t∂4yx,t∂4x+m¯x,t∂2yx,t∂2t+cx,t∂yx,t∂t=δx−xptPt,where EIx,t and m¯x,t are the flexural rigidity and mass per unit length of the beam, respectively; cx,t is the damping coefficient; yx,t, xpt, and Pt 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](16)yx,t=−∑i=1∞2Pm¯ωisiniπ/Lxiπv/L2−ωi2ωisiniπvLt−iπvLsinωit,0≤t≤Lv,ωi=iπL2EIm¯.The beam is divided into 20 equal units. All calculation parameters are assumed to be dimensionless (among whichL and l are the beam length and unit length, respectively).(17)EI=100,c=0.0,L=10,m¯=1.0,P=1.0,v=0.1.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η¯=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.Figure 2
Optimization time steps variation over time. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)Figure 3
Relative errors after time step optimization. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)When Node 5 is chosen as the research object, as shown in Figure2(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 Figure2(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 Figure2(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 Figure3 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 η¯=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 Figure4. The maximum relative error in Figure 4(a) is 0.0118, 0.0124 in Figure 4(b), and 0.0135 in Figure 4(c).Figure 4
Relative errors before time step optimization. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)


### 3.4. Comparison of Calculation Efficiency

It can be seen from the comparison between Figures3 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 Table1.Table 1
Comparison of computation time.
First computation time (s)Second computation time (s)Third computation time (s)Average computation time (s)Node 5Before optimization14.37514.38814.36614.376After optimization6.5926.5956.5796.587Node 11Before optimization14.00313.99013.87613.956After optimization3.8323.8733.8713.859Node 17Before optimization14.16914.14914.18414.167After optimization5.7695.7585.7665.764The 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.


## 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. Figure1 shows the corresponding numerical model.Figure 1
Simply supported beam under moving load and node distribution.Under a moving load, the differential equation of movement for the Bernoulli–Euler beam with homogeneous uniform sections can be written as follows [17]:(15)EIx,t∂4yx,t∂4x+m¯x,t∂2yx,t∂2t+cx,t∂yx,t∂t=δx−xptPt,where EIx,t and m¯x,t are the flexural rigidity and mass per unit length of the beam, respectively; cx,t is the damping coefficient; yx,t, xpt, and Pt 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](16)yx,t=−∑i=1∞2Pm¯ωisiniπ/Lxiπv/L2−ωi2ωisiniπvLt−iπvLsinωit,0≤t≤Lv,ωi=iπL2EIm¯.The beam is divided into 20 equal units. All calculation parameters are assumed to be dimensionless (among whichL and l are the beam length and unit length, respectively).(17)EI=100,c=0.0,L=10,m¯=1.0,P=1.0,v=0.1.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η¯=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.Figure 2
Optimization time steps variation over time. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)Figure 3
Relative errors after time step optimization. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)When Node 5 is chosen as the research object, as shown in Figure2(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 Figure2(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 Figure2(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 Figure3 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 η¯=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 Figure4. The maximum relative error in Figure 4(a) is 0.0118, 0.0124 in Figure 4(b), and 0.0135 in Figure 4(c).Figure 4
Relative errors before time step optimization. (a) Node 5. (b) Node 11. (c) Node 17.
(a)(b)(c)


## 3.4. Comparison of Calculation Efficiency

It can be seen from the comparison between Figures3 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 Table1.Table 1
Comparison of computation time.
First computation time (s)Second computation time (s)Third computation time (s)Average computation time (s)Node 5Before optimization14.37514.38814.36614.376After optimization6.5926.5956.5796.587Node 11Before optimization14.00313.99013.87613.956After optimization3.8323.8733.8713.859Node 17Before optimization14.16914.14914.18414.167After optimization5.7695.7585.7665.764The 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(18)∂σij∂xj+ρbi=ρu¨is,where ρ is the mixture’s density; σij is the total stress tensor; bi is the body force vector; and u¨is is the displacement vector of the soil skeleton.The continuity equation for pore water is derived as(19)nKwdpwdt−nkγw∂2pw∂xi2−ρwε¨iis−ε˙iis=0,where n is the porosity; Kw is the fluid bulk modulus; pw is the pore water pressure; k is the coefficient of permeability; εiis is the volumetric strain of solid; and ρ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:(20)dσij′=Dijklepdεkls,where dεkls is the strain increment; Dijklep is the elastoplastic stiffness matrix of the material; and dσij′ 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(21)M⋅u¨N+C⋅u˙N+K⋅ΔuN+Kv⋅pw=F−Rt,ρwKvT⋅u¨N−γwnkKvT⋅u˙N+A⋅p˙w−∑i=1mαipiw+αpW=0,where M,C, and K are the mass, damping, and stiffness matrix; Rt is the internal force vector; F is the external force vector; Kv is the modulus vector for the pore water pressure; and α 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 timet+Δt are represented according to the variables at time t by the Newmark-β method as follows:(22)u¨t+Δt=u¨t+Δu¨,(23)u˙t+Δt=u˙t+1−γu¨t+γu¨t+Δt×Δt,(24)ut+Δt=ut+u˙t×Δt+121−2βu¨t+2βu¨t+Δt×Δt2,where γ and β are controlling parameters for the Newmark-β method.The pore water pressure at timet+Δt is represented in terms of the variables at time t by the backward difference method as follows:(25)pt+Δtw=ptw+p˙t+Δtw⋅Δt.Upon substituting equations (22)–(25) into governing equation (21) of the u−p form, the final governing equation is obtained as(26)M+ΔtγCt+Δt+Δt2βKt+ΔtΔtKvρw−γwΔtγnkKvTA+αΔtu¨t+Δtp˙t+Δtw+0Δt∑i=1mαip˙i,t+Δtw=Ft+Δt−Rt−Ct+Δtu˙t+Δt1−γu¨t−Kt+ΔtΔtu˙t+Δt212−βu¨t−KvptwγwnkKvTu˙t+Δt1−γu¨t+∑i=1mαipi,tw−αptw.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 Figures5 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 5
2D earth dam.Figure 6
Finite element mesh.The cyclic elastoplastic constitutive model is adopted for soil mass. Soil parameters are shown in Table2. The input horizontal seismic wave is shown in Figure 7.Table 2
Soil properties of 2D earth dam.
Soil parametersSaturated soilDry soilInitial void ratio,e00.800.75Compression index,λ0.0150.018Density,Δt(t/m3)1.9511.85Swelling index0.00450.0055Initial shear modulus ration,G0/Gm0′930.0930.0Permeability coefficient,km/s2.6×10−5Phase transformation stress ratio,Mm0.910.91Failure stress ratio,Mf1.121.12Hardening parameterB025003000Hardening parameterB13060Bulk modulus of pore water,Kw(N/m2)2.0×106Figure 7
Seismic wave in horizontal direction.


### 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Δt0=1.0×10−4s and calculate the first two steps to obtain the initial variable value (ut2,u˙t2,u¨t2, and u¨t1); the rate of acceleration change is determined by the following equation:(27)u…t2=u¨t2−u¨t1Δt0,set the relative error limit η¯=1.0×10−4 and calculate the optimization time step Δt1 by equation (14); then substitute Δt1 into the optimized program to determine ut3,u˙t3,u¨t3, and u⃛t3; then calculate the optimization time step Δt2 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 Figure8: there are a total of 3026 points, the minimum time step is 1.5×10−3, the maximum time step is 0.144s, and the maximum time step is 960 times the minimum time step.Figure 8
Optimization time step for all nodes.The vertical displacement time history curve of the dam crest is shown in Figure9, 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 9
Vertical displacement of earth dam top.Figure 10
Horizontal displacement of earth dam toe.Figure 11
Distribution of excess pore water pressure after time step optimization. (a)t = 2 s. (b) t = 22 s.
(a)(b)Figure 12
Distribution of vertical displacement after time step optimization (t = 22s).


### 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 stepm¯=1.0 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.Figure 13
Distribution of excess pore water pressure before time step optimization (t = 22s).Figure 14
Distribution of vertical displacement before time step optimization (t = 22s).As shown in Figures9 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(28)2716−6752716×100%=75.15%.Due to the drastic change in the external load, the time steps vary in a very large range: the minimum time step is1.5×10−3s, 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.


## 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(18)∂σij∂xj+ρbi=ρu¨is,where ρ is the mixture’s density; σij is the total stress tensor; bi is the body force vector; and u¨is is the displacement vector of the soil skeleton.The continuity equation for pore water is derived as(19)nKwdpwdt−nkγw∂2pw∂xi2−ρwε¨iis−ε˙iis=0,where n is the porosity; Kw is the fluid bulk modulus; pw is the pore water pressure; k is the coefficient of permeability; εiis is the volumetric strain of solid; and ρ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:(20)dσij′=Dijklepdεkls,where dεkls is the strain increment; Dijklep is the elastoplastic stiffness matrix of the material; and dσij′ 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(21)M⋅u¨N+C⋅u˙N+K⋅ΔuN+Kv⋅pw=F−Rt,ρwKvT⋅u¨N−γwnkKvT⋅u˙N+A⋅p˙w−∑i=1mαipiw+αpW=0,where M,C, and K are the mass, damping, and stiffness matrix; Rt is the internal force vector; F is the external force vector; Kv is the modulus vector for the pore water pressure; and α 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 timet+Δt are represented according to the variables at time t by the Newmark-β method as follows:(22)u¨t+Δt=u¨t+Δu¨,(23)u˙t+Δt=u˙t+1−γu¨t+γu¨t+Δt×Δt,(24)ut+Δt=ut+u˙t×Δt+121−2βu¨t+2βu¨t+Δt×Δt2,where γ and β are controlling parameters for the Newmark-β method.The pore water pressure at timet+Δt is represented in terms of the variables at time t by the backward difference method as follows:(25)pt+Δtw=ptw+p˙t+Δtw⋅Δt.Upon substituting equations (22)–(25) into governing equation (21) of the u−p form, the final governing equation is obtained as(26)M+ΔtγCt+Δt+Δt2βKt+ΔtΔtKvρw−γwΔtγnkKvTA+αΔtu¨t+Δtp˙t+Δtw+0Δt∑i=1mαip˙i,t+Δtw=Ft+Δt−Rt−Ct+Δtu˙t+Δt1−γu¨t−Kt+ΔtΔtu˙t+Δt212−βu¨t−KvptwγwnkKvTu˙t+Δt1−γu¨t+∑i=1mαipi,tw−αptw.In equation (26), variable parameters are shown in literature [13].


## 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(18)∂σij∂xj+ρbi=ρu¨is,where ρ is the mixture’s density; σij is the total stress tensor; bi is the body force vector; and u¨is is the displacement vector of the soil skeleton.The continuity equation for pore water is derived as(19)nKwdpwdt−nkγw∂2pw∂xi2−ρwε¨iis−ε˙iis=0,where n is the porosity; Kw is the fluid bulk modulus; pw is the pore water pressure; k is the coefficient of permeability; εiis is the volumetric strain of solid; and ρ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:(20)dσij′=Dijklepdεkls,where dεkls is the strain increment; Dijklep is the elastoplastic stiffness matrix of the material; and dσij′ 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(21)M⋅u¨N+C⋅u˙N+K⋅ΔuN+Kv⋅pw=F−Rt,ρwKvT⋅u¨N−γwnkKvT⋅u˙N+A⋅p˙w−∑i=1mαipiw+αpW=0,where M,C, and K are the mass, damping, and stiffness matrix; Rt is the internal force vector; F is the external force vector; Kv is the modulus vector for the pore water pressure; and α 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 timet+Δt are represented according to the variables at time t by the Newmark-β method as follows:(22)u¨t+Δt=u¨t+Δu¨,(23)u˙t+Δt=u˙t+1−γu¨t+γu¨t+Δt×Δt,(24)ut+Δt=ut+u˙t×Δt+121−2βu¨t+2βu¨t+Δt×Δt2,where γ and β are controlling parameters for the Newmark-β method.The pore water pressure at timet+Δt is represented in terms of the variables at time t by the backward difference method as follows:(25)pt+Δtw=ptw+p˙t+Δtw⋅Δt.Upon substituting equations (22)–(25) into governing equation (21) of the u−p form, the final governing equation is obtained as(26)M+ΔtγCt+Δt+Δt2βKt+ΔtΔtKvρw−γwΔtγnkKvTA+αΔtu¨t+Δtp˙t+Δtw+0Δt∑i=1mαip˙i,t+Δtw=Ft+Δt−Rt−Ct+Δtu˙t+Δt1−γu¨t−Kt+ΔtΔtu˙t+Δt212−βu¨t−KvptwγwnkKvTu˙t+Δt1−γu¨t+∑i=1mαipi,tw−αptw.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 Figures5 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 5
2D earth dam.Figure 6
Finite element mesh.The cyclic elastoplastic constitutive model is adopted for soil mass. Soil parameters are shown in Table2. The input horizontal seismic wave is shown in Figure 7.Table 2
Soil properties of 2D earth dam.
Soil parametersSaturated soilDry soilInitial void ratio,e00.800.75Compression index,λ0.0150.018Density,Δt(t/m3)1.9511.85Swelling index0.00450.0055Initial shear modulus ration,G0/Gm0′930.0930.0Permeability coefficient,km/s2.6×10−5Phase transformation stress ratio,Mm0.910.91Failure stress ratio,Mf1.121.12Hardening parameterB025003000Hardening parameterB13060Bulk modulus of pore water,Kw(N/m2)2.0×106Figure 7
Seismic wave in horizontal direction.


## 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Δt0=1.0×10−4s and calculate the first two steps to obtain the initial variable value (ut2,u˙t2,u¨t2, and u¨t1); the rate of acceleration change is determined by the following equation:(27)u…t2=u¨t2−u¨t1Δt0,set the relative error limit η¯=1.0×10−4 and calculate the optimization time step Δt1 by equation (14); then substitute Δt1 into the optimized program to determine ut3,u˙t3,u¨t3, and u⃛t3; then calculate the optimization time step Δt2 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 Figure8: there are a total of 3026 points, the minimum time step is 1.5×10−3, the maximum time step is 0.144s, and the maximum time step is 960 times the minimum time step.Figure 8
Optimization time step for all nodes.The vertical displacement time history curve of the dam crest is shown in Figure9, 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 9
Vertical displacement of earth dam top.Figure 10
Horizontal displacement of earth dam toe.Figure 11
Distribution of excess pore water pressure after time step optimization. (a)t = 2 s. (b) t = 22 s.
(a)(b)Figure 12
Distribution of vertical displacement after time step optimization (t = 22s).


## 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 stepm¯=1.0 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.Figure 13
Distribution of excess pore water pressure before time step optimization (t = 22s).Figure 14
Distribution of vertical displacement before time step optimization (t = 22s).As shown in Figures9 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(28)2716−6752716×100%=75.15%.Due to the drastic change in the external load, the time steps vary in a very large range: the minimum time step is1.5×10−3s, 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 Figures7 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.



---

*Source: 8863211-2020-09-29.xml*