可重復(fù)使用液體火箭發(fā)動機設(shè)計參數(shù)對推力室身部棘輪應(yīng)變的影響

羽中豪，金 平，蔡國飆

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

1 引言

可重復(fù)使用運載器是未來航天領(lǐng)域的重要發(fā)展方向，而能夠長期穩(wěn)定、可重復(fù)工作的液體火箭發(fā)動機，是研制可重復(fù)使用運載器的關(guān)鍵技術(shù)。在可重復(fù)使用發(fā)動機循環(huán)工作過程中，發(fā)動機推力室身部處于高低溫交變載荷之下，會產(chǎn)生極大的塑形應(yīng)變，并隨著循環(huán)次數(shù)逐漸積累，形成棘輪效應(yīng)，直至發(fā)生破壞［1］。為了在設(shè)計階段延緩熱疲勞與棘輪效應(yīng)并提高發(fā)動機推力室壽命，研究發(fā)動機推力室設(shè)計參數(shù)對推力室在工作過程中應(yīng)變發(fā)展情況的影響十分必要。

國內(nèi)外的許多學(xué)者都對液體火箭發(fā)動機推力室的工作狀態(tài)和熱疲勞現(xiàn)象通過仿真實驗進(jìn)行了研究［2-4］，并且對推力室進(jìn)行了基于棘輪應(yīng)變與蠕變應(yīng)變的疲勞壽命預(yù)估［5-6］。德國的Riccius等開展了大量熱-機械疲勞平板實驗研究［7-9］，以驗證彈塑性模型在發(fā)動機推力室疲勞計算中的準(zhǔn)確性［10］；程誠研究了不同冷卻通道的排布形式對推力室疲勞壽命的影響［11］。

但這些研究多是針對已成型的發(fā)動機，較少學(xué)者關(guān)注從設(shè)計層面如何避免或減輕熱棘輪現(xiàn)象的產(chǎn)生。為此，本文在分析火箭發(fā)動機推力室在單循環(huán)、多循環(huán)工作過程中的溫度場及結(jié)構(gòu)場的基礎(chǔ)上，針對推力室的設(shè)計過程，研究室壓、推力和混合比的變化對工作過程中推力室產(chǎn)生的棘輪應(yīng)變及其發(fā)展的影響。

2 物理模型和計算方法

2.1 物理模型

物理模型采用液體火箭發(fā)動機設(shè)計方法進(jìn)行設(shè)計，推力室型面與冷卻通道的具體設(shè)計流程與準(zhǔn)則主要參考文獻(xiàn)［12］和［13］，并參考國內(nèi)某型號氫氧發(fā)動機。銑槽式再生冷卻推力室包括圓筒段、收斂段以及擴張段短噴管三部分，由帶槽道的NARloy-Z合金內(nèi)壁與電鍍鎳合金外壁釬焊而成。以液氫作為冷卻劑，采用單向逆流布局，短噴管采用最大推力噴管設(shè)計。實際計算中取二分之一個冷卻通道作為計算的物理模型，如圖1所示。

圖1 推力室冷卻通道物理模型Fig．1 Physical model of thrust chamber and cooling channel

其中，模型邊界上A點為外壁中心點，B點為槽底中心點，C點為槽的氣壁中心點，D點為筋的氣壁中心點，F(xiàn)、E點分別為槽的上下角點，G點為槽頂中心點。

2.2 計算模型

2.2.1 傳熱分析部分

采用準(zhǔn)二維的方法對冷卻通道工作階段的傳熱過程進(jìn)行分析。將冷卻通道沿軸向分為若干段，每一段的傳熱過程都視為一維傳熱模型，段與段之間依靠冷卻劑傳熱與沿程損失相聯(lián)系，得到冷卻通道沿軸線的溫度與壓強分布。準(zhǔn)二維的方法在保證了精度的同時極大地提高了計算效率，同時易于進(jìn)行集成計算，計算采用的控制方程與具體的實施過程參見文獻(xiàn)［14］。

2.2.2 熱-結(jié)構(gòu)耦合分析部分

采用順序耦合法，首先分析推力室的溫度場。對于無熱源的穩(wěn)態(tài)溫度場，熱傳導(dǎo)方程為式（1）［11］：

其中，W為熱源函數(shù)；λ為內(nèi)壁導(dǎo)熱系數(shù)；ΔT為溫度梯度。燃?xì)鈧?cè)與冷卻劑側(cè)的邊界條件為第三類邊界條件［14］，即式（2）：

其中，q為熱流密度；hf為流體對壁面的對流換熱系數(shù)；Tf為流體主流溫度；Tw為內(nèi)壁溫度；n為壁面外法線方向矢量。外壁為絕熱邊界條件［2］，即式（3）：

以溫度場為基礎(chǔ)進(jìn)行結(jié)構(gòu)分析，在準(zhǔn)靜態(tài)過程中結(jié)構(gòu)場的平衡方程為式（4）［11］：

其中，［KT（δ）］為切線剛度矩陣；δ為節(jié)點位移向量，與總應(yīng)變ε的關(guān)系由柯西公式確定；R為節(jié)點機械載荷向量；Rth為節(jié)點熱載荷向量?？倯?yīng)變由熱應(yīng)變εth、彈性應(yīng)變εe、塑性應(yīng)變εp以及蠕變應(yīng)變 εcr四部分組成，即式（5）［2］：

對于各向同性的材料，熱應(yīng)變是溫度變化率的函數(shù)，具體為式（6）［6］：

其中，a為材料的線膨脹系數(shù)，ΔT為相對于參考溫度的溫度變化，I2為二階單元張量。彈性和塑性應(yīng)變的計算依靠滿足Von Mises屈服準(zhǔn)則與隨動強化準(zhǔn)則的本構(gòu)關(guān)系，本構(gòu)方程采用Chaboche 模型［15］，即式（7）：

其中，α為背應(yīng)力；M為階數(shù)，本文取3；dp為塑性應(yīng)變變化率的絕對值；Ci與γi為材料參數(shù)，由式（8）確定［16］：

因為推力室工作過程中的蠕變行為多是在較高應(yīng)力作用下發(fā)生的，所以蠕變準(zhǔn)則采用牛頓準(zhǔn)則［6］的一個改型，即指數(shù)形式的蠕變模型，如式（9）：

其中，σ為等效應(yīng)力，由應(yīng)力張量確定；C1、C2和C3為材料參數(shù)。

3 計算結(jié)果與討論

3.1 推力室單循環(huán)熱-結(jié)構(gòu)計算結(jié)果分析

以現(xiàn)有液體火箭發(fā)動機推力、室壓和混合比設(shè)計范圍［17］內(nèi)適中位置的推力室為例，對其在循環(huán)工作過程中溫度、應(yīng)變情況進(jìn)行分析，設(shè)計推力150 t、室壓16 MPa、混合比為6的發(fā)動機推力室開展熱-結(jié)構(gòu)仿真。推力室在工作過程中分為預(yù)冷、熱試、后冷、松弛四個階段，各階段時間持續(xù)如表1，計算得到的溫度場分布如圖2所示。

表1 各工作階段時間Table 1 Time of each work state

圖2 各工作階段溫度分布Fig．2 Temperature distribution of different states

圖2 表明推力室壁在工作過程中溫度變化較大，并且存在溫度分布不均的現(xiàn)象。結(jié)構(gòu)邊界上各點溫度隨時間變化的曲線如圖3所示。

圖3 單次循環(huán)工作溫度變化曲線Fig．3 Temperature during one cycle work

從圖3中可以看到，由于內(nèi)外壁材料傳熱性質(zhì)的差異，外壁溫度的變化總是滯后于內(nèi)壁溫度的變化；而內(nèi)壁材料傳熱性質(zhì)較好，內(nèi)壁各點溫度時序上保持一致。由于冷卻劑的對流換熱系數(shù)很大，冷卻通道上各點（B、E、F和G）的溫度在整個工作過程中保持一致；燃?xì)鈧?cè)C、D兩點在熱試階段相差5 K，其他階段溫度完全相同。

在整個工作過程中，推力室內(nèi)壁的溫度變化最高超過700 K，并且是發(fā)生在不到1 s的時間內(nèi)，所以內(nèi)壁將產(chǎn)生嚴(yán)重的熱變形；同時因為內(nèi)外壁溫度、材料屬性差異較大，膨脹（或收縮）程度不同，內(nèi)外壁彼此擠壓（或拉伸），導(dǎo)致內(nèi)壁材料屈服，產(chǎn)生塑性應(yīng)變。內(nèi)壁燃?xì)鈧?cè)冷卻通道中點的C點在整個工作過程中徑向應(yīng)變、切向應(yīng)變以及等效應(yīng)變?nèi)鐖D4所示。

圖4 C點的等效應(yīng)變、切向應(yīng)變與徑向應(yīng)變Fig．4 Equivalent strain， tangential strain and radial strain of point C

推力室進(jìn)入預(yù)冷階段，溫度從室溫降到冷卻劑溫度時，內(nèi)壁材料收縮，C點切向產(chǎn)生拉應(yīng)變，徑向產(chǎn)生壓應(yīng)變。隨著內(nèi)壁溫度穩(wěn)定，外壁溫度才開始下降，收縮并對內(nèi)壁產(chǎn)生徑向拉應(yīng)力，內(nèi)壁徑向壓應(yīng)變、切向拉應(yīng)變逐漸減小，直至外壁溫度也逐漸穩(wěn)定，推力室達(dá)到預(yù)冷溫度的穩(wěn)態(tài)。隨后熱試階段開始，內(nèi)壁的最高溫度在0.6 s內(nèi)上升至825 K，內(nèi)壁材料大幅膨脹，C點切向被擠壓產(chǎn)生壓應(yīng)變，徑向產(chǎn)生拉應(yīng)變，同時依舊因為外壁溫度變化延遲，內(nèi)壁切向與徑向的應(yīng)變會小范圍恢復(fù)。25 s時熱試結(jié)束，進(jìn)入后冷階段，內(nèi)壁的應(yīng)變變化與預(yù)冷類似，但由于在熱試階段產(chǎn)生了無法恢復(fù)的塑性應(yīng)變，所以后冷階段C點的切向應(yīng)變與徑向應(yīng)變均比預(yù)冷階段的大。最后推力室進(jìn)入松弛階段，不再產(chǎn)生新的應(yīng)變。

圖5為推力室各點徑向應(yīng)變在工作過程中的變化曲線，內(nèi)壁上B點與D點應(yīng)變的趨勢與C點表現(xiàn)一致。但因為D點距冷卻通道較遠(yuǎn)，整體的應(yīng)變幅值大于C點；而B點在冷卻通道上，所以熱載荷小于C點與D點，應(yīng)變的幅值也較小。推力室外壁熱載荷較小，所以在熱試階段產(chǎn)生的應(yīng)變主要由于內(nèi)壁的擠壓，與內(nèi)壁應(yīng)變方向相反。

圖5 單次循環(huán)工作各點徑向應(yīng)變Fig．5 Radial strain evaluation during one cycle work

圖6 為C點熱應(yīng)變、蠕變應(yīng)變、彈性應(yīng)變和塑形應(yīng)變隨時間變化的曲線，可以看到，熱應(yīng)變與機械應(yīng)變的方向相反，并且都遠(yuǎn)大于蠕變應(yīng)變。

3.2 推力室多循環(huán)熱-結(jié)構(gòu)計算結(jié)果分析

可重復(fù)使用發(fā)動機在使用中會進(jìn)行多次循環(huán)工作，工作過程中推力室最先出現(xiàn)破壞的點是冷卻通道中點C［1］。 因此研究推力150 t、室壓16 MPa、混合比為6的發(fā)動機的推力室冷卻通道中點C點在循環(huán)過程中的切向應(yīng)變隨時間變化。由計算結(jié)果知推力室在多次循環(huán)后應(yīng)變變化趨于穩(wěn)定，所以取前10次的循環(huán)結(jié)果進(jìn)行分析，如圖7所示。

圖6 C點的熱應(yīng)變、機械應(yīng)變與蠕變應(yīng)變Fig．6 Thermal strain， mechanic strain and creep strain of point C

圖7 循環(huán)工作時C點的切向應(yīng)變Fig．7 Tangential strain of point C during cycle works

從圖7中可以看到，因為每次循環(huán)工作的載荷相同，所以C點的應(yīng)變幅值也相同；但因為每次循環(huán)推力室內(nèi)壁都屈服、產(chǎn)生塑形應(yīng)變而形成棘輪效應(yīng)，隨著循環(huán)次數(shù)增加塑形應(yīng)變累積，每次松弛結(jié)束C點的應(yīng)變逐漸增大，推力室結(jié)構(gòu)變形逐漸嚴(yán)重并最終發(fā)生破壞。

棘輪應(yīng)變是衡量推力室棘輪效應(yīng)的重要指標(biāo)［10］，圖8為C點的棘輪應(yīng)變隨循環(huán)次數(shù)變化的曲線，圖中中點法如式（10）：

其中，εr，n為第n次循環(huán)后推力室產(chǎn)生的棘輪應(yīng)變；εmiddle，n－1、εmiddle，n分別為第 n － 1 次和第 n次循環(huán)中最大應(yīng)變與最小應(yīng)變的平均值。

圖8 循環(huán)工作時的棘輪應(yīng)變Fig．8 Ratchet strain during cycle works

由殘余應(yīng)變確定的棘輪應(yīng)變?nèi)缡剑?1）：

其中，εremain，n－1，εremain，n分別為第 n － 1 次和第n次循環(huán)后的殘余應(yīng)變。

從圖8中可以看到兩種方法計算得到的棘輪應(yīng)變在前幾次循環(huán)有所差異，而之后的差異較小，可以認(rèn)為兩種方法是等效的，因此后文中全部采用殘余應(yīng)變計算棘輪應(yīng)變。

為研究工作方式對推力室身部應(yīng)變的影響，對兩種不同工作方式的發(fā)動機推力室進(jìn)行分析，如圖9和表2所示。

圖9 兩種不同工作方式下推力室的應(yīng)變Fig．9 Tangential strain of chambers with two different work plans

表2 兩種不同的工作方式Table 2 Two different work plans

圖9中，方式A為循環(huán)工作10次、但每次熱試階段只持續(xù)10 s的發(fā)動機推力室身部的切向應(yīng)變情況，方式B為進(jìn)行一次工作、熱試階段持續(xù)100 s的發(fā)動機推力室身部的切向應(yīng)變情況。從圖中可以看到，由于存在蠕變現(xiàn)象，以方式B工作結(jié)束后推力室的應(yīng)變略大于以方式A工作第一次循環(huán)后的應(yīng)變，但因為蠕變應(yīng)變較小，所以差異不大；當(dāng)以方式A工作兩次循環(huán)后，推力室產(chǎn)生的應(yīng)變就大于方式B，其后因為棘輪效應(yīng)以方式A工作10次的推力室應(yīng)變將比方式B大4.8×10－4。綜上，雖然兩種工作方式的總熱試時間一致，但相比于一次工作，多次循環(huán)的工作方式將使推力室產(chǎn)生更大的變形。

3.3 推力室設(shè)計參數(shù)對推力室棘輪效應(yīng)的影響

為在設(shè)計階段就避免或減弱上述棘輪效應(yīng)，開展設(shè)計參數(shù)對棘輪應(yīng)變影響的研究。推力室的設(shè)計過程中，室壓、推力與混合比決定了推力室的構(gòu)型與工作狀態(tài)，因此以下將研究這三個主要設(shè)計參數(shù)與棘輪應(yīng)變之間的關(guān)系。設(shè)計參數(shù)的研究范圍主要參考國內(nèi)外現(xiàn)有液體火箭發(fā)動機的設(shè)計參數(shù)［17］?？芍貜?fù)使用火箭發(fā)動機作為一級發(fā)動機，推力通常在50 t至200 t之間。室壓的選取綜合考慮燃燒效率和冷卻、強度的影響，變化范圍取10～22 MPa。對于混合比，氫氧推進(jìn)劑的最佳混合比為8，但考慮燃燒室的冷卻問題，氫氧火箭發(fā)動機混合比通常取6左右，并且變化范圍不大［13］。在研究混合比對棘輪應(yīng)變的影響時，混合比的變化上限不應(yīng)超過最佳混合比，因此取7.5；混合比過低則會損失較多推力，取5.5作為混合比的變化的下限。

3.3.1 室壓對推力室棘輪應(yīng)變的影響

推力110 t、混合比6、室壓范圍10～22 MPa的7臺發(fā)動機推力室的棘輪應(yīng)變隨循環(huán)次數(shù)變化的曲線如圖10所示。從圖中可以看到，隨著室壓升高，推力室的棘輪應(yīng)變增大，一方面是內(nèi)壁兩側(cè)流體壓力升高的結(jié)果，另一方面是熱載荷發(fā)生變化的結(jié)果。

表3為熱試階段不同室壓推力室的熱載荷情況，其中Q為熱流密度，Twg為推力室內(nèi)壁溫度?？梢钥吹剑覊焊蟮耐屏κ覠彷d荷更重，這是因為，一方面，推力不變時提升設(shè)計室壓會使推力室整體尺寸減小，不利于再生冷卻；另一方面，更高的室壓會使燃料更充分地燃燒，推力室內(nèi)的燃?xì)鉁囟纫簿碗S之升高。

圖10 不同室壓推力室循環(huán)工作時的棘輪應(yīng)變Fig．10 Ratchet strain evaluation of chambers with different pressure during cycle works

表3 不同室壓推力室的熱載荷Table 3 Thermal load of chambers with different pressures

其次，從圖10中還能看到，較高室壓的推力室棘輪應(yīng)變隨循環(huán)次數(shù)增加呈減小趨勢，而較低室壓的推力室棘輪應(yīng)變隨循環(huán)次數(shù)增加保持不變，甚至略微呈增加趨勢。

當(dāng)循環(huán)工作的熱載荷完全相同時，每次循環(huán)產(chǎn)生的棘輪應(yīng)變也應(yīng)相同。但是對于室壓較高的推力室而言，第一次循環(huán)就使結(jié)構(gòu)產(chǎn)生了較大的變形，所以盡管其后每次循環(huán)的載荷相同，但總變形的增大會使其后每次變形的積累越來越小，即棘輪應(yīng)變呈現(xiàn)減小趨勢。反之，對于室壓較低的推力室而言，結(jié)構(gòu)殘余應(yīng)變較小甚至總應(yīng)變?yōu)榉捶较颍啈?yīng)變即會維持不變或是呈現(xiàn)上升趨勢。3.3.2 推力對推力室棘輪應(yīng)變的影響

圖11、圖12分別為室壓16 MPa、混合比6、推力范圍50～200 t的7臺發(fā)動機推力室的殘余應(yīng)變和棘輪應(yīng)變隨循環(huán)次數(shù)變化的曲線?？梢钥吹?，推力對推力室身部應(yīng)變的影響較為復(fù)雜：推力較小的推力室在開始工作時殘余應(yīng)變較小，但隨著循環(huán)次數(shù)的增加，殘余應(yīng)變增長較快，即棘輪應(yīng)變一直較大；推力較大的推力室在開始工作時殘余應(yīng)變較大，但其后增長較慢，即棘輪應(yīng)變迅速減小。

圖11 不同推力推力室循環(huán)工作時的殘余應(yīng)變Fig．11 Remain strain evaluation of chambers with different thrust during cycle works

圖12 不同推力推力室循環(huán)工作時的棘輪應(yīng)變Fig．12 Ratchet strain evaluation of chambers with different thrust during cycle works

從圖11和圖12中還可以看到，較高推力的發(fā)動機推力室第一次循環(huán)應(yīng)變較大，200 t時殘余應(yīng)變達(dá)到了5.2×10－4，因此之后循環(huán)積累的應(yīng)變大幅減少，棘輪應(yīng)變由1.0×10－4降至0.4×10－4。而推力較低的發(fā)動機推力室第一次循環(huán)殘余應(yīng)變較低，50 t時殘余應(yīng)變?yōu)?.0×10－4，之后每次循環(huán)積累的應(yīng)變基本不變，棘輪應(yīng)變維持在0.6 ×10－4左右。

推力對推力室熱環(huán)境的影響見表4。室壓相同的情況下，更大推力的推力室結(jié)構(gòu)尺寸更大，冷卻劑的流量更大，冷卻效果較好，熱流較大，推力室壁的最高溫度較低。而因為推力室身部的機械載荷來源于燃?xì)馀c冷卻劑的擠壓，保持室壓不變、改變推力對流體的壓強影響不大，所以推力對棘輪應(yīng)變的影響主要是由影響熱環(huán)境實現(xiàn)的。

表4 不同推力推力室的熱載荷Table 4 Thermal load of chambers with different thrusts

3.3.3 混合比對推力室棘輪應(yīng)變的影響

圖13為推力110 t、室壓16 MPa、混合比范圍5.5～7.5的7臺發(fā)動機推力室的棘輪應(yīng)變隨循環(huán)次數(shù)變化的曲線，表5為熱試階段不同混合比推力室的熱載荷情況。

圖13 不同混合比推力室循環(huán)工作時的棘輪應(yīng)變Fig．13 Ratchet strain evaluation of chambers with different mixture ratio during cycle works

表5 不同混合比推力室的熱載荷Table 5 Thermal load of chambers with different mixture ratios

實際發(fā)動機為保證工作時燃燒處于富燃狀態(tài)，所取混合比往往小于最佳混合比，因此較高混合比的發(fā)動機燃?xì)馊紵浞?，推力室?nèi)溫度也更高。同時，采用再生冷卻的發(fā)動機往往以燃料作為冷卻劑，較高混合比的發(fā)動機冷卻劑的流量較小，冷卻熱流較小，不利于推力室的冷卻。所以從圖13中可以看到，混合比較小的發(fā)動機推力室的棘輪效應(yīng)較輕。

其次，不同混合比的發(fā)動機推力室在工作過程中棘輪的發(fā)展情況與室壓對棘輪發(fā)展的影響是類似的：混合比較大時，推力室的溫度變化范圍較大，棘輪應(yīng)變隨循環(huán)次數(shù)增加的下降趨勢更明顯；當(dāng)混合比較小時，溫度變化范圍較小，棘輪應(yīng)變即會維持不變或是呈現(xiàn)上升趨勢。

4 結(jié)論

本文通過設(shè)計不同室壓、不同推力以及不同混合比的多臺發(fā)動機，開展熱-結(jié)構(gòu)仿真分析，對比研究了室壓、推力及混合比對推力室循環(huán)工作過程中產(chǎn)生的棘輪應(yīng)變的影響，具體如下：

1）總的熱試時間相同，多次循環(huán)、單次工作熱試時間短的發(fā)動機比一次循環(huán)、單次工作熱試時間長的發(fā)動機的推力室身部產(chǎn)生更嚴(yán)重的變形。

2）室壓、推力及混合比對推力室棘輪應(yīng)變的影響主要是由于改變了推力室的熱環(huán)境而造成的，在設(shè)計時選取更高的室壓、更小的推力或是更高的混合比會使推力室冷卻效果變差，進(jìn)而使推力室在循環(huán)工作過程中產(chǎn)生更嚴(yán)重的棘輪應(yīng)變。

3）推力室的棘輪應(yīng)變總是隨著循環(huán)工作次數(shù)的增加而趨于穩(wěn)定的；對于高室壓、大推力或是高混合比的發(fā)動機，推力室的棘輪應(yīng)變隨循環(huán)工作次數(shù)的增加會有明顯的下降趨勢；對于低室壓、小推力或是低混合比的發(fā)動機，推力室的棘輪應(yīng)變隨循環(huán)工作次數(shù)的增加變化較小。

參考文獻(xiàn)（References）

［1］ Hannum N，Kasper H，Pavli A.Experimental and theoretical investigation of fatigue life in reusable rocket thrust chambers［C］ ／／Aiaa／sae Propulsion Conference， New York， 1976：685-709.

［2］ 楊進(jìn)慧，陳濤，金平，等.可重復(fù)使用火箭發(fā)動機再生冷卻槽失效分析［J］.北京航空航天大學(xué)學(xué)報，2013，39（9）：1187-1191.Yang J， Chen T， Jin P， et al.Failure analysis of reusable rocket engine coolant passage［J］.Journal of Beijing University of Aeronautics ＆ Astronautics， 2013， 39（9）：1187-1191.（in Chinese）

［3］ 孫冰，宋佳文.液體火箭發(fā)動機推力室壁瞬態(tài)加載三維熱結(jié)構(gòu)分析［J］. 推進(jìn)技術(shù)，2016，37（7）：1328-1333.Sun B，Song J W.Three dimensional transient loading thermomechanical analysis of LRE thrust chamber wall［J］.Journal of Propulsion Technology， 2016， 37（7）：1328-1333. （in Chinese）

［4］ Amakawa1 H，Negishi H，Nishimoto M，et al.Numerical investigation of longer life combustion chambers of liquid rocket engines based on coupled thermal-fluid-structure simulation［C］ ／／Aiaa／asme／sae／asee Structures， Structural Dynamics， and Materials Conference， Grapevine， Texas， 2017：1989-2004.

［5］ 孫冰，丁兆波，康玉東.液體火箭發(fā)動機推力室內(nèi)壁壽命預(yù)估［J］. 航空動力學(xué)報， 2014，29（12）：2980-2986.Sun B，Ding Z B，Kang Y D.Life prediction of liquid rocket engine thrust chamber liner wall［J］.Journal of Aerospace Power， 2014， 29（12）：2980-2986. （in Chinese）

［6］ Asraff A K，Sunil S，Muthukumar R，et al.New concepts in structural analysis and design of double walled LPRE thrust chambers［C］ ／／Aiaa／asme／sae／asee Joint Propulsion Conference ＆ Exhibit， Sacramento， California， 2013：4368-4381.

［7］ Riccius J R，Jayaganesan B.A TMF panel optimization strategy applied to the FLPP storable engine hot gas wall geometry［C］ ／／Aiaa／sae／asee Joint Propulsion Conference， Orlando，F(xiàn)lorida， 2013：4071-4081.

［8］ Riccius J R，Hilsenbeck M，Haidn O.Optimization of geometric parameters of cryogenic liquid rocket combustion chambers［C］ ／／Aiaa／asme／sae／asee Joint Propulsion Conference＆ Exhibit， Salt Lake City， Ut， 2013：3408-3420.

［9］ Thiede G，Zametaev E B，Riccius J R，et al.Comparison of damage parameter based finite element fatigue life analysis results to combustion chamber type TMF panel test results［C］ ／／Aiaa／asme／sae／asee Joint Propulsion Conference， Orlando， Florida， 2015：4070-4082.

［10］ Riccius J R， Bouajila W， Zametaev E B.Comparison of Finite Element analysis and experimental results of a combustion chamber type TMF panel test［C］ ／／Aiaa／asme／sae／asee Joint Propulsion Conference， Son Jose， California， 2013： 3846-3860.

［11］ 程誠，王一白，劉宇，等.冷卻流道布局對銑槽噴管低周疲勞壽命的影響［J］.推進(jìn)技術(shù)，2013，34（9）：1257-1265.Cheng C，Wang Y B，Liu Y，et al.Effects of coolant passage layout on low cycle fatigue life of milled channel nozzle［J］.Journal of Propulsion Technology， 2013， 34（9）：1257-1265.（in Chinese）

［12］ 張其陽.液體火箭發(fā)動機推力室結(jié)構(gòu)與冷卻設(shè)計［D］.北京：清華大學(xué)，2012.Zhang Q Y.Structure and Cooling Design of Liquid Rocket Engine Thrust Chamber［D］.Beijing： Tsinghua University，2012. （in Chinese）

［13］ Huang D H， Huzel D K.Modern Engineering For Design of Liquid-Propellant Rocket Engines［M］.New York： American Institute of Aeronautics＆ Astronautics，1992：81-108.

［14］ Jason C W， William E A， Philip A H， et al.Supercritical flows in high aspect ratio cooling channels［R］.AIAA-2005-4302，2000.

［15］ Chaboche J， Jung O.Application of a kinematic hardening viscoplasticity model with thresholds to the residual stress relaxation［J］.International Journal of Plasticity， 1997， 13（10）：785-807.

［16］ Shafiqul B.Anatomy of coupled constitutive models for ratcheting simulation ［ J］.International Journal of Plasticity，2000， 16（3-4）：381-409.

［17］ 朱寧昌.液體火箭發(fā)動機設(shè)計［M］.北京：中國宇航出版社，1994：25-46.Zhu N C.Design of Liquid Rocket Engine［M］.Beijing：China Space Navigation Press， 1994：25-46.（in Chinese）