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

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

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

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

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

1 原始DDA位移模式

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

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

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

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

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

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

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

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

2 原始DDA存在的若干問題

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

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

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

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

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

將剛體轉(zhuǎn)動位移增量精確解式(14)、式(15)和近似解式(16)、式(17)進(jìn)行幾何對比，如圖1所示，線段P0P旋轉(zhuǎn)后長度為P0P′，且P0P=P0P′。而小角度假定條件下的近似解P0P′導(dǎo)致線段沿著r向量發(fā)生偏移，此時P0P＜P0P′。當(dāng)每時步內(nèi)的轉(zhuǎn)角增量較大時，誤差逐步累積，最終導(dǎo)致塊體體積產(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 應(yīng)變增量的直接疊加導(dǎo)致塊體應(yīng)力、應(yīng)變場畸變

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

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

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

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

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

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

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

圖3 塊體應(yīng)變示意圖Fig.3 Sketch of block strain

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

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

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

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

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

3 算例驗(yàn)證

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

3.1 正八邊形繞定點(diǎn)擺動問題

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

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

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

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

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

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

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

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的水平初始應(yīng)力。采用一個塊體模擬該矩形薄板的運(yùn)動，假定塊體不受重力作用，為消除應(yīng)變對塊體內(nèi)應(yīng)力場、速度場的影響，假定每時步內(nèi)的應(yīng)變分量增量

表2 原始DDA計算的主應(yīng)變Table 2 Principal strains by original DDA

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

圖7 原始DDA計算出的塊體應(yīng)力Fig.7 Stresses calculated by original DDA

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

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

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

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

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

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

圖9 不同方法應(yīng)變對比Fig.9 Strains calculated by different methods

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

4 結(jié)論

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

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

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

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

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

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