，，

(長江科學(xué)院 水利部巖土力學(xué)與工程重點實驗室, 武漢 430010 )

1 研究背景

石根華博士提出的非連續(xù)變形方法(DDA)，是繼關(guān)鍵塊體理論后研究塊體系統(tǒng)的一種新的計算方法[1-2]。該方法通過滿足塊體系統(tǒng)力系的平衡、開-閉迭代及動力求解收斂，來獲得塊體系統(tǒng)力與變形真解。其自提出以來，一直作為巖土工程數(shù)值模擬方法研究中的前沿問題，為國內(nèi)外所重視[3-11]。

非連續(xù)變形分析方法是一種隱式求解的動力學(xué)計算方法，其中塊體接觸理論是DDA的核心問題。由于塊體系統(tǒng)的接觸是一個高度非線性問題，涉及參數(shù)眾多，如接觸面上的彈簧剛度、塊體的彈性模量、時間步長等。關(guān)于塊體系統(tǒng)接觸以及合適參數(shù)的取值等問題的研究，長期以來，主要借助DDA程序運行和具體問題分析，通過參數(shù)組合和試算等方式來獲取一些經(jīng)驗性認識。在某種條件下，建立一個解析的理論模型，來反映DDA計算過程中塊體接觸位移與上述參數(shù)間的關(guān)系，無疑是非常有意義的。本文試圖通過建立2塊體接觸分析模型，在一定的力學(xué)條件下，推導(dǎo)塊體間接觸位移的理論公式。根據(jù)接觸位移解析解，對塊體系統(tǒng)間塊體接觸開-閉迭代的過程以及參數(shù)間影響規(guī)律進行研究，以揭示DDA中塊體系統(tǒng)開-閉迭代收斂的本質(zhì)要求。

由于剛度矩陣組裝和聯(lián)立方程組求解較復(fù)雜，本文塊體接觸分析模型中，僅考慮2個規(guī)則的矩形塊體接觸，下部塊體固定，上部塊體受均勻分布的鉛直向載荷作用。因此，在該模型中，可只考慮法向嵌入這種情況，這樣僅涉及到法向剛度彈簧，整個計算過程可以手工完成，而且所得到的結(jié)果是解析的。建立了位移解析計算式后，便可以從量化的角度去研究時間步長、塊體彈性模量以及彈簧剛度等控制性參數(shù)之間的關(guān)系及其對計算結(jié)果的影響規(guī)律。根據(jù)本文的解析，證明了石根華博士在DDA計算中采用時間步長來調(diào)整和加快開-閉迭代收斂的策略是合理、高效的。該研究也揭示DDA各參數(shù)取值應(yīng)依據(jù)的基本規(guī)律，可以推動DDA的實際工程應(yīng)用。

2 DDA塊體接觸開-閉迭代分析模型

為定量研究DDA塊體接觸開-閉迭代的本質(zhì)問題，這里建立了如圖1所示的2塊體接觸開-閉迭代的塊體接觸分析模型。在該模型中，假定塊體N固定，塊體M的邊長為a,上部邊界受到均布荷載σ的作用，其變形模量為E。

圖1 塊體接觸開-閉迭代分析模型

應(yīng)用該模型推導(dǎo)各項剛度子矩陣和載荷子矩陣，從而形成DDA總體方程，來求解2個塊體系統(tǒng)的接觸位移。另外，從解析角度研究塊體接觸部位的變形與塊體系統(tǒng)各控制參數(shù)間的關(guān)系，包括塊體的變形模量E、時間步長Δ、彈簧剛度p等參數(shù)。

3 DDA塊體接觸位移理論公式推導(dǎo)

非連續(xù)變形分析(DDA)是以嚴格遵循經(jīng)典力學(xué)規(guī)則為基礎(chǔ)的。該方法是在假定塊體位移模式條件下，利用勢能最小原理建立塊體系統(tǒng)總體平衡方程。

3.1 二維DDA方法的一階位移形式

DDA方法做了如下假設(shè)：塊體的大位移和大變形是由塊體每步的小位移、小變形積累得到的。在每個時步內(nèi)，各塊體的位移足夠小，塊體的物理方程滿足彈性力學(xué)假定。在一階位移模式條件下，每個塊體內(nèi)應(yīng)力、應(yīng)變均為常應(yīng)力、常應(yīng)變。塊體的未知量列陣D由以下6個變形變量組成，即D=(u0v0γ0εxεyγxy)T。其中：(u0,v0)為塊體內(nèi)指定點(x0,y0)處的剛體位移；γ0為塊體運動時以點(x0,y0)為中心的轉(zhuǎn)動角(弧度)；(εx,εy,γxy)分別為塊體的正應(yīng)變和剪應(yīng)變。點(x0,y0)一般設(shè)為塊體的重心。

根據(jù)塊體運動學(xué)方程，塊體i內(nèi)任意點(x,y)的位移(u,v)可由位移轉(zhuǎn)置矩陣Ti和變形變量矩陣Di表示，即

(1)

式中：

Ti=

Di=[d1id2id3id4id5id6i]T=

[u0v0γ0εxεyγxy]T。

3.2 塊體系統(tǒng)子矩陣

由于塊體N為固定塊體，則塊體N變形變量的6個分量為0，即

D2=

(u2v2r2εx2εy2γxy2)T=(0 0 0 0 0 0)T

綜上，總體方程中只需要形成剛度系數(shù)矩陣K11、載荷矩陣F1、變形變量矩陣D1即可。而根據(jù)上述已知條件，將塊體的慣性子矩陣、彈性子矩陣、彈簧剛度子矩陣和載荷子矩陣總體集成，形成總體方程。其各子矩陣推導(dǎo)如下。

3.2.1 慣性力子矩陣

記(u(t)，v(t))為塊體i中任意點(x，y)的隨時間變化的位移，m為塊體單位面積質(zhì)量，則塊體單位面積慣性力為

(2)

塊體i的慣性力勢能為

(3)

假定D(tk)為第k時間步開始時位移，Δ為時間步長，Dk為該時步的位移增量，則根據(jù)時間積分可以求得

(4)

由此得

(5)

(6)

根據(jù)總勢能最小，對∏i求導(dǎo)，可得到相應(yīng)的慣性力子矩陣：

(7)

(8)

經(jīng)計算分析可得

?TiTTidxdy=

(9)

式中：s為塊體面積；s1,s2,s3為中間變量，s1=sxx-x0sx,s2=syy-y0sy,s3=sxy-x0sy;sx為對x軸的面積矩；sy為對y軸的面積矩；sxx為對x軸的面積慣性矩；syy為對y軸的面積慣性矩；sxy為對x和y軸的面積慣性矩。

對于塊體M有

s=a2,

這里只對靜力條件下的慣性子矩陣進行推導(dǎo)。即每個時步的開始速度為0，故V0=0，0表示零矩陣，公式(8)為“0”，該處0也表示零矩陣。故經(jīng)計算，只形成加到K11的慣性子矩陣，其慣性子矩陣為

(10)

3.2.2 彈性子矩陣

塊體i的彈性應(yīng)力所產(chǎn)生的應(yīng)變能為

(11)

根據(jù)總勢能最小，對∏e求導(dǎo)，可得到相應(yīng)的彈性子矩陣為

sEi→Kii。

(12)

對于塊體M，形成的彈性子矩陣為

(13)

3.2.3 彈簧子矩陣

由于彈簧子矩陣推導(dǎo)很復(fù)雜，涉及的變量很多，這里就沒有寫出詳細推導(dǎo)過程。針對該塊體接觸模型，直接引用文獻[1]中的公式：

(14)

(15)

(16)

第2個彈簧矩陣(x1=a，y1=0)為

(17)

3.2.4 載荷子矩陣

假設(shè)荷載是分布在從點(x1,y1)到點(x2,y2)的線段上，受荷線段的直線方程為

(18)

(19)

此處荷載沿作用線變化。線荷載(Fx(k),Fy(k))的勢能為

(20)

計算∏l的導(dǎo)數(shù)，即求勢能∏l的極小值：

r=1,…,6 。

(21)

其中

(22)

對于塊體M，F(xiàn)x=0，F(xiàn)y=-σ，(x1,y1)=(0,a)，(x2,y2)=(a,a)，l=a。則加到F1的載荷子矩陣為

=

(23)

3.3 接觸位移理論公式

對于靜力問題，由式(10)、式(13)、式(16)和式(17)形成剛度系數(shù)子矩陣K11，式(23)形成F1，由總方程KD=F求解變形變量。經(jīng)計算分析，最終得到3個方程：

經(jīng)計算求解，獲得接觸部位的變形變量D1=(0,v1,0,εx1,εy1,0)T。各相關(guān)分量計算公式為

(24)

表2 p=1E下的接觸位移

表3 p=10E下的接觸位移

由式(24)可知，2塊體的接觸位移與塊體的變形模量E、彈簧剛度p和時間步長Δ相關(guān)。在變形模量E和彈簧剛度p一定條件下，2塊體的接觸位移與時間步長直接有關(guān)。而且DDA中慣性力的引入直接與時間步長相關(guān)，因此，將慣性問題引入塊體接觸收斂求解的本質(zhì)問題是通過調(diào)整時間步長，實現(xiàn)塊體接觸開-閉迭代收斂。

4 DDA塊體接觸力及位移理論解分析

取塊體M邊長a=1 m,重度為2.7 T/m3，均布荷載σ=2.7 T/m2。由公式(24)計算2接觸彈簧處的變形變量。

4.1 接觸力計算

取E=10 GPa，p=10E,不同時間步長下,法向接觸力如表1所示。

表1 不同的時間步長下的法向接觸力

由于塊體M的外荷載σ=26.5 kN/m2，表1中法向接觸力與外荷載的相對誤差小于0.000 05%，則根據(jù)式(24)理論公式計算的法向接觸力可認為是與外力平衡的。由此，接觸部位滿足力的平衡，該公式的正確性得到了驗證。

4.2 接觸位移理論解計算及分析

分別取不同的塊體變形模量E、彈簧剛度p和時間步長Δ，由公式(1)計算塊體M接觸位移。經(jīng)計算分析，2彈簧接觸處的位移完全相同，現(xiàn)得到如表2至表4計算結(jié)果。

表4 p=100E下的接觸位移

分析表2至表4，當(dāng)塊體變形模量E、彈簧剛度系數(shù)一定時，時間步長Δ越小，接觸部位的侵入位移越??；當(dāng)時間步長Δ、彈簧剛度系數(shù)一定時，塊體變形模量E越大，接觸部位的侵入位移越?。划?dāng)塊體變形模量E、時間步長Δ一定時，彈簧剛度系數(shù)越大，接觸部位的侵入位移越小。

取E=10 GPa，p=10E, 接觸部位的侵入位移隨時間步長的關(guān)系曲線如圖2所示。

圖2 接觸位移與時間步長的關(guān)系曲線

圖2表明，在變形模量和彈簧剛度一定條件下，時間步長越小，接觸部位的侵入位移越小。為了保證塊體系統(tǒng)的開-閉迭代收斂，應(yīng)合理地調(diào)整時間步長。

5 結(jié) 論

非連續(xù)變形分析是一種隱式求解的動力學(xué)計算方法，且采用在塊體界面加減剛硬彈簧的方式來滿足塊體界面無張拉和無嵌入的接觸準則。本文基于塊體接觸分析模型，推導(dǎo)了靜力條件下塊體接觸位移理論公式，并通過接觸力的計算驗證了該公式的正確性。進而探討了塊體變形模量、彈簧剛度和時間步長與塊體接觸位移的相互關(guān)系，在變形模量和彈簧剛度一定條件下，2塊體在接觸部位的位移將直接與時間步長有關(guān)，時間步長越小，接觸部位的侵入位移越小。該結(jié)論也證明了石根華博士在DDA計算中采用時間步長來調(diào)整和加快開-閉迭代收斂的策略是合理、高效的。而且，DDA中慣性力的引入直接與時間步相關(guān)，因此，將慣性問題引入塊體接觸收斂求解的本質(zhì)問題是通過調(diào)整時間步長，實現(xiàn)塊體接觸開-閉迭代收斂。

參考文獻：

[1] SHI Gen-hua. Discontinuous Deformation Analysis: A New Numerical Model for the Statics and Dynamics of Block System[D]. Berkeley: Department of Civil Engineering, University of California, 1988.

[2] 石根華. 數(shù)值流形方法與非連續(xù)變形分析[M]. 裴覺民，譯. 北京：清華大學(xué)出版社, 1997. (SHI Gen-hua. Numerical Manifold Method and Discontinuous Deformation Analysis [M]. Translated by PEI Jue-min. Beijing: Tsinghua University Press, 1997.(in Chinese))

[3] MACLAUGHLIN M M, BERGER E A. A Decade of DDA Validation[C]∥Proceedings of the 6th International Conference on Analysis of Discontinuous Deformation, Norway, October 5-8, 2003:13-32.

[4] OHNISHI Y, NISHIYAMA S, SAKAKI T,etal. The Application of DDA to Practical Rock Engineering Problems: Issues and Recent insights[C]∥Proceedings of the 7th International Conference on Analysis of Discontinuous Deformation, Hawaii, December 10-12, 2005: 277-287.

[5] WANG C Y, SHENG J, CHEN M H. Time Integration Theories for the DDA Method with Finite Element Meshes[C]∥Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, Albuquerque: TSI Press, June 12-14, 1996: 263-287.

[6] 鄔愛清，丁秀麗，陳勝宏,等. DDA方法在復(fù)雜地質(zhì)條件下地下廠房圍巖變形與破壞特征分析中的應(yīng)用研究[J]. 巖石力學(xué)與工程學(xué)報，2006,25(1):1-8. (WU Ai-qing, DING Xiu-li, CHEN Sheng-hong,etal. Researches on Deformation and Failure Characteristics of an Underground Powerhouse with Complicated Geological Conditions by DDA Method[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(1):1-8.(in Chinese))

[7] MACLAUGHLIN M, SITAR N, DOOLIN D,etal. Investigation of Slope-stability Kinematics Using Discontinuous Deformation Analysis[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(5): 753-762.

[8] HATZOR Y H, FEINTUCH A. The Validity of Dynamic Block Displacement Prediction Using DDA[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38: 599-606.

[9] 鄔愛清. DDA方法塊體穩(wěn)定性驗證及其在巖質(zhì)邊坡穩(wěn)定性分析中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報，2008,27(4): 664-672. (WU Ai-qing. Validation for Rock Block Stability and Its Application to Rock Slope Stability Evaluation Using DDA Method[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(4): 664-672.(in Chinese))

[10] 邱寬紅，林紹忠，黃 斌. 基于DDA的膨脹土邊坡破壞過程模擬[J]. 長江科學(xué)院院報，2009,26(11)：58-61. (QIU Kuan-hong, LIN Shao-zhong, HUANG Bin. Failure Simulation of Expansive Soil Slope Based on DDA [J]. Journal of Yangtze River Scientific Research Institute, 2009, 26(11): 58-61. (in Chinese))

[11] 劉 軍，李仲奎. 非連續(xù)變形分析方法中一些控制參數(shù)的設(shè)置[J]. 成都理工大學(xué)學(xué)報(自然科學(xué)版)，2004, 31(5):522-526.(LIU Jun, LI Zhong-kui. Setting Some Control Parameters in the Discontinuous Deformation Analysis Method[J]. Journal of Chengdu University of Technology (Science & Technology Edition), 2004, 31(5): 522-526.(in Chinese))