固體結構損傷破壞統一相場理論、算法和應用1)

吳建營

(華南理工大學亞熱帶建筑科學國家重點實驗室,廣州 510641)

引言

自1920 年Griffith[1]創立線彈性斷裂力學一個世紀以來,工程材料和結構的損傷破壞問題一直是國際固體力學領域方興未艾的研究熱點和難點.一方面,伴隨著研究人員在應力強度因子[2]、內聚裂縫模型[3-5]、J 積分[6]等方面取得的重要進展,斷裂力學理論漸趨成熟[7].另一方面,Kachanov[8]于1958年提出了金屬高溫蠕變的損傷模型,直至20 世紀80年代引入裂縫帶理論[9]處理含應變軟化段的本構關系,損傷力學理論也獲得了長足發展[10-11].將斷裂力學理論或損傷力學理論與有限元等數值方法相結合,研究人員分別發展了損傷破壞分析的非連續方法[12]或連續方法[13],損傷與破壞力學逐漸演變為固體力學的重要分支,并在工程結構安全分析和優化設計等方面發揮了重要作用[14-15].

然而,無論是斷裂力學或損傷力學理論,還是相應的非連續或連續數值方法,均存在一系列長期難以解決的難題.斷裂力學理論自身缺乏裂縫起裂和分叉準則,研究人員對最大環向應力、最大能量釋放率、最大應變能密度等裂縫擴展方向判據的認識也有明顯分歧;數值上,直接處理位移跳躍的非連續方法需要繁瑣的網格重劃分(對于界面單元)或復雜的裂縫幾何表征、裂縫路徑跟蹤等(對于擴展有限元、內嵌裂縫方法等)[16].損傷力學模型可以自動處理裂縫起裂、擴展和分叉等復雜行為,但損傷變量缺乏明確的物理意義,與裂縫之間也不存在一一對應關系,故而難以準確預測裂縫位置、裂縫寬度等工程中極為關心的局部信息;將含應變軟化段的經典損傷力學模型應用于結構損傷破壞的數值分析時,會遭遇網格敏感性、應力鎖死等難題[17];即便是非局部損傷模型(無論積分型或微分型),也存在材料內部尺度隨意選取、邊界效應處理困難、裂縫帶無限變寬等眾多缺點[18].事實上,正如2009 年底美國三院院士Bazˇant[19]在鐵木辛柯獎章頒獎禮上的致辭中指出:“固體和結構的損傷破壞問題,與流體力學中的湍流,并稱為21 世紀工程科學的兩大難題”;20 世紀90 年代中期,中國科學院院士楊衛[20]也有類似的表述.

20 世紀末、21 世紀初,研究人員提出了兩類新的固體損傷破壞分析方法,即近場動力學[21]和相場斷裂模型[22-23].近場動力學方法引入非局部微分算子,將固體力學中的變形協調和力平衡偏微分方程改寫為積分方程,無需對固體損傷破壞時的不連續位移場進行特別考慮,能夠處理裂縫起裂、擴展、分叉等問題[24-26].然而,該方法難以與經典材料本構關系直接對應[27]、存在波色散[28]和零能模式[29]等缺點,在數值上也與有限元方法不兼容,因此這里不做進一步討論.

相場斷裂模型根植于Francfort 和Marigo[22]提出的線彈性斷裂能量變分原理,即實際裂縫路徑使得固體系統的總能量最小,完美地解決了Griffith 能量方法需要預設裂縫路徑的難題[30].在理論上,該模型借鑒了計算機圖像分割中Mumford-Shah 泛函[31]的橢圓型正則化方法[32](因此,該模型也被稱為AT2相場模型[33]),將尖銳裂縫正則化為具有一定尺度的裂縫帶,并引入在[0,1]范圍內連續分布的裂縫相場變量及其梯度描述裂縫狀態.由能量變分原理給出的裂縫相場控制方程在形式上與凝聚態物理中的Ginzberg-Landau 相變方程[34]、多相流體力學中的Allen-Cahn 擴散界面方程[35]以及計算機圖像分割中的Ambrosio-Tortorelli 方程[32]均有異曲同工之處.在數值上,該模型可以非常方便地通過有限元、無網格、物質點甚至近場動力學等多種數值方法加以實現.可以嚴格證明[36]:當裂縫尺度趨近于零時,上述相場斷裂模型Γ?收斂于Griffith 線彈性斷裂力學.2010 年以來,德國的Miehe 等[37-38]引入經典損傷力學中的能量拉?壓分解[39]、歷史最大損傷能釋放率[40]等概念,在熱力學理論框架內重新推導了上述脆性斷裂相場模型,并將其應用于固體結構的靜力、動態和多場耦合破壞分析.至此以后,相場斷裂模型迅速得到了學術界的廣泛重視,研究人員開展了大量的跟蹤研究和應用工作.例如,法國的Marigo團隊[41]提出了脆性斷裂相場模型AT1,解決了AT2模型缺乏線彈性段的問題;在國內,近年來清華大學莊茁團隊[42-44]將相場斷裂模型AT2 應用于動態熱?力耦合條件下的絕熱剪切帶分析,中國科學技術大學李良彬團隊[45-46]采用AT2 模型開展了高速沖擊下脆性和超彈性材料的動態破壞分析.

令人遺憾的是,上述相場斷裂模型AT1/2 均僅適用于脆性破壞,無法應用于工程中普遍存在的準脆性材料和結構.更糟糕的是,模型的分析結果存在嚴重的裂縫尺度敏感性問題:裂縫尺度越小,結構的極限承載力(峰值載荷)越大.因此,若將裂縫尺度視為數值參數以保證模型的Γ?收斂特性,則此類模型與Griffith 線彈性斷裂力學類似,同樣無法描述完好或僅有弱奇異性固體的裂縫起裂.迫不得已,研究人員只好放棄相場斷裂模型的Γ?收斂性質,并類似非局部損傷模型[47-48],將裂縫尺度視為材料屬性[33].然而,最新研究表明:上述處理方法仍然無法完全解決裂縫起裂問題[49],這一致命缺陷也被相場斷裂模型的創建者所證實[50].

為了解決上述問題,荷蘭的de Borst 團隊[51-53]以及意大利的Iurlano 團隊[54-55]等分別嘗試從內聚裂縫模型出發建立相場模型.然而,前者僅適用于界面裂縫而無法考慮任意裂縫擴展路徑,后者則同樣存在裂縫尺度敏感性問題.

2017 年以來,筆者將斷裂力學和損傷力學有機結合,引入兩個參數化的特征函數即裂縫幾何函數和能量退化函數,分別表征固體破壞時的裂縫相場分布特征以及開裂固體應變能的退化規律,基于熱力學基本原理建立了同時適用于脆性斷裂和準脆性破壞的統一相場理論[56].從該理論出發,不僅脆性斷裂相場模型AT1/2 可以作為其特例直接給出,還在國際上首次提出了一類相場正則化內聚裂縫模型PF-CZM[57-58].理論分析證明:當裂縫尺度趨近于零時,該模型收斂為一類混合型破壞的內聚裂縫模型[59-60],能夠準確預測線性軟化以及指數軟化、雙曲型軟化、混凝土科內列森軟化[61]等常用的非線性軟化曲線.由于同時考慮了基于強度的裂縫起裂準則、基于能量的裂縫擴展準則以及基于變分原理的裂縫擴展方向判據,該模型僅需少量標準材料參數即可定量預測工程結構的損傷破壞行為;特別是,當裂縫尺度小于某一上限時,其取值對脆性斷裂和準脆性破壞的分析結果幾乎不產生影響.歸功于上述優點,統一相場理論提出后迅速得到了國內外學者的廣泛認可和跟蹤研究,并被成功應用于靜力、動力和多場耦合條件下混凝土[62]、巖石[63]、冰[64]、玻璃[65]、PMMA[58]、復合材料[66-67]、橡膠[68-69]等多種固體材料和結構的損傷破壞分析.

前已述及,相場模型可以非常方便地通過多種常用數值方法加以實現.然而,由于裂縫相場梯度的引入,為保證計算結果的數值精度,裂縫帶附近的空間離散必須具有足夠的數值分辨率,往往導致系統總自由度數目較大,對非線性方程數值求解算法提出了很高要求.已有研究大多遵循文獻[23,30]的建議,采用裂縫相場?位移場交錯迭代求解器.該求解器具有優異的數值健壯性,但收斂速度極慢(對于某些裂縫快速擴展的關鍵步,往往需要成百甚至上千次全局迭代才能收斂).因此,已有固體結構損傷破壞的相場分析往往局限于簡單二維算例,復雜三維問題的相場模擬極為少見[70].近年來,筆者采用整體隱式BFGS 擬牛頓迭代算法求解裂縫相場?位移場耦合方程,計算效率較交錯迭代求解器大幅度提升[71-73],成功地在個人計算機上開展了復雜三維問題的相場分析[74].結合并行計算技術和大型超算平臺,有望應用于實際工程結構的損傷破壞分析.

考慮到統一相場理論的相關研究工作分散發表于若干外文期刊,本工作將對該理論及其數值算法和應用算例加以系統總結、整理和介紹,以期為工程結構的損傷破壞分析提供一種可行的方法,并對該方向今后的研究工作進行展望.

1 統一相場理論

不失一般性,考慮如圖1 所示內嵌裂縫/界面的固體? ?Rndim(ndim=1,2,3),其外邊界和外法向矢量分別記為?? ?Rndim?1和n.雖然統一相場理論同樣適用于有限變形[75],這里僅考慮小應變情況,故采用位移場u(x)和線性應變場?(x):=?su(x)描述固體的變形狀態,這里x為空間坐標,?s(·)表示對稱梯度算子.為保證邊值問題的適定性,外邊界?? 分成互補的兩部分??u和??t,并分別施加給定的位移邊界u?(x)和力邊界t?(x).

圖1 固體內嵌裂縫/界面及其幾何正則化[60]Fig.1 A cracking solid and its geometric regularization[60]

1.1 裂縫/界面的相場正則化

相場模型將如圖1 所示的尖銳裂縫/界面S 彌散為有限尺度b>0 的裂縫帶B ??,并引入裂縫相場d(x):B →[0,1]描述其狀態;裂縫帶B 的外邊界記為?B,其外法向矢量表示為nB.需要指出的是,在固體損傷破壞過程中裂縫帶B 并非預先設定或保持固定,而是遵循自身特定的本構行為進行擴展、演化.類似于位移場,裂縫相場也可以施加某種強迫邊界條件.例如,對于彈性區域有d(x)=0;對于預制的初始裂縫有d(x)=1 等.

根據圖像分割理論[36],裂縫面積AS可以表示為如下體積分[37]

式(1)中,左側第一個大括號代表尖銳裂縫,右側大括號代表正則化裂縫.意大利數學家Ambrosio 和Tortorelli[32]在對圖像分割Mumford-Shah 泛函[31]進行正則化時,最早給出了形如式(2)的裂縫面積密度泛函,其裂縫幾何函數取為α(d)=d2.可以嚴格證明[36]:當裂縫尺度b趨近于零時,式(1)給出的裂縫面積AS滿足Γ ?收斂即AS=limb→0Ad(d).裂縫面積密度γ(d;?d)是裂縫相場d及其梯度?d的如下函數

式中b>0 為裂縫尺度,裂縫幾何函數α(d)∈[0,1]是裂縫相場d(x)的遞增函數,并滿足如下條件[56]

為了保證固體破壞時局部能量耗散的遞增性質,這里引入了條件(3)1(公式號的下標1 代表公式中的第1 式;以此類推),在物理學領域,美國Karma教授團隊最早采用雙勢阱函數α(d)=d2(1?d)2并基于Ginzberg-Landau 相變方程[34]建立脆性斷裂相場模型[76-77].雙勢阱函數滿足性質α(0)=α(1)=0,即材料處于中間過渡狀態時的勢函數比純相時更大,這一特點與材料相變、多相流或結構拓撲優化等問題[78]較為契合.雖然雙勢阱函數同樣滿足Γ ?收斂特性,但該函數在[0,1]之間不具有單調遞增性,也無法區分d=0 和d=1 這兩種截然不同的裂縫狀態.因此,若將其作為裂縫幾何函數,將會導致裂縫面積和能量耗散隨裂縫擴展先增大后降低,這與固體斷裂時Griffith 能量原理(即應變能減小、裂縫表面能增大、最終總能量達到最小值)相矛盾.本文僅討論滿足條件(3)的裂縫幾何函數.

受裂縫影響,固體Helmholtz 自由能ψ 發生退化,并表示為如下一般形式

表征裂縫相場影響的能量退化函數ω(d)∈[0,1]應滿足以下條件[37]

特別指出,條件(5)4保證固體損傷破壞過程中裂縫帶的寬度不會持續增大.

1.2 裂縫相場演化的能量原理

熱力學第一定律和第二定律要求:等溫絕熱條件下,任意容許的變形過程必須滿足如下能量耗散不等式

式中,應力張量σ與應變張量?功共軛;為時間導數.代入Helmholtz 自由能函數(4),并考慮應變率的任意性,可以給出固體的應力?應變關系并定義裂縫相場的能量釋放率

其中,:=?ψ/?ω 為有效能量釋放率,ω′(d):=?ω/?d為能量退化函數的導數.相應的,能量耗散不等式(6)簡化為

式(8)表征了固體損傷破壞過程中的應變能釋放.

由此,可以給出如下裂縫相場演化的能量原理[56-57].

(1)裂縫相場演化的不可逆性:在封閉系統中,裂縫相場演化滿足不可逆性,即

因此,裂縫相場的變分δd≥0 非負(這里不考慮裂縫愈合).

(2)裂縫相場演化的能量準則:裂縫相場的準靜態演化滿足如下關系式

式中,Gf為固體的斷裂能常數,第一個大括號代表應變能釋放,第二個大括號代表裂縫表面能.

利用高斯?斯托克斯散度定理并根據裂縫相場變分的非負特性δd≥0,可以給出如下裂縫相場演化的能量準則

裂縫面積密度函數γ 的變分導數δdγ 表示為

式中,?d=?· ?d為裂縫相場的拉普拉斯運算符;α′(d):=?α/?d為裂縫幾何函數的導數.此外,還有如下邊界條件

其中,通量q與裂縫相場梯度?d之間滿足線性本構關系,比例系數與裂縫尺度b>0 成正比.

(3)裂縫相場演化的能量守恒:任意時刻,裂縫相場演化還需要滿足能量守恒原理,即

當g(Y,d)<0 時,裂縫相場發生卸載即=0;當裂縫相場發生進一步演化即>0 時,必須滿足能量準則g(Y,d)=0 和邊界條件q·nB=0.

式(9)、式(10)和式(13)可以統一寫成

上述裂縫相場演化的能量原理將斷裂力學和經典損傷力學緊密聯系,構建了兩種理論之間的橋梁.

1.3 裂縫相場?位移場控制方程

可以將裂縫相場演化能量準則(10)和邊界條件(12)改寫為如下形式

即裂縫相場通量q的散度與相場源Q(?,d)平衡

同時,位移場u應滿足經典的力平衡方程

式中,b?為體積均布力.

可以看出,上述偏微分方程在形式上一致且相互耦合,構成了統一相場理論的裂縫相場?位移場耦合控制方程.

可以證明[59,79],上述裂縫相場?位移場耦合控制方程與系統總能量最小原理[23,30]等效,即

其中,系統總能量E 表示為

式中,第一、二、三個大括號分別代表固體應變能、裂縫表面能、系統總勢能.換句話說,在所有可能的裂縫路徑中,實際裂縫使得系統總能量最小.與經典Griffith 能量理論不同的是,這里不要求預設或已知裂縫路徑.

統一相場理論同樣適用于固體結構的動態損傷破壞分析.此時,平衡方程(15c)1需考慮慣性力的影響,即[59,80]

1.4 應力?應變本構關系

控制方程(15)還需補充Helmholtz 自由能ψ的具體表達式.顯然,其定義與分析對象和分析目的直接相關.特別是,對于脆性和準脆性固體,與Helmholtz 自由能對應的應力?應變本構關系和裂縫相場能量釋放率(7)必須考慮材料受拉?受壓力學行為的不同特性.為此,類似于損傷力學的思路[39,82],研究人員分別提出了應變球?偏分解[83-84]、應變拉?壓譜分解[38]、有效應力拉?壓能量正交分解[60,85]等多種方法,對材料Helmholtz 自由能ψ 進行分解,詳見綜述[59]中的討論.

為簡單起見,這里采用各向同性Helmholtz 自由能定義

有效應力張量由如下線彈性關系給出

式中,E0和C0分別為材料的彈性剛度張量和柔度張量.于是,由式(7)得到

有效能量釋放率表示為

其中,E0和ν0分別為材料楊氏模量和泊松比.如圖2 所示,上述應力?應變本構關系將給出完全對稱的受拉和受壓力學行為,這與脆性和準脆性材料的實際力學特性不符.

圖2 不同有效能量釋放率定義給出的強度包絡線:平面應力狀態Fig.2 Strength envelopes given by different effective energy release rates:Plane stress state

為考慮拉、壓應力狀態下的不同力學行為,以單軸受拉強度為基準,將有效能量釋放率重新定義為

需要說明的是,上述修正后的應力?應變關系僅適用于裂縫處于受拉張開狀態的情況,而無法考慮反復加載條件下裂縫閉合引起的剛度恢復.對于后者,可以采用有效應力張量的能量正交分解法[60,85]加以描述;如果需要進一步考慮裂縫擴展對材料受壓力學行為的影響,類似筆者在雙標量損傷模型方面的研究工作[85-86],可以引入另一個裂縫相場變量,此處略過不表.

1.5 應變局部化分析

筆者之前的研究工作[87-89]表明:非彈性材料發生應變局部化并最終形成尖銳裂縫應滿足麥克斯韋變形協調條件,即

這里[|u|]為裂縫的位移跳躍.上述麥克斯韋變形協調條件如圖3 所示.

圖3 統一相場理論的應變局部化[60]Fig.3 Strain localization of the unified phase-field theory[60]

對于式(21)1給出的應力?應變關系,非彈性應變?in表示為

式中引入了單調遞增的函數φ(d)

由式(5)可知,函數φ(d)需滿足如下條件

當裂縫尺度趨近于零即b→0 時,裂縫帶B 內的應力幾乎均勻分布.于是,由式(25)和式(26)給出

相應的內聚力確定為

其中,局部化損傷變量?(d)表示為

二階對稱張量G0:=nS·E0·nS為裂縫的初始剛度.

上述分析表明:裂縫尺度b→0 時,統一相場理論收斂為一類同時考慮拉伸?剪切混合破壞行為的內聚裂縫模型[90-91].這一事實,一方面展示了統一相場理論的合理性;另一方面,也為后文確定裂縫幾何函數和能量退化函數提供了理論依據.

2 脆性斷裂和準脆性破壞相場模型

上述統一相場理論中,裂縫幾何函數α(d)和能量退化函數ω(d)這兩個特征函數決定了其分析結果.本節針對脆性和準脆性這兩種典型的固體結構破壞行為,確定相應的特征函數,并由此建立脆性斷裂和準脆性破壞相場模型.

2.1 解析解:一般情況

為了確定統一相場理論的特征函數,分析其單軸受拉情況下的解析解.假設長度足夠長的桿x∈[?L,L],其左、右兩端作用單調遞增且大小相等、方向相反的受拉位移.不考慮桿的均布體積力.

在開始階段,應變場沿桿長均勻分布,直至桿內應力σ 達到如下峰值σc之前[56,59]根據裂縫幾何函數α(d)和能量退化函數ω(d)的性質,峰值應力對應的裂縫相場dc∈(0,1)和系數βc按以下三種情況給出[59]

即最終破壞時,裂縫相場分布僅與函數α(d)有關.

2.2 解析解:特殊情況

為了更好地說明上述解析解,本節考慮如下滿足條件(3)和(28)的最簡函數形式[56]

式中,P(d)=1 +a2d+a3d2+··· 為多項式;系數ξ ∈[0,2],a1>0,a2,a3,··· 以及指數m>0 均為待確定的模型參數.需要指出的是,為了保證裂縫幾何函數(38a)滿足α(d)∈[0,1]單調性遞增條件,參數ξ應在[0,2]范圍內取值,如圖4(a)所示.

圖4 參數ξ 取值對裂縫幾何函數及最終裂縫相場分布的影響Fig.4 Effects of the parameter ξ on the crack geometric function and the ultimate distribution of the crack phase-field

對于裂縫幾何函數(38a),由式(37a)可以給出最終破壞時裂縫相場分布du(x)和裂縫帶半寬度Du的解析表達式[56,59].特別對于典型參數ξ,其最終裂縫相場分布du(x)如圖4(b)所示.可見,隨著參數ξ ∈[0,2]增大,裂縫帶半寬度Du(ξ)逐漸減小.

此時,由特征函數(38a)和(38b)的性質可以給出

峰值應力對應的裂縫相場dc與參數ξ ∈[0,2]有關.

根據參數a1>0 的取值,考慮以下兩種情況:

(1)參數a1與裂縫尺度b無關(記為a1=c1):此時,峰值應力σc與裂縫尺度b的平方根成反比,即

當裂縫尺度b→0 時,峰值應力σc無窮大,這一特性與經典線彈性斷裂力學一致.相應的,統一相場理論退化為一類脆性斷裂相場模型,見2.2.1 節.

與經典Griffith 或Irwin 線彈性斷裂力學類似,脆性斷裂相場模型難以考慮裂縫起裂[22].為了解決這一難題,峰值應力必須為有限大小.最常用的方法是假定裂縫尺度參數b按下式確定[33,41,84]

式中,ft為材料的破壞強度(一般取為單軸抗拉強度);lch:=為表征固體材料脆性程度的Griffith 或Irwin 特征長度:其值越小,材料越脆.對于脆性材料,這一方法雖然能夠描述裂縫起裂,但存在以下主要缺點[49,58]:①裂縫尺度b為材料屬性而并非數值參數,導致裂縫面積正則化公式(1)和整個模型失去Γ收斂性質;②對于某些脆性材料,由此確定的裂縫尺度b過大(如砂漿b≈50 mm),會導致裂縫面積計算精度不足、裂縫帶過寬、破壞模式失真等問題,從而高估裂縫能量耗散、結構承載力和延性等關鍵指標,詳見文獻[50]中的討論.

(2)參數a1與裂縫尺度b成反比:為了保證峰值應力σc為有限值(即材料破壞強度ft)且與裂縫尺度b無關,由式(39)有

換句話說,a1并非材料屬性,而是與裂縫尺度b成反比的數值參數.此時,所給出的模型仍然滿足Γ 收斂性質:其值越小,裂縫面積(式(1))的精度越高,同時能量退化函數ω(d)越陡峭(如圖5 所示),裂縫對周邊固體的影響范圍也越小,正則化裂縫越接近真實裂縫.

圖5 參數a1>0 對能量退化函數ω(d)的影響[58]Fig.5 Effects of the parameter a1on the energetic degradation function[58]

此種情況下,當其他模型參數滿足一定條件時,統一相場理論可以退化為一類可以描述應變軟化行為的準脆性破壞相場模型,見2.2.2 節所述.

2.2.1 脆性斷裂

對于a1=c1為常數的脆性斷裂相場模型,可以考慮以下參數

即[23,30]

根據裂縫幾何函數(38a)的性質,有以下兩種情況.

(1)參數ξ=0,即α(d)=d2和cα=2.峰值應力對應的裂縫相場由式(33)確定

上述AT2 模型由文獻[23,30]提出并經Miehe 等[37]完善,是最早應用于脆性斷裂分析的相場模型.

(2)ξ ∈(0,2]:此時,由式(34)可以得到

常用模型有以下兩種

這兩種模型給出的應力?應變曲線均為理想彈脆性,且裂縫相場為有限分布,與實際情況較為相符,因此在固體脆性斷裂中的應用較為廣泛[33].

圖6 比較了上述不同脆性斷裂相場給出的單軸受拉結果.可見,AT2 模型給出的應力?應變關系缺乏線彈性階段,這與線彈性斷裂力學不一致,而且最終破壞時裂縫相場的分布范圍為[?∞,+∞],也與實際情況不符.按式(42)確定裂縫尺度b后,AT1 模型和WN 模型雖然能夠地描述脆性材料的起裂和裂縫擴展,但無法考慮準脆性材料的損傷破壞行為.

圖6 不同脆性斷裂相場模型給出的應力?應變曲線Fig.6 Stress-strain curves given by different phase-field models for brittle fracture

2.2.2 準脆性破壞

為了更好地解決裂縫起裂和準脆性破壞問題,裂縫幾何函數(38a)的參數取為ξ ∈(0,2],以保證裂縫相場為有限分布.

此時,有α′(0)>0 且ω′(0)<0.相應的,應力達到峰值σc之前,應變場沿桿全長均勻分布,裂縫相場為零,應力?應變關系為線彈性.峰值應力σc和對應的裂縫相場dc表示為

為了保證峰值應力σc為有限值(即材料破壞強度ft)且與裂縫尺度b無關,根據式(43)有

式中,裂縫尺度b為數值參數,其取值越小,裂縫面積(1)的精度越高;其取值小于某一上限值后,對計算結果幾乎不產生影響.

隨著施加的位移進一步增大,桿內會出現應變局部化和裂縫相場非均勻分布,應力?名義裂縫位移跳躍曲線σ(w)由式(36)給出,如圖7 所示.

圖7 準脆性破壞相場模型給出的應力?名義位移跳躍曲線[56]Fig.7 Stress-displacement(apparent)jump curves given by different phase-field models for quasi-brittle fracture[56]

軟化曲線σ(w)的初始斜率k0和極限張開位移wc確定為[56]

此外,裂縫帶半寬度的初始值D0表示為

其中,P(1)=1+a2+a3+··· 為多項式P(d)在d=1時的取值.從式(37)和式(53)可以看出,裂縫帶寬度與裂縫尺度b成正比.

由圖7 和式(52)可知:m<2 給出的極限裂縫張開位移為零,即軟化曲線出現回跳,這與準脆性破壞特征不符.因此,能量退化函數(38b)中的指數參數m的取值范圍為

具體而言,m>2 給出的極限裂縫張開位移wc無窮大,而m=2 給出的wc為有限值.

上述情況下,統一相場理論可以給出合理的軟化曲線,其破壞強度ft、初始斜率k0和斷裂能Gf均為有限值,表現出典型的準脆性破壞特征.

2.3 相場正則化內聚裂縫模型

上述結果表明,為了合理反映損傷破壞過程中的應變軟化段,應采用如下裂縫幾何函數和能量退化函數

式中,參數ξ ∈(0,2]以及m≥2,且a1>0,a2和a3確定為[56]

此時,式(36)給出的應力?名義位移跳躍曲線σ(w)與裂縫尺度b無關,而僅取決于彈性模量E0、破壞強度ft、斷裂能Gf以及軟化段參數(如初始斜率k0、裂縫極限位移wc)等材料屬性.由此,對于裂縫幾何函數和能量退化函數(55),統一相場理論給出了一類相場正則化內聚裂縫模型.例如,文獻[64]選用ξ=1/2即PF1/2-CZM,并應用于飛機除冰優化設計;文獻[94-96]選用ξ=1 即PF1-CZM,并應用于混凝土和巖石等準脆性材料結構.然而,這些模型均無法應用于線性軟化.

理論上,對于所有的參數ξ ∈(0,2]上述相場正則化內聚裂縫模型都成立.然而,不同參數ξ 對應的模型適用范圍有所不同.特別是,裂縫相場不可逆≥0 要求:裂縫帶半寬度初始值D0不能超過其最終值Du,即[56]

基于上述考慮,筆者建議采用如下裂縫幾何函數

參數a1>0,a2和a3仍由關系式(56)給出,即

這里βk和βw分別表示為

對于線性軟化,有βk=βw=1.此時,式(36)給出的典型軟化曲線如圖8 所示.可以看出,盡管P(d)只采用了二次多項式,線性軟化、指數軟化和雙曲線軟化等均得到了高精度的描述,混凝土科內列森軟化曲線[61]也相當吻合;其他如雙線性等軟化曲線也可以類似得到,見文獻[62].因此,參數ξ=2 對應的相場正則化內聚裂縫模型PF2-CZM 最為完備,且能夠廣泛適用于各類常用軟化曲線,在混凝土[62]、玻璃[65]、PMMA[58]、復合材料[66-67]、橡膠[68-69]、巖石[63]等多種固體材料結構中得到成功應用.

圖8 相場正則化內聚裂縫模型PF2-CZM 給出的軟化曲線[56]Fig.8 Softening curves given by the phase-field regularized cohesive zone model(PF2-CZM)[56]

圖8 相場正則化內聚裂縫模型PF2-CZM 給出的軟化曲線[56](續)Fig.8 Softening curves given by the phase-field regularized cohesive zone model(PF2-CZM)[56] (continued)

在計算量允許的條件下,裂縫尺度b的取值越小，裂縫面積式(1)以及相應能量耗散的計算精度越高.極限情況下,當裂縫尺度b→0 時,裂縫帶寬度也趨近于零,統一相場理論退化為1.5 節的內聚裂縫模型.因此,相場正則化內聚裂縫模型不僅適用于準脆性破壞,還可采用線性軟化描述脆性斷裂[58],實現了兩種典型破壞特征的統一反映.

3 數值實現和求解算法

統一相場理論數值實現時,可以將偏微分控制方程(15)轉化為其弱形式:求解滿足邊界條件u(x∈??u)=u?和不可逆條件d˙(x)≥0 的位移場和裂縫相場(u,d),使得下式對于任意容許的位移場和裂縫相場(δu,δd)均成立

與位移邊界條件類似,裂縫相場d(x)也可以施加諸如d(x∈?B)=0 或d(x∈S0)=1 等強制邊界條件,這里S0為初始預制裂縫.

3.1 多場有限元數值實現

式(60)是關于位移場u(x)和裂縫相場d(x)的耦合方程,通常采用多場有限元方法進行數值求解.

圖9 計算區域的有限元空間離散[97]Fig.9 Finite element discretization of the computational domain[97]

在有限元空間離散時,可以將整個計算區域內所有的單元節點均賦予位移和裂縫相場自由度,也可以如圖9 所示將整個計算域分為兩部分[97]以提高計算效率,即:可能出現裂縫的相場子區域Bh以及剩余彈性部分?hBh,對于前者,單元節點同時具有位移自由度和裂縫相場自由度,而后者的單元節點僅有普通位移自由度.已有研究表明[56-58,62,79],為保證計算精度,一般相場子區域Bh內的單元大小一般取h≤b/5.由于固體開裂后會發生劇烈的變形局部化,相場模型的數值實現必須反映裂縫帶內相場變量的梯度變化,這對網格大小提出了要求.Γ ?收斂分析表明[30]:由于有限元空間離散的數值誤差,正則化關系式(1)將高估裂縫面積約h/(cαb).為減少計算量,已有相場模擬大多采用h=b/2,此時將高估裂縫面積16%～25%左右.出于計算精度的考慮,筆者[56]建議采取h≤b/5,相應裂縫面積或能量耗散的誤差控制在6%～10%以內.這一取值也得到了數值模擬結果的驗證[62,79,97]和其他學者[33,50,98]的采納.在實際工程應用中,若需要采用較粗的有限元網格,可通過減小材料斷裂能等方式[70],部分補償有限元空間離散引起的數值誤差.

有限元空間離散后,位移場和應變場、裂縫相場及其梯度分別表示為(在數值實現中,張量均采用Voigt 標記)

式中,fext為標準外力列向量[99].由于裂縫相場演化的不可逆性質,其平衡方程為不等式,需要在數值求解時加以處理.

方程組(62)通常采用增量迭代法進行求解,即將時間步[0,T]劃分為N個子步[tn,tn+1],這里n=0,1,···,N?1.對于時間增量為?t:=tn+1?tn的典型子步[tn,tn+1],tn時刻所有的狀態變量已知,在此基礎上采用迭代方法求解tn+1時刻的狀態變量.

對于考慮慣性力的動力平衡方程,有限元離散后的控制方程表示為

式中:=d2a/dt2為單元節點的加速度列向量,M:=為一致質量矩陣.此時,可以采用隱式算法(如Newmark-β 法、HHT-α 法等),將上述常微分方程轉換為節點位移a的非線性代數方程進行求解,也可以采用顯式動力積分方法,直接計算節點位移.具體步驟詳見文獻[59].

3.2 裂縫相場演化不可逆性

為了處理裂縫相場演化的不可逆性,研究人員提出了增廣拉格朗日[100]、罰函數[101]等多種方法,但前者會引入新的自由度,而后者存在罰剛度選用難題.這里介紹另外兩種較為常用的方法,即約束優化方法和歷史最大變量方法.

(1)在有限元方法和增量求解算法的框架內,相場子區域Bh內每個單元節點J均滿足不可逆條件,即,這里下標n和n+1 分別表示tn和tn+1時刻.相應的,裂縫相場子問題(62)2可改寫為如下約束優化問題[102-104]

并采用減縮空間激活集牛頓方法[105]等約束優化求解器進行求解,此時迭代求解過程僅需考慮等式(64)1.具體過程見開源工具箱PETSC[106]中的非線性求解器SNES.

(2)受損傷力學[40]的啟發,將式(15b)的有效能量釋放率替換為其歷史最大值H,即[37]

上述兩種方法的共同點是將裂縫相場子問題不等式(62)2變為等式,即

3.3 控制方程的數值求解

為了求解非線性代數方程組(66),可以考慮以下三種數值算法:整體牛頓迭代算法、子問題交錯迭代算法以及整體BFGS 擬牛頓迭代算法.

3.3.1 整體牛頓迭代算法

對于牛頓迭代算法,將平衡方程(66)進行整體線性化后給出

子剛度矩陣分別表示為

由于所有被積分函數均連續,可采用標準高斯數值積分方法進行數值計算.

需要指出:由于系統總能量(17)是整體未知量z:={a,}的非外凸函數[23,30],上述整體牛頓迭代算法的收斂性難以保證.為此,研究人員提出了基于裂縫面積的弧長法[53,107]、反常線性搜索[108]、自適應修正牛頓方法[109-110]等,但效果均不十分理想.

對于弱耦合非線性問題,通?？梢院雎苑菍ΨQ項,即

已有研究表明[110],對于具有強非線性和強耦合特征的相場斷裂模型,這種忽略耦合項的修正牛頓迭代算法效果也很差.

3.3.2 子問題交錯迭代算法

雖然系統總能量(17)是整體未知量(u,d)的非外凸函數,但卻分別單獨是位移場u或裂縫相場d的外凸函數[23,30].因此,可以通過交錯迭代算法進行求解式(62).對于某一時間增量步[tn,tn+1],以第k次迭代過程為例,有

上述位移子問題可以通過牛頓迭代方法進行求解,其線性化方程表示為

剛度矩陣Kuu由式(68)給出.

(2)利用上述求得的節點位移a(k),通過相場子問題(66)2求解相場自由度,即

其線性化方程表示為

剛度矩陣Kdd由式(68c)給出.

上述求解過程交替往復,直至最后達到收斂,具體過程詳見[59,79].收斂判據可以采用以下幾種.

(1)不迭代一次過[38]:將相場模型視為裂縫相場?位移場的順序耦合問題,不檢查整個控制方程是否收斂.這種方法雖然得到了廣泛應用[111-113],但必須采取極小的時間步長(通常為10?5～10?7量級)才能保證計算精度;更糟糕的是,這種方法會延緩裂縫擴展,從而可能導致不準確甚至錯誤的數值結果.

(2)基于裂縫相場的收斂準則[23,30,84]:當連續兩次整體迭代后的裂縫相場足夠接近時判斷為收斂,即,這里|·|2為L2范數,容差可取為?=1.0×10?5.

(3)基于能量的收斂準則[108,114]:類似于上述裂縫相場收斂準則,檢查連續兩次整體迭代后的系統總能量(17),當二者足夠接近時即判斷為達到收斂標準.由于系統總能量對裂縫相場不敏感,采用該收斂準則需非常小心,詳見文獻[114].

(4)基于殘差的收斂準則[66,108,110,115]:每次整體迭代結束后檢查平衡方程(66)的殘差r和ˉr判斷是否收斂.這種收斂準則非常便于在商業有限元軟件中實現.

與整體牛頓迭代或修正牛頓迭代算法相比,子問題交錯迭代算法具有很好的數值健壯性,特別是與局部弧長法[116]聯合使用其收斂性還可以進一步提高[79].然而,該算法的收斂速度較慢,某些關鍵增量步甚至需要上千次整體迭代,求解效率無法滿足大規模計算要求.

為了克服子問題交錯迭代算法計算效率低的問題,文獻[104]提出了一種復合算法,即在某一增量步的開始階段采用子問題交錯迭代算法,一旦判斷求解進入收斂半徑后換用整體牛頓迭代算法.然而,這種復合算法的轉換時刻缺乏理論依據,因而只能針對具體問題進行試算,難以實現通用化[56-57],其數值實現也非常復雜.

3.3.3 整體BFGS 擬牛頓迭代算法

由于系統總能量泛函非外凸,整體牛頓迭代或修正牛頓算法的收斂性非常差;子問題交錯迭代算法具有很好的收斂性,但其收斂速度慢、計算效率低.近年來,筆者首次將整體BFGS(Broyden-Fletcher-Goldfarb-Shanno)擬牛頓迭代算法[117]應用于統一相場理論控制方程的求解,在保證算法收斂性的同時,收斂速度和計算效率大幅度提升.

在求解非線性方程組g(z)=0 時,BFGS 擬牛頓迭代算法通過如下遞推公式更新計算剛度矩陣

與整體修正牛頓迭代算法不同,雖然上述初始剛度忽略了耦合項的影響,更新后的剛度矩陣仍然能夠考慮裂縫相場?位移場之間的耦合.此外,為了更加合理地反映耦合項的影響,可以構造更準確的初始剛度矩陣,進一步提升計算效率[143].

在實際計算中,BFGS 擬牛頓迭代算法往往與線性搜索聯合使用[117],并根據需要設置間隔一定迭代次數(通常為5～10 次)才對剛度矩陣進行更新計算.因此,所需迭代次數可能比牛頓迭代算法更多(在后者能夠收斂的前提下),但由于BFGS 算法無需每次迭代都更新剛度矩陣且計算過程較為簡單,計算效率通常反而更高.較之廣泛采用的子問題交錯迭代算法,整體BFGS 算法的求解效率可以提高5～8倍[71-73].

3.4 數值實現平臺

從上述數值算法可以看出,統一相場理論非常便于通過有限元方法加以實現,且其二維和三維數值實現幾乎沒有區別.比較常用的開源和商業有限元計算平臺有:

(1)Matlab:借助Matlab 軟件友好的開發環境和強大的圖形處理功能,文獻[45,119]等實現了脆性斷裂相場模型AT2,并應用于固體的靜力和動態裂縫擴展分析.

(2)Feap:該有限元軟件采用Fortran 語言編程,被廣泛應用于學術研究.文獻[120-121]基于該平臺實現了脆性斷裂相場模型AT2.

(3)Mef90:文獻[70]針對脆性斷裂相場模型AT1/2,采用Fortran90 語言開發了該軟件.特別是,為保證系統總能量取得全局最小值,該軟件采用了回溯搜索算法.

(4)Deal.II:是一種C++編程的開源有限元平臺,支持有限元網格的自適應細分,且通過PETSC 工具包[106]提供了約束優化求解器.文獻[109,122]在該平臺上實現了脆性斷裂相場模型AT2,并通過網格自適應提升計算效率.

(5)FENICS:是一種C++和PYTHON 混合編程的開源有限元平臺,用戶僅需提供系統的總能能量泛函或控制方程的弱形式,非常便于相場模型的初步驗證.利用這一平臺,文獻[104,123]等實現了脆性斷裂相場模型AT1,并進行了靜力和動力裂縫擴展分析.

(6)MOOSE:采用C++編程,是一種應用廣泛的多場耦合分析平臺.基于該平臺，文獻[124-125]針對脆性斷裂相場模型AT2、文獻[126]針對準脆性破壞PF1-CZM 分別進行了數值實現.

(7)Jive:是一種C++編程的開源有限元軟件平臺,非常便于自定義材料和自定義單元二次開發.筆者和合作者[49,59,65]基于該平臺實現了統一相場理論,開發了約束優化求解器和BFGS 擬牛頓算法,并將其應用于靜力、動態以及多場耦合條件下的裂縫擴展分析.

由于大多針對研究目的,上述數值實現平臺通常采用GMSH[127]、PARAVIEW[128]等開源軟件進行前后處理,不便于工程應用.

為此,部分學者基于ABAQUS[127],COMSOL[131]等商業軟件進行二次開發,實現了相場斷裂模型AT2.然而,受軟件接口限制,這些數值實現大多僅適用于脆性斷裂相場模型,且通常采用一次過不迭代或整體牛頓迭代算法,計算效率不高且計算精度難以保證.

最近,筆者通過用戶自定義材料接口UMAT 和用戶自定義單元接口UEL,分別采用子問題交錯迭代算法和整體BFGS 擬牛頓算法,首次在ABAQUS軟件平臺實現了包括AT1/2 和PF-CZM 等在內的各種相場模型[71].相關源代碼已在網上公開,可供研究人員免費下載使用(https://github.com/jianyingwu/pfczm-abaqus).

4 應用算例

本節給出統一相場理論在脆性和準脆性結構損傷破壞分析方面的若干驗證和應用算例.所有數值模擬均采用2.3 節給出的相場正則化內聚裂縫模型(PF2-CZM,對應于ξ=2).對于脆性破壞,采用線性軟化曲線,此時僅PF2-CZM 能嚴格保證裂縫相場演化的不可逆性質;對于混凝土,采用科內列森軟化曲線.所有二維問題假定為平面應力,網格劃分采用四節點雙線性單元;三維問題采用八節點實體單元.所有數值模擬均在小型高性能計算工作站(配置:Intel(R)Core(TM)i7-7700HQ CPU@2.80 GHz 中央處理器、512 GB 內存)上完成.

4.1 脆性破壞分析[58]

首先考慮一組非對稱缺口開孔梁三點彎試驗,該實驗由國際著名斷裂力學學者Ingraffea 教授課題組[132-133]開展,已成為線彈性斷裂力學的標準算例之一.特別是,該系列試驗受預制裂縫與孔洞之間的距離影響較大,因此對分析模型的要求較高.

試件幾何尺寸如圖10 所示,預制缺口的寬度為0.05 in;小孔直徑均為0.5 in (1 in=2.54 cm),試件由同一批PMMA 脆性材料加工制作.文獻[132]給出的材料參數分別為:彈性模量E0=4.75 × 105psi (1 psi=6.895 kPa),泊松比ν0=0.35,斷裂能Gf=1.8 lbf/in (1 lbf/in=0.113 N·m).數值模擬采用b={0.025,0.050} in 兩種裂縫尺度,相場子區域的有限元網格大小h=0.005 in;總單元數約67 萬,計算時間約為20.3 h.

圖10 非對稱缺口開孔梁三點彎試驗:試件尺寸(長度單位:in,1 in=2.54 cm)、載荷和邊界條件[132]Fig.10 Asymmetrically notched three-point bending beams:Specimen dimensions(unit:in,1 in=2.54 cm),loading and boundary conditions[132]

試驗并未測定材料的破壞強度ft,但給出了編號0 的未開孔試件(預制裂縫偏心距e1=6 in,裂縫高度e2=1 in)的載荷?裂縫開口位移曲線,因此先對該試件進行分析.由圖11 可以看出,材料破壞強度取為ft=2500 psi(對應的材料特征長度lch=0.14 in),計算給出的載荷?裂縫開口位移曲線和裂縫路徑均與試驗結果吻合良好.此外,只要保證裂縫面積(1)和能量耗散(17)的計算精度(滿足條件h≤b/5),裂縫尺度對數值模擬結果幾乎沒有影響,這一結論對于固體損傷破壞分析至關重要.

接下來考慮其他3 種開孔試件(每個試件均帶有3 個直徑為0.05 in 的小孔).(1)試件a:e1=6 in,e2=1 in;(2)試件b:e1=5.1 in,e2=1.5 in;(3)試件c:e1=4.75 in,e2=1.5 in.

除試驗給出的材料參數外,數值模擬中材料破壞強度取為ft=2500 psi,采用基于裂縫開口位移(CMOD)的局部弧長法進行加載,以跟蹤梁的損傷破壞全過程.

由于試驗并未給出這些試件的載荷?位移曲線結果,這里僅討論裂縫路徑.試件a 的數值模擬結果與試驗結果的對比如圖12(a)所示.可以看出,二者的吻合結果十分良好,特別是能夠再現底部小孔附近局部應力集中引起裂縫路徑的細微變化.隨著預制裂縫與底部開孔的距離縮短,局部應力集中引起的裂縫路徑擾動更加明顯,但裂縫仍能繞過底部開孔繼續向上擴展直到到達中間第二個開孔,見圖12(b);當二者的距離進一步縮短時,局部應力集中的影響進一步加強,裂縫直接與底部開孔交匯,如圖12(c)所示.同樣,裂縫尺度僅影響裂縫帶的寬度,對數值模擬結果沒有影響.

4.2 準脆性破壞分析[97]

1967 年,位于印度的Koyna 混凝土重力壩遭受里氏6.5 級地震作用,導致壩體迎水面(高度在背水面變坡度處)出現一條裂縫.該混凝土重力壩的幾何尺寸如圖13 所示,裂縫長度為1.93 m.為評估該受損結構的安全性能,這里對大壩在洪水溢流工況下的承載力和損傷破壞過程進行分析.

圖13 Koyna 混凝土重力壩溢流破壞:幾何尺寸、邊界和載荷條件Fig.13 Overflow of Koyna dam:Dimensions(unit:m),loading and boundary conditions

根據文獻[134],數值模擬中混凝土材料參數分別為:密度ρ=2450 kg/m3,楊氏模量E0=2.5 ×104MPa,泊松比ν0=0.20,破壞強度ft=1.0 MPa,斷裂能Gf=100 J/m2,對應的材料特征長度lch=2.5 m.荷載施加分為三個荷載步:重力載荷、三角形分布的滿水庫靜水壓力以及不斷增加的矩形均布載荷(模擬溢流水壓力);前兩步均直接采用載荷控制,第三步采用基于左上角壩頂水平位移的間接位移加載控制,以便跟蹤整個裂縫擴展全過程.

數值模擬考慮b={0.3,0.2} m 兩種裂縫尺度參數,相場子區域的有限元網格大小取為h=b/5;總單元數分別約16.5 萬和36.1 萬,計算時間分別約為3.2 h 和11.8 h.從圖14 可以看出,裂縫帶寬度隨著裂縫尺度b的減小而變窄,但裂縫擴展路徑不受影響:從初始裂縫尖端開始,裂縫沿著大致平行于背水面的方向斜向下擴展.有趣的是,由于重力壩的特點,壩體下部壓應力逐漸變大,因此裂縫達到一定的深度后不再繼續向下擴展,而是在距背水面變坡度最近的位置發生分叉;分叉裂縫發展到一定程度后,由于背水面變坡度位置局部壓應力的影響,裂縫會發生二次分叉等復雜的損傷破壞行為,使得整個結構的總能量最小.

圖14 Koyna 混凝土重力壩溢流破壞:裂縫擴展路徑[97]Fig.14 Overflow of Koyna dam:Predicted crack patterns[97]

當左上角壩頂水平位移達到75 mm 時,大壩結構變形如圖15 所示.這一數值模擬結果與內聚裂縫模型分析給出的結果一致(見美國Bazˇant 教授撰寫的專著[15]封面圖片).然而,內聚裂縫模型需要預設裂縫路徑并進行復雜繁瑣的網格重劃分,而廣泛采用的擴展有限元方法則無法預測上述裂縫分叉現象[134].截然不同的是,無需其他任何處理手段,相場正則化內聚裂縫模型即可自然地預測上述復雜的裂縫擴展引起的結構損傷破壞行為.

圖15 Koyna 混凝土重力壩溢流破壞:左側壩頂水平位移75 mm 時的結構變形圖Fig.15 Overflow of Koyna dam:Deformations when the horizontal displacement of the crest node reaches 75 mm

不同裂縫尺度給出的溢流高度?壩頂水平位移曲線如圖16 所示.可以看出,裂縫尺度參數對計算結果幾乎沒有影響:溢流高度較小(低于5.5 m)時,壩頂水平位移基本呈線性變化;裂縫快速擴展時,結構所能抵抗的溢流高度會出現短時降低;裂縫進入穩定擴展階段,壩頂水平位移隨溢流高度增大而非線性增長,并逐漸趨于某一極限值.

圖16 Koyna 混凝土重力壩溢流破壞:溢流高度?壩頂水平位移[97]Fig.16 Overflow of Koyna dam:Curves of the overflow height-horizontal displacement of the crest[97]

4.3 動力破壞分析[135]

接下來考慮如圖17 所示的內部沖擊載荷作用下厚壁圓柱的動態破碎問題[136-137].

圖17 厚壁圓柱動態破碎問題:幾何尺寸(左)和載荷時程(右).圓柱內外半徑分別為rin=80 mm 和rout=150 mm,平面外厚度為70 mm.沖擊載荷的幅值為p0=400 MPa,參數t0=100μsFig.17 Fragmentation of a thick cylinder:Geometry(left)and the loading history(right).The innner and outer radii are rin=80 mm and rout=150 mm,and the thickness is 70 mm.For the impact pressure,the amplitude is p0=400 MPa and the parameter is t0=100μs

根據文獻[137]中的討論,考慮材料斷裂能Gf服從如下Weibull 分布

式中,平均斷裂能Gf0=20 N/mm;Rand為[0,1]之間的隨機數,指數取為m=2.基于類似的考慮,文獻[95,138]考慮彈性模量的隨機分布.其他材料參數均由文獻[139]給出:密度ρ=7850 kg/m3,楊氏模量E0=2.1×105MPa,泊松比ν0=0.30,破壞強度ft=1000.0 MPa.

數值模擬采用裂縫尺度b=1.0 mm,有限元網格大小統一取為h=0.2 mm;總單元數約162 萬,計算時間約為64.2 h.動力平衡方程求解采用顯式求解器并行計算,模擬給出的破碎演化過程如圖18 所示:在t≈55μs 時刻,內邊緣的裂縫相場幾乎均勻分布,隨后裂縫繼續沿徑向向外邊緣擴展;在t=70μs 時刻,受有限元網格劃分和材料參數隨機性以及兩條主裂縫間壓應力的影響,有的裂縫停止擴展,部分裂縫則繼續擴展直至圓柱外邊緣;裂縫擴展過程中還會出現如圖18 所示的分叉.

為了分析有限元網格大小的影響,裂縫尺度參數同樣取為b=1 mm,但考慮h={0.2,0.15}mm 兩種網格尺寸;總單元數分別約162 萬和288.5 萬,計算時間分別約為64.2 h 和121.6 h.數值模擬給出的結果如圖19,兩種有限元網格尺寸給出的碎片數量均為15(這里,碎片指兩條從內邊緣貫穿至外邊緣裂縫之間的部分),且計算得到的應變能和裂縫表面能也幾乎重合.

接下來考慮裂縫尺度參數的影響.裂縫尺度參數分別取為b={1.0,1.5,2.0}mm,有限元網格大小均取為h=b/5;總單元數分別約162 萬、72.3 萬和40.8萬,計算時間分別約為64.2 h,31.7 h 和15.1 h.如圖20所示,由于材料隨機性的影響,數值模擬給出的裂縫路徑雖然有所差異,但碎片數量同樣均為15 片,而應變能和裂縫表面能則幾乎完全一樣.上述結果表明,即便對于動力損傷破壞問題,相場正則化內聚裂縫模型的計算結果也與裂縫尺度參數無關,而脆性斷裂相場模型AT1/2 完全無法做到這一點[135].

4.4 多場耦合破壞分析[140]

統一相場理論的控制方程是一組裂縫相場?位移場耦合的偏微分方程,因此特別合適于多場耦合(包括力學載荷在內)條件下固體結構的破壞分析.

為說明此點,這里考慮如圖21 所示的金屬氫脆斷裂問題,即氫離子濃度聚集會降低金屬材料的破壞強度和斷裂能,加速結構損傷破壞,這一病害在金屬管道、壓力容器等實際工程中經常出現.

圖18 厚壁圓柱體動態破碎問題:破碎演化過程[135]Fig.18 Fragmentation of a thick cylinder:Cracks evolution process[135]

圖19 厚壁圓柱動態破碎問題:有限元網格大小敏感性分析[135]Fig.19 Fragmentation of a thick cylinder:Mesh dependence analyses[135]

圖20 厚壁圓柱動態破碎問題:裂縫尺度參數敏感性分析[135]Fig.20 Fragmentation of a thick cylinder:Length scale sensitivity[135]

圖21 金屬氫脆腐蝕問題:幾何尺寸(單位:mm; x=25 mm)、邊界條件和載荷Fig.21 Hydrogen induced corrosion pits problem:Geometry(unit:mm; x=25 mm),loading and boundary conditions

為了描述氫離子的影響,引入如下耦合關系[140]

式中,ft0和Gf0分別為氫離子濃度為零時材料的破壞強度和斷裂能;氫離子表面濃度θ 與體積濃度C之間滿足Langmuir-McLean 關系式[141],后者則根據質量擴散方程和Fick 擴散定律給出,具體詳見文獻[140].

表征氫離子濃度對材料破壞強度和斷裂能退化的影響函數φ1(θ)和φ2(θ)根據試驗給出.為簡單起見,這里假定氫離子濃度不影響材料特征長度lch,即

退化函數φ(θ)=1 ?χθ,系數χ 由試驗或第一性原理計算得到.材料參數由文獻[142]給出:彈性模量E0=2.0×105MPa,泊松比ν0=0.3,破壞強度ft0=1778 MPa,斷裂能Gf0=90 N/mm,退化系數χ=0.89.這里考慮兩種氫離子穩態濃度C={0.0,1.0×10?6}對結構損傷破壞的影響.

數值模擬考慮b={0.4,0.2}mm 兩種裂縫尺度參數,相場子區域的有限元網格大小取為h=b/5;總單元數分別約7 萬和29.5 萬,計算時間分別約為4.2 h和20.7 h.從圖22 和23 可以看出,由于假定氫離子濃度不改變材料的特征長度lch,兩種濃度下裂縫擴展路徑和結構破壞模式基本相同,模擬給出的載荷?位移響應也基本相似,但氫離子的存在會顯著降低極限承載力(見圖24).與純力學問題一樣,無論有無氫離子濃度的影響,裂縫尺度參數b僅改變裂縫帶寬度,而不影響裂縫擴展路徑、破壞模式以及載荷?位移響應等結構損傷破壞行為.

圖22 金屬氫脆腐蝕問題:裂縫尺度b 對裂縫擴展路徑的影響(氫離子濃度C?=0.0)[140]Fig.22 Hydrogen induced corrosion pits problem:Length scale analysis for the hydrogen concentration C?=0.0[140]

圖23 金屬氫脆腐蝕問題:裂縫尺度b 對裂縫擴展路徑的影響(氫離子濃度C?=1.0×10?6)[140]Fig.23 Hydrogen induced corrosion pits problem:Length scale analysis for the hydrogen concentration C?=1.0×10?6[140]

圖24 金屬氫脆腐蝕問題:不同裂縫尺度參數和有限元網格大小對載荷?位移曲線的影響[140]Fig.24 Hydrogen induced corrosion pits problem:Effect of the length scale parameter and mesh size on the load-displacement curves[140]

4.5 三維破壞分析[143]

統一相場理論不受空間維度限制、無需顯式表征裂縫面積并跟蹤裂縫路徑,因此特別適用于非光滑、非規則裂縫擴展引起的三維結構損傷破壞問題.

為此,考慮如圖25 所示I+III 混合型破壞的斜切口三點受彎梁問題[144-145].試件總長度260 mm,跨度240 mm,高度60 mm,平面外厚度10 mm,45?豎向斜切口的寬度和高度分別為2 mm 和20 mm.為了觀測裂縫擴展過程,試驗采用透明PMMA 脆性材料制作.試驗發現,隨著加載位移的增大,開始平行于斜向切口的裂縫會慢慢轉動,直至裂縫角度與梁平面垂直后繼續豎直向上擴展.這一過程中,梁的破壞模式從平面外III 型破壞逐漸轉為平面內I 型.

圖25 斜切口三點受彎梁I+III 混合型破壞問題:幾何尺寸、邊界和載荷條件Fig.25 Mixed mode I+III failure of a skewly notched beam:Geometry,loading and boundary conditions

數值模擬采用文獻[146-147]給出的材料參數:彈性模量E0=2800 MPa,泊松比ν0=0.38,破壞強度ft=20 MPa 以及斷裂能Gf=0.5 N/mm,相應材料特征長度lch=3.5 mm.裂縫尺度取為b=1.0 mm,相場子區域的有限元網格大小h=0.2 mm;總單元數約135 萬,計算時間約為42.6 h.

模擬給出的裂縫相場分布等值面圖以及不同高度截面的裂縫相場云圖分別如圖26 和圖27 所示.可以看出,相場正則化內聚裂縫模型能夠再現試驗觀測到的非光滑裂縫面從斜切口處逐漸向梁頂部發生轉動、然后豎直向上擴展這一復雜的損傷破壞過程.

圖26 斜切口三點受彎梁I+III 混合型破壞問題:裂縫模式對比[143]Fig.26 Mixed mode I+III failure of a skewly notched beam:Geometry,loading and boundary conditions[143]

圖27 斜切口三點受彎梁I+III 混合型破壞問題:不同截面(以斜切口處為高度基準)的裂縫相場云圖[143]Fig.27 Mixed mode I+III failure of a skewly notched beam:Numerical damage profiles at various heights above the notch[143]

圖28 進一步對比了試驗觀測和數值模擬給出的梁側面裂縫路徑.可以看出,數值模擬給出的裂縫路徑完全落在試驗值范圍內[145].

圖28 斜切口三點受彎梁I+III 混合型破壞問題:裂縫路徑與試驗結果對比[143]Fig.28 Mixed mode I+III failure of a skewly notched beam:Comparison of crack patterns[143]

圖28 斜切口三點受彎梁I+III 混合型破壞問題:裂縫路徑與試驗結果對比[143](續)Fig.28 Mixed mode I+III failure of a skewly notched beam:Comparison of crack patterns[143] (continued)

數值模擬給出的載荷?位移曲線如圖29 所示.正如所預料的情況,該梁表現出明顯的彈脆性特征:峰值載荷之前梁幾乎為線彈性;伴隨著裂縫快速擴展,載荷陡降至100 kN 左右,裂縫擴展進入穩定階段,承載力穩步下降,直至最終裂縫擴展到梁頂部結構完全破壞.上述結果表明,相場正則化內聚裂縫模型不僅適用于準脆性材料,同樣適用于脆性材料.

圖29 斜切口三點受彎梁I+III 混合型破壞問題:載荷?位移曲線[143]Fig.29 Mixed mode I+III failure of a skewly notched beam:Load-displacement curves[143]

5 結論與展望

聚焦于固體結構損傷破壞這一世紀難題,本文重點介紹了筆者近年來提出的統一相場理論、算法和若干應用算例.

統一相場理論從裂縫的幾何正則化和熱力學基本原理出發,給出了裂縫相場演化的不可逆性、能量準則和能量守恒條件,由此建立了固體損傷破壞分析的裂縫相場?位移場耦合控制方程.該理論將斷裂力學和經典損傷力學有機結合,同時考慮了基于強度的裂縫起裂準則、基于能量的裂縫擴展準則以及基于變分原理的裂縫擴展方向判據,為固體材料和結構的損傷破壞分析提供了新的理論框架.應變局部化分析表明:裂縫尺度趨近于零時,該理論Γ ?收斂為一類混合型的內聚裂縫模型.

將該理論應用于單軸受拉問題,給出了相應的解析解.特別是,針對一組參數化的裂縫幾何函數和能量退化函數,分別得到了適用于脆性斷裂和準脆性破壞的兩類相場模型:前者的宏觀響應嚴重依賴于裂縫尺度參數,而后者則與其無關.在此基礎上,進一步提出了同時適用于線性軟化和非線性軟化的相場正則化內聚裂縫模型,實現了脆性斷裂和準脆性破壞的統一反映.

統一相場理論的裂縫相場?位移場耦合控制方程非常方便通過有限元等數值方法加以實現.相應單元節點的自由度除節點位移外,還包括裂縫相場自由度.有限元空間離散后得到的非線性代數方程雖然可以采用整體牛頓迭代、子問題交錯迭代以及整體BFGS 擬牛頓迭代等算法進行求解,但其性能差異明顯:雖然整體牛頓迭代算法在收斂半徑內的收斂速度較快,但由于系統能量泛函非外凸,其數值健壯性較差;子問題交錯迭代算法雖然具有優異的收斂性,但收斂速度極慢;整體BFGS 擬牛頓迭代方法同時具有收斂性好、收斂速度快的優點,是目前求解裂縫相場?位移場耦合控制方程最合適的算法.

數值算例表明,統一相場理論、特別是相場正則化內聚裂縫模型僅需少量材料力學參數(彈性模量、泊松比、抗拉強度、斷裂能、軟化曲線類型等),即可直接反映靜力、動力和多物理場耦合等條件下裂縫萌生、擴展、分叉等過程,還能夠有效解決三維非光滑裂縫、I/II 或I/III 型破壞模式轉換等復雜問題.與其他模型/方法相比,相場正則化內聚裂縫模型消除了經典損傷力學等連續方法數值結果的網格敏感性,還消除了非局部損傷模型錯誤預測裂縫起裂、裂縫帶偽加載等問題;同時,克服了傳統內聚裂縫模型網格重劃分復雜繁瑣、彈性罰剛度難以選取等缺點,而且還解決了裂縫幾何表征、裂縫路徑跟蹤、裂縫分叉準則等非連續方法(如內嵌裂縫或擴展有限元等)長期存在的難題;此外,相場正則化內聚裂縫不僅解決了脆性斷裂相場模型存在的裂縫尺度敏感性問題,還可廣泛適用于線性軟化、指數軟化、混凝土科內列森軟化等準脆性破壞.歸功于上述優點,統一相場理論提出后迅速被國內外學者應用于混凝土、巖石、冰、玻璃、PMMA、復合材料、橡膠等多種固體材料和結構的損傷破壞分析.

盡管統一相場理論獲得了較為成功的應用,但在以下方面仍然存在亟需進一步開展的研究工作.

(1)考慮剪切滑移影響的相場正則化內聚裂縫模型.通過修正有效能量釋放率,目前相場正則化內聚裂縫模型可以近似描述裂縫剪應力的影響.然而,單一裂縫相場變量僅能考慮I 型、II 型或以某種破壞形式為主的混合型破壞,尚無法解決巖石力學或地殼斷層分析中復雜的混合型破壞問題.因此,從固體損傷破壞的物理機理出發,更加合理地考慮描述裂縫的法向?切向耦合行為,進一步拓展統一相場理論的邊界,是值得深入研究的課題.

(2)考慮裂縫擴展微慣性和微阻尼的相場動力模型.對于動態破壞問題,裂縫相場與位移場之間的耦合遠較靜力情況復雜.現有動力相場模型大多僅考慮與位移加速度相關的宏觀慣性力,并假定宏觀慣性力與裂縫相場無關,而忽略了結構動態破壞過程中材料應變率效應的影響.文獻[81,148-149]雖然考慮了微慣性和微阻尼對裂縫相場演化的影響,但其參數缺乏物理機理,導致模型模擬結果受參數取值影響極大而失去預測能力.因此,如何從固體動態破壞的物理機理出發,合理地描述裂縫動態擴展過程中材料的應變率效應,值得進一步加以探討.

(3)考慮亞臨界狀態和時間效應的相場疲勞?徐變模型.正常使用條件下工程結構通常處于亞臨界狀態,其應力水平低于材料強度.然而,疲勞載荷或長期載荷作用下,工程結構仍然會出現裂縫起裂和裂縫擴展,嚴重者甚至引發結構災變破壞.工程結構的時變損傷破壞分析一直也是結構分析的重要難點.因此十分有必要考慮加載歷史和能量累積對亞臨界裂縫起裂和擴展的影響,建立疲勞或長期加載條件下工程材料的相場疲勞?徐變模型,并應用于工程結構的全生命周期性能分析和壽命預測.

(4)考慮復雜多物理場耦合效應的相場分析方法.復雜服役環境下工程結構除了受到力學載荷作用外,往往還會遭受溫度場、化學場、電磁場等多物理場耦合作用,例如混凝土早齡期抗裂性能和高溫損傷破壞分析會涉及濕?化?熱?力多場耦合,頁巖氣、可燃冰、地熱等能源開采利用過程中也需要處理濕?熱?力多場耦合問題,電子元器件(如壓電陶瓷、鋰電池、半導體芯片等)研發和制造則涉及復雜的力?電耦合效應.作為一種裂縫相場?位移場耦合分析方法,統一相場理論在復雜多物理場耦合作用下的結構損傷破壞分析方面有著廣闊的前景.

(5)基體裂縫?界面裂縫耦合破壞的相場分析方法.工程結構往往并非由單一材料構成,而是由若干具有不同性能或功能的材料復合而成,典型實例包括鋼筋混凝土結構、新舊混凝土結構、FRP 加固混凝土結構、高鐵無砟軌道板結構、復合材料結構等.這些工程結構破壞時,不僅基體會出現裂縫起裂和擴展,不同材料層間還會出現界面裂縫,兩種破壞模式相互耦合、互相競爭,導致破壞模式具有多樣性.將相場正則化內聚裂縫模型與非連續方法聯合使用,結合兩種方法各自的優點,有望在復合材料結構損傷破壞分析方面發揮重要作用.

(6)細觀非均質材料的相場分析和材料?結構優化設計.在細觀尺度,工程材料的非均質性不僅會直接影響到裂縫擴展路徑和能量耗散,還會對工程結構的力學性能和耗能能力產生明顯影響.將統一相場理論與X 射線計算機斷層掃描(XCT)方法或隨機介質模型相結合,可以分析細觀非均質材料的損傷破壞行為[150].在此基礎上,結合拓撲優化設計方法,對工程材料的細觀結構進行優化設計和性能調控,有望獲取更高的材料強度和斷裂韌性.類似的思想也可以應用于工程結構,對可能出現裂縫的關鍵區域進行局部加強或對其它非重點區域進行局部削弱,誘導甚至主動控制裂縫路徑,實現工程結構的優化設計和性能調控.

(7)統一相場理論的高效數值實現方法.為了保證相場分析時裂縫面積和能量耗散的計算精度,有限元空間離散時破壞區域內的單元網格必須足夠細密,導致整個系統的自由度數量和求解規模均較大,即便是效率最高的整體BFGS 擬牛頓迭代算法也需要耗費較多的計算時間.因此,采用多尺度有限元、各向異性自適應網格劃分等方法降低計算量,進一步發展更為高效的數值求解算法,并采用GPU 并行計算等方法提高計算效率,也是將統一相場理論應用于大規模實際工程結構損傷破壞分析的重要方面.