，，

(長(zhǎng)江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室, 武漢 430010 )

1 研究背景

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

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

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

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

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

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

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

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

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

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

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

根據(jù)塊體運(yùn)動(dòng)學(xué)方程，塊體i內(nèi)任意點(diǎn)(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個(gè)分量為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中任意點(diǎn)(x，y)的隨時(shí)間變化的位移，m為塊體單位面積質(zhì)量，則塊體單位面積慣性力為

(2)

塊體i的慣性力勢(shì)能為

(3)

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

(4)

由此得

(5)

(6)

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

(7)

(8)

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

?TiTTidxdy=

(9)

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

對(duì)于塊體M有

s=a2,

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

(10)

3.2.2 彈性子矩陣

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

(11)

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

sEi→Kii。

(12)

對(duì)于塊體M，形成的彈性子矩陣為

(13)

3.2.3 彈簧子矩陣

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

(14)

(15)

(16)

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

(17)

3.2.4 載荷子矩陣

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

(18)

(19)

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

(20)

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

r=1,…,6 。

(21)

其中

(22)

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

=

(23)

3.3 接觸位移理論公式

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

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

(24)

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

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

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

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

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

4.1 接觸力計(jì)算

取E=10 GPa，p=10E,不同時(shí)間步長(zhǎng)下,法向接觸力如表1所示。

表1 不同的時(shí)間步長(zhǎng)下的法向接觸力

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

4.2 接觸位移理論解計(jì)算及分析

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

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

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

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

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

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

5 結(jié) 論

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

參考文獻(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é)報(bào)，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)定性驗(yàn)證及其在巖質(zhì)邊坡穩(wěn)定性分析中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào)，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的膨脹土邊坡破壞過(guò)程模擬[J]. 長(zhǎng)江科學(xué)院院報(bào)，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é)報(bào)(自然科學(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))