臨近空間太陽(yáng)能飛行器著陸狀態(tài)陣風(fēng)響應(yīng)

吳金洋，吳光輝，魏小輝，聶 宏，房興波，陳 虎

（1.南京航空航天大學(xué)航空學(xué)院，江蘇 南京 210016；2.中國(guó)商用飛機(jī)有限責(zé)任公司，上海 200126；3.南京航空航天大學(xué)通用航空與飛行學(xué)院，江蘇 溧陽(yáng) 213300）

0 引言

臨近空間太陽(yáng)能飛行器通常巡航飛行在20～100 km 的臨近空間區(qū)域。與常規(guī)布局飛行器相比，該類飛行器往往具有尺寸大、質(zhì)量輕、翼展長(zhǎng)、機(jī)體柔性大等幾何非線性結(jié)構(gòu)特點(diǎn)，以及高升力、低阻力等氣動(dòng)性能，且著陸過(guò)程中飛行速度極低。這些因素極易在風(fēng)切變、湍流等突發(fā)氣流載荷條件下使機(jī)翼產(chǎn)生大的彎曲和扭轉(zhuǎn)形變，引起機(jī)翼氣動(dòng)載荷分布的改變，進(jìn)而對(duì)整個(gè)飛行器的氣動(dòng)特性產(chǎn)生較大影響，甚至導(dǎo)致飛行器解體墜毀事故。2016 年，美國(guó)Facebook 公司的Aquila 太陽(yáng)能無(wú)人機(jī)在降落初期由于遭遇陣風(fēng)載荷，最后導(dǎo)致墜毀解體。基于此，研究臨近空間太陽(yáng)能飛行器起降過(guò)程中的氣動(dòng)彈性陣風(fēng)響應(yīng)問(wèn)題對(duì)其安全性具有重要意義。

針對(duì)臨近空間太陽(yáng)能飛行器的陣風(fēng)響應(yīng)問(wèn)題，國(guó)內(nèi)外研究人員開(kāi)展了大量的研究工作，如使用非線性梁結(jié)構(gòu)和大迎角氣動(dòng)力模型[1-2]，準(zhǔn)模態(tài)法和非定常渦格法[3-4]，準(zhǔn)模態(tài)法和修正Theodorsen 方法結(jié)合片條理論[5]等進(jìn)行結(jié)構(gòu)與氣動(dòng)建模，探討了大柔性飛行器在不同陣風(fēng)載荷下結(jié)構(gòu)和氣動(dòng)特性的變化。近年來(lái)，基于歐拉方程或N-S 方程的CFD 與CSD 相互耦合解決非線性氣動(dòng)彈性問(wèn)題的精確計(jì)算發(fā)展迅速[6]。在研究幾何非線性大柔性飛行器陣風(fēng)響應(yīng)時(shí)，CFD/CSD耦合展現(xiàn)出良好的性能，能夠準(zhǔn)確反映飛行器在大位移及結(jié)構(gòu)大變形條件下的動(dòng)態(tài)特性[7]。

為了對(duì)飛行器陣風(fēng)響應(yīng)進(jìn)行理論分析并設(shè)計(jì)相應(yīng)的減緩控制系統(tǒng)，陣風(fēng)模型也是1 個(gè)關(guān)鍵因素。早期研究者主要采用離散陣風(fēng)模型[8]，如“銳邊”陣風(fēng)和“斜坡”陣風(fēng)，它們可描述大氣紊流中的峰值、飛行器尾流區(qū)的流動(dòng)以及地形誘導(dǎo)的氣流等[9]。隨著陣風(fēng)研究的深入，“銳邊”和“斜坡”陣風(fēng)逐漸被1-cos 型離散陣風(fēng)取代，該類型陣風(fēng)至今仍在使用，并被FAR-25、CCAR-25 等適航法規(guī)作為分析飛行器陣風(fēng)載荷的理論模型。由于空間尺度上非均勻陣風(fēng)的復(fù)雜性，上述陣風(fēng)模型都基于一維均勻陣風(fēng)場(chǎng)假設(shè)進(jìn)行了簡(jiǎn)化，即認(rèn)為飛行器陣風(fēng)速度僅在飛行方向上變化，沿展向方向保持不變。然而，對(duì)于擁有較大展向尺寸的大展弦比柔性飛行器，采用一維均勻陣風(fēng)場(chǎng)模型進(jìn)行分析存在缺陷，需要進(jìn)一步考慮展向陣風(fēng)變化[10-11]。

為此，本文基于CFD/CSD 松耦合分析法，利用網(wǎng)格速度法引入陣風(fēng)載荷，研究了二維陣風(fēng)對(duì)臨近空間太陽(yáng)能飛行器著陸狀態(tài)時(shí)陣風(fēng)響應(yīng)的影響，并將二維陣風(fēng)響應(yīng)與一維陣風(fēng)響應(yīng)進(jìn)行對(duì)比。研究獲得了臨近空間太陽(yáng)能飛行器著陸時(shí)，二維陣風(fēng)沿飛行器翼展方向的變化對(duì)其結(jié)構(gòu)和氣動(dòng)性能的影響規(guī)律。

1 臨近空間太陽(yáng)能飛行器模型

參考Facebook Aquila 無(wú)人機(jī)的設(shè)計(jì)，本文所計(jì)算的臨近空間太陽(yáng)能飛行器模型如圖1所示。

圖1 臨近空間太陽(yáng)能飛行器模型Fig.1 Model of near-space solar-powered aircraft

該太陽(yáng)能飛行器采用內(nèi)、外翼飛翼式布局；飛行器的展長(zhǎng)為48.00 m，平均氣動(dòng)弦長(zhǎng)為1.47 m，展弦比為32.65，總質(zhì)量為280 kg。飛行器所用材料及參數(shù)見(jiàn)文獻(xiàn)[12]。利用Ansys 軟件進(jìn)行模型構(gòu)建。建模時(shí)，針對(duì)著陸狀態(tài)，飛行器模型有一定迎角，且來(lái)流速度很低，空氣密度接近海平面。鑒于計(jì)算對(duì)象具有對(duì)稱特性，為提升計(jì)算效率，僅對(duì)一側(cè)流場(chǎng)進(jìn)行分析，為此選擇半實(shí)體模型。為降低流體網(wǎng)格數(shù)量，構(gòu)建流體域時(shí)使前后平面與飛行器前緣平行。借助Ansys Mechanical 進(jìn)行流體和固體網(wǎng)格劃分，對(duì)流固耦合面網(wǎng)格進(jìn)行細(xì)化處理。圖2展示了模型部分網(wǎng)格示意圖。

圖2 CFD局部網(wǎng)格示意圖Fig.2 Schematic diagram of CFD local grid

2 計(jì)算方法

盡管所研究的臨近空間太陽(yáng)能飛行器巡航階段低雷諾數(shù)效應(yīng)顯著，但其在起降過(guò)程中雷諾數(shù)超過(guò)百萬(wàn)，屬于高雷諾數(shù)范圍。據(jù)此，采用CFD/CSD 松耦合方法進(jìn)行求解。在CFD 部分，利用基于SST 湍流模型的有限體積法對(duì)可壓縮積分形式的雷諾平均Navier-Stokes（RANS）方程組進(jìn)行計(jì)算。在空間離散過(guò)程中：對(duì)流通量項(xiàng)采用Roe 格式的二階精度迎風(fēng)格式；黏性通量項(xiàng)選用二階精度的中心離散方式；時(shí)間推進(jìn)則使用LU-SGS隱式時(shí)間推進(jìn)技術(shù)。

2.1 氣動(dòng)模型

SST 湍流模型是將標(biāo)準(zhǔn)k-ω模型與標(biāo)準(zhǔn)k-ε模型通過(guò)混合函數(shù)相結(jié)合而形成[13]。其在繼承k-ε模型的穩(wěn)定性、計(jì)算時(shí)間少以及k-ω模型無(wú)需壁面函數(shù)即可進(jìn)行壁面積分等優(yōu)勢(shì)的基礎(chǔ)上，剔除了k-ε模型對(duì)逆壓力梯度和邊界層分離不敏感以及k-ω模型收斂困難等不足。因此，該模型在近壁區(qū)域和自由剪切層中展現(xiàn)出優(yōu)良的數(shù)值特性，并且能夠較好地模擬具有較大逆壓梯度特征的流動(dòng)。SSTk-ω湍流模型的通用控制方程為：

式（1）（2）中：t為計(jì)算時(shí)間（飛行時(shí)間）；ρ為密度；k為湍流動(dòng)能；ω為湍流特殊耗散；uˉi、uˉj為湍流速度平均值；xi、xj為坐標(biāo)分量；Γk、Γω為有效擴(kuò)散系數(shù)；Fk、Fω為湍流生成項(xiàng)；Yk、Yω為k、ω的耗散項(xiàng)；Dω為擴(kuò)散項(xiàng)。

2.2 網(wǎng)格速度法

陣風(fēng)模型使用網(wǎng)格速度法引入，流體區(qū)域內(nèi)的流場(chǎng)速度可表示為：

式（3）中：u、v和w為直角坐標(biāo)系中的速度分量；xt、yt和zt分別為3個(gè)方向上的網(wǎng)格速度分量；i、j和k代表直角坐標(biāo)系在各個(gè)方向上的單位向量。

當(dāng)流體區(qū)域受到1 個(gè)向上、速度為wi的陣風(fēng)影響時(shí)，z方向上的流體區(qū)域網(wǎng)格會(huì)出現(xiàn)1個(gè)速度為wi的突增，而實(shí)際上流體域網(wǎng)格未發(fā)生運(yùn)動(dòng)。此時(shí)，流體域相當(dāng)于在網(wǎng)格不動(dòng)的情況下整體都受到-wi的來(lái)流作用。據(jù)此，在陣風(fēng)作用下，流體域中流場(chǎng)的速度可以表示為：

2.3 陣風(fēng)計(jì)算模型

一維陣風(fēng)模型采用1-cos工程陣風(fēng)模型[14]，整個(gè)飛行器在豎直方向上受到如式（5）所示的陣風(fēng)影響。

式（5）中：wg為陣風(fēng)速度；Ude為陣風(fēng)幅值；Lg為陣風(fēng)長(zhǎng)度；V為飛行器速度；t為計(jì)算時(shí)間（飛行時(shí)間）。

圖3 為飛行器所遭遇的1-cos 型陣風(fēng)速度型[15]，設(shè)飛行器平均氣動(dòng)弦長(zhǎng)為c，無(wú)窮遠(yuǎn)處來(lái)流速度為V∞，定義無(wú)量綱時(shí)間S=2V∞t c。

圖3 一維1-cos型陣風(fēng)速度型[15]Fig.3 One-dimension 1-cos gust speed pattern

以只有豎直方向速度的1-cos陣風(fēng)模型為基礎(chǔ)，二維陣風(fēng)模型考慮了陣風(fēng)幅值沿展向的變化。因此，二維陣風(fēng)速度型為[16-17]：

式（6）中：Lspan為飛行器翼展；a為陣風(fēng)沿展向分布的因子；y為沿展向坐標(biāo)。

圖4 為a=1.0 時(shí)的陣風(fēng)模型；圖5 為不同a值下陣風(fēng)沿展向的分布。由圖可看出，隨著a的增大，二維陣風(fēng)向一維陣風(fēng)逼近，當(dāng)a接近無(wú)窮大時(shí)，二維陣風(fēng)退化為一維陣風(fēng)。

圖4 二維陣風(fēng)模型[17]Fig.4 Two-dimensional gust model

圖5 二維陣風(fēng)的展向分布[17]Fig.5 Two-dimensional gust spanwise distribution

2.4 計(jì)算流程

本研究中，即便降落時(shí)飛行器未遭遇陣風(fēng)，機(jī)翼本身也處在一定變形下的平衡狀態(tài)。為此，可建立如圖6 所示的CFD/CSD 松耦合方式的計(jì)算流程。首先，通過(guò)CFD/CSD 松耦合迭代收斂求出飛行器降落時(shí)的靜平衡狀態(tài)；然后，在平衡后某一時(shí)刻添加陣風(fēng)載荷；最后，通過(guò)CFD/CSD松耦合迭代，直到收斂為止。

圖6 松耦合計(jì)算流程圖Fig.6 Flowchart of loose coupling calculation

3 計(jì)算結(jié)果及分析

3.1 陣風(fēng)維數(shù)的影響

取飛行器飛行高度為500 m，飛行速度為40 km/h（11.1 m/s），大氣密度為1.167 kg/m3，降落時(shí)初始迎角為4°，先算出飛行器的靜平衡狀態(tài)，然后以此為基礎(chǔ)，計(jì)算到10 s 時(shí)分別施加一維和二維陣風(fēng)載荷，研究不同陣風(fēng)維數(shù)對(duì)飛行器著陸時(shí)陣風(fēng)響應(yīng)的影響。其中，陣風(fēng)幅值Ude1為6 m/s，陣風(fēng)長(zhǎng)度Lg為25 倍飛行器平均氣動(dòng)弦長(zhǎng)[18]，二維陣風(fēng)沿展向分布因子設(shè)為1.0。參照?qǐng)D4，在-0.25～0.25 飛行器翼展的跨度范圍內(nèi)，陣風(fēng)速度方向豎直向上，而其余跨度范圍內(nèi)的陣風(fēng)速度方向豎直向下。

圖7 給出了在一維和二維陣風(fēng)載荷作用下，飛行器達(dá)到靜平衡狀態(tài)后，受到陣風(fēng)載荷時(shí)翼尖位移和扭轉(zhuǎn)角的對(duì)比?？煽闯觯涸谝痪S和二維陣風(fēng)載荷作用下，翼尖變形均出現(xiàn)較大的響應(yīng)峰值，但翼尖z向位移響應(yīng)的峰值方向相反，且二維陣風(fēng)載荷作用下得到的峰值比一維陣風(fēng)載荷作用下的小，出現(xiàn)時(shí)間相對(duì)滯后；此外，達(dá)到峰值后，一維陣風(fēng)載荷作用下得到的翼尖位移和扭轉(zhuǎn)角恢復(fù)到平衡位置，而二維陣風(fēng)載荷作用下得到的翼尖z向位移未恢復(fù)到靜平衡狀態(tài)，只是恢復(fù)到-0.09 m，翼尖扭轉(zhuǎn)角則減小到-1.99°。

圖7 翼尖位移和扭轉(zhuǎn)角陣風(fēng)響應(yīng)計(jì)算結(jié)果對(duì)比Fig.7 Gust response calculation results comparison for wingtip displacement and torsional angle

圖8 給出了在一維和二維陣風(fēng)載荷作用下，飛行器達(dá)到靜平衡狀態(tài)后，受到陣風(fēng)載荷時(shí)翼根所受的結(jié)構(gòu)整體彎矩的對(duì)比。

圖8 翼根所受結(jié)構(gòu)彎矩陣風(fēng)響應(yīng)計(jì)算結(jié)果對(duì)比Fig.8 Gust response calculation results comparison for structural moment on the wing root

可看出：二維陣風(fēng)載荷作用下得到的結(jié)構(gòu)彎矩響應(yīng)結(jié)果峰值比一維陣風(fēng)載荷作用下的小，且方向相反，出現(xiàn)時(shí)間相對(duì)滯后；此外，二維陣風(fēng)載荷作用下得到的彎矩在達(dá)到峰值后未恢復(fù)到靜平衡狀態(tài)，而是振蕩恢復(fù)到-900 N·m左右，導(dǎo)致二維陣風(fēng)載荷作用下得到的翼尖變形減小到負(fù)值。這是由于展向分布因子為1.0 時(shí)，翼尖部分所受的二維陣風(fēng)速度方向豎直向下，與一維陣風(fēng)速度方向相反所致。

圖9 給出了在一維和二維陣風(fēng)載荷作用下，飛行器達(dá)到靜平衡狀態(tài)后，受到陣風(fēng)載荷時(shí)整體升力系數(shù)的對(duì)比。

圖9 飛行器整體升力系數(shù)陣風(fēng)響應(yīng)計(jì)算結(jié)果對(duì)比Fig.9 Gust response calculation results comparison for lift coefficient of aircraft

可看出：二維陣風(fēng)載荷作用下得到的結(jié)構(gòu)整體升力系數(shù)響應(yīng)峰值比一維陣風(fēng)載荷作用下的小，且出現(xiàn)時(shí)間相對(duì)滯后；此外，二維陣風(fēng)載荷作用下得到的整體升力系數(shù)在達(dá)到峰值后未恢復(fù)到靜平衡狀態(tài)，而是振蕩恢復(fù)到0.43 左右。這是由于飛行器受到展向分布因子為1.0 的二維陣風(fēng)載荷后，整體未恢復(fù)到靜平衡狀態(tài)，而是出現(xiàn)反方向扭轉(zhuǎn)，使得飛行器整體迎角減小所致。

3.2 展向分布因子的影響

圖10給出了飛行器達(dá)到靜平衡狀態(tài)后，展向分布因子分別取0.5、1.0、2.0 和5.0 時(shí)，二維陣風(fēng)載荷作用下翼尖z向位移、扭轉(zhuǎn)角、翼根所受結(jié)構(gòu)彎矩和整體升力系數(shù)的對(duì)比。結(jié)合圖7～9 可看出：隨著展向分布因子的增加，飛行器在二維陣風(fēng)載荷作用下的陣風(fēng)響應(yīng)逐漸接近于一維陣風(fēng)作用下的陣風(fēng)響應(yīng)。當(dāng)a≥2.0時(shí)，響應(yīng)峰值隨a的增加逐漸增大，接近受到一維陣風(fēng)作用下的響應(yīng)峰值；尤其當(dāng)a=5.0 時(shí)，二維陣風(fēng)的響應(yīng)峰值和響應(yīng)趨勢(shì)與一維陣風(fēng)的具有較好的一致性；而當(dāng)a ＜2.0時(shí)，由于陣風(fēng)存在速度方向豎直向下的部分，飛行器的陣風(fēng)響應(yīng)過(guò)程與受一維陣風(fēng)載荷作用下的響應(yīng)過(guò)程差異較大；當(dāng)a=0.5時(shí)，所有的響應(yīng)結(jié)果峰值較小，除扭轉(zhuǎn)角響應(yīng)結(jié)果出現(xiàn)2 個(gè)峰值外，翼尖位移、翼根結(jié)構(gòu)彎矩以及升力系數(shù)響應(yīng)結(jié)果峰值方向均與受一維陣風(fēng)載荷作用下的方向相反；而當(dāng)a=1.0時(shí)，所有的響應(yīng)結(jié)果達(dá)到峰值后未恢復(fù)到平衡狀態(tài)，且翼尖z向位移和翼根所受結(jié)構(gòu)彎矩的響應(yīng)結(jié)果的峰值方向均與受一維陣風(fēng)載荷作用下的響應(yīng)方向相反，并達(dá)到最大。

圖10 不同展向分布因子下飛行器陣風(fēng)響應(yīng)結(jié)果Fig.10 Aircraft gust response results for different spanwise distribution factor

4 結(jié)論

針對(duì)臨近空間太陽(yáng)能飛行器著陸時(shí)可能遭遇的二維陣風(fēng)載荷，基于CFD/CSD 松耦合分析法，利用網(wǎng)格速度法引入陣風(fēng)載荷，以1-cos 陣風(fēng)模型為基礎(chǔ)，探討了二維陣風(fēng)載荷對(duì)臨近空間太陽(yáng)能飛行器著陸狀態(tài)的陣風(fēng)響應(yīng)的影響，并將二維響應(yīng)與一維響應(yīng)進(jìn)行了對(duì)比，得出結(jié)論如下。

1） 相較于傳統(tǒng)的一維陣風(fēng)，飛行器受到二維陣風(fēng)作用下的陣風(fēng)響應(yīng)更為復(fù)雜，且受展向分布因子的影響較大。當(dāng)a≥2.0 時(shí)，二維陣風(fēng)響應(yīng)隨a的增加逐漸趨近于一維陣風(fēng)響應(yīng)；當(dāng)a ＜2.0時(shí)，由于二維陣風(fēng)速度存在方向豎直向下的部分，飛行器的響應(yīng)與受到一維陣風(fēng)作用下的響應(yīng)相差較大。

2） 當(dāng)二維陣風(fēng)的展向分布因子為1.0時(shí)，翼尖z向位移和翼根所受結(jié)構(gòu)彎矩的響應(yīng)結(jié)果峰值方向均與受一維陣風(fēng)載荷作用下的響應(yīng)方向相反，并達(dá)到最大，且所有響應(yīng)結(jié)果達(dá)到峰值后未恢復(fù)到平衡狀態(tài)。為此，研究臨近空間太陽(yáng)能飛行器著陸狀態(tài)陣風(fēng)響應(yīng)時(shí)，需重點(diǎn)關(guān)注a=1.0的二維陣風(fēng)模型。

3） 臨近空間太陽(yáng)能飛行器結(jié)構(gòu)質(zhì)量輕、機(jī)體柔性大，降落速度低，在著陸過(guò)程中易受到陣風(fēng)載荷的影響，對(duì)整個(gè)飛行器的結(jié)構(gòu)和飛行軌跡具有重大影響。建議應(yīng)對(duì)該類飛行器著陸過(guò)程中的陣風(fēng)響應(yīng)和減緩技術(shù)進(jìn)行廣泛深入的研究。