基于Newmark隱式時間積分方案的裂紋動態擴展的數值計算方法

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

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

隨著21世紀我國基礎設施大規模建設的進行，西部大開發戰略的實施以及世界經濟危機以來國家對基礎設施建設的投資，我國的鐵路、公路、大中型水電站建設以及南水北調、西氣東輸等工程將有大量的長大隧道需要建設.因此，隧道掘進機(TBM)在我國的應用前景非常廣闊，我國對掘進機及其技術的需求猛增.

在工程建設和運營過程中，結構或者裂隙巖體會受到多種形式的動力作用(如爆破、地震等)，在動荷載的作用下更容產生失穩破壞，故針對裂紋動態擴展的研究引起了廣泛的關注[1-3].Sun等[4]基于有限元方法比較了顯式時間積分和隱式時間積分在分析線性和非線性動力問題中的區別.2012年，Nilsson等[5]在黏結單元中嵌入裂紋，分析了動荷載作用于彈塑性厚板的動力響應.

擴展有限單元法(XFEM)[6-9]是近十五年發展起來的一種在常規有限元框架內求解不連續問題，特別是裂紋擴展問題的數值方法.其原理是在裂紋影響區域的單元結點上用裂尖漸近位移場函數及跳躍函數加強傳統有限元的基，以考慮由于裂紋存在引起的位移不連續性，繼承了標準有限元的所有優點，但其所使用的網格與結構內部的幾何和物理界面無關，從而避免了常規有限元方法中的網格重構，不需要裂紋面與有限元網格一致，不需要在裂縫周圍布置高密度網格，大大簡化了裂紋擴展的分析過程.能夠很容易刻畫出裂紋面上位移的強不連續性質和裂紋尖端的應力奇異性，并且無需重新劃分網格.2001年，Stolarska等[10]將水平集函數引入XFEM來描述裂紋的位置和裂紋擴展之后的更新位置.裂紋的幾何位置能夠很容易用兩個零水平集函數來表達，這兩個零水平集函數在裂紋尖端處彼此正交.同時，隨著裂紋的擴展，裂紋面和裂紋尖端處需要富集的節點能夠實時更新.

Zhou等[11]在推導擴展有限元算法的基礎上，分析了應力強度因子的J積分計算方法及積分區域的選取.基于擴展有限元法對I型裂紋進行了計算，有限元網格獨立于裂紋面，無需在裂紋尖端加密網格.Chen等[12]引入裂紋交叉匯合加強函數以分析多裂紋交叉匯合過程，并且在裂紋附近區域使用廣義形函數，并引入線增函數消除混合單元，可有效提高裂紋附近的精度.黃宏偉等[13]在這些研究的基礎上，采用擴展有限元研究了襯砌在主要影響因素作用下的裂縫分布規律、裂縫擴展過程、裂縫外觀表現形式及發生機制.阮濱等[14]基于擴展有限元法對均質土壩壩頂的初始裂紋擴展路徑進行了模擬，研究結果表明擴展有限元法對網格劃分的要求比較高，網格須均勻，網格的疏密程度對計算的結果影響不大.Menouillard等[15]提出利用無網格近似方法的XFEM富集方案，該方法采用顯式時間積分方案來分析動力過程.Wen等[16]基于改進的擴展有限單元法研究了裂紋動態擴展問題的計算精度和算法穩定性問題.

但是標準有限元在處理時間積分時，在裂紋不斷擴展的過程中整體剛度矩陣的自由度也會不斷的增大，從而導致迭代計算無法進行.本文基于擴展有限單元法模擬動態裂紋擴展的方法，提出了新的時間積分方案.將所有節點都富集Heaviside函數和裂紋尖端的漸近位移場函數，即每個節點都有12個自由度，從而使得總體剛度矩陣式保持一致，避免迭代計算式無法進行的情況.提出了一種稀疏矩陣技術來解決矩陣所占內存大和計算時間長的問題.

1 動力擴展有限元的控制方程

1.1 控制方程的強形式

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

動力平衡方程的強形式為

(1)

該二維區域的應力邊界條件和位移邊界條件為

σ·nt=T, 在Γt上

(2)

σ·nc=0, 在Γc上

(3)

U·nu=0, 在Γu上

(4)

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

在小應變和小位移的條件下，幾何方程可以表達為

(5)

根據廣義胡克定律，σ與應變張量的之間的關系可以表達為

σ=C∶ε

(6)

式中：C為Hooke張量.

1.2 控制方程的等效弱形式

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

(7)

再根據分部積分，式(7)可以化簡為

(8)

由變分原理將式(8)的第一項進行變換為

(9)

因此，該二維區域的能量系統的泛函數可以用下式表達

(10)

對式(10)進行變分，并結合式(8)可以得到：

(11)

因此有，

(12)

1.3 區域的空間離散

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

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

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

(13)

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

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

ψi+1(x)=

(14)

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

如圖3所示，區域Ω被離散化成ne個單元，I表示該區域所有節點，J表示裂紋面上被Heaviside函數富集的節點集合，由藍色正方形標識，K表示裂紋尖端上被漸近位移場函數富集的節點集合，由紅色圓圈形標識.

圖3 擴展有限單元法的節點富集方案Fig.3 Enrichment scheme of XFEM

根據文獻[6，18]的研究，單元位移場可以表達為

(15)

(16)

(17)

c4x1c4y1c4x2c4y2c4x3c4y3c4x4c4y4]T

(18)

式中：aix和aiy分別為第i個節點在x和y方向的自由度位移值；bix和biy分別為第i個Heaviside函數富集節點在x和y方向的附加自由度位移值；cixj和ciyj分別為第i個節點的第j個附加自由度沿x和y方向的附加自由度位移值.

根據Wu等[19]的研究，線彈性材料中的I型、II 型及 III 型裂紋尖端漸近位移場均可由幾個特定的基函數組成的函數形式來表達，為了能體現出裂紋尖端漸近位移場的奇異性，將該組基函數引入裂紋尖端的位移場計算，如下：

(19)

(20)

(21)

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

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

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

(22)

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

(23)

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

(24)

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

(25)

式中：

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

(26)

根據式(26)，擴展有限單元法的運動學方程為

(27)

1.4 時間積分方案

在擴展有限單元法的計算中，時間積分會遇到困難.時間積分基于迭代的算法，而在裂紋不斷擴展的過程中整體剛度矩陣的自由度也會不斷的增大，從而導致迭代計算無法進行.在每次迭代中，時間積分會涉及當前步的位移場Ut和下一步的位移場Ut+Δt，具體表達式為

(28)

(29)

本文提出的時間積分方案是將所有節點都富集Heaviside函數和裂紋尖端的漸近位移場函數，即每個節點都有12個自由度，從而使得總體剛度矩陣式保持一致，而避免迭代計算式無法進行.但是，將所有節點都富集Heaviside函數和裂紋尖端的漸近位移場函數的代價是剛度矩陣的階數增多，從而使得參與運算的矩陣所占內存和計算時間急劇增加.為此，本文提出了一種稀疏矩陣技術來解決矩陣所占內存大和計算時間長的問題.

對于兩個常規有限元的自由度，若該節點的某個方向被約束(約束狀態)，那么所對應的自由度值被定義為0，若該節點的某個方向沒有被約束(激活狀態)，那么所對應的自由度值作為未知數參與計算.對于兩個Heaviside函數富集的自由度，若該節點不屬于集合J時(約束狀態)，這兩個自由度位移值被定義為0，若該節點屬于集合J時(激活狀態)，這兩個自由度位移值作為未知數參與計算.對于另外8個裂紋尖端漸近位移場函數富集的自由度，若該節點不屬于集合K時(約束狀態)，這8個自由度位移值被定義為0，若該節點屬于集合K時(激活狀態)，這8個自由度位移值作為未知數參與計算.

根據上述定義，可以得到一個將每個節點有兩個自由度的體系(簡稱二維體系)映射到每個節點有12個自由度的體系(簡稱十二維體系)中的轉換矩陣，該轉換矩陣的構造方法如下：

(30)

式中：變量n和m根據模型的節點個數來定.二維體系中的第i個自由度映射于12維體系中的第j個自由度，當該自由度是激活狀態時，該矩陣中所有元素為1，當該自由度是約束狀態時，該矩陣中所有元素為0.

(31)

(32)

(33)

(34)

那么，擴展有限單元法在12維體系中的運動方程為

(35)

對于時間積分，本文采用Newmark隱式時間積分方案如下：

(36)

(37)

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

在Newmark隱式時間積分方案中，t+Δt時刻的位移場為

(38)

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

(39)

該方程的初始條件是：

(40)

(41)

(42)

在12維體系中，t+Δt時刻的加速度場和速度場為

(43)

(44)

根據Newmark時間積分方法，式(39)求解穩定的條件是ξ和η應滿足[20]：

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

(45)

2 實例與分析

2.1 實例1

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

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

σ(t)=σ0f0(t)

(46)

根據Wang等[22-23]的研究，裂紋尖端的應力強度因子KI的解析解為

(47)

(48)

(49)

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

圖6 KI的歷時曲線圖Fig.6 Time history of KI

圖7 不同時間步KI的震蕩歷時曲線圖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 實例2

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

σ(t)=σ0f0(t)

(50)

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

圖10 不同網格密度下KI的歷時曲線圖Fig.10 Time history of KI at different mesh densities

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

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

3 結語

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