郭德平，李 錚，彭森林，曾志凱，吳岱峰

(1. 敘鎮(zhèn)鐵路有限責(zé)任公司， 云南 昭通 657900； 2. 武漢大學(xué) 土木建筑工程學(xué)院， 武漢 430072； 3. 重慶市城市建設(shè)投資(集團(tuán))有限公司， 重慶 400023； 4. 重慶大學(xué) 土木工程學(xué)院， 重慶 400045)

隨著21世紀(jì)我國(guó)基礎(chǔ)設(shè)施大規(guī)模建設(shè)的進(jìn)行，西部大開(kāi)發(fā)戰(zhàn)略的實(shí)施以及世界經(jīng)濟(jì)危機(jī)以來(lái)國(guó)家對(duì)基礎(chǔ)設(shè)施建設(shè)的投資，我國(guó)的鐵路、公路、大中型水電站建設(shè)以及南水北調(diào)、西氣東輸?shù)裙こ虒⒂写罅康拈L(zhǎng)大隧道需要建設(shè).因此，隧道掘進(jìn)機(jī)(TBM)在我國(guó)的應(yīng)用前景非常廣闊，我國(guó)對(duì)掘進(jìn)機(jī)及其技術(shù)的需求猛增.

在工程建設(shè)和運(yùn)營(yíng)過(guò)程中，結(jié)構(gòu)或者裂隙巖體會(huì)受到多種形式的動(dòng)力作用(如爆破、地震等)，在動(dòng)荷載的作用下更容產(chǎn)生失穩(wěn)破壞，故針對(duì)裂紋動(dòng)態(tài)擴(kuò)展的研究引起了廣泛的關(guān)注[1-3].Sun等[4]基于有限元方法比較了顯式時(shí)間積分和隱式時(shí)間積分在分析線(xiàn)性和非線(xiàn)性動(dòng)力問(wèn)題中的區(qū)別.2012年，Nilsson等[5]在黏結(jié)單元中嵌入裂紋，分析了動(dòng)荷載作用于彈塑性厚板的動(dòng)力響應(yīng).

擴(kuò)展有限單元法(XFEM)[6-9]是近十五年發(fā)展起來(lái)的一種在常規(guī)有限元框架內(nèi)求解不連續(xù)問(wèn)題，特別是裂紋擴(kuò)展問(wèn)題的數(shù)值方法.其原理是在裂紋影響區(qū)域的單元結(jié)點(diǎn)上用裂尖漸近位移場(chǎng)函數(shù)及跳躍函數(shù)加強(qiáng)傳統(tǒng)有限元的基，以考慮由于裂紋存在引起的位移不連續(xù)性，繼承了標(biāo)準(zhǔn)有限元的所有優(yōu)點(diǎn)，但其所使用的網(wǎng)格與結(jié)構(gòu)內(nèi)部的幾何和物理界面無(wú)關(guān)，從而避免了常規(guī)有限元方法中的網(wǎng)格重構(gòu)，不需要裂紋面與有限元網(wǎng)格一致，不需要在裂縫周?chē)贾酶呙芏染W(wǎng)格，大大簡(jiǎn)化了裂紋擴(kuò)展的分析過(guò)程.能夠很容易刻畫(huà)出裂紋面上位移的強(qiáng)不連續(xù)性質(zhì)和裂紋尖端的應(yīng)力奇異性，并且無(wú)需重新劃分網(wǎng)格.2001年，Stolarska等[10]將水平集函數(shù)引入XFEM來(lái)描述裂紋的位置和裂紋擴(kuò)展之后的更新位置.裂紋的幾何位置能夠很容易用兩個(gè)零水平集函數(shù)來(lái)表達(dá)，這兩個(gè)零水平集函數(shù)在裂紋尖端處彼此正交.同時(shí)，隨著裂紋的擴(kuò)展，裂紋面和裂紋尖端處需要富集的節(jié)點(diǎn)能夠?qū)崟r(shí)更新.

Zhou等[11]在推導(dǎo)擴(kuò)展有限元算法的基礎(chǔ)上，分析了應(yīng)力強(qiáng)度因子的J積分計(jì)算方法及積分區(qū)域的選取.基于擴(kuò)展有限元法對(duì)I型裂紋進(jìn)行了計(jì)算，有限元網(wǎng)格獨(dú)立于裂紋面，無(wú)需在裂紋尖端加密網(wǎng)格.Chen等[12]引入裂紋交叉匯合加強(qiáng)函數(shù)以分析多裂紋交叉匯合過(guò)程，并且在裂紋附近區(qū)域使用廣義形函數(shù)，并引入線(xiàn)增函數(shù)消除混合單元，可有效提高裂紋附近的精度.黃宏偉等[13]在這些研究的基礎(chǔ)上，采用擴(kuò)展有限元研究了襯砌在主要影響因素作用下的裂縫分布規(guī)律、裂縫擴(kuò)展過(guò)程、裂縫外觀(guān)表現(xiàn)形式及發(fā)生機(jī)制.阮濱等[14]基于擴(kuò)展有限元法對(duì)均質(zhì)土壩壩頂?shù)某跏剂鸭y擴(kuò)展路徑進(jìn)行了模擬，研究結(jié)果表明擴(kuò)展有限元法對(duì)網(wǎng)格劃分的要求比較高，網(wǎng)格須均勻，網(wǎng)格的疏密程度對(duì)計(jì)算的結(jié)果影響不大.Menouillard等[15]提出利用無(wú)網(wǎng)格近似方法的XFEM富集方案，該方法采用顯式時(shí)間積分方案來(lái)分析動(dòng)力過(guò)程.Wen等[16]基于改進(jìn)的擴(kuò)展有限單元法研究了裂紋動(dòng)態(tài)擴(kuò)展問(wèn)題的計(jì)算精度和算法穩(wěn)定性問(wèn)題.

但是標(biāo)準(zhǔn)有限元在處理時(shí)間積分時(shí)，在裂紋不斷擴(kuò)展的過(guò)程中整體剛度矩陣的自由度也會(huì)不斷的增大，從而導(dǎo)致迭代計(jì)算無(wú)法進(jìn)行.本文基于擴(kuò)展有限單元法模擬動(dòng)態(tài)裂紋擴(kuò)展的方法，提出了新的時(shí)間積分方案.將所有節(jié)點(diǎn)都富集Heaviside函數(shù)和裂紋尖端的漸近位移場(chǎng)函數(shù)，即每個(gè)節(jié)點(diǎn)都有12個(gè)自由度，從而使得總體剛度矩陣式保持一致，避免迭代計(jì)算式無(wú)法進(jìn)行的情況.提出了一種稀疏矩陣技術(shù)來(lái)解決矩陣所占內(nèi)存大和計(jì)算時(shí)間長(zhǎng)的問(wèn)題.

1 動(dòng)力擴(kuò)展有限元的控制方程

1.1 控制方程的強(qiáng)形式

圖1 在動(dòng)荷載作用下含有裂紋的二維區(qū)域示意圖Fig.1 Diagram of two-dimensional domain containing a crack subjected to dynamic loads

動(dòng)力平衡方程的強(qiáng)形式為

(1)

該二維區(qū)域的應(yīng)力邊界條件和位移邊界條件為

σ·nt=T, 在Γt上

(2)

σ·nc=0, 在Γc上

(3)

U·nu=0, 在Γu上

(4)

式中：nt為邊界Γt的單位外法向量；nc為邊界Γc的單位外法向量；nu為邊界Γu的單位外法向量.

在小應(yīng)變和小位移的條件下，幾何方程可以表達(dá)為

(5)

根據(jù)廣義胡克定律，σ與應(yīng)變張量的之間的關(guān)系可以表達(dá)為

σ=C∶ε

(6)

式中：C為Hooke張量.

1.2 控制方程的等效弱形式

根據(jù)虛位移原理，引入虛位移δU，式(1)可以改寫(xiě)為

(7)

再根據(jù)分部積分，式(7)可以化簡(jiǎn)為

(8)

由變分原理將式(8)的第一項(xiàng)進(jìn)行變換為

(9)

因此，該二維區(qū)域的能量系統(tǒng)的泛函數(shù)可以用下式表達(dá)

(10)

對(duì)式(10)進(jìn)行變分，并結(jié)合式(8)可以得到：

(11)

因此有，

(12)

1.3 區(qū)域的空間離散

水平集法是裂紋和非連續(xù)面追蹤強(qiáng)有力的工具.Zhou等[17]首次引入水平集法來(lái)研究材料內(nèi)部界面的追蹤定位和形狀描述.如圖2所示，一根裂紋可以用2個(gè)在裂紋尖端彼此相互正交的零水平集函數(shù)φ(x)和ψ(x)表示.

圖2 用于追蹤裂紋的兩個(gè)零水平集函數(shù)原理圖Fig.2 Schematic diagram of two zero-level set functions used to track cracks

Stolarska等[10]給出了水平集法的更新算法，水平集函數(shù)φ(x)的更新算法為

(13)

式中：φi+1(x)為水平集函數(shù)φ(x)的更新值；Δt為時(shí)間增量；ui和vi分別為裂紋尖端沿x和y方向擴(kuò)展的速度.

對(duì)于水平集函數(shù)ψ(x)，其更新算法為

ψi+1(x)=

(14)

式中：F=[FxFy]為裂紋尖端在裂紋面外法線(xiàn)方向的速度矢量；(xi,yi)為裂紋尖端的坐標(biāo).

如圖3所示，區(qū)域Ω被離散化成ne個(gè)單元，I表示該區(qū)域所有節(jié)點(diǎn)，J表示裂紋面上被Heaviside函數(shù)富集的節(jié)點(diǎn)集合，由藍(lán)色正方形標(biāo)識(shí)，K表示裂紋尖端上被漸近位移場(chǎng)函數(shù)富集的節(jié)點(diǎn)集合，由紅色圓圈形標(biāo)識(shí).

圖3 擴(kuò)展有限單元法的節(jié)點(diǎn)富集方案Fig.3 Enrichment scheme of XFEM

根據(jù)文獻(xiàn)[6，18]的研究，單元位移場(chǎng)可以表達(dá)為

(15)

(16)

(17)

c4x1c4y1c4x2c4y2c4x3c4y3c4x4c4y4]T

(18)

式中：aix和aiy分別為第i個(gè)節(jié)點(diǎn)在x和y方向的自由度位移值；bix和biy分別為第i個(gè)Heaviside函數(shù)富集節(jié)點(diǎn)在x和y方向的附加自由度位移值；cixj和ciyj分別為第i個(gè)節(jié)點(diǎn)的第j個(gè)附加自由度沿x和y方向的附加自由度位移值.

根據(jù)Wu等[19]的研究，線(xiàn)彈性材料中的I型、II 型及 III 型裂紋尖端漸近位移場(chǎng)均可由幾個(gè)特定的基函數(shù)組成的函數(shù)形式來(lái)表達(dá)，為了能體現(xiàn)出裂紋尖端漸近位移場(chǎng)的奇異性，將該組基函數(shù)引入裂紋尖端的位移場(chǎng)計(jì)算，如下：

(19)

(20)

(21)

式中：(r,θ)是以裂紋尖端為原點(diǎn)建立的極坐標(biāo)系的坐標(biāo)值；(xtip,ytip)是絕對(duì)坐標(biāo)系中裂紋尖端的坐標(biāo)值，坐標(biāo)系的建立如圖4所示，x′O′y′為以裂紋尖端建立的局部坐標(biāo)系，其中α為裂紋中心線(xiàn)與絕對(duì)坐標(biāo)系x軸的夾角.

圖4 裂紋面上的兩套坐標(biāo)系Fig.4 Two coordinate systems of crack surface

裂紋尖端富集后的形函數(shù)Nφ(x)可以表達(dá)為

(22)

φj為尖端位移場(chǎng).將式(15)代入式(5)，可以得到單元的應(yīng)變張量εe(x)如下：

(23)

再將式(15)、(23)代入式(10)中可以得到：

(24)

將式(24)代入式(11)并化簡(jiǎn)，可以得到：

(25)

式中：

考慮到3種類(lèi)型的自由度位移值Ua、Ub及Uc式彼此相互獨(dú)立的，結(jié)合變分原理，將式(25)代入式(12)，可以得到：

(26)

根據(jù)式(26)，擴(kuò)展有限單元法的運(yùn)動(dòng)學(xué)方程為

(27)

1.4 時(shí)間積分方案

在擴(kuò)展有限單元法的計(jì)算中，時(shí)間積分會(huì)遇到困難.時(shí)間積分基于迭代的算法，而在裂紋不斷擴(kuò)展的過(guò)程中整體剛度矩陣的自由度也會(huì)不斷的增大，從而導(dǎo)致迭代計(jì)算無(wú)法進(jìn)行.在每次迭代中，時(shí)間積分會(huì)涉及當(dāng)前步的位移場(chǎng)Ut和下一步的位移場(chǎng)Ut+Δt，具體表達(dá)式為

(28)

(29)

本文提出的時(shí)間積分方案是將所有節(jié)點(diǎn)都富集Heaviside函數(shù)和裂紋尖端的漸近位移場(chǎng)函數(shù)，即每個(gè)節(jié)點(diǎn)都有12個(gè)自由度，從而使得總體剛度矩陣式保持一致，而避免迭代計(jì)算式無(wú)法進(jìn)行.但是，將所有節(jié)點(diǎn)都富集Heaviside函數(shù)和裂紋尖端的漸近位移場(chǎng)函數(shù)的代價(jià)是剛度矩陣的階數(shù)增多，從而使得參與運(yùn)算的矩陣所占內(nèi)存和計(jì)算時(shí)間急劇增加.為此，本文提出了一種稀疏矩陣技術(shù)來(lái)解決矩陣所占內(nèi)存大和計(jì)算時(shí)間長(zhǎng)的問(wèn)題.

對(duì)于兩個(gè)常規(guī)有限元的自由度，若該節(jié)點(diǎn)的某個(gè)方向被約束(約束狀態(tài))，那么所對(duì)應(yīng)的自由度值被定義為0，若該節(jié)點(diǎn)的某個(gè)方向沒(méi)有被約束(激活狀態(tài))，那么所對(duì)應(yīng)的自由度值作為未知數(shù)參與計(jì)算.對(duì)于兩個(gè)Heaviside函數(shù)富集的自由度，若該節(jié)點(diǎn)不屬于集合J時(shí)(約束狀態(tài))，這兩個(gè)自由度位移值被定義為0，若該節(jié)點(diǎn)屬于集合J時(shí)(激活狀態(tài))，這兩個(gè)自由度位移值作為未知數(shù)參與計(jì)算.對(duì)于另外8個(gè)裂紋尖端漸近位移場(chǎng)函數(shù)富集的自由度，若該節(jié)點(diǎn)不屬于集合K時(shí)(約束狀態(tài))，這8個(gè)自由度位移值被定義為0，若該節(jié)點(diǎn)屬于集合K時(shí)(激活狀態(tài))，這8個(gè)自由度位移值作為未知數(shù)參與計(jì)算.

根據(jù)上述定義，可以得到一個(gè)將每個(gè)節(jié)點(diǎn)有兩個(gè)自由度的體系(簡(jiǎn)稱(chēng)二維體系)映射到每個(gè)節(jié)點(diǎn)有12個(gè)自由度的體系(簡(jiǎn)稱(chēng)十二維體系)中的轉(zhuǎn)換矩陣，該轉(zhuǎn)換矩陣的構(gòu)造方法如下：

(30)

式中：變量n和m根據(jù)模型的節(jié)點(diǎn)個(gè)數(shù)來(lái)定.二維體系中的第i個(gè)自由度映射于12維體系中的第j個(gè)自由度，當(dāng)該自由度是激活狀態(tài)時(shí)，該矩陣中所有元素為1，當(dāng)該自由度是約束狀態(tài)時(shí)，該矩陣中所有元素為0.

(31)

(32)

(33)

(34)

那么，擴(kuò)展有限單元法在12維體系中的運(yùn)動(dòng)方程為

(35)

對(duì)于時(shí)間積分，本文采用Newmark隱式時(shí)間積分方案如下：

(36)

(37)

式中：ξ和η分別為Newmark隱式時(shí)間積分方案中的參數(shù).

在Newmark隱式時(shí)間積分方案中，t+Δt時(shí)刻的位移場(chǎng)為

(38)

聯(lián)合式(36)～(38)，可以得到12維體系的位移場(chǎng)的求解方程式：

(39)

該方程的初始條件是：

(40)

(41)

(42)

在12維體系中，t+Δt時(shí)刻的加速度場(chǎng)和速度場(chǎng)為

(43)

(44)

根據(jù)Newmark時(shí)間積分方法，式(39)求解穩(wěn)定的條件是ξ和η應(yīng)滿(mǎn)足[20]：

ξ≥0.5，η≥0.25(0.5+ξ)2

(45)

2 實(shí)例與分析

2.1 實(shí)例1

如圖5所示，平面板的長(zhǎng)L=10 m，高H=10 m.預(yù)制裂紋在板的左側(cè)中部，其長(zhǎng)度l0=5 m，距離上邊界高度h=2 m.板的下邊界豎向被約束，左下角被水平方向約束.板的材料屬性為：彈性模量E=210 GPa，泊松比ν=0.3，密度ρ=8 000 kg/m3，阻尼系數(shù)μ=0.05.上邊界受到動(dòng)荷載σ(t)的作用，作用力的函數(shù)表達(dá)式為

圖5 實(shí)例1的幾何布置和荷載條件(m)Fig.5 Geometric layout and load conditions of example 1 (m)

σ(t)=σ0f0(t)

(46)

根據(jù)Wang等[22-23]的研究，裂紋尖端的應(yīng)力強(qiáng)度因子KI的解析解為

(47)

(48)

(49)

采用本文提出的方法計(jì)算不同時(shí)刻裂紋尖端的動(dòng)態(tài)應(yīng)力強(qiáng)度因子，再將其歸一化.采用了3種不同的時(shí)間步Δt=10 μs，Δt=20 μs及Δt=50 μs，所得數(shù)值結(jié)果如圖6所示.裂紋尖端動(dòng)態(tài)應(yīng)力強(qiáng)度因子在時(shí)間t=0～τc內(nèi)幾乎為0，因?yàn)檫@段時(shí)間應(yīng)力波還沒(méi)有到達(dá)裂紋的尖端.當(dāng)t>τc時(shí)，裂紋尖端的動(dòng)態(tài)應(yīng)力強(qiáng)度因子開(kāi)始逐漸增大.從如圖6可以看出，本方法計(jì)算的結(jié)果與解析解的結(jié)果吻合.本方法計(jì)算的的動(dòng)態(tài)應(yīng)力強(qiáng)度因子具有一定的震蕩性，但是震蕩性較小，如圖7所示(KI,S為震蕩值)，并且震蕩性隨著時(shí)間的增長(zhǎng)而逐漸衰弱.本方法計(jì)算的x方向應(yīng)力σx的分布如圖8所示，在裂紋尖端的應(yīng)力集中比較明顯.

圖6 KI的歷時(shí)曲線(xiàn)圖Fig.6 Time history of KI

圖7 不同時(shí)間步KI的震蕩歷時(shí)曲線(xiàn)圖Fig.7 Time history of oscillation of KI at different time steps

圖8 σx分布云圖(Δt=10 μs，t=3τc)Fig.8 Contour of σx (Δt=10 μs，t=3τc)

2.2 實(shí)例2

圖9 實(shí)例2的幾何布置和荷載條件(m)Fig.9 Geometric layout and load conditions of example 2 (m)

σ(t)=σ0f0(t)

(50)

在本實(shí)例中，研究了25×100，50×200和70×280共3種不同的網(wǎng)格密度(行數(shù)×列數(shù)).不同網(wǎng)格密度下裂紋動(dòng)態(tài)應(yīng)力強(qiáng)度因子的歷時(shí)曲線(xiàn)如圖10所示，網(wǎng)格密度為25×100和50×200的裂紋動(dòng)態(tài)應(yīng)力強(qiáng)度因子的相關(guān)系數(shù)為 0.999 3，網(wǎng)格密度為50×200和70×280的裂紋動(dòng)態(tài)應(yīng)力強(qiáng)度因子的相關(guān)系數(shù)為 0.999 4，網(wǎng)格密度為25×100和70×280的裂紋動(dòng)態(tài)應(yīng)力強(qiáng)度因子的相關(guān)系數(shù)為 0.998 5.當(dāng)t<0.5 ms時(shí)，裂紋尖端的應(yīng)力強(qiáng)度因子幾乎為0；當(dāng)t≥0.5 ms時(shí)，裂紋尖端的應(yīng)力強(qiáng)度因子隨著時(shí)間的增加而逐漸增大，并呈現(xiàn)很小的震蕩特性.變形后的簡(jiǎn)支梁如圖11所示，裂紋基本沿著直線(xiàn)向上擴(kuò)展.該簡(jiǎn)支梁變形后的應(yīng)力分布如圖12所示，裂紋尖端應(yīng)力集中突出，水平方向應(yīng)力值達(dá)2.09×103MPa.

圖10 不同網(wǎng)格密度下KI的歷時(shí)曲線(xiàn)圖Fig.10 Time history of KI at different mesh densities

圖11 變形后的簡(jiǎn)支梁(網(wǎng)格密度為50×200)Fig.11 Deformed mesh for mesh density (at a mesh density of 50×200)

圖12 t=2 ms時(shí)刻的σx云圖(網(wǎng)格密度為50×200)Fig.12 Contour of σx at t=2 ms (at a mesh density of 50×200)

3 結(jié)語(yǔ)

標(biāo)準(zhǔn)有限元在處理時(shí)間積分時(shí)，在裂紋不斷擴(kuò)展的過(guò)程中整體剛度矩陣的自由度也會(huì)不斷增大，從而導(dǎo)致迭代計(jì)算無(wú)法進(jìn)行.本文提出的基于擴(kuò)展有限單元法模擬動(dòng)態(tài)裂紋擴(kuò)展的方法提出了新的時(shí)間積分方案，將所有節(jié)點(diǎn)都富集Heaviside函數(shù)和裂紋尖端的漸近位移場(chǎng)函數(shù)，即每個(gè)節(jié)點(diǎn)都有12個(gè)自由度，從而使得總體剛度矩陣式保持一致，避免迭代計(jì)算式無(wú)法進(jìn)行.并且提出了一種稀疏矩陣技術(shù)來(lái)解決矩陣所占內(nèi)存大和計(jì)算時(shí)間長(zhǎng)的問(wèn)題.數(shù)值計(jì)算的結(jié)果表明，利用空間變換理論計(jì)算動(dòng)力問(wèn)題時(shí)，施加的荷載以膨脹波的形式傳遞，裂紋尖端的應(yīng)力強(qiáng)度因子比荷載施加的時(shí)間稍有延遲，數(shù)值解的結(jié)果能夠與解析的結(jié)果很好地吻合.線(xiàn)性荷載的數(shù)值結(jié)果比瞬時(shí)脈沖荷載的數(shù)值結(jié)果的震蕩性更小，應(yīng)力分布更加平穩(wěn)光滑.不同網(wǎng)格密度條件下的數(shù)值計(jì)算結(jié)果相差不大，說(shuō)明該方法計(jì)算裂紋擴(kuò)展對(duì)網(wǎng)格的依賴(lài)性較小.