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

湯健華,錢林方,徐亞棟

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

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

湯健華,錢林方,徐亞棟

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

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

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

DOI:10.3969/j.issn.1000-1093.2016.01.021

0 引言

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

圖1 空投貨臺及吊掛結構示意圖Fig.1 Schematic diagram of Airdrop sling system

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

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

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

1 貨臺多體動力學分析

1.1 貨臺動力學問題描述

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

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

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

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

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

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

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

1.2 平衡點法統一迭代格式

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

忽略連接節點自由度,多體動力學運動方程變為

1)固定連接節點的位置,求解其他部件的的響應。

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

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

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

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

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

式中:變化的位移為

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

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

由關系式

可以得到如下關系:

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

式中:I為單位矩陣。吊掛關于連接節點的雅克比矩陣為

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

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

圖5 可能導致奇異的情形Fig.5 Singularity phenomenon

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

2 非定常空投貨臺流體-固體耦合分析

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

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

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

3 算例

3.1 平衡點法分析

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

圖7 貨臺的初始狀態Fig.7 Initial state of platform

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

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

圖8 連接節點χ位置隨時間變化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]曾經通過試驗手段(見圖 11和圖12)獲得定常流動與非定常流動的升阻力系數,并指出定常流動下得到的升阻力系數與非定常流動下得到的測量結果相差較大。試驗的結構如圖11所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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 結論

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

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

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

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

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"型吊掛系統的物傘組合體動力分析[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—),男,教授,博士生導師。E-mail:lfqian@vip.163.com