基于相場模型的一維拉桿脆性斷裂分析

袁集海 陳昌萍

摘? ?要：基于相場模型分析了受拉結構脆性斷裂破壞，介紹了裂縫分布函數和能量耗散函數的構造，根據能量守恒定律推導了相場損傷模型控制方程，并采用交錯數值計算方法對方程進行求解.通過對比一維無缺陷受拉桿數值解和解析解論證了數值解的正確性.隨后模擬了含幾何缺陷受拉結構和含損傷受拉結構斷裂破壞，數值結果表明有限元數值模擬相場分析方法能夠較準確地模擬破壞區應變局部化行為和脆性斷裂.

關鍵詞：一維拉桿;相場模型;脆性斷裂;數值模擬

Abstract：Based on phase-field model，brittle fracture of the structures under tension is studied numerically. An introduction for construction of crack distribution function and energy dissipation function is presented. The governing equations for the phase field model are derived according to the conservation of energy and further solved by alternating minimization （staggered） algorithm. Good agreements between analytical solutions and numerical solutions of one-dimensional bar under tension validate the accuracy of analytical solutions. Numerical simulations of structures with geometric imperfection and predefined damage under tension are also presented. The numerical results show that deformation localization and brittle fracture in failure areas can be modeled accurately.

Key words：one-dimensional pull rod;phase-field model;brittle fracture;numerical simulation

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

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

1? ?相場損傷模型

2? ?結果與討論

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

2.1? ?無缺陷受拉桿

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

2.2? ?含幾何缺陷受拉桿

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

圖5給出了含缺陷桿（實線）和無缺陷桿（虛線）荷載位移曲線，從圖中可以看出曲線上升階段二者保持一致，下降階段無缺陷桿最大內力和臨界位移荷載均大于含缺陷桿最大內力和臨界位移荷載. 圖6給出了不同荷載步n相場、應變和位移分布，其對應荷載位移點如圖5所示. 當位移荷載較小（n = 500）時，相場和應變均勻分布，位移線性分布;當位移荷載接近臨界位移荷載時，相場和應變開始局部化，

其分布集中在弱化區，且在該區位移不再線性分布.當位移荷載超過臨界荷載（n = 1 889）時，在桿中點處發生斷裂（相場在該點達到1，應變趨于無窮大，位移均勻分布：左側位移為0，右側位移為所施加位移荷載）.

2.3? ?損傷受拉桿

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

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

3? ?結? ?論

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

參考文獻

[1]? ? FRANCFORT G A，MARIGO J J. Revisiting brittle fracture as an energy minimization problem[J]. Journal of the Mechanics and Physics of Solids，1998，46（8）：1319—1342.

[2]? ? GRIFFITH A A. The phenomena of rupture and flow in solids[J]. Philosophical Transactions of the Royal Society of London，1921，221：163—198.

[3]? ? BOURDIN B，FRANCFORT G A，MARIGO J J. Numerical experiments in revisited brittle fracture[J]. Journal of the Mechanics and Physics of Solids，2000，48（4）：797—826.

[4]? ? MIEHE C，WELSCHINGER F，HOFACKER M. Thermodynamically consistent phase-field models of fracture：Variational principles and multi-field FE implementations[J]. International Journal for Numerical Methods in Engineering，2010，83（10）：1273—1311.

[5]? ? BORDEN M J，HUGHES T J R，LANDIS C M，et al. A higher-order phase-field model for brittle fracture：Formulation and analysis within the isogeometric analysis framework[J].Computer Methods in Applied Mechanics and Engineering，2014，273：100—118.

[6]? ? BORDEN M J，VERHOOSEL C V，SCOTT M A，et al. A phase-field description of dynamic brittle fracture[J].Computer Methods in Applied Mechanics and Engineering，2012，217/218/219/220：77—95.

[7]? ? BORDEN M J，HUGHES T J R，LANDIS C M，et al. A phase-field formulation for fracture in ductile materials：Finite deformation balance law derivation，plastic degradation，and stress triaxiality effects[J].Computer Methods in Applied Mechanics and Engineering，2016，312：130—166.

[8]? ? 鄧俊俊. 基于相場模型的無網格損傷模擬[D]. 廈門：廈門大學，2017.

DENG J J. Meshfree failure simulations with phase field models[D].Xiamen：Xiamen University，2017. （In Chinese）

[9]? ? 張飛.基于自適應移動網格及相場逼近的水力裂縫延伸模擬[D]. 北京：中國石油大學（北京），2018.

ZHANG F. A study on adaptive finite element solution of phase-field models for hydraulic fracturing[D]. Beijing：China University of Petroleum （Beijing），2018. （In Chinese）

[10]? 付禹銘.基于相場法的脆性斷裂數值模擬[D].成都：西南交通大學，2019.

FU Y M.Numerical simulation of brittlefracture based on phase field method[D].Chengdu：Southwest Jiaotong University，2019.（In Chinese）

[11]? 仇杰峰.混凝土相場損傷模型與擴展有限元方法對比研究[D]. 廣州：華南理工大學，2019：1—87.

QIU J F. Comparative study of phase-field regularized cohesive zone model and extended finite element method in concrete [D]. Guangzhou：South China University of Technology，2019：1—87. （In Chinese）

[12]? 莊洛嘉.基于統一相場損傷理論的混凝土破壞全過程分析應用研究[D]. 廣州：華南理工大學，2019：1—93.

ZHUANG L J Application study of unified phase field theory for modeling localized failure in concrete [D]. Guangzhou：South China University of Technology，2019：1—93. （In Chinese）

[13]? WU J Y，NGUYEN V P，NGUYEN C T，et al. Phase-field modeling of fracture[J]. Advances in Applied Mechanics. https：//www.researchgate.net/publication/326258921.

[14]? WU J Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure[J]. Journal of the Mechanics and Physics of Solids，2017，103：72—99.

[15]? MIEHE C，HOFACKER M，WELSCHINGER F. A phase field model for rate-independent crack propagation：Robust algorithmic implementation based on operator splits[J]. Computer Methods in Applied Mechanics and Engineering，2010，199（45/46/47/48）：2765—2778.