廖 鑫 鄧 強(qiáng) 張 超

（成都飛機(jī)工業(yè)（集團(tuán)）有限責(zé)任公司，四川 成都 610092）

脆性材料（例如陶瓷、巖石等）在日常生活中隨處可見，在溫度載荷作用下，產(chǎn)生的內(nèi)應(yīng)力會(huì)導(dǎo)致斷裂破壞，如何準(zhǔn)確模擬脆性材料的復(fù)雜裂紋萌生和擴(kuò)展一直是國(guó)內(nèi)外研究的難題。

PD 理論假設(shè)連續(xù)體中粒子之間存在相互作用，不需要指定額外的失效準(zhǔn)則，材料的斷裂損傷自然發(fā)生，常用于研究裂紋擴(kuò)展問(wèn)題。Silling S A 等[1]根據(jù)PD 理論成功地捕捉了帶有初始缺陷厚板的裂紋萌生和擴(kuò)展現(xiàn)象，證明了近場(chǎng)動(dòng)力學(xué)可以用來(lái)模擬裂紋擴(kuò)展問(wèn)題。秦洪遠(yuǎn)等[2]構(gòu)建了可以消除泊松比限制的雙參數(shù)微彈脆性本構(gòu)模型，模擬了含雙裂紋巴西圓盤的劈裂破壞過(guò)程。孫杰等[3]同樣使用雙參數(shù)微彈脆性本構(gòu)模型研究了初始裂紋角度對(duì)多裂紋擴(kuò)展路徑和臨界破壞載荷的影響。Huang D 等[4]分析了混凝土結(jié)構(gòu)的靜態(tài)彈性變形和開裂行為。針對(duì)脆性材料的熱沖擊損傷，鄭州大學(xué)的王振宇[5]采用PD 方法研究了巖石材料的動(dòng)態(tài)斷裂與瞬時(shí)熱傳導(dǎo)行為。Zhang H 等[6]提出了一種考慮熱力載荷的常規(guī)態(tài)基PD 模型，以預(yù)測(cè)金屬和陶瓷雙材料結(jié)構(gòu)的斷裂損傷問(wèn)題。PD 理論適合處理復(fù)雜的裂紋萌生和擴(kuò)展問(wèn)題，根據(jù)該理論可以較好地研究脆性材料的熱沖擊損傷問(wèn)題。

該文基于PD 理論建立一種分析脆性材料熱沖擊損傷問(wèn)題的方法，推導(dǎo)近場(chǎng)動(dòng)力學(xué)的傳熱方程和熱力耦合本構(gòu)方程，模擬脆性材料在熱沖擊載荷作用下發(fā)生多裂紋隨機(jī)萌生和擴(kuò)展的現(xiàn)象。

1 鍵基PD 熱力耦合模型

鍵基PD 的運(yùn)動(dòng)方程由虛功原理推導(dǎo)而來(lái)，滿足動(dòng)量守恒定律和角動(dòng)量守恒定律，如公式（1）所示。

式中：ρ為質(zhì)點(diǎn)的密度；ü為加速度矢量；H為視界范圍；b為體積力密度矢量；f為鍵力密度矢量，與質(zhì)點(diǎn)的材料屬性、質(zhì)點(diǎn)間的拉伸情況和質(zhì)點(diǎn)的視界半徑有統(tǒng)計(jì)學(xué)意義；x為質(zhì)點(diǎn)的坐標(biāo)矢量；t為時(shí)間；u為質(zhì)點(diǎn)的位移矢量；u'為x視界內(nèi)其他質(zhì)點(diǎn)的位移矢量；x'為視界內(nèi)其他質(zhì)點(diǎn)的坐標(biāo)矢量。

鍵力密度矢量f與鍵的伸長(zhǎng)率s呈線性關(guān)系，如公式（2）所示。

式中：ξ、η分別為2 個(gè)質(zhì)點(diǎn)之間的相對(duì)位置、相對(duì)位移；c為鍵合常數(shù)；s為鍵的伸長(zhǎng)率；μ為質(zhì)點(diǎn)的位移矢量。

對(duì)各向同性材料來(lái)說(shuō)，可以通過(guò)WPD=WCM（WPD為基于近場(chǎng)動(dòng)力學(xué)的應(yīng)變能；WCM為基于連續(xù)介質(zhì)力學(xué)的應(yīng)變能）來(lái)獲得PD 參數(shù)與連續(xù)介質(zhì)力學(xué)參數(shù)間的關(guān)系[3]，鍵合常數(shù)c在平面應(yīng)變狀態(tài)下如公式（3）所示。

式中：h為厚度；δ為近場(chǎng)半徑；E為彈性模量。

引入一個(gè)標(biāo)量值函數(shù)g，以控制質(zhì)點(diǎn)間鍵力的存在和消失，如公式（4）所示。

式中：sc為鍵失效時(shí)的臨界拉伸；α為熱膨脹系數(shù)；T為溫差。

令消除所有穿過(guò)新裂紋表面的鍵所做的功等于臨界能量釋放率G，可以推導(dǎo)二維情況臨界拉伸與臨界能量釋放率間的關(guān)系，如公式（5）所示。

式中：Gc為臨界能力釋放率。

溫度是標(biāo)量，但是熱傳導(dǎo)具有方向性，質(zhì)點(diǎn)之間的溫差是熱傳導(dǎo)的驅(qū)動(dòng)力，由傅里葉熱傳導(dǎo)定律可知，熱流密度q如公式（6）所示。

式中：k為材料的熱導(dǎo)率；T為溫度梯度。

鍵基PD 傳熱方程可以用溫差、熱流密度等參數(shù)來(lái)表示，如公式（7）所示。

式中：ρ為質(zhì)點(diǎn)的密度；cv為質(zhì)點(diǎn)的比熱容；T為 溫度對(duì)時(shí)間的一階導(dǎo)數(shù)；T'為質(zhì)點(diǎn)x'的溫度；x'為x視界內(nèi)其他質(zhì)點(diǎn)的坐標(biāo)矢量；Vx'為x視界內(nèi)所有質(zhì)點(diǎn)；sb為熱源；fh為熱響應(yīng)函數(shù)。

熱響應(yīng)函數(shù)與質(zhì)點(diǎn)的溫差、熱導(dǎo)率和質(zhì)點(diǎn)的相對(duì)位置有統(tǒng)計(jì)學(xué)意義，因此熱響應(yīng)函數(shù)如公式（8）所示。

式中：τ為質(zhì)點(diǎn)的溫差；κ為微熱導(dǎo)率，可以通過(guò)假設(shè)簡(jiǎn)單線性溫度場(chǎng)下PD 熱勢(shì)能和經(jīng)典連續(xù)介質(zhì)力學(xué)熱勢(shì)能是等價(jià)的來(lái)獲得。

二維情況下近場(chǎng)動(dòng)力學(xué)的微熱導(dǎo)率如公式（9）所示。

式中：s為厚度。

對(duì)流邊界條件下的單位體積產(chǎn)熱率如公式（10）所示。

式中：n為垂直于邊界Sf的法線方向單位矢量；Sf為實(shí)體熱流邊界。

包括標(biāo)量函數(shù)g的傳熱方程如公式（11）所示。

式中：hs為單位體積的產(chǎn)熱率。

通過(guò)時(shí)間的積分可以獲得所有質(zhì)點(diǎn)的當(dāng)前溫度，然后在鍵力密度矢量中引入熱膨脹項(xiàng)c'T，以表征溫度變化與鍵力間的關(guān)系，并計(jì)算鍵力密度矢量，如公式（12）所示。

式中：c'為PD 熱力耦合參數(shù)。

c'與鍵合參數(shù)c的關(guān)系如公式（13）所示。

式中：α為熱膨脹系數(shù)。

獲得特定溫度場(chǎng)影響下的所有質(zhì)點(diǎn)鍵力密度矢量后，將其代入鍵基PD 運(yùn)動(dòng)方程，如公式（14）所示。

該方程可計(jì)算所有質(zhì)點(diǎn)的加速度ü，再積分求得所有質(zhì)點(diǎn)的位移結(jié)果。

2 數(shù)值模擬

2.1 PD 模型設(shè)置

PD模型取直徑為l3 mm、厚度為1 mm的圓形薄板建立二維模型，薄板初始溫度為300 ℃~500 ℃，溫度分布均勻，并在圓盤外周施加1 圈對(duì)流換熱邊界，環(huán)境溫度設(shè)置為T0=15 ℃。材料參數(shù)如下：E=370 GPa，υ=1/3，p=3 980 kg/m3，Gc=12.16 J·m2。

初步設(shè)置模型的質(zhì)點(diǎn)間距Δ=50 μm，視界尺寸為δ=3.015Δ。熱沖擊溫度為500℃，溫度場(chǎng)計(jì)算步長(zhǎng)設(shè)置為dt0=1×10-5s，位移場(chǎng)計(jì)算步長(zhǎng)為dt=5×10-3s，為了保證ADR收斂，設(shè)置最大收斂步數(shù)為200、收斂安全系數(shù)為5。

2.2 PD 參數(shù)

進(jìn)行δ收斂分析，以確定視界最優(yōu)值，在較細(xì)的網(wǎng)格尺寸（Δ=50 μm）下進(jìn)行，δ如公式（15）所示。

式中：δ為視界尺寸；m為視界尺寸與網(wǎng)格大小的比值，m分別取2、3、5、7、10 和15。

保持其他條件完全一致，給出水平中線上的質(zhì)點(diǎn)，在不同m值下的PD 解與FEM 解的相對(duì)誤差歐幾里德范數(shù)的均值如公式（16）所示。

式中：ai為質(zhì)點(diǎn)的待分析參數(shù)值；為對(duì)照組的參數(shù)值；h為質(zhì)點(diǎn)數(shù)量。

由表1 可知，當(dāng)m=3 時(shí)，溫度場(chǎng)和位移場(chǎng)的相對(duì)誤差最小分別為0.052 3%和0.547 5%，計(jì)算精度較高。

表1 不同m 值的PD 解相對(duì)誤差

分析網(wǎng)格大小的最優(yōu)值，將m值固定為3，以保證結(jié)果精度，令Δ為13 μm、25 μm、50 μm、100 μm、200 μm 和500 μm，其他條件完全一致，不同Δ值的PD 解與FEM 解的相對(duì)誤差歐幾里德范數(shù)的均值見表2。當(dāng)Δ分別為50 μm和100 μm 時(shí)，相對(duì)誤差較小，分別為0.547 5%和0.156 9%?？紤]計(jì)算精度和效率，認(rèn)為Δ為50 μm 或100 μm 是較好的選擇。

表2 不同Δ 值下PD 解與FEM 解的相對(duì)誤差

綜上所述，當(dāng)T=500 ℃時(shí)，m=3、Δ=50 μm 的結(jié)果精度較高。

3 計(jì)算結(jié)果分析

采用近場(chǎng)動(dòng)力學(xué)方法計(jì)算不同溫度下的熱沖擊裂紋擴(kuò)展情況，在裂紋數(shù)量和長(zhǎng)度方面，與試驗(yàn)結(jié)果[7]進(jìn)行對(duì)比，如圖1 所示。在熱沖擊載荷的作用下，陶瓷薄板發(fā)生了垂直于周向的裂紋萌生和擴(kuò)展現(xiàn)象，裂紋具有周期性和層次性。熱沖擊裂紋沿陶瓷薄板四周均勻分布，周期性地出現(xiàn)長(zhǎng)裂紋、中裂紋和短裂紋，裂紋層次鮮明，均由圓板四周萌生并向圓心處擴(kuò)展。當(dāng)溫度為300 ℃時(shí)，PD 模擬結(jié)果出現(xiàn)4 條長(zhǎng)裂紋，分布在0°、90°、180°和270°方向，分布較均勻，而試驗(yàn)結(jié)果同樣明顯出現(xiàn)4 條長(zhǎng)裂紋，數(shù)量一致，大致分布在0°、70°、190°和270°方向。當(dāng)溫度為350 ℃時(shí)，PD 模擬結(jié)果出現(xiàn)8 條長(zhǎng)裂紋，分布在0°、30°、60° 、90°、120°、150°、180°和270°方向，而試驗(yàn)結(jié)果明顯出現(xiàn)6 條長(zhǎng)裂紋，分布在0°、60° 、120°、150°、210°和270°方向，相比之下，PD模擬結(jié)果多產(chǎn)生了2 條長(zhǎng)裂紋，二者有5 條長(zhǎng)裂紋方向一致。當(dāng)溫度為400 ℃時(shí)，PD 模擬結(jié)果出現(xiàn)8 條長(zhǎng)裂紋，分布在0°、30°、60° 、90°、120°、150°、180°和270°方向，而試驗(yàn)結(jié)果明顯出現(xiàn)6 條長(zhǎng)裂紋，分布在30°、60° 、120°、190°、270°和330°，相比之下，PD 模擬結(jié)果多產(chǎn)生了2 條長(zhǎng)裂紋，二者有4 條長(zhǎng)裂紋方向一致。當(dāng)溫度為500 ℃時(shí)，PD 模擬結(jié)果出現(xiàn)8 條長(zhǎng)裂紋，分布在0°、30°、60°、90°、180°、210°、240°和270°方向，而試驗(yàn)結(jié)果同樣明顯出現(xiàn)8 條長(zhǎng)裂紋，數(shù)量一致，分布在30°、60° 、140°、190°、210°、240°、270°和340°方向，二者有5 條長(zhǎng)裂紋方向一致。綜上所述，PD 模擬結(jié)果與試驗(yàn)結(jié)果的熱沖擊裂紋趨勢(shì)較一致，PD 方法成功捕捉到了垂直于圓形薄板周向的裂紋萌生和擴(kuò)展，裂紋的周期性和層次性規(guī)律也十分相似，主裂紋的長(zhǎng)度也比較接近。但是裂紋分布情況存在一定差別，推測(cè)可能是由試驗(yàn)所用陶瓷板材料分布不均勻或存在初始缺陷導(dǎo)致，試驗(yàn)誤差引起的淬火不均勻性也會(huì)影響試驗(yàn)結(jié)果。

圖1 熱沖擊裂紋對(duì)比（左為PD 模擬結(jié)果，右為試驗(yàn)結(jié)果）

圖2 統(tǒng)計(jì)了300 ℃~500 ℃，PD 模擬結(jié)果與試驗(yàn)結(jié)果的裂紋數(shù)量對(duì)比，根據(jù)裂紋長(zhǎng)度與圓盤半徑R的比值歸類，將裂紋長(zhǎng)度在0.1R~0.5R（一種定量的對(duì)比手段，將裂紋長(zhǎng)度與圓盤半徑進(jìn)行對(duì)比，當(dāng)裂紋長(zhǎng)度等于1R時(shí)，表示裂紋貫穿；當(dāng)裂紋長(zhǎng)度0.5R~1.0R時(shí)，表示裂紋較長(zhǎng)接近貫穿；當(dāng)裂紋長(zhǎng)度0.1R~0.5R時(shí)，表示裂紋較短；當(dāng)裂紋長(zhǎng)度小于0.1R時(shí)，屬于微損傷，裂紋小且數(shù)量龐大，統(tǒng)計(jì)難度較大，因此忽略）范圍內(nèi)視為短裂紋，大于0.5R為長(zhǎng)裂紋，忽略小于0. 1R的微小裂紋。模擬結(jié)果的長(zhǎng)裂紋和短裂紋交錯(cuò)出現(xiàn)，規(guī)律性較強(qiáng)。當(dāng)熱沖擊溫度為300 ℃時(shí)，PD 模擬結(jié)果和試驗(yàn)結(jié)果的長(zhǎng)裂紋數(shù)量均為4 條，短裂紋數(shù)量均為20 條，因此總裂紋數(shù)量相等。當(dāng)熱沖擊溫度為350 ℃時(shí)，PD 模擬結(jié)果和試驗(yàn)結(jié)果的長(zhǎng)裂紋數(shù)量相差2 條，短裂紋數(shù)量均為20條，總裂紋數(shù)量較接近。同樣，當(dāng)熱沖擊溫度為400 ℃時(shí)，PD 模擬結(jié)果和試驗(yàn)結(jié)果的長(zhǎng)裂紋數(shù)量相差2 條，短裂紋數(shù)量均為23 條，總裂紋數(shù)量較接近。當(dāng)熱沖擊溫度為500 ℃時(shí)，PD 模擬結(jié)果和試驗(yàn)結(jié)果的長(zhǎng)裂紋數(shù)量均為8 條，短裂紋數(shù)量均為28 條，因此總裂紋數(shù)量相等。綜上所述，在不同熱沖擊溫度下，對(duì)陶瓷薄板熱沖擊裂紋進(jìn)行模擬，PD 方法較準(zhǔn)確，不僅短裂紋數(shù)量高度一致，而且長(zhǎng)裂紋數(shù)量偏差較小。

圖2 試驗(yàn)與PD 模擬裂紋數(shù)量

最后，對(duì)比不同溫度下的PD 模擬結(jié)果發(fā)現(xiàn)，長(zhǎng)裂紋數(shù)量由300 ℃的4 條逐漸增至500 ℃的8 條，而短裂紋數(shù)量由300 ℃的20 條逐漸增至500 ℃的28 條，隨著熱沖擊溫度的升高，長(zhǎng)裂紋和短裂紋的數(shù)量逐漸增加，主裂紋長(zhǎng)度也逐漸增加，陶瓷板的損傷程度加劇。這一現(xiàn)象說(shuō)明，熱沖擊載荷不僅會(huì)影響裂紋的數(shù)量，而且還會(huì)影響裂紋的整體長(zhǎng)度。

綜上所述，該文的PD 方法較好地模擬了陶瓷圓形薄板的熱沖擊裂紋擴(kuò)展問(wèn)題，僅使用少量計(jì)算資源就得到與試驗(yàn)較一致的現(xiàn)象，證明了PD 方法對(duì)該問(wèn)題有較高的適用性和經(jīng)濟(jì)性。

4 結(jié)語(yǔ)

該文基于鍵基PD 理論建立了熱沖擊損傷模型，并對(duì)陶瓷薄板的熱沖擊過(guò)程進(jìn)行模擬以及參數(shù)收斂性分析，保證了計(jì)算精度和效率，最后將模擬結(jié)果與試驗(yàn)結(jié)果進(jìn)行比較。

模型的視界尺寸δ和網(wǎng)格尺寸Δ對(duì)計(jì)算精度和效率的影響較大。將m值分別取為2、3、5、7、10 和15，以改變視界尺寸。當(dāng)m=3 時(shí)，PD 方法與FEM 方法的溫度和位移結(jié)果相對(duì)誤差均較小，計(jì)算時(shí)間短。網(wǎng)格尺寸Δ值分別為13 μm、25 μm、50 μm、100 μm、200 μm 和500 μm，當(dāng)Δ為50 μm或100 μm 時(shí)，相對(duì)誤差較小。因此認(rèn)為m=3、Δ=50 μm 是該模型精度高、效率高的組合。

當(dāng)熱沖擊溫度為300 ℃~500 ℃時(shí)，分別對(duì)比了PD 結(jié)果與試驗(yàn)結(jié)果的裂紋情況，發(fā)現(xiàn)無(wú)論是裂紋分布規(guī)律，還是裂紋數(shù)量、長(zhǎng)度和形式都較為一致，現(xiàn)象吻合較好，說(shuō)明PD 模型較好地模擬了陶瓷薄板的熱沖擊損傷問(wèn)題。此外，還發(fā)現(xiàn)隨著溫度升高，長(zhǎng)裂紋和短裂紋的數(shù)量均逐漸增加，長(zhǎng)裂紋在較高的熱沖擊溫度下變得更長(zhǎng)，陶瓷板的損傷程度加劇。