(1.東北石油大學 機械科學與工程學院，黑龍江大慶 163318；2.中石化石油機械股份有限公司，武漢 430205)

0 引言

工業(yè)生產(chǎn)的高參數(shù)化使得換熱器的服役環(huán)境日漸苛刻，換熱器的管殼程溫差的逐漸加大，使得單一材料的管板無法滿足工藝要求[1]，使用層合的復合材料管板[2-6]也可能會出現(xiàn)層間開裂等失效問題。相比之下，功能梯度材料在制備過程中，結(jié)構(gòu)兩側(cè)由不同性能的材料組成，但這兩種材料并不直接相連，而是在它們之間形成一個在成分、組織及性能上呈現(xiàn)連續(xù)梯度變化的過渡區(qū)域，從而達到內(nèi)部界面弱化、甚至完全消失，材料的宏觀特性(如彈性模量、熱膨脹系數(shù)等)在空間位置上呈現(xiàn)梯度變化，進而減小和克服結(jié)合部位性能的不匹配，以滿足不同結(jié)構(gòu)元件對材料使用性能的不同要求，達到提高承載結(jié)構(gòu)能力并合理使用材料的目的[7]，因此，功能梯度材料(Functionally Graded Materials，簡稱FGM)在換熱器管板上具有廣闊的應用前景[8-9]。目前，Ganczarski等[10]推導出了FGM厚板在熱載荷作用下的平面應力狀態(tài)；Ashrafi等[11-12]研究了厚功能梯度材料板的非線性動力響應和振動;李世榮等[13]研究了經(jīng)典板理論下功能梯度板的自由振動響應問題和一階剪切變形理論下功能梯度圓板與相應均勻圓板軸對稱彎曲解的線性相似轉(zhuǎn)換關(guān)系;Sui等[14]研究了四邊簡支的功能梯度層合板彎曲的一階剪切變形理論解與相應均勻板彎曲的經(jīng)典解的關(guān)系，并且給出了相應的轉(zhuǎn)換關(guān)系的表達式?；谠摬牧系牧W行為研究使得功能梯度板力學問題得到簡化，對工程上功能梯度材料管板的設計與制造具有重要意義?；诖?，本文利用物理中面的概念，將功能梯度管板中材料的各向異性問題轉(zhuǎn)化到經(jīng)典版理論中，從而求得功能梯度管板的應力分布。

1 彈性理論下功能梯度管板應力計算

當換熱工藝存在較大溫差時，通常選用U形管式換熱器來減小溫差應力，因此,功能梯度材料制備的換熱器管板基本上都用在U形管式換熱器上[15]。同時FGM管板作為復合材料管板，可能出現(xiàn)與筒體及管箱材料不同的情況，為便于連接，采用通過墊片與管箱法蘭及殼體法蘭連接形式的管板，即圖1所示的a型，這樣既便于設計，也便于加工制造。因此，以U形管式換熱器中法蘭連接且不兼做法蘭的功能梯度管板為例，進行內(nèi)壓載荷下FGM管板的應力分析。

圖1 管板的連接形式

1.1 功能梯度管板應力計算公式推導

1.1.1 冪函數(shù)模型

考慮到成分分布范圍廣泛以及材料加工時的難易程度，文中選用冪函數(shù)模型[16]。

(1)

式中fij——FGM成分分布模型函數(shù);

xc,xp——兩組分邊界區(qū)位置;

β——梯度因子。

取FGM板z向為厚度方向，板的物性參數(shù)沿厚度方向呈梯度變化，將兩種材料分別定義為材料1和材料2，根據(jù)冪函數(shù)模型表達式來表達功能梯度材料制備的換熱器管板(以下簡稱FGM管板)的物性參數(shù)[16]彈性模量、線膨脹系數(shù)如下。因目前的功能梯度材料大多是陶瓷-金屬或金屬-金屬，泊松比μ都相差不大，故可認為泊松比是與空間域無關(guān)的一個常量，即功能梯度材料的泊松比是一個常數(shù)。

彈性模量E(hz)：

(2)

線膨脹系數(shù)α(hz)：

(3)

式中hz——FGM板厚度方向的坐標，0≤hz≤h；

h——板厚度，mm;

c——下角標，材料1的物性參數(shù)；

m——下角標，材料2的物性參數(shù)。

1.1.2 基本求解方程

管板作為管殼式換熱器的重要零部件，是一塊開有多個密集小孔的圓平板?；谝粋€厚度h、半徑R的功能梯度材料圓平板，并建立O-rθz柱坐標系，原點O與板的圓心重合，O-rθ面位于全金屬相所在的平面，z軸垂直于該平面，材料沿z軸呈連續(xù)梯度變化，并作出以下假設：

(1)在FGM圓板中垂直于中性面的平面在變形后仍然垂直于中性面，即Kirchhoft假設；

(2)FGM圓板為小變形；

(3)FGM圓板的彈性模量及泊松比沿厚度方向按照體分比函數(shù)連續(xù)變化。

取FGM圓板的物理中面為hz=h0，則：

(4)

FGM圓板內(nèi)任意一點的位移在r,θ,z方向上的位移分別為u,v,w。對于周邊簡支的FGM圓板，其半徑為R，厚度為h，受到均布載荷q0，在平板上取單元體，如圖2所示。

圖2 FGM圓平板單元體

z方向所受的力的平衡方程為[17]：

-Vrrdθ+(Vr+dVr)(r+dr)dθ+q0rdrdθ=0

(5)

式中Vr——單位長度上的剪力,N;

q0——均布載荷，MPa。

根據(jù)微元體的力矩平衡條件，所有內(nèi)力與外力對圓柱體外壁面切線的力矩代數(shù)和為零，即：

(6)

式中Mr——徑向彎矩，N·m；

Mθ——周向彎矩，N·m。

根據(jù)前文假設在橫截面上的切向應變γrθ,γθz及軸向應力σz均為0，可以得出下式：

(7)

式中σr——徑向應力，MPa；

μ——泊松比；

εr——徑向應變，mm；

εθ——周向應變，mm；

σθ——周向應力，MPa；

τrθ——切向應力，MPa；

γrθ——切向應變，mm。

(8)

式中r——任意一點處的半徑，mm。

因彎矩Mr,Mθ是應力σr,σθ沿板厚方向的積分，建立彎矩Mr與Mθ與撓度w的關(guān)系如下：

(9)

對于軸對稱FGM圓平板，利用彈性理論中的位移法即可建立可解的平衡方程。聯(lián)立式(9)與式(6)可得：

(10)

聯(lián)立式(5)，式(10),可得周邊簡支的FGM圓板位移函數(shù)：

(11)

式(11)中有4個未知系數(shù)A,B,C,D，需補充如下邊界條件。

(1)因r=0時lnr無意義，所以為滿足w在r=0處存在，所以D=0。

將系數(shù)A,B,C,D代入式(11)，求得周邊簡支的FGM圓板位移函數(shù)表達式：

(12)

將式(12)代入式(9)得到彎矩方程，求得周邊簡支的FGM圓板應力表達式為：

(13)

1.1.3 FGM管板應力分布計算公式

GB/T 151—2014《熱交換器》中對管板計算的基本假設是把管板視為承受均布載荷、放置在彈性基礎上的且受管孔均勻削弱的當量圓平板，當量圓平板應力分析的主要依據(jù)是彈性力學板殼理論中的Kirchhoff薄板理論。在計算中考慮管板的整體彎曲變形和面內(nèi)拉伸變形，而忽略由于介質(zhì)壓力作用于管孔帶上產(chǎn)生的局部變形，同時將管孔周邊和管孔帶上的應力集中按應力分類原則劃分為峰值應力，在強度校核中不予考慮[18]。因此，針對U形管式換熱器中法蘭連接且不兼做法蘭的功能梯度管板的應力分析中，首先對管板進行模型簡化，簡化后的力學模型見圖3。

圖3 管板簡化力學模型

對a型管板來說，R即為墊片壓緊力作用中心圓半徑，Rt為布管區(qū)當量半徑，pt為管程內(nèi)壓，ps為殼程內(nèi)壓。從圖3中可知，簡化后的a型管板相當于周邊簡支的當量圓平板，當量圓平板由一個不布管區(qū)的圓環(huán)與一個將布管區(qū)通過管孔削弱系數(shù)簡化成的實心圓板組成。計算布管區(qū)域應力時，利用剛度削弱系數(shù)進行計算[19]。

由于梯度材料兩個相的泊松比相差不大，可以認為泊松比為常量，不隨坐標變化，令其為μ；令管板的設計壓力為q0。并且管板簡化后的力學模型中圓環(huán)和環(huán)內(nèi)的圓盤彈性模量不同，結(jié)合式(11)可知，管板的位移函數(shù)是由兩段函數(shù)構(gòu)成的一個連續(xù)函數(shù)，即：

(14)

式中η——彎曲剛度系數(shù)；

Rt——管板布管區(qū)當量半徑,mm。

位移函數(shù)除了滿足前面的邊界條件外，還應滿足下面的變形協(xié)調(diào)方程：

(15)

根據(jù)位移函數(shù)、邊界條件及變形協(xié)調(diào)方程，求出管板的應力分布：

(16)

(17)

式中μ*——強度削弱系數(shù)。

由于簡支約束在r=R處彎矩為0，因此可能存在應力最大值的點只剩下r=0和r=Rt處，兩處應力值的計算方法如下。

當r=0時：

σr=σθ

(18)

當r=Rt時：

(19)

1.2 功能梯度管板應力理論值計算

以某廠的冷卻器功能梯度管板為研究對象，采用上述推導的理論計算公式進行計算，獲取管板的理論計算應力值。該大型冷卻器管板采用功能梯度材料制造，設計參數(shù)如表1所示。

表1 管板設計參數(shù)

管板采用功能梯度材料制造，連接方式為a型連接方式，即管箱與殼體皆通過螺栓法蘭墊片與管板進行連接，殼程側(cè)采用耐高溫且隔熱性能良好的ZrO2陶瓷，管程側(cè)為022Cr17Ni12Mo2不銹鋼，中間采用梯度方式連續(xù)過渡，梯度因子β=1，管板厚度110 mm。

定義Orθ位于管程側(cè)表面，O位于管板管程側(cè)表面中心，從管程側(cè)指向殼程側(cè)方向為z軸正方向，計算均布載荷q0在管殼程壓力共同作用時危險截面應力值。將相關(guān)參數(shù)帶入式(1)中，得到管板的中性面位置h0=52.762 mm，計算系數(shù)K=1.898×1010。

根據(jù)ASME規(guī)范相應的計算得到剛度削弱系數(shù)η=0.471，強度削弱系數(shù)μ*=0.53，管板布管區(qū)當量半徑Rt=448.6 mm。

在管殼程壓力共同作用下，根據(jù)式(18),(19)計算危險點處的應力如下。

(1)在r=0處時。

σr=σθ

=-6.423×10-6×(193000-381.818z)

×(z-52.762)

經(jīng)計算，在管程側(cè)表面z=0處，σr=σθ=65.40 MPa；在殼程側(cè)表面z=110處，σr=σθ=-55.51 MPa。

(2)在r=Rt處時。

σr=-1.279×10-6×(193000-381.818z)×(z-52.762)

σθ=-1.989×10-6×(193000-381.818z)×(z-52.762)

經(jīng)計算，在管程側(cè)表面z=0處，σr=13.02 MPa，σθ=20.25 MPa;在殼程側(cè)表面z=110處，σr=-11.051 MPa,σθ=-16.76 MPa。

2 壓力作用下功能梯度管板應力有限元分析

功能梯度材料的材料參數(shù)隨位置變化的特性給其有限元模型帶來了很大困難。文中選用ABAQUS、并利用ABAQUS的USDFLD子程序來進行功能梯度管板的建模及有限元分析，通過自定義場變量來控制材料屬性設置，使材料屬性成為空間域的函數(shù)，進而實現(xiàn)功能梯度材料有限元模型的建立，使其滿足有限元分析要求[20]。根據(jù)第1.2節(jié)以及表2中參數(shù)，進行功能梯度材料管板與相關(guān)USDFLD子程序編寫。

表2 FGM管板金屬相與陶瓷相材料屬性

首先進行FGM管板的幾何建模，考慮管板的對稱性，建立管板的1/4模型，根據(jù)圣維南定理，計算選取換熱管外伸長度為100 mm，建立幾何模型如圖4所示。

在property模塊中根據(jù)表2中數(shù)據(jù)定義材料屬性，并在材料定義模塊中添加場變量及用戶自定義場變量選項。ABAQUS不能直接進行映射網(wǎng)格劃分，只能通過Use mapped meshing where appropriate選項讓程序自動判斷。在可選的劃分方式中，結(jié)構(gòu)化網(wǎng)格的質(zhì)量最高，后續(xù)計算也容易收斂。管板本身的模型復雜，不能采用結(jié)構(gòu)化網(wǎng)格劃分，但可以運用partition工具將管板沿著換熱管孔切割，分成形狀簡單的區(qū)域，使其能夠進行結(jié)構(gòu)化網(wǎng)格劃分。對管板進行合理切割，采用結(jié)構(gòu)化網(wǎng)格進行劃分，所得單元都是六面體單元，共劃分的單元數(shù)量為601 662個，節(jié)點數(shù)量為721 463個。FGM管板的有限元模型如圖5所示。

圖4 FGM管板幾何模型

圖5 管板有限元模型

在第三強度理論下分析管殼層壓力共同作用下FGM管板的等效應力、徑向應力與周向應力。選擇管板有限元模型的單元類型為C3D8，在step模塊中選擇靜力分析模塊，在load模塊中施加載荷，在管板的管程側(cè)和換熱管內(nèi)表面施加管程設計壓力0.78 MPa，在殼程側(cè)及換熱管外表面施加殼程設計壓力1.58 MPa，在管板的兩個對稱面施加對稱約束，在管板與法蘭的連接面上施加夾緊約束，在換熱管的截面上施加換熱管軸向力。進入job模塊，加載編好的USDFLD子程序并進行求解。文中主要是進行管板應力分析，因此需要略去換熱管，取管板為局部模型進行詳細的應力分析。管板的第三強度理論應力云圖見圖6，可以看出，最大值位于管板中心開孔處，為68.56 MPa。

圖6 第三強度理論管板等效應力云圖

以全局坐標系為基礎建立柱坐標系，查看管板的徑向應力和周向應力，如圖7,8所示。

圖7 管板徑向應力分布云圖

圖8 管板周向應力分布云圖

從圖中可以看出，管板的徑向應力最大值與管板周向應力最大值均位于管板中心附近的開孔處。忽略計算誤差，管板中心處的徑向應力與周向應力相等，在管程側(cè)的值為63.59 MPa，在殼程側(cè)的值為-55.15 MPa。

3 FGM管板理論解與有限元解對比分析

第1節(jié)中求出了彈性理論下FGM管板在危險點處沿z軸(即厚度方向)的應力分布函數(shù)，并計算出在管程側(cè)與殼程側(cè)表面危險點處的應力值，根據(jù)第2節(jié)的有限元分析可以提取出危險點處沿厚度方向的應力分布的有限元解及管程側(cè)與殼程側(cè)表面危險點處應力的有限元解，其結(jié)果如表3所示。

表3 管板危險點處應力

在r=0及r=Rt處分別列出管板沿厚度分布的徑向應力及周向應力的理論解析解及有限元數(shù)值解曲線。考慮到開孔削弱及隔板槽的影響，在繪制有限元解曲線時選取最大值及最小值所在處的應力分布進行繪制。圖9～11中h=0處為管程側(cè)表面，h=110處為殼程側(cè)表面。

圖9 管殼程壓力共同作用下管板中心r=0處徑向應力及周向應力分布

圖10 管殼程壓力共同作用下布管區(qū)邊緣r=Rt處徑向應力分布

圖11 管殼程壓力共同作用下布管區(qū)邊緣r=Rt處周向應力分布

圖9示出管殼程壓力共同作用下FGM管板在r=0處的徑向應力及周向應力沿管板厚度分布曲線的理論解與數(shù)值解對比圖，可以看出，在管程側(cè)表面應力值大于0，為拉應力，在殼程側(cè)應力值小于0，為壓應力，在管板的兩側(cè)表面分別取最大值和最小值。

圖10為管殼程壓力共同作用下FGM管板在r=Rt處的徑向應力沿管板厚度分布的理論解與數(shù)值解曲線對比圖，圖11為管殼程壓力共同作用下FGM管板在r=Rt處的徑向應力沿管板厚度分布的理論解與數(shù)值解曲線對比圖，可以看出,在管程側(cè)表面應力值大于0，為拉應力，在殼程側(cè)應力值小于0，為壓應力，在管板的兩側(cè)表面分別取最大值和最小值。

分析圖9,10可知，利用第1節(jié)中所推導的FGM管板應力分布計算公式所得到的理論解與有限元分析所獲得的數(shù)值解曲線幾乎重合；分析圖11可知,在r=Rt處的周向應力曲線擬合差一些，但誤差在允許范圍內(nèi)，充分說明了理論解的正確性。

4 結(jié)語

基于彈性理論，研究了功能梯度管板在壓力載荷下的應力分布情況，并且推導了相關(guān)應力分布計算公式。之后，利用ABAQUS及其USDFLD子程序建立了功能梯度管板的有限元模型，對壓力載荷下的功能梯度管板進行了靜力分析，對比分析了理論解曲線與有限元解曲線的擬合效果，驗證所推導應力分布計算公式的正確性，為功能梯度管板的設計校核提供了理論基礎。文中雖然對功能梯度管進行了彈性理論解的推導及有限元數(shù)值分析，但由于受到客觀條件的限制，筆者僅僅給出了靜力載荷下的功能梯度管板的彈性理論應力計算公式，后續(xù)可以繼續(xù)研究彈性理論下的熱應力理論解及熱力耦合載荷下的應力分布。