楊鵬飛， 陳 剛， 薛 杰， 王 珺

(1.西安交通大學(xué) 航天航空學(xué)院，西安 710049； 2.氣動(dòng)噪聲與控制重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000；3.西安航天動(dòng)力研究所，西安 710100)

液體晃動(dòng)問題是指液體自由表面由于受到外加激勵(lì)而產(chǎn)生的運(yùn)動(dòng)，這種問題廣泛存在于工程實(shí)際應(yīng)用的各個(gè)領(lǐng)域，如航空航天、石油化工等。就航天器貯箱來說，其受力及振動(dòng)環(huán)境復(fù)雜，微小的外激勵(lì)可能會(huì)造成內(nèi)部液體劇烈晃動(dòng)，引起附加慣性質(zhì)量分布和流-固耦合剛度改變，嚴(yán)重時(shí)會(huì)影響結(jié)構(gòu)振動(dòng)模態(tài)。特別是當(dāng)儲(chǔ)液容器結(jié)構(gòu)很薄時(shí)，內(nèi)部液體晃動(dòng)可能與貯箱彈性模態(tài)耦合，產(chǎn)生不可預(yù)料和控制的振動(dòng)，這對航天器的安全性和穩(wěn)定性產(chǎn)生了極大影響。

目前對液體晃動(dòng)問題的研究也越來越受到研究者們重視。近幾十年來，國內(nèi)外學(xué)者對貯箱液體晃動(dòng)問題已展開了深入的研究。Biswal等[1]在不同的圓環(huán)形阻尼板的尺寸、位置及數(shù)量情況下，通過有限元方法研究了液體晃動(dòng)特性。Cho等[2]同樣在有阻尼板作用下采用了有限元方法研究了二維矩形容器的液體晃動(dòng)特性，探討了不同阻尼板的尺寸及安裝位置對晃動(dòng)特性的影響。郝亞娟等[3]采用單一拉格朗日法,對理想流體貯箱中矩形隔層板的耦合振動(dòng)問題進(jìn)行了研究,探討了上下層流體的高度、隔層板的幾何尺寸、流體密度、流體晃動(dòng)頻率等因素對隔層板動(dòng)力特性的影響。王瓊瑤等[4]針對部分充液罐車轉(zhuǎn)向或制動(dòng)過程中罐體內(nèi)液體的晃動(dòng)問題，研究了防波板的幾何參數(shù)(開孔的大小、形狀及位置等)對液體瞬態(tài)晃動(dòng)時(shí)的載荷轉(zhuǎn)移及晃動(dòng)力的影響。Chen等[5]利用改進(jìn) 的SPH法模擬了二維水箱液面晃動(dòng)問題，得到了液體自由液面形狀和壁面處最大壓力，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對比驗(yàn)證。程選生等[6]基于勢流理論在微幅晃動(dòng)范圍內(nèi)研究矩形貯液結(jié)構(gòu)的動(dòng)力響應(yīng)，結(jié)果表明在微幅晃動(dòng)范圍內(nèi)，液體最大晃動(dòng)波高和地震加速度幅值成線性關(guān)系。張展博等[7]研究了低頻大幅晃動(dòng)條件下水平多孔擋板不同浸沒深度對液面晃動(dòng)的影響，發(fā)現(xiàn)了波峰和波谷變化的規(guī)律以及在特殊情況下?lián)醢鍖?dǎo)致波峰高度增加的現(xiàn)象。

對于考慮流固耦合效應(yīng)的液體晃動(dòng)問題，Strand等[8]研究了具有柔性壁面的二維矩形水箱，分析了液體晃動(dòng)與柔性壁面間的耦合關(guān)系。劉小民等[9]針對鐵路罐車剎車過程中的液體晃動(dòng)進(jìn)行了研究，分析了不同材料下防波板受液體晃動(dòng)沖擊時(shí)的受力情況 ，同時(shí)揭示了充液比對罐體壁面及防波板受力的影響。得出了罐體端面受力隨充液比的增加而增大。吳作偉等[10]以高速動(dòng)車組車下水箱為研究對象，采用雙向流固耦合方法，對車載水箱箱體內(nèi)液體沖擊晃動(dòng)問題進(jìn)行了建模與數(shù)值分析，結(jié)果表明:箱體所受最大應(yīng)力基本與充液比成正比 ,且最大應(yīng)力位置基本相同。姚燦等[11]針對薄壁貯箱的流固耦合問題，在理論分析基礎(chǔ)上，采用有限元數(shù)值方法進(jìn)行貯箱在充液不加壓和充液加壓兩種工況下的模態(tài)分析，得到流固耦合下貯箱的含阻尼非線性固有頻率。Yang等[12]采用雙向流固耦合方法，研究了三維儲(chǔ)罐的液面分布特征和容器的結(jié)構(gòu)響應(yīng)，結(jié)果表明了儲(chǔ)罐結(jié)構(gòu)的響應(yīng)振幅隨著地震的增加而增大，并主要受液體晃動(dòng)的影響。

上述研究主要探討了對擋板的幾何參數(shù)與液體晃動(dòng)的相互影響關(guān)系，大多沒有采用雙向流固耦合方法，僅有的少數(shù)采用雙向流固耦合方法的研究也沒有分析結(jié)構(gòu)柔性對液體晃動(dòng)和結(jié)構(gòu)應(yīng)力分布的影響。本文采用雙向流固耦合方法，獲取液貯箱液體晃動(dòng)波高以及貯箱結(jié)構(gòu)應(yīng)力隨時(shí)間的變化歷程，研究了結(jié)構(gòu)柔性對液體晃動(dòng)和結(jié)構(gòu)應(yīng)力分布的瞬態(tài)響應(yīng)的影響。分析了貯箱的材料、厚度以及箱內(nèi)液體的充液比對彈性貯箱液體晃動(dòng)及結(jié)構(gòu)受力特性的影響規(guī)律和流固耦合效應(yīng)，從而為利用結(jié)構(gòu)柔性進(jìn)行彈性貯箱液體晃動(dòng)抑制提供設(shè)計(jì)依據(jù)和方法支撐。

1 數(shù)值方法

1.1 流動(dòng)控制方程

假設(shè)箱內(nèi)液體無旋、不可壓縮，流體模型采用層流模型，忽略溫度的影響，則貯箱內(nèi)液體流動(dòng)滿足三個(gè)基本方程：

1.1.1 質(zhì)量守恒方程

單位時(shí)間內(nèi)流體微元中質(zhì)量的增量，等于同一時(shí)間間隔內(nèi)流入該微元的凈質(zhì)量。將質(zhì)量守恒方程寫成微分形式，表示如下：

(1)

式中:散度divu=?u1/?x+?u2/?y+?u3/?z,t、u、ρ分別為時(shí)間、速度矢量、流體密度，u1、u2和u3是速度矢量u在各方向上的分量。

1.1.2 動(dòng)量守恒方程

在單位時(shí)間內(nèi)流場域流體微元體中動(dòng)量增加的量，必須等于同一單位時(shí)間內(nèi)出入該微元體流體的凈動(dòng)量和凈沖量的和，其中邊界表面上的表面力和整個(gè)體積區(qū)域上的體積力會(huì)對微元體產(chǎn)生沖量。寫成微分形式，并利用連續(xù)性方程化簡得：

(2)

1.1.3 能量守恒方程

單位時(shí)間內(nèi)流體微元中能量的增量，等于同一時(shí)間間隔內(nèi)流入該微元的能量與外力做功的和。其微分形式為：

(3)

e表示單位質(zhì)量流體的內(nèi)能，體積力F=(Fx,Fy,Fz)。

1.2 結(jié)構(gòu)動(dòng)力學(xué)方程

由經(jīng)典力學(xué)理論，動(dòng)力學(xué)平衡方程表示如下：

(4)

1.3 流固耦合基本方程

綜合固體與流體控制方程，流固耦合問題的有限元控制方程可以整理為下式:

(5)

式中，

式中：q表示來自流體的變量，fs表示流體結(jié)構(gòu)的相互作用。

1.4 雙向流固耦合策略

在很多耦合問題中，流體的作用力影響結(jié)構(gòu)的變形，同時(shí)結(jié)構(gòu)的位移又影響流場的形態(tài)，我們把這種類型的分析叫做“雙向耦合”。雙向流固耦合計(jì)算分為流體域和固體域兩部分，流體和固體計(jì)算通過耦合界面?zhèn)鬟f數(shù)據(jù)，不考慮溫度的影響，流固耦合界面滿足運(yùn)動(dòng)學(xué)條件和動(dòng)力學(xué)條件，即：

(6)

式中：df和ds分別表示流體和結(jié)構(gòu)的位移，τf和τs分別為流體和結(jié)構(gòu)的應(yīng)力。

對于上述方程，當(dāng)給定初始條件及邊界條件后，可以通過建立耦合方程的通用形式進(jìn)行統(tǒng)一求解，目前，主要有兩種求解方法：直接耦合解法和分離法。本文采用直接耦合解法，將流體和固體的控制方程耦合至同一方程矩陣求解，方程表述如下：

(7)

求解過程概括如下：首先，將總計(jì)算時(shí)間劃分為多個(gè)時(shí)間步，在每個(gè)時(shí)間步內(nèi)，求解耦合系統(tǒng)的線性方程組并更新得到的解。計(jì)算應(yīng)力和位移殘量并和給定的迭代容差比較。如果解還不收斂，且沒有到達(dá)流固耦合迭代的最大次數(shù)，則回到第1步繼續(xù)下一個(gè)迭代。隨后，該時(shí)間步的結(jié)果將作為下一時(shí)間步的初值進(jìn)行計(jì)算，繼而完成所有時(shí)間步的求解，如圖 1 所示。

圖1 雙向流固耦合方法

2 計(jì)算模型及條件

2.1 計(jì)算模型

本文研究對象的結(jié)構(gòu)尺寸如圖2所示，貯箱的圓柱部分的長度為142 mm，箱體球形端半徑為213 mm ，隔板位于圓柱端與球形端的交接處，圓柱端外部設(shè)有法蘭結(jié)構(gòu)，具體位置和尺寸如圖所示。箱體厚度2 mm,法蘭結(jié)構(gòu)厚度5 mm。

圖2 計(jì)算模型及結(jié)構(gòu)尺寸示意圖

對箱體和流體分開建模，建模時(shí)保持殼單元和流體單元接觸面處節(jié)點(diǎn)一一對應(yīng)。在幾何空間位置上，流場和結(jié)構(gòu)場的流固耦合邊界應(yīng)該重合，由于建模時(shí)需要考慮薄壁特點(diǎn)，選用殼單元來建立罐壁模型；設(shè)罐體外一點(diǎn)與罐體進(jìn)行剛性連接，其方向?yàn)閤軸正方向，罐體重力加載到x軸負(fù)方向；沿罐體軸向?yàn)閥方向。流體域的網(wǎng)格數(shù)為45 080，固體域的網(wǎng)格數(shù)為5 958。貯箱整體結(jié)構(gòu)(左)和流體(右)的有限元模型如圖3。

圖3 貯箱結(jié)構(gòu)及流體網(wǎng)格

2.2 載荷加載

采用輸入加速度時(shí)程曲線的方法對貯箱結(jié)構(gòu)施加振動(dòng)，該方法加載簡單，且易于實(shí)現(xiàn)。選取的振動(dòng)時(shí)程曲線如圖4。定義時(shí)間步長為0.002 s,共計(jì)算500步，計(jì)算總時(shí)間1 s，加速度加載到y(tǒng)方向。

圖4 加速度時(shí)程曲線

3 結(jié)果分析

3.1 不同貯箱材料對流固耦合響應(yīng)的影響

本部分在充液比為0.5,箱內(nèi)液體為酒精，箱體厚度為2 mm的情況下，采用雙向流固耦合方法研究了不同材料的貯箱模型對箱內(nèi)液體微幅晃動(dòng)及結(jié)構(gòu)應(yīng)力變化的影響。箱體材料及其相應(yīng)的材料屬性，如表1所示。

表1 貯箱材料及其屬性

圖5為不同材料的貯箱模型液體晃動(dòng)波高時(shí)程曲線。從圖中可以看出，貯箱內(nèi)液體晃動(dòng)呈周期性變化，隨著時(shí)間的增加，液體晃動(dòng)愈來愈劇烈，在0.688 s時(shí)達(dá)到峰值。

圖5 液體晃動(dòng)波高時(shí)程曲線

圖6為不同材料的貯箱在液體晃動(dòng)波高峰值時(shí)刻云圖，在0.688 s時(shí)刻，結(jié)構(gòu)鋼材料貯箱模型晃動(dòng)波高最大，為18.23 mm；鈦合金材料貯箱模型晃動(dòng)波高次之，為14.28 mm;鋁合金材料貯箱模型晃動(dòng)波高最小，為10.85 mm。

圖6 液體晃動(dòng)波高峰值時(shí)刻云圖

可以發(fā)現(xiàn)，隨著材料楊氏彈性模量和密度的減小，貯箱內(nèi)液體晃動(dòng)波高峰值減小，在外加激勵(lì)不變的情況下，晃動(dòng)幅度逐漸趨于平緩。

圖7為不同材料的貯箱模型等效應(yīng)力隨時(shí)間變化曲線，從圖中可以看出，箱體等效應(yīng)力呈不規(guī)則的周期振蕩，隨著材料楊氏模量和密度的增加，箱體所受的等效應(yīng)力增大，且振蕩幅度增加。

圖7 貯箱等效應(yīng)力時(shí)程曲線

圖8為不同材料的貯箱模型應(yīng)力峰值時(shí)刻云圖，從圖中可以看出，等效應(yīng)力最大值均位于箱體靠近法蘭位置，在0.004 s時(shí)刻，鋁合金材料貯箱最大等效應(yīng)力為523 599 Pa,在0.216 s時(shí)，鈦合金材料貯箱和結(jié)構(gòu)鋼材料貯箱等效應(yīng)力達(dá)到最大值，分別為584 793 Pa和698 608 Pa?？梢园l(fā)現(xiàn)，隨著材料楊氏模量和密度的增加，貯箱模型等效應(yīng)力增加。

圖8 貯箱等效應(yīng)力峰值時(shí)刻云圖

3.2 貯箱厚度對流固耦合響應(yīng)的影響

本部分采用雙向流固耦合方法對2 mm、3 mm和4 mm厚度的貯箱液體晃動(dòng)及其箱體結(jié)構(gòu)所受應(yīng)力的影響進(jìn)行了對比研究。貯箱箱體材料為鈦合金，箱內(nèi)充酒精。

圖9為不同厚度的貯箱模型液體晃動(dòng)波高時(shí)程曲線。如圖所示，貯箱內(nèi)液體晃動(dòng)呈周期性變化，隨著時(shí)間的增加，液體晃動(dòng)愈來愈劇烈，在達(dá)到峰值后晃動(dòng)幅度逐漸減小。在0.688 s時(shí)刻，3種貯箱模型晃動(dòng)波高均達(dá)到峰值，2 mm貯箱模型晃動(dòng)波高最小，為14.28 mm; 3 mm貯箱模型晃動(dòng)波高較大，為17.12 mm;在0.686 s時(shí)，4 mm貯箱模型晃動(dòng)波高最大，為19.03 mm?？梢园l(fā)現(xiàn)，隨著貯箱厚度的增加，液體晃動(dòng)幅度增大。

圖9 液體晃動(dòng)波高時(shí)程曲線

圖10為不同材料的貯箱模型應(yīng)力峰值時(shí)刻云圖。從圖中可以看出，在0.216 s時(shí)刻，充2 mm厚度貯箱最大等效應(yīng)力為584 793 Pa；在0.268 s時(shí)，3 mm厚度貯箱最大等效應(yīng)力為418 996 Pa，4 mm厚度貯箱最大等效應(yīng)力為342 455 Pa?？梢园l(fā)現(xiàn)，隨著貯箱厚度的增加，貯箱模型等效應(yīng)力減小。

圖10 貯箱等效應(yīng)力時(shí)程曲線

3.3 不同充液比對貯箱流固耦合響應(yīng)的影響

本部分采用雙向流固耦合方法研究不同充液比對箱內(nèi)液體晃動(dòng)和箱體的影響。仍然選取鈦合金作為貯箱箱體材料，箱體厚度2 mm,貯箱內(nèi)充入液體為酒精。

圖11為不同充液比的貯箱模型液體晃動(dòng)波高時(shí)程曲線。從圖中可以看出，貯箱內(nèi)液體晃動(dòng)呈周期性變化，隨著時(shí)間的增加，液體晃動(dòng)愈來愈劇烈，在達(dá)到峰值后晃動(dòng)幅度逐漸減??；在0.716 s時(shí)刻，充液比為0.2貯箱峰值波高最大，為18.91 mm；在0.688 s時(shí)，充液比為0.5貯箱峰值波高次之，為14.28 mm;在0.6 s時(shí)，充液比為0.8貯箱峰值波高最小，為9.743 mm。隨著貯箱充液比的增加，液體晃動(dòng)幅度減小，且晃動(dòng)峰值時(shí)刻提前。

圖11 液體晃動(dòng)波高時(shí)程曲線

圖12為不同材料的貯箱模型應(yīng)力峰值時(shí)刻云圖。從圖中可以看出，在0.216 s時(shí)刻，充液比為0.2貯箱最大等效應(yīng)力為352 640 Pa；在0.216 s時(shí)，充液比為0.5貯箱最大等效應(yīng)力為584 793 Pa；在0.004 s時(shí)，充液比為0.8貯箱最大等效應(yīng)力為806 876 Pa?？梢园l(fā)現(xiàn)，隨著充液比的增加，貯箱模型等效應(yīng)力增大。

4 結(jié) 論

本文針對薄壁貯箱的流固耦合問題，采用雙向流固耦合方法對考慮流固耦合效應(yīng)的充液貯箱進(jìn)行了瞬態(tài)響應(yīng)分析。通過改變貯箱的材料、充液比以及貯箱厚度，研究了貯箱結(jié)構(gòu)柔性對液體晃動(dòng)及結(jié)構(gòu)受力特性的影響。得到以下結(jié)論：

(1)減小貯箱材料的楊氏模量和密度，箱內(nèi)液體的晃動(dòng)幅度減小，同時(shí)貯箱結(jié)構(gòu)等效應(yīng)力減小?？梢园l(fā)現(xiàn)，適當(dāng)增加模型的柔性，可以有效抑制箱內(nèi)液體晃動(dòng)，同時(shí)減小液體晃動(dòng)對箱體造成的強(qiáng)度破壞。

(2)增加貯箱結(jié)構(gòu)的厚度，箱內(nèi)液體晃動(dòng)幅度增大，但貯箱結(jié)構(gòu)受到的應(yīng)力減小。

(3)通過改變貯箱內(nèi)液體的充液比，可以發(fā)現(xiàn)，對本文所研究的貯箱模型，隨著充液比的增加，箱內(nèi)液體的晃動(dòng)幅度減小，但箱體結(jié)構(gòu)所受的應(yīng)力增大。