陳 晨, 郭琪磊

(中國(guó)民用航空飛行學(xué)院，廣漢 618307)

降落傘的充氣展開過程歷時(shí)短、變形大，是典型的復(fù)雜非線性流固耦合問題，是降落傘研究的重點(diǎn)和難點(diǎn)。近年來，流固耦合ALE方法成為解決該問題的主流方式。

文獻(xiàn)[1-6]首先將LS-DYNA軟件應(yīng)用于降落傘FSI問題中，基于歐拉-拉格朗日罰函數(shù)算法，得到了初始形狀為傘衣平鋪展開的情況下，開傘過程中的阻力系數(shù)、擺動(dòng)頻率，擺動(dòng)角等參數(shù)，并觀察到降落傘的過度充氣和循環(huán)呼吸現(xiàn)象，證明了該方法具有很高的可靠性和實(shí)用前景。在此基礎(chǔ)上，文獻(xiàn)[7-10]改進(jìn)了歐拉-拉格朗日耦合算法，把傘衣的透氣性參數(shù)引入降落傘的開傘模擬計(jì)算中，取得較大成功。Cheng等[11]、高興龍等[12]和賈賀等[13-14]采用多孔介質(zhì)模型研究了無限質(zhì)量、不計(jì)傘重情況下，基于ALE方法的開傘充氣過程，得到傘衣投影面積和有效透氣性的變化規(guī)律曲線，該傘衣模型具有較真實(shí)的結(jié)構(gòu)細(xì)節(jié)，如傘衣幅的突起，得到的流場(chǎng)數(shù)據(jù)也更加符合實(shí)際情況。馬曉冬等[15]對(duì)渦環(huán)旋轉(zhuǎn)傘進(jìn)行了流固耦合分析與傘塔實(shí)驗(yàn)，得到渦環(huán)旋轉(zhuǎn)傘的充氣時(shí)間與開傘動(dòng)載。隨著研究的深入，發(fā)現(xiàn)基于有限質(zhì)量方式建立的開傘模型更為貼近空投實(shí)際情況，其減速性能與開傘力也更受關(guān)注，但此時(shí)流入傘衣的氣流速度和傘衣結(jié)構(gòu)都在短時(shí)間內(nèi)急劇變化，流場(chǎng)非穩(wěn)態(tài)、流固耦合模型與降落傘系統(tǒng)彈道方程相耦合，求解難度較大，相應(yīng)的研究結(jié)果較少。Tutt等[16]首先提出用有限質(zhì)量與動(dòng)網(wǎng)格方式建立拉直-展開下落模型，并通過實(shí)驗(yàn)對(duì)比，證明了該模型的正確性。高興龍等[17]基于任意拉格朗日-歐拉(arbitrary Lagrange Euler，ALE)方法和動(dòng)網(wǎng)格技術(shù)，建立了平面縱向開縫傘低速空投的有限質(zhì)量開傘模型，得到傘衣阻力面積變化曲線，并驗(yàn)證了開傘動(dòng)載的變化規(guī)律，最終分析了不同拉直速度對(duì)開傘時(shí)間的影響規(guī)律。Cheng等[18]、Yu等[19]基于有限元方法計(jì)算了有限質(zhì)量情況下小速度開傘動(dòng)載變化特征，并分析了傘衣上受載嚴(yán)重位置。

前人研究為數(shù)值計(jì)算降落傘空投流構(gòu)耦合(fluid structure interaction，F(xiàn)SI)問題奠定了堅(jiān)實(shí)的基礎(chǔ)，并積累了豐富的經(jīng)驗(yàn)，但以上計(jì)算中大部分存在一些模型假設(shè)的不合理性，如未考慮傘衣織物孔隙率、透氣性影響，用無限質(zhì)量模型模擬空投工況等。更為普遍的是：在降落傘拉直-展開建模計(jì)算中，大部分文獻(xiàn)假設(shè)傘衣初始為拉直狀態(tài)，并設(shè)拉直速度為初始速度。但實(shí)際工況中，初始投放瞬間，傘繩上并未加載初始應(yīng)力，初始不處于張緊狀態(tài)(即尚未拉直)，這與實(shí)際工況不符，導(dǎo)致模型計(jì)算的拉直速度、下落距離和拉直/開傘時(shí)間有一定誤差。

為完善數(shù)值計(jì)算模型，擬基于ALE方法和動(dòng)網(wǎng)格技術(shù)，考慮織物透氣性并預(yù)設(shè)初始投放速度，建立降落傘充氣展開的有限元模型。以C-9型傘作為算例，計(jì)算充氣過程中，有限質(zhì)量情況下的開傘力，并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，擬通過實(shí)驗(yàn)驗(yàn)證該模型的正確性與合理性，隨后分析其開傘過程中應(yīng)力分布特征與氣動(dòng)減速特性。預(yù)期為降落傘數(shù)值分析和設(shè)計(jì)提供參考依據(jù)。

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

1.1 模型假設(shè)

(1)降落傘傘衣初始為折疊狀態(tài)，傘繩為近似拉直狀態(tài)。

(2)忽略傘衣織物摩擦力的影響。

(3)由于降落傘投放區(qū)域海拔不高，且運(yùn)動(dòng)速度小于0.3馬赫數(shù)，可忽略空氣的壓縮性，用時(shí)變非穩(wěn)態(tài)不可壓流模擬流場(chǎng)特性。

1.2 材料本構(gòu)模型

1.2.1 傘衣模型

柔性織物在展開過程中可看作線彈性各向同性材料，則根據(jù)彈性力學(xué)理論，其本構(gòu)方程為

σij=λεkkδij+2?εij

(1)

式(1)中：σ為應(yīng)力張量，σij為σ的分量；ε為Green應(yīng)變張量，εkk、εij為ε的分量；δij為Kronecker記號(hào)；λ和?為拉梅常數(shù)，可以由材料的楊氏模量(E)和泊松比(υ)計(jì)算，其計(jì)算公式為

(2)

(3)

用三節(jié)點(diǎn)殼單元模擬傘衣織物，則根據(jù)平面應(yīng)力假設(shè)和材料對(duì)稱性假設(shè)，式(1)可簡(jiǎn)化為

(4)

1.2.2 空氣模型

理想氣體狀態(tài)方程結(jié)合流體本構(gòu)方程模型可建立流體域計(jì)算模型。其中理想氣體狀態(tài)方程為

p=(1-γ)ρωint

(5)

式(5)中：p為靜水壓力；γ為氣體比熱；ρ為氣體密度；ωint為氣體單位質(zhì)量?jī)?nèi)能。

(6)

式(6)中：v為運(yùn)動(dòng)速度，vi,j、vj,i為v分量；μ為動(dòng)力黏性系數(shù)。

1.2.3 傘繩模型

用一維線彈性本構(gòu)方程計(jì)算傘繩拉力：

F=Kmax(Δl,0)

(7)

式(7)中：Δl為繩長(zhǎng)變化量；K為材料剛度。

Δl=l-(l0-loff)

(8)

(9)

式中：l為傘繩當(dāng)前長(zhǎng)度；l0為傘繩原長(zhǎng)；E為繩索彈性模量；A為繩索截面積；loff為偏移量，模型中傘繩初始無預(yù)應(yīng)力，即為松弛狀態(tài)，故設(shè)偏移量為負(fù)值。

1.3 控制方程

用ALE方程求解自由界面流動(dòng)及典型的流固耦合問題。具有物質(zhì)運(yùn)輸功能的單元在兩種坐標(biāo)體系的網(wǎng)格間平動(dòng)，在耦合界面處流體節(jié)點(diǎn)隨結(jié)構(gòu)變形運(yùn)動(dòng)。將流體質(zhì)點(diǎn)速度通過質(zhì)量守恒和動(dòng)量守恒定律在黏性不可壓流體域內(nèi)離散。

1.3.1 結(jié)構(gòu)動(dòng)力學(xué)模型

傘衣織物是柔性大變形體，具有強(qiáng)非線性動(dòng)力學(xué)特征。令Ωs為空間結(jié)構(gòu)域，固體邊界用?Ωs表示，可得結(jié)構(gòu)控制方程，在固體邊界?Ωs上[20]有

(10)

式(10)中:y為位移矢量；ρs為結(jié)構(gòu)材料密度；f為作用在結(jié)構(gòu)上的體積力；σs為結(jié)構(gòu)體Cauchy應(yīng)力張量;t為時(shí)間。

1.3.2 流場(chǎng)動(dòng)力學(xué)模型

ALE方程描述中，在參考坐標(biāo)系下，不可壓流體的納維-斯托克斯方程(Navier-Stokes,N-S)方程為

(11)

(12)

(13)

式中：u為流體速度；w為參考構(gòu)型下網(wǎng)格點(diǎn)的移動(dòng)速度；ρf為流體密度；e為物質(zhì)內(nèi)能。對(duì)于牛頓流體，應(yīng)力張量(σf)定義為

σf=-pI+μ[gradu+(gradu)T]

(14)

式(14)中：p為壓力；I為二階單位張量；μ為動(dòng)力黏度系數(shù)。

其中，Dirichlet和Neumann邊界條件為

(15)

(16)

u(x,0)=u0

(17)

式(17)中:u0為初始來流速度。

1.4 耦合方式-罰函數(shù)法

傘衣與流體域的流構(gòu)耦合要求在建立有限元模型時(shí)，結(jié)構(gòu)與流體的網(wǎng)格可以重疊在一起，但又能實(shí)現(xiàn)力學(xué)參量的傳遞，使用罰函數(shù)法。追蹤在任意tn時(shí)間步內(nèi)，拉格朗日節(jié)點(diǎn)(結(jié)構(gòu)節(jié)點(diǎn)，從節(jié)點(diǎn))和歐拉流體物質(zhì)(流體節(jié)點(diǎn)，主節(jié)點(diǎn))位置間的相對(duì)位移dn，主從節(jié)點(diǎn)的相對(duì)速度為vrel，則dn的迭代更新方程為

(18)

耦合界面上的力與發(fā)生的貫穿數(shù)量成正比：

Fn=kidn

(19)

式(19)中：ki表示基于主、從節(jié)點(diǎn)質(zhì)量模型特征的剛度系數(shù)。

1.5 透氣性

傘衣為柔性透氣織物，具有滲透性，其耦合力可通過薄膜狀Shell單元體積的Ergun方程導(dǎo)出得[12]

(20)

式(20)中：r為傘衣單元的法向；a(μ,ε)為滲透性Shell單元的滲透系數(shù)(黏性系數(shù))；b(ρ,ε)為慣性系數(shù)；a、b組合反映了織物的透氣性特征[14]。

1.6 初始投放速度

在降落傘初始空投時(shí)，初始投放速度并不等于拉直速度。初始投放時(shí)，傘衣和重物速度相等，重物受的空氣阻力小，運(yùn)動(dòng)速度較快，到達(dá)某個(gè)位置后具有相對(duì)傘衣的最大速度差Δvmax，在此位置，傘繩被繃緊并開始拉伸。當(dāng)傘繩拉直時(shí)，傘衣被重新加速，以致在此瞬間傘衣與重物的速度差Δv=0，而傘繩的黏彈性伸長(zhǎng)度已吸收了減速器的全部動(dòng)能后，軸向張力載荷出現(xiàn)瞬時(shí)峰值，即拉直力，而此時(shí)傘衣和重物的速度才是拉直速度。

根據(jù)Wolf提出的最大速度差估算方法預(yù)估投放初始速度[21]：

(21)

(22)

(23)

式中：Δvmax為傘衣相對(duì)回收體的最大速度差；v0為初始投放速度；f為擬合曲線斜率；ρ為空氣密度；(CDS)p為傘衣的阻力面積；l1為傘繩初始長(zhǎng)度；mp為減速器質(zhì)量；mc為傘衣質(zhì)量；m1為傘繩質(zhì)量；Kb為系數(shù)，表示不同工況，可從文獻(xiàn)[21]中讀??；CDA為重物阻力面積；mb為重物質(zhì)量；vb為傘繩拉伸時(shí)重物速度。

根據(jù)式(21)～式(23)預(yù)估的初始投放速度可能由于傘型、有無吊帶、降落傘質(zhì)量分布等因素影響與實(shí)際情況稍有偏差，需通過試算調(diào)整初始投放速度，獲取真實(shí)值。

2 計(jì)算模型

以C-9型傘作為算例，建立全尺寸數(shù)值模型，計(jì)算其空投自由下落→充氣→開傘(折疊展開)→穩(wěn)定下落的整個(gè)過程，重點(diǎn)關(guān)注充氣展開階段。模型參數(shù)如表1所示。

表1 降落傘模型參數(shù)

降落傘下懸掛重物(載荷)，以90°的軌跡傾斜角空投，具體工況如表2所示。

表2 空投實(shí)驗(yàn)工況

由于不考慮懸掛載荷的氣動(dòng)特性，為簡(jiǎn)化計(jì)算，用一個(gè)長(zhǎng)方體板殼結(jié)構(gòu)模擬重物。帶有傘頂孔的降落傘折疊模型網(wǎng)格(圖1)，其中包括傘衣、傘繩和重物，傘繩處于未完全拉直狀態(tài)，無初始預(yù)應(yīng)力。為模擬有限質(zhì)量情況下的空投過程，用動(dòng)網(wǎng)格技術(shù)追蹤降落傘，所建立的流體網(wǎng)格域如圖1所示。傘衣位于流體域的中下部，在傘衣附近加密流場(chǎng)網(wǎng)格，且等比放大填充整個(gè)流體域，以優(yōu)化計(jì)算效率。考慮傘衣織物透氣性，基于顯式有限元方法，采用二階精度的Van Leer MUSCL格式求解控制方程。

圖1 降落傘及流場(chǎng)域模型

3 結(jié)果驗(yàn)證

計(jì)算結(jié)果與空投實(shí)驗(yàn)結(jié)果[21]對(duì)比如圖2所示。圖2中，F(xiàn)為開傘力，Wb為重物重力，t為時(shí)間，tf為傘衣充氣時(shí)間。從圖2可知，計(jì)算的開傘力與實(shí)驗(yàn)結(jié)果吻合較好：開傘時(shí)間基本一致；開傘力變化趨勢(shì)基本相同；開傘力有兩個(gè)峰值，且其出現(xiàn)時(shí)機(jī)基本一致。但計(jì)算值的開傘過載整體偏大，且比實(shí)驗(yàn)稍有提前，這主要是由于傘衣結(jié)構(gòu)模型未能充分考慮阻尼耗散與摩擦力的影響造成的。兩條曲線中比較明顯的區(qū)別是，計(jì)算曲線在t/tf=0.25出現(xiàn)明顯沖擊載荷，而實(shí)驗(yàn)值沒有，這正是數(shù)值模擬計(jì)算中體現(xiàn)了拉直力的作用。在設(shè)計(jì)和實(shí)際投放降落傘時(shí)，都已采用了引導(dǎo)傘或開傘袋來抵消拉直時(shí)的猛烈沖擊，而本模型計(jì)算中并未考慮此項(xiàng)，而是采用直接拉直的方式，從而反映未作任何設(shè)計(jì)時(shí)，拉直力的大小與拉直時(shí)間。

圖2 開傘力隨時(shí)間的變化曲線

參考文獻(xiàn)[21]中的另一算例中，給出了拉直力和開傘力間的關(guān)系，由于實(shí)驗(yàn)工況未明，只有拉直力與開傘力的峰值與變化趨勢(shì)可參考。將計(jì)算結(jié)果中拉直瞬間作為時(shí)間零點(diǎn)繪制圖3。

圖3 拉直力與開傘力間的關(guān)系

Fs-test表示實(shí)驗(yàn)中拉直力峰值，F(xiàn)1-test表示開傘力峰值，F(xiàn)s-cal表示計(jì)算拉直力峰值，F(xiàn)1-cal表示開傘力峰值。從圖4可得，實(shí)驗(yàn)結(jié)果中，F(xiàn)s-test/F1-test=0.83，而計(jì)算結(jié)果中，F(xiàn)s-cal/F1-cal=0.8，相對(duì)誤差為3.6%，屬于工程可接受誤差范圍，這也證明了本文模型的正確性。

圖4 折疊開傘過程中的應(yīng)力分布

4 結(jié)果分析

4.1 折疊開傘過程

隨著時(shí)間的變化，開傘過程與應(yīng)力分布如圖4所示。其中零時(shí)刻為投放時(shí)刻。

0時(shí)刻，投放傘體，傘衣折疊，與重物一起自由落體下降。0.015 s傘衣底口張開，空氣質(zhì)量開始以氣團(tuán)的形式進(jìn)入傘衣內(nèi)部，并沿傘衣運(yùn)動(dòng)。0.02 s時(shí)少量氣團(tuán)首次到達(dá)傘衣頂部。0.045～0.075 s時(shí)，氣團(tuán)繼續(xù)匯聚，擴(kuò)大傘衣內(nèi)部，傘衣因徑向動(dòng)量而膨脹，而傘衣底邊面積基本保持不變。最大應(yīng)力出現(xiàn)在傘繩連接和傘頂處，前者是受傘繩拖曳力的影響，后者是受氣團(tuán)附加質(zhì)量力的影響。0.075 s時(shí)，更多氣團(tuán)到達(dá)傘頂，一小部分從傘頂孔逸出，傘衣從底邊充滿到傘頂，大部分形成空氣附加質(zhì)量，拉伸傘繩，傘繩受力伸長(zhǎng)，傘衣受傘繩、重物和空氣附加質(zhì)量的影響，形成瞬間拉直力的峰值，此時(shí)拉直速度為76.2 m/s，最大應(yīng)力在傘衣與傘繩連接處。0.095 s時(shí)，更多氣流匯聚到傘頂孔處，傘頂膨脹，應(yīng)力大多集中在頂孔附近。隨后，在0.095～0.185 s，空氣質(zhì)量開始沿傘衣流管充入傘衣內(nèi)部，除傘頂外的傘衣也逐漸膨脹，此時(shí)，傘衣周向受力基本均勻分布，整個(gè)降落傘的下落速度也趨于穩(wěn)定，最大應(yīng)力仍集中在頂孔和傘繩連接處。0.185 s時(shí)，大量空氣質(zhì)量團(tuán)首次匯聚到傘頂孔處，形成“烏賊狀”傘衣，代表著充氣過程中最初階段的終結(jié)，也形成了開傘作用力的第一個(gè)峰值，頂孔褶皺處有少許應(yīng)力集中。0.185 s后，降落傘減速效果明顯，傘衣徑向膨脹速度加快，頂孔和傘衣上部應(yīng)力不斷增大，受氣流與織物耦合力的影響出現(xiàn)應(yīng)力不對(duì)稱、形狀不對(duì)稱等開傘過程(0.285 s)，但隨著時(shí)間推移，受氣動(dòng)力影響，這種不對(duì)稱被減弱乃至消除(0.335 s)，傘衣趨向于均勻膨脹，頂孔處傘衣完全展開，應(yīng)力逐漸減小，傘衣上部應(yīng)力逐漸增大，應(yīng)力集中部位也由頂孔處轉(zhuǎn)移至傘衣上部褶皺處。0.4 s時(shí)，傘衣投影直徑(Dp)/傘衣名義直徑(D0)=0.58，接近完全張滿，此時(shí)，出現(xiàn)開傘力的第2個(gè)峰值。隨著傘衣膨脹，傘衣中部應(yīng)力逐漸增大而頂部減小，柔性織物不規(guī)則變形，褶皺處出現(xiàn)最大應(yīng)力(參考0.41 s時(shí))。在0.435～0.475 s時(shí)，Dp/D0為0.67～0.7傘衣充滿[21]，傘衣周向應(yīng)力趨于平均，僅褶皺處存在應(yīng)力集中，隨后隨著傘衣膨脹逐漸消失。之后，傘衣繼續(xù)穩(wěn)定下降，空氣附加質(zhì)量趨于恒定，出現(xiàn)輕微的呼吸現(xiàn)象，傘衣提供的穩(wěn)態(tài)氣動(dòng)阻力使重物減速下落(圖5)，圖5中，橫軸為時(shí)間t(s)，縱軸為傘衣和重物的下落速度V(m/s)。隨后，周圍空氣質(zhì)量的動(dòng)量造成傘衣底邊過度擴(kuò)張膨脹和傘頂塌陷(0.55 s后)。

圖5 傘衣及重物運(yùn)動(dòng)速度變化曲線

4.2 傘衣周圍流場(chǎng)

開傘過程中，值得注意的是5個(gè)典型階段，即拉直→充氣→烏賊狀→膨脹→充滿，此時(shí)傘衣周圍流場(chǎng)變化如圖6所示。

圖6 傘衣周圍流場(chǎng)變化

0.045 s時(shí)，氣流由傘底流向傘頂，并從頂孔加速逸出，在頂孔處形成扇形低壓區(qū)，誘導(dǎo)更多氣流從傘底進(jìn)入傘衣流管內(nèi)部，傘底處形成入流渦。0.075 s時(shí)，傘繩拉直，更多氣團(tuán)到達(dá)傘衣頂部，并在頂孔周圍匯聚、加速，形成傘衣內(nèi)外較大壓差，頂孔處的氣流穩(wěn)定流出，形成尾跡，而此前形成的扇形低壓區(qū)逐漸減弱消失。有少量氣體從傘衣織物表面逸出，體現(xiàn)了傘衣的透氣性特征。0.125 s時(shí)，傘衣內(nèi)部氣流加速，形成徑向動(dòng)量，膨脹整個(gè)傘衣，傘衣周圍氣流在結(jié)構(gòu)的作用下形成數(shù)個(gè)渦，引起傘衣結(jié)構(gòu)瞬間非線性大變形，而結(jié)構(gòu)同時(shí)影響渦的位置的大小。相對(duì)的，頂孔處的尾流更加流暢和穩(wěn)定。0.185 s時(shí)，傘衣呈“烏賊狀”，結(jié)束充氣過程的最初階段，傘頂?shù)臏u在傘衣外形和壓力梯度變化的作用下脫落，形成非對(duì)稱尾跡，這與文獻(xiàn)[20]結(jié)論相同，也符合鈍體擾流尾跡特征。0.285 s時(shí)，傘衣在氣動(dòng)力和慣性力的作用下繼續(xù)膨脹并左右輕微搖晃，傘衣周圍氣流在傘衣結(jié)構(gòu)作用下被加速，產(chǎn)生附加質(zhì)量力，尾渦加速形成和脫落。0.4 s時(shí)，傘衣周圍氣流擾動(dòng)加強(qiáng)，阻力面積增大，開傘力達(dá)到峰值，尾跡也形成規(guī)律的周期性脫落渦。0.455 s，傘衣充滿，傘衣內(nèi)部徑向動(dòng)量趨于消失，周圍氣流也趨于穩(wěn)定，形成有效阻流體，尾跡渦繼續(xù)周期性形成并脫落，并且更加穩(wěn)定和強(qiáng)烈。0.55 s后，傘衣過度擴(kuò)張，傘頂附近的渦破碎分離，氣流動(dòng)量引起傘頂塌陷，但整個(gè)系統(tǒng)依舊保有減速回收效果。

5 結(jié)論

基于ALE方法，考慮織物透氣性并預(yù)設(shè)初始投放速度，建立了有限質(zhì)量情況下，降落傘折疊展開的有限元模型。以C-9傘為算例，計(jì)算其拉直-充氣-展開過程，并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了該模型的正確性，得到如下結(jié)論。

(1)計(jì)入織物透氣性并預(yù)設(shè)初始投放速度后，該模型能較為準(zhǔn)確地計(jì)算折疊-開傘過程，其開傘時(shí)間與開傘應(yīng)力都滿足工程精度要求，且能體現(xiàn)傘繩上的拉直力。

(2)傘衣充氣初期，氣流進(jìn)入傘衣這個(gè)細(xì)長(zhǎng)的流管，受氣流和織物耦合影響，其周向應(yīng)力分布并不均勻，容易造成傘頂變形、傘頂甩打等異?，F(xiàn)象而導(dǎo)致開傘失敗。

(3)開傘過程中，傘衣周向應(yīng)力分布逐漸趨于均勻，應(yīng)力較大處也由頂孔周圍轉(zhuǎn)為傘衣上部，再到傘衣中部。不同開傘時(shí)刻時(shí)的應(yīng)力危險(xiǎn)截面各有不同。

(4)開傘初期，頂孔附近形成低壓區(qū)，傘衣逐漸膨脹，其周圍氣流不斷加速。隨著開傘過程，傘頂處形成附著渦，隨后周期性脫落。