金屬膜片充壓翻轉的動態仿真方法研究

宋冠華，周仕明，李道奎

(國防科技大學 空天科學學院，湖南長沙 410073)

0 引言

推進劑貯箱作為航天器的重要組成部分，具有管理控制推進劑的功能[1]。金屬隔膜貯箱因具有工作原理簡單、可靠性高、排空效率高,以及不受干擾加速度影響等特點，而被廣泛應用于需要頻繁調姿、經常變軌的航天器。此外，金屬隔膜貯箱具有強度高、耐腐蝕性好、變形過程可控、抗滲透性好等特點，可解決液體燃料晃動、精確控制和長期貯存的難題，顯著提高航天器的壽命和機動性，具有廣闊的應用前景[2]。

金屬隔膜貯箱工作時，增壓氣體擠壓膜片，將推進劑輸送到管路，其間膜片經歷了從上半球翻轉到下半球的彈塑性大變形過程，工作原理如圖1所示。金屬膜片在翻轉過程中易產生褶皺、偏心、破裂等缺陷，進而導致翻轉失效[3]。金屬膜片的翻轉特性直接影響到整個航天器能否正常運行。因此，對其翻轉規律及失效原因進行系統的研究是非常必要的。對于金屬膜片翻轉過程的研究，試驗費用過于昂貴，尚未有好的解決方法，因此利用有限元法進行數值模擬是一種有效途徑。

圖1 金屬隔膜貯箱工作原理Fig.1 Operational principle of metal diaphragm tank

國內外學者對膜片翻轉特性做了大量的研究工作。RADTKE等[4-5]利用翻轉實驗和數值模擬等手段，對金屬膜片的結構設計，選材和制造工藝進行了研究。強洪夫等[6-7]采用數值手段，研究了徑厚比和典型結構參數對錐柱形膜片翻轉特性的影響。然而，現階段對膜片翻轉過程的仿真分析大多是在膜片上施加恒定的壓力，然后采用弧長法進行計算[8-10]。這種方法的缺點是計算效率低，經常出現計算結果不收斂的問題。此外，該方法還將膜片翻轉過程作為準靜態過程處理。而實際的膜片翻轉是一個動態過程，氣瓶按一定流速向氣腔內充氣，當氣腔內壓強達到膜片翻轉臨界載荷時，膜片開始翻轉；此后，氣腔體積增大，氣腔內氣體的量也不斷增加，氣腔內壓強隨時間變化。

本文采用ABAQUS有限元軟件建立了金屬膜片充壓翻轉的流固耦合有限元模型，采用顯式動態分析法進行計算。氣腔中氣體的量、氣腔體積和壓強是隨時間變化的量。該方法可更好地模擬膜片真實翻轉的整個過程，解決了現有方法存在的計算收斂性差、計算效率低等問題，為金屬膜片的設計提供了一種重要的仿真手段。

1 仿真方法介紹

1.1 弧長法

膜片翻轉過程屬于薄殼結構產生大變形的力學問題，涉及材料非線性、幾何非線性的耦合[11]。傳統的非線性平衡方程的求解方法包括Newton-Raphson法、修正的Newton-Raphson法、Quasi-Newton法等。在求解非線性屈曲問題時，當載荷增加至臨界載荷，由于結構的切線剛度K逼近于零，求解會非常困難，因此很難得到收斂解，求解式為

(1)

WEMPNER和RIKS提出的弧長法，經過CRISFIELD，RAMN，FORDE等人的修正和完善，可用于非線性屈曲極值點附近的分析，通過追蹤整個翻轉過程中實際載荷、位移關系，獲得全部狀態信息[12]。當載荷參數為λk時，位移xk已知；當載荷參數增量為Δλk時，相應的位移增量為Δxk，則

ψ(xk+Δxk,λk+Δλk)=0

(2)

式中：Δxk和Δλk均是未知的，需要補充一個附加條件，即

f(Δxk,Δλk)=Δlk

(3)

式中：Δlk為增量弧長。聯立求解式(2)，(3)可求得既滿足平衡條件，又滿足輔助方程的Δxk和Δλk。

弧長法是目前膜片翻轉仿真分析中使用最多的方法，仿真結果與實驗結果大致吻合[13-15]，但弧長法計算效率很低，即使構型合理的膜片有時也會出現計算不收斂的問題，且該方法還將膜片翻轉過程近似作為準靜態過程處理。因此，尋找一種計算效率高、收斂性好、能模擬膜片真實翻轉過程的仿真分析方法是非常必要的。

1.2 顯式分析方法

(4)

式中：P為所施加的外力；I為單元內力。

在當前增量步開始時(t時刻)，計算加速度為

(5)

采用中心差分法對加速度在時間上進行積分，用速度變化值加上前一個增量步中點的速度來確定當前增量步中點的速度，有

(6)

速度對時間的積分加上在增量步開始時的位移可確定增量步結束時的位移，有

(7)

在增量步開始時提供滿足動力學平衡條件的加速度，然后在時間上“顯式”前推速度和位移。為了使該方法的計算結果精確，時間增量必須足夠小，這樣在增量步中加速度幾乎為常數。由于時間增量步很小，因此一個分析需要成千上萬個增量步。但不必同時求解聯立方程組，所以每一個增量步的計算成本很低，總體計算效率相對于弧長法有很大程度的提高。

顯式分析法最顯著的特點是不需要弧長法中所需的整體切線剛度矩陣。由于其為“顯式”前推模型的狀態，因此不需要迭代和收斂準則，不存在收斂性問題。另外，通過建立膜片流體腔流固耦合有限元模型，結合顯式分析法，可真實模擬膜片充壓翻轉的過程。

2 有限元仿真分析

2.1 金屬膜片充壓翻轉流固耦合有限元模型的建立

本文研究的膜片構型為錐柱形，包括下部預翻邊(I部分)、錐線段(II部分)和上部圓弧段(III部分)，如圖2所示。其中，錐線段與預翻邊圓弧和頂部圓弧均為相切幾何關系。

圖2 錐柱形膜片結構Fig.2 Structure of cone-bar diaphragm

膜片幾何參數如下：上部圓弧段半徑R1為200 mm，錐線段長度L為20 mm，錐角α為78°，預翻邊半徑r為3 mm。膜片為變厚度分布，預翻邊厚度為1.0 mm，錐線段最低點到膜片頂點，膜片厚度從1.0 mm到1.3 mm隨高度線性變化。

整體有限元模型包括膜片和外層剛體，外層剛體與膜片圍成氣腔，如圖3所示。圖3(c)為整體有限元模型的1/2。膜片采用S4R殼單元，外層剛體采用R3D4單元，膜片材料為純鈦，此處看作理想彈塑形材料，彈性模量為113 GPa，泊松比為0.32，屈服極限為250 MPa，密度為4 500 kg/m3。氣腔內的壓強根據氣腔內氣體的量與氣腔體積通過理想氣體方程計算，采用顯式分析法進行仿真分析。

圖3 流固耦合有限元模型Fig.3 Fluid-solid coupling finite element model

2.2 流固耦合有限元模型的計算原理

氣瓶以一定速度向氣腔內充氣，氣腔壓力逐漸增大。在氣腔壓力達到膜片翻轉臨界載荷之前，膜片主要發生彈性變形，膜片位移很小，氣腔體積幾乎不變，氣腔壓力隨時間線性增加。當氣腔壓力達到膜片翻轉臨界載荷時，膜片開始翻轉，此后氣腔體積隨著膜片翻轉逐漸增大。此時，氣腔中的壓強為

(8)

式中：P(t)為氣腔內氣體壓強；n(t)為氣腔內氣體的物質的量，n(t)=k·t，k為充壓氣體流速；R為理想氣體常量；T為溫度；V(t)為氣腔體積，可通過膜片翻轉過程中的節點坐標計算確定[17]。

在有限元分析中，顯式動態分析采用時間增量步的方式。首先計算t時刻氣腔內氣體的量，并根據t時刻膜片和外層剛體的節點信息計算膜片和外層剛體圍成氣腔的體積，然后根據理想氣體狀態方程計算氣腔內氣體壓強，給定時間增量Δt，最后計算膜片在氣腔壓力作用下在Δt時間增量內發生的位移，得到t+Δt時刻膜片節點信息，開始下個時間增量步的計算。圖4為流固耦合有限元模型計算流程圖。

2.3 顯式動態分析法與弧長法仿真結果的對比

分別采用顯式動態分析法和弧長法對圖3中膜片進行翻轉仿真分析，顯式動態分析法采用流固耦合有限元模型，而弧長法采用單獨的膜片有限元模型，通過在膜片上施加恒定壓力載荷完成計算。將顯式動態法與弧長法的分析結果進行對比，兩種仿真分析方法計算得到的膜片翻轉過程如圖5，6所示。其中，U2表示膜片頂點軸向位移，顯式動態分析中，充氣速度為100 g/s。可以看出，膜片從預翻邊處開始翻轉，整個過程膜片未出現褶皺、凹坑，但都出現一定程度的偏心現象，兩種仿真分析方法都能成功模擬膜片的整個翻轉過程。

整個翻轉過程中膜片翻轉壓差隨膜片頂點軸向位移的變化關系如圖7所示。兩種方法計算得到的膜片翻轉壓差曲線基本重合，驗證了顯式動態分析仿真結果的合理性。

2.4 膜片翻轉失效的動態仿真

建立與文獻[3]中完全相同的膜片有限元模型和流體腔，膜片幾何參數如下：上部圓弧段半徑R1為384 mm，錐線段長度L為144 mm，錐角α為90°，預翻邊半徑r為14 mm。膜片為變厚度分布，預翻邊厚度為2.0 mm，錐線段到膜片頂點膜片厚度從2.0 mm到3.5 mm隨高度線性變化。膜片材料為鋁，彈性模量為69 GPa，屈服強度為30 MPa，抗拉強度為70 MPa，延伸率為39%。采用顯式動態分析法進行翻轉仿真分析，得到膜片變形過程如圖8所示。可以看出：0.3 s內，膜片正常向下翻轉一段距離；0.4 s時，膜片在最下端出現屈曲現象；0.7 s時，膜片上出現大面積凹陷，翻轉失效。

通過與文獻[3]中實驗結果進行對比，如圖9所示，可以看出，顯式動態仿真結果與實驗結果基本吻合，膜片上均出現6個凹坑，屈曲模態相同。二者主要的區別是膜片變形的幅度不同，這可能是因為膜片在制造過程中存在幾何缺陷，并且膜片在屈曲之后存在接觸、擠壓、摩擦等現象，而動態仿真中忽略了這些因素。由此得出，采用顯式動態分析法能仿真出膜片翻轉失效的整個過程，并且吻合度很高。

3 充氣速度對膜片翻轉特性的影響

圖6 弧長法仿真結果Fig.6 Simulation result of arc length method

圖7 膜片翻轉壓差曲線Fig.7 Relationship between differential pressure and top node axial displacement

膜片實際翻轉過程與氣瓶充氣速度有關，為研究充氣速度對膜片翻轉特性的影響，分別以100， 200， 500， 1 000 g/s這4種充氣速度，用顯式動態分析法對圖3中膜片進行翻轉仿真分析。4種充氣速度下膜片均成功完成翻轉，未出現褶皺、凹坑。

圖8 膜片翻轉失效過程Fig.8 Failure process of diaphragm overturning

圖9 動態仿真結果與實驗結果Fig.9 Dynamic simulation result and experimental result

4種充氣速度下，膜片翻轉壓差隨膜片頂點軸向位移的變化關系如圖10所示。可以看出：4種充氣速度下膜片翻轉壓差曲線與弧長法計算得到的壓差曲線基本吻合；唯一不同的是，顯式動態分析法得到的壓差曲線在中間翻轉過程有一定程度的上下波動，并且充氣速度越快，翻轉壓差波動的幅值越大，這一特點與膜片實際翻轉過程相吻合。

膜片頂點橫向位移可反映出膜片翻轉過程中的偏心程度[18]，膜片翻轉過程中頂點橫向位移隨頂點軸向位移的變化關系如圖11所示。4種充氣速度下膜片頂點最大橫向位移分別為59，31，12，3 mm。可以看出：充氣速度越快，膜片翻轉過程中的偏心越小；1 000 g/s的充氣速度下，膜片翻轉過程如圖12所示。與圖5比較可發現，在1 000 g/s速度下，膜片翻轉過程中基本無偏心。

圖10 膜片翻轉壓差曲線Fig.10 Relationship between differential pressure and top node axial displacement

圖11 膜片偏心程度Fig.11 Relationship between top node transversal displacement and top node axial displacement

圖12 顯式分析仿真結果(充氣速度1 000 g/s)Fig.12 Simulation result of dynamic explicit method (Mass flux: 1 000 g/s)

4 結束語

本文建立膜片充壓翻轉的膜片流體腔流固耦合有限元模型進行顯式動態仿真，通過與弧長法進行分析結果、實驗結果的對比，驗證了顯式動態分析仿真結果的正確性；通過改變充氣速度，研究了充氣速度對膜片翻轉過程的影響。主要結論如下：

1) 與弧長法相比，顯式分析法計算效率更高，不存在收斂性問題。

2) 通過建立流固耦合有限元模型，采用顯式動態分析法，將膜片翻轉過程作為動態過程處理，氣體按一定流速向氣腔內充氣，膜片上所受壓力依賴于氣腔體積和氣腔內氣體的量，這是隨時間變化的量，與膜片實際的充壓翻轉過程更接近。

3) 膜片實際翻轉過程受氣瓶充氣速度影響，通過建立流固耦合有限元模型，結合顯式動態分析法，可以研究充氣速度對膜片翻轉特性的影響，得出不同充氣速度下翻轉壓差曲線與偏心程度。

4) 充氣速度越快，膜片中間平穩翻轉過程的壓力波動越大，偏心程度越小，這一特點與膜片實際翻轉過程相吻合,而弧長法無法仿真出這一特點。

本文提出的顯式動態分析方法法也可應用于球形、頂部凹陷形等其他類型的金屬膜片翻轉仿真，以及其他基于理想氣體狀態方程的有限元仿真分析。本文提出的仿真方法，可用于構型、幾何參數、材料等對膜片翻轉效果影響的研究，解決了目前仿真方法中存在的構型合理的膜片計算結果不收斂的問題，為金屬膜片的優化設計提供了參考。