空投貨臺擺蕩階段流體-固體耦合建模分析

湯健華,錢林方,徐亞棟

(南京理工大學機械工程學院,江蘇南京210094)

空投貨臺擺蕩階段流體-固體耦合建模分析

湯健華,錢林方,徐亞棟

(南京理工大學機械工程學院,江蘇南京210094)

為獲得空投貨臺擺蕩階段的精確動力學響應(yīng),考慮貨臺大幅度擺蕩對空氣產(chǎn)生的非定常效應(yīng),基于平衡點法構(gòu)建了空投吊掛系統(tǒng)多體動力學模型,并采用基于ALE格式的有限體積法建立了空氣流場模型,在此基礎(chǔ)上推導(dǎo)建立了非定常流動作用下的貨臺-空氣流體-固體耦合分析模型。通過數(shù)值算例研究:驗證了該多體動力學模型的魯棒性和正確性;明確了層流非定常模型能較好地描述流體對平板的作用力,可用于空投貨臺擺蕩階段的流場分析;對比指出通過流體-固體耦合分析得到的非定常流動氣動力顯示出震蕩特性,而定常流動假設(shè)模型會給貨臺帶來類似粘滯阻尼的特性,其減速效果要比非定常流動模型強烈。

流體力學;非定常流動;平衡點法;空投貨臺;流體-固體耦合

DOI:10.3969/j.issn.1000-1093.2016.01.021

0 引言

當空投貨臺從機艙拉出以后,空投貨臺在自身重力、初速度作用下進行落體運動,同時降落傘隨著貨臺及牽引傘的作用拉出傘包并開始充氣,空氣對降落傘的氣動阻力通過吊掛連接繩索傳遞到貨臺上,使貨臺瞬時減速,并在空中進行大范圍的擺蕩,其結(jié)構(gòu)示意圖如圖1所示。在周圍空氣的作用下,貨臺動能以加速貨臺周圍空氣運動的形式耗散,最終貨臺平穩(wěn)下降。這個過程稱為空投系統(tǒng)的擺蕩階段,在此期間伴隨有強烈的空氣流動,貨臺經(jīng)歷大幅度的位移、翻轉(zhuǎn)和速度突變,可能對空投物資產(chǎn)生損害。因此,建立精確的數(shù)學模型預(yù)測貨臺在擺蕩階段的響應(yīng),對于空投系統(tǒng)的改進、優(yōu)化具有實際的工程意義。

圖1 空投貨臺及吊掛結(jié)構(gòu)示意圖Fig.1 Schematic diagram of Airdrop sling system

研究貨臺在擺蕩階段的響應(yīng),關(guān)鍵問題在于構(gòu)建精確的空氣對貨臺的氣動力模型[1-2]。早期研究工作主要通過試驗或數(shù)值的方法獲得貨臺在定常流動下的氣動力系數(shù)。Kenneth[3]通過試驗的方法測定了定常流動下鈍體平板的升阻力系數(shù)與迎角的關(guān)系,并羅列出了影響流場及升力的因素[4]。Mcquilling等[4]使用計算流體力學方法,計算了定常流動下貨臺在不同迎角所受的升力和阻力,并與試驗結(jié)果對比吻合較好。然而,基于定常流動的氣動力模型與貨臺擺蕩階段的實際情況有較大的差別。文獻[2]進行空投貨臺試驗發(fā)現(xiàn),采用定常流動的氣動力模型作為驅(qū)動,計算得到的貨臺響應(yīng)與試驗有明顯區(qū)別。表現(xiàn)為,當貨臺仍然做較大范圍的擺蕩運動時,仿真計算的結(jié)果已經(jīng)趨于穩(wěn)定。Kenneth[3,5]通過試驗研究了剛性平板在液體中的升阻力,指出定常流動下得到的升阻力系數(shù)與非定常流動下得到的測量結(jié)果相差較大。貨臺在擺蕩過程中迫使周圍的空氣運動改變空氣流場分布,呈現(xiàn)出明顯的非定常流動特性,因此應(yīng)基于非定常流動理論建立氣動力模型。

貨臺因為擺蕩運動改變了周圍流場的分布,同時也受到氣動力的驅(qū)動作用,運動姿態(tài)發(fā)生調(diào)整,二者之間相互耦合作用。為了獲得貨臺擺蕩階段精確的動態(tài)響應(yīng),必須要建立多體動力學和流場計算相耦合的分析模型。然而由于問題的復(fù)雜性,關(guān)于貨臺的流固耦合分析模型,尤其是非定常流體-貨臺耦合模型的研究鮮見報道。

本文基于平衡點法和ALE格式有限體積法推導(dǎo)建立了考慮非定常流動作用的貨臺-空氣流體-固體耦合分析模型。首先,推導(dǎo)基于平衡點法的Newton-Raphson統(tǒng)一迭代格式,得到能夠精確捕捉吊掛松弛-張緊的貨臺多體動力學計算模型。然后,對時間域進行離散,在每一個時間步內(nèi)通過建立的多體模型獲得貨臺位置響應(yīng),更新貨臺流體位置界面。通過基于ALE格式有限體積法求解經(jīng)過界面更新的流場,得到空氣流體作用力。將獲得的氣動力參數(shù)作為下一時間步多體模型的輸入載荷。最后,通過算例驗證說明了考慮非定常流動作用的耦合模型與采用定常流動氣動力模型的區(qū)別。

1 貨臺多體動力學分析

1.1 貨臺動力學問題描述

吊掛結(jié)構(gòu)是空投系統(tǒng)最常用的一種結(jié)構(gòu)。典型的吊掛結(jié)構(gòu)通過數(shù)根連接繩索把貨臺與承載結(jié)構(gòu)(如降落傘等)連接起來,這種結(jié)構(gòu)在工程中已有廣泛的應(yīng)用[6-7]。在擺蕩過程中,貨臺由于自身慣性,其頭部或尾部翹起。使吊掛結(jié)構(gòu)中的吊帶經(jīng)歷松弛-繃緊-松弛的交替狀態(tài),如圖2所示。

圖2 貨臺在空中擺蕩時繩索的松弛Fig.2 Slackness of sling lines

當繩索兩連接點間的距離小于繩索的原始長度時,繩索剛度為0.相反地,當繩索兩連接點間的距離大于繩索的原始長度時,繩索近似服從胡克定律,其材料特性如圖3所示。張力-應(yīng)變可以通過分段函數(shù)表示[8-9],如 (1)式所示。

式中:l是加載前未變形的原始長度;p是繩索變形后的長度;s是繩索剛度的斜率,很明顯,繩索的剛度斜率不滿足1階連續(xù)性。

圖3 繩索的材料特性Fig.3 Constitutive relations of sling lines

由于繩索只能傳遞力而不能承受彎矩,因此,把連接節(jié)點考慮為質(zhì)點,不考慮連接節(jié)點的轉(zhuǎn)動對繩索的影響。通過這一假設(shè),多體動力學運動方程可以描述為

式中:M是質(zhì)量矩陣;K是剛度矩陣;U是繩索節(jié)點的位置向量;x是廣義位移向量。由于繩索是典型的非線性材料,且剛度斜率不滿足1階連續(xù)性,使得吊掛系統(tǒng)的微分方程組具有較大的剛性,迭代求解時有可能在根附近來回震蕩。采用平衡點法可以有效的解決這個問題,使吊掛系統(tǒng)的多體動力學模型更加魯棒穩(wěn)定。

1.2 平衡點法統(tǒng)一迭代格式

平衡點法最早由Poole提出并用于求解剛性方程組[10],經(jīng)過文獻[11-13]的發(fā)展,已經(jīng)證實能夠應(yīng)用于多根繩索的模型以及吊掛系統(tǒng)計算中。平衡點法就是通過求解連接節(jié)點的力平衡方程組得到連接節(jié)點的位移,本文采用Newton-Raphson方法推導(dǎo)建立了基于平衡點法的統(tǒng)一迭代格式。

忽略連接節(jié)點自由度,多體動力學運動方程變?yōu)?/p>

1)固定連接節(jié)點的位置,求解其他部件的的響應(yīng)。

2)固定其他的部件,反復(fù)迭代求出連接節(jié)點的位置,直到收斂。

圖4 尋找連接點坐標的迭代過程Fig.4 Iteration process of connective nodes

不失一般性,如(3)式所示,可以得到在連接點處的合力為

式中:Fi是第i根繩索的阻力。

為了能夠得到數(shù)值求解連接點的坐標信息,這里使用牛頓方法進行迭代,迭代的過程表示為

式中:變化的位移為

Jk是第k個Newton迭代步的雅克比矩陣。對于三維情況而言,矩陣求解規(guī)模為3階。

式中:Pi是連接點到繩索另一端的向量;pi是該向量的模;ki是繩索的彈性模量;li是沒有加載變形前繩索的原始長度。對該式的彈性力向量求微分得到

由關(guān)系式

可以得到如下關(guān)系:

第i根繩索的雅克比矩陣可以修改為

式中:I為單位矩陣。吊掛關(guān)于連接節(jié)點的雅克比矩陣為

考慮到繩索的剛度可能為0,在迭代的過程中,有可能只有兩個或單個繩索受力,導(dǎo)致奇異,如圖5所示。對雅克比矩陣J的行列式進行判斷,當行列式小于容差時,需要對雅克比矩陣J進行滿秩分解:

式中:r為張緊繩索的條數(shù)。使用廣義逆公式,得到雅克比矩陣的廣義逆來替代雅克比矩陣的逆:

圖5 可能導(dǎo)致奇異的情形Fig.5 Singularity phenomenon

因為連接點的維數(shù)限制,迭代矩陣的規(guī)模為3階,大大降低了求解的難度,提高求解效率。而且對時間步長也沒有了穩(wěn)定性限制。

2 非定?？胀敦浥_流體-固體耦合分析

為了考察貨臺大幅度擺蕩對空氣產(chǎn)生的非定常效應(yīng),需要考慮隨時間不斷變化的貨臺對周圍空氣流場分布的影響。本文通過基于ALE格式的有限體積法對貨臺周圍空氣流場分布進行求解,流體質(zhì)量方程和動量方程為

應(yīng)當注意的是,基于定常假設(shè)的流體-固體耦合模型通過貨臺角度與氣動力系數(shù)計算作用在貨臺上的氣動力。因此需要額外獲取定常氣動力系數(shù)以及力矩系數(shù)隨角度變化關(guān)系。而采用上述流體-固體耦合分析方法由于已經(jīng)通過求解NS方程獲得作用在貨臺上的升阻力及作用在貨臺上的力矩,不需要額外獲取氣動力系數(shù)。

圖6 流體-固體耦合分析過程Fig.6 Flow chart of FSI Process

3 算例

3.1 平衡點法分析

為了驗證所建立動力學模型的正確性,本文計算了貨臺在真空中的響應(yīng)。將降落傘的連接點固定,建立的模型如圖7所示,將計算結(jié)果與商業(yè)多體動力學程序RECURDYN結(jié)果進行對比。

圖7 貨臺的初始狀態(tài)Fig.7 Initial state of platform

驗證的模型參數(shù)為:平板的尺寸為4.864 1 m× 2.0574 m×0.101 6 m,繩索連接于平板的4個角,4根連接繩索的原始長度為l1=4.750 71 m.初始速度幅值為5 m/s,方向沿著χ軸正方向,連接節(jié)點與固定點繩索的原始長度 l2=2 m.慣性矩陣為104IKg·m2,繩索的本構(gòu)關(guān)系為

通過控制較小的積分步長得到精確的解,RECURDYN程序的最大時間步長取10-6s,平衡點法的時間步長取為10-3s.分別計算了連接點位移,平板速度與角速度,如圖8～圖10所示,對比顯示兩者計算結(jié)果一致。說明基于平衡點方法建立的多體動力學模型能夠在較大的時間步長下準確地描述貨臺吊掛結(jié)構(gòu)的響應(yīng),具有良好的魯棒性。

圖8 連接節(jié)點χ位置隨時間變化Fig.8 Positon of connective node vs.time in direction χ

圖9 χ方向速度隨時間變化曲線Fig.9 Velocity of platform vs.time in direction χ

圖10 角速度隨時間變化曲線Fig.10 Angular velocity vs.time

3.2 非定常流動模型對氣動力的影響

Kenneth[5]曾經(jīng)通過試驗手段(見圖 11和圖12)獲得定常流動與非定常流動的升阻力系數(shù),并指出定常流動下得到的升阻力系數(shù)與非定常流動下得到的測量結(jié)果相差較大。試驗的結(jié)構(gòu)如圖11所示。

圖11 平板的主視圖與俯視圖Fig.11 Front and vertical views of plate

圖12 平板角度α隨時間變化曲線Fig.12 Angle α vs.time

在水槽中的平板與水槽來流方向成角度α,使用的流體是水,來流速度為0.2 m/s,雷諾數(shù)為3.8×104.試驗的具體尺寸參數(shù)如表1所示。

表1 定常模型尺寸參數(shù)Tab.1 Sizes of water tunnel and plate

Kenneth[5]通過固定平板攻角,待來流趨于穩(wěn)定情況下測量流體對平板的作用力,獲得定常流動下平板的升阻力與角度α的關(guān)系,通過(19)式進一步得到升阻力系數(shù)與角度α的關(guān)系。

式中:C為升阻力系數(shù)向量;F為平板所受的力向量;v∞為來流速度;Sp為貨臺的面積。

本文通過求解定常N-S方程得到升阻力系數(shù)隨角度的變化關(guān)系,并與試驗測量數(shù)據(jù)進行對比,如圖13、圖14所示。從對比結(jié)果可以看到,雖然最大阻力系數(shù)與試驗結(jié)果有一定誤差,但是變化趨勢與試驗結(jié)果吻合較好。

Kenneth[5]通過驅(qū)動平板按預(yù)定規(guī)律運動,攻角α隨時間的變化關(guān)系如圖12所示,實時測量流場對平板的作用力,獲得非定常流動下的升阻力系數(shù)(如圖15～圖18中的實驗數(shù)據(jù))。在已知平板運動規(guī)律(圖12所示)和升阻力系數(shù)與攻角關(guān)系(圖13和圖14)的前提下,確定對應(yīng)時刻平板的攻角值,并對該攻角值所對應(yīng)的升阻力系數(shù)進行插值,就可以近似得到定常流動假設(shè)下升阻力系數(shù)的時間歷程。圖15和圖16對比了基于定常流動假設(shè)計算的升阻力系數(shù)與試驗測量結(jié)果,由于通過定常流動計算得到的升阻力系數(shù)并沒有充分考慮平板運動對流場的影響,因此與試驗結(jié)果有較大差別。

圖13 定常阻力系數(shù)隨角度變化曲線Fig.13 Steady drag coefficient vs.α

圖14 定常升力系數(shù)隨角度變化曲線Fig.14 Steady lift coefficient vs.α

圖15 插值得到的阻力系數(shù)時間歷程Fig.15 Time history of interpolated steady flow drag coefficient

圖16 插值得到的升力系數(shù)時間歷程Fig.16 Time history of interpolated steady flow lift efficient

考慮到固體運動與流場之間的相互耦合作用,非定常流動模型更加符合實際情況。采用本文提出的非定常流固耦合分析方法,平板運動參數(shù)預(yù)先設(shè)定,平板轉(zhuǎn)動會導(dǎo)致周圍流場改變,氣動力也發(fā)生相應(yīng)的變化。基于非定常流動模型的升阻力系數(shù)計算結(jié)果如圖17和圖18所示,可以看到計算結(jié)果與試驗結(jié)果吻合較好。因此可以相信,在低雷諾數(shù)流動中,使用層流非定常模型能比較好地獲得空氣對平板的作用力。

圖17 考慮平板運動的非定常阻力系數(shù)時間歷程Fig.17 Time history drag coefficient of moving plate

3.3 流體-固體耦合算例分析

本算例旨在模擬貨臺在空氣中自由擺蕩,觀察氣動力對貨臺的影響,吊掛結(jié)構(gòu)和貨臺參數(shù)與3.1節(jié)算例相同,空氣來流速度為0.1 m/s沿著z方向。為了保證貨臺在有效流場網(wǎng)格范圍內(nèi)運動,構(gòu)建的流場尺寸為36 m×36 m×36 m,劃分的網(wǎng)格數(shù)目為499 381個,時間歷程為10 s.在本算例中,多體力學方程通過編寫基于平衡點法的多體力學程序進行求解。流體求解通過商業(yè)CFD求解器FLUENT進行。編寫的多體程序以用戶自定義函數(shù)(UDF)的形式經(jīng)FLUENT調(diào)用,使用聯(lián)想服務(wù)器(Think Station D30),通過30線程并行計算,計算時間約為48 h.

圖18 考慮平板運動的非定常升力系數(shù)時間歷程Fig.18 Time history of lift coefficient of moving plate

截取貨臺對稱平面 Oχz平面(如圖20所示),直觀展示貨臺的氣動力分布,貨臺在Oχz平面的壓力分布隨時間變化關(guān)系如圖19所示,在半個擺蕩周期內(nèi)貨臺上下表面的壓力分布不斷交替,貨臺χ方向氣動力合力隨時間變化曲線如圖21所示。由于定常流模型下氣動力僅是攻角的函數(shù),而考慮非定常流模型后氣動力除了與貨臺的攻角有關(guān),還與上一時刻流場參數(shù)和當前貨臺的運動狀態(tài)有關(guān),需要考慮固體和流體之間的耦合作用逐步計算獲得。因此,在圖21中,通過非定常流動計算的氣動力變化劇烈,而通過定常流動升阻力系數(shù)預(yù)測的氣動力卻只有較少的震蕩,整個氣動力曲線顯得光滑規(guī)整。對比顯示出,非定常流動模型下得到的氣動力小于定常流動模型的計算值,二者在多體系統(tǒng)中都起到粘滯阻力的作用。

貨臺在擺蕩過程中除了氣動力的影響以外,還受到自身重力以及吊掛結(jié)構(gòu)提供的牽引力的影響。由于貨臺重量及慣量相對較大,貨臺在自身重力、慣性力及牽引力作用下來回擺蕩,而氣動力則主要起到阻尼作用。為了考察耦合作用下空氣對貨臺的耗能效果,觀察采用不同的流動模型得到的貨臺速度幅值變化情況籍此了解能量的衰減,對比結(jié)果如圖22所示。在貨臺擺蕩過程中,定常流動模型與非定常流動模型都展示了良好的阻尼性質(zhì),兩種模型下系統(tǒng)的運動規(guī)律比較一致,但運動衰減幅度不同,這可以說明在空中擺蕩的貨臺有大部分能量通過加速流體的形式耗散在空氣中。其中,定常流動模型顯示出較大的阻尼效果,平板速度幅值衰減快于非定常流動模型。這在一定程度上可以解釋Peter在貨臺空投檢驗試驗中的發(fā)現(xiàn)[2]:貨臺仍然做較大范圍的擺蕩時,通過定常升阻力系數(shù)得到的仿真結(jié)果已經(jīng)趨于穩(wěn)定。

圖19 貨臺在Oχz平面的壓力分布Fig.19 Pressure distribution of platform on Oχz plane

圖20 貨臺的運動平面Fig.20 Oχz plane of moving platform

圖21 χ方向氣動力時間歷程Fig.21 Time history of aerodynamic force along χ direction

圖22 貨臺速度幅值時間變化曲線Fig.22 Platform velocity amplitude vs.time

4 結(jié)論

本文考慮貨臺大幅度擺蕩對空氣產(chǎn)生的非定常動效應(yīng),基于平衡點法和ALE格式有限體積法推導(dǎo)建立了非定常流動作用下的貨臺-空氣流體-固體耦合分析模型。通過數(shù)值算例研究得到了以下的結(jié)論:

1)基于平衡點法建立的繩索吊掛系統(tǒng)多體動力學模型具有迭代格式簡單、魯棒性較好的特點,能夠精確捕捉連接繩索的松弛-張緊關(guān)系。

2)在低雷諾數(shù)流動中,使用非定常層流模型能比較好地獲得流體對平板的作用力,非定常流動得到的氣動力更加符合實際,可作為空投貨臺擺蕩階段的流場分析方法。

3)通過流體-固體耦合分析得到非定常流動計算的氣動力變化劇烈,而通過定常流動升阻力系數(shù)計算得到的結(jié)果卻只有較少的震蕩,但其顯示出類似粘滯阻尼的特性耗能效果要比非定常模型強烈,這有可能是造成文獻中數(shù)值分析與試驗不相符的原因之一。

References)

[1]Peter C.A software simulation of cargo drop tests[C]//17th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar.Monterey,CA:AIAA,2003:2003-2132.

[2]Peter C,Kenneth D.Validation of a cargo airdrop software simulator [C]//17th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar.Monterey,CA:AIAA,2003:2003-2133.

[3]Kenneth D.The motion and aerodynamics of an airdrop platform [C]//22nd Applled Aerodynamics Conference and Exhibit.Provldence,RI:AIAA,2004:2004-4845.

[4]Mcqulling M,Potvin J,Riley J.Simulating the flows about cargo containers used during parachute airdrop operations[J].Journal of Aircraft,2011,48(4):1405-1411.

[5]Kenneth D.Aerodynamic forces on an airdrop platform[C]//18th Aerodynamic Decelerator Systems Technology Conference and Seminar.Munich,GE:AIAA,2005:2005-1634.

[6]Vaughan J,Kim D,Singhose W.Control of tower cranes with double-pendulum payload dynamics[J].IEEE Transactions on Control Systems Technology,2010,18(6):1345-1358.

[7]Masoud Z N,Nayfeh A H,Mook D T.Cargo pendulation reduction of ship-mounted cranes[J].Nonlinear Dynamics,2004,35(3):299-311.

[8]Track H H S T,Holloman N M.Designed experiments for nylon band characterization[C]//US Air Force T&E Days 2010.Nashville,TN:AIAA,2010:2010-1711.

[9]Thomas L,K R.Biaxially oriented nylon-6 as a long duration material[C]//International Technology Conference.Albuquerque,NM:AIAA,1991:1991-3659.

[10]Poole L R,Huckins E K.Evaluation of massless-spring modeling of suspension-line elasticity during the parachute unfurling process,NASA-TN-D-6671[R].WA,US:NASA Langley Research Center,1972.

[11]John F,Robert T,Behzad R.Multibody parachute flight simulations using singular perturbation theory[C]//20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar.Seattle,WA:AIAA,2009:2009-2920.

[12]Raiszadeh B.Mutibody parachute flight simulations for planetary entry trajectories using“equilibrium points"[J].Advances in the Astronautical Sciences,2003,114:913-923.

[13]程文科,秦子增,張曉今.具有倒“Y"型吊掛系統(tǒng)的物傘組合體動力分析[J].彈道學報,1998,10(2):10-14.CHENG Wen-ke,QIN Zi-zeng,ZHANG Xiao-jin.The dynamic analysis of a parachute and payload assembly with an inverted“Y"suspension system[J].Journal of Ballistics,1998,10(2): 10-14.(in Chinese)

Fluid-structure Interaction Modelling of Airdrop Cargo Platform Swinging

TANG Jian-hua,QIAN Lin-fang,XU Ya-dong
(School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing 210094,Jiangsu,China)

A platform-sling model is established to study the dynamic response of airdrop platform.An uniform equivalent point Newton-Rapshon iteration method is presented.In order to study the unsteady fluid behavior,the time domain is discretized,and an ALE base finite volume method is used to solve the NS equations.The aero-force is obtained and imposed on the platform.The equivalent point method is validated through commercial software RECURDYN.The results show that the equivalent point method can well represent the slack-taut cases of the sling system.It can also be found that the laminar flow model is used to describe the acting force of fluid on plate well at low Reynolds number.Finally,the comparisons are made between the steady and unsteady FSI models.The time history of unsteady fluid force displays severe oscillation while the steady fluid force varies evenly.

fluid mechanics;unsteady flow;equivalent point;airdrop cargo platform;fluid-structure interaction

TJ011

A

1000-1093(2016)01-0141-08

2015-03-22

國家自然科學基金項目(11472137)

湯健華(1984—),男,博士研究生。E-mail:jian_hua_tang@126.com;錢林方(1961—),男,教授,博士生導(dǎo)師。E-mail:lfqian@vip.163.com