砂土地震液化分析中Newmark時(shí)域離散的誤差評(píng)估

張西文，唐小微，渦岡良介，胡記磊，張帥芳，付培帥

（1.大連理工大學(xué)巖土工程研究所，遼寧大連116024；2.德島大學(xué)巖土工程研究所，德島市日本770?8506；3.大連理工大學(xué)工程力學(xué)研究所，遼寧大連116024）

砂土地震液化分析中Newmark時(shí)域離散的誤差評(píng)估

張西文1，2，唐小微1，渦岡良介2，胡記磊1，張帥芳3，付培帥1

（1.大連理工大學(xué)巖土工程研究所，遼寧大連116024；2.德島大學(xué)巖土工程研究所，德島市日本770?8506；3.大連理工大學(xué)工程力學(xué)研究所，遼寧大連116024）

顯式有限元存在計(jì)算精度低，對(duì)計(jì)算時(shí)間步長(zhǎng)較敏感等缺點(diǎn)?；诤篁?yàn)誤差評(píng)估的方法，給出了顯式算法下Newmark時(shí)域離散誤差的來(lái)源和評(píng)估方法。通過(guò)飽和砂土地震液化響應(yīng)的數(shù)值算例評(píng)估了時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響。數(shù)值分析結(jié)果表明：時(shí)間步長(zhǎng)不同，結(jié)點(diǎn)位移和單元孔壓的時(shí)程曲線明顯不同，同時(shí)計(jì)算耗時(shí)也呈雙曲線關(guān)系；相對(duì)誤差主要分布在變形較大的區(qū)域，全域平均相對(duì)誤差在動(dòng)力響應(yīng)劇烈的時(shí)間段內(nèi)較大。通過(guò)對(duì)計(jì)算時(shí)間步長(zhǎng)和離散誤差的評(píng)估，可有效恰當(dāng)?shù)膶?duì)計(jì)算時(shí)間步長(zhǎng)進(jìn)行取值，也為自動(dòng)步長(zhǎng)調(diào)整提供了依據(jù)。

誤差評(píng)估；砂土液化；堤壩；地震；Newmark法；計(jì)算成本

對(duì)微分方程的求解中常用到隱式算法和顯式算法。隱式方法通過(guò)迭代計(jì)算，降低了時(shí)間步長(zhǎng)大小對(duì)結(jié)果的影響，但計(jì)算過(guò)程非常耗時(shí)，在靜力計(jì)算中應(yīng)用較多。顯式方法則直接根據(jù)當(dāng)前時(shí)間步的狀態(tài)推導(dǎo)出下一時(shí)間步的狀態(tài)，具有快捷性和高效性而被廣泛應(yīng)用于動(dòng)力問(wèn)題的求解中。但是顯式算法的計(jì)算結(jié)果對(duì)時(shí)間步長(zhǎng)的敏感性較高，如何權(quán)衡計(jì)算精度和計(jì)算效率成為使用者較為頭疼的問(wèn)題。因此，對(duì)時(shí)域離散誤差進(jìn)行正確評(píng)估和進(jìn)行時(shí)間步長(zhǎng)自動(dòng)調(diào)整是顯式計(jì)算中的關(guān)鍵技術(shù)。在時(shí)域積分中常用的方法有中心差分法、Newmark法［1］、Houbolt法［2］和Wilson?θ法［3］等。其中Newmark?β時(shí)域離散的方法最為實(shí)用，在結(jié)構(gòu)動(dòng)力學(xué)和砂土地震液化的計(jì)算中被廣泛使用。關(guān)于Newmark法很多學(xué)者曾進(jìn)行了改進(jìn)來(lái)提高計(jì)算精度［4?6］，本文的工作是對(duì)顯式二相耦合計(jì)算中的Newmark?β離散方法進(jìn)行了誤差評(píng)估，并通過(guò)砂土液化的算例進(jìn)行了驗(yàn)證。關(guān)于時(shí)域離散誤差的評(píng)估前人做了很多工作，如Zienkiewicz等［7］在單自由度的動(dòng)力分析中，對(duì)比Taylor級(jí)數(shù)展開(kāi)式獲得了時(shí)域離散的截?cái)嗾`差。Zeng等［8?9］改進(jìn)了Zienkiewicz的方法，提出了一種后驗(yàn)的評(píng)估方法。在巖土工程領(lǐng)域，Sloan等［10?11］探討了時(shí)間步長(zhǎng)對(duì)固結(jié)計(jì)算的影響。

土體在地震作用下的液化流動(dòng)特性及堤壩的地震破壞目前被廣泛研究［12?17］，常用方法也是Ne?wmark顯式計(jì)算方法。本文首先是基于有限元-有限差分（FEM?FDM）的分析方法［18?19］和循環(huán)彈塑性的本構(gòu)模型［20?21］建立了u?p格式的控制方程。然后采用Zeng等［8］提出的后驗(yàn)誤差評(píng)估方法u?p格式顯式計(jì)算進(jìn)行了誤差評(píng)估。同時(shí)也采用堤壩作為數(shù)值算例，對(duì)誤差評(píng)估結(jié)果進(jìn)行驗(yàn)證。

1 控制方程及Newmark時(shí)域離散方法

對(duì)飽和砂土，將加速度的影響作為慣性力，忽略孔隙流體的加速度，可得水-土二相混合體的平衡方程：

式中：σij為應(yīng)力張量，u··si為土顆粒的加速度，bi為體力向量。

對(duì)孔隙水，建立連續(xù)性方程：

式中：n為孔隙率，Kw為水的體積模量，pw為孔隙水壓力，k為滲透系數(shù)，為土體應(yīng)變張量對(duì)時(shí)間的導(dǎo)數(shù)。

在空間域進(jìn)行有限元離散，考慮超孔隙水壓力和阻尼的影響，可建立在動(dòng)力荷載作用下u?p格式的控制方程：

為了簡(jiǎn)便，采用了Rayleigh阻尼，阻尼C取如下形式：

在時(shí)間域采用Newmark?β法對(duì)速度和位移進(jìn)行離散，即獲得t+Δt時(shí)刻與t時(shí)刻變量的關(guān)系：

將式（5）、（6）表示的t+Δt時(shí)刻位移、速度代入式（3）中，得控制方程的最終表達(dá)形式為

2 后驗(yàn)時(shí)域離散誤差評(píng)估

根據(jù)Zeng等［8］的方法，假定t時(shí)刻的變量為已知量，t+Δt時(shí)刻的變量為未知量。且存在τ，τ∈［t，t+Δt］。由式（6），Newmark法在［t，t+Δt］內(nèi)計(jì)算的平均加速度為

克深區(qū)塊單井試油期間開(kāi)展測(cè)試放空氣回收，在鉆井顯示發(fā)現(xiàn)良好后，開(kāi)始組織測(cè)試氣回收，由承包商鋪設(shè)回收管線，并在井場(chǎng)連接好回收氣設(shè)備，試油放噴開(kāi)始后直接開(kāi)展放空氣回收，大大減少天然氣放空，降低耗電量，實(shí)現(xiàn)節(jié)能減排，保護(hù)環(huán)境，取得了良好的效果。

在［t，t+Δt］內(nèi)，更精確的加速度可描述為

在τ時(shí)刻的加速度誤差可表示為u··NM與u··?的差值：

在t+Δt時(shí)刻的位移誤差可用加速度誤差的兩次積分表示：

式（11）給出的為絕對(duì)誤差值，為了表示誤差對(duì)計(jì)算結(jié)果的影響，本文后面提到的誤差都是以相對(duì)誤差的形式表示：

3 堤壩的地震響應(yīng)分析

3.1 有限元模型

本算例為一個(gè)堤壩在地震荷載作用下的動(dòng)力響應(yīng)，計(jì)算模型如圖1所示。底面和左右兩側(cè)為固定位移邊界條件，頂面為排水邊界。結(jié)點(diǎn)N1、N2和單元E1作為計(jì)算的輸出結(jié)點(diǎn)和單元。

圖1 堤壩計(jì)算模型Fig.1 The model of embankment

土體分為干土和飽和砂土2個(gè)部分，采用循環(huán)彈塑性本構(gòu)模型，土體的數(shù)值參數(shù)如表1。

輸入的地震荷載時(shí)程如圖2所示，輸入地震波的步長(zhǎng)為0.01 s，地震波時(shí)長(zhǎng)為22 s。

表1 土體參數(shù)Table1 Soil properties in numerical example

圖2 輸入地震波Fig.2 Input seismic wave

3.2 不同時(shí)間步長(zhǎng)時(shí)堤壩的動(dòng)力響應(yīng)

圖5表示了單元E1在不同計(jì)算時(shí)間步長(zhǎng)時(shí)的超孔隙水壓力比（EPWPR）的時(shí)程曲線。從圖中可以看出，在振動(dòng)過(guò)程中，超孔隙水壓力比從零變?yōu)?，即土體由初始狀態(tài)逐漸變?yōu)橐夯癄顟B(tài)。當(dāng)取不同時(shí)間步長(zhǎng)時(shí)，超孔隙水壓力的發(fā)展是不一致的，初始液化的時(shí)間也不相同，也就是說(shuō)時(shí)間步長(zhǎng)的大小對(duì)孔隙水壓力的變化也有影響。

圖3 不同時(shí)間步長(zhǎng)時(shí)N1點(diǎn)豎向位移時(shí)程曲線Fig.3 Vertical displacement history of point N1in dif?ferent time step

圖4 不同時(shí)間步長(zhǎng)時(shí)N2點(diǎn)水平位移時(shí)程曲線Fig.4 Horizontal displacement history of point N1in different time step

圖5 不同時(shí)間步長(zhǎng)時(shí)單元E1的超孔隙水壓力比時(shí)程曲線Fig.5 EPWPR history of E1in different time step

圖6 表示了地震結(jié)束之后，超孔隙水壓力比EPW?PR的云圖，可見(jiàn)在水平地震作用下，堤壩周圍的飽和砂土層已經(jīng)基本達(dá)到了液化狀態(tài)。圖7表示了堤壩最終的變形圖，與初始網(wǎng)格對(duì)比可見(jiàn)周圍土體的液化使堤壩產(chǎn)生了較大的豎向沉陷和水平擴(kuò)展。

圖6 超孔隙水壓力比云圖Fig.6 Distribution of EPWPR

圖7 網(wǎng)格變形圖（Δt＝1×10-5s）Fig.7 Final deformation of the mesh（Δt＝1×10-5s））

為了對(duì)比不同時(shí)間步長(zhǎng)時(shí)的計(jì)算成本，本文統(tǒng)一采用IBM Intel（R）Xeon（R）3.00 GHz處理器進(jìn)行了計(jì)算，圖8表示了CPU計(jì)算時(shí)間和時(shí)間步長(zhǎng)大小的關(guān)系曲線?？梢钥闯鲭S著時(shí)間步長(zhǎng)的減小，計(jì)算成本成雙曲線式增長(zhǎng)。通過(guò)擬合，得計(jì)算時(shí)間y與時(shí)間步長(zhǎng)x的關(guān)系式如下

可見(jiàn)對(duì)于水土二相耦合顯式計(jì)算，如果采用較大步長(zhǎng)，則會(huì)產(chǎn)生較大的計(jì)算誤差，采用較小的步長(zhǎng)則會(huì)耗費(fèi)計(jì)算時(shí)間。

圖8 不同時(shí)間步長(zhǎng)時(shí)CPU計(jì)算時(shí)間Fig.8 Computation time in different time steps

3.3 時(shí)域離散誤差評(píng)估

為了分析上述不同步長(zhǎng)下的計(jì)算差異，根據(jù)式（11）、（12），計(jì)算了節(jié)點(diǎn)位移的相對(duì)誤差，并繪制了分布云圖，如圖9所示?？梢钥闯?，誤差較大的區(qū)域分布在變形較大的堤壩下部和兩側(cè)。另外，隨著時(shí)間步長(zhǎng)的減小，時(shí)域離散的誤差也隨之減小，且呈現(xiàn)出不同的數(shù)量級(jí)。

為了進(jìn)一步對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行誤差評(píng)估，下式表示了計(jì)算域的平均相對(duì)誤差：

式中：n表示結(jié)點(diǎn)個(gè)數(shù)；j＝1，2表示x和y方向的相對(duì)誤差。

圖10表示了2個(gè)不同時(shí)間步長(zhǎng)時(shí)，平均相對(duì)誤差的時(shí)程曲線?？梢园l(fā)現(xiàn)，在前10 s地震響應(yīng)較劇烈的時(shí)間段內(nèi)，全域平均相對(duì)誤差較大。在不同的時(shí)間步長(zhǎng)下誤差時(shí)程曲線的形狀類似，誤差的大小卻有數(shù)量級(jí)上的差別。

圖10 全域平均相對(duì)誤差時(shí)程曲線Fig.10 Time history of average relative error

4 結(jié)論

本文采用有限元后驗(yàn)誤差計(jì)算的方法對(duì)顯式Newmark?β時(shí)域離散誤差進(jìn)行了評(píng)估，并應(yīng)用于堤壩砂土液化的動(dòng)力數(shù)值分析中。從分析結(jié)果中可得以下結(jié)論：

1）顯式有限元的計(jì)算精度與時(shí)間步長(zhǎng)的大小密切相關(guān)，時(shí)間步長(zhǎng)越小越精確，但花費(fèi)的計(jì)算成本隨時(shí)間步長(zhǎng)的減小成雙曲線增長(zhǎng)。

2）通過(guò)后驗(yàn)誤差評(píng)估的方法，給出了Newmark?β時(shí)域離散誤差的來(lái)源和評(píng)估方法，在砂土液化的數(shù)值算例中給出了相對(duì)誤差的分布和時(shí)程曲線?？梢园l(fā)現(xiàn)時(shí)間步長(zhǎng)越大，相對(duì)誤差則越大；變形較大的區(qū)域和地震響應(yīng)劇烈的時(shí)間段內(nèi)相對(duì)誤差較大。

3）通過(guò)時(shí)域離散相對(duì)誤差的評(píng)估可以幫助我們更客觀地選擇時(shí)間步長(zhǎng)的大小，同時(shí)為自動(dòng)時(shí)間步長(zhǎng)的調(diào)整提供了依據(jù)。

［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］鐘萬(wàn)勰.結(jié)構(gòu)動(dòng)力方程的精細(xì)時(shí)程積分法［J］.大連理工大學(xué)學(xué)報(bào)，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］郭澤英，李青寧，張守軍.結(jié)構(gòu)地震反應(yīng)分析的一種新精細(xì)積分法［J］.工程力學(xué)，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］李爽，翟長(zhǎng)海，劉洪波，等.Newmark積分方法在負(fù)剛度時(shí)的數(shù)值穩(wěn)定性分析［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，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］阮光雄，傅少君，陳勝宏.固結(jié)問(wèn)題有限元法時(shí)步自適應(yīng)研究［J］.巖土力學(xué)，2005，26（4）：591?599.NGUYEN Quanghung，F(xiàn)U 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］汪明武.可液化場(chǎng)地土工抗震性能動(dòng)態(tài)離心機(jī)試驗(yàn)與數(shù)值模擬［M］.北京：科學(xué)出版社，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］蘭景巖，劉紅帥，呂悅軍.渤海土類動(dòng)力非線性參數(shù)及合理性［J］.哈爾濱工程大學(xué)學(xué)報(bào)，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，F(xiàn)U 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.網(wǎng)絡(luò)出版時(shí)間：2015?01?09.

國(guó)家863計(jì)劃資助項(xiàng)目（2012AA112510）.

張西文（1987?），男，博士研究生；

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

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