宋冠華，周仕明，李道奎

(國(guó)防科技大學(xué) 空天科學(xué)學(xué)院，湖南長(zhǎng)沙 410073)

0 引言

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

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

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

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

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

1 仿真方法介紹

1.1 弧長(zhǎng)法

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

(1)

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

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

(2)

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

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

(3)

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

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

1.2 顯式分析方法

(4)

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

在當(dāng)前增量步開(kāi)始時(shí)(t時(shí)刻)，計(jì)算加速度為

(5)

采用中心差分法對(duì)加速度在時(shí)間上進(jìn)行積分，用速度變化值加上前一個(gè)增量步中點(diǎn)的速度來(lái)確定當(dāng)前增量步中點(diǎn)的速度，有

(6)

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

(7)

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

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

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，錐線段長(zhǎng)度L為20 mm，錐角α為78°，預(yù)翻邊半徑r為3 mm。膜片為變厚度分布，預(yù)翻邊厚度為1.0 mm，錐線段最低點(diǎn)到膜片頂點(diǎn)，膜片厚度從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)的壓強(qiáng)根據(jù)氣腔內(nèi)氣體的量與氣腔體積通過(guò)理想氣體方程計(jì)算，采用顯式分析法進(jìn)行仿真分析。

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

2.2 流固耦合有限元模型的計(jì)算原理

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

(8)

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

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

2.3 顯式動(dòng)態(tài)分析法與弧長(zhǎng)法仿真結(jié)果的對(duì)比

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

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

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

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

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

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

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

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

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

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

圖9 動(dòng)態(tài)仿真結(jié)果與實(shí)驗(yàn)結(jié)果Fig.9 Dynamic simulation result and experimental result

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

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

圖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é)束語(yǔ)

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

1) 與弧長(zhǎng)法相比，顯式分析法計(jì)算效率更高，不存在收斂性問(wèn)題。

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

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

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

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