高空氣球發(fā)放過程中的流固耦合動(dòng)力學(xué)研究

蘇 一,張向強(qiáng),劉 強(qiáng)

(1. 中國科學(xué)院空天信息創(chuàng)新研究院,北京 100094;2. 中國科學(xué)院大學(xué),北京 100049)

1 引言

臨近空間(Near Space)通常被定義為距離地面18-100km高度的高空區(qū)域,其位于地球大氣中的平流層和中層,是一個(gè)待開拓的具有重要戰(zhàn)略意義的空間圈層。在這一空間圈層的中低層(18-50km)運(yùn)用最為廣泛的平臺(tái)是高空氣球,高空氣球根據(jù)球膜是否承受內(nèi)外壓差通常被劃分為零壓氣球和超壓氣球,從上世紀(jì)五六十年代開始,零壓氣球由于其成本低、載重大、發(fā)放簡易等優(yōu)勢,在空間科學(xué)領(lǐng)域方面發(fā)揮了重要的作用[5-7];1997年,NASA提出了ULDB計(jì)劃,其目的是完成超壓氣球載重1.6噸、駐空33.5km,續(xù)航100天的飛行任務(wù),目前已實(shí)現(xiàn)駐空54天的飛行[8]。由于臨近空間的大氣密度較小,海拔20km的大氣密度是地面的1/14,海拔35km的大氣密度是地面的1/150,即為維持浮力重力平衡,高空氣球的體積很大,同時(shí)材料面密度要小或厚度薄。目前國內(nèi)外常用的球膜材料一般為聚乙烯,體積達(dá)幾萬至幾十萬立方米,球膜厚度在20～80微米之間。若設(shè)計(jì)升限為40～50km,材料的厚度降為幾微米左右[13]。因此,基于氣球的形變和應(yīng)力分析的球形設(shè)計(jì)始終是高空氣球的難點(diǎn)[14]。

高空氣球發(fā)放方式大體上可分為靜態(tài)發(fā)放、半動(dòng)態(tài)發(fā)放和動(dòng)態(tài)發(fā)放[2-4],不管采用哪種發(fā)放方式,解除約束后氣球受力和形態(tài)變化劇烈,典型的氣球發(fā)放過程如圖1所示。

圖1 高空氣球發(fā)放過程

2 發(fā)放過程系纜分析

2.1 數(shù)學(xué)模型

假設(shè)氣泡(球體充氣部分)剛性,結(jié)纜系統(tǒng)和氣泡以下的球體部分視為柔性系纜,簡化后的模型受力如下圖所示。

考慮在X,Y坐標(biāo)下氣球某位置上球體微段的受力,它包括張力T和T′=T-dT,重量dW.g,風(fēng)對該線段的阻力dDr,在瞬時(shí)位置要求每一個(gè)微段受力平衡,在系纜沒有觸地的時(shí)候滿足以下平衡方程[1]:

X方向有,

T·cosθ+dDr·cosα=(T-dT)·cos(θ-dθ)

(1)

Y方向有,

T·sinθ=dW·g+dDr·sinα+(T-dT)·sin(θ-dθ)

(2)

化簡得,

dT=dDr·cosα·cosθ+(dW·g+dDr·sinα)·sinθ

(3)

另外有

dx=cosθ·dl

dy=sinθ·dl

(4)

(5)

ODE45是一種自適應(yīng)步長的采用Runge-Kutta算法常微分方程數(shù)值解法,給定微分方程、積分范圍和參數(shù)初值可以得到積分范圍內(nèi)的解析解。積分的初始值為[T0θ0X0Y0],X0=0,Y0由高空氣球所處的位置決定。通過打靶法,并以二分法優(yōu)化初值的搜索過程,初值T0和θ0建立以下的關(guān)聯(lián)

T0=F·sinθ0

(6)

其中,F為高空氣球產(chǎn)生的凈浮力,這樣T0的初始值由θ0決定,初始值變?yōu)閇F·sinθ0θ0X0Y0],變量只有一個(gè),確定θ0的取值范圍,通過打靶法找到確定的θ0,就可以完成曲線的擬合。

假設(shè)有

dDr=fDr·dl

dW=ρ·dl

(7)

其中fDr系纜單位長度所受風(fēng)阻,ρ為線密度,根據(jù)(7)化簡(3),這樣積分時(shí)以dl為自變量,積分區(qū)間為系纜長度,根據(jù)公式,并考慮初值的正負(fù)性,則ODE45初值公式為

=-fDr·cosα·cosθ+(ρ·g-fDr·sinα)·sinθ

(8)

其中α和θ之間有

(9)

當(dāng)系纜觸地后每一個(gè)微段會(huì)有

T′+Ff=T

(10)

Ff=f·FN

(11)

根據(jù)圖2可知,觸地后系纜微段的微分方程為

圖2 系纜受力圖

(12)

確定系纜的曲線后,考慮系纜是否觸地,將曲線做出后和地面的高度比較,判斷系纜曲線是否與地面直線相交,若不相交,則系纜曲線為所求;若相交,則分段打靶求曲線。第一段打靶根據(jù)初始條件找出滿足與地面相切的點(diǎn)C1,相切條件為目標(biāo)點(diǎn)C1(y=0,θ=0),然后根據(jù)觸地后的微分方程,根據(jù)區(qū)間設(shè)定點(diǎn)C2的范圍為(C2,O),作為第三段打靶的初始條件,打靶最后收斂于目標(biāo)點(diǎn)O,反過來確定C2的位置。

高空氣球在發(fā)放過程中會(huì)受到來自空氣的阻力

(13)

其中,Fair為空氣阻力,ρa(bǔ)ir為空氣密度,V為高空氣球體積,Cd空氣阻力系數(shù),Sd為高空氣球受力面積。由上文可知,每次可以通過確定的位置得到當(dāng)時(shí)球體的受力情況,如圖3所示,確定系纜的曲線初值為[F·sinθ0θ0X0Y0],而只需要給定Y0就可以通過打靶法確定曲線,故只需分析Y方向的受力情況,根據(jù)基本的運(yùn)動(dòng)學(xué)公式

圖3 高空氣球受力圖

(14)

Y方向受力情況為

F=T·sinθ+Fair+(MHe+Mballoon)·a

(15)

本次數(shù)學(xué)模型采用VBalloon=7m3的小型高空氣球,充入30%的氦氣作為浮升氣體,經(jīng)過

FHe=ρa(bǔ)ir·VHe·g

GHe=ρHe·VHe·g

(16)

F=FHe-GHe-MBalloon·g

(17)

其中,VHe=0.3·VBalloon,經(jīng)過計(jì)算可知F=22.12N,因?yàn)槌錃饬看嬖谡`差并且球膜的質(zhì)量不可控,經(jīng)測量凈浮力為21.52N。

通過MATLAB中ODE45函數(shù)求解后,高空氣球在高度Y方向的加速度、速度、位移和系纜末端O′點(diǎn)拉力隨時(shí)間的變化如圖4所示。

圖4 高空氣球運(yùn)動(dòng)曲線

由高空氣球運(yùn)動(dòng)軌跡可知,氣球在初始時(shí)加速度較大,速度在到達(dá)峰值后逐漸變小。

2.2 有限元仿真

通過MATLAB中的計(jì)算結(jié)果可知,在初始時(shí)刻發(fā)放車和高空球底部的直線距離為20.76m,在DYNA的前處理中將系纜未與高空氣球相連的一端通過添加初始位移的方法將距離調(diào)至相同,并通過添加重力達(dá)到發(fā)放的初始狀態(tài)的參數(shù)設(shè)置相同距離,除此之外,高空氣球高度、發(fā)放車高度、系纜的粗細(xì)、長度和線密度以及系纜所受的重力和高空氣球產(chǎn)生的浮力均與數(shù)學(xué)模型一致,在LS-DYNA中建立模型如圖5所示。

圖5 高空氣球系纜隨時(shí)間變化(間隔1s)

圖6 LS-DYNA模型初始狀態(tài)

圖7 LS-DYNA仿真系纜形狀變化(間隔1s)

高空氣球在上升中受浮力和重力相互影響,在LS-DYNA中取全局阻尼為3.5,計(jì)算結(jié)果如圖8～11所示。

圖8 氣球加速度隨時(shí)間變化

圖9 氣球速度隨時(shí)間變化

圖10 氣球高度隨時(shí)間變化

圖11 系纜拉力隨時(shí)間變化

圖12 高空氣球受力分析

圖13 高空氣球初始形狀

通過LS-DYNA仿真可以看出,MATLAB和LS-DYNA計(jì)算結(jié)果說明高空氣球在發(fā)放過程中確實(shí)會(huì)受到系纜拉力作用,但LS-DYNA中無法做到根據(jù)用戶要求定義阻尼大小,故速度峰值會(huì)有所不同。

3 發(fā)放仿真

在LS-DYNA中進(jìn)行動(dòng)力學(xué)仿真,觀察球膜的形狀以及應(yīng)力變化,因此球膜需要轉(zhuǎn)化成柔性體,通過在球膜里面加入SPH粒子來模擬氦氣,通過SPH和球膜的相互作用,對發(fā)生沖擊時(shí)球膜和氣體進(jìn)行耦合分析。

3.1 高空氣球球型模型

氣球球膜一般由厚度為20μm左右的薄膜制成,球膜的厚度和球體尺度相比很小,符合無矩薄殼理論使用的條件[12]。

由回轉(zhuǎn)薄殼無力矩理論可知:

(18)

其中

(19)

(20)

將式(19)(20)代入(18),并化簡得

(21)

在截面圓上的總負(fù)載T=2πrσm在弧長方向的改變程度和周向應(yīng)力在弧長方向的分力以及膜重的和應(yīng)保持平衡,如式(22)所示。

(22)

化簡式(22)得

(23)

球形參數(shù)λ和∑的定義如式(24)所示,λ和Σ表示了氣球負(fù)載L及飛行高度對球形的影響程度

(24)

自然形高空氣球應(yīng)當(dāng)保證周向應(yīng)力σc=0,上式共同組成球形母線的求解微分方程組,在NASA高空氣球手冊中有Σ表,不同的Σ值對應(yīng)不同的母線形狀,本文采用的Σ=0.3,成型后的形狀圖14所示。

圖14 高空氣球柔性球膜發(fā)放過程

3.2 SPH方法

SPH方法(Smoothed Particle Hydrodynamics)全稱為光滑粒子流體動(dòng)力學(xué)方法[10],是一種為求解流體的動(dòng)力學(xué)問題而提出的一種無網(wǎng)格方法。一般流體動(dòng)力學(xué)問題需要基于密度、速度、能量等變量場的偏微分方程組,而SPH方法是通過采用插值方法將連續(xù)的流體用相互作用的質(zhì)點(diǎn)組來描述,每個(gè)質(zhì)點(diǎn)上包含各種物理量的信息,這樣再過連續(xù)介質(zhì)力學(xué)的守恒定律可以將整個(gè)流體過程描述出來,從這個(gè)方面來說,只要采用的質(zhì)點(diǎn)粒子數(shù)夠多,就能準(zhǔn)確的描述流體的力學(xué)行為。

在有限元計(jì)算中,SPH是一種純Lagrange方法,SPH由于不存在網(wǎng)格關(guān)系,所以不存在由于變形過大導(dǎo)致的網(wǎng)格畸變問題,因此適合求解高速碰撞過程中的大變形問題。

3.3 接觸算法

LS-dyna在處理接觸問題的過程中主要采用以下三種不同的算法:

1)動(dòng)力約束法

2)分配參數(shù)法

3)對稱罰函數(shù)法。

LS-DYNA的一般默認(rèn)采用對稱罰函數(shù)法,是計(jì)算中最常用的一種方法。

其首先要確定接觸面的主從關(guān)系,一般從接觸面比主接觸面的網(wǎng)格更密,然后在計(jì)算過程中的每一個(gè)時(shí)間步都判斷從節(jié)點(diǎn)是否會(huì)穿透主表面,沒穿透則不操作;若穿透,則對從節(jié)點(diǎn)和接觸面之間引入一種限制穿透作用的類似法向彈簧的接觸力,接觸力的大小和穿透深度、接觸剛度有關(guān),具體如式(25)所示

F=K·δ

(25)

其中,K為接觸界面剛度,δ為穿透量。這種接觸算法應(yīng)用簡單,不需要添加碰撞和釋放條件,而且計(jì)算過程中保持動(dòng)量守恒。在球膜與氦氣耦合過程中,氦氣為從節(jié)點(diǎn),球膜為主接觸面[9]。

3.4 沖擊仿真

3.4.1 材料參數(shù)

表1 仿真參數(shù)表

3.4.2 高空氣球發(fā)放沖擊模擬

本文研究的高空氣球的發(fā)放過程是高空氣球在充氣完畢后在滾筒的約束下保持穩(wěn)定為發(fā)放零時(shí)刻,之后打開發(fā)放滾筒到氣球升起并在升到系纜的極限長度后在下方發(fā)放車的約束下停止并保持穩(wěn)定的一個(gè)過程。

3.4.2.1 建立模型

首先根據(jù)Σ=0.3做出球型的母線,將母線導(dǎo)入U(xiǎn)G建立模型,再用ICEM為模型畫網(wǎng)格,用LS-prepost進(jìn)行前處理,并加入SPH粒子,建立好的網(wǎng)格模型如下圖所示。其中高空氣球體積為7m3,沖入的氦氣為氣球體積的30%,然后對氦氣施加浮力,對球膜和系纜施加重力,獲得初始穩(wěn)定的形狀。

利用文章2.2章節(jié)中所采用的仿真方法,如圖5所示,將球膜轉(zhuǎn)換成柔性體并沖入30%氦氣,會(huì)得到如圖15所示的球膜以及系纜曲線。

圖15 高空氣球充氦氣后形狀與實(shí)際對比圖

發(fā)放過程中會(huì)出現(xiàn)沖擊是在系纜拉直,球膜會(huì)向系纜方向收緊,和豎直向上所受沖擊模式類似,故可以采用豎直向上發(fā)放的方式研究受到系纜和氣體沖擊時(shí)球膜所受應(yīng)力變化。

3.4.2.2 添加場力

因?yàn)榻⒌哪Ｐ屯獠坎]有氣體,所以要讓氦氣產(chǎn)生向上的浮力氣體的密度并非氦氣的密度,而是通過浮力公式進(jìn)行換算,算出的密度應(yīng)當(dāng)將外部氣體和氦氣的密度都考慮進(jìn)去,根據(jù)式(16),有

F凈浮力=FHe-GHe=ρa(bǔ)ir·VHe·g-ρHe·VHe·g

=ρHe′·VHe·g

(26)

故有

ρHe′=ρa(bǔ)ir-ρHe

(27)

這樣,預(yù)先充入的氣體在受到向上的加速度g時(shí),可以產(chǎn)生向上的且符合實(shí)際的浮力,球膜和系纜通過關(guān)鍵字*LOAD＿GRAVITY＿PART施加重力,獲得初始形狀如圖16所示。

圖16 高空氣球沖擊受力圖

高空氣球成型后的形狀是內(nèi)部流場和外部流場對球膜共同作用的結(jié)果,由于仿真后的高空氣球模型沒有加入外部流場,故來自于外部流場的壓強(qiáng)和風(fēng)場并沒有作用在球膜上,圖15中球膜的形狀是在內(nèi)部流場作用下穩(wěn)定的結(jié)果。

3.4.2.3 后處理

獲得初始穩(wěn)定形狀后,仿真模型在穩(wěn)定形狀的基礎(chǔ)上給其中的SPH粒子和球膜節(jié)點(diǎn)通過關(guān)鍵字*INITIAL＿VELOSITY施加初速度4m/s,在仿真時(shí)不施加阻尼。 以下為沖擊過程中兩次沖擊間隔過程中球膜的應(yīng)力以及形狀變化。

通過球膜的受力情況可知,開始時(shí)球頂中心有比較大的應(yīng)力,之后隨著氣體擴(kuò)散,球頂部受力比較均勻,但較大的應(yīng)力還是集中在球頂和氣球底部。當(dāng)氣體和球膜沖擊后會(huì)造成球膜有比較大形狀變化,并由于氣體運(yùn)動(dòng)的關(guān)系和球膜的接觸變少,所以此時(shí)球膜應(yīng)力并不大。在浮力的作用下,氣體再次上升,會(huì)和球膜再一次發(fā)生碰撞,這時(shí)球頂部和球底部應(yīng)力較大,通過應(yīng)力云圖發(fā)現(xiàn)球底部應(yīng)力較頂部更大。

取高空氣球頂部中心的單元,做出其速度和應(yīng)力曲線如圖17所示。

圖17 球頂?shù)乃俣扰c應(yīng)力隨時(shí)間變圖

由圖中應(yīng)力和速度曲線可知,在速度下降的同時(shí)伴隨著球膜應(yīng)力的上升,說明氣體沖擊會(huì)造成球膜應(yīng)力增加,同時(shí)球膜會(huì)降低氣體的速度,說明球膜和氣體產(chǎn)生了耦合作用。

4 結(jié)論

本文以7m3高空氣球?yàn)槔?通過MATLAB和LS-DYNA對高空氣球的發(fā)放過程中系纜隨時(shí)間的變化以及受到?jīng)_擊時(shí)球膜形狀及應(yīng)力的變化的進(jìn)行了模擬。通過MATLAB的數(shù)值仿真可以得到系纜的拉力,并可以對高空氣球做受力分析,通過運(yùn)動(dòng)學(xué)方程得到每一時(shí)刻高空氣球的位置,并得到高空氣球每個(gè)位置的受力情況,通過LS-DYNA將簡化后的模型進(jìn)行了有限元仿真,并對受到?jīng)_擊時(shí)的球膜的形狀及應(yīng)力變化進(jìn)行了仿真,通過對比可以得到以下結(jié)論:

1)高空氣球發(fā)放時(shí)會(huì)產(chǎn)生2～3g的加速度,這時(shí)候由于氦氣和球膜的初始狀態(tài)均是靜止的,相對速度不大,并且系纜并沒有拉緊,初始狀態(tài)更接近于高空氣球做自由上升運(yùn)動(dòng),系纜對球膜的拉力隨著高空氣球高度的增加逐漸增大,阻力也隨著球膜的速度的增大而增大,這樣高空氣球在阻力和系纜拉力的共同作用下速度上升不會(huì)太快,能有效減小到達(dá)沖擊點(diǎn)時(shí)的速度,達(dá)到減小沖擊的效果。

2)通過LS-DYNA將發(fā)放時(shí)球膜所受沖擊進(jìn)行仿真,通過應(yīng)力云圖可知?dú)馇蝽敳亢蜌馇虻撞繛檩^大的應(yīng)力集中區(qū)域,而由于球頂部受力面積較大,應(yīng)力相對較小,而受力更加集中的底部反而是更容易遭到破壞的地方。在以后球頂部中心要加強(qiáng)材料強(qiáng)度,而球底部要加強(qiáng)對球膜的保護(hù),避免應(yīng)力集中。

在以后計(jì)算高空氣球發(fā)放過程中,可以采用MATLAB數(shù)值仿真的方法,改變系纜長度,高空氣氣球的大小,和充氣量的不同,計(jì)算不同情況下高空氣球的速度峰值。通過LS-DYNA可以仿真小型高空氣球球膜應(yīng)力及形狀變化,但難以采用SPH方法對大型高空氣球進(jìn)行模擬,因?yàn)殡S著高空氣球的體積增大,網(wǎng)格數(shù)量和SPH粒子數(shù)量也在增加,導(dǎo)致計(jì)算量呈幾何級增長,需要找到更加合適的計(jì)算方法。