基于有限元-近場動力學(xué)耦合方法的擋土墻變形破損分析*

牟太平，孫偉，邱志雄，肖紅萍，張嘎，林沛元

1. 廣東省高速公路有限公司，廣東 廣州 510623

2. 中山大學(xué)土木工程學(xué)院，廣東 廣州 510275

3. 北京城投地下空間開發(fā)建設(shè)有限公司，北京 100044

4. 清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084

高速公路改擴建中的路堤常采用擋土墻進行支擋處理。在背后填土土壓力和外部荷載（如車輛荷載）作用下，擋土墻可能發(fā)生整體傾斜失穩(wěn)甚至倒塌等變形破壞。特別是對于剛性墻，可能會由于應(yīng)力過大而導(dǎo)致墻身局部出現(xiàn)裂縫等。

目前，有學(xué)者采用模型試驗、理論分析或數(shù)值模擬等方法，對擋土墻應(yīng)力變形、土壓力計算方法，加筋擋土墻性能等進行了研究。章為民等［1］進行了加筋擋土墻離心模型試驗；宋飛等［2］設(shè)計了可考慮剛性擋墻平動位移模式的土壓力離心模型設(shè)備；李永剛等［3］引入主應(yīng)力拱的概念，利用水平層分析方法，研究了墻背豎直且填土無黏性條件下的土壓力分布。通過理論分析，劉洋等［4］得到了滑動土體任一點的應(yīng)力解，以及擋土墻土壓力的大小、分布及作用點位置等?？紤]地基土和填土的流變以及格柵的蠕變性，欒茂田等［5］采用非線性有限元法（FEM，finite element method）分析了眾多因素對擋土墻應(yīng)力變形性能的影響；采用有限差分法，鄒維列等［6］研究了重力式加筋土擋土墻的工作性能和土壓力分布特點。Rao 等［7］提出了一種平動模式下路肩擋土墻主動土壓力的簡化計算方法。并與試驗數(shù)據(jù)、庫侖和朗肯理論等的預(yù)測結(jié)果進行對比，驗證了方法的有效性。

針對擋土墻的失穩(wěn)破壞，也有學(xué)者進行了研究。張波等［8］采用有限元法分析了黏土加筋擋土墻的破壞過程，獲得的滑裂面與實際破壞形式較為相符。采用有限差分法，吳順川等［9］研究了車輛動荷載作用下路基加筋擋土墻的失穩(wěn)機理和破壞模式。孟云偉等［10］揭示了墻背階梯式石籠擋土墻的破壞模式主要是傾覆破壞。Wang 等［11］采用ABAQUS 建立有限元分析模型，開展土釘和預(yù)應(yīng)力錨桿加固復(fù)合支護結(jié)構(gòu)的失效機理研究，研究發(fā)現(xiàn)，預(yù)應(yīng)力錨桿的施加使基坑變形減小，提高了高坡度的整體穩(wěn)定性。

然而，關(guān)于剛性擋土墻在外部荷載作用下墻身可能發(fā)生局部破損的問題，目前的研究還相當(dāng)不足。本文以某改擴建工程為例，提出了一種有限元-近場動力學(xué)（FEM-NOSBPD）耦合方法，并將其應(yīng)用于算例的懸臂式L 型擋土墻的變形破損分析。

1 有限元-近場動力學(xué)耦合方法

有限元-近場動力學(xué)耦合方法［12-13］基于疊加多尺度理論［14］，在耦合區(qū)域內(nèi)將總位移分解為有限元表達的位移和近場動力學(xué)表達的位移，并在耦合區(qū)域邊界施加齊次邊界條件，通過變分原理實現(xiàn)兩者的有效耦合。在該方法中，僅需要在結(jié)構(gòu)破損的較小區(qū)域采用近場動力學(xué)離散，而在其他大部分區(qū)域采用計算量較小的有限元方法離散，可實現(xiàn)大型結(jié)構(gòu)變形破損的全過程精確模擬。

近場動力學(xué)（PD，peridynamics）是一種新的非局部計算力學(xué)理論［15］，可以自發(fā)地判斷裂紋擴展方向，不需要記錄和分析裂紋網(wǎng)格的幾何關(guān)系；并可同時描述多條裂紋，能夠方便地模擬材料的破壞過程，并可有效避免破損模擬中的單元尺寸效應(yīng)問題，從而獲得與離散密度無關(guān)的客觀解。2007 年，Silling 提出了常規(guī)（ordinary state-based）和非常規(guī)態(tài)型（non-ordinary state-based） PD 理論［16］。本文采用NOSBPD理論，其運動方程為

式中ρ是質(zhì)量密度，u是位移，b是給定的體積力，Hx表示點x的鄰域，它是指初始參考構(gòu)型中截斷半徑δ所包圍的區(qū)域，f表示對點力函數(shù)，[x，t]＜x' -x＞-[x'，t]＜x-x' ＞表示作用在點x上由點x'給的力密度向量；[ ]x，t和[x'，t]表示力向量態(tài)。初始位置參考構(gòu)形的向量態(tài)為

某一時刻，變形向量態(tài)為

給定某一經(jīng)典連續(xù)介質(zhì)本構(gòu)關(guān)系，其應(yīng)力張量可通過式（5）定義的非局部變形梯度求得。經(jīng)過推導(dǎo)，力向量態(tài)可以用應(yīng)力σ和變形梯度F表示為

式（6）建立了NOSBPD 中的力態(tài)與經(jīng)典連續(xù)介質(zhì)力學(xué)理論中的應(yīng)力、變形梯度等物理量之間的關(guān)系。

孫偉從變分角度嚴格證明了有限元和非常規(guī)態(tài)型近場動力學(xué)理論（NOSBPD，non-ordinary state-based peridynamics）［16］兩者內(nèi)力虛功不同形式之間的等效關(guān)系，然后推導(dǎo)了可以保持兩者各自求解格式不變的總的控制方程和耦合切向剛度矩陣（即有限元采用網(wǎng)格離散和弱形式求解；近場動力學(xué)采用粒子離散和強形式求解），從而解決了兩者離散和求解方式不同的問題，提出了一種新的耦合界面處理技術(shù)以消除或減小界面處的虛假的相互作用力，從而解決了兩者耦合界面的問題［17］，建立了適用于結(jié)構(gòu)破損分析的有限元-近場動力學(xué)耦合分析方法。

本文采用FEM-NOSBPD 耦合分析方法進行破損計算。首先，進行裂縫萌生起裂判斷，而后進行裂縫擴展分析。裂縫的萌生起裂可由材料的抗拉、抗壓或抗剪強度來判斷，這里采用抗拉強度進行起裂判斷。對于裂縫的擴展，則將其視為斷裂力學(xué)問題，采用基于能量的破損準(zhǔn)則［18-19］，即當(dāng)近場動力學(xué)中的單個鍵儲存的能量密度達到某一閾值（與材料的斷裂能有關(guān)），則該鍵發(fā)生斷裂。本文采用隱式方法進行耦合方程的求解。隱式計算過程中，每一荷載步的每一增量步應(yīng)迭代至該步完全沒有新的鍵發(fā)生斷裂為止，方可進入下一增量步。

接下來，采用混凝土偏心三點彎曲梁算例［20］來驗證FEM-NOSBPD 耦合方法在混凝土結(jié)構(gòu)破損模擬方面的有效性。 混凝土計算參數(shù)為E= 38 GPa，υ= 0.2，ft= 3.0 MPa，GI= 69 N/m，近場動力學(xué)離散參數(shù)為Δx= 1 mm，δ= 3Δx。圖1為試驗布置和幾何尺寸。圖2為耦合方法的計算模型和最終裂縫形態(tài)。圖3 為加載點承載力-位移曲線。分析可知，無論是裂縫形態(tài)還是承載力曲線，耦合方法的計算結(jié)果都能與試驗值吻合較好。因此，本文建立的FEM-NOSBPD 耦合方法能夠有效地模擬混凝土結(jié)構(gòu)破損。

圖1 混凝土偏心三點彎曲梁試驗［20］（單位：mm）Fig.1 Three-point bending eccentrically concrete beam test[20](unit:mm)

圖2 FEM-PD耦合方法計算模型和計算結(jié)果Fig.2 The computational model of the FEM-PD coupling method and computational results

圖3 加載點承載力-位移曲線Fig.3 The curve of bearing capacity versus displacement

2 計算模型

以某高速公路改擴建工程為例進行分析，計算模型如圖4 所示。懸臂式L 型擋土墻高5 m，底寬3 m，懸臂寬度0.27 m。根據(jù)相關(guān)規(guī)范［21］，設(shè)定車輛等效荷載為20 kPa。土體采用摩爾-庫倫模型模擬。計算參數(shù)基于現(xiàn)場地勘資料確定，如表1所示。擋土墻為C30混凝土，本文計算暫不考慮擋土墻的配筋問題，計算參數(shù)如表2 所示。其中，ft，fc和GI分別表示混凝土單軸抗拉強度、抗壓強度和I 型斷裂能。斷裂能參數(shù)取自文獻［22］，其他參數(shù)取值來自混凝土結(jié)構(gòu)設(shè)計規(guī)范［23］。擋土墻和地基土及填土之間設(shè)置了無厚度的Goodman 單元，接觸本構(gòu)采用Clough-Duncan 雙曲線模型，計算參數(shù)如表3所示（表中參數(shù)取自文獻［17］）。近場動力學(xué)粒子離散參數(shù)為粒子間距5 mm，鄰域半徑3 倍粒子間距，粒子僅需在破損位置附近布置。計算邊界條件為底邊雙向固定、側(cè)邊固定法向，計算按照平面應(yīng)變假設(shè)進行。計算中，應(yīng)力均以壓為正、拉為負，主應(yīng)力滿足σ1＞σ2＞σ3。

圖4 計算模型和邊界條件Fig.4 Calculation model and boundary conditions

表1 土體計算參數(shù)Table 1 Calculation parameters of the soil

表2 擋土墻變形破損計算參數(shù)Table 2 Parameters for the deformation and failure analysis of the retaining wall

表3 接觸面計算參數(shù)Table 3 Parameters of the interface

3 擋土墻應(yīng)力變形及破損分析

計算得到的土體水平位移、豎直位移及塑性區(qū)分布如圖5 所示。由圖5 可知，由于回填土是無黏性砂土，因此在擋土墻墻頂附近位移較大，最大水平位移約為1 cm，最大豎直位移約為2.8 cm，同時在該位置還出現(xiàn)了一定的塑性區(qū)，墻體和土體在該位置的接觸面也有一定的滑動或脫開。

圖5 土體計算結(jié)果（變形與塑性區(qū)分布）Fig.5 Computational results of the soil deformation and plastic zone distribution

擋土墻的變形圖和位移如圖6～7 所示。由圖可知，在本文的計算條件下，墻體處于繞墻底轉(zhuǎn)動模式。沿著擋土墻高度方向，豎直位移變化不大，水平位移則是隨著高程增加而增大，墻頂達到最大。擋土墻的應(yīng)力分布如圖8 所示。由圖8 可知，1）墻體在豎直方向主要承受墻體和土體不均勻沉降所帶來的摩擦力以及墻體自身的重力，而水平方向則是墻體一側(cè)填土帶來的土壓力。從應(yīng)力分布來看，墻底壓應(yīng)力并不是很大，最大僅2.24 MPa；但墻底由于水平荷載帶來的彎矩所引起的拉應(yīng)力卻較大，最大約5.47 MPa，出現(xiàn)在L 型轉(zhuǎn)折交叉內(nèi)側(cè)處。因此，懸臂式擋土墻結(jié)構(gòu)以承受水平荷載為主，為受彎構(gòu)件。2）墻體應(yīng)力較大的位置主要分布在L 型轉(zhuǎn)折處，內(nèi)側(cè)受拉、外側(cè)受壓；該位置是易破損位置，實際施工中可在該位置考慮配置一定的鋼筋。

圖6 變形前后的擋土墻Fig.6 The retaining wall before and after deformation

圖7 擋土墻位移分布Fig.7 Distribution of displacement of the retaining wall

圖8 擋土墻應(yīng)力分布Fig.8 Distribution of stress of the retaining wall

圖9給出了擋土墻的破損分布。結(jié)果表明，在L型擋墻轉(zhuǎn)折交叉內(nèi)側(cè)處出現(xiàn)了一條裂縫，裂縫走向為水平向下約20°，裂縫長度約為6 cm?？傮w上看裂縫擴展長度和破損區(qū)域并不是很大，這可能與拉應(yīng)力分布區(qū)域較小，僅集中在轉(zhuǎn)折交叉位置有關(guān)。

圖9 擋土墻破損分布（單位：m）Fig.9 Distribution of damage of the retaining wall(unit:m)

采用耦合模型進行分析時，近場動力學(xué)某粒子點的鍵達到設(shè)定閾值后將發(fā)生斷裂，從而形成裂縫。在裂縫尖端附近存在著明顯的應(yīng)力集中和重分布，使得該點其他的鍵或者周邊點的鍵可能進入破壞狀態(tài)，從而促進裂縫向前擴展。圖10 為加載結(jié)束時刻破損區(qū)裂縫附近的大小主應(yīng)力分布，裂縫尖端存在奇異性，應(yīng)力較大，而其他位置的超額應(yīng)力自動發(fā)生了轉(zhuǎn)移，應(yīng)力較小，即擋土墻破損區(qū)存在明顯的應(yīng)力集中和應(yīng)力重分布。

圖10 擋土墻破損區(qū)主應(yīng)力分布Fig.10 Distribution of principle stress in the damage zone of the retaining wall

4 扶壁影響分析

扶壁式擋土墻是路堤擋土、碼頭岸壁的一種主要形式［24］。通過計算，本文探討了扶壁對于擋墻應(yīng)力和破損的影響。扶壁的具體尺寸如圖11 所示，其他計算條件與前述懸臂式擋土墻相同。圖12 給出了扶壁式擋土墻的應(yīng)力分布。對比圖8 可知：（1）在應(yīng)力分布規(guī)律上，兩者基本相似。應(yīng)力較大的位置仍主要分布在扶壁底部的轉(zhuǎn)折處，內(nèi)側(cè)受拉、外側(cè)受壓，即擋墻易破損位置沒有變化。（2）在應(yīng)力大小上，扶壁式擋土墻相比于懸臂式減小很多，如：扶壁式擋土墻拉應(yīng)力最大不超過0.6 MPa，與懸臂式最大拉應(yīng)力5.47 MPa 相比，減小超過80%；扶壁式擋土墻壓應(yīng)力最大不超過0.8 MPa，與懸臂式最大壓應(yīng)力2.24 MPa 相比，減小約60%。扶壁式擋土墻沒有出現(xiàn)裂縫，應(yīng)力的改善是與增設(shè)了扶壁、截面的抗彎剛度得到增強有關(guān)。因此，增加扶壁是改善擋土墻受力條件、防止發(fā)生破損的有效措施。

圖11 扶壁式擋土墻的幾何尺寸（單位：m）Fig.11 Dimensions of the counterfort retaining wall(unit:m)

圖12 扶壁式擋土墻應(yīng)力分布Fig.12 Distribution of stress of the counterfort retaining wall

5 結(jié) 論

采用有限元-近場動力學(xué)耦合方法，分析了某路基懸臂式L型擋土墻的變形破損狀況。主要結(jié)論如下：

1）通過混凝土偏心三點彎曲梁算例，發(fā)現(xiàn)耦合方法的計算結(jié)果能較好地再現(xiàn)試驗裂縫形態(tài)，且承載力曲線與試驗結(jié)果吻合，這說明本文采用的FEM-NOSBPD 耦合方法可有效地模擬混凝土結(jié)構(gòu)損傷問題。

2）懸臂式擋土墻處于繞墻底轉(zhuǎn)動模式。水平方向撓度隨著墻高的增高而增大，墻頂達到最大；擋墻墻底的壓應(yīng)力較小，但墻底由于水平荷載帶來的彎矩所引起的拉應(yīng)力卻較大，出現(xiàn)在L型轉(zhuǎn)折交叉內(nèi)側(cè)處。因此，懸臂式擋土墻結(jié)構(gòu)以承受水平荷載為主，為受彎構(gòu)件。

3）懸臂式擋土墻易破損位置是L 型轉(zhuǎn)折處，其內(nèi)側(cè)受拉、外側(cè)受壓，實際施工中可在該位置考慮配置一定的鋼筋。破損分析表明，在L型轉(zhuǎn)折交叉內(nèi)側(cè)處出現(xiàn)了一條裂縫，裂縫走向為水平向下約20°，裂縫長度約為6 cm?？傮w上看，裂縫擴展長度和破損區(qū)域并不是很大，破損區(qū)存在明顯的應(yīng)力集中和應(yīng)力重分布。

4）同等計算條件下，扶壁式擋土墻與懸臂式擋土墻應(yīng)力分布規(guī)律基本相似，擋墻易破損位置沒有變化。然而，扶壁式擋土墻由于截面抗彎剛度的增大，應(yīng)力得到極大改善，拉應(yīng)力最大值減小超過80%，壓應(yīng)力最大值減小約60%，擋墻沒有出現(xiàn)裂縫。因此，增加扶壁是改善擋土墻受力條件、防止發(fā)生破損的有效措施。

5）本文采用的分析方法僅適用于二維靜態(tài)或準(zhǔn)靜態(tài)分析，對更復(fù)雜的三維或動力分析，仍需要加以研究。