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

(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)用于地震引起的山體滑坡分析、巖石爆破過(guò)程模擬、落石軌跡追蹤、隧道及采礦、巖石邊坡穩(wěn)定性評(píng)價(jià)及其他方面[5－12]。

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

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

1 原始DDA位移模式

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

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

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

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

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

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

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

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

2 原始DDA存在的若干問(wèn)題

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

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

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

式(1)、式(2)可簡(jiǎn)化為：

式(5)、式(6)簡(jiǎn)化為：

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

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

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

2.2 應(yīng)變?cè)隽康闹苯盈B加導(dǎo)致塊體應(yīng)力、應(yīng)變場(chǎng)畸變

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

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

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

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

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

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

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

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

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

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

DDA通過(guò)選取不同動(dòng)力系數(shù)gg(gg=0～1)，即塊體當(dāng)前時(shí)步的初速度與上一時(shí)步末速度的比值，可同時(shí)處理靜力學(xué)與動(dòng)力學(xué)問(wèn)題。對(duì)于動(dòng)力學(xué)問(wèn)題，塊體運(yùn)動(dòng)過(guò)程中的慣性項(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)過(guò)程會(huì)產(chǎn)生離心力及科氏力，原始DDA采用一階近似后的線性表達(dá)式推導(dǎo)塊體的加速度表達(dá)式，而塊體加速度應(yīng)是從線性化前的位移增量表達(dá)式導(dǎo)出后再線性化，直接采用線性化后的位移增量表達(dá)式忽略了塊體旋轉(zhuǎn)時(shí)離心力和科氏力的作用。

3 算例驗(yàn)證

設(shè)計(jì)兩個(gè)典型算例來(lái)說(shuō)明原始 DDA采用一階近似后的位移增量簡(jiǎn)單疊加并推導(dǎo)塊體相應(yīng)的加速度表達(dá)式，及開(kāi)閉迭代后應(yīng)變?cè)隽康闹苯盈B加，在模擬塊體大轉(zhuǎn)動(dòng)時(shí)帶來(lái)的相關(guān)問(wèn)題。

3.1 正八邊形繞定點(diǎn)擺動(dòng)問(wèn)題

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

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

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

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

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

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

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

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

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

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

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

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

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

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

如圖8所示為以角速度ω=0.2rad/s 繞質(zhì)心勻速旋轉(zhuǎn)的等厚度矩形薄板，長(zhǎng)2L=6 m，寬2b=4 m。材料質(zhì)量密度為ρ=2800kg/m3，彈性模量為E=2 GPa。為方便起見(jiàn)，假定矩形薄板處于平面應(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

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

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

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

4 結(jié)論

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

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

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

(2) 應(yīng)變分量增量是相對(duì)當(dāng)前時(shí)刻定義的，不同時(shí)步應(yīng)變?cè)隽繀⒄諛?gòu)型不同，簡(jiǎn)單的累加導(dǎo)致應(yīng)變分量物理意義不明確；

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

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