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

（1.南京航空航天大學航空學院，江蘇 南京 210016；2.中國商用飛機有限責任公司，上海 200126；3.南京航空航天大學通用航空與飛行學院，江蘇 溧陽 213300）

0 引言

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

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

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

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

1 臨近空間太陽能飛行器模型

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

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

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

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

2 計算方法

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

2.1 氣動模型

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

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

2.2 網(wǎng)格速度法

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

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

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

2.3 陣風計算模型

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

式（5）中：wg為陣風速度；Ude為陣風幅值；Lg為陣風長度；V為飛行器速度；t為計算時間（飛行時間）。

圖3 為飛行器所遭遇的1-cos 型陣風速度型[15]，設(shè)飛行器平均氣動弦長為c，無窮遠處來流速度為V∞，定義無量綱時間S=2V∞t c。

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

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

式（6）中：Lspan為飛行器翼展；a為陣風沿展向分布的因子；y為沿展向坐標。

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

圖4 二維陣風模型[17]Fig.4 Two-dimensional gust model

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

2.4 計算流程

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 展向分布因子的影響

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

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

4 結(jié)論

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

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

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

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