含隔板球形貯箱液體晃動激勵頻率影響分析

馬 亮， 劉 昊， 魏 承， 湯 亮， 趙 陽

(1.哈爾濱工業(yè)大學(xué) 航天工程系，哈爾濱 150001； 2.北京控制工程研究所，北京 100190)

光滑粒子流體動力學(xué)方法(Smoothed Particle Hydrodynamics,SPH)自1977年由Lucy等[1-2]首次提出，在經(jīng)歷了40年的發(fā)展和改進(jìn)過程后，現(xiàn)已在航空航天[3-5]、船舶海洋[6-8]、公路運(yùn)輸[9-11]等多個領(lǐng)域取得了廣泛應(yīng)用。

在航天領(lǐng)域，隨著現(xiàn)代航天器對高精度載荷及長壽命飛行的深入需求，航天器貯箱內(nèi)液體晃動所產(chǎn)生的晃動力對系統(tǒng)動力學(xué)行為的影響已經(jīng)成為總體設(shè)計(jì)重點(diǎn)考慮的問題。

基于SPH方法的流體動力學(xué)問題研究經(jīng)歷了逐步完善的過程。Monaghan[12]于1994年首次將SPH方法引入流體動力學(xué)問題的求解，隨后多位學(xué)者基于此方法對鈍體繞流[13]、剪切流動[14]、潰壩[15]等問題進(jìn)行了仿真分析。同時，SPH方法在計(jì)算精度[16]、計(jì)算穩(wěn)定性[17]及邊界處理[18]等方面經(jīng)過不斷地改進(jìn)和修正，使其成為一種解決流體動力學(xué)問題的成熟建模方法。

現(xiàn)代航天器貯箱多為球形結(jié)構(gòu)，而針對球形貯箱的液體晃動問題，現(xiàn)有研究多采用等效模型法或有限元法：黃華等[19]通過建立三維質(zhì)心面等效模型，利用質(zhì)心動力學(xué)方程及接觸碰撞模型對橢球形貯箱的晃動力進(jìn)行求解；趙陽等[20]采用任意拉格朗日-歐拉法(Arbitrary Lagrange-Euler,ALE)，對常重力及微重力環(huán)境下球形貯箱內(nèi)的液體晃動現(xiàn)象進(jìn)行了仿真分析；王天舒等[21]提出適用于微重力環(huán)境的球形貯箱液體晃動彈簧-質(zhì)量等效模型，同時利用虛功率原理對晃動力進(jìn)行了求解；苗楠等[22]建立了靜平衡面垂直于等效重力的等效模型，從而將液體大幅晃動分解為等效重力的整體運(yùn)動及在此基礎(chǔ)上的小幅晃動；岳寶增等[23]提出網(wǎng)格移動算法，在ALE方法基礎(chǔ)上對球形貯箱大幅晃動進(jìn)行了數(shù)值模擬。

由于航天器在執(zhí)行任務(wù)過程中，伴隨著頻繁的軌道及姿態(tài)機(jī)動，使得航天器貯箱存在大范圍運(yùn)動或大載荷受力，從而導(dǎo)致貯箱內(nèi)液體出現(xiàn)大幅度晃動?，F(xiàn)有研究所采用的等效模型方法[24]或有限元方法均無法適用于大幅晃動下液體出現(xiàn)破碎的現(xiàn)象，而針對航天器球形貯箱的液體大幅晃動問題又亟待解決。

通過以上分析，本文開展了基于SPH方法的球形貯箱液體晃動問題研究：首先根據(jù)SPH方法的基本方程，給出N-S方程的求解形式，同時確定人工黏度項(xiàng)及壓力項(xiàng)方程，并采用動態(tài)邊界以提高邊界計(jì)算精度；然后基于SPH方法建立球形貯箱液體晃動模型，采用兩種激勵模式及三種激勵頻率進(jìn)行仿真分析，分別記錄艙壁受力和觀測點(diǎn)壓強(qiáng)情況；最終設(shè)計(jì)一種帶隔板的球形貯箱結(jié)構(gòu)，并分析隔板對液體晃動的抑制作用。

1 基于SPH方法的N-S方程求解

光滑粒子流體動力學(xué)方法的基本思想是通過兩步近似，將空間函數(shù)及空間函數(shù)導(dǎo)數(shù)表示為核函數(shù)支持域內(nèi)各離散點(diǎn)處函數(shù)值與核函數(shù)及核函數(shù)導(dǎo)數(shù)乘積求和的形式。有利于簡化空間函數(shù)及空間函數(shù)導(dǎo)數(shù)求解的難度，使得該方法適用于求解N-S方程。

由Gingold等的研究可知，在經(jīng)歷核近似、粒子近似兩步近似過程后，空間函數(shù)f(x)在粒子i處的函數(shù)值可以表示為

(1)

式中:W為核函數(shù);h為光滑長度;N為核函數(shù)支持域內(nèi)粒子總數(shù)。

空間函數(shù)f(x)的導(dǎo)數(shù)在粒子i處的函數(shù)值可以表示為

(2)

核函數(shù)的支持域由光滑長度h及光滑因子κ決定，支持域半徑為κh，如圖1所示。本文采用三次樣條核函數(shù)[25]，其具體表達(dá)式為

式中:R=|x-x′|/h;光滑因子取為κ=2。

圖1 核函數(shù)支持域及粒子近似Fig.1 Support domain of kernel function and the particle approximation

采用SPH方法的基本方程，可將N-S方程組中的連續(xù)性方程、動量方程表示為

(4)

(5)

式中:vij=vi-vj;Wij=W(xi-xj,h);fe為單位質(zhì)量外力;α和β為坐標(biāo)軸方向。

(6)

在弱可壓縮模型下[27]，壓強(qiáng)項(xiàng)可表示為

(7)

式中:ρ0為流體參考密度；壓強(qiáng)-密度剛性系數(shù)γ=7;聲速cs取決于流體的最大速度，取為cs=10Vmax。

SPH方法存在邊界問題：在距邊界一個κh長度的范圍內(nèi)(包括邊界處)會出現(xiàn)核函數(shù)支持域與問題域交錯的現(xiàn)象。如圖2所示，方形粒子代表邊界粒子，三角形粒子代表發(fā)生截?cái)嗟牧黧w粒子，圓形粒子代表未發(fā)生截?cái)嗟牧黧w粒子。為使這一范圍內(nèi)的粒子計(jì)算保持精度，需對邊界問題進(jìn)行特殊處理。

為解決SPH方法的邊界問題，本文采用動態(tài)邊界法[28]，該方法的思想是直接構(gòu)造與流體實(shí)粒子性質(zhì)相同的邊界虛粒子，在計(jì)算過程中將邊界虛粒子代入計(jì)算，解決邊界附近(一個κh長度范圍內(nèi))流體粒子核函數(shù)支持域截?cái)鄦栴}。但邊界虛粒子自身參數(shù)不隨計(jì)算結(jié)果發(fā)生變化，而是保持固定或隨規(guī)劃運(yùn)動軌跡進(jìn)行運(yùn)動，該方法下虛粒子厚度與核函數(shù)支持域半徑相關(guān)。

圖2 支持域與問題域截?cái)鄬?dǎo)致邊界問題Fig.2 Boundary problem caused by the overlap of support domain and computational domain

2 仿真工況設(shè)計(jì)

采用SPH方法對N-S方程進(jìn)行離散化，并給出人工黏度項(xiàng)、壓力項(xiàng)的經(jīng)驗(yàn)方程，針對邊界問題提出動態(tài)邊界法，從而實(shí)現(xiàn)對N-S方程的完整求解。在此基礎(chǔ)上完成球形貯箱液體晃動仿真模型的建立，主要包括貯箱及液體建模、激勵形式及工況設(shè)計(jì)。

2.1 貯箱及液體

貯箱設(shè)計(jì)分為無隔板模型及有隔板模型，兩種結(jié)構(gòu)的唯一區(qū)別在于隔板結(jié)構(gòu)，其他參數(shù)均相同。球形貯箱半徑r=0.125 m，液體深度h=0.1 m，隔板高度hb=0.05 m，隔板位于xoz平面內(nèi)。為測量液體晃動對艙壁產(chǎn)生的壓力，設(shè)置了4個測量點(diǎn)，測量點(diǎn)均位于yoz平面的艙壁面上，不同測量點(diǎn)的高度不同。其中，p2點(diǎn)高度與隔板高度相同，p3點(diǎn)高度與液面高度相同，具體位置如圖3所示，貯箱與隔板均設(shè)計(jì)為剛體。

圖3 貯箱及液體建模(m)Fig.3 Model of spherical tank and fluid(m)

2.2 激 勵

由于航天器貯箱在實(shí)際飛行過程中存在大范圍運(yùn)動或大載荷受力，仿真設(shè)計(jì)兩種類型激勵模式：力模式、運(yùn)動模式。力模式的加載對象為液體，控制量為液體所受加速度；運(yùn)動模式的加載對象為貯箱，控制量為貯箱位移。兩種模式下均采用簡諧激勵形式：x=Asin(2πf·t+φ)，力模式下x表示液體y方向所受加速度，運(yùn)動模式下x表示貯箱y方向位移，仿真的重力環(huán)境選為常重力。設(shè)計(jì)6組仿真工況，激勵模式及簡諧運(yùn)動參數(shù)如表1所示。

表1 仿真工況設(shè)計(jì)Tab.1 Design of simulation and parameters

3 激勵頻率影響分析

圖4所示為兩種貯箱模型下，4個時刻點(diǎn)處晃動液體的側(cè)視圖。每個時刻下有隔板貯箱的液體晃動幅度均小于無隔板情況，為了定量分析隔板影響，分別從艙壁受力、測點(diǎn)壓強(qiáng)兩方面進(jìn)行對比分析。

3.1 艙壁受力

力模式激勵下液體的加載方向?yàn)閥方向，液體在yoz平面的晃動行為明顯，因此主要分析艙壁y,z兩方向的受力情況。

艙壁y方向初始受力為零，隨著液體受到簡諧激勵，開始出現(xiàn)周期性晃動力。不同的激勵頻率造成艙壁受力不同的曲線形式，當(dāng)激勵頻率f=1.5 Hz時晃動力幅度達(dá)到最大，而三種激勵頻率下隔板的存在均使得晃動幅度減小，如圖5(a)和圖5(b)所示。

艙壁z方向初始受力為負(fù)值，這是由液體重力造成的，隨著液體受到簡諧激勵，同樣出現(xiàn)周期性晃動力。不同的激勵頻率對艙壁z方向晃動力影響巨大，且同樣在激勵頻率f=1.5 Hz時晃動力幅度達(dá)到最大，隔板的存在同樣使得晃動幅度減小，如圖5(c)和圖5(d)所示。可以看出，由于存在液體晃動現(xiàn)象，艙壁在y,z兩方向的受力均增大，且為周期性受力。

表2給出了各工況下，艙壁受力穩(wěn)定后的平均晃動幅值，體現(xiàn)出的規(guī)律與曲線圖一致，同時可以看出隔板對激勵頻率大的工況晃動抑制效果更明顯。

圖4 晃動過程中各時刻下液面形態(tài)Fig.4 Shape of the fluid at different moment in sloshing

圖5 艙壁受力隨振動頻率變化情況Fig.5 Force applied on bulkhead varies with vibrational frequency

表2 艙壁受力平均晃動幅值Tab.2 Average force amplitude applied on bulkhead

3.2 測點(diǎn)壓強(qiáng)

運(yùn)動模式激勵下貯箱的加載方向同樣為y方向，液體在yoz平面的晃動行為明顯，因此在yoz平面設(shè)置了4個壓力測量點(diǎn)：p1點(diǎn)高度低于隔板高度，p2點(diǎn)高度與隔板高度相同，p3點(diǎn)高度與液面高度相同，p4點(diǎn)高度高于液面。

對于p1點(diǎn)，在無隔板工況下，不同的激勵頻率使得壓力曲線呈現(xiàn)出截然不同的形式：當(dāng)f=1 Hz時，晃動行為與激勵形式一致，表現(xiàn)為簡諧波動；當(dāng)f=1.5 Hz時，出現(xiàn)“雙波峰”現(xiàn)象，前壓力峰值略高于后壓力峰值，這是由液體撞擊艙壁后沿球形內(nèi)輪廓上升后回落所形成；當(dāng)f=2 Hz時，也會出現(xiàn)“雙波峰”現(xiàn)象，但兩次壓力峰值相差較大，這是由于液體撞擊艙壁后的上升幅度較小，沒有形成翻卷、破碎，如圖6(a)所示。在有隔板工況下，三種激勵頻率下的測點(diǎn)壓力值均略有下降。當(dāng)f=1 Hz時，壓力曲線形式同樣表現(xiàn)為簡諧波動；當(dāng)f=1.5 Hz時，“雙波峰”的前后壓力峰值更加接近；當(dāng)f=2 Hz時，后壓力峰值反高于前壓力峰值，如圖6(b)所示。這是由于p1點(diǎn)位于隔板下方，隔板的存在會對液體首次沖擊艙壁存在阻礙作用，而當(dāng)液體回落時隔板同樣阻礙了液體的回流，正是這種作用使得前壓力峰值有所減小而后壓力峰值明顯增大。但隔板的存在使得壓力峰值達(dá)到“平均化”效果，這樣不會使艙壁局部在某時刻的壓力過大，而導(dǎo)致結(jié)構(gòu)變形或損壞。

對于p2點(diǎn)，其高度高于p1點(diǎn)，但仍位于液面之下，因此在有無隔板及三種激勵頻率下的壓力曲線形式與p1點(diǎn)都相似，但局部現(xiàn)象存在不同之處。在無隔板及激勵頻率為1.5 Hz,2 Hz的工況下，沒有出現(xiàn)明顯的“雙波峰”，而變?yōu)榍安ǚ寮印昂髩毫ζ矫妗钡慕M合,如圖6(c)所示。這是由于p2點(diǎn)高度增加，液體不存在回落或回落幅度極??；在有隔板及激勵頻率為1.5 Hz,2 Hz的工況下，前壓力峰值高于后壓力峰值，如圖6(d)所示。這是由于p2點(diǎn)高度與隔板高度相同，隔板對液體首次沖擊艙壁及液體回落回流的阻礙作用都減小，但隔板的存在同樣使得測點(diǎn)壓力值有所下降。

對于p3點(diǎn)，其高度位于液面處，壓力曲線形式與前兩個位于液面下的測點(diǎn)存在較大差異。在三種激勵頻率下壓力曲線均呈現(xiàn)“單波峰”加零壓力組合。當(dāng)不存在隔板時，激勵頻率越大壓力峰值越大,如圖6(e)所示；當(dāng)存在隔板時，隨著激勵頻率的增加隔板對測點(diǎn)處壓力的抑制效果越明顯，如圖6(f)所示。

對于p4點(diǎn)，其高度位于液面之上，壓力曲線形式與p3點(diǎn)相似。在激勵頻率為1 Hz工況下已檢測不到壓力值，而在激勵頻率為1.5 Hz,2 Hz的工況下，壓力曲線也變?yōu)椤懊}沖波”加零壓力組合。隔板的存在同樣對激勵頻率越大的工況測點(diǎn)處壓力抑制越明顯,如圖6(g)和圖6(h)所示。

圖6 測點(diǎn)壓強(qiáng)隨振動頻率變化情況Fig.6 Pressure of measuring point varies with vibrational frequency

4 結(jié) 論

通過建立有無隔板兩種球形貯箱模型，針對兩種激勵模式及三種激勵頻率進(jìn)行仿真分析，最終得到以下結(jié)論：

(1) 力模式激勵下，無論貯箱有無隔板均在激勵頻率f=1.5 Hz時，得到艙壁受力的最大峰值，而隔板的存在會減弱艙壁所受晃動力。

(2) 力模式激勵下，由于液體晃動現(xiàn)象的產(chǎn)生，艙壁在y,z兩方向的受力均增大，且為周期性受力。

(3) 運(yùn)動模式激勵下，位于液面下的艙壁測點(diǎn)p1,p2表現(xiàn)出一致性，液體沖擊艙壁及回落作用導(dǎo)致“雙波峰”現(xiàn)象的出現(xiàn)，測點(diǎn)高度的不同又會使得該現(xiàn)象在兩測點(diǎn)處存在差異。而隔板的存在對艙壁所受壓力具有“平均化”效果。

(4) 運(yùn)動模式激勵下，位于液面處及液面以上的艙壁測點(diǎn)p3,p4表現(xiàn)出一致性，出現(xiàn)“單波峰”或“脈沖波”加零壓力組合形式，且隔板的存在對激勵頻率越大的工況抑制效果越明顯。