羽中豪，金 平，蔡國飆

（北京航空航天大學宇航學院，北京100191）

1 引言

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

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

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

2 物理模型和計算方法

2.1 物理模型

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

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

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

2.2 計算模型

2.2.1 傳熱分析部分

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

2.2.2 熱-結構耦合分析部分

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

其中，W為熱源函數；λ為內壁導熱系數；ΔT為溫度梯度。燃氣側與冷卻劑側的邊界條件為第三類邊界條件［14］，即式（2）：

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

以溫度場為基礎進行結構分析，在準靜態(tài)過程中結構場的平衡方程為式（4）［11］：

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

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

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

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

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

其中，σ為等效應力，由應力張量確定；C1、C2和C3為材料參數。

3 計算結果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由殘余應變確定的棘輪應變如式（11）：

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

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

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

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

3.3 推力室設計參數對推力室棘輪效應的影響

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

3.3.1 室壓對推力室棘輪應變的影響

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

表3為熱試階段不同室壓推力室的熱載荷情況，其中Q為熱流密度，Twg為推力室內壁溫度?？梢钥吹?，室壓更大的推力室熱載荷更重，這是因為，一方面，推力不變時提升設計室壓會使推力室整體尺寸減小，不利于再生冷卻；另一方面，更高的室壓會使燃料更充分地燃燒，推力室內的燃氣溫度也就隨之升高。

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

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

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

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

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

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

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

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

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

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

3.3.3 混合比對推力室棘輪應變的影響

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

圖13 不同混合比推力室循環(huán)工作時的棘輪應變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ā)動機燃氣燃燒更充分，推力室內溫度也更高。同時，采用再生冷卻的發(fā)動機往往以燃料作為冷卻劑，較高混合比的發(fā)動機冷卻劑的流量較小，冷卻熱流較小，不利于推力室的冷卻。所以從圖13中可以看到，混合比較小的發(fā)動機推力室的棘輪效應較輕。

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

4 結論

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

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

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

3）推力室的棘輪應變總是隨著循環(huán)工作次數的增加而趨于穩(wěn)定的；對于高室壓、大推力或是高混合比的發(fā)動機，推力室的棘輪應變隨循環(huán)工作次數的增加會有明顯的下降趨勢；對于低室壓、小推力或是低混合比的發(fā)動機，推力室的棘輪應變隨循環(huá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］ 楊進慧，陳濤，金平，等.可重復使用火箭發(fā)動機再生冷卻槽失效分析［J］.北京航空航天大學學報，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)加載三維熱結構分析［J］. 推進技術，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ā)動機推力室內壁壽命預估［J］. 航空動力學報， 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，Florida， 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］.推進技術，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ā)動機推力室結構與冷卻設計［D］.北京：清華大學，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ā)動機設計［M］.北京：中國宇航出版社，1994：25-46.Zhu N C.Design of Liquid Rocket Engine［M］.Beijing：China Space Navigation Press， 1994：25-46.（in Chinese）