袁集海，陳昌萍?

（1.廈門大學建筑與土木工程學院，福建廈門 361005；2.廈門理工學院土木工程與建筑學院，福建廈門 361024）

上世紀末，F(xiàn)rancfort 和Marigo[1]基于能量最小化原理，提出了解決固體結構斷裂損傷的變分方法.該方法認為裂紋的出現(xiàn)和演化使總勢能最小，其實質是從能量變分的角度對Griffith[2]脆性斷裂理論的發(fā)展和推廣，然而在連續(xù)介質力學框架下求解過程中仍然存在位移不連續(xù)帶來的求解困難，該工作為二階相場模型的建立奠定了理論基礎.Bourdin[3]等在斷裂變分方法的基礎上，通過引入在0 到1 之間連續(xù)變化的相場變量d-來表征裂縫的起裂和延伸，同時引入裂縫寬度參數(shù)l 來控制裂縫的彌散程度.隨后，Miehe[4]等建立了與熱力學自洽的二階相場損傷模型.二階相場模型的完善迅速吸引了一大批國內外學者的關注，為后續(xù)相場損傷模型的發(fā)展做出了重要的貢獻.Borden[5]等提出了四階相場損傷模型并采用等幾何法對耦合方程進行求解，另外，還研究了動態(tài)脆性斷裂[6]和塑性材料[7]斷裂問題.鄧俊俊[8]采用無網(wǎng)格相場模型研究了復雜邊坡破壞問題.張飛[9]基于自適應移動網(wǎng)格方法優(yōu)化了相場模型并模擬了水力裂縫延伸過程.付禹銘[10]針對各向異性損傷相場模型提出了新的投影算子算法.仇杰峰[11]和莊洛嘉[12]采用相場損傷模型研究了混凝土破壞問題.

相場損傷模型采用連續(xù)方法來描述不連續(xù)問題，與傳統(tǒng)求解斷裂問題方法相比具有一定的優(yōu)勢[13-14].本文以受拉結構為研究對象，采用相場損傷模型數(shù)值模擬了無缺陷桿、含幾何缺陷桿以及含損傷桿應變局部化行為和脆性斷裂行為.

1 相場損傷模型

考慮一維桿斷裂問題，其長度和截面分別為L和A.假設桿在x∈Γd處斷裂，其他地方完好無損，并引入相場變量d（x）來描述桿的破壞狀態(tài)，則該變量分布（如圖1（a）所示）為

為了便于數(shù)值計算，引入如下形式的指數(shù)函數(shù)來逼近上述不連續(xù)函數(shù)分布，即

式中：l 為正則化參數(shù)，該參數(shù)可以控制相場變量的分布區(qū)域，如圖1（b）所示，正則化參數(shù)越小，相場變量分布區(qū)域也越小，當l 趨近于0 時，該連續(xù)函數(shù)分布趨近于上述不連續(xù)函數(shù)分布.顯然，式（2）是如下微分方程的解[4].

且滿足下列邊界條件

對方程（3）運用變分原理再在整個域內積分可構建相應泛函

且在整個積分域內滿足dV=Γdx，則

圖1 相場變量分布Fig.1 Distribution of phase field

因此，斷裂面Γ 可用裂縫分布函數(shù)近似表示為

該函數(shù)描述了斷裂面處裂縫的分布情況.引入裂縫分布密度函數(shù)

對于多維斷裂問題，在求解域Ω 內相應裂縫分布函數(shù)和裂縫分布密度函數(shù)分別為

基于格里菲斯斷裂理論，彈性體在單位面積上耗散的能量為Gc，則在斷裂面處Γl（d）結構斷裂耗散能量[8].

在不考慮動能的情況下，根據(jù)能量守恒定律，彈性體總能量變化率為0，即

將式（13）～（15）代入能量守恒方程即可得到相場損傷模型控制方程

采用有限元方法對方程（16）進行數(shù)值求解，平衡方程伽遼金弱形式

采用線性單元對相場和位移場及其相應增量進行離散化

式中：N、B 分別為形函數(shù)和梯度矩陣；a、d 為節(jié)點位移和相場.將式（24）代入弱形式和增量方程，并將所得結果代入方程（23）可得增量迭代方程組.本文采用交錯最小算法[15]，故忽略耦合剛度后的迭代方程組為

式中：剛度矩陣和節(jié)點力向量分別為

2 結果與討論

考慮如圖2 所示一維受拉桿，其長度L=100 mm，截面面積A=4 mm×4 mm，左端固定，右端施加位移荷載u（L），不計體力.若無特別說明，彈性模量E=210 GPa；泊松比υ=0.3；正則化參數(shù)l=0.2 mm；斷裂能釋放率Gc=2.7 × 10-3kN/mm.本節(jié)通過三個數(shù)值算例分析一維受拉桿結構脆性破壞.

圖2 一維受拉桿Fig.2 One dimensional bar under tension

2.1 無缺陷受拉桿

對于如圖2 所示無缺陷受拉桿，由于應變和相場均勻分布，故可忽略控制方程中相場變量梯度項，即可得到解析解.圖3 給出了荷載位移曲線有限元解和解析解，其中，虛線為解析解，實線為有限元解.對于有限元解，當位移荷載超過臨界荷載后，曲線垂直下降到0，這是因為桿發(fā)生了脆性斷裂，開裂后，內力立即減小為0.而對于解析解，曲線上升階段與有限元解保持一致.另外，值得指出的是解析解曲線只有上升階段才有意義，脆性斷裂無軟化階段.

圖3 有限元解與解析解對比Fig.3 Comparision beween numerical solution and analytical solution

2.2 含幾何缺陷受拉桿

本節(jié)第二個數(shù)值算例為含缺陷受拉桿，其長度為L=10 mm，截面面積為A=1 mm×1 mm；中間部分為弱化區(qū)，長度為l0=1 mm，截面積為A′=0.9A，如圖4 所示.由解析解可以得到臨界位移荷載uc=0.032 7 mm；加載過程荷載步長Δu=uc/2 000.

圖4 含幾何缺陷受拉桿Fig.4 One dimensional bar with geometric imperfection under tension

圖5 給出了含缺陷桿（實線）和無缺陷桿（虛線）荷載位移曲線，從圖中可以看出曲線上升階段二者保持一致，下降階段無缺陷桿最大內力和臨界位移荷載均大于含缺陷桿最大內力和臨界位移荷載.圖6給出了不同荷載步n 相場、應變和位移分布，其對應荷載位移點如圖5 所示.當位移荷載較小（n=500）時，相場和應變均勻分布，位移線性分布；當位移荷載接近臨界位移荷載時，相場和應變開始局部化，其分布集中在弱化區(qū)，且在該區(qū)位移不再線性分布.當位移荷載超過臨界荷載（n=1 889）時，在桿中點處發(fā)生斷裂（相場在該點達到1，應變趨于無窮大，位移均勻分布：左側位移為0，右側位移為所施加位移荷載）.

圖5 荷載位移曲線Fig.5 Load deflection curve

圖6 相場、應變和位移演化Fig.6 Evolution of phase field，strain and displacement

2.3 損傷受拉桿

類似損傷力學中損傷變量的概念，可以將相場變量視作量化結構損傷指標.本節(jié)最后一個算例在桿中點預設相場變量，如圖7 所示.其他材料參數(shù)和幾何參數(shù)同算例2，荷載步Δu=uc/4 000.

圖7 損傷受拉桿Fig.7 One dimensional bar with damage under tension

圖8 給出了損傷受拉桿荷載位移曲線，由于桿在中點處發(fā)生了嚴重的損傷，損傷桿臨界位移荷載和最大內力明顯遠小于無損傷桿臨界位移荷載和最大內力.圖9 給出了不同荷載步n 相場、應變和位移分布，其對應荷載位移點如圖8 所示.與含幾何缺陷桿不同的是，位移荷載較小時，預設損傷變量會導致桿結構發(fā)生局部化現(xiàn)象，如圖9 所示.

圖8 荷載位移曲線Fig.8 Load deflection curve

圖9 相場、應變和位移演化Fig.9 Evolution of phase field，strain and displacement

3 結論

基于相場模型，本文采用有限元數(shù)值方法通過三個數(shù)值算例分析了一維受拉桿在位移荷載下脆性斷裂破壞.算例一模擬了無缺陷桿受拉斷裂破壞，并將數(shù)值解與解析解對比，論證了數(shù)值解的正確性.算例二和算例三分別分析了含幾何缺陷桿和損傷桿受拉破壞.計算結果表明有限元相場分析方法能夠較準確地模擬破壞區(qū)應變局部化行為和破壞演化過程.