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

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

1 引言

臨近空間(Near Space)通常被定義為距離地面18-100km高度的高空區(qū)域,其位于地球大氣中的平流層和中層,是一個(gè)待開(kāi)拓的具有重要戰(zhàn)略意義的空間圈層。在這一空間圈層的中低層(18-50km)運(yùn)用最為廣泛的平臺(tái)是高空氣球,高空氣球根據(jù)球膜是否承受內(nèi)外壓差通常被劃分為零壓氣球和超壓氣球,從上世紀(jì)五六十年代開(kāi)始,零壓氣球由于其成本低、載重大、發(fā)放簡(jiǎn)易等優(yōu)勢(shì),在空間科學(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í)材料面密度要小或厚度薄。目前國(guó)內(nèi)外常用的球膜材料一般為聚乙烯,體積達(dá)幾萬(wàn)至幾十萬(wàn)立方米,球膜厚度在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ā)放過(guò)程如圖1所示。

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

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

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

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

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

X方向有,

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

(1)

Y方向有,

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

(2)

化簡(jiǎn)得,

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

(3)

另外有

dx=cosθ·dl

dy=sinθ·dl

(4)

(5)

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

T0=F·sinθ0

(6)

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

假設(shè)有

dDr=fDr·dl

dW=ρ·dl

(7)

其中fDr系纜單位長(zhǎng)度所受風(fēng)阻,ρ為線密度,根據(jù)(7)化簡(jiǎn)(3),這樣積分時(shí)以dl為自變量,積分區(qū)間為系纜長(zhǎng)度,根據(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,反過(guò)來(lái)確定C2的位置。

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

(13)

其中,Fair為空氣阻力,ρa(bǔ)ir為空氣密度,V為高空氣球體積,Cd空氣阻力系數(shù),Sd為高空氣球受力面積。由上文可知,每次可以通過(guò)確定的位置得到當(dāng)時(shí)球體的受力情況,如圖3所示,確定系纜的曲線初值為[F·sinθ0θ0X0Y0],而只需要給定Y0就可以通過(guò)打靶法確定曲線,故只需分析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)過(guò)

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

GHe=ρHe·VHe·g

(16)

F=FHe-GHe-MBalloon·g

(17)

其中,VHe=0.3·VBalloon,經(jīng)過(guò)計(jì)算可知F=22.12N,因?yàn)槌錃饬看嬖谡`差并且球膜的質(zhì)量不可控,經(jīng)測(cè)量?jī)舾×?1.52N。

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

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

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

2.2 有限元仿真

通過(guò)MATLAB中的計(jì)算結(jié)果可知,在初始時(shí)刻發(fā)放車和高空球底部的直線距離為20.76m,在DYNA的前處理中將系纜未與高空氣球相連的一端通過(guò)添加初始位移的方法將距離調(diào)至相同,并通過(guò)添加重力達(dá)到發(fā)放的初始狀態(tài)的參數(shù)設(shè)置相同距離,除此之外,高空氣球高度、發(fā)放車高度、系纜的粗細(xì)、長(zhǎng)度和線密度以及系纜所受的重力和高空氣球產(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 高空氣球初始形狀

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

3 發(fā)放仿真

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

3.1 高空氣球球型模型

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

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

(18)

其中

(19)

(20)

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

(21)

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

(22)

化簡(jiǎn)式(22)得

(23)

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

(24)

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

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

3.2 SPH方法

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

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

3.3 接觸算法

LS-dyna在處理接觸問(wèn)題的過(guò)程中主要采用以下三種不同的算法:

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

2)分配參數(shù)法

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

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

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

F=K·δ

(25)

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

3.4 沖擊仿真

3.4.1 材料參數(shù)

表1 仿真參數(shù)表

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

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

3.4.2.1 建立模型

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

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

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

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

3.4.2.2 添加場(chǎng)力

因?yàn)榻⒌哪Ｐ屯獠坎](méi)有氣體,所以要讓氦氣產(chǎn)生向上的浮力氣體的密度并非氦氣的密度,而是通過(guò)浮力公式進(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ò)關(guān)鍵字*LOAD＿GRAVITY＿PART施加重力,獲得初始形狀如圖16所示。

圖16 高空氣球沖擊受力圖

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

3.4.2.3 后處理

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

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

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

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

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

4 結(jié)論

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

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

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

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