宋冠華，周仕明，李道奎

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

0 引言

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

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

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

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

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

1 仿真方法介紹

1.1 弧長法

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

(1)

WEMPNER和RIKS提出的弧長法，經(jīng)過CRISFIELD，RAMN，F(xiàn)ORDE等人的修正和完善，可用于非線性屈曲極值點附近的分析，通過追蹤整個翻轉(zhuǎn)過程中實際載荷、位移關(guān)系，獲得全部狀態(tài)信息[12]。當(dāng)載荷參數(shù)為λk時，位移xk已知；當(dāng)載荷參數(shù)增量為Δλk時，相應(yīng)的位移增量為Δxk，則

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

(2)

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

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

(3)

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

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

1.2 顯式分析方法

(4)

式中：P為所施加的外力；I為單元內(nèi)力。

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

(5)

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

(6)

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

(7)

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

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

2 有限元仿真分析

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

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

圖2 錐柱形膜片結(jié)構(gòu)Fig.2 Structure of cone-bar diaphragm

膜片幾何參數(shù)如下：上部圓弧段半徑R1為200 mm，錐線段長度L為20 mm，錐角α為78°，預(yù)翻邊半徑r為3 mm。膜片為變厚度分布，預(yù)翻邊厚度為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。氣腔內(nèi)的壓強根據(jù)氣腔內(nèi)氣體的量與氣腔體積通過理想氣體方程計算，采用顯式分析法進行仿真分析。

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

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

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

(8)

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

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

2.3 顯式動態(tài)分析法與弧長法仿真結(jié)果的對比

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

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

2.4 膜片翻轉(zhuǎn)失效的動態(tài)仿真

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

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

3 充氣速度對膜片翻轉(zhuǎn)特性的影響

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

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

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

圖8 膜片翻轉(zhuǎn)失效過程Fig.8 Failure process of diaphragm overturning

圖9 動態(tài)仿真結(jié)果與實驗結(jié)果Fig.9 Dynamic simulation result and experimental result

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

膜片頂點橫向位移可反映出膜片翻轉(zhuǎn)過程中的偏心程度[18]，膜片翻轉(zhuǎn)過程中頂點橫向位移隨頂點軸向位移的變化關(guān)系如圖11所示。4種充氣速度下膜片頂點最大橫向位移分別為59，31，12，3 mm?？梢钥闯觯撼錃馑俣仍娇欤てD(zhuǎn)過程中的偏心越??；1 000 g/s的充氣速度下，膜片翻轉(zhuǎn)過程如圖12所示。與圖5比較可發(fā)現(xiàn)，在1 000 g/s速度下，膜片翻轉(zhuǎn)過程中基本無偏心。

圖10 膜片翻轉(zhuǎn)壓差曲線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 顯式分析仿真結(jié)果(充氣速度1 000 g/s)Fig.12 Simulation result of dynamic explicit method (Mass flux: 1 000 g/s)

4 結(jié)束語

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

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

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

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

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

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