砂土地震液化分析中Newmark时域离散的误差评估

张西文，唐小微，渦岡良介，胡记磊，张帅芳，付培帅

（1.大连理工大学岩土工程研究所，辽宁大连116024；2.德岛大学岩土工程研究所，德岛市日本770⁃8506；3.大连理工大学工程力学研究所，辽宁大连116024）

砂土地震液化分析中Newmark时域离散的误差评估

张西文1，2，唐小微1，渦岡良介2，胡记磊1，张帅芳3，付培帅1

（1.大连理工大学岩土工程研究所，辽宁大连116024；2.德岛大学岩土工程研究所，德岛市日本770⁃8506；3.大连理工大学工程力学研究所，辽宁大连116024）

显式有限元存在计算精度低，对计算时间步长较敏感等缺点。基于后验误差评估的方法，给出了显式算法下Newmark时域离散误差的来源和评估方法。通过饱和砂土地震液化响应的数值算例评估了时间步长对计算结果的影响。数值分析结果表明：时间步长不同，结点位移和单元孔压的时程曲线明显不同，同时计算耗时也呈双曲线关系；相对误差主要分布在变形较大的区域，全域平均相对误差在动力响应剧烈的时间段内较大。通过对计算时间步长和离散误差的评估，可有效恰当的对计算时间步长进行取值，也为自动步长调整提供了依据。

误差评估；砂土液化；堤坝；地震；Newmark法；计算成本

对微分方程的求解中常用到隐式算法和显式算法。隐式方法通过迭代计算，降低了时间步长大小对结果的影响，但计算过程非常耗时，在静力计算中应用较多。显式方法则直接根据当前时间步的状态推导出下一时间步的状态，具有快捷性和高效性而被广泛应用于动力问题的求解中。但是显式算法的计算结果对时间步长的敏感性较高，如何权衡计算精度和计算效率成为使用者较为头疼的问题。因此，对时域离散误差进行正确评估和进行时间步长自动调整是显式计算中的关键技术。在时域积分中常用的方法有中心差分法、Newmark法［1］、Houbolt法［2］和Wilson⁃θ法［3］等。其中Newmark⁃β时域离散的方法最为实用，在结构动力学和砂土地震液化的计算中被广泛使用。关于Newmark法很多学者曾进行了改进来提高计算精度［4⁃6］，本文的工作是对显式二相耦合计算中的Newmark⁃β离散方法进行了误差评估，并通过砂土液化的算例进行了验证。关于时域离散误差的评估前人做了很多工作，如Zienkiewicz等［7］在单自由度的动力分析中，对比Taylor级数展开式获得了时域离散的截断误差。Zeng等［8⁃9］改进了Zienkiewicz的方法，提出了一种后验的评估方法。在岩土工程领域，Sloan等［10⁃11］探讨了时间步长对固结计算的影响。

土体在地震作用下的液化流动特性及堤坝的地震破坏目前被广泛研究［12⁃17］，常用方法也是Ne⁃wmark显式计算方法。本文首先是基于有限元-有限差分（FEM⁃FDM）的分析方法［18⁃19］和循环弹塑性的本构模型［20⁃21］建立了u⁃p格式的控制方程。然后采用Zeng等［8］提出的后验误差评估方法u⁃p格式显式计算进行了误差评估。同时也采用堤坝作为数值算例，对误差评估结果进行验证。

1 控制方程及Newmark时域离散方法

对饱和砂土，将加速度的影响作为惯性力，忽略孔隙流体的加速度，可得水-土二相混合体的平衡方程：

式中：σij为应力张量，u··si为土颗粒的加速度，bi为体力向量。

对孔隙水，建立连续性方程：

式中：n为孔隙率，Kw为水的体积模量，pw为孔隙水压力，k为渗透系数，为土体应变张量对时间的导数。

在空间域进行有限元离散，考虑超孔隙水压力和阻尼的影响，可建立在动力荷载作用下u⁃p格式的控制方程：

为了简便，采用了Rayleigh阻尼，阻尼C取如下形式：

在时间域采用Newmark⁃β法对速度和位移进行离散，即获得t+Δt时刻与t时刻变量的关系：

将式（5）、（6）表示的t+Δt时刻位移、速度代入式（3）中，得控制方程的最终表达形式为

2 后验时域离散误差评估

根据Zeng等［8］的方法，假定t时刻的变量为已知量，t+Δt时刻的变量为未知量。且存在τ，τ∈［t，t+Δt］。由式（6），Newmark法在［t，t+Δt］内计算的平均加速度为

克深区块单井试油期间开展测试放空气回收，在钻井显示发现良好后，开始组织测试气回收，由承包商铺设回收管线，并在井场连接好回收气设备，试油放喷开始后直接开展放空气回收，大大减少天然气放空，降低耗电量，实现节能减排，保护环境，取得了良好的效果。

在［t，t+Δt］内，更精确的加速度可描述为

在τ时刻的加速度误差可表示为u··NM与u··∗的差值：

在t+Δt时刻的位移误差可用加速度误差的两次积分表示：

式（11）给出的为绝对误差值，为了表示误差对计算结果的影响，本文后面提到的误差都是以相对误差的形式表示：

3 堤坝的地震响应分析

3.1 有限元模型

本算例为一个堤坝在地震荷载作用下的动力响应，计算模型如图1所示。底面和左右两侧为固定位移边界条件，顶面为排水边界。结点N1、N2和单元E1作为计算的输出结点和单元。

图1 堤坝计算模型Fig.1 The model of embankment

土体分为干土和饱和砂土2个部分，采用循环弹塑性本构模型，土体的数值参数如表1。

输入的地震荷载时程如图2所示，输入地震波的步长为0.01 s，地震波时长为22 s。

表1 土体参数Table1 Soil properties in numerical example

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

3.2 不同时间步长时堤坝的动力响应

图5表示了单元E1在不同计算时间步长时的超孔隙水压力比（EPWPR）的时程曲线。从图中可以看出，在振动过程中，超孔隙水压力比从零变为1，即土体由初始状态逐渐变为液化状态。当取不同时间步长时，超孔隙水压力的发展是不一致的，初始液化的时间也不相同，也就是说时间步长的大小对孔隙水压力的变化也有影响。

图3 不同时间步长时N1点竖向位移时程曲线Fig.3 Vertical displacement history of point N1in dif⁃ferent time step

图4 不同时间步长时N2点水平位移时程曲线Fig.4 Horizontal displacement history of point N1in different time step

图5 不同时间步长时单元E1的超孔隙水压力比时程曲线Fig.5 EPWPR history of E1in different time step

图6 表示了地震结束之后，超孔隙水压力比EPW⁃PR的云图，可见在水平地震作用下，堤坝周围的饱和砂土层已经基本达到了液化状态。图7表示了堤坝最终的变形图，与初始网格对比可见周围土体的液化使堤坝产生了较大的竖向沉陷和水平扩展。

图6 超孔隙水压力比云图Fig.6 Distribution of EPWPR

图7 网格变形图（Δt＝1×10-5s）Fig.7 Final deformation of the mesh（Δt＝1×10-5s））

为了对比不同时间步长时的计算成本，本文统一采用IBM Intel（R）Xeon（R）3.00 GHz处理器进行了计算，图8表示了CPU计算时间和时间步长大小的关系曲线。可以看出随着时间步长的减小，计算成本成双曲线式增长。通过拟合，得计算时间y与时间步长x的关系式如下

可见对于水土二相耦合显式计算，如果采用较大步长，则会产生较大的计算误差，采用较小的步长则会耗费计算时间。

图8 不同时间步长时CPU计算时间Fig.8 Computation time in different time steps

3.3 时域离散误差评估

为了分析上述不同步长下的计算差异，根据式（11）、（12），计算了节点位移的相对误差，并绘制了分布云图，如图9所示。可以看出，误差较大的区域分布在变形较大的堤坝下部和两侧。另外，随着时间步长的减小，时域离散的误差也随之减小，且呈现出不同的数量级。

为了进一步对整个计算区域进行误差评估，下式表示了计算域的平均相对误差：

式中：n表示结点个数；j＝1，2表示x和y方向的相对误差。

图10表示了2个不同时间步长时，平均相对误差的时程曲线。可以发现，在前10 s地震响应较剧烈的时间段内，全域平均相对误差较大。在不同的时间步长下误差时程曲线的形状类似，误差的大小却有数量级上的差别。

图10 全域平均相对误差时程曲线Fig.10 Time history of average relative error

4 结论

本文采用有限元后验误差计算的方法对显式Newmark⁃β时域离散误差进行了评估，并应用于堤坝砂土液化的动力数值分析中。从分析结果中可得以下结论：

1）显式有限元的计算精度与时间步长的大小密切相关，时间步长越小越精确，但花费的计算成本随时间步长的减小成双曲线增长。

2）通过后验误差评估的方法，给出了Newmark⁃β时域离散误差的来源和评估方法，在砂土液化的数值算例中给出了相对误差的分布和时程曲线。可以发现时间步长越大，相对误差则越大；变形较大的区域和地震响应剧烈的时间段内相对误差较大。

3）通过时域离散相对误差的评估可以帮助我们更客观地选择时间步长的大小，同时为自动时间步长的调整提供了依据。

［1］NEWMARK N M.A method of computation for structure dy⁃namics［J］.Journal of the Engineering Mechanics Division，1959，85（3）：67⁃94.

［2］HOUBOLT J C.A recurrence matrix solution for the dynamic response of elastic aircraft［J］.Journal of the Aeronautical Sciences，1950，17（9）：540⁃550.

［3］WILSON E L.Nonlinear dynamic analysis of complex struc⁃tures［J］.Earthquake Engineering and Structural Dynamics，1973，1：241⁃252.

［4］钟万勰.结构动力方程的精细时程积分法［J］.大连理工大学学报，1994，34（2）：131⁃136.ZHONG Wanxie.On precise time⁃integration method for structural dynamics［J］.Journal of Dalian University of Technology，1994，34（2）：131⁃136.

［5］郭泽英，李青宁，张守军.结构地震反应分析的一种新精细积分法［J］.工程力学，2007，24（4）：35⁃40.GUO Zeying，LI Qingning，ZHANG Shoujun.A new preciseintegration method for structural seismic response analysis［J］.Engineering Mechanics，2007，24（4）：35⁃40.

［6］李爽，翟长海，刘洪波，等.Newmark积分方法在负刚度时的数值稳定性分析［J］.哈尔滨工业大学学报，2011，43（8）：1⁃5.LI Shuang，ZHAI Changhai，LIU Hongbo et al.On the sta⁃bility of Newmark integration method for structure with nega⁃tive stiffness［J］.Journal of Harbin Institute of Technology，2011，43（8）：1⁃5.

［7］ZIENKIEWICZ O C，XIE Y M.A simple error estimator and adaptive time stepping procedure for dynamic analysis［J］.Earthquake Engineering＆Structural Dynamics，1991，20（9）：871⁃887.

［8］ZENG L F，WIBERG N E，LI X D.A posteriori local error estimation and adaptive time⁃stepping for Newmark integra⁃tion in dynamic analysis［J］.Earthquake Engineering＆Structural Dynamics，1992，21（7）：555⁃571.

［9］WIBERG N E，LI X D.A post⁃processing technique and an a posteriori error estimate for the newmark method in dynam⁃ic analysis［J］.Earthquake Engineering＆Structural Dy⁃namics，1993，22（6）：465⁃489.

［10］SLOAN S W，ABBO A J.Biot consolidation analysis with automatic time stepping and error control Part 2：Applica⁃tions［J］.International Journal for Numerical and Analyti⁃cal Methods in Geomechanics，1999，23（6）：493⁃529.

［11］阮光雄，傅少君，陈胜宏.固结问题有限元法时步自适应研究［J］.岩土力学，2005，26（4）：591⁃599.NGUYEN Quanghung，FU Shaojun，CHEN Shenghong.Study on adaptive time step of consolidation problems by fi⁃nite element method［J］.Rock and Soil Mechanics，2005，26（4）：591⁃599.

［12］IAI S.Three dimensional formulation and objectivity of a strain place multiple mechanism model for sand［J］.Soils and Foundations，1993，33（1）：192⁃199.

［13］汪明武.可液化场地土工抗震性能动态离心机试验与数值模拟［M］.北京：科学出版社，2012：110⁃130.WANG Mingwu.Dynamic centrifuge tests and numerical modelling for geotechnical earthquake engineering in lique⁃fiable soils［M］.Beijing：Science Press，2012：110⁃130.

［14］TANG Xiaowei，SHAO Qi.Numerical simulation on seismic liquefaction by adaptive mesh ref i nement due to two recov⁃ered f i elds in error estimation［J］.Soil Dynamics and Earthquake Engineering，2013，49：109⁃121.

［15］WANG Mingwu，CHEN Guangyi，IAI S.Seismic perform⁃ances of dyke on liquefiable soils［J］.Journal of Rock Me⁃chanics and Geotechnical Engineering，2013，5（4）：294⁃305.

［16］ZHANG Xiwen，TANG Xiaowei，SHAO Qi，et al.The up⁃lift behavior of large underground structures in liquefied field［J］.Applied Mechanics and Materials，2011，90⁃93：2112⁃2118.

［17］兰景岩，刘红帅，吕悦军.渤海土类动力非线性参数及合理性［J］.哈尔滨工程大学学报，2012，33（9）：1079⁃1085.LAN Jingyan，LIU Hongshuai，LYU Yuejun.Dynamic non⁃linear parameters of soil in the Bohai Sea and their rational⁃ity［J］.Journal of Harbin Engineering University，2012，33（9）：1079⁃1085.

［18］AKAI K，TAMURA T.Numerical analysis of multi⁃dimen⁃sional consolidation accompanied with elastic⁃plastic con⁃stitutive equation［C］／／Proceedings of the Society of Civil Engineering，1978，269（3）：95⁃104.

［19］OKA F，YASHIMA A，SHIBATA T，et al.FEM⁃FDM cou⁃pled liquefaction analysis of a porous soil using an elastic⁃plastic model［J］.Applied Scientific Reasearch 1994，52：209⁃245.

［20］OKA F，YASHIMA A，TATEISHI A，et al.A cyclic elas⁃to⁃plastic constitutive model for sand considering a plastic⁃strain dependence of the shear modulus［J］.Geotech⁃nique，1999，49（5）：661⁃680.

［21］UZUOKA R.Analytical study on the mechanical behavior and prediction of soil liquefaction and f l ow［D］.Kyoto：Kyoto University，2000：47⁃77.

Temporal discretization error estimation for the Newmark scheme in sand liquefaction analysis

ZHANG Xiwen1，2，TANG Xiaowei1，UZUOKA Ryosuke2，HU Jilei1，ZHANG Shuaifang3，FU Peishuai1
（1.Institute of Geotechnical Engineering，Dalian University of Technology，Dalian 116024，China；2.Institute of Geotechnical Engi⁃neering，The University of Tokushima，Tokushima 770⁃8506，Japan；3.Institute of Engineering Mechanics，Dalian University of Tech⁃nology，Dalian 116024，China）

The disadvantage of the explicit finite element method is low accuracy and sensitive to the time step size.In this paper，a posteriori error evaluation method was introduced for the explicit Newmark scheme，giving the source and estimation method of the Newmark temporal discretization error.Then，a numerical example of saturated sand liquefaction in earthquake was conducted to evaluate the influence of time step size on the calculation result.The results showed that the time step size has effects on the time histories of node displacements and element pore 's water pressures.Moreover，the relation between the computation cost and the time step size is a hyperbolic curve.The relative error is mainly generated in the area with large deformation and the period with rapid dynamic response.The mean relative error in the entire area is larger in the time range with violent dynamic response.Through the es⁃timated discrete error，a proper time step size can be determined to meet the desired accuracy.For automatic time step adjustment，the error estimation can provide a criterion to control the time step size.

error estimation；sand liquefaction；embankment；earthquake；Newmark method；computation time

10.3969／j.issn.1006⁃7043.201311001

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1657.020.html

TU43

A

1006⁃7043（2015）03⁃0322⁃05

2013⁃11⁃04.网络出版时间：2015⁃01⁃09.

国家863计划资助项目（2012AA112510）.

张西文（1987⁃），男，博士研究生；

唐小微（1968⁃），男，副教授，博士.

唐小微，E⁃mail：tangxw＠dlut.edu.cn.