自适应时间步长法及三维堤坝地震液化数值模拟①

自适应时间步长法及三维堤坝地震液化数值模拟①

通信作者：唐小微(1968-)，男，博士生导师，主要从事土动力学、海洋土动力学等方面研究。E-mail：tangxw@dlut.edu.cn。

张西文1， 唐小微2， 白旭2， 唐晓成3， 韩小凯2

(1.济南大学土木建筑学院，山东 济南 250022； 2.大连理工大学海岸与近海工程国家重点实验室，辽宁 大连 116024；

3.东北电力大学建筑工程学院，吉林 吉林 132000)

摘要：三维大模型数值计算因巨大的单元和结点数目而非常耗时，在地震响应分析中受计算时间步长的限值则更加耗时。在饱和砂土动力液化计算平台上开发时域离散误差评估方法和时间步长自适应调整的计算程序，并成功应用于三维堤坝地震液化响应分析。时域离散误差包括土骨架的位移误差和单元孔压误差，通过定义孔压误差影响系数计算出混合误差，根据混合误差和设定的误差允许值进行计算步长的自适应调整。在三维堤坝地震液化数值模拟中，采用自适应时间步长法有效避免小步长精确但耗时、大步长省时而不精确的缺点。在大模型和超大模型计算中，最优调整每一步的计算时间步长，完美实现既节省时间又不失精度的时域离散策略。

关键词：自适应时间步长法； 混合误差； 堤坝； 地震液化； 计算时间

收稿日期：①2014-08-20

基金项目：国家自然科学

作者简介：：张西文(1987-)，男，讲师，主要从事地震液化数值模拟、自适应等方法研究。E-mail：zhangxw_11@163.com。

中图分类号:TU43 文献标志码：A

DOI:10.3969/j.issn.1000-0844.2015.03.0672

Three-Dimensional Numerical Simulation of Embankment Seismic

Liquefaction Using Adaptive Time-Stepping Method

ZHANG Xi-wen1, TANG Xiao-wei2, BAI Xu2, TANG Xiao-cheng3, HAN Xiao-kai2

(1.SchoolofCivilEngineeingandArchitecture,TheUniversityofJinan,Jinan250022,Shangdong,China;

2.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,Liaoning，China;

3.SchoolofCivilEngineeringandArchitecture,NortheastDianliUniversity,Jilin132000,Jilin,China)

Abstract：Three-dimensional (3D) numerical simulation of a large model is a time-consuming process due to the huge number of elements and nodes. This is particularly true for seismic behavior analysis in which the accuracy is limited by the size of the time increment. In this study, the temporal discretization error estimator and time increment adaptive adjustment method is proposed on the basis of the numerical platform for seismic liquefaction analysis of saturated soil. The proposed method is successfully applied in 3D embankment seismic response analysis. The temporal discretization error estimation includes errors in soil skeleton displacement and pore water pressure. These two types of errors are combined to form a mixed error as an effect of pore water pressure error. The time increment is adaptively adjusted by the relationship of mixed error and the given error tolerance. Through a numerical example of 3D embankment, it is determined that the adaptive time-stepping method can be used to improve the accuracy of small time increments. The results of this method are comparable to those obtained using highly efficient but time-consuming large time increments.

Key words： adaptive-time stepping method; mixing error; embankment; seismic liquefaction; computation time

0引言

饱和砂土地基的液化是常见地震灾害之一，在中国的汶川地震(2008)、日本阪神地震(1995)和3·11大地震(2011)中均发现了液化现象和破坏的实例[1-2]。对于堤坝尤其是处于可液化地基上的堤坝来讲，由地震液化造成的竖向沉陷和横向侧移是最重要的破坏模式。在过去的几十年里众多学者对堤坝进行了研究[3-7]，其研究方法可以归结为经验理论法、振动台试验法和数值模拟法。随着计算机科学技术的发展，数值计算方法逐渐受到众多科研工作者的青睐。为了达到数值计算的精确度，开始采用更细密的单元划分，更小的计算时间步长，但同时计算效率和计算精度的矛盾也显得更为突出。自适应有限元的出现开创了自动调整网格尺寸和自动调整计算时间步长的策略，自动进行最佳的离散方式，有效实现了效率和精度的统一。网格自适应方面，Babuska[8]最早提出了残差法为依据的网格调整方法，由于形式复杂没有得到广泛应用。随后被广泛应用的则是由Zienkiewicz等[9]提出的后验误差评估方法(也称为z-z误差评估)。在时间离散方面，起步和发展都较网格自适应晚，Bergan[10]在1985年提出了一种时间步调整的方法，Zienkiewicz[11]将Newmark解与Taylor级数解的差值作为截断误差。Zeng[12]对Zienkiewicz的方法进行了改进，提出了后验的误差计算方法，误差评估的对象是单元的结点位移，在不考虑孔隙水的情况下是适用的。此后，Zhang Zihua、王怀忠和冯领香等都采用了Zeng的误差评估方法[13-15]。

时间自适应算法在结构和流体工程中应用较广，如Hulbert[16]在结构分析中采用了自适应时间步长法，Tang[17]、Hirthe[18]在流体的计算中都采用了变步长法和自动步长法。在岩土工程领域，Sloan[19-20]、阮光雄[21]在岩土固结分析中采用自适应时间步长法，极大节省了计算时间。在岩土动力分析中，自适应时间步长法的应用和可查阅的文献非常少，在水土耦合分析和砂土动力液化中研究更少。然而动力计算常采用显式有限元的方法，计算的精度和时间步长的选取存在必然的联系；同时，由于计算时间长，计算步数多，模型巨大而耗费大量的时间。因此对于多单元大模型的动力计算问题，开发一套误差评估方法和时间步长自适应调整程序是非常有必要的。

在本文的研究中，首先采用有限元-有限差分(FEM-FDM)的方法进行水土耦合分析。除了结点的位移误差需要考虑之外，单元孔隙水压力的误差也是不可忽略的。因此，借鉴自适应时间步长法的研究进展，对结点位移误差和单元孔压误差同时进行评估，并将二者的混合误差作为时间步长调整的依据，开发适合于饱和砂土液化分析的自适应步长法程序。

1水土二相耦合计算理论

饱和砂土可看做是由土和水组成的二相混合物。假定土颗粒本身和水本身的密度为ρs和ρw，混合体的孔隙率为n，则混合体的密度可以表示为：

采用有效应力原理建立有效应力和总应力之间的关系为：

在动力荷载作用下，将加速度对混合体的影响作为惯性力，并忽略孔隙流体的加速度，可得水-土二相混合体的平衡方程。

式中:usi为土骨架的位移向量；ρ为混合体的密度；bi为体力向量。

根据达西定律和质量守恒条件，可建立孔隙水连续性方程：

式中:n为土体的孔隙率；pw为土体孔隙水压力；k为土体的渗透系数；γw为孔隙水的重度；εsii为土体骨架的应变张量。

通过有限单元方法进行空间离散，对和孔压有关的相用有限差分法离散，获得u-p格式的控制方程：

式中:[C]为阻尼阵，采用Rayleigh阻尼，表示为：

采用循环弹塑性本构建立应力增量和应变增量的关系：

式中:Depijkl为弹塑性矩阵。关于此本构的具体方程和本构参数设置可参考Oka等[22]的文献。

对结点位移采用Newnark-β法进行时间域离散。

对单元孔压，采用向后差分的方法进行离散，

将式(8)和式(9)表示的t+Δt时刻的速度和位移、式(10)表示的t+Δt时刻的孔压代入式(5)中，经过整理后得到最终的控制方程式为

2时域离散误差评估及自适应调整

2.1位移误差评估

图1　Newmark-β法位移时间离散误差 Fig.1　Temporal discretization error for node displacement 　　　 using Newmark-β method

由式(9)可知当采用Newmark-β法时，在[t,t+Δt]内的加速度近似值为：

按照Zeng的线性化假设，存在一个变量τ∈[t,t+Δt]，更精确的加速度描述可表示为：

在τ时刻的速度误差可表示为加速度误差在[t,τ]区间上的积分，即：

在t+Δt时刻位移误差可表示为速度误差在[t,t+Δt]区间上的积分：

位移相对误差表示为：

2.2孔压误差评估

在差分法中，向前和向后差分为1阶精度，中心差分则为2阶精度。在孔压误差的评估中，将中心差分的解和向后差分的解之间的差值作为孔压的截断误差。向后差分的解见式(18)，中心差分的解见式(19)。

孔压的误差表示为：

孔压的相对误差表示为：

2.3时间步长自适应方法

为了综合考虑位移误差和孔压误差，采用混合误差作为时间步长调整的依据。混合误差的定义为：

经过自适应调整后的时间步长为：

3三维堤坝数值算例

3.1三维堤坝模型

图2　堤坝数值算例(单位：m) Fig.2　Numerical model of the embankment (Unit:m)

堤坝土体分为干土和饱和砂土两部分，土体的参数见表1。

3.2输入地震波荷载

输入的地震波见图 3，为水平方向振动。地震波加速度峰值为0.3g。地震波时间步长为0.01s，共2 200步，持续振动22s。

3.3位移分析

(1) 坝顶竖向位移时程

图 4给出了坝顶N1点的竖向位移时程曲线。可以发现，位移曲线与时间步长的大小密切相关，时间步长越小则计算出的竖向位移越大。采用自适应时间步长法时，可以获得与固定步长法一致的变化规律，且当误差允许值越小时越接近小步长的结果，即越精确。

表 1　土体参数

图3　输入地震波 Fig.3　Input seismic wave

图4　坝顶竖向位移时程曲线 Fig.4　Time history of vertical displacement on the 　　　 embankment crest

(2) 坝趾水平位移时程

图 5 给出了坝趾N2点的水平位移时程曲线。采用自适应时间步长法，通过设定误差允许值可以获得与小步长同样的精度。同时，也可以发现在地震荷载作用下堤坝会发生竖向和横向的较大位移。

图5　坝趾水平位移时程曲线 Fig.5　Time history of horizontal displacement 　　　 at the embankment toe

(3) 位移云图

图 6 给出了地震结束后的竖向位移分布图。可以发现整个堤坝发生了较大的沉陷变形。

图6　位移云图(t=22 s) Fig.6　Distribution of the vertical displacement (t=22 s)

3.4超孔隙水压力比分布

图 7(a)和(b)给出了地震发生过程中t=2 s和地震结束时t=22 s的超孔隙水压力比(epwpr)的云图。可以发现液化范围在地震过程中不断扩大，也可以推断由于发生了地基液化而造成了堤坝的较大变形。

图7　超孔隙水压力比云图 Fig.7　Distribution of the excess pore water pressure

3.5计算时间步长

图 8表示采用自适应步长法时的计算时间步长时程曲线。通过自适应步长的算法实现了计算时间步长实时自动的调整，可以合理地分配计算资源。当采用较小的混合误差允许值时，从图 4和图 5可见计算更加精确，同时从图 8中可以发现时间步长的取值会相应变小。

图8　计算时间步长的时程曲线 Fig.8　Time history of computation time step

3.6计算时间对比

为了对计算效率作对比，本文所有算例统一在DELL Precision T5610 Intel®Xeon®3.4 GHz 的工作站上计算。计算时间列于表2中。可见在相同计算精度条件下，自适应时间步长法更节省计算时间。

表 2　计算时间对比

4结论

采用自适应时间步长动力有限元法对三维堤坝进行地震液化灾害模拟。通过对结果的分析，可以得到以下结论:

(1) 计算时间步长的取值大小决定计算的精度，但获得高精度的代价则是花费巨大的计算资源。自适应步长法可以实时地调整计算时间步长大小，合理地分配计算资源，与固定步长法相比，在获得相同精度的条件下可最大地节省计算时间。

(2) 通过数值模拟，研究三维堤坝在地震液化情况下的破坏模式，包括竖向的沉陷变形和水平的流动扩展变形，为揭示堤坝地震液化灾害机理提供依据和参考。

参考文献(References)

[1]黄雨,于淼,Bhattacharya S.2011年日本东北地区太平洋近海地震地基液化灾害综述[J].岩土工程学报:2013，35(5)：834-840.

HUANG Yu,YU Miao,Bhattacharya S.Review on Liquefaction Induced Damages of Soils and Foundations During the 2011 off the Pacific Coast of Tohoku Earthquake (Japan)[J].Chinese Journal of Geotechnical Engineering,2013，35(5)：834-840.(in Chinese)

[2]Chen L W,Yuan X M,Cao Z Z,et al.Liquefaction Macrophenomena in the Great Wenchuan Earthquake[J].Earthquake Engineering and Engineering Vibration, 2009,8 (2):219-229.

[3]Ishihara K,Yasuda S,Yoshida Y.Liquefaction Induced Flow Failure of Embankment and Residual Strength of Silty Sand[J].Soils and Foundations, 1990, 30 (3): 69-80.

[4]Manzari M T.Finite Deformation Analysis of Liquefaction Induced Flow Failure in Soil Embankments[J].Geotechnical and Geological Engineering, 1996, 14: 83-110.

[5]Elgamal A,Parra E,Yang Z,et al.Numerical Analysis of Embankment Foundation Liquefaction Countermeasures[J].Journal of Earthquake Engineering, 2002, 6(4): 447-471.

[6]汪明武,Iai S,Tobita T.液化场地堤坝地震响应动态土工离心试验及模拟[J].水利学报,2008,39: 1346-1352.

WANG Ming-wu,Iai S,Tobita T.Centrifuge Test and Numerical Analysis of Seismic Responses of Dyke on Liquefiable Soils Foundation[J].Journal of Hydraulic Engineering,2008,39:1346-1352. (in Chinese)

[7]HuangY,Yashima A,Sawada K,et al.A Case Study of Seismic Response of Earth Embankment Foundation on Liquefiable Soils[J].Journal of Central South University of Technology,2009,16(6):994-1000.

[8]Babuska I,Rheinboloat W C.Error Estimates for Adaptive FiniteElement Computations[J].Siam Journal on Numerical Analysis,1978,15(4):736-754.

[9]Zienkiewicz O C,Zhu J Z.A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis[J].International Journal for Numerical Methods in Engineering,1987,24:337-357.

[10]Bergan P G, Mollestad E.An Automatic Time Stepping Algorithm for Dynamic Problems[J].Computer Methods in Applied Mechanics and Engineering, 1985, 49(3):299-318.

[11]Zienkiewicz O C,Xie Y M.A Simple Error Estimator and Adaptive Time Stepping Procedure for Dynamic Analysis[J].Earthquake Engineering and Structure Dynamics, 1991,20:871-887.

[12]Zeng L F,Wiberg N E,Li X D.A Posteriori Local Error Estimation and Adaptive Time-Stepping for Newmark Integration in Dynamic Analysis[J].Earthquake Engineering and Structure Dynamics,1992,21:555-571.

[13]Zhang Z H,Yang Z J,Liu G H.An Adaptive Time-Stepping Procedure Based on the Scaled Boundary Finite Element Method for Elastodynamics[J].International Journal of Computational Methods,2012,9(1):13.

[14]王怀忠,孙钧.能量的局部误差估计和直接积分法的时步自适应算法[J].岩土工程学报,1995,17(5):32-41.

WANG Huai-zhong,SUN Jun.Local Error Estimator of Energy and Adaptive Time Stepping Procedure for Direct Integration Method[J].Chinese Journal of Geotechnical Engineering, 1995,17(5):32-41.(in Chinese)

[15]冯领香,魏建国,王森林,等.一种可自调步长的改进newmark算法[J].河北农业大学学报,2004,27(3):111-114.

FENG Ling-xiang,WEI Jian-guo,WANG Seng-lin.The Improved Newmark Algorithm With Adaptive Step Size[J].Journal of Agricultural University of Hebei,2004,27(3):111-114.(in Chinese)

[16]Hulbert G M,Jang I.Automatic Time Step Control Algorithms for Structural Dynamics[J].Computer Methods in Applied Mechanics and Engineering,1995,126:155-178.

[17]Tang G P,Alshawabkeh A N,Mayes M A.Automatic Time Stepping with Global Error Control for Groundwater Flow Models[J].Journal of Hydrologic Engineering,2008,13(9):803-810.

[18]Hirthe E M,Graf T.Non-Iterative Adaptive Time Stepping Scheme with Temporal Truncation Error Control for Simulating Variable-density Flow[J].Advances in Water Resources,2012,49:46-55.

[19]Sloan S W,Abbo A J.Biot Consolidation Analysis with Automatic Time Stepping and Error Control,Part 1:Theory and Implementation[J].International Journal for Numerical and Analytical Methods in Geomechanics,1999,23(6):467-492.

[20]Sloan S W,Abbo A J.Biot Consolidation Analysis with Automatic Time Stepping and Error Control,Part 2: Applications[J].International Journal for Numerical and Analytical Methods in Geomechanics,1999,23(6):493-529.

[21]阮光雄,傅少君,陈胜宏.固结问题有限元法时步自适应研究[J].岩土力学,2005, 26(4):591-599.

Ruan Guang-xiong,FU Shao-jun,CHEN Sheng-hong.Study on Adaptive Time Step of Consolidation Problems by Finite Element Method[J].Rock and Soil Mechanics,2005,26(4):591-599.(in Chinese)

[22]Oka F,Yashima A,Tateishi A,et al.A Cyclic Elastic-plastic Constitutive Model for Sand Considering a Plastic-strain Dependence of the Shear Modulus[J].Geotechnique,1999,49(5):661-680.(in Chinese)