基于相場(chǎng)模型的一維拉桿脆性斷裂分析

袁集海，陳昌萍?

（1.廈門(mén)大學(xué)建筑與土木工程學(xué)院，福建廈門(mén) 361005；2.廈門(mén)理工學(xué)院土木工程與建筑學(xué)院，福建廈門(mén) 361024）

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

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

1 相場(chǎng)損傷模型

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

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

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

且滿足下列邊界條件

對(duì)方程（3）運(yùn)用變分原理再在整個(gè)域內(nèi)積分可構(gòu)建相應(yīng)泛函

且在整個(gè)積分域內(nèi)滿足dV=Γdx，則

圖1 相場(chǎng)變量分布Fig.1 Distribution of phase field

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

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

對(duì)于多維斷裂問(wèn)題，在求解域Ω 內(nèi)相應(yīng)裂縫分布函數(shù)和裂縫分布密度函數(shù)分別為

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

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

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

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

采用線性單元對(duì)相場(chǎng)和位移場(chǎng)及其相應(yīng)增量進(jìn)行離散化

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

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

2 結(jié)果與討論

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

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

2.1 無(wú)缺陷受拉桿

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

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

2.2 含幾何缺陷受拉桿

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

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

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

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

圖6 相場(chǎng)、應(yīng)變和位移演化Fig.6 Evolution of phase field，strain and displacement

2.3 損傷受拉桿

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

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

圖8 給出了損傷受拉桿荷載位移曲線，由于桿在中點(diǎn)處發(fā)生了嚴(yán)重的損傷，損傷桿臨界位移荷載和最大內(nèi)力明顯遠(yuǎn)小于無(wú)損傷桿臨界位移荷載和最大內(nèi)力.圖9 給出了不同荷載步n 相場(chǎng)、應(yīng)變和位移分布，其對(duì)應(yīng)荷載位移點(diǎn)如圖8 所示.與含幾何缺陷桿不同的是，位移荷載較小時(shí)，預(yù)設(shè)損傷變量會(huì)導(dǎo)致桿結(jié)構(gòu)發(fā)生局部化現(xiàn)象，如圖9 所示.

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

圖9 相場(chǎng)、應(yīng)變和位移演化Fig.9 Evolution of phase field，strain and displacement

3 結(jié)論

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