熱沖擊下帶裂紋功能梯度厚壁圓筒的分析

彭亦鵬， 陳愛軍

(南京理工大學 理學院，南京 210094)

熱沖擊下帶裂紋功能梯度厚壁圓筒的分析

彭亦鵬， 陳愛軍

(南京理工大學 理學院，南京 210094)

針對功能梯度材料熱沖擊下的斷裂問題，對功能梯度材料的材料參數分類進行加權處理，建立帶裂紋功能梯度厚壁圓筒的模型；給出熱量平衡方程與邊界條件，建立相關問題的有限元分析格式；改進Calahan算法以適用于帶狀稀疏矩陣的一階非齊次常微分方程組，完善溫度場的求解過程，計算了裂紋尖端的溫度應力強度因子，并且分析功能梯度材料組分構成對溫度應力強度因子的影響。以上研究為帶裂紋功能梯度厚壁圓筒的可靠性分析及結構優化設計提供了參考。

功能梯度材料；厚壁圓筒；熱沖擊；溫度應力強度因子；有限元分析；裂紋

功能梯度材料最早是作為熱格柵提出的，故而此類結構在熱環境中的響應分析一直是研究的重點。過去數年內，各國學者開展了大量的研究工作，Hui等[1]用無網格法對功能梯度材料進行了熱力學分析；陳健橋等[2]將無網格Petrov-Galerkin方法應用于功能梯度材料的三維熱傳導問題；Shen等[3]分析了功能梯度材料板在熱環境中和橫向壓力作用下的非線性彎曲；Sutradhar等[4]用邊界元法分析了功能梯度材料的瞬態熱傳導問題；Sladek等[5]利用無網格局部邊界元法，分析了功能梯度材料的瞬態熱傳導問題； Chen等[6]用精細積分法對功能梯度材料的熱傳導問題進行了靈敏度分析。Liu等[7]基于二維各向異性熱彈性理論，給出了均勻熱流下帶裂紋功能梯度材料的應力強度因子的解決辦法。Lee等[8]采用有限單元法分析了梯度板在熱彈性加載下的非線性響應問題。田建輝等[9]采用混合數值法分析了功能梯度材料板中瞬態熱響應問題。曹蕾蕾等[10]針對功能梯度材料物性參數隨坐標變化的特點，基于Trefftz完備解提出了求解功能梯度材料熱傳導問題的有限梯度元方法。針對功能梯度材料的二維瞬態熱傳導問題，藍林華等[11]提出了一種有效的降維精細積分法。李世榮等[12]基于Timoshenko梁理論研究了功能梯度材料梁在一維熱沖擊載荷作用下的瞬態動力響應。高效偉等[13]采用徑向積分邊界元法分析了功能梯度材料動態斷裂力學問題。唐雪松[14]研究了含中心裂紋的正交各向異性板，在遠場均勻熱流作用下溫度場的分布。Wang等[15]通過FE/FD方法計算了功能梯度材料的瞬態溫度場和相關的熱應力。趙軍等[16]采用攝動法推導出無限大對稱型梯度功能材料平板的一維非定常溫度場及非定常熱應力場的解析解，并計算了瞬態熱應力強度因子。Zhang等[17]通過DQM對功能梯度殼在瞬態溫度場中的熱應力進行分析。同樣采用DQM數值解法，蒲育等[18]得到了功能梯度梁自由振動的二維彈性解。杜長城等[19]研究了熱環境中功能梯度圓柱殼的內共振非線性模態。

基于現有的文獻的不足，本文首先希望能將關于功能梯度材料研究由解析解引入到試件的數值分析當中，建立一種相對自由的功能梯度厚壁圓筒力學性能參數函數，在此基礎上的編寫相應模塊的程序。與此同時，參考了呂和祥等[20]對一階齊次常微分方程組的解法，改進了Calahan算法，提供了一種行之有效的算法來計算有限元離散后的熱沖擊下功能梯度厚壁圓筒微分方程。

本文的研究首先建立功能梯度厚壁圓筒的微分方程；然后推導了改進的Calahan算法以求解一階非齊次常微分方程；之后進行數值計算，通過計算熱沖擊下帶裂紋功能梯度厚壁圓筒的溫度應力強度因子，驗證了算法的有效性與適用性；提出熱沖擊下帶裂紋功能梯度厚壁圓筒相關問題的結論。

1 熱沖擊下功能梯度厚壁圓筒的微分方程

1.1 均勻材料瞬態熱傳導問題的微分方程

首先，瞬態溫度場的溫度變量φ在直角坐標系中應滿足微分方程

(1)

即熱量平衡方程，式中的第一項是微體升溫需要的熱量；第2,3,4項是由x,y和z方向傳入微體的熱量；最后一項是微體內熱源產生的熱量。微分方程表明，微體升溫所需的熱量應與傳入微體的熱量以及微體內熱源產生的熱量相平衡。

邊界條件如下：

(2)

q(在Γ2邊界上)

(3)

h(φa-φ)(在Γ3邊界上)

(4)

1.2 功能梯度厚壁圓筒的材料特性

對于如圖1功能梯度厚壁圓筒模型的功能梯度厚壁圓筒，長度為L，中心半徑為r，厚度為2h，建立坐標系(x,θ,z)，其中，坐標原點于圓心，θ軸在圓周方向，z軸在圓筒中心方向。考慮圓筒的內表面承受熱沖擊載荷，外表面與環境做熱對流的邊界條件。

圖1 功能梯度厚壁圓筒模型

復合材料的材料特性是用混合物的線性規則來描述的。引入參數c(x)∈[0,1]。基于圖1給出的模型，可令

(5)

式中：x是半徑，c0、k和n是材料參數。在式(5)的描述下，通過選取適當的n和k，可以實現功能梯度材料特性范圍較大的非線性連續變化。

對于材料的溫度特性，即熱傳導系數kx、kθ，比熱容Cx、Cθ，密度ρx、ρθ，用P來表示。以兩種材料為例，P可以表示為

P=Pic(r)+Pj[1-c(r)]

(6)

式中，Pi、Pj表示這兩種材料的溫度特性。

在小變形假設下，應變與位移的關系是

(7)

式中：u是徑向位移；εx,εθ是應變。

依據混合物的線性規則以及兩種線彈性材料的均勻假設，功能梯度厚壁圓筒的應力由下式給出

(8)

其中，

(9)

下標0,1表示兩種不同的材料，λi、μi是拉梅系數。式(8)表示混合物的線性規則應用于拉梅系數得出的結果。

由Xin[21]可以得出，如此描述功能梯度材料的等效性能，對試件模型的計算精度有很大提升，也更加符合實際情況。

1.3 邊界條件的形成

在如圖1所示模型中，假設厚壁圓筒的外部與外界熱傳導，內部受到熱沖擊荷載。

溫度的初始條件和外部、內部邊界條件是

(10)

(11)

式中：ΔT是熱載荷振幅;a是載荷變化的常數;hr是厚壁圓筒與外界的換熱系數。

1.4 有限元格式的建立

由于功能梯度材料斷裂問題十分復雜，從理論上給出問題的解析解耗時耗力，不具有通用的意義。而實際工程應用當中，功能梯度材料制備的多樣性，使材料模型與荷載的形成各式各樣，增加了理論分析的難度。通過運用有限單元法編寫程序，求解功能梯度材料的斷裂問題，能夠在工程中得到廣泛的應用。厚壁圓筒模型可以簡化為平面問題，在程序的編寫過程當中應用四邊形八節點等參單元。

對于形如式(1)的微分方程，建立等效積分形式

(12)

w=w2=w3=δφ

(13)

將式(13)代入式(12)，并對其中第1個域Ω內積分的第2-4項進行分布積分，則可得到

(14)

利用式(12)可以建立瞬態溫度場有限元的一般格式，首先將空間域Ω離散為有限個單元體，在典型單元內溫度φ仍可以近似地用節點溫度φi插值得到，但要注意此時的結點溫度是時間的函數，即

(15)

插值函數Ni只是空間域的函數，應具有插值函數的基本性質，將式(15)代入式(14)，并考慮到dφi的任意性，就可以得到用來確定n個結點溫度φi的有限元求解方程

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

至此，已將時間域和空間域的偏微分問題在空間域內離散為N個結點溫度φi(t)的常微分方程的初值問題，對于給定溫度值的邊界Γ1上n1個結點，方程組始終的相應項應引入的條件是

(24)

式中，i是Γ1上ni個結點的編號。

2 一階非齊次常微分方程組具有三階精度的無條件穩定的隱式算法

將式(6)所確定的材料系數代入式(16)中，可以得到以時間為變量的一階非齊次常微分方程組

(25)

本節中提出一個新算法，使Calahan算法適用于矩陣K為稀疏帶狀矩陣，也就是對一般的非齊次常微分方程組提出一個具有三階精度，無條件穩定自起步的算法，并且能夠不破壞原稀疏矩陣的存儲方式及其運算規則。

當C是非奇異矩陣時，可分解為

C=LU

(26)

其中L是下三角陣，U是上三角陣。當C是對稱矩陣時U=LT。將式(26)代入式(25)，并前乘L-1，同時令

ψ=Uφ

(27)

則方程(25)化為標準形式

(28)

其中

(29)

2.1 自起步、三階精度的逐步積分法

令ψ1=ψ(tk+f),ψ0=ψ(tk)，由(28)有

(30)

其中，

(31)

將式(31)代入(29)，有

(32)

由式(32)解出ψ1，有

(33)

調整T使得式(33)左端的系數為平方的形式，因此令

(34)

由式(34)，可以求得

(35)

由式(35)，式(34)變為

(36)

通過移項與化簡，得到在原來帶狀陣上運算的遞推公式

φ1=(C+fUK)-1C(C+fUK)-1×

(37)

式(37)所有的運算都可以在稀疏矩陣存儲的方式上運行，其中所有的逆矩陣只是一個形式上的寫法，因為所要求的并不是逆矩陣本身，而是它與列向量的乘積，即可以用稀疏矩陣和列向量想成運算規則，用解線性方程組的方法得到。

2.2 穩定性分析

把方程(37)用于單自由度方程

則有

x1=Hx0

(38)

其中

(39)

顯然只有

(40)

由式(39)計算得

3 熱沖擊下帶裂紋功能梯度厚壁圓筒的數值分析

3.1 應力強度因子

應力強度因子作為斷裂力學的重要參量，其量綱是[力]×[長度]-3/2,是表征裂尖應力場強度的物理量，一般以符號K表示。由于非均質線彈性體裂紋尖端漸進場與均質線彈性體裂紋尖端漸進場相同，因此由裂紋應力強度因子與裂紋張開位移之間的關系可以得出應力強度因子的表達式，為

(41)

(42)

式(41)中，KI，KII和KIII分別稱為Ⅰ型、Ⅱ型和Ⅲ型裂紋的應力強度因子，μtip為裂紋尖端剪切模量，r為裂紋面上考慮的點到裂紋尖端的距離，Δui(r),i=1,2,3為局部坐標系下裂紋張開位移，其表達式為

(43)

3.2 程序的設計

對于程序的設計，旨在通過有限元的工具，對模型進行計算，從而得出較為精確的結果。

首先是通過商業軟件Ansys進行建模，劃分網格。導出單元信息Elements.txt和單元信息Nodes.txt。模型不考慮機械荷載。Materials計算功能梯度材料中的材料參數，并編入與單元等容的數組中。Boundary TemperatureConditions用以輸入溫度邊界條件，并離散為列陣。Ske和Stiff函數用以計算單元熱傳導矩陣與總體熱傳導矩陣。在單元剛度矩陣的計算中，通過InverseMatrix來計算Jacobi矩陣的逆陣。CaMatrix用以計算熱容矩陣。在形成熱容矩陣C，熱傳導矩陣K和溫度載荷列陣P之后，用如Caption2的方法計算結點溫度列陣φ。隨后，設計Strain函數來計算結點應變，StressFactorsIntense函數以計算溫度應力強度因子，輸出為Result.txt。并用線性擬合方法進行后處理，并在Matlab得到圖形形式的結果。程序設計框圖,見圖2。

3.3 溫度場的計算

在下文中，考慮由SiC和Ni構成的功能梯度厚壁圓筒。成分的材料參數在表1中給出，泊松比為0.2。

在這一節中，功能梯度厚壁圓筒內部邊界為全陶瓷，外部邊界為全金屬，組分比N1/N2=k，功能梯度厚壁圓筒的長度L=10 m，R=6 m，r=2 m。熱沖擊載荷為T0=300 K，ΔT=100 K以及a=10，熱交換系數是hr=100。

使用ANSYS劃分網格，網格與裂紋尖端顯示于圖3。圖4中是功能梯度厚壁圓筒的內壁，外壁與中心的溫度場。

表1 材料參數

圖2 程序框圖

圖3 模型網格

圖4 功能梯度厚壁圓筒內壁、外壁與中心溫度

顯然，內壁溫度隨著熱載荷而迅速增加。系統溫度隨響應時間而趨于穩態，這與李世榮等[12]的研究保持一致。

3.5 應力與應力強度因子的計算

由上節求得的溫度場求得應變，依照彈性力學的本構關系進一步可以得到功能梯度厚壁圓筒的應力。

為了驗證上一部分中改進的Calahan算法(Improved Calahan Method,ICM)的可靠性，先做一組對比分析。分析各向同性的均勻材料Cu作為厚壁圓筒材料，模型尺寸、溫度荷載與上一節相同。

從圖5和圖6中可以看出，對于以Cu為材料的厚壁圓筒，ICM求得的結果與有限元分析軟件Ansys得出的結果相近。當系統達到穩態后，誤差在0.03%以內。改進的Calahan算法也將應用到接下來的計算中。

圖5 Cu為材料的厚壁圓筒瞬態熱應力σxx

圖6 Cu為材料的厚壁圓筒瞬態熱應力σθθ

圖7和圖8展示了上節所示的功能梯度厚壁圓筒的徑向和切向應力。可以看出，隨著厚壁圓筒的溫度場趨于穩態，不同位置的應力也趨于定值。

圖7 不同表面的瞬態熱應力σxx

圖8 不同表面的瞬態熱應力σθθ

圖9展示了在熱載荷變化時候的裂紋尖端應力強度因子。

圖9 T=0.25ΔT,0.5ΔT,ΔT時應力強度因子

在T=0.5ΔT時，變化式(5)中n的值，通過運算得到了不同組分下，溫度應力強度因子隨時間變化的結果展示于圖10。

圖10 不同組分下的溫度應力強度因子

4 結 論

本文基于線彈性理論，對內壁受熱沖擊載荷的厚壁圓筒進行了分析。首先，基于微分方程等效積分形式的伽遼金提法在空間域進行有限元離散，得到一階常微分方程組，并引入改進的Calahan算法進行溫度場的求解。研究的主要結論是：

(1) 熱應力隨著熱沖擊激增，峰值熱應力發生在功能梯度厚壁圓筒的內、外表面。與溫度場不同的是，應力在不同表面收斂到不同的定值。

(2) 當功能梯度厚壁圓筒受到熱沖擊時，功能梯度材料組分連續變化的特性使其在很短的時間內緩和溫度場的峰值。然而，當溫度趨于穩定時，由于功能梯度材料不均勻的特性，殘余環向應力仍較大。因此殘余環向應力不能忽略。

(3) 隨著熱沖擊溫度的變化，即ΔT/T逐漸增大，溫度應力強度因子依次變大，并逐漸趨于穩定，在達到穩定時比較接近。

(4) 組分變化對溫度應力強度因子也有影響，隨著組分變化，溫度應力強度因子也隨之變化。

[1] HUI W, QIN H Meshless approach for thermomechanical analysis of functionally graded materials[J]. Engineering Analysis with Boundary Elements, 2008, 32: 704-712.

[2] 陳健橋，丁亮.功能梯度材料瞬態熱傳導問題的MLPG方法[J]. 華中科技大學學報(自然科學版)，2007, 35(4): 119-121.

CHEN Jianqiao, DING Liang, A MLPG method of transient heat transference in FGMs[J]. Journal of Huazhong University of Science and Technology(Nature Science), 2007, 35(4): 119-121.

[3] SHEN H S. Nonlinear bending response of functionally graded plates subjected to transverse loads and in thermal environments[J]International Journal of Mechanical Sciences 2002, 44: 561-584.

[4] SUTRADHAR A, PAULINO G H. The simple boundary element method for transient heat conduction in functionally graded materials[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193: 4511-4539.

[5] SLADEK J, SLADEK V, ZHANG C. Transient heat conduction analysis in functionally graded materials by the meshless local boundary integral equation method[J]. Computational Materials Science, 2003,28: 494-504.

[6] CHEN B S, TONG L Y. Sensitivity analysis of heat conduction for functionally graded materials[J]. Materials and Design, 2004,25: 663-672.

[7] LIU L, KARDOMATEAS G A. A dislocation approach for the thermal stress intensity factors of a crack in an infinite anisotropic medium under uniform heat flow[J]. Composites:Part A Appled Science and Manufacturing, 2006,37: 989-996.

[8] LEE S J, REDDY J N. Non-linear response of laminated composite plates under thermomechanical loading[J]. International Journal of Non-Linear Mechanics, 2005, 40(7): 971-985.

[9] 田建輝,韓旭,孫小衛. 基于混合數值法的功能梯度材料板瞬態熱響應分析[J]. 固體力學學報, 2008,29(4): 396-401.

TIAN Jianhui, HAN Xu, SUN Xiaowei. Transient thermal response of functionally graded material plates based on the hybrid numerical method[J]. Chinese Journal of Solid Mechanics, 2008, 29(4): 396-401.

[10] 曹蕾蕾, 裴建中, 張春國,等. 基于混合Trefftz有限梯度元法的功能梯度材料熱傳導數值模擬[J]. 中國科技論文, 2015,10(2): 222-226.

CAO Leilei, PEI Jianzhong, ZHANG Chunguo, et al. A hybrid Trefftz graded finite element method to steady-state heat conduction in functionally graded materia[J]. China Sciencepaper, 2015,10(2): 222-226.

[11] 藍林華，富明慧，程正陽.功能梯度材料瞬態熱傳導問題的降維精細積分法[J].固體力學學報，2010,31(4): 406-410.

LAN Linhua, FU Minghui, CHENG Zhengyang. Decrement-dimensional precise time integration of 2-D transient heat conduction equation for functionally graded materials[J]. Acta Mechanica Solida Sinica, 2010, 31(4): 406-410.

[12] 李世榮，范亮亮.功能梯度梁在熱沖擊下的動態響應[J].振動工程學報，2009, 22(4): 371-378.

LI Shirong, FAN Liangliang. Dynamic responses of functionally graded material beams under thermal shock[J]. Journal of Vibration Engineering, 2009, 22(4): 371-378.

[13] 高效偉, 鄭保敬, 劉健. 功能梯度材料動態斷裂力學的徑向積分邊界元法[J]. 力學學報, 2015, 47(5): 868-873.

GAO Xiaowei, ZHENG Baojing, LIU Jian. Dynamic fracture analysis of functionally graded materials by radial integration bem[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015,47(5): 868-873.

[14] 唐雪松. 含中心裂紋正交各向異性板均勻熱流作用下的溫度場解析解[J]. 工程力學, 2007,24(3): 28-41.

TANG Xuesong. Analytical solutions of temperature fields for an orthotropic plate with a central crack under remote uniform heat flows[J]. Engineering Mechanics, 2007,24(3): 28-41.

[15] WANG B L, MAI Y W, ZHANG X H. Thermal shock resistance of functionally graded materials[J]. Acta Materialia, 2004,52: 4961-4973.

[16] 趙軍, 王朝霞, 董一芬,等. 對稱型梯度功能材料的瞬態熱應力強度因子[J]. 山東大學學報(工學版), 2004, 34(3): 9-13.

ZHAO Jun, WANG Chaoxia, DONG Yifen, et al. Transient thermal stress intnsity factors for a functionally gradient material with symmetrical structure[J]. Journal of Shandong University(Engineering Science), 2004, 34(3): 9-13.

[17] ZHANG J H, LI G Z, LI S R, et al. DQM based thermal stresses analysis of a functionally graded cylindrical shell under thermal shock[J]. Journal of Thermal Stresses, 2015, 38(9): 959-982.

[18] 蒲育，滕兆春.Winkler-Pasternak彈性地基FGM梁自由振動二維彈性解[J]. 振動與沖擊, 2015, 34(20): 74-79.

PU Yu, TENG Zhaochun. Two-dimensional elasticity solutions for free vibration of FGM beams resting on Winkler-Pasternak elastic foundations[J]. Journal of Vibration and Shock, 2015, 34(20): 74-79.

[19] 杜長城，李映輝，金學松.熱環境中功能梯度圓柱殼的內共振非線性模態[J]. 振動與沖擊, 2014, 33(6): 161-165.

DU Changcheng, LI Yinghui, JIN Xuesong. Nonlinear normal modes of functionally graded cylindrical shells with internal resonance in thermal environment[J]. Journal of Vibration and Shock, 2014, 33(6): 161-165.

[20] 呂和祥，董志強，陳建峰. 一階常微分方程組的一個有效解法[J]. 大連理工大學學報, 2006, 46(5): 629-632.

Lü Hexiang, DONG Zhiqiang, CHEN Jianfeng. Effective solution method for the first-order ordinary differential equation system[J]. Journal of dalian university of Technology, 2006, 46(5): 629-632.

[21] XIN L B, DUI G S, YANG S Y, et al. An elasticity solution for functionally graded thick-walled tube subjected to internal pressure[J]. International Journal of Mechanical Sciences, 2014(89): 334-349.

Numerical simulation for a functionally graded thick-walled cylinder with cracks under thermal shock

PENG Yipeng, CHEN Aijun

(School of Sciences, Nanjing University of Science & Technology, Nanjing 210094, China)

Aiming at the fracture problem of functionally graded materials under thermal shock, physical parameters of functionally graded materials were classified and weighted to establish a model of a functionally graded thick-walled cylinder with cracks. Its heat balance equations and boundary conditions were deduced to build the relevant problem’s finite element analysis procedure. Calahan algorithm was improved to solve the first order inhomogeneous ordinary differential equations of the corresponding banded sparse matrix and obtain the cylinder’s temperature field and stress field. Then the thermal stress intensity factor of the crack tip was calculated, the effects of functionally graded material’s components on the thermal stress intensity factor were analyzed. The results provided a reference for reliability analysis and structural optimal design of a functionally graded thick-walled cylinder with cracks under thermal shock.

functionally graded materials; thick-walled cylinder; thermal shock; thermal stress intensity factor; finite element analysis; crack

江蘇省自然科學基金(BK20150766)

2016-02-19 修改稿收到日期：2016-06-11

彭亦鵬 男，碩士生，1991年3月生

陳愛軍 男，博士，副教授，1972年12月生

O316

A

10.13465/j.cnki.jvs.2017.15.010