固體結(jié)構(gòu)損傷破壞統(tǒng)一相場理論、算法和應(yīng)用1)

吳建營

(華南理工大學(xué)亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,廣州 510641)

引言

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

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

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

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

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

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

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

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

考慮到統(tǒng)一相場理論的相關(guān)研究工作分散發(fā)表于若干外文期刊,本工作將對該理論及其數(shù)值算法和應(yīng)用算例加以系統(tǒng)總結(jié)、整理和介紹,以期為工程結(jié)構(gòu)的損傷破壞分析提供一種可行的方法,并對該方向今后的研究工作進(jìn)行展望.

1 統(tǒng)一相場理論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(2)裂縫相場演化的能量準(zhǔn)則:裂縫相場的準(zhǔn)靜態(tài)演化滿足如下關(guān)系式

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

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

裂縫面積密度函數(shù)γ 的變分導(dǎo)數(shù)δdγ 表示為

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

其中,通量q與裂縫相場梯度?d之間滿足線性本構(gòu)關(guān)系,比例系數(shù)與裂縫尺度b>0 成正比.

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

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

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

上述裂縫相場演化的能量原理將斷裂力學(xué)和經(jīng)典損傷力學(xué)緊密聯(lián)系,構(gòu)建了兩種理論之間的橋梁.

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

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

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

同時,位移場u應(yīng)滿足經(jīng)典的力平衡方程

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

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

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

其中,系統(tǒng)總能量E 表示為

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

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

1.4 應(yīng)力?應(yīng)變本構(gòu)關(guān)系

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

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

有效應(yīng)力張量由如下線彈性關(guān)系給出

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

有效能量釋放率表示為

其中,E0和ν0分別為材料楊氏模量和泊松比.如圖2 所示,上述應(yīng)力?應(yīng)變本構(gòu)關(guān)系將給出完全對稱的受拉和受壓力學(xué)行為,這與脆性和準(zhǔn)脆性材料的實(shí)際力學(xué)特性不符.

圖2 不同有效能量釋放率定義給出的強(qiáng)度包絡(luò)線:平面應(yīng)力狀態(tài)Fig.2 Strength envelopes given by different effective energy release rates:Plane stress state

為考慮拉、壓應(yīng)力狀態(tài)下的不同力學(xué)行為,以單軸受拉強(qiáng)度為基準(zhǔn),將有效能量釋放率重新定義為

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

1.5 應(yīng)變局部化分析

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

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

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

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

式中引入了單調(diào)遞增的函數(shù)φ(d)

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

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

相應(yīng)的內(nèi)聚力確定為

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

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

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

2 脆性斷裂和準(zhǔn)脆性破壞相場模型

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

2.1 解析解:一般情況

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

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

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

2.2 解析解:特殊情況

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

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

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

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

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

峰值應(yīng)力對應(yīng)的裂縫相場dc與參數(shù)ξ ∈[0,2]有關(guān).

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

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

當(dāng)裂縫尺度b→0 時,峰值應(yīng)力σc無窮大,這一特性與經(jīng)典線彈性斷裂力學(xué)一致.相應(yīng)的,統(tǒng)一相場理論退化為一類脆性斷裂相場模型,見2.2.1 節(jié).

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

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

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

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

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

此種情況下,當(dāng)其他模型參數(shù)滿足一定條件時,統(tǒng)一相場理論可以退化為一類可以描述應(yīng)變軟化行為的準(zhǔn)脆性破壞相場模型,見2.2.2 節(jié)所述.

2.2.1 脆性斷裂

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

即[23,30]

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

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

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

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

常用模型有以下兩種

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

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

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

2.2.2 準(zhǔn)脆性破壞

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

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

為了保證峰值應(yīng)力σc為有限值(即材料破壞強(qiáng)度ft)且與裂縫尺度b無關(guān),根據(jù)式(43)有

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

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

圖7 準(zhǔn)脆性破壞相場模型給出的應(yīng)力?名義位移跳躍曲線[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+··· 為多項(xiàng)式P(d)在d=1時的取值.從式(37)和式(53)可以看出,裂縫帶寬度與裂縫尺度b成正比.

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

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

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

2.3 相場正則化內(nèi)聚裂縫模型

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值實(shí)現(xiàn)和求解算法

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

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

3.1 多場有限元數(shù)值實(shí)現(xiàn)

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

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

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

有限元空間離散后,位移場和應(yīng)變場、裂縫相場及其梯度分別表示為(在數(shù)值實(shí)現(xiàn)中,張量均采用Voigt 標(biāo)記)

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

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

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

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

3.2 裂縫相場演化不可逆性

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

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

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

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

上述兩種方法的共同點(diǎn)是將裂縫相場子問題不等式(62)2變?yōu)榈仁?即

3.3 控制方程的數(shù)值求解

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

3.3.1 整體牛頓迭代算法

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

子剛度矩陣分別表示為

由于所有被積分函數(shù)均連續(xù),可采用標(biāo)準(zhǔn)高斯數(shù)值積分方法進(jìn)行數(shù)值計算.

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

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

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

3.3.2 子問題交錯迭代算法

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

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

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

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

其線性化方程表示為

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

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

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

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

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

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

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

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

3.3.3 整體BFGS 擬牛頓迭代算法

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

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

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

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

3.4 數(shù)值實(shí)現(xiàn)平臺

從上述數(shù)值算法可以看出,統(tǒng)一相場理論非常便于通過有限元方法加以實(shí)現(xiàn),且其二維和三維數(shù)值實(shí)現(xiàn)幾乎沒有區(qū)別.比較常用的開源和商業(yè)有限元計算平臺有:

(1)Matlab:借助Matlab 軟件友好的開發(fā)環(huán)境和強(qiáng)大的圖形處理功能,文獻(xiàn)[45,119]等實(shí)現(xiàn)了脆性斷裂相場模型AT2,并應(yīng)用于固體的靜力和動態(tài)裂縫擴(kuò)展分析.

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

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

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

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

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

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

由于大多針對研究目的,上述數(shù)值實(shí)現(xiàn)平臺通常采用GMSH[127]、PARAVIEW[128]等開源軟件進(jìn)行前后處理,不便于工程應(yīng)用.

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

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

4 應(yīng)用算例

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

4.1 脆性破壞分析[58]

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

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

圖10 非對稱缺口開孔梁三點(diǎn)彎試驗(yàn):試件尺寸(長度單位: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]

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

接下來考慮其他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.

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

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

4.2 準(zhǔn)脆性破壞分析[97]

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

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

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

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

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

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

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

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

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

4.3 動力破壞分析[135]

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

圖17 厚壁圓柱動態(tài)破碎問題:幾何尺寸(左)和載荷時程(右).圓柱內(nèi)外半徑分別為rin=80 mm 和rout=150 mm,平面外厚度為70 mm.沖擊載荷的幅值為p0=400 MPa,參數(shù)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

根據(jù)文獻(xiàn)[137]中的討論,考慮材料斷裂能Gf服從如下Weibull 分布

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

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

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

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

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

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

為說明此點(diǎn),這里考慮如圖21 所示的金屬氫脆斷裂問題,即氫離子濃度聚集會降低金屬材料的破壞強(qiáng)度和斷裂能,加速結(jié)構(gòu)損傷破壞,這一病害在金屬管道、壓力容器等實(shí)際工程中經(jīng)常出現(xiàn).

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

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

圖20 厚壁圓柱動態(tài)破碎問題:裂縫尺度參數(shù)敏感性分析[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

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

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

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

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

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

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

圖23 金屬氫脆腐蝕問題:裂縫尺度b 對裂縫擴(kuò)展路徑的影響(氫離子濃度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 金屬氫脆腐蝕問題:不同裂縫尺度參數(shù)和有限元網(wǎng)格大小對載荷?位移曲線的影響[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]

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

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

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

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

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

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

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

圖28 進(jìn)一步對比了試驗(yàn)觀測和數(shù)值模擬給出的梁側(cè)面裂縫路徑.可以看出,數(shù)值模擬給出的裂縫路徑完全落在試驗(yàn)值范圍內(nèi)[145].

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

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

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

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

5 結(jié)論與展望

聚焦于固體結(jié)構(gòu)損傷破壞這一世紀(jì)難題,本文重點(diǎn)介紹了筆者近年來提出的統(tǒng)一相場理論、算法和若干應(yīng)用算例.

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

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

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

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

盡管統(tǒng)一相場理論獲得了較為成功的應(yīng)用,但在以下方面仍然存在亟需進(jìn)一步開展的研究工作.

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

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

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

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

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

(6)細(xì)觀非均質(zhì)材料的相場分析和材料?結(jié)構(gòu)優(yōu)化設(shè)計.在細(xì)觀尺度,工程材料的非均質(zhì)性不僅會直接影響到裂縫擴(kuò)展路徑和能量耗散,還會對工程結(jié)構(gòu)的力學(xué)性能和耗能能力產(chǎn)生明顯影響.將統(tǒng)一相場理論與X 射線計算機(jī)斷層掃描(XCT)方法或隨機(jī)介質(zhì)模型相結(jié)合,可以分析細(xì)觀非均質(zhì)材料的損傷破壞行為[150].在此基礎(chǔ)上,結(jié)合拓?fù)鋬?yōu)化設(shè)計方法,對工程材料的細(xì)觀結(jié)構(gòu)進(jìn)行優(yōu)化設(shè)計和性能調(diào)控,有望獲取更高的材料強(qiáng)度和斷裂韌性.類似的思想也可以應(yīng)用于工程結(jié)構(gòu),對可能出現(xiàn)裂縫的關(guān)鍵區(qū)域進(jìn)行局部加強(qiáng)或?qū)ζ渌侵攸c(diǎn)區(qū)域進(jìn)行局部削弱,誘導(dǎo)甚至主動控制裂縫路徑,實(shí)現(xiàn)工程結(jié)構(gòu)的優(yōu)化設(shè)計和性能調(diào)控.

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