有限元與離散元混合法在裂紋擴(kuò)展中的應(yīng)用

孫 翔,劉傳奇,薛世峰

(中國石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院,山東青島 266580)

有限元與離散元混合法在裂紋擴(kuò)展中的應(yīng)用

孫 翔,劉傳奇,薛世峰

(中國石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院,山東青島 266580)

基于有限元與離散元混合方法研究裂紋擴(kuò)展模擬問題。對(duì)含原生裂紋的結(jié)構(gòu)進(jìn)行單元離散,用有限元計(jì)算單元內(nèi)部,采用離散元計(jì)算單元界面,通過單元連接形式的轉(zhuǎn)變實(shí)現(xiàn)連續(xù)到非連續(xù)的轉(zhuǎn)化；采用平面半彈簧法進(jìn)行接觸判斷,通過顯示迭代求解運(yùn)動(dòng)方程,不斷更新單元坐標(biāo),實(shí)現(xiàn)裂紋擴(kuò)展的數(shù)值模擬。以單裂紋與雁形裂紋在單向位移載荷作用下的擴(kuò)展為例,對(duì)比數(shù)模結(jié)果與試驗(yàn)結(jié)果。結(jié)果表明:采用有限元與離散元混合方法可有效模擬單、多裂紋的擴(kuò)展過程；巖橋?yàn)?0°的雁形裂紋受對(duì)向擠壓載荷作用發(fā)生翼裂-翼裂貫通。

有限單元法；離散單元法；裂紋擴(kuò)展；數(shù)值模擬

擴(kuò)展過程的本質(zhì)是連續(xù)面到非連續(xù)面的轉(zhuǎn)化,如何對(duì)不連續(xù)界面準(zhǔn)確刻畫是問題的關(guān)鍵。近幾十年,裂紋擴(kuò)展研究發(fā)展迅速[1-2]。經(jīng)典有限單元法通過尖端設(shè)置奇異單元,可用來模擬裂紋擴(kuò)展。困難在于裂紋擴(kuò)展后須進(jìn)行網(wǎng)格重構(gòu)[3]。無網(wǎng)格方法可高精度地模擬裂紋擴(kuò)展,但邊界條件以及數(shù)值積分難以確定[4]。黏聚力模型通過設(shè)置界面單元可實(shí)現(xiàn)結(jié)構(gòu)的斷裂,但須預(yù)先估計(jì)擴(kuò)展路徑[5]。擴(kuò)展有限元與廣義有限元是基于單位分解框架的單位分解法,難點(diǎn)在于確定復(fù)雜裂紋的附加函數(shù),由于沒有接觸算法,亦無法考慮開裂紋的閉合影響[6]。離散單元法作為一種非連續(xù)方法,在巖土工程領(lǐng)域得到廣泛應(yīng)用[7]。離散單元法中求解單元內(nèi)部變形常采用高斯定理和Wilkins的大變形有限差分法[8]。近年來,由Munjiza[9]提出的離散元與有限元混合方法發(fā)展迅速,李世海等[10]提出基于連續(xù)介質(zhì)離散元法,其本質(zhì)亦為有限元與離散元的混合方法?，F(xiàn)有文獻(xiàn)中尚未有人將離散元與有限元混合方法應(yīng)用于裂紋擴(kuò)展的數(shù)值模擬中。筆者編制二維程序,模擬單裂紋、多裂紋在壓剪作用下的擴(kuò)展情況,通過與試驗(yàn)結(jié)果的對(duì)比,驗(yàn)證此方法模擬裂紋擴(kuò)展時(shí)的有效性。

1 基本原理

采用有限元與離散元混合方法模擬裂紋擴(kuò)展的基本過程可概括為:采用網(wǎng)格剖分對(duì)含原生裂紋的結(jié)構(gòu)進(jìn)行單元離散,采用有限元計(jì)算單元內(nèi)部,采用離散元計(jì)算單元邊界,以修正的Mohr-Coulomb準(zhǔn)則作為斷裂判據(jù),通過彈簧連結(jié)的斷裂產(chǎn)生單元分離,采用半彈簧法進(jìn)行接觸處理,通過顯示迭代求解運(yùn)動(dòng)方程,不斷更新單元坐標(biāo),實(shí)現(xiàn)裂紋擴(kuò)展的數(shù)值模擬。

1.1 基本方程

選用三角形常應(yīng)變單元,單元與單元之間通過法向與切向彈簧相聯(lián)(圖1),計(jì)算過程須滿足運(yùn)動(dòng)方程與本構(gòu)方程。其中,運(yùn)動(dòng)方程為核心方程,在單元內(nèi)部滿足本構(gòu)方程,彈簧僅作為數(shù)學(xué)處理,通過增大其剛度以保證結(jié)構(gòu)變形主要發(fā)生在單元內(nèi)部,而非單元界面上,采用類似于有限元法處理接觸問題時(shí)的罰函數(shù)法。

圖1 單元間聯(lián)接關(guān)系Fig.1 Connections between elements

圖1中,i單元1號(hào)節(jié)點(diǎn)的運(yùn)動(dòng)方程為

式中,mi1為i單元1號(hào)節(jié)點(diǎn)的集中質(zhì)量(在三角形常應(yīng)變單元中,mi1=ρA/3)；ρ為材料密度；A為單元面積；Fd為單元變形力；Fj為界面作用力；Fc為阻尼力；Fe為外力。

(1)單元內(nèi)部。單元為彈性,采用有限元計(jì)算變形力,由于顯式求解,無須形成總剛,僅需單元?jiǎng)偠?便可由節(jié)點(diǎn)位移求得節(jié)點(diǎn)力?；舅悸啡缦?

由最小勢能原理可推得單元的剛度矩陣為

式中,[B]為幾何矩陣；[D]為彈性矩陣。

當(dāng)單元節(jié)點(diǎn)位移已知時(shí),由

可求得節(jié)點(diǎn)力,即為單元變形力。單元內(nèi)部應(yīng)變、應(yīng)力可分別求得為

(2)界面作用。將單元界面的相互作用分為連接型和接觸型,單元界面作用力即可分為彈簧作用力與單元接觸力。單元破壞之前均為連接型,提供彈簧作用力,此時(shí)有

式中,F為彈簧作用力；K為彈簧剛度；Δd為彈簧位移；上標(biāo)j表示第j根彈簧,下標(biāo)n,s分別表示法向與切向。

進(jìn)行破壞計(jì)算時(shí),采用帶拉伸破壞的修正Mohr-Coulomb準(zhǔn)則[11]對(duì)彈簧力進(jìn)行修正。

式中,T為抗拉強(qiáng)度；φ為內(nèi)摩擦角；C為黏聚力。

彈簧無論發(fā)生何種模式的斷裂,連接型均變?yōu)榻佑|型,此時(shí)界面力為接觸力。單元不發(fā)生接觸時(shí),相互作用力為零,單元相互擠壓時(shí),作用力如下:

式中,Δdjn＜0表示單元相互擠壓。

設(shè)定彈簧有兩個(gè)作用:①求解連續(xù)過程中起到連接單元的作用；②通過彈簧的斷裂產(chǎn)生破壞,起到黏聚力模型[5]中界面單元的作用。

(3)阻尼力與外力。阻尼采用局部無黏性阻尼[12],節(jié)點(diǎn)上的阻尼力與非平衡力成比例,阻尼力方向確保能量總是能夠消散,阻尼力表示為

1.2 求解方法

運(yùn)動(dòng)方程的求解分為顯式解法和隱式解法,采用動(dòng)態(tài)松弛方法求解。動(dòng)態(tài)松弛法是將線性或非線性靜力問題轉(zhuǎn)化為動(dòng)力問題,通過設(shè)置阻尼消耗能量,最終使位移收斂到平衡位置(即勢能最小位置)的顯式求解方法[13]。

設(shè)時(shí)間步長為Δt,采用中心差分格式,節(jié)點(diǎn)速度與節(jié)點(diǎn)新坐標(biāo)為

在每個(gè)計(jì)算時(shí)步后,需更新節(jié)點(diǎn)坐標(biāo)。整個(gè)計(jì)算流程如圖2所示。

圖2 計(jì)算流程Fig.2 Flow chart of calculation process

2 關(guān)鍵問題的處理

按照Cundall的理論[14],離散元程序必須能夠自動(dòng)識(shí)別接觸,接觸處理是離散元程序中最重要的部分。此外,選擇時(shí)步和加速計(jì)算等問題是編程的關(guān)鍵。

2.1 接觸判斷

目前,離散元程序?qū)佑|的判斷普遍采用公共面法(common plane method),但此方法需進(jìn)行接觸類型(點(diǎn)點(diǎn),點(diǎn)線,線線)的判斷。將半彈簧法[15]簡化至二維情況,以此進(jìn)行接觸判斷。如圖3所示,將判斷單元的各邊頂點(diǎn)向中心縮進(jìn)小段距離(如5%),生成半彈簧,當(dāng)半彈簧與目標(biāo)單元的距離在一定范圍內(nèi)時(shí),生成接觸。同樣過程應(yīng)用到各單元中,通過對(duì)若干特殊情況的處理(如半彈簧生成點(diǎn)又是目標(biāo)單元的插值點(diǎn)時(shí)須進(jìn)行接觸力減半處理),接觸力可直接根據(jù)半彈簧與目標(biāo)面間的相對(duì)位移確定,而無須了解母塊體與子塊體間的具體接觸類型。

圖3 接觸判斷模型Fig.3 Contact judgement model

2.2 時(shí)步選擇

離散元的求解是條件穩(wěn)定的,時(shí)間步長的選擇應(yīng)以計(jì)算收斂為前提,提高收斂速度為目的。時(shí)間步長的設(shè)定須滿足單元內(nèi)部計(jì)算與單元間相對(duì)位移的計(jì)算均收斂。本例參照離散元軟件UDEC,設(shè)定步長。

單元內(nèi)部計(jì)算收斂的時(shí)間步長為

式中,mi為單元節(jié)點(diǎn)i的集中質(zhì)量；ki為節(jié)點(diǎn)周圍的彈簧剛度。

單元間相對(duì)位移計(jì)算收斂確定的時(shí)間步長為

式中,Mmin為最小單元質(zhì)量；Kmax為最大接觸剛度； frac為修正系數(shù),通常取為0.1。

時(shí)間步長設(shè)定為

2.3 加速處理

程序設(shè)計(jì)中,很大部分機(jī)時(shí)消耗在判斷接觸上。接觸檢測算法[9]主要有:基于二叉樹結(jié)構(gòu)的接觸檢測算法；直接映射算法；篩選接觸檢測算法；排序算法以及Munjiza提出的非二叉搜索算法。

針對(duì)裂紋擴(kuò)展模擬中,單元并不完全分離這一特點(diǎn),僅須將各單元共節(jié)點(diǎn)的單元設(shè)定為其潛在接觸單元,若研究開裂紋受壓破壞過程,還須將含裂紋壁的所有單元互相定義為潛在接觸單元。此外,僅當(dāng)彈簧發(fā)生斷裂時(shí),在斷裂彈簧兩節(jié)點(diǎn)相應(yīng)位置,生成半彈簧,而其余彈簧仍為連接型。采用上述加速處理,計(jì)算受壓條件下,10 000個(gè)單元的結(jié)構(gòu)破壞,使用Pentium E5800處理器,僅須計(jì)算30 min。

3 裂紋擴(kuò)展模擬

分別對(duì)初始裂紋為傾角為45°的單裂紋以及巖橋?yàn)?0°的雁形裂紋在單向位移載荷作用下,裂紋的擴(kuò)展過程進(jìn)行模擬。

3.1 單裂紋擴(kuò)展模擬

圖4為計(jì)算模型示意圖(單位:cm),位移載荷按0.001 cm/s施加。

圖4 單裂紋計(jì)算模型Fig.4 Single crack computation model

巖體彈性模量E為20 GPa；摩擦系數(shù)μ為0. 2；材料的內(nèi)聚力C為100 kPa；材料抗拉強(qiáng)度T為20 kPa。彈簧剛度Kn、Ks均為100 GN/m。

若界面單元的剛度過小,則材料表現(xiàn)軟弱；若界面剛度充分大,但會(huì)顯著減小顯式計(jì)算的穩(wěn)定時(shí)間步長,降低計(jì)算效率[16],通過試算,將彈簧剛度設(shè)為巖石彈性模量的50倍,以減小彈簧帶來的誤差。

模擬結(jié)果顯示,計(jì)算前2000步,裂紋并未發(fā)生擴(kuò)展,僅產(chǎn)生了較明顯的應(yīng)力集中,原生裂紋寬度有所減小,但未發(fā)生裂紋壁的接觸。繼續(xù)加載,當(dāng)計(jì)算到4000步時(shí),可較清楚地觀察到由于拉伸破壞產(chǎn)生的翼型裂紋,裂紋進(jìn)入曲線擴(kuò)展過程,次生裂紋不斷向最大主應(yīng)力方向趨近。在裂紋寬度方面,越靠近原生裂紋尖端,擴(kuò)展裂紋的寬度越大。圖5為不同時(shí)步下的結(jié)構(gòu)水平向位移云圖,圖6為數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果[17]的對(duì)比?？梢钥闯?采用有限元與離散法混合方法模擬裂紋擴(kuò)展是十分有效的。

圖5 單裂紋擴(kuò)展水平位移Fig.5 Horizontal displacements of single crack propagation

圖6 數(shù)模與試驗(yàn)結(jié)果對(duì)比Fig.6 Contrast of simulation and experiment results

3.2 多裂紋擴(kuò)展模擬

研究巖橋傾角為90°的雙裂紋在按0.001 cm/s施加的對(duì)向擠壓位移載荷作用下的裂紋擴(kuò)展情況。圖7為結(jié)構(gòu)尺寸(單位:cm),且將內(nèi)端定義為巖橋所在區(qū)域,外端定義為另外兩個(gè)裂紋尖端影響區(qū)域。

圖7 多裂紋計(jì)算模型Fig.7 Multi-crack computation model

類似于單一裂紋,兩裂紋尖端起初均發(fā)生較明顯的應(yīng)力集中,隨著載荷的繼續(xù)施加,首先出現(xiàn)I型翼型裂紋,當(dāng)計(jì)算6000～8000步時(shí),內(nèi)端次生裂紋迅速發(fā)展以致貫穿,外端次生裂紋發(fā)育相對(duì)緩慢,均發(fā)生拉伸破壞。圖8為不同時(shí)步下,多裂紋擴(kuò)展水平位移云圖?？梢钥闯?貫穿路徑并非一條連續(xù)曲線,而是在裂紋尖端發(fā)生復(fù)雜破壞,產(chǎn)生多條細(xì)小的次生裂紋,因而采用以往路徑假定的方式處理裂紋擴(kuò)展有本質(zhì)的缺陷。

圖8 多裂紋擴(kuò)展水平位移Fig.8 Horizontal displacements of multi-crack propagation

試驗(yàn)研究[18]表明,巖橋?yàn)?0°時(shí)的擴(kuò)展過程為:預(yù)制裂隙尖端產(chǎn)生翼形裂紋,翼形裂紋拐折向最大壓應(yīng)力方向穩(wěn)定擴(kuò)展,與另一條預(yù)制原生裂紋的翼裂紋發(fā)生相互連接(翼裂—翼裂貫通),最后形成貫通性的宏觀破裂帶。此描述與數(shù)模過程相匹配。圖9為最終數(shù)模結(jié)果與試驗(yàn)結(jié)果對(duì)比。由圖9可以看出采用有限元與離散法混合方法模擬多裂紋擴(kuò)展是可行的。

圖9 數(shù)模與試驗(yàn)結(jié)果對(duì)比Fig.9 Contrast of simulation and experiment results

4 結(jié) 論

(1)巖橋?yàn)?0°的雁行裂紋受對(duì)向擠壓位移載荷作用發(fā)生翼裂—翼裂貫通。

(2)采用有限元與離散法混合方法可有效模擬巖土材料下單、多裂紋的擴(kuò)展過程。有限元與離散法混合方法模擬裂紋擴(kuò)展具有模擬范圍廣,計(jì)算迅速,無須重構(gòu)網(wǎng)格,同步計(jì)算裂紋寬度等優(yōu)點(diǎn)。

[1] 趙益忠,曲連忠,王幸尊,等.不同巖性地層水力壓裂裂縫擴(kuò)展規(guī)律的模擬試驗(yàn)[J].中國石油大學(xué)學(xué)報(bào):自然科學(xué)版,2007,31(3):63-66.

ZHAO Yi-zhong,QU Lian-zhong,WANG Xing-zun,et al.Simulation experiment on prolongation law of hydraulic fracture for different lithologic formations[J].Journal of China University of Petroleum(Edition of Natural Science),2007,31(3):63-66.

[2] 李瑋,閆鐵,畢雪亮.基于分形方法的水力壓裂裂縫起裂擴(kuò)展機(jī)理[J].中國石油大學(xué)學(xué)報(bào):自然科學(xué)版, 2008,32(5):87-91.

LI Wei,YAN Tie,BI Xue-liang.Mechanism of hydraulically created fracture breakdown and propagation based on fractal method[J].Journal of China University of Petroleum(Edition of Natural Science),2008,32(5):87-91.

[3] BOUCHARD P O,BAY F,CHASTEL Y.Numerical modelling of crack propagation:automatic remeshing and comparison of different criteria[J].Computer Methods in Applied Mechanics and Engineering,2003,192(35): 3887-3908.

[4] BELYTSCHKO T,LU Y Y,GU L.Element-free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2):229-256.

[5] RENé de Borst.Numerical aspects of cohesive-zone models[J].Engineering Fracture Mechanics,2003,70(14): 1743-1757.

[6] DUARTE C A,HAMZEH O N,LISZKA T J,et al.A generalized finite element method for the simulation of three-dimensional dynamic crack propagation[J].Computer Methods in Applied Mechanics and Engineering, 2001,190(15):2227-2262.

[7] MISHRA B K,RAJ K Rajamani.The discrete element method for the simulation of ball mills[J].Applied Mathematical Modelling,1992,16(11):598-604.

[8] LANRU Jing,STAPHANSSON Ove.Fundamentals of discrete element methods for rock engineering[M]. Elesevier,2007.

[9] MUNJIZA Ante.The combined finite-discrete element method[M].Queen Mary,London:University of London,2004.

[10] LI Shi-hai,LIU Xiao-yu,LIU Tian-ping,et al.Continuum-based discrete element method and its applications: UK-China Summer School/International Symposium on DEM[C].Beijing:c2008.

[11] BRADY B H G,BROWN E T.Rock mechanics for underground mining[M].2nd ed.London:Chapama& Hall,1993:106-108.

[12] CUNDALL P A,BROWN E T.Distinct element models of soil and rock structure[J].Analytical and Computational Methods in Engineering Rock Mechanics,1987, 4:129-163.

[13] TONG P.An adaptive dynamic relaxation method for static problems[J].Computational Mechanics,1986,1, 127-140.

[14] CUNDALL P A,HART R D.Numerical modelling of discontinua[J].Engineering Computations,1992,9 (2):101-113.

[15] 馮春,李世海,劉曉宇.半彈簧接觸模型及其在邊坡破壞計(jì)算中的應(yīng)用[J].力學(xué)學(xué)報(bào),2011,43(1):184-192.

FENG Chun,LI Shi-hai,LIU Xiao-yu.Semi-spring contact model and its application to failure simulation of slope[J].Acta Mechanica Sinica,2011,43(1):184-192.

[16] 常曉林,胡超,馬剛,等.模擬巖體失效全過程的連續(xù)-非連續(xù)變形體離散元方法及應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2011,30(10):2004-2011.

CHANG Xiao-lin,HU Chao,MA Gang,et al.Continuous-discontinuous deformable discrete element method to simulate the whole failure process of rock masses and application[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(10):2004-2011.

[17] 李強(qiáng),楊慶,欒茂田,等.曲線翼型裂紋擴(kuò)展路徑理論分析及試驗(yàn)驗(yàn)證[J].巖土力學(xué),2010,31(2):345-349.

LI Qiang,YANG Qing,LUAN Mao-tian,et al.Study of curved wing crack path by theory and testing methods [J].Rock and Soil Mechanics,2010,31(2):345-349.

[18] 趙延林.裂隙巖體滲流-損傷-斷裂耦合理論及應(yīng)用研究[D].長沙:中南大學(xué)資源與安全工程學(xué)院, 2009.

ZHAO Yan-lin.Coupling theory of seepage-damagefracture in fracture rock masses and its application[D]. School of Resources and Safety Engineering in Central South University,2009.

(編輯 沈玉英)

Application of combined finite-discrete element method for crack propagation

SUN Xiang,LIU Chuan-qi,XUE Shi-feng
(College of Pipeline and Civil Engineering in China University of Petroleum,Qingdao 266580,China)

The process of crack growth was simulated based on the combined finite-discrete element method.The media existing primary cracks was divided by elements.The unit and unit interface were calculated using the finite element and discrete element method,respectively.The continuum was converted into discontinuum by transformation of connection types between elements.The contacts were detected using the half-spring method in plane.The motion equation was solved in an explicit iterative scheme.Updating the coordinates of elements,the process of crack propagation was simulated.The simulation results of a single crack and en echelon cracks under uniaxial compression were contrasted with experimental results.The results show that the process of crack propagation can be effectively simulated using the combined finite-discrete element method.The en echelon cracks with 90°rock bridge under uniaxial compression were linked by wing cracks.

finite element method；discrete element method；crack propagation；numerical simulation

TU 45

A

1673-5005(2013)03-0126-05

10.3969/j.issn.1673-5005.2013.02.022

2012-11-22

中國石油重大科技專項(xiàng)(2012-ZG-009)

孫翔(1973-),男,高級(jí)工程師,博士研究生,主要從事工程力學(xué)研究。E-mail:sunx839@126.com。