砂土地震液化分析中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.