高興龍，陳 欽，*，張青斌，李志輝

(1. 中國空氣動力研究與發(fā)展中心，綿陽 621000；2. 國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073)

0 引 言

翼傘是一種雙層結(jié)構(gòu)的柔性矩形翼，上、下翼面用翼型的肋幅分隔成若干氣室，翼型前緣開口，在前進(jìn)飛行中形成“沖壓空氣”，維持若干個(gè)氣室的內(nèi)壓以保持翼型[1]。當(dāng)翼傘系統(tǒng)需要進(jìn)行機(jī)動轉(zhuǎn)彎和雀降等操縱動作時(shí)，會對翼型后緣進(jìn)行下拉偏轉(zhuǎn)操作來實(shí)現(xiàn)。

翼傘后緣偏轉(zhuǎn)的操縱過程會顯著改變翼面的整體氣動布局，同時(shí)需要多根操縱繩精確協(xié)同控制，是典型的氣動與結(jié)構(gòu)緊耦合問題，涉及到的動力學(xué)問題復(fù)雜多變。對于翼傘系統(tǒng)操縱過程的動力學(xué)機(jī)理問題研究一直是降落傘領(lǐng)域的關(guān)鍵技術(shù)和熱點(diǎn)問題。著名的美國X-38 太空成員返回飛行器利用翼傘后緣下偏進(jìn)行減速著陸的操縱，成功實(shí)現(xiàn)了太空人員返回安全著陸，是迄今為止世界上最大最先進(jìn)的翼傘研究項(xiàng)目[2]。除此之外，比較典型的應(yīng)用案例有美國研制的聯(lián)合精確空投系統(tǒng)（joint precision airdrop system,JPADS），該系統(tǒng)是一個(gè)高空、全天候、GPS 制導(dǎo)的精確空投系統(tǒng)，包括超輕型和輕型兩種布局，前后兩者的最大載重量分別約為998 kg 和4536 kg，空投精度優(yōu)于50～100 m。其他案例還有加拿大的“雪雁”（SnowGoose）偵察系統(tǒng)和“夏爾巴人”（Sherpa）精確空投系統(tǒng)、歐洲的SLG-Sys 自主滑翔傘降系統(tǒng)等[2]。

目前對于翼傘后緣偏轉(zhuǎn)的動力學(xué)行為研究主要依靠數(shù)值仿真和試驗(yàn)，其中數(shù)值仿真主要采用流固耦合仿真技術(shù)對沖壓翼傘動力學(xué)行為進(jìn)行預(yù)測。目前國內(nèi)外專門針對翼傘流固耦合過程的動力學(xué)機(jī)理研究已經(jīng)逐漸由傳統(tǒng)的完全試驗(yàn)法、半試驗(yàn)半理論法和完全理論法三種手段，轉(zhuǎn)向了基于先進(jìn)CFD 和CSD 方法的大規(guī)模并行數(shù)值仿真技術(shù)。Tang 等[3]基于LS-DYNA 軟件中的ICFD 流體求解器和結(jié)構(gòu)求解器仿真計(jì)算了翼傘后緣下偏過程的結(jié)構(gòu)變形行為。Tezduyar 所領(lǐng)導(dǎo)的團(tuán)隊(duì)基于DSD/SST 有限元方法開展了沖壓翼傘系統(tǒng)的流固耦合計(jì)算仿真[4]。Fogell等[5]基于松耦合、利用CFD 和CSD 方法對翼傘穩(wěn)態(tài)飛行的充氣外形進(jìn)行流固耦合仿真，研究了翼傘結(jié)構(gòu)變形和周圍的流場特性。Altmann[6]使用勢流理論和簡化有限元法分析了翼傘的流固耦合問題。近年來，國內(nèi)專門針對翼傘系統(tǒng)的流固耦合研究也逐漸增多，Nie 和Cao[7]基于Dirchlet/Neumann 迭代方法研究了翼傘充氣過程的流固耦合問題，騰海山團(tuán)隊(duì)[8,9]基于LS-DYNA 求解器對翼傘穩(wěn)態(tài)構(gòu)型進(jìn)行了流固耦合仿真，同時(shí)對大型翼傘操縱轉(zhuǎn)彎過程的氣動特性進(jìn)行仿真分析。汪龍芳等[10]基于索膜有限元模型和CFD 方法對翼傘穩(wěn)態(tài)飛行的氣動變形進(jìn)行了三維數(shù)值模擬，并分析了氣室的“鼓包”現(xiàn)象。孫青林等[11-12]研究了襟翼偏轉(zhuǎn)過程翼傘系統(tǒng)的氣動性能變化。郭一鳴等[13]基于松耦合方法和參數(shù)辨識開展了翼傘變形和下偏情況的動力學(xué)和氣動性能仿真。但是針對翼傘后緣下偏過程的氣動與結(jié)構(gòu)耦合動力學(xué)特性研究較少，本文分別考慮翼傘后緣偏轉(zhuǎn)過程的非線性氣動特性和結(jié)構(gòu)響應(yīng)行為，采用數(shù)值仿真與試驗(yàn)對比的方法深入分析了翼傘單側(cè)和雙側(cè)后緣偏轉(zhuǎn)過程的流固耦合動力學(xué)特性。

本文基于Structured ALE（S-ALE）流固耦合方法對翼傘后緣偏轉(zhuǎn)過程進(jìn)行動力學(xué)建模和仿真分析。研究翼傘三維模型后緣偏轉(zhuǎn)過程、傘衣結(jié)構(gòu)場和周圍流場的時(shí)變演化規(guī)律及分布特性，為進(jìn)一步指導(dǎo)大型翼傘精確空投系統(tǒng)的飛控系統(tǒng)設(shè)計(jì)和技術(shù)應(yīng)用提供參考。

1 流固耦合建模

1.1 問題描述及假設(shè)

本文所研究的翼傘后緣偏轉(zhuǎn)過程是針對充滿鼓包狀態(tài)的翼傘三維模型進(jìn)行的。翼傘系統(tǒng)包括傘衣、傘繩和掛重載荷，幾何模型如圖1 所示。實(shí)際流固耦合仿真過程只考慮傘衣結(jié)構(gòu)與流場的雙向耦合作用；傘繩在翼傘偏轉(zhuǎn)過程承受拉力，且通過傘繩施加后緣下拉過程的作用力載荷；忽略傘繩與周圍流體的耦合作用和繩索的阻尼效應(yīng)。

1.2 結(jié)構(gòu)動力學(xué)

翼傘傘衣為柔性織物結(jié)構(gòu)，一般由透氣量極低的抗撕錦絲綢材料制成，本文忽略傘衣織物結(jié)構(gòu)的透氣量，則結(jié)構(gòu)域的控制方程可以表示為[14]：

其中，ρs為材料密度，u為結(jié)構(gòu)質(zhì)點(diǎn)的速度矢量，σs為柯氏應(yīng)力張量，fs為作用在結(jié)構(gòu)的外部力， ?s代表結(jié)構(gòu)域的速度邊界。

1.3 流體動力學(xué)

翼傘系統(tǒng)周圍流體域的動力學(xué)方程為[15]：

其中，ρf、uf分別為密度和速度矢量，ff和 σf分別為外部力和應(yīng)力張量， ?f代表流體域的速度邊界。通過引入ALE 格式，網(wǎng)格可以自由移動。在參考構(gòu)型下不可壓縮流的ALE 格式控制方程為：

其中，流體速度uf=v-w，v和w分別為參考構(gòu)型下的流體質(zhì)點(diǎn)速度和材料網(wǎng)格點(diǎn)速度，g是單位密度的外力矢量。 當(dāng)v=w時(shí)，公式為拉格朗日形式；當(dāng)v≠w時(shí)，公式為歐拉格式。同時(shí)給出Dirichlet 和Neuman形式的邊界條件如下：

其中，n是法向矢量，h為流固耦合邊界接觸力矢量。

1.4 罰函數(shù)法

在翼傘系統(tǒng)的整個(gè)后緣偏轉(zhuǎn)過程中，傘衣結(jié)構(gòu)變形與周圍流場變化相互作用，通過流固耦合界面信息的傳遞保證相互作用過程的能量守恒。通常有限元模型很難實(shí)現(xiàn)流體和結(jié)構(gòu)網(wǎng)格的完全匹配，本文基于歐拉-拉格朗日描述下的罰函數(shù)法進(jìn)行流體與傘衣結(jié)構(gòu)間節(jié)點(diǎn)力信息的傳遞，該方法允許非匹配的流體和結(jié)構(gòu)網(wǎng)格[16]。

流固耦合界面為傘衣織物結(jié)構(gòu)表面，罰函數(shù)的穿透力fc計(jì)算公式為：

其中，ks代表罰函數(shù)的剛度系數(shù)；m為積分時(shí)間步；dm代表結(jié)構(gòu)節(jié)點(diǎn)在第m個(gè)時(shí)間步進(jìn)入流體域的穿透深度，與結(jié)構(gòu)和流場網(wǎng)格節(jié)點(diǎn)的相對運(yùn)動速度有關(guān)，可以在積分過程進(jìn)行迭代計(jì)算。

利用穿透深度和相對運(yùn)動速度可以在每個(gè)時(shí)間步內(nèi)判斷穿透是否發(fā)生，若穿透發(fā)生則在流固耦合界面上對流場和結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)施加大小相等、方向相反的作用力，以滿足流固耦合界面的受力平衡：

假設(shè)流體為滲透性介質(zhì)，采用顯式動力學(xué)積分方法對仿真模型進(jìn)行迭代計(jì)算，從而完成流固耦合信息的關(guān)聯(lián)和傳遞。

2 仿真方法驗(yàn)證

首先針對翼傘氣室單元后緣偏轉(zhuǎn)過程進(jìn)行仿真計(jì)算，初步驗(yàn)證本文所采用方法的計(jì)算收斂性和可行性。建立翼傘氣室結(jié)構(gòu)的三維有限元模型，并在后緣沿操縱繩方向施加下拉載荷。流場域完全包裹結(jié)構(gòu)模型，流場網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格。為避免因流體和結(jié)構(gòu)單元之間尺寸差異過大而導(dǎo)致顯式動力學(xué)積分過程可能出現(xiàn)的非物理特征“沙漏現(xiàn)象”，進(jìn)而引起計(jì)算發(fā)散，流場網(wǎng)格尺寸與結(jié)構(gòu)網(wǎng)格尺寸盡量接近1∶1，如圖2 所示。

圖2 翼傘氣室流固耦合仿真網(wǎng)格模型Fig. 2 Mesh model of the parafoil cell for the FSI simulation

本文采用S-ALE 求解方法對流固耦合模型進(jìn)行仿真計(jì)算，S-ALE 方法與傳統(tǒng)ALE 方法的基本理論相同，均包括了映射過程的對流輸運(yùn)、界面重構(gòu)和歐拉流場與拉格朗日結(jié)構(gòu)相互作用的流固耦合過程[17]。不同的是，在網(wǎng)格的處理方法上，S-ALE 方法采用自動生成網(wǎng)格技術(shù)，即流場網(wǎng)格根據(jù)控制點(diǎn)設(shè)定的方向、增長率、網(wǎng)格尺寸、網(wǎng)格密度等參數(shù)在仿真過程中隨著時(shí)間步的推進(jìn)逐漸產(chǎn)生，仿真前無需單獨(dú)建立流場網(wǎng)格。這可以極大減小網(wǎng)格處理時(shí)間并提高計(jì)算效率。經(jīng)過仿真測算，與傳統(tǒng)ALE 方法相比，S-ALE 方法的計(jì)算效率可以提高60%[18]。

計(jì)算的來流速度為15 m/s，來流攻角為0°，流場和結(jié)構(gòu)域的計(jì)算結(jié)果如圖3～圖5 所示。

圖3 翼傘氣室后緣偏轉(zhuǎn)過程流場速度變化云圖Fig. 3 Velocity contour variations of the parafoil cell under trailing-edge deflection

圖5 翼傘氣室后緣偏轉(zhuǎn)過程傘衣表面壓力分布云圖Fig. 5 Pressure contours of the parafoil cell under trailing-edge deflection

從圖3 流場演化結(jié)果中的速度云圖可以明顯看到渦流由上翼面向后緣逐漸過渡和破碎的過程。隨著后緣向下偏轉(zhuǎn)程度的增大，在后緣尾部的上端出現(xiàn)明顯的低速區(qū)域，形成高壓的逆壓梯度，從而引起流動分離現(xiàn)象。從圖4 和圖5 的傘衣結(jié)構(gòu)應(yīng)力和壓力分布云圖中可以看出，氣室中間鼓包區(qū)域?yàn)槭芰^為集中且強(qiáng)度高的區(qū)域。

將采用本文流固耦合方法仿真計(jì)算得到的下偏狀態(tài)翼傘氣室單元結(jié)構(gòu)有限元模型導(dǎo)出，作為幾何邊界模型，采用Fluent 對該導(dǎo)出的后緣偏轉(zhuǎn)氣室單元進(jìn)行氣動特性仿真計(jì)算驗(yàn)證，來流速度同樣取為15 m/s，湍流模型選擇標(biāo)準(zhǔn)k-ω模型，計(jì)算得到的翼傘氣室單元后緣下偏狀態(tài)翼型表面的壓力系數(shù)為0.65。采用本文流固耦合仿真方法計(jì)算得到的同樣狀態(tài)翼型表面壓力系數(shù)為0.62。兩種方法的結(jié)果較為接近，初步驗(yàn)證了本文方法的有效性和可靠性。

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

翼傘氣室單元仿真計(jì)算初步驗(yàn)證了本文方法的可行性和收斂性，將本文方法進(jìn)一步應(yīng)用于整個(gè)翼傘系統(tǒng)多氣室后緣偏轉(zhuǎn)過程的流固耦合仿真。為了模擬翼傘機(jī)動轉(zhuǎn)彎和雀降著陸兩個(gè)經(jīng)典操縱動作，分別對整個(gè)翼傘系統(tǒng)模型的后緣單側(cè)下偏和雙側(cè)下偏過程進(jìn)行仿真計(jì)算。

3.1 流固耦合結(jié)果

圖6 和圖7 為翼傘后緣單側(cè)和雙側(cè)下偏操縱過程的傘衣表面結(jié)構(gòu)Von Mises 應(yīng)力分布云圖。從圖中可以明顯看出傘衣的上翼面整體應(yīng)力分布較為均勻，但當(dāng)后緣下偏時(shí)，在偏轉(zhuǎn)與未偏轉(zhuǎn)的轉(zhuǎn)折交接區(qū)域會出面明顯的應(yīng)力集中且受力增大；同時(shí)翼型前緣也相應(yīng)出現(xiàn)類似現(xiàn)象，且與后緣應(yīng)力集中區(qū)域的受力水平接近。說明翼傘后緣偏轉(zhuǎn)發(fā)生的轉(zhuǎn)折交接區(qū)域和偏轉(zhuǎn)側(cè)前緣切口區(qū)域是在翼傘操縱過程比較容易發(fā)生撕裂現(xiàn)象的區(qū)域，應(yīng)考慮在這兩處增加結(jié)構(gòu)強(qiáng)度。從傘衣結(jié)構(gòu)受力水平來講，在下偏量相同的情況下，雙側(cè)下偏的傘衣結(jié)構(gòu)受力明顯高于單側(cè)下偏情況。

圖6 不同時(shí)刻單側(cè)下偏傘衣表面結(jié)構(gòu)Von Mises 應(yīng)力云圖Fig. 6 Von Mises stress contours of the canopy under single deflection at different time instances

圖7 不同時(shí)刻雙側(cè)下偏傘衣表面結(jié)構(gòu)Von Mises 應(yīng)力云圖Fig. 7 Von Mises stress contours of the canopy under double deflection at different time instances

圖8 和圖9 為翼傘后緣雙側(cè)下偏過程不同時(shí)刻的傘衣周圍流場速度變化云圖。從單側(cè)下偏過程的速度變化分布云圖可以看出，在翼傘的尾緣頂部緊貼上翼面附近出現(xiàn)近似圓形的高速區(qū)域，當(dāng)下偏時(shí)該區(qū)域被拉長，下偏完成后此區(qū)域減小。該高速區(qū)域會使上翼面頂部的氣流自后緣向前緣倒流，產(chǎn)生流動分離，這也是翼傘轉(zhuǎn)彎操縱出現(xiàn)失穩(wěn)現(xiàn)象的主要原因。從雙側(cè)下偏過程的橫向剖面流場速度分布云圖可以看出，該高速區(qū)域主要存在于翼傘下偏發(fā)生的轉(zhuǎn)折交界區(qū)域附近。

圖8 不同時(shí)刻單側(cè)下偏傘衣周圍流場速度變化云圖Fig. 8 Velocity contours around the canopy under single deflection at different time instances

圖9 不同時(shí)刻雙側(cè)下偏過程周圍流場速度變化云圖Fig. 9 Velocity contours around the canopy under double deflection at different time instances

3.2 性能曲線結(jié)果

圖10 為翼傘后緣單側(cè)和雙側(cè)下偏過程所施加的下偏量隨時(shí)間的變化曲線。兩種情況的下偏量相同。圖11～圖13 分別為兩種情況的翼傘氣動阻力面積、升力系數(shù)、阻力系數(shù)三者隨時(shí)間變化曲線。后緣單側(cè)下偏對應(yīng)翼傘的操縱轉(zhuǎn)彎動作，下拉量開始時(shí)間為0.2 s，此時(shí)升力系數(shù)出現(xiàn)一個(gè)峰值，下偏量達(dá)到最大的時(shí)刻為0.3 s，此時(shí)阻力系數(shù)最大，阻力面積的變化與下偏量變化結(jié)果趨勢較為一致。

圖10 翼傘后緣下偏過程下偏量隨時(shí)間變化Fig. 10 Time variation of the deflection for the parafoil under trailing-edge deflection

圖11 翼傘后緣下偏過程氣動阻力面積隨時(shí)間變化Fig. 11 Time variation of the drag area for the parafoil under trailing-edge deflection

圖12 翼傘后緣下偏過程氣動升力系數(shù)隨時(shí)間變化Fig. 12 Time variation of the lift coefficient of the parafoil under trailing-edge deflection

圖13 翼傘后緣下偏過程氣動阻力系數(shù)隨時(shí)間變化Fig. 13 Time variation of the drag coefficient of the parafoil under trailing-edge deflection

雙側(cè)下偏過程對應(yīng)翼傘的雀降著陸操縱動作，翼傘下偏量達(dá)到最大的時(shí)刻為0.2 s，但翼傘最大阻力系數(shù)出現(xiàn)的時(shí)間約為0.44 s，最大升力系數(shù)出現(xiàn)的時(shí)間約為0.38 s，均晚于下拉操縱量達(dá)到最大的時(shí)刻，表明沖壓翼傘后緣雙側(cè)下偏操縱過程具有延遲響應(yīng)的特性。

3.3 試驗(yàn)驗(yàn)證

為進(jìn)一步驗(yàn)證本文仿真結(jié)果，將流固耦合仿真計(jì)算獲得的后緣雙側(cè)下偏過程達(dá)到不同下偏量的翼傘系統(tǒng)結(jié)構(gòu)有限元模型輸出保存，利用ANSYS 軟件將有限元模型轉(zhuǎn)換為翼傘結(jié)構(gòu)外輪廓的實(shí)體幾何模型，并采用3D 打印技術(shù)打印出實(shí)體模型。

在0.55 m×0.4 m 低噪聲航空聲學(xué)風(fēng)洞對翼傘模型進(jìn)行靜態(tài)測壓吹風(fēng)試驗(yàn)，風(fēng)速同樣為15 m/s。在翼傘每個(gè)氣室的翼型模型中剖面布置一定數(shù)量的測壓孔，對測得的壓力分布在升力方向積分，得到該狀態(tài)翼傘系統(tǒng)的升力系數(shù)。

表1 為翼傘后緣偏轉(zhuǎn)模型升力系數(shù)仿真與試驗(yàn)的對比結(jié)果，其中仿真結(jié)果為利用本文方法對不同下偏量的翼傘結(jié)構(gòu)外輪廓模型進(jìn)行的流固耦合仿真。從對比可以看出仿真結(jié)果與試驗(yàn)結(jié)果較為一致，且不同下偏量對應(yīng)的升力系數(shù)變化趨勢相同。

表1 翼傘后緣偏轉(zhuǎn)模型升力系數(shù)對比結(jié)果Table 1 Lift coefficient comparison of the parafoil under trailing-edge deflection

4 結(jié) 論

本文研究了翼傘后緣下偏過程的流固耦合動力學(xué)問題。在單個(gè)氣室翼型的后緣偏轉(zhuǎn)過程仿真研究中發(fā)現(xiàn)偏轉(zhuǎn)過程氣室周圍流場出現(xiàn)分離流動現(xiàn)象。之后對整個(gè)翼傘的后緣下偏過程進(jìn)行流固耦合仿真，并利用風(fēng)洞試驗(yàn)對仿真結(jié)果進(jìn)行驗(yàn)證，研究結(jié)果表明基于S-ALE 和罰函數(shù)法的流固耦合技術(shù)可以進(jìn)行翼傘后緣偏轉(zhuǎn)過程的流固耦合動力學(xué)行為研究：1）傘衣結(jié)構(gòu)計(jì)算結(jié)果表明后緣下偏過程應(yīng)力集中區(qū)域位于翼型前緣切口和后緣下偏轉(zhuǎn)折發(fā)生區(qū)域；2）通過分析氣動參數(shù)的時(shí)間歷程演化趨勢，發(fā)現(xiàn)沖壓翼傘的后緣下偏過程會出現(xiàn)明顯的延遲響應(yīng)現(xiàn)象?？傊帽疚乃捎玫牧鞴恬詈戏椒梢詫σ韨憧胀断到y(tǒng)的氣動結(jié)構(gòu)耦合動力學(xué)行為進(jìn)行較好的預(yù)測，有助于翼傘精確空投系統(tǒng)的技術(shù)驗(yàn)證和設(shè)計(jì)分析。