基于一類非局部宏?微觀損傷模型的裂紋模擬1)

盧廣達 陳建兵

(同濟大學土木工程防災國家重點實驗室,土木工程學院,上海 200092)

引言

工程結構的損傷、破壞和解體往往是微裂紋形核、萌生、發展和傳播擴散的結果.因此,準確把握固體材料與結構中裂紋的產生和發展過程至關重要[1].自20 世紀20 年代初Griffith 的開創工作[2]以來,人們對此進行了大量的研究,并先后發展了線彈性斷裂力學[3-6]、彈塑性斷裂力學[7-11]和損傷力學[12-18].然而,固體中的裂紋萌生與傳播的模擬,迄今仍面臨一系列挑戰性問題.例如,經典斷裂力學難以解決裂紋形核以及擴展方向等問題[19].20 世紀末,通過引入能量釋放率最小化等準則和裂紋面密度連續化思想,基于變分原理發展了脆性材料的相場理論,為裂紋擴展模擬開辟了新的道路[20-21].最近,結合損傷力學和相場理論,Wu 等[22-25]提出了一類統一相場理論框架,對準脆性材料的靜力與動力裂紋擴展分析問題取得了重要突破.這一途徑,從計算角度來看,屬于彌撒型裂紋模型的范疇.在傳統上,人們認為斷裂力學與損傷力學是分別處理和分析裂紋出現之后與產生之前的兩種力學理論.統一相場模型的成功,從一個側面表明,這一認識存在局限性.事實上,二者之間是可以相容而相得益彰的.與此同時,人們意識到,僅僅在宏觀連續介質力學尺寸不足以把握材料損傷萌生與裂紋擴展的物理?力學本質,要從根本上解決上述問題,需要結合微細觀機理的分析,由此發展了宏微觀斷裂力學[1].另一方面,人們逐步認識到非局部化性質對克服經典局部損傷力學中網格敏感性問題的根本性意義[26].鑒于裂紋本質上是不連續的,20 世紀末以來人們提出了一系列新的途徑,包括內聚裂紋模型[27]、擴展有限元方法[28]和無網格方法[29].自2000 年以來,Silling 提出peridynamics(近場動力學[30]或毗域動力學[31])理論[32],將固體介質的相互作用采用積分和離散形式表達,給出了一類新的、具有非局部化性質的離散裂紋萌生和擴展模擬方法,并在一系列問題中取得了相當的成功[33-37].

雖然如此,統一相場理論需要額外求解相場演化方程,從而使得理論和求解的復雜程度較高.另一方面,近場動力學方法中本構模型的確定以及對準脆性裂紋的模擬尚存在困難.為此,結合統一相場理論和近場動力學的基本思想,針對準脆性材料的裂紋模擬,最近提出了一類新的非局部宏?微觀損傷模型[38].在此基礎上,本文進一步修正微觀損傷準則,并對具有強非線性回彈特性的典型裂紋擴展問題進行分析,驗證了所改進模型的有效性.

1 宏?微觀損傷模型的基本理論

從微細觀層次來看,裂紋的萌生和發展從幾何上表現為出現了不連續性,從宏觀上雖依然可認為是連續體,但與原初固體材料相比,出現了宏觀損傷.這一損傷將引起材料自由能的變化,使得本構關系中出現非線性,因而導致材料表現出非線性行為.本節將給出這一基本路徑的理論刻畫[38].

1.1 拓撲損傷

考慮連續固體B,其邊界為?B.該連續體中的物質點記為x=(x1,x2,x3)∈B.考察物質點對(x,x′),其中∈B,ξ=x′?x為兩點之間的空間差,則為相應的距離,這里為向量的長度.連續體發生變形后,物質點x和x′的位移分別記為u=u(x)和u′=u(x′).

根據變形幾何(如圖1 所示),可以定義物質點對(x,x′)之間的變形為

其中,ν=為物質點對的單位方向向量.由此,可定義物質點對的伸長量為

式中,H(·)為Heaviside 階躍函數,當λ0 時,H(λ)=0,否則H(λ)=1.

圖1 物質點對的變形Fig.1 Deformation of material point pair

可以合理地假設,當伸長量超過給定閾值ˉλ 時,物質點對之間的聯系(物質鍵) 開始弱化,直至完全分離,即物質鍵斷裂.這一過程顯然與單調非減的加載歷史參數有關.為此,可定義加載歷史參數為即物質鍵變形歷史的最大超越伸長量,其中t為時間變量(對靜力問題應理解為虛擬時間).

注記1物質鍵[38]的觀念受啟發于近場動力學方法[32-34],類似的概念也見于分子動力學[39]和虛內鍵模型[40]中.需要指出的是,在本文中是臨界伸長量,而不是文獻[38]中的臨界伸長率.

不難驗證,歷史最大超越伸長量有dtκ0,這里,后續符號類似.因此,對于任意物質鍵(x,x′) 可以定義一個刻畫其聯系強弱程度的量ω(x,x′,t).不失一般性,可以對其進行歸一化,即ω ∈[0,1].它是歷史最大超越伸長量κ 的單調非減函數,即

注記2從物理意義上看,ω 是微細觀損傷的度量,其具體形式可由分子動力學的相互作用勢函數通過計算加以確定,或者根據理性假設給出.在本文中,假設為如下形式

其中γ>0 為模型參數,它表征微細觀損傷的發展速度.不難注意到,是以γ 為參數的、Heaviside 階躍函數H(κ)的極限函數列,而Heaviside 階躍函數表征了不連續性.因此,從這個意義上說,微細觀損傷刻畫了固體在物質鍵層次的不連續程度,這一點與李杰等[18]認為損傷變量應打破物質空間連續性的觀點是一致的.

事實上,連續體中某物質點處的損傷,往往不僅取決于該點本身,還必然與其周邊一定鄰域范圍內的物質點相互作用有關.換言之,其損傷的變化,本質上是非局部化的.記存在影響的鄰域(影響域)半徑為?,即?=diamD?(x)/2,這里D?(x)表示物質點x的影響域.根據上述表述,可以進一步提煉出損傷的數學語言:

定義1設連續體B?Rd(d=1,2,3)及其位移場u?Rd,若給定影響域半徑? > 0,存在閾值ˉλ > 0,使得物質點x∈B有

則稱連續體在x處發生損傷,此時微細觀損傷ω>0.它從數學上刻畫了位移場u在空間上關于方向ν的不連續程度,因此不妨稱這一定義為ˉλ ?? 語言.

作為對照,這里給出位移場u作為d元d維向量值函數其連續性的? ?δ 語言[41].

定義2設u:B?Rd→Rd,若對任意? > 0,存在δ>0,使得x∈B有

其中Dδ(x)=,則稱向量值函數u在x處連續.

注記3定義1 和定義2 恰好構成關于函數不連續性和連續性的一組對偶表述,而這要得益于本文采用物質鍵伸長量表征微細觀損傷的改進.需要說明的是,這一損傷的?? 語言是文獻[38]中尚未認識到的,并區別于經典連續損傷力學理論[13,17]中損傷變量只表征材料剛度退化的觀念.

對?x∈B物質點給定的影響域D?(x),設任意兩物質點之間的影響函數為?(x,x′),并滿足如下條件

從宏觀上看,對于固體B內任意物質點x,其損傷應是該物質點所連接的物質鍵之微細觀損傷的加權平均.因此,對?x∈B定義宏觀層次損傷為

即根據微細觀層次的物質鍵(x,x′) 的斷裂程度在局部空間上進行幾何上的加權平均,它連續地刻畫了位移場在一點附近的整體不連續性.不難注意到,損傷區域J?={x∈B:?(x)>0}構成了裂紋拓撲J(其中dimJ=d?1)的非局部化鄰域表達.因此,從這個意義上,不妨稱之為拓撲損傷.

注記4拓撲損傷? 雖然在形式上與光滑處理的非局部損傷[42]相同,但二者存在著本質的區別:前者積分項內的損傷變量是建立在點對(兩點)意義上的,而后者則是基于一點的.需要指出的是,基于一點的非局部損傷已經被證明存在應力鎖死的問題[43],而拓撲損傷則不存在此問題[38].特別地,當微細觀損傷ω 取為Heaviside 階躍函數時,則式(10)退化為近場動力學的損傷表示[33].

注記5事實上,式(10) 以積分形式定義的拓撲損傷與統一相場理論[22]中含梯度算子的橢圓泛函[44]

的作用相同,即采用非局部化思路將裂紋拓撲表示為含一定帶寬的彌散型界面,從而使得分析過程無需人為預設裂紋擴展路徑.其中,式(11)的參數φ 和α 分別為相場函數和裂紋幾何函數.

進一步地,可以證明以下命題:

命題1設裂紋拓撲J?Rd,以及? 和Θ 分別如式(10)和式(11)所定義,則

式中,?(d)=Υ(d)/2,其中Υ(d) 表示d維單位球域D?=1(0)的表面積測度.

需要說明的是,上式第一個等式(即橢圓泛函逼近)已經由Braides[44]給出了證明,而第二個等式(拓撲損傷逼近) 的具體證明因限于本文的篇幅而考慮另文給出.

1.2 能量退化因子

注意到,連續體的自由能密度必然隨著拓撲損傷的發展(即裂紋萌生和擴展)而退化.設能量退化因子為g,則有損連續體的自由能為

式中小應變張量ε=?Su,其中?S(·) 表示對稱梯度算子,ψ0為無損材料的自由能密度函數,即

其中E為四階彈性張量.

為簡化計,在本文中暫不考慮塑性的影響.因此,這里彈性應變和總應變相同,后文對二者將不再加以區分.

顯然,能量退化因子與拓撲損傷有關,因而它是拓撲損傷的函數g=g(?),表征了物質點儲能能力的退化,稱為能量退化因子[45]且滿足如下要求[22]

其中,g(0)=1 表示物質點處于線彈性無損狀態;g(1)=0 則是完全損傷的刻畫,即喪失了全部的儲能能力; 此外,d?g< 0 表示損傷總是耗能的; 而d?g(1)=0 則保證當?=1 時,損傷驅動力??ψ=d?g(?)ψ0趨于零,即當一點的損傷發展完全后,其驅動力應隨之消失.

若從能量意義上說,不妨稱Φ=1 ?g(?) 為基于能量的損傷.事實上,基于拓撲損傷的能量退化因子還可以進一步從物理上論證,g是[0,1] 上的凸函數[38],其定性的物理解釋如下.

由于拓撲損傷? 考慮的是影響域內所有方向上物質鍵損傷的幾何“平均”,而各個方向上的物質鍵并非同時地發生損傷,因為物質鍵的伸長量在各個方向上的分布一般是不均勻的.在初始階段,只有位移場在空間變化劇烈方向上那些少數(相對于所有方向上的物質鍵總數而言) 超過閾值的物質鍵率先進入損傷狀態,此時拓撲損傷? 在數值上非常小,但這些率先發生損傷的物質鍵因變形劇烈而儲存了較大比重的變形能.因此,在?=0 附近處基于能量的損傷Φ 要大于拓撲損傷?,即隨著其他方向物質鍵損傷的發展及其能量的釋放,逐漸地一點鄰域內所能釋放的變形能比重越來越小,即能量退化因子g隨著拓撲損傷? 的增大而趨于平緩.

根據上述解釋,可以清楚地看到:在本模型中能量退化因子是拓撲損傷的非線性函數,不同于經典連續損傷力學[13,17]1 ?? 的線性形式.借鑒于統一相場模型[22]的結果并加以適當修正,可取能量退化因子為如下凸函數[38]

其中,p和q均為退化函數.一些代表性取值的函數曲線如圖2 所示.

圖2 能量退化因子Fig.2 The energetic degradation factor

2 控制方程

2.1 本構方程

根據Clausius-Duhem 不等式[17],有能量耗散率dtΞ 不等式

其中,σ為Cauchy 應力張量.將式(13) 代入上述不等式并稍作整理,得

由前述拓撲損傷即能量退化因子的定義和性質可知,上述積分項內第二項必然不大于零.因此,該不等式成立的充要條件是

式中,σ0為有效應力張量.

注記6由此可見,能量退化因子表征了材料的物性關系,式(20)與傳統損傷力學理論[18]的本構關系在形式上完全相同.其中,本文式(17)退化因子給出的是準脆性材料本構,p和q的取值刻畫了材料的宏觀脆性程度,可認為是物性參數,應根據標準試件進行標定.

2.2 平衡方程與邊界條件

平衡方程和邊界條件與常規連續介質力學的基本表達相同.考慮圖3 所示固體B的參考構形,容易列出平衡方程為

其中,t為內力向量,為體力密度向量(當不計體力時,=0).

圖3 固體的構形及其邊界Fig.3 The solid configuration and its boundary

根據Cauchy 假定,即t=σ·nP,其中nP為子集P的邊界外法向向量.由式(20)所表述的本構關系,利用散度定理并注意到P的任意性,可得

而在邊界上則分別有

其中,nB為固體的邊界外法向向量,和分別為給定的邊界面力和邊界位移.

注記7得益于式(10)所定義的拓撲損傷關于位移場是Lipschitz 連續的(具體證明見文獻[38]),本文模型仍然可以采用連續介質力學的Cauchy 假定建立平衡方程.這一點區別于近場動力學積分形式的平衡方程中關于應力散度項的低階逼近[46]

其中,f和T分別為近場動力學內力密度和內力狀態向量.盡管通過引入高階微分算子[47]可以一定程度提高其逼近精度以及規避因內力作用方式重構而帶來的零能模式問題[37],但也增加了問題求解的復雜性和計算量.此外,近場動力學方法本來的計算效率就低于有限元方法[48].

3 數值求解方法

3.1 有限元離散方程

上述問題,可以直接采用常規有限單元法進行離散和求解.式(22)和式(23)這一邊值問題經過弱形式逼近后,可得單元的平衡方程為

式中,Ue為單元結點位移向量,Ke為單元剛度矩陣,Re為單元結點力向量,其中

這里Ae為單元區域,?tAe為邊界單元的力邊界,Be和Ψe分別為幾何矩陣和形函數矩陣,其表達式可在一般有限元專著中找到,茲不贅述.此外,De為彈性矩陣,對于平面應力問題有

其中,E為彈性模量,μ為泊松比.

根據單元組集規則,可得結構整體平衡方程為

其中,整體剛度矩陣K、整體等效載荷向量R和整體位移向量U分別為

式中,Λe為單元轉換矩陣,Ne為單元數目.

3.2 單元剛度矩陣的計算

在式(26)中,為了計算單元剛度矩陣,需要計算給定積分點的能量退化因子g.這時,應通過離散化方式計算式(10) 中的拓撲損傷?,然后將其代入能量退化因子式(17)中.此時,拓撲損傷可根據離散化后的各單元變形計算.對于常應變單元,有

其中,xe和xe′分別為單元Ae和Ae′的形心坐標,而=vol[D?(xe)∩Ae′].若考慮采用等參單元進行離散,則二者為相應的高斯點及其權重體積.

注意到,能量退化因子是位移的函數.由式(26)可知,單元剛度矩陣也是位移的函數,因而由式(30)給出的整體剛度矩陣亦為位移的函數.因此,平衡方程(29)是一個非線性方程.在本文中,采用改進局部弧長法進行求解.限于篇幅,茲不贅述.

注記8在求解中,拓撲損傷? 是整體節點位移向量U的“顯式”函數,即一旦求得整體節點位移向量U已知,直接將其代入式(31)就可以得到拓撲損傷?,整個過程自始至終只需要求解(29) 的平衡方程,無需額外的待求解方程.這一點比相場模型簡單而高效,因為相場模型中相場函數和位移場均是待求未知量,需要采用聯立或交替的策略求解平衡方程和額外的相場演化方程[49].

4 典型實例分析

以Ingraffea 和Grigoriu 所做的有機玻璃剪切梁試件[50]為例,來驗證和分析本模型的有效性.注意到,該試件以其強非線性回彈特性而著稱,是檢驗數值模型優良性的經典算例.近些年來,經過該算例檢驗的數值模型有基于混合(應變/位移)元的各向同性損傷模型[51]和統一相場模型[23]等.分析中試件按平面應力問題處理,其材料參數按文獻[51]取:彈性模量E=3.102 GPa,泊松比μ=0.35.宏?微觀損傷模型中的參數取:影響域半徑?=4 mm,臨界伸長量=1.38×10?2mm,γ=2.0×105,以及p=4 和q=4.

4.1 不含圓孔情況

注意到,經典的局部損傷力學存在網格敏感性[18]:一旦所分析的固體存在或出現初始裂紋時,當計算點的間距加密時,基于經典局部損傷力學的分析結果,將會出現即使是很微小的外力作用就可以引起嚴重的損傷(當網格非常密集時,更甚者會出現不需要任何外力就可以使固體發生破壞的問題) 這一違背物理真實的現象.

為了考察本文損傷模型的網格敏感性問題,在有限元離散時采用3 種不同疏密程度的三角形網格進行劃分.為此,首先考慮不含圓孔的情況,試件形狀和尺寸如圖4 所示.

圖4 剪切梁基本尺寸(切口寬度為2 mm)Fig.4 Size of the shear beam(The notch width=2 mm)

具體地,考慮在切口一定范圍的區域進行單元加密處理,3 種不同疏密程度的單元總數目分別為17 308,23 564 和55 262,單元分布情況見圖5.

圖5 有限元網格劃分Fig.5 Finite element meshes

圖6 裂紋發展模式Fig.6 The crack development modes

圖6 裂紋發展模式(續)Fig.6 The crack development modes(continued)

圖7 不同網格的載荷?位移曲線Fig.7 Load-deformation curves for different meshes

采用本文方法進行分析和求解.3 種不同網格劃分情況下的裂紋擴展模式如圖6 所示.與圖6(d)的試驗結果對比可見,3 種不同網格的分析結果均與試驗結果相符,且隨著網格加密裂紋擴展路徑趨于試驗結果.圖7 進一步給出了外加載荷?切口張開位移和外加載荷?加載點位移曲線.特別值得指出的是,圖7(b)表現出一般數值方法很難處理的強非線性回彈性質.從圖中可見,3 種不同網格的數值分析結果具有良好的一致性.可見采用本文建議方法分析的結果克服了經典局部損傷力學所面臨的網格敏感性問題.

圖8 是載荷?切口張開位移和載荷?加載點位移曲線與試驗結果的對比.從圖8(a)可見,數值模擬的載荷?切口張開位移曲線與試驗曲線的吻合良好.在圖8(b)中,雖然未能完美吻合,但載荷?加載點位移曲線的最大載荷與最終加載點位移與實驗值的接近,且總體的定性特征一致.

圖8 數值與試驗結果比較Fig.8 Comparisons between the numerical and experimental results

圖9 和圖10 分別給出了本文建議方法與基于混合元的各向同性損傷模型[51]和統一相場模型[23]分析所獲得的載荷?變形曲線與裂紋模式的對比,可見三者較為一致.

圖9 與其他數值模型的載荷?位移曲線比較Fig.9 The load-deformation curves compared with other numerical models

圖10 與其他數值模型破壞模式的比較Fig.10 The comparisons of failure modes with other numerical models

圖11 典型載荷步損傷與位移云圖Fig.11 Damage and displacement field at typical loading steps

圖11 典型載荷步損傷與位移云圖(續)Fig.11 Damage and displacement field at typical loading steps(continued)

圖11 是圖8 的載荷?變形曲線中幾個典型階段的損傷云圖和位移云圖演化情況.從圖中可見,裂紋隨著加載過程而不斷發展,同時位移場在裂紋兩側出現顯著的不連續性.特別值得注意的是,當外加載荷已經達到峰值的時候,裂紋才剛剛具有肉眼可見的發展(圖11(a)),從圖11(b)亦可見,此時位移場尚未出現顯著的不連續性.當裂紋明顯萌生和發展時(圖11(c),圖11(e),圖11(g)),試件的載荷?變形曲線已經進入顯著的下降段(圖8(a))或回彈段(圖8(b)).換言之,對于此類問題,當出現肉眼可見的裂紋時,裂紋已經進入不穩定發展階段.這與人們從實際工程的損傷和斷裂分析中獲得的經驗完全相符.

4.2 含圓孔情況

考慮含圓孔情況,試件形狀與尺寸如圖12 所示.類似地采用3 種不同疏密程度的三角形網格進行劃分,單元數目和網格劃分具體情況見圖13.

圖12 剪切梁基本尺寸(切口寬度為2 mm)Fig.12 Size of the shear beam(the notch width=2 mm)

圖13 有限元網格劃分Fig.13 Finite element meshes

圖14 裂紋發展模式Fig.14 The crack development modes

采用本文建議的方法,得到3 種不同的網格劃分獲得的裂紋最終發展模式如圖14 所示.與圖14(d)的試驗結果對比表明,3 種網格均能良好地捕捉裂紋發展過程,并給出了試驗未展示的裂紋后續擴展部分.圖15 給出了載荷?變形曲線的分析結果.同樣地,在圖15(b) 中也表現出強非線性回彈特性.由此可見,采用不同的網格時,無論是載荷?裂紋張開位移曲線、還是載荷?加載點位移曲線都具有很好的一致性,不存在經典局部損傷力學中的網格敏感性問題.

圖15 不同網格的載荷?位移曲線Fig.15 Load-deformation curves for different meshes

圖16 是采用網格C 時的載荷?變形曲線.與圖8 類似,此時載荷?加載點位移曲線出現了明顯的回彈現象.

圖17 是與圖16 中標注的典型階段相應的損傷云圖與位移云圖.與前述情況類似地可見,當外載荷達到峰值時,試件中才剛剛出現可見裂紋的擴展.而且,裂紋兩側位移場出現顯著的不連續性.有意思的是,當存在圓孔時,裂紋的發展會被圓孔所誘導(圖17(c)),并最終穿透圓孔繼續發展(圖17(e)).值得注意的是,在被圓孔俘虜后,當裂紋進一步穿透圓孔時,出現了載荷的小幅上升(見圖16),然后繼續下降,直至完全破壞并基本喪失承載力.

圖16 載荷?位移曲線(網格C)Fig.16 The load-deformation curve(mesh C)

圖18 是本文模型與統一相場模型分析結果[23]的對比.可見二者的裂紋擴展模式是基本一致的.

5 結論

在連續損傷力學的框架下,結合相場理論與近場動力學的基本思想,基于一類新的非局部宏?微觀損傷模型,實現了典型試件的裂紋發展過程模擬和分析.主要結論如下.

(1)本文對非局部宏?微觀損傷模型進行了改進.該模型可以方便地在連續介質力學與有限元離散的框架下求解.與統一相場模型相比,不需要求解額外的相場演化方程.與近場動力學相比,不需要根據連續介質力學重新確定本構關系.

(2)應用該損傷模型分析固體破壞問題,不需要人為預設裂紋擴展路徑,可自動描述和跟蹤裂紋擴展.不僅可以獲得正確的裂紋模式,而且可以定量地獲得載荷?變形曲線.

圖17 典型載荷步損傷與位移云圖Fig.17 Damage and displacement field at typical loading steps

圖18 與統一相場模型破壞模式的比較Fig.18 The failure mode compared to that obtained by the unified phase field model

(3)該損傷模型是一類非局部化模型,不存在經典局部損傷力學所面臨的網格敏感性問題.

本文的研究工作初步驗證了非局部宏?微觀損傷模型的有效性.在未來工作中,尚有大量研究工作需要深入推進,例如基本參數與斷裂能之間的內在關系、塑性變形的考慮、各向異性問題、動態裂紋擴展與分叉等挑戰性問題.