孫 燊 易 敏

(南京航空航天大學(xué) 航空航天結(jié)構(gòu)力學(xué)及控制全國(guó)重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210016)

0 引言

可重復(fù)使用運(yùn)載器具有高效、經(jīng)濟(jì)、環(huán)保等特點(diǎn),是目前航天技術(shù)發(fā)展的重要方向[1]。2015年以來(lái),美國(guó)私人航天公司SpaceX回收并重復(fù)使用獵鷹9號(hào)運(yùn)載火箭,國(guó)外研究機(jī)構(gòu)不斷突破,以求擺脫單次使用運(yùn)載火箭高昂費(fèi)用對(duì)空間探索的制約。我國(guó)一方面密切追蹤國(guó)外先進(jìn)技術(shù)發(fā)展[2-3],另一方面努力推進(jìn)自身關(guān)鍵技術(shù)攻關(guān)。可重復(fù)使用運(yùn)載器技術(shù)既能降低空間運(yùn)輸成本,貫徹可持續(xù)發(fā)展戰(zhàn)略,又是彰顯綜合國(guó)力、搶占空間優(yōu)勢(shì)的必要工具。液體火箭發(fā)動(dòng)機(jī)作為動(dòng)力裝置,其耐久性和可靠性是保障運(yùn)載器可重復(fù)使用的關(guān)鍵。推力室喉部結(jié)構(gòu)是液體火箭發(fā)動(dòng)機(jī)的核心部件,在服役過(guò)程中會(huì)經(jīng)歷復(fù)雜、交變的熱機(jī)載荷作用,不均勻的溫度分布、局部的塑性應(yīng)變更會(huì)加快結(jié)構(gòu)的疲勞失效。因此,對(duì)推力室喉部結(jié)構(gòu)進(jìn)行傳熱分析和疲勞壽命預(yù)測(cè)具有重要意義。

推力室在工作時(shí)會(huì)經(jīng)歷急速的升溫與降溫,在局部會(huì)發(fā)生較大的熱變形,分析結(jié)構(gòu)承受的熱載荷需要計(jì)算結(jié)構(gòu)的傳熱情況和溫度分布。推力室內(nèi)壁需要承受高溫燃?xì)獾臎_刷作用,冷卻劑在冷卻通道中逆流降溫。多種材料分別具有不同的物性,對(duì)傳熱分析的計(jì)算結(jié)果提出很高的要求。吳峰等應(yīng)用湍流模型進(jìn)行了三維數(shù)值模擬,計(jì)算了再生冷卻推力室通道的流動(dòng)與傳熱[4]。韓煒對(duì)推力室的再生冷卻情況進(jìn)行了三維數(shù)值模擬,得到了內(nèi)壁面上溫度分布一般規(guī)律[5]。吳有亮等提出了推力室準(zhǔn)二維傳熱計(jì)算的通用方法,相較于三維計(jì)算,降低了計(jì)算時(shí)間,極大地提升了計(jì)算效率[6]。根據(jù)傳熱分析得到的溫度分布可以計(jì)算結(jié)構(gòu)中的熱本征應(yīng)變,可仿真模擬燃?xì)鉁囟群腿細(xì)鈮簭?qiáng)載荷的共同作用效果,并預(yù)估推力室使用壽命。孫冰等對(duì)推力室結(jié)構(gòu)進(jìn)行了三維瞬態(tài)熱分析,并根據(jù)應(yīng)力應(yīng)變分布預(yù)測(cè)了結(jié)構(gòu)危險(xiǎn)點(diǎn)的失效情況[7-9]。徐紹桐等建立了推力室壁在燃?xì)?冷卻液溫度和燃?xì)鈮簭?qiáng)載荷下的二維彈塑性計(jì)算模型,提高了失效行為模擬的計(jì)算速度[10]。張明等通過(guò)數(shù)值計(jì)算研究了變截面冷卻通道的傳熱情況[11]。目前,推力室疲勞壽命預(yù)測(cè)的研究主要是基于損傷力學(xué)計(jì)算損傷的累積,或者根據(jù)應(yīng)力-應(yīng)變法計(jì)算材料的疲勞壽命,較少?gòu)臄嗔蚜W(xué)出發(fā)預(yù)測(cè)壽命。

傳統(tǒng)方法處理斷裂問(wèn)題往往將界面作為一個(gè)突變的量來(lái)處理,突變問(wèn)題的數(shù)值奇異性增加了問(wèn)題的求解難度。在材料斷裂的界面,數(shù)值上難以進(jìn)行微分算法,材料參數(shù)需要考慮因材料破壞而損失的剛度。相場(chǎng)模型通過(guò)引入在界面處急劇變化但連續(xù)的相場(chǎng)變量——序參量來(lái)描述不同的相[12-14],斷裂序參量在整個(gè)計(jì)算區(qū)域上是連續(xù)變化的,因此相場(chǎng)模型中不再需要采用特殊的數(shù)學(xué)方法追蹤界面的幾何形態(tài),也避免了追蹤斷裂面所引起的誤差。相場(chǎng)斷裂方法可以反映加載過(guò)程中斷裂序參量的演化,直觀地觀察到裂紋的演化。Francfort等[15]通過(guò)對(duì)Griffith斷裂理論[16]及變分公式的研究,提出了準(zhǔn)靜態(tài)脆性斷裂的相場(chǎng)模型。Miehe等提出了應(yīng)變張量的譜分解方法,用以區(qū)分壓縮和拉伸對(duì)應(yīng)的應(yīng)變能,并且推廣到熱-彈-塑性耦合的斷裂問(wèn)題[17-19]。文獻(xiàn)[20-21]通過(guò)修正斷裂相場(chǎng)模型中的斷裂韌性,實(shí)現(xiàn)對(duì)疲勞斷裂問(wèn)題的相場(chǎng)模擬。

本文使用相場(chǎng)疲勞斷裂模型對(duì)再生冷卻推力室喉部結(jié)構(gòu)進(jìn)行有限元模擬,通過(guò)熱傳導(dǎo)方程考慮溫度載荷引入的熱本征應(yīng)變,得到結(jié)構(gòu)的溫度和應(yīng)力/應(yīng)變分布情況,根據(jù)相場(chǎng)序參量的變化分析結(jié)構(gòu)的疲勞斷裂過(guò)程,研究推力室結(jié)構(gòu)在服役過(guò)程中的失效行為和熱機(jī)疲勞壽命。

1 熱力耦合疲勞相場(chǎng)模型

1.1 熱力耦合模型

銑槽式再生冷卻推力室采用等肋寬的薄壁結(jié)構(gòu),其喉部自?xún)?nèi)向外分別為內(nèi)壁、冷卻通道和外壁,冷卻通道使用肋片進(jìn)行強(qiáng)化傳熱。推力室中實(shí)際發(fā)生的熱量傳遞包括:燃?xì)鈱?duì)內(nèi)壁的加熱、冷卻劑對(duì)冷卻通道的對(duì)流換熱。考慮熱載荷對(duì)受力的單向耦合,采用穩(wěn)態(tài)熱平衡方程對(duì)溫度場(chǎng)進(jìn)行傳熱分析,熱傳導(dǎo)的控制方程為無(wú)熱源的導(dǎo)熱微分方程,即

(1)

式中:cp為定壓比熱容;ρ為材料密度;T為結(jié)構(gòu)溫度;t為時(shí)間;?為哈密頓算子;λ為材料熱導(dǎo)率。

發(fā)動(dòng)機(jī)推力室在工作過(guò)程中同時(shí)受到溫度載荷和機(jī)械載荷的作用。結(jié)構(gòu)的總應(yīng)變?chǔ)舤otal由熱應(yīng)變?chǔ)舤h、彈性應(yīng)變?chǔ)舉la、塑性應(yīng)變?chǔ)舙la和蠕變應(yīng)變?chǔ)與r構(gòu)成,即

εtotal=εth+εela+εpla+εcr

(2)

熱應(yīng)變與溫度的變化成線性關(guān)系,即

εth=αthΔTI

(3)

式中:αth為線膨脹系數(shù);ΔT為相對(duì)參考溫度的溫度變化量;I為二階單位張量。

1.2 塑性模型

在分離彈性應(yīng)變和塑性應(yīng)變時(shí),使用von Mises屈服準(zhǔn)則來(lái)判斷材料是否達(dá)到屈服,屈服條件表達(dá)式為

(4)

式中:f為屈服函數(shù);s為偏應(yīng)力張量;a為背應(yīng)力張量,即屈服面中心;σy為屈服應(yīng)力。

達(dá)到屈服面后塑性應(yīng)變產(chǎn)生,塑性應(yīng)變沿著屈服面梯度方向增加,即

(5)

式中:dλ為塑性乘子,即等效塑性應(yīng)變?cè)隽?σ為應(yīng)力張量。

1.3 蠕變應(yīng)變模型

計(jì)算蠕變應(yīng)變時(shí),考慮陳化理論,即Norton定律。在蠕變過(guò)程中,時(shí)效、擴(kuò)散和恢復(fù)等因素影響蠕變行為,其中時(shí)間是最大的影響因素。Norton模型的表達(dá)式為

(6)

1.4 斷裂相場(chǎng)模型

定義相場(chǎng)序參量c表示材料的斷裂情況,如果材料完好未損壞,則c為 1;如果存在裂紋,則c為 0。因此,變量c可以看作是損傷力學(xué)模型中的損傷參數(shù)。斷裂模型中總自由能密度由彈性應(yīng)變能、塑性應(yīng)變能和斷裂界面能組成,表達(dá)式為

ψ=ψela+(1-c)2ψpla+ψc

(7)

假定彈性應(yīng)變中只有拉伸分量驅(qū)動(dòng)裂紋擴(kuò)展,將彈性應(yīng)變譜分解為

εela=PΛPT

(8)

式中:Λ=diag(λ1,λ2,λ3)是與彈性應(yīng)變張量相似的對(duì)角陣,λ1、λ2、λ3為彈性應(yīng)變的特征值;P為特征矩陣。

根據(jù)譜分解,彈性應(yīng)變的正負(fù)分量為

εela+=PΛ+PT

εela-=PΛ-PT

(9)

彈性應(yīng)變能密度分解為驅(qū)動(dòng)裂紋擴(kuò)展的拉伸部分ψela+和不可驅(qū)動(dòng)裂紋擴(kuò)展的壓縮部分ψela-,表達(dá)式為

(10)

根據(jù)Borden的塑性斷裂相場(chǎng)模型[22],塑性應(yīng)變能密度表達(dá)式為

(11)

裂紋的相場(chǎng)演化方程表示為

(12)

式中:η為動(dòng)力學(xué)系數(shù);l0為界面寬度尺寸;G為臨界斷裂能密度;ψh為驅(qū)動(dòng)裂紋擴(kuò)展的歷史變量,即

(13)

1.5 考慮疲勞的相場(chǎng)模型

推力室結(jié)構(gòu)在工作過(guò)程中會(huì)經(jīng)歷周期性的熱應(yīng)力或熱應(yīng)變,單次循環(huán)載荷作用下,應(yīng)力水平可能遠(yuǎn)低于使材料斷裂的臨界應(yīng)力,但是仍會(huì)經(jīng)歷往復(fù)循環(huán)的塑性變形。因此,需要考慮循環(huán)載荷下的相場(chǎng)疲勞斷裂模型。通過(guò)引入疲勞退化函數(shù),對(duì)材料的斷裂韌性進(jìn)行修正,使其在循環(huán)載荷下不斷衰減,模擬結(jié)構(gòu)在疲勞載荷下的斷裂失效。采用的退化函數(shù)為

(14)

(15)

(16)

疲勞斷裂的相場(chǎng)演化方程表示為

(17)

2 計(jì)算結(jié)果與討論

推力室喉部為環(huán)形對(duì)稱(chēng)結(jié)構(gòu),選取1/2冷卻再生通道截面作為計(jì)算模型,結(jié)構(gòu)的幾何特征與邊界條件如圖1所示。

圖1 模擬結(jié)構(gòu)幾何示意圖Fig.1 Geometry of the structure

本文在 MOOSE(Multiphysics Object-Oriented Simulation Environment)框架下進(jìn)行模擬,通過(guò)編寫(xiě)有限元程序,利用PETSc求解器和Libmesh庫(kù)進(jìn)行并行計(jì)算。為了提高計(jì)算效率,使用二維平面應(yīng)變結(jié)構(gòu)進(jìn)行模擬。本文采用四邊形網(wǎng)格進(jìn)行網(wǎng)格劃分。模型網(wǎng)格基礎(chǔ)尺寸為0.1 mm,在肋片和內(nèi)壁進(jìn)行加密,尺寸最小為0.01 mm,網(wǎng)格數(shù)量為16 305。對(duì)于網(wǎng)格無(wú)關(guān)性驗(yàn)證,進(jìn)行了6 495、16 305和29 888不同單元數(shù)量的測(cè)試,3種網(wǎng)格的溫度分布基本一致,第一種網(wǎng)格與后兩種網(wǎng)格所獲得的應(yīng)力/應(yīng)變分布有明顯差別,后兩種網(wǎng)格所獲得的結(jié)果基本一致,故本文結(jié)果皆是基于第二種網(wǎng)格(16 305單元)得到的。計(jì)算采用單向耦合求解,在每一個(gè)時(shí)間步,先進(jìn)行傳熱計(jì)算,得到溫度分布后,將溫度轉(zhuǎn)化為熱應(yīng)變,代入力學(xué)方程求解應(yīng)力/應(yīng)變和斷裂相場(chǎng)值,然后進(jìn)行下一時(shí)間步的計(jì)算。模型的兩邊設(shè)置為對(duì)稱(chēng)邊界。發(fā)動(dòng)機(jī)工作過(guò)程分為啟動(dòng)-熱試車(chē)-后冷-關(guān)機(jī)4個(gè)階段,各階段時(shí)長(zhǎng)分別為5、300、20、600 s,環(huán)境溫度為295.15 K。結(jié)構(gòu)的溫度邊界條件和壓力邊界條件參見(jiàn)文獻(xiàn)[23]。

結(jié)構(gòu)的內(nèi)壁和肋片材料為NARloy-Z銅合金,外壁材料為1Cr21Ni5Ti。NARloy-Z銅合金性能參數(shù)參照Esposito 在高低溫環(huán)境下的單軸拉伸試驗(yàn)數(shù)據(jù)[24],蠕變模型使用Ellis得到的不同溫度、不同應(yīng)力水平下的蠕變?cè)囼?yàn)數(shù)據(jù)[25]。相場(chǎng)疲勞斷裂模型的參數(shù)見(jiàn)表1。為了獲得相場(chǎng)模型參數(shù),根據(jù)材料的屈服強(qiáng)度、抗拉強(qiáng)度和伸長(zhǎng)率等參數(shù)先進(jìn)行不同參數(shù)下的單軸拉伸計(jì)算測(cè)試,選取能實(shí)現(xiàn)吻合實(shí)驗(yàn)結(jié)果的模型參數(shù),再進(jìn)行推力室受載計(jì)算。

表1 NARloy-Z合金參數(shù)Tab.1 Parameters of NARloy-Z alloy

2.1 溫度分布結(jié)果

不同階段下結(jié)構(gòu)的溫度分布如圖2所示。在熱試車(chē)階段,結(jié)構(gòu)的升溫比較明顯,最高溫度出現(xiàn)在內(nèi)壁下表面,G點(diǎn)溫度最高可達(dá)881.17 K。

圖2 各階段溫度分布Fig.2 Temperature distribution at different stages

選取結(jié)構(gòu)上7個(gè)參考點(diǎn)的溫度變化曲線,如圖3所示?？梢钥闯?自下而上,內(nèi)壁溫度最高,肋片次之,外壁溫度最低。由于內(nèi)壁與肋片的材料熱導(dǎo)率高于外壁,外壁上參考點(diǎn)的升溫速率也明顯慢于其他兩部分。外壁上B、C兩點(diǎn)的溫度基本一致。

圖3 各參考點(diǎn)溫度變化Fig.3 Temperature evolution of reference points

2.2 單次工作結(jié)構(gòu)應(yīng)力/應(yīng)變分析

推力室單次工作各階段結(jié)構(gòu)的應(yīng)力/應(yīng)變分布如圖4和圖5所示。在開(kāi)機(jī)階段,結(jié)構(gòu)升溫幅度較小,整體的壓應(yīng)力和壓應(yīng)變的數(shù)值較低。由于肋片膨脹受到內(nèi)、外壁的限制,整體處于壓應(yīng)力水平。在熱試車(chē)階段,結(jié)構(gòu)整體升溫,受熱膨脹,由于受到對(duì)稱(chēng)邊界的限制,結(jié)構(gòu)整體表現(xiàn)為受壓。肋片膨脹受到外壁限制,C處會(huì)產(chǎn)生較高的壓力。該階段下,結(jié)構(gòu)內(nèi)壁下表面溫度最高,G點(diǎn)受到左側(cè)邊界的限制,受到的壓應(yīng)變最大,最高達(dá)到1.34%。在后冷階段,結(jié)構(gòu)彈性變形恢復(fù),內(nèi)壁由于壓載較大,降溫后有殘余的塑性變形。內(nèi)壁下表面F點(diǎn)和G點(diǎn)分別表現(xiàn)為殘余拉應(yīng)變和殘余壓應(yīng)變。F點(diǎn)受到拉伸載荷,殘余應(yīng)變約為1%。在關(guān)機(jī)階段,結(jié)構(gòu)的溫度分布與后冷階段相差不大,結(jié)構(gòu)最終的受載情況和殘余應(yīng)變與后冷階段基本一致。

圖4 各階段應(yīng)力分布Fig.4 Stress distribution at different stages

圖5 各階段應(yīng)變分布Fig.5 Strain distribution at different stages

根據(jù)應(yīng)力和應(yīng)變分布圖,確定D、E、F、G這4點(diǎn)在工作過(guò)程中受到的局部應(yīng)力-應(yīng)變較大,為服役危險(xiǎn)點(diǎn)。這4點(diǎn)的應(yīng)力、應(yīng)變曲線如圖6所示。4個(gè)點(diǎn)在循環(huán)中都經(jīng)歷壓應(yīng)力變?yōu)槔瓚?yīng)力,D、E、G點(diǎn)始終受到壓應(yīng)變,F點(diǎn)最終由壓應(yīng)變轉(zhuǎn)變?yōu)槔瓚?yīng)變。

圖6 單次工作中4個(gè)特征點(diǎn)的應(yīng)力-應(yīng)變曲線Fig.6 Stress-strain curves of four points during single cycle

2.3 10次循環(huán)結(jié)構(gòu)應(yīng)力/應(yīng)變分析

對(duì)10個(gè)工作循環(huán)下的推力室喉部結(jié)構(gòu)進(jìn)行仿真。根據(jù)單次循環(huán)應(yīng)力應(yīng)變分析,結(jié)構(gòu)的危險(xiǎn)點(diǎn)位于內(nèi)壁下表面的F和G點(diǎn),循環(huán)仿真后重點(diǎn)關(guān)注F和G點(diǎn)的應(yīng)力應(yīng)變情況。F點(diǎn)和G點(diǎn)10次循環(huán)的應(yīng)力應(yīng)變曲線分別如圖7和圖8所示。隨著循環(huán)數(shù)的增加,F、G兩點(diǎn)表現(xiàn)為應(yīng)力控制的非對(duì)稱(chēng)循環(huán)加載,并且每次循環(huán)結(jié)束殘余應(yīng)變值不斷增加。其中,F點(diǎn)為拉伸應(yīng)變,G點(diǎn)為壓縮應(yīng)變。F、G兩點(diǎn)殘余應(yīng)變發(fā)展曲線如圖9所示。

圖7 F點(diǎn)10次循環(huán)應(yīng)力-應(yīng)變響應(yīng)Fig.7 Stress-strain response of point F in 10 cycles

圖8 G點(diǎn)10次循環(huán)應(yīng)力-應(yīng)變響應(yīng)Fig.8 Stress-strain response of point G in 10 cycles

圖9 F、G點(diǎn)10次循環(huán)殘余應(yīng)變的變化Fig.9 Residual strain evolution at F and G in 10 cycles

2.4 結(jié)構(gòu)疲勞失效分析

由于單次循環(huán)結(jié)構(gòu)承受的最大應(yīng)力小于材料斷裂破壞的臨界應(yīng)力,結(jié)構(gòu)中的斷裂序參量增長(zhǎng)較小。循環(huán)加載并計(jì)算結(jié)構(gòu)的斷裂序參量,直至結(jié)構(gòu)中出現(xiàn)序參量值達(dá)到0.95則停止計(jì)算。最終發(fā)現(xiàn),91個(gè)循環(huán)后結(jié)構(gòu)失效,計(jì)算停止,F點(diǎn)斷裂序參量最先達(dá)到臨界值c=0.95,表明結(jié)構(gòu)的疲勞壽命為91次循環(huán),失效點(diǎn)位于內(nèi)壁下表面F點(diǎn)。

單次循環(huán)后的序參量分布如圖10所示。圖10表明,結(jié)構(gòu)中內(nèi)壁下表面中心點(diǎn)F序參量最大,在單次工作結(jié)束后,F點(diǎn)的損傷最大,但此時(shí)序參量的值還處于很低的水平,對(duì)結(jié)構(gòu)的承載能力影響不大,不會(huì)使應(yīng)力衰減,仍可以進(jìn)行多次循環(huán)工作。結(jié)構(gòu)在內(nèi)壁下表面受到的溫度載荷最高,并且F點(diǎn)由于對(duì)稱(chēng)邊界受到約束,在溫度恢復(fù)常溫后,受到的拉應(yīng)力最大,所以序參量演化最快。根據(jù)應(yīng)力分布云圖,在結(jié)構(gòu)變截面處,即C、D點(diǎn)會(huì)出現(xiàn)應(yīng)力集中。C點(diǎn)受到的溫度載荷較低,所以溫度載荷恢復(fù)后受到的載荷水平較低,沒(méi)有出現(xiàn)損傷。D點(diǎn)受到的溫度載荷較高,相場(chǎng)值為6×10-7,但與F點(diǎn)相比仍較低,因此結(jié)構(gòu)并未在變截面處出現(xiàn)裂紋。

圖10 單次工作后相場(chǎng)序參量分布Fig.10 Distribution of phase-field order parameter after 1 cycle

10次循環(huán)后的序參量分布如圖11所示,結(jié)果與單次循環(huán)相似,仍在F點(diǎn)序參量值最高。

圖11 10次工作后相場(chǎng)序參量分布Fig.11 Distribution of phase-field order parameter after 10 cycles

結(jié)構(gòu)中的斷裂序參量演化如圖12所示。隨著循環(huán)數(shù)增加,下表面的斷裂序參量值不斷上升,并且EF邊界明顯快于內(nèi)壁其他區(qū)域,從F點(diǎn)到E點(diǎn)從大到小基本呈線性分布。

圖12 相場(chǎng)序參量演化Fig.12 Evolution of phase-field order parameter

結(jié)構(gòu)失效破壞時(shí)的徑向位移如圖13所示。疲勞破壞時(shí),EF邊界明顯向下凹陷,內(nèi)壁不斷變薄,最終在F點(diǎn)處斷裂失效。

圖13 結(jié)構(gòu)破壞時(shí)的徑向位移Fig.13 Radial displacement at failure

3 結(jié)論

本文采用熱-彈-塑性耦合的疲勞斷裂相場(chǎng)模型,分析了火箭發(fā)動(dòng)機(jī)再生冷卻推力室結(jié)構(gòu)在循環(huán)運(yùn)行過(guò)程中的溫度、應(yīng)力和應(yīng)變分布及變化,根據(jù)斷裂序參量的演化,分析了熱機(jī)疲勞破壞行為并計(jì)算了結(jié)構(gòu)的疲勞壽命,主要結(jié)論如下。

1)推力室結(jié)構(gòu)在工作過(guò)程中熱量自下而上傳導(dǎo),溫度逐漸降低。結(jié)構(gòu)升溫降溫產(chǎn)生的塑性變形是導(dǎo)致結(jié)構(gòu)疲勞失效的主要因素。

2)結(jié)構(gòu)內(nèi)壁部分變形較大,隨著循環(huán)次數(shù)增加,F點(diǎn)和G點(diǎn)分別不斷累積殘余正應(yīng)變和殘余壓應(yīng)變。

3)F點(diǎn)受到殘余正應(yīng)變積累,驅(qū)動(dòng)裂紋擴(kuò)張的拉伸應(yīng)變能不斷增加,最終導(dǎo)致結(jié)構(gòu)發(fā)生疲勞破壞。內(nèi)壁中心受到殘余正應(yīng)變,使內(nèi)壁向下表面塌陷,壁面變薄。