凌道盛，鞏師林，胡成寶，鈕家軍

(1.浙江大學巖土工程研究所，浙江，杭州 310058；2.浙江大學軟弱土與環(huán)境土工教育部重點實驗室，浙江，杭州 310058；3.浙江大學寧波理工學院土木建筑工程學院，浙江，寧波 315100)

包含大量孔隙、節(jié)理甚至斷層的非連續(xù)巖土體結(jié)構(gòu)在自然界中廣泛存在，巖土體的非連續(xù)變形分析日益引起人們的關(guān)注。石根華[1－2]提出的非連續(xù)變形分析方法(discontinuous deformation analysis,DDA)是非連續(xù)變形分析中最有效的方法之一[3－4]，已經(jīng)廣泛應用于地震引起的山體滑坡分析、巖石爆破過程模擬、落石軌跡追蹤、隧道及采礦、巖石邊坡穩(wěn)定性評價及其他方面[5－12]。

原始DDA假設(shè)塊體滿足小應變、大轉(zhuǎn)動，將塊體位移增量分解為常應變對應的變形、剛體平動和剛體轉(zhuǎn)動三部分之和，即位移增量分量直接相加。求解時以上一時步構(gòu)型作為參照構(gòu)型，在時間步內(nèi)將塊體位移增量線性化，并利用求得的線性化位移增量更新塊體構(gòu)型。研究表明，原始DDA的位移模式采用全一階近似和小角度假定，導致塊體體積自由膨脹、應變分量物理意義不明確、應力場畸變等問題[13]，給數(shù)值計算帶來較大誤差。同時，采用一階近似后的位移模式推導塊體相關(guān)加速度表達式，忽略了塊體轉(zhuǎn)動時離心力與科氏力的作用。針對上述問題，MacLaughlin等[14]將一階近似前位移公式中余弦值二階泰勒級數(shù)展開，當每時步轉(zhuǎn)角小于 0.4 rad時塊體體積自由膨脹產(chǎn)生的誤差可以接受。Cheng等[15]以轉(zhuǎn)角增量的正弦作為基本未知量，基于三角函數(shù)恒等式改寫了有限轉(zhuǎn)動位移增量公式，并利用上一時步轉(zhuǎn)角的增量簡化形函數(shù)矩陣，一定程度上克服了塊體體積不合理膨脹現(xiàn)象，在塊體勻速轉(zhuǎn)動時可獲得精確解。Ke[16]及Koo等[17]在每一時步計算完成后采用有限轉(zhuǎn)動位移公式對塊體位移進行后修正(post-adjustment)。高亞楠等[18]采用有限變形理論將塊體運動過程分解為變形和轉(zhuǎn)動兩部分，并根據(jù)有限變形幾何場更新塊體構(gòu)型，可消除塊體轉(zhuǎn)動帶來的自由膨脹。Jiang等[19]在每一塊體中固定隨動坐標系，應變分量增量僅在隨動坐標系下疊加，利用有限轉(zhuǎn)動位移全量公式計算塊體全量位移，更新當前構(gòu)型，避免了塊體大轉(zhuǎn)動時的虛假體積膨脹及應力場畸變。Lin等[20]基于Green應變和第二類 Piola-Kirchhoff應力，給出采用二階多項式導出的塊體有限變形公式，可模擬大應變、有限轉(zhuǎn)動問題。Fan等[21]將塊體的運動進行S-R(strain-rotation)分解以同時獲取塊體的變形與轉(zhuǎn)動，采用原有DDA的形函數(shù)構(gòu)造控制方程，克服了塊體體積自由膨脹問題，并賦予其處理大變形問題的能力。

綜上可見，現(xiàn)有DDA及其改進方法大都基于原始DDA的位移模式及構(gòu)型更新方法，在此基礎(chǔ)上，本文分析了原始DDA存在的若干問題，包括小角度假設(shè)及線性化位移增量導致的塊體體積膨脹、應變分量直接疊加導致的應變場畸變及其物理意義不明確以及采用線性化后的位移增量公式推導塊體的加速度導致忽略塊體旋轉(zhuǎn)的離心力和科氏力作用。通過算例證實了上述若干問題的存在，為DDA的進一步發(fā)展，從而為真實準確地模擬巖土工程中的具體問題提供參考。

1 原始DDA位移模式

石根華提出的 DDA假設(shè)塊體為常應變，且滿足小變形假設(shè)。不失一般性，假定tn時刻塊體的構(gòu)型已知，物質(zhì)點 {X,Y}T的空間坐標為 {xn,yn}T。以tn時刻為參照構(gòu)型，從tn時刻到tn+1時刻的位移增量{Δun, Δvn}T可分解為三部分，即剛體平動、繞參考點的剛體轉(zhuǎn)動和常應變對應的變形。根據(jù)原始DDA的表述，位移增量由上述三部分直接疊加，表示為：

DDA的研究對象以巖石為主，在小應變、小轉(zhuǎn)動假設(shè)的前提下，塊體的幾何方程表示為：

塊體在滿足小變形、小轉(zhuǎn)動假定后，式(1)、式(2)可線性化為：

將上式，采用矢量形式表示：

式中，Tn、 Δdn分別為塊體在tn時刻的形函數(shù)矩陣與廣義位移矢量增量。分別表示為：

類似地，tn+1時刻塊體內(nèi)任意物質(zhì)點的加速度表示為：

式中，圓點頂標“??”表示物理量關(guān)于時間的二階導數(shù)。

開閉迭代求得塊體廣義位移增量后，原始DDA通過簡單累加求得tn+1時刻塊體的廣義位移矢量總量dn+1和塊體任意點位移un+1：

2 原始DDA存在的若干問題

2.1 轉(zhuǎn)角增量的一階近似導致大轉(zhuǎn)動時塊體體積自由膨脹

建立在小變形、小轉(zhuǎn)動假設(shè)之上的原始 DDA在模擬邊坡滑動、巖石崩塌等實際工程災害時，作為研究對象的塊體常涉及大轉(zhuǎn)動問題，此時會觀察到塊體體積出現(xiàn)不合理膨脹現(xiàn)象。

在tn時刻，塊體中任意點P{xn,yn}T繞參考點沿著x軸、y軸位移增量的精確解{Δun, Δvn}T，見式(1)、式(2)。假定研究對象為剛體，此時：

式(1)、式(2)可簡化為：

式(5)、式(6)簡化為：

將剛體轉(zhuǎn)動位移增量精確解式(14)、式(15)和近似解式(16)、式(17)進行幾何對比，如圖1所示，線段P0P旋轉(zhuǎn)后長度為P0P′，且P0P=P0P′。而小角度假定條件下的近似解P0P′導致線段沿著r向量發(fā)生偏移，此時P0P＜P0P′。當每時步內(nèi)的轉(zhuǎn)角增量較大時，誤差逐步累積，最終導致塊體體積產(chǎn)生自由膨脹。

圖1 剛體轉(zhuǎn)動位移精確解與近似解幾何圖Fig.1 Gemetry of approximation and the exact solution to rigid block rotation

圖2 小角度假設(shè)的相對誤差Fig.2 Relative error induced by small rotation assumption

2.2 應變增量的直接疊加導致塊體應力、應變場畸變

現(xiàn)階段，針對原始DDA采用線性化位移增量更新塊體構(gòu)型引起的塊體體積自由膨脹現(xiàn)象，一些學者提出了不同的修正方法。雖然不同修正方法可在一定程度上克服塊體體積的自由膨脹，但是模擬塊體大轉(zhuǎn)動問題時，塊體應力及應變計算仍會引入較大誤差。

原始DDA模擬的塊體發(fā)生有限轉(zhuǎn)動時，開閉迭代結(jié)束后得到的變形變量增量Δd可分為兩類。一類是當塊體發(fā)生有限轉(zhuǎn)動時可直接疊加的水平和豎向位移增量以及剛體轉(zhuǎn)角增量以上三種塊體的變形變量分量的增量當塊體發(fā)生有限轉(zhuǎn)動時可直接疊加。另一類為小應變假設(shè)前提下，塊體另外三項應變相關(guān)的變形變量分量的增量在塊體發(fā)生有限轉(zhuǎn)動時，應變分量增量的直接疊加會帶來明顯的誤差。

而原始DDA及大多數(shù)修正后的DDA程序在開閉迭代后，應變分量增量直接疊加獲得應變分量總量，即：

假定塊體為線彈性體，原始DDA根據(jù)塊體應變求得應力，并將求得的應力作為塊體初始應力加到總體平衡方程中去。在平面應力條件下，nt時刻應力-應變關(guān)系滿足：

式中：E、μ分別為塊體的彈性模量及泊松比；{σx,σy,τxy}nT及 {σx,σy,τxy}n-1T分別為tn、tn-1時刻的塊體應力分量總量。

同樣，原始DDA中每時步的初速度為上一時步結(jié)束后的末速度，在時間步長為Δt條件下，塊體在tn時步內(nèi)的速度與tn-1時步內(nèi)的速度有如下關(guān)系：

然而，開-閉迭代后求得的應變分量增量均是基于局部坐標系，即圖 3(b)中的X′OY′坐標系。若應變分量增量直接進行疊加而不作坐標系旋轉(zhuǎn)，不同時步應變增量的參照構(gòu)型不同，塊體內(nèi)應變物理意義不明確。且若將塊體的應變限制為0，由式(19)、式(20)，當塊體旋轉(zhuǎn)時，塊體內(nèi)應力場及速度場將產(chǎn)生畸變。

圖3 塊體應變示意圖Fig.3 Sketch of block strain

2.3 采用加法分解的位移增量公式忽略塊體旋轉(zhuǎn)過程中的離心力與科氏力

Jiang在原始DDA基礎(chǔ)上，在每個塊體中固定隨動局部坐標系，開閉迭代后求得的應變增量分量旋轉(zhuǎn)至局部坐標系下疊加，并采用有限轉(zhuǎn)動位移總量公式計算塊體總位移，更新塊體構(gòu)型。不僅克服了大轉(zhuǎn)動時塊體體積的自由膨脹，而且避免了塊體應變的累計誤差，計算出的應變更為合理。

DDA通過選取不同動力系數(shù)gg(gg=0～1)，即塊體當前時步的初速度與上一時步末速度的比值，可同時處理靜力學與動力學問題。對于動力學問題，塊體運動過程中的慣性項不可忽略。DDA將加速度引起的慣性力考慮在內(nèi)，并表示為：

在塊體位移增量精確表示式(1)、式(2)的基礎(chǔ)上，原始DDA基于小角度假設(shè)，將式(1)、式(2)一階近似為式(5)、式(6)，并以一階近似后的式(5)、式(6)推導塊體相關(guān)子矩陣。慣性力子矩陣的推導與塊體加速度有關(guān)，而原始DDA中塊體加速度由式(5)、式(6)直接導出，即：

塊體旋轉(zhuǎn)過程會產(chǎn)生離心力及科氏力，原始DDA采用一階近似后的線性表達式推導塊體的加速度表達式，而塊體加速度應是從線性化前的位移增量表達式導出后再線性化，直接采用線性化后的位移增量表達式忽略了塊體旋轉(zhuǎn)時離心力和科氏力的作用。

3 算例驗證

設(shè)計兩個典型算例來說明原始 DDA采用一階近似后的位移增量簡單疊加并推導塊體相應的加速度表達式，及開閉迭代后應變增量的直接疊加，在模擬塊體大轉(zhuǎn)動時帶來的相關(guān)問題。

3.1 正八邊形繞定點擺動問題

如圖4所示，對角線長10 m、繞頂點E擺動的正八邊形，質(zhì)量密度為 2800 kg/m3，彈性模量為2 GPa，泊松比為0.25。以八邊形質(zhì)心為坐標原點O，如圖示建立坐標系XOY。

采用一個塊體模擬重力作用下八邊形繞頂點E的擺動，取時間步長Δt=0.002 s。圖 5給出了原始DDA程序計算的該八邊形面積隨計算時步的變化曲線。由圖 5可見，隨著計算時步的增加，原始DDA計算出正八邊形面積發(fā)生較大膨脹，該正八邊形的初始面積為70.7 m2，計算到5000步時，面積已經(jīng)達到72.4 m2，膨脹了2.45%，這一誤差在工程分析中是不可接受的。

表1給出了正八邊形塊體在200步、2000步及5000步時，四條對角線長度變化及相應的線應變，表2給出相應時步塊體應變及主應變計算的結(jié)果。

圖4 正八邊形塊體繞固定點擺動Fig.4 A swing regular octagon fixed at one vertex point

圖5 八邊形面積隨計算時步變化曲線Fig.5 Change of octagon area with time step

表1 原始DDA計算的各對角線長度及應變Table 1 The length changes of four segments on the swing octagon by original DDA

由表1可得，隨著計算時步的增加，當計算至5000步時，四條對角線伸長量均較大，對應的線應變均在1.2×10-2。而5000步時，原始DDA計算的該正八邊形塊體最大、最小主應變分別 3.32×10-5及-2.60×10-5，此時四條對角線的線應變均未處于塊體最大和最小主應變范圍內(nèi)的情況。這一現(xiàn)象隨著計算時步的增大愈加明顯。說明原始DDA在處理大轉(zhuǎn)動問題時不僅會帶來塊體體積的自由膨脹，而且應變的直接疊加導致應變分量物理意義不明確，無法正確反映塊體的實際變形。

3.2 矩形塊體受均布荷載勻速旋轉(zhuǎn)問題

如圖6所示，以角速度ω= 0 .5π r ad/s 繞質(zhì)心勻速旋轉(zhuǎn)的等厚矩形薄板，長2L=6 m，寬2b=4 m。材料質(zhì)量密度為ρ=2400kg/m3，彈性模量E=100 MPa，塊體受大小為1 MPa的水平初始應力。采用一個塊體模擬該矩形薄板的運動，假定塊體不受重力作用，為消除應變對塊體內(nèi)應力場、速度場的影響，假定每時步內(nèi)的應變分量增量

表2 原始DDA計算的主應變Table 2 Principal strains by original DDA

圖6 矩形塊體的幾何形狀及受力示意圖Fig.6 Geometry and loading of the rectangular block

圖7 原始DDA計算出的塊體應力Fig.7 Stresses calculated by original DDA

取時間步長Δt=1 s，原始DDA計算出的塊體應力如圖7。由圖7可見，塊體的水平應力、豎向應力及剪應力與旋轉(zhuǎn)的累積角度無關(guān)。表明原始DDA在計算轉(zhuǎn)動塊體時，計算出的錯誤應力方向直接影響塊體內(nèi)應力分布及彈性形變。

3.3 矩形塊體繞質(zhì)心勻速旋轉(zhuǎn)問題

如圖8所示為以角速度ω=0.2rad/s 繞質(zhì)心勻速旋轉(zhuǎn)的等厚度矩形薄板，長2L=6 m，寬2b=4 m。材料質(zhì)量密度為ρ=2800kg/m3，彈性模量為E=2 GPa。為方便起見，假定矩形薄板處于平面應力狀態(tài)，泊松比μ=0。以薄板質(zhì)心為坐標原點O，如圖示建立坐標系XOY。矩形薄板內(nèi)任意點的正應變分量為：

定義薄板平均正應變?yōu)椋?/p>

圖8 繞質(zhì)心勻速旋轉(zhuǎn)矩形塊體Fig.8 A rectangular block rotateing uniformly around the gravity center

采用一個塊體模擬該矩形薄板的運動。圖9給出時間步長 Δt= 0 .002s 時，原始 DDA及式(25)計算出的應變。

圖9 不同方法應變對比Fig.9 Strains calculated by different methods

由圖9可見，原始DDA不能反映離心力作用下薄板的應變，不能計算出由式(25)給出的塊體平均應變解析解，給出了錯誤的應變。

4 結(jié)論

原始DDA的位移增量公式由塊體的平動、轉(zhuǎn)動及變形直接疊加而得，并假定塊體每時步的轉(zhuǎn)角增量足夠小，由此得到線性化后的位移增量公式及相應形函數(shù)矩陣。而當模擬的塊體產(chǎn)生轉(zhuǎn)動時，塊體體積出現(xiàn)錯誤膨脹，且塊體內(nèi)應變(應力)場發(fā)生畸變，從而影響DDA的模擬準確性及精度。

通過對原始DDA的位移模式、構(gòu)型更新方式等分析，并通過數(shù)值算例驗證原始DDA方法存在如下不足：

(1) 位移增量的簡單疊加僅適用小轉(zhuǎn)動條件；

(2) 應變分量增量是相對當前時刻定義的，不同時步應變增量參照構(gòu)型不同，簡單的累加導致應變分量物理意義不明確；

(3) 位移增量求解采用線性化后的近似表達式，隨著計算時步的增加，引起誤差累積，導致塊體出現(xiàn)體積自由膨脹；

(4) 塊體加速度應是從線性化前的位移增量表達式導出后再線性化，采用線性化后的位移增量公式忽略了塊體轉(zhuǎn)動引起的離心力和科氏力作用。