基于Ansys-Workbench EPKM點火過程流固耦合仿真分析①

劉迎賓，李 卓，江曉瑞

(內蒙古工業大學 理學院，呼和浩特 010051)

0 引言

固體火箭發動機的點火過程具有持續時間短、升壓梯度大的特點。點火過程中，藥柱內表面容易產生升壓不均勻的現象，導致藥柱結構完整性發生破壞，這屬于發動機內流場與固體推進劑相互作用的流固耦合問題。

關于固體火箭發動機流固耦合問題的研究始于1991年，美國PQM-1固體火箭助推器在地面靜態試驗時點火階段發生爆炸，CHANG等通過流固耦合仿真計算，對爆炸原因進行了分析，為其改進提供了可靠依據。近年來，國外一些學者在流固耦合計算模型及數值方法方面做了研究，如JOHNSTON等采用一種特殊的時間迭代技術將流場求解器、結構場求解器進行耦合，用于分析運載火箭的推力振蕩現象，TOLA等將內彈道流場與結構場相互耦合分析，對槽形藥柱的結構進行優化。國內關于固體火箭發動機點火過程流固耦合的研究起步較晚。2000年，國防科技大學張為華等利用N-S方程，新型的ENO差分格式進行數值計算，定性研究了點火藥噴流與主裝藥形狀的相互影響，提出固體火箭發動機裝藥設計中應該注意的問題。之后，周紅梅、曹琪等指出了對固體火箭發動機點火過程流固耦合研究的重要性，并利用MPCCI軟件作為流場與結構場數據的交換平臺，分別對含有裂紋裝藥及帶有徑向翼槽大長徑比裝藥的固體火箭發動機點火過程進行了仿真分析。黃永恒利用流固單向耦合方法對管形及星形裝藥發動機點火過程中藥柱的受力情況及影響因素進行分析。近年來，一些學者通過冷流沖擊試驗來模擬發動機藥柱在點火過程中受高壓燃氣的沖擊作用，利用流固耦合分析方法，對藥柱結構場在沖擊過程中的變化規律進行研究。

通過對前人工作的學習與總結，發現有關固體火箭發動機點火過程流固耦合的研究成果相對較少，而且由于計算機軟硬件的限制，前人大多只能做一些相對簡單的計算，卻需要大而繁瑣的工作。近年來，隨著Ansys商用軟件的發展，其Workbench可以通過調用其內部的各個模塊，進行耦合分析。本文通過Workbench調用Fluent及Mechanical APDL模塊分別用于對發動機流場及結構場的計算，并通過調用System Coupling模塊對其進行耦合分析，以EPKM型固體火箭發動機為算例，研究了其點火升壓階段燃燒室內流場及結構場的變化規律，為固體火箭發動機的設計及故障診斷提供參考。

1 仿真計算模型

1.1 幾何模型

EPKM型固體火箭發動機用于將載荷從近地點軌道推進到地球同步轉移軌道，其主裝藥頭部含有12個沿周向均勻分布的翼槽。為了節約計算資源，只對從翼槽的中間對稱面，繞發動機軸線旋轉15°的區域進行建模，模型尺寸見表1。由于在點火計算過程中，固體推進劑點燃后，需要利用Fluent的源項法，通過UDF(User-defined function)從靠近藥柱內表面的一薄層流場區域內加入高溫燃氣，故在建模過程中，需對該區域進行單獨建模，即圖1中Source area。采用四面體非結構化網格對模型進行網格劃分，如圖1所示。

表1 EPKM模型尺寸

圖1 EPKM幾何模型及其網格劃分

1.2 控制方程

燃氣流動與傳熱的連續方程，動量方程，能量方程可用共同的形式表達，使用通用變量，通用微分方程為

(1)

式中 第一項為瞬變項或時間項；第二項為對流項；第三項為擴散項；為對應于變量的擴散系數；末項為源項。

結構場單元基本平衡代數方程為

[]{}=[]{ü}+{}+{}

(2)

1.3 計算模型

湍流模型選用Spalart-Allmaras單方程模型，其輸運方程：

(3)

固體推進劑的燃速模型：

=1044×100278

(4)

式中為燃燒室壓強。

溫度為3300 K的燃氣在相應流率下所具有的焓可用下式計算：

(5)

式中為燃氣的比定壓熱容，J/(kg·K)。

固體推進劑為粘彈性材料，其體積松弛模量用1～9參數的Prony級數給定：

(6)

其中，各參數如表2所示，其瞬時彈性模量(=0)=18.87 MPa，泊松比為0.498，密度1750 kg/m。

EPKM的點火裝置為點火發動機，采用HTPB系統配方推進劑藥柱，計算中假設點火裝置產生的燃氣與主裝藥產生的燃氣具有相同物理性質，如表3所示。

表2 Prony級數參數

表3 燃氣的物理性質

1.4 流固耦合方法

文中采用的流固耦合方法包括兩種，如圖2所示。其中，圖2(a)表示流固雙向耦合，主要通過System coupling對耦合迭代的時間步長及流場與結構場對應面上壓力與變形量數據的傳遞進行控制。圖2(b)表示一種簡化的流固單向耦合方法，其計算過程是先完成流場部分的計算，再將不同時刻流場對結構場的壓力直接傳遞到結構場進行受力分析，該方法不需要在每個時間步進行數據的雙向傳遞。因此，計算速度相對于第一種會非常快，但只能計算流場流動狀態基本不受結構場變形影響的情況，且只能是某些間斷時間點處結構場的狀態。

(a)Bidirectional coupling (b)Unidirectional and simplified coupling

1.5 初始和邊界條件

仿真計算中流場區域入口位置如圖3所示，將其設置為質量流率入口，流率按圖示曲線變化，為保證主裝藥的穩定燃燒，點火發動機在整個點火過程一直處于工作狀態。流場區域出口，即噴管出口截面設置為壓強出口，大氣背壓約1000 Pa，溫度約1200 K。堵蓋打開前，將堵蓋處邊界設置為Wall，堵蓋打開后，將其改為Interface，堵蓋的打開壓強為2.5 MPa。流場模型的側面設置為Symmetry對稱邊界。初始狀態，燃燒室內整個流場速度為零，壓強和溫度為在地面對燃燒室進行密封時的壓強和溫度，即1 atm和300 K。

圖3 點火燃氣入口位置及其流量曲線

將藥柱內表面與流場外表面都設置為流固耦合面，并在流場中設置移動邊界。計算過程中，將流場的壓力數據傳遞到結構場，將結構場產生的變形量數據傳遞到流場。將藥柱的剩余面都設置為光滑面約束。

只有在高溫燃氣流經藥面時，藥面才會被加熱引燃。因此，用藥柱內表面附近流體單元的溫度和壓強來判斷藥柱是否被點燃。通過對Fluent進行編程嵌入，計算過程中獲取藥面附近流體單元的溫度和壓強并進行判斷，當其同時達到一定值時，認為該單元附近藥柱被點燃。然后，從該單元處往流場中加入3300 K的高溫燃氣。流率由燃速決定：

(7)

計算過程中忽略侵蝕燃燒對推進劑燃速的影響。由于點火過程所用時間極短，因此忽略藥柱燃面的徑向位移。不考慮藥柱由于固化降溫引起的溫度載荷，只對點火升壓過程中的內壓載荷進行計算。

1.6 網格無關性驗證

為分別對流場和結構場網格進行無關性驗證，先增加流場網格數，計算0.1 ms時發動機頭部某點處壓強，再增加結構場網格數，計算藥柱在該時刻的最大應力，如表4所示。

表4 不同網格數計算結果

隨著網格數增加，流場某點壓強變化最大約 0.75%，結構場最大應力變化最大約0.26%。網格數的增加對計算結果影響很小。最終采用流場網格數230 392，結構場網格數26 246。

2 仿真結果分析

藥柱在初始幾個毫秒內，由于點火燃氣的沖擊作用，會產生相對較大的變形，因此采用圖2(a)所示流固雙向耦合計算方法對點火過程計算5 ms，藥柱最大變形量隨時間變化曲線如圖4所示。可以發現，其最大變形量不到1 mm，對流場的影響基本可以忽略。因此，為了加快計算速度，假設在整個點火過程中藥柱的變形量很小，基本不會對流場產生影響，采用圖2(b)所示簡化的流固單向耦合方法完成后續的計算，并根據燃燒室壓強達到穩定工作壓強時藥柱的最大變形量大小對假設進行驗證。

圖4 結構場最大變形量隨時間變化曲線

2.1 燃燒室內流場分析

主裝藥點燃前，燃燒室內燃氣完全來源于點火發動機。由圖5可知，點火燃氣由點火發動機噴出后，首先不會直接流過藥柱的頭部內表面，而是先從燃燒室軸線附近流動到其尾部，經堵蓋作用，再回流到頭部翼槽區域。在翼槽內的流動情況如圖5中=26 ms、=34 ms和=57 ms所示，先流動到翼槽的前緣，然后沿著前緣逐漸向其底部和后緣流動，最后填充滿整個翼槽區域，而且發動機燃燒室內初始時刻的冷空氣會被點火燃氣壓縮到翼槽內，使該區域的溫度相對較低，藥柱的點燃延時相對較長。

(a)t=10 ms (b)t=26 ms

(c)t=34 ms (d)t=57 ms

點火燃氣剛噴入燃燒室后會形成復雜的多個渦流，其后它們會逐漸進行融合并變得穩定，見圖6。在=12 ms時，管型藥柱段內的流場已基本變得穩定，只存在一個相對較大的渦流，在=18 ms時，翼槽內流場也變得穩定，到=22 ms時，整個燃燒室流場基本達到穩定的狀態。此時，燃燒室管型藥柱段和翼槽內各存在一個相對較大且穩定的渦流，直到堵蓋打開，該流動狀態基本不會再發生變化。

(a)t=2.5 ms (b)t=12ms

(c)t=18 ms (d)t=22 ms

對比圖7中=57 ms和=57.8 ms時刻主裝藥內表面沿軸向溫度分布曲線，顯然在=57.8 ms時，在距燃燒室頭部約0.6 m處溫度突然升高，即在=57.8 ms時，主裝藥已經被點燃，而且由圖8可知，首先被點燃的位置位于翼槽到管型藥柱段的過渡區域。

(a)t=57 ms

(b)t=57.8 ms

圖8 t=57.8 ms時翼槽附近溫度云圖

主裝藥點燃后，燃面從首先被點燃處分別向頭尾兩端擴展，且在管型藥柱段的擴展速度遠大于向翼槽內的擴展速度，如圖9所示。其原因如前所述，翼槽內含有大量壓縮的初始時刻燃燒室內冷空氣，相對溫度較低。燃面在翼槽內的擴展過程為：首先，在翼槽到管型藥柱段的過渡區被點燃，燃面以非常慢的速度向翼槽內擴展，大量推進劑燃氣會隨著從燃燒室尾部回流的點火燃氣向翼槽的前緣流動，前緣被點燃；然后，燃面逐漸向翼槽底部和后緣方向擴展。

當燃燒室堵蓋面平均壓強達到2.5 MPa(約67 ms時)，堵蓋打開，燃燒室內渦流逐漸消失，如圖10所示。最終主裝藥產生的高溫燃氣直接流向噴管，由噴管噴出，不再產生渦流，燃燒室內達到新的穩定流動狀態。

(a)t=64 ms (b)t=68 ms

(c)t=90 ms (d)t=117 ms

(a)t=67.5 ms (b)t=69 ms

(c)t=72.5 ms (d)t=80 ms

2.2 藥柱結構場分析

點火初期，發動機燃燒室堵蓋打開前，由圖11中的=13 ms和=64 ms時刻藥柱變形量云圖可知，藥柱尾部的變形大于其頭部，原因是由于噴管堵蓋的作用，使得燃燒室內流場尾部的壓強始終大于其頭部，即藥柱在尾部受到的壓力始終大于其頭部，這從圖12與圖13也可以看出。

(a)t=13 ms (b)t=64 ms

(c)t=69 ms (d)t=136 ms

(a)t=13 ms (b)t=64 ms

(c)t=69 ms (d)t=136 ms

(a)t=13 ms (b)t=64 ms

(c)t=69 ms (d)t=136 ms

堵蓋打開后，由圖11中的=69 ms和=136 ms時刻云圖可知，藥柱最大變形量發生的位置逐漸從尾部向頭部移動，當流場基本達到穩定的流動狀態時，即=136 ms時刻所對應的狀態，藥柱最大變形量約5.4 mm，位置在翼槽的側面，這也證明了前面關于結構場變形很小，基本不會對流場流動狀態產生影響的假設的正確性。堵蓋打開后，使得燃燒室堵蓋附近壓強突然下降，但由于此時燃燒室內整體壓強較高，整個藥柱內表面壓強受堵蓋打開影響很小。因此，整個藥柱內表面應力分布比較均勻，且沿徑向逐漸降低，如圖12、圖13中=69 ms及=136 ms時刻云圖所示。

在=136 ms時，燃燒室內壓強基本已達到工作壓強，該時刻沿翼槽對稱面與藥柱內表面交線的應力、應變分布曲線如圖14所示，最大應力約1.16 MPa，對應最大應變約6.2%，發生在圖中點所對應的位置。對于無缺陷的固體推進劑在內壓載荷作用下，其破壞準則一般采用八面體剪切應變準則：

<()

(8)

式中為臨界值；為安全系數。

式中為該固體推進劑的泊松比,=0.498；為圖13與圖14所示的等效應變。

該推進劑單向拉伸最大延伸率≥40%，當安全系數取為5時，等效應變的最大值6.2%<8%，即<(15)。因此，藥柱的結構完整性滿足要求。

圖14 t=136 ms沿圖中特征線應力、應變曲線

3 結論

(1)點火燃氣噴入燃燒室后，形成的復雜的多個渦流會很快的進行融合，最終在藥柱的管型裝藥段及各個翼槽內各融合為一個穩定的渦流，且在噴管堵蓋打開之前，該流動狀態會一直保持穩定。因此，在發動機研制中，需要重點考慮點火燃氣剛噴入燃燒室的幾個毫秒。

(2)翼槽內燃面的擴展方向主要取決于其內部燃氣的流動方向，與燃氣流動方向一致。

(3)點火過程中，由于噴管堵蓋的作用，會使藥柱尾部所受應力大于其頭部，當噴管堵蓋打開后，整個結構場應力沿軸向分布比較均勻，沿徑向逐漸降低，且應力及應變的最大值發生在翼槽與管形藥柱段的過渡處。因此，在實際發動機藥柱設計時，應注意堵蓋附近及藥面形狀過渡處藥柱的強度。

由于計算時間及計算機硬件的限制，本文沒有對整個點火過程進行流固雙向耦合計算，有待于進一步研究。