臨近空間太陽能飛行器著陸狀態陣風響應

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

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

0 引言

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

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

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

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

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

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

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

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

圖2 CFD局部網格示意圖Fig.2 Schematic diagram of CFD local grid

2 計算方法

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

2.1 氣動模型

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

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

2.2 網格速度法

陣風模型使用網格速度法引入，流體區域內的流場速度可表示為：

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

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

2.3 陣風計算模型

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

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

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

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

以只有豎直方向速度的1-cos陣風模型為基礎，二維陣風模型考慮了陣風幅值沿展向的變化。因此，二維陣風速度型為[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 計算流程

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

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

3 計算結果及分析

3.1 陣風維數的影響

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

圖7 給出了在一維和二維陣風載荷作用下，飛行器達到靜平衡狀態后，受到陣風載荷時翼尖位移和扭轉角的對比。可看出：在一維和二維陣風載荷作用下，翼尖變形均出現較大的響應峰值，但翼尖z向位移響應的峰值方向相反，且二維陣風載荷作用下得到的峰值比一維陣風載荷作用下的小，出現時間相對滯后；此外，達到峰值后，一維陣風載荷作用下得到的翼尖位移和扭轉角恢復到平衡位置，而二維陣風載荷作用下得到的翼尖z向位移未恢復到靜平衡狀態，只是恢復到-0.09 m，翼尖扭轉角則減小到-1.99°。

圖7 翼尖位移和扭轉角陣風響應計算結果對比Fig.7 Gust response calculation results comparison for wingtip displacement and torsional angle

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

圖8 翼根所受結構彎矩陣風響應計算結果對比Fig.8 Gust response calculation results comparison for structural moment on the wing root

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

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

圖9 飛行器整體升力系數陣風響應計算結果對比Fig.9 Gust response calculation results comparison for lift coefficient of aircraft

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

3.2 展向分布因子的影響

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

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

4 結論

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

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

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

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