基于CEM的加錨節(jié)理巖體動(dòng)力學(xué)算法研究*

23

(1.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究有限責(zé)任公司，湖北 武漢 430010；2.武漢大學(xué)土木建筑工程學(xué)院，湖北 武漢 430072；3.陜西省混凝土結(jié)構(gòu)安全與耐久性重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710123；4.長(zhǎng)江巖土工程總公司(武漢)，湖北 武漢 430010)

大量巖土邊坡地震崩塌與滑動(dòng)主要發(fā)生在多震的西部地區(qū)，這些地區(qū)地形地質(zhì)條件非常復(fù)雜，巖體中存在斷層、節(jié)理、裂隙等各種結(jié)構(gòu)面。這些地區(qū)復(fù)雜的地質(zhì)構(gòu)造和較高的地震烈度使得錨固巖體的抗震穩(wěn)定也成為工程的關(guān)鍵問題之一。

數(shù)值計(jì)算方法[1]能夠體現(xiàn)出復(fù)雜巖體的力學(xué)特性，已被成功地應(yīng)用于巖質(zhì)邊坡的靜力穩(wěn)定性分析中，在巖質(zhì)邊坡的動(dòng)力計(jì)算研究方面也取得了較大的進(jìn)展，眾多學(xué)者采用有限單元法(FEM)[2]、非連續(xù)變形分析法(DDA)[3]、離散單元法[4]等多種計(jì)算方法對(duì)節(jié)理巖體進(jìn)行動(dòng)力學(xué)計(jì)算分析，但在模擬加錨節(jié)理巖體方面大多處于初步發(fā)展階段，考慮的形式比較簡(jiǎn)單，還需諸多努力。目前仍然缺少一種既能同時(shí)實(shí)現(xiàn)對(duì)結(jié)構(gòu)面、加固措施的精細(xì)模擬，又能滿足動(dòng)態(tài)設(shè)計(jì)對(duì)于建模效率要求的數(shù)值模擬方法[5]。

為解決這個(gè)難題，相關(guān)專家提出了復(fù)合單元法(CEM)，該方法主要有以下三個(gè)優(yōu)點(diǎn)：?準(zhǔn)確：能夠?qū)崿F(xiàn)結(jié)構(gòu)面和加固措施的離散模擬，具有細(xì)致精確的優(yōu)點(diǎn)；?快速：在建模時(shí)無(wú)須考慮結(jié)構(gòu)面和加固措施，通過(guò)引入復(fù)合單元，用復(fù)雜的計(jì)算問題來(lái)代替復(fù)雜的結(jié)構(gòu)面和加固措施建模問題，拓?fù)溆?jì)算得到整合后的復(fù)合單元均由程序自動(dòng)實(shí)現(xiàn)；?動(dòng)態(tài)設(shè)計(jì)：當(dāng)設(shè)計(jì)方案優(yōu)化調(diào)整時(shí)，不需要重新建模，根據(jù)新的支護(hù)措施，程序重新拓?fù)溆?jì)算得到新的復(fù)合單元即可。

目前為止，靜力學(xué)計(jì)算方面，復(fù)合單元法已經(jīng)實(shí)現(xiàn)了對(duì)結(jié)構(gòu)面(如節(jié)理、斷層等)[5]、多物理場(chǎng)模擬(如溫度、滲流)等和加固措施(如錨桿、錨索等)的模擬，但是，在動(dòng)力學(xué)計(jì)算方面，僅實(shí)現(xiàn)了節(jié)理巖體模擬[6]，還未實(shí)現(xiàn)對(duì)加錨節(jié)理巖體動(dòng)力學(xué)計(jì)算。

本文進(jìn)一步發(fā)展了復(fù)合單元法，建立了加錨節(jié)理巖體復(fù)合單元法動(dòng)力學(xué)算法，并首次將粘彈性邊界在復(fù)合單元上進(jìn)行了初步應(yīng)用，通過(guò)算例證明了算法的合理性，為加錨節(jié)理巖體的動(dòng)力學(xué)計(jì)算研究提供了新的思路。

1 復(fù)合單元法動(dòng)力學(xué)算法步驟

復(fù)合單元法動(dòng)力學(xué)模型的建立和有限元?jiǎng)恿Ψ治龇椒ㄊ穷愃频?，通過(guò)位移空間離散化，可以最終得到系統(tǒng)求解方程為：

(1)

M——質(zhì)量矩陣；

C——阻尼矩陣；

K——?jiǎng)偠染仃嚕?/p>

P——節(jié)點(diǎn)載荷向量。

在得到上述系統(tǒng)運(yùn)動(dòng)方程之后，可以選擇用直接積分法或者模態(tài)分解法來(lái)求解系統(tǒng)響應(yīng)。

復(fù)合單元法動(dòng)力學(xué)方法可以分為以下幾步：

a.生成復(fù)合單元網(wǎng)格信息，一個(gè)復(fù)合單元可能會(huì)包括若干子單元和接觸面。

d.根據(jù)b和c生成的復(fù)合單元?jiǎng)偠染仃嚭蛷?fù)合單元質(zhì)量矩陣進(jìn)行模態(tài)分析，根據(jù)其自振頻率和阻尼比組裝形成Rayleigh阻尼矩陣Ce。

e.將復(fù)合單元區(qū)域的積分結(jié)果和有限元積分結(jié)果求和，即組裝系統(tǒng)整體剛度矩陣K=∑Ke、整體質(zhì)量矩陣M=∑Me、整體阻尼矩陣C=∑Ce和整體節(jié)點(diǎn)載荷向量P=∑Pe。

f.求解系統(tǒng)方程。

2 復(fù)合單元質(zhì)量矩陣和積分方法

對(duì)復(fù)合單元采用一致質(zhì)量矩陣為：

(2)

在復(fù)合單元法中，對(duì)于含各種不連續(xù)面的復(fù)合單元，對(duì)于子單元塊體積分方法是先將組成該塊體的頂點(diǎn)按照一定順序和規(guī)則連接形成若干四面體，再將每個(gè)四面體用行波法自動(dòng)剖分網(wǎng)格的方式分成4個(gè)六面體，具體方法見圖1。

圖1 子單元塊體積分方法示意

在復(fù)合單元法中，對(duì)于錨桿桿體和砂漿子單元塊體積分方法為是建立局部笛卡兒坐標(biāo)系和局部柱坐標(biāo)系，對(duì)其進(jìn)行積分，坐標(biāo)系見圖2。

圖2 局部笛卡爾坐標(biāo)系和錨桿局部柱坐標(biāo)系

3 復(fù)合單元質(zhì)量剛度矩陣

加錨節(jié)理巖體的復(fù)合單元示意見圖3。

圖3 加錨節(jié)理巖體的復(fù)合單元示意

如圖3所示，在復(fù)合單元內(nèi)，分別定義不同材料區(qū)域?yàn)樽訂卧?，在?fù)合單元引入5套虛節(jié)點(diǎn)加上實(shí)節(jié)點(diǎn)共6套節(jié)點(diǎn)，用于描述兩個(gè)塊體中的巖石、錨桿、砂漿等子單元，相應(yīng)地定義6套位移分量分別用來(lái)插值6個(gè)子單元內(nèi)部的位移。

根據(jù)虛功原理，可得

(3)

經(jīng)過(guò)推導(dǎo)得到以下平衡方程組：

(4)

(5)

4 復(fù)合單元阻尼矩陣

阻尼矩陣采用Rayleigh阻尼假定，其阻尼矩陣是質(zhì)量矩陣和剛度矩陣的組合：

[C]=α[M]+β[K]

(6)

式中α——Rayleigh阻尼的質(zhì)量系數(shù)；

β——Rayleigh阻尼的剛度系數(shù)。

將上述復(fù)合單元的矩陣[M]、[C]和[K]代入式(1)，然后采用Wilson-θ法進(jìn)行每一個(gè)時(shí)間步的求解。

5 復(fù)合單元黏彈性人工邊界

進(jìn)行邊坡動(dòng)力分析需要考慮邊坡與周圍基巖的動(dòng)力相互作用，無(wú)限地基如何模擬是常常要面對(duì)的問題。黏彈性人工邊界是一種十分理想的動(dòng)力邊界[7]，基本思想是在人工邊界節(jié)點(diǎn)的法向和切向設(shè)置并聯(lián)的彈簧和阻尼器來(lái)實(shí)現(xiàn)。

三維人工邊界上法向和切向的彈簧剛度和阻尼系數(shù)取值如下：

(7)

(8)

式中kbN——單位面積的法向彈簧剛度；

kbT——單位面積的切向彈簧剛度；

CbN——單位面積的法向阻尼器的阻尼系數(shù)；

CbT——單位面積的切向阻尼器的阻尼系數(shù)；

r——波源至人工邊界點(diǎn)的距離；

λ、G——介質(zhì)Lame常數(shù)；

ρ——介質(zhì)質(zhì)量密度；

cP——介質(zhì)的P波波速；

cS——介質(zhì)的S波波速；

A、B——參數(shù)，取建議值[7]為A=0.8，B=1.1。

當(dāng)黏彈性邊界單元含有復(fù)合單元時(shí)，以復(fù)合單元內(nèi)子單元為基本單位，依次得到每個(gè)子單元對(duì)應(yīng)的虛實(shí)節(jié)點(diǎn)的法向與切向彈簧剛度和法向與切向阻尼器的阻尼系數(shù)，取值為單位面積的法向與切向彈簧剛度和單位面積的法向與切向阻尼器的阻尼系數(shù)均分別乘以其有效面積。

6 技術(shù)路線

技術(shù)路線見圖4。

圖4 技術(shù)路線

7 算例考證

7.1 算例一

取一個(gè)尺寸為1m×1m×10m(長(zhǎng)×寬×高)的巖體，巖體中心內(nèi)含有一條長(zhǎng)度為10m砂漿錨桿，鋼筋的半徑16mm，砂漿的半徑80mm，5m高處含有一個(gè)1m×1m水平節(jié)理。復(fù)合單元法計(jì)算網(wǎng)格見圖5，單元數(shù)為480，節(jié)點(diǎn)數(shù)為820。模型底部全約束。表1為實(shí)體材料力學(xué)參數(shù)。表2為接觸面力學(xué)參數(shù)。模型Y方向施加線荷載F(t)=100kN/m(0≤t≤1)，作用位置見圖5。采用一階振型的瑞利阻尼，分別考慮無(wú)阻尼、阻尼比0.05和阻尼比0.01的計(jì)算工況。

圖5 復(fù)合單元法網(wǎng)格和荷載示意

表1 彈性材料參數(shù)

表2 接觸面材料參數(shù)

7.2 算例一結(jié)果分析

圖6～圖8是以加錨節(jié)理巖體為對(duì)象不同阻尼比下線荷載作用中點(diǎn)位移歷時(shí)曲線圖，算例對(duì)比了復(fù)合單元法和有限單元法的計(jì)算結(jié)果，從圖中可以看出，隨著阻尼系數(shù)的增大，達(dá)到平衡位置的時(shí)間變短，并且，復(fù)合單元法和有限單元法這兩種方法計(jì)算的結(jié)果變化規(guī)律相似，誤差較小，量值也較接近，說(shuō)明復(fù)合單元法的計(jì)算結(jié)果是合理的。

圖6 無(wú)阻尼時(shí)線荷載作用中點(diǎn)位移歷時(shí)曲線

圖7 阻尼比0.01線荷載作用中點(diǎn)位移歷時(shí)曲線

圖8 阻尼比0.05線荷載作用中點(diǎn)位移歷時(shí)曲線

7.3 算例二

從半無(wú)限空間截取一個(gè)尺寸為4m×1m×3m(長(zhǎng)×寬×高)的彈性材料作為有限計(jì)算區(qū)域，3m高處還有一水平節(jié)理，頂端自由，法向?yàn)閅的面法向約束，法向?yàn)閄和Z的面施加粘彈性邊界。材料彈性模量為2.5Pa，泊松比為0.25，材料密度為1kg/m3，節(jié)理法向剛度為10Pa/m，切向剛度為5Pa/m。時(shí)間步長(zhǎng)為0.01s，總計(jì)算時(shí)間為15s。動(dòng)荷載為作用在頂面的壓力荷載F(t)，數(shù)值大小見圖9。網(wǎng)格模型和測(cè)點(diǎn)見圖10，單元數(shù)為300，節(jié)點(diǎn)數(shù)為672。

圖9 自由表面的壓力荷載F(t)歷時(shí)曲線

圖10 網(wǎng)格模型和測(cè)點(diǎn)示意

為了研究半無(wú)限空間粘彈性邊界計(jì)算結(jié)果的合理性，增加兩組計(jì)算工況。工況一：采用固定邊界，保持模型和荷載不變，將粘彈性邊界改為法向約束；工況二：在工況一的基礎(chǔ)上，將模型尺寸變?yōu)?00m×1m×300m(長(zhǎng)×寬×高)，頂面自由，其余為法向約束，此工況的結(jié)果可以作為準(zhǔn)理論解。

7.4 算例二結(jié)果分析

圖11和圖12為測(cè)點(diǎn)1和測(cè)點(diǎn)2不同工況下的Z向位移歷時(shí)曲線，從計(jì)算結(jié)果可以看出，粘彈性人工邊界能夠很有效地吸收散射波在邊界處產(chǎn)生的能量，消除邊界的反射，與準(zhǔn)理論解從量值和規(guī)律上都比較吻合。固定邊界由于不能吸收邊界處散射波的反射，夸大了能量的輸入，造成了振蕩現(xiàn)象。由此，也可以證明粘彈性邊界在復(fù)合單元法中的應(yīng)用是合理的。

圖11 測(cè)點(diǎn)1不同工況下的Z向位移歷時(shí)曲線

圖12 測(cè)點(diǎn)2不同工況下的Z向位移歷時(shí)曲線

8 結(jié)論與展望

本文開展了加錨節(jié)理巖體復(fù)合單元法動(dòng)力學(xué)研究，推導(dǎo)了加錨節(jié)理巖體復(fù)合單元法質(zhì)量矩陣、剛度矩陣和阻尼矩陣，建立了復(fù)合單元粘彈性人工邊界，建立了慣性力、阻尼力、動(dòng)力荷載及彈性力作用下的復(fù)合單元?jiǎng)恿ο到y(tǒng)的控制方程組，編制相應(yīng)程序，并進(jìn)行算例驗(yàn)證。

a.通過(guò)算例結(jié)果的對(duì)比可知，復(fù)合單元法的計(jì)算精度與有限單元法基本一致，粘彈性邊界在復(fù)合單元法中的應(yīng)用是合理的，驗(yàn)證了復(fù)合單元法能有效地模擬加錨節(jié)理巖體，進(jìn)行動(dòng)力學(xué)計(jì)算。

b.提供了一套完整的錨固巖體復(fù)合單元法動(dòng)力學(xué)分析的技術(shù)路線，為加錨節(jié)理巖體動(dòng)力學(xué)計(jì)算提供了一套高效的數(shù)值模擬方法。

c.本文建立的加錨節(jié)理巖體復(fù)合單元法動(dòng)力學(xué)模型僅是初步研究成果，只能處理一些簡(jiǎn)單的算例，僅考慮了彈性狀態(tài)，后續(xù)研究工作需要合理的錨固巖體動(dòng)力學(xué)本構(gòu)模型作為技術(shù)支撐。

d.復(fù)合單元法進(jìn)行動(dòng)力分析時(shí)，復(fù)合單元粘彈性邊界條件下地震波的輸入，垂直從地基底部入射或者考慮波動(dòng)斜入射等問題，都值得未來(lái)更加深入地研究。