楊燕初，張航悅，趙 榮

(1. 中國科學(xué)院大學(xué)， 北京 100190； 2. 中國科學(xué)院 光電研究院， 北京 100094)

高空氣球是利用密度小于空氣的氣體產(chǎn)生的浮力進(jìn)行飛行，屬于臨近空間低動態(tài)飛行器，具有飛行高度高、速度低、載重大的特點(diǎn)，為高空科學(xué)實驗研究提供了可靠的平臺。目前，高空氣球根據(jù)內(nèi)外壓差大小可分為零壓氣球和超壓氣球，零壓氣球底部有排氣管與大氣相通，超壓氣球則與外界封閉，可以認(rèn)為超壓氣球是零壓氣球在達(dá)到升限高度后進(jìn)一步上升膨脹形成的[1]。零壓氣球外形為軸對稱球形，且目前大多采用了自然形球形的設(shè)計。從1782年氣球第一次升空之后的150多年里，并沒有進(jìn)行過系統(tǒng)理論性的研究，工程師們依靠直覺選擇了諸如正球形、圓柱形及四面體形等形狀[1]。直到1939年，Upson提出了自然形氣球的概念[2]，規(guī)定載荷完全靠子午線方向的球膜應(yīng)力承擔(dān)，周向應(yīng)力始終為零，并且初步推導(dǎo)出氣球形狀方程。20世紀(jì)50年代，研究人員對軸對稱自然形高空氣球球形系統(tǒng)性地建立了數(shù)學(xué)模型并數(shù)值求解得到一組球形數(shù)據(jù)，推動了氣球球形的研究[3]。1965年，Smalley進(jìn)一步使用計算機(jī)數(shù)值求解得到了較為精確和完整的球形數(shù)據(jù)，并在氣球設(shè)計中得到了實際運(yùn)用[4]。之后的球形設(shè)計都是在自然形氣球的基礎(chǔ)上進(jìn)行改進(jìn)的，包括Li等在20世紀(jì)80年代提出的“混合形”球形[5]以及Yajima等提出的“三維球幅”設(shè)計[6]。

自然形球形方程最終可以化簡為一組微分方程，在求解方法上，通常是選擇標(biāo)準(zhǔn)的Runge-Kutta方法，姜魯華從提高計算精度方面選擇了Gill法[7]，Baginski等利用打靶法得到氣球外形數(shù)據(jù)并指出使用多段打靶法求解部分膨脹狀態(tài)下球形的可行性[8]。在此基礎(chǔ)上，本文提出將多段打靶法與最優(yōu)化方法中的等式二次規(guī)劃(Equality constrained Quadratic Programs, EQP)算法結(jié)合，進(jìn)而求解出完全膨脹及部分膨脹狀態(tài)下的氣球母線形狀。

高空氣球球膜在漂浮高度完全張開，不存在褶皺。在上升過程中，由于球體體積小于漂浮高度時的體積，球膜未完全張開，存在褶皺現(xiàn)象，較小體積時氣球底部球膜合攏形成聚束。對球膜褶皺及聚攏的分析涉及改變材料本構(gòu)關(guān)系及褶皺判斷等研究[9]，球形設(shè)計上通常選擇對其進(jìn)行一定的簡化。Smalley引入了“ad hoc”假設(shè)，認(rèn)為未完全膨脹狀態(tài)下多余的氣球球膜均勻地進(jìn)行軸對稱分布，并將球底球膜形成的聚束附加到負(fù)載質(zhì)量中，只計算上方升力氣體體積部分氣球形狀[4]。Baginski等則以同一母線長度處完全膨脹球形與部分膨脹球形截面半徑的比值作為球膜密度變化因子，改變部分膨脹狀態(tài)氣球膜密度以使得球體質(zhì)量在兩種狀態(tài)下都保持恒定，進(jìn)而忽略褶皺及多余材料的影響，得到氣球在部分膨脹狀態(tài)下球膜同樣完全張開的球形[8]。本文亦采用了Baginski的處理方法，將球膜面密度作為變量處理，對部分膨脹狀態(tài)下氣球整體形狀進(jìn)行分析。

1 自然形球形方程推導(dǎo)及數(shù)值求解

1.1 自然形球形母線微分方程組推導(dǎo)

氣球球膜一般由厚度為20 μm左右的薄膜制成，厚度遠(yuǎn)小于球體尺度，符合無矩薄殼理論。球形方程推導(dǎo)需使用的相關(guān)參數(shù)定義見表1，并且取球膜單元進(jìn)行受力分析，如圖1所示，其中單元經(jīng)放大處理。

表1 相關(guān)參數(shù)定義

圖1 氣球球膜單元Fig.1 Balloon film unit

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

(1)

將其代入式(1)并化簡得到：

(2)

在截面圓上的總負(fù)載T滿足T=2πrσm，在弧長方向的改變量d(2πrσm)/ds與周向應(yīng)力在弧長方向的分量以及膜重存在如式(3)所示關(guān)系。

(3)

化簡式(3)可得：

(4)

球形參數(shù)λ和Σ的定義如式(5)所示，λ、Σ表示了氣球負(fù)載L及飛行高度(決定bd值)對球形的影響程度。

(5)

自然形氣球規(guī)定周向應(yīng)力σc=0。以上各式共同組成球形微分方程組。并使用Simulink搭建仿真程序，如圖2所示。

圖2 Simulink結(jié)構(gòu)圖Fig.2 Simulink structure diagram

1.2 球形方程組數(shù)值求解

對于完全膨脹狀態(tài)的氣球，底部起始角度較大，且有邊界條件r0=z0=0，rs=0，r0、z0為氣球底點(diǎn)半徑及高度，rs為氣球頂點(diǎn)半徑，運(yùn)用一般打靶法將兩點(diǎn)邊值問題轉(zhuǎn)換為初值問題進(jìn)行求解[10]。此時猜測變量為氣球母線弧長s及底部起始角度θ0。

(6)

球形底點(diǎn)與頂點(diǎn)兩點(diǎn)邊界條件表示為g(z(0),z(l))=0，f、g為非線性函數(shù)。將邊值問題轉(zhuǎn)換為初值問題，設(shè)猜測初值為c0。

(7)

打靶法為尋找滿足頂點(diǎn)邊界條件的c0，通過不斷改變初值c0使得式(8)得以滿足。

g(c0,y(l;c0))=0

(8)

則z=y(s,c0)為非線性微分方程組的解。

(9)

在以上多段打靶法的基礎(chǔ)上，引入EQP算法，即將原始問題轉(zhuǎn)換為一系列等式約束二次規(guī)劃問題。將多段打靶法融入EQP二次規(guī)劃問題中，目標(biāo)函數(shù)設(shè)為常數(shù)，尋找滿足打靶法分段約束以及球形邊界條件約束的各段起始猜測點(diǎn)。

EQP計算過程中當(dāng)?shù)玫揭粋€迭代點(diǎn)ck，使用ck處的二次規(guī)劃模型確定ck+1，約束條件為約束函數(shù)G(c)在ck處Taylor展開的線性部分。

(10)

s.t. [G(ck)]Td+G(ck)=0

(11)

其中,d=c-ck，目標(biāo)函數(shù)O(c)=0；Hk為增廣Lagrange函數(shù)的Hesse矩陣；步長固定為1，即ck+1=ck+d。

對于上升過程中的氣球，按照式(12)改變對應(yīng)弧長截面的膜密度，以滿足在兩種狀態(tài)下氣球膜總質(zhì)量相等的要求。

(12)

式中，r、R分別為完全膨脹狀態(tài)及部分膨脹狀態(tài)球形截面圓半徑。求解部分膨脹狀態(tài)時的猜測變量為氣球零壓點(diǎn)高度a及起始角度θ0。

圖3 程序流程圖Fig.3 Program flow chart

2 氣球設(shè)計算例

2.1 工程算例

假定某次飛行任務(wù)中，設(shè)計指標(biāo)如表2所示。

表2 氣球設(shè)計要求

利用上述算法和流程得到部分參數(shù)計算結(jié)果，見表3。圖4展示了上升階段及最后穩(wěn)定階段氣球球形，圖5為制作氣球時聚乙烯膜裁制尺寸(裁膜曲線縱坐標(biāo)為氣球母線長度，橫坐標(biāo)為由分膜數(shù)確定的角度區(qū)間所劃出的氣球截面圓圓弧長度)。

表3 部分參數(shù)計算結(jié)果

圖4 氣球母線 圖5 裁膜曲線Fig.5 Cutting curve Fig.4 Balloon generatrix

以往工程上根據(jù)∑表[1]設(shè)計氣球，球膜質(zhì)量按經(jīng)驗系數(shù)確定。本文計算結(jié)果與∑表[1]設(shè)計及文獻(xiàn)[7]采用的Gill法計算結(jié)果相關(guān)參數(shù)對比見表4。

表4不同設(shè)計方法計算結(jié)果對比

Tab.4 Comparison of calculation results

結(jié)果體積/m3表面積/m2母線長度/m∑表203 76117 014115.57文獻(xiàn)[7]210 88117 410116.61本文213 45217 541116.78

將氣球附件質(zhì)量附加到載荷質(zhì)量上重新迭代計算得到的氣球體積比文獻(xiàn)[7]的大1.2%，較以往設(shè)計方法更為精確。

上升過程中氣球體積、球體高度、零壓點(diǎn)高度及單位體積浮力變化曲線如圖6所示。

圖6 上升過程球形相關(guān)參數(shù)對比Fig.6 Comparison of shape parameters in the rise process

由圖6可以得出，隨著氣球飛行高度升高，單位體積浮力降低，氣球體積膨脹，體積增大速率隨高度增加而增加，氣球球體高度降低，零壓點(diǎn)逐漸下移至底點(diǎn)。

氣球上升過程中各高度下氣球表面徑向應(yīng)力的計算結(jié)果如圖7所示。由圖7可以得出，球膜徑向應(yīng)力呈現(xiàn)中間小、兩邊大的分布趨勢，在兩極趨向于無窮大。這是由于氣球兩極截面圓半徑趨向于零的緣故。在氣球結(jié)構(gòu)處理上，常采用雙層或是多層頭部以減弱氣球頭部應(yīng)力集中[6]。文獻(xiàn)[8]所提出的“混合形”球形則是假設(shè)球膜周向應(yīng)力并不始終為零，而是在氣球的下半部(球底至最大半徑處)周向應(yīng)力為零，在氣球的上半部(最大半徑處至球頂)周向應(yīng)力從零單調(diào)增加，以至在球頂處周向應(yīng)力與徑向應(yīng)力相等[7]。據(jù)此約束條件求解得到的氣球球形能夠滿足在氣球頂部應(yīng)力趨于有限值[5]。

圖7 各高度下氣球表面徑向應(yīng)力分布Fig.7 Radial stress distribution on the surface of a balloon at various heights

2.2 氣球有限元分析及驗證

有限元方法也是一種為偏微分方程的邊界值問題尋找近似解的數(shù)值計算方法，能夠?qū)碚撚嬎愕玫降那蛐卧诘撞繎覓燧d荷、膜重力及壓力差梯度作用下的應(yīng)力應(yīng)變分布進(jìn)行計算。導(dǎo)入計算所得氣球模型，網(wǎng)格選擇三維減縮積分四節(jié)點(diǎn)M3D4R單元，采用掃掠算法生成結(jié)構(gòu)網(wǎng)格。按照氣球所在飛行高度計算(ρa(bǔ)-ρh)·g得到壓強(qiáng)差隨氣球高度變化的梯度系數(shù)。計算結(jié)果轉(zhuǎn)換在柱坐標(biāo)下查看。

圖8 徑向應(yīng)力分布云圖(36 km)Fig.8 Radial stress distribution cloud chart(36 km)

結(jié)合上述算例，選取典型狀態(tài)，分別對36 km處(完全膨脹)及30 km處(部分膨脹)氣球進(jìn)行有限元分析，如圖8～11所示。結(jié)果顯示，氣球能夠維持外形，應(yīng)力值均保持正值，氣球膜處于張緊狀態(tài)，應(yīng)力沿高度方向呈環(huán)帶分布(即同一高度處應(yīng)力值相同)，據(jù)此可判斷氣球沒有褶皺產(chǎn)生。載荷基本靠球膜徑向應(yīng)力承擔(dān)，并且符合應(yīng)力從中間向兩極逐漸增大的趨勢，周向不承受應(yīng)力(在氣球頂部局部區(qū)域出現(xiàn)周向應(yīng)力的原因是仿真中球膜會出現(xiàn)變形導(dǎo)致局部曲率改變，在理論推導(dǎo)中將球膜視作剛體)。因此，有限元分析結(jié)果能較好說明所設(shè)計氣球外形的正確性。

圖9 周向應(yīng)力分布云圖(36 km)Fig.9 Circumferential stress distributioncloud chart(36 km)

圖10 徑向應(yīng)力分布云圖(30 km)Fig.10 Radial stress distribution cloud chart(30 km)

圖11 周向應(yīng)力分布云圖(30 km)Fig.11 Circumferential stress distribution cloud chart(30 km)

3 關(guān)鍵參數(shù)敏感性分析

在球形設(shè)計的基礎(chǔ)上，選取在高空氣球設(shè)計時最重要的若干個關(guān)鍵參數(shù)，即飛行高度、載荷質(zhì)量及晝夜溫差，研究了這幾個量在變化時對氣球形狀的影響。

3.1 高度變化的影響分析

對氣球在攜帶相同質(zhì)量載荷(1000 kg) 的情況下，進(jìn)行了不同飛行高度(10 km,15 km,20 km,25 km,30 km,35 km,40 km)球形的計算，結(jié)果如圖12～13所示。

由圖可以得出，隨著高度的增加，氣球體積增加逐漸加快，而球體高度、截面圓半徑和母線長度增加相對平緩。并且截面圓半徑比球體高度增加稍慢，氣球形狀有向圓柱形發(fā)展的趨勢。

圖12 相同載荷不同高度球形對比Fig.12 Shape contrast at different height under the same load

圖13 相同載荷不同高度球形參數(shù)對比Fig.13 Comparison of shape parameters at different height under the same load

3.2 載荷變化的影響分析

對氣球在相同飛行高度(36 km)、攜帶不同質(zhì)量載荷(50 kg,100 kg,200 kg,400 kg,600 kg,800 kg,1000 kg)時的球形進(jìn)行計算，結(jié)果如圖14～15所示。

由圖可以得出，在同一飛行高度下，隨著載荷的增加，氣球體積幾乎保持線性增加，對比圖8的結(jié)果，說明高度(壓強(qiáng))是影響氣球體積大小的主要原因。球體高度、截面圓半徑和母線長度增加逐漸平緩，有線性增加的趨勢，這說明氣球載荷對同一飛行高度氣球的影響基本呈線性關(guān)系，符合以往的設(shè)計經(jīng)驗。

圖14 相同高度不同載荷球形對比Fig.14 Shape contrast of different load under the same height

圖15 相同高度不同載荷球形參數(shù)對比Fig.15 Comparison of shape parameters at different load under the same height

3.3 晝夜溫差的影響分析

另外，零壓氣球底部有排氣管與外界大氣相通，使得底部壓強(qiáng)差為零。在晝夜溫差影響下，球內(nèi)氦氣溫度發(fā)生變化，導(dǎo)致氣球體積變化而引起氣球飛行高度改變[12]。設(shè)計參數(shù)為飛行高度20 km，載荷質(zhì)量1000 kg，膜密度20 g/m2的零壓氣球在北緯30°，春分時期，高度20 km處24 h內(nèi)氦氣溫度變化引起的氣球飛行高度、體積、球體高度及零壓點(diǎn)高度變化曲線如圖16所示。

可以得出，氦氣晝夜溫度有50 K的變化，引起氣球體積在夜間縮小2500 m3，飛行高度偏離設(shè)計高度(20 km)達(dá)1.2 km，球體高度從30.5 m增高至32.75 m，零壓點(diǎn)高度則有近10 m的變化。因此需要進(jìn)行高度控制，工程上常采用投擲壓艙物以減輕負(fù)載質(zhì)量的手段[1]。超壓氣球則是利用超壓量大于溫度差引起的壓強(qiáng)變化，使氣球晝夜始終完全膨脹，能夠在穩(wěn)定高度飛行較長時間。

圖16 晝夜溫差引起零壓氣球參數(shù)變化Fig.16 Variation of zero pressure balloon parameters caused by temperature difference

4 總結(jié)

本文對零壓式高空氣球球形設(shè)計進(jìn)行了系統(tǒng)分析，包括理論推導(dǎo)，設(shè)計數(shù)值方法進(jìn)行求解，編寫整套氣球球體設(shè)計程序，以及對影響氣球形狀的關(guān)鍵參數(shù)進(jìn)行敏感性分析及應(yīng)用。首先借助薄膜無矩理論推導(dǎo)了氣球球形的控制方程組，選擇打靶法求解完全膨脹狀態(tài)球形，并將多段打靶法與等式二次規(guī)劃法相結(jié)合求解上升過程中的氣球外形。隨后，結(jié)合工程實際進(jìn)行了算例分析：一方面與以往工程計算及文獻(xiàn)結(jié)果進(jìn)行了對比；另一方面結(jié)合有限元分析方法對計算得到的球形進(jìn)行仿真，氣球外形變化、受力平衡及應(yīng)力分布情況表明球形符合零壓氣球應(yīng)力假設(shè)，驗證了設(shè)計方法的正確性。最后，對影響氣球球形的飛行高度、攜帶載荷以及晝夜溫差這些關(guān)鍵參數(shù)進(jìn)行了靈敏度分析，給出了定量的計算結(jié)果，獲得了有參考性的結(jié)論?？偟膩碚f，本文在之前學(xué)者的基礎(chǔ)上對其方法做了進(jìn)一步優(yōu)化與改進(jìn)，能夠為后續(xù)研究人員提供相應(yīng)的參考和借鑒。