趙旋，彭緒浩，鄧子辰，張偉偉

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

0 引 言

現(xiàn)有風(fēng)洞試驗中的載荷測量主要包括直接力測量和表面壓力測量，飛行器風(fēng)洞縮比模型也因此分為測力模型和測壓模型?？傮w而言，風(fēng)洞試驗的可信度較高，獲取的氣動力/載荷結(jié)果往往作為考核數(shù)值仿真方法精度的標(biāo)模。然而，由于風(fēng)洞試驗的周期較長，對試驗人員的經(jīng)驗性依賴較高，試驗方案的合理性直接影響氣動載荷的獲取效率和效果。現(xiàn)有的工程實踐普遍認(rèn)為翼型表面至少需要布置50～100 個測壓孔，經(jīng)過壓力分布積分得到的升力和俯仰力矩系數(shù)精度才較為可信。為獲得翼型表面完整的流場信息，傳統(tǒng)方法通常在翼型表面布置足夠多的測壓孔進行風(fēng)洞試驗，通過簡單的插值重構(gòu)獲得翼型全表面的壓力分布，這往往需要較多的測壓孔。對于復(fù)雜飛行器，受限于空間位置和試驗成本，測壓數(shù)據(jù)獲得不充分，使得這種傳統(tǒng)的方法精度不夠；且復(fù)雜飛行器跨聲速風(fēng)洞試驗的氣動載荷受參數(shù)影響敏感，精細(xì)化的氣動力/載荷的測量難度更大、周期更長，因此亟需發(fā)展一種利用有限測量點流場數(shù)據(jù)重構(gòu)復(fù)雜飛行器表面流場信息的方法。

數(shù)據(jù)融合技術(shù)針對多種來源信息（同質(zhì)或異質(zhì)），根據(jù)某種標(biāo)準(zhǔn)在空間或時間上進行組合，獲得被測對象的一致性描述，并使得該融合系統(tǒng)具有更好的性能。目前，在空氣動力學(xué)領(lǐng)域中，已經(jīng)開展了多種多源信息融合方面的研究。Belyaev 等利用CFD 計算數(shù)據(jù)與風(fēng)洞試驗數(shù)據(jù)融合構(gòu)建了非均勻、非完整試驗狀態(tài)空間中的代理模型，并利用風(fēng)洞試驗數(shù)據(jù)驗證了模型的精度。Ghoreyshi 等利用試驗數(shù)據(jù)與 CFD 計算結(jié)果構(gòu)建了飛行器的高維變精度氣動力模型。王文正等提出了基于數(shù)學(xué)模型的氣動力數(shù)據(jù)融合準(zhǔn)則和方法，以氣動力數(shù)據(jù)滿足的準(zhǔn)則為依據(jù)，將不同類型的數(shù)據(jù)進行融合。Wang 等利用本征正交分解（POD）技術(shù)分別得到低、高精度流場數(shù)據(jù)的低維表示，通過 Kriging 模型建立兩者的映射關(guān)系，實現(xiàn)了不同狀態(tài)的定常流場融合預(yù)測。Kou 等采用多核神經(jīng)網(wǎng)絡(luò)模型將變精度模型推廣到非定常氣動力模型中，實現(xiàn)了利用低精度歐拉結(jié)果對 N–S 數(shù)值結(jié)果的逼近。He 等采用Kriging 模型構(gòu)建氣動熱流分布數(shù)據(jù)融合模型，獲得了較好的效果。Misfud 等采用Gappy POD 結(jié)合風(fēng)洞試驗數(shù)據(jù)和CFD 計算數(shù)據(jù)構(gòu)建了針對二維翼型和三維翼身組合體的數(shù)據(jù)融合模型，利用正則化方法提升了模型的預(yù)測精度。Perron針對不同構(gòu)型表面壓力分布場，采用POD 和流形對齊構(gòu)建了高維數(shù)據(jù)融合模型，但是建模效果有待提高。Renganathan等發(fā)展了基于貝葉斯理論框架和帶約束的本征正交分解技術(shù)，用于解決帶噪聲的不完整風(fēng)洞測量場與確定但有偏的數(shù)值模擬場的融合問題。Sun 等開發(fā)了一種快速風(fēng)場重建方法，提出一種傳感器布局優(yōu)化求解過程，基于傳感器獲得的有限測量數(shù)據(jù)，通過POD 技術(shù)對風(fēng)場進行逆過程快速重構(gòu)。Zhao 等提出了基于壓縮感知算法的針對翼型表面進行壓力分布重構(gòu)的方法。Li 等基于深度神經(jīng)網(wǎng)絡(luò)，融合風(fēng)洞試驗數(shù)據(jù)和CFD 計算數(shù)據(jù)，預(yù)測了機翼不同狀態(tài)、不同截面的壓力系數(shù)。

現(xiàn)代數(shù)值模擬規(guī)模和分辨率的提高，提供了大量仿真數(shù)據(jù)，因此本文提出了一種融合有限的風(fēng)洞試驗數(shù)據(jù)和數(shù)值計算（CFD）數(shù)據(jù)的新方法，通過本征正交分解技術(shù)對仿真得到的壓力分布數(shù)據(jù)進行特征提取，然后融合稀疏的風(fēng)洞試驗測壓數(shù)據(jù)，基于壓縮感知（Compressed Sensing，CS）算法進行分布載荷精細(xì)化重構(gòu)研究。與Misfud、Renganathan等工作的不同之處在于：本文采用稀疏的試驗測壓數(shù)據(jù)進行壓力分布重構(gòu)，通過壓縮感知算法求解基函數(shù)對應(yīng)的坐標(biāo)，提高了方法的魯棒性，對于變幾何變來流狀態(tài)算例仍具有較好的適應(yīng)性，對不同翼型具有較好的泛化性，且積分獲得的升力系數(shù)和俯仰力矩系數(shù)精度較高。

1 研究模型及方法

1.1 研究模型

算法框架如圖1所示，主要包括以下4 個方面：

圖1 算法框架流程Fig.1 Algorithm framework process

1）使用CFD 方法計算不同工況下（幾何改變、狀態(tài)改變）的壓力分布數(shù)據(jù)。

2）采用本征正交分解技術(shù)提取CFD 計算的壓力分布數(shù)據(jù)的低維特征，即POD 基函數(shù)。

3）優(yōu)化翼型表面壓力傳感器位置：采用粒子群優(yōu)化算法（PSO）獲得最佳測量點位置。

4）利用有限試驗測壓數(shù)據(jù)，基于壓縮感知算法進行壓力分布重構(gòu)，對壓力分布進行積分獲得升力系數(shù)和俯仰力矩系數(shù)，并通過試驗數(shù)據(jù)進行驗證。

1.2 研究方法

本文中算例均為定常流動，CFD 求解器基于定常雷諾平均Navier-Stokes（N-S）方程，求解方法采用二階有限體積法，空間離散格式采用Roe 格式，采用偽時間推進法，通過隱式Gauss-Seidel 迭代求解。

本征正交分解技術(shù)是一種常用的降維手段，能夠?qū)?fù)雜動力學(xué)從高維離散化系統(tǒng)投影到低維系統(tǒng)。該技術(shù)的本質(zhì)是通過對流場樣本進行矩陣變換及正交分解，得到使樣本殘差最小的若干正交基函數(shù)，用于描述流場的主要規(guī)律?；赑OD 技術(shù)，根據(jù)特征值的大小對模態(tài)按照能量進行排序，可以提取出主要的流動模態(tài)。通常采取降維方式簡化模型：按能量截取一個維數(shù)為M（M＜＜N）的低維空間（式（1））使（）≥；一般取99%（略小于1），以捕獲絕大部分能量。

其中，N 為樣本空間維度，為特征值大小，為能量大小。

粒子群優(yōu)化算法是一種進化計算技術(shù)，其基本思想是通過群體中個體之間的協(xié)作和信息共享來尋找最優(yōu)解。本文使用粒子群優(yōu)化算法獲得測壓點位置，基本步驟如下：

1）將CFD 所計算的壓力分布數(shù)據(jù)分為兩部分，一部分為訓(xùn)練集，另一部分為測試集，基函數(shù)在訓(xùn)練集中通過POD 獲取。

2）在測試集中選取一部分壓力分布數(shù)據(jù)，通過粒子群優(yōu)化算法，測量點遍歷翼型上所有點，通過壓縮感知算法重構(gòu)流場，以重構(gòu)誤差作為判斷標(biāo)準(zhǔn)。

3）輸出一組使重構(gòu)誤差較小的測量點位置，將輸出的測量點位置用于壓力分布重構(gòu)。

壓縮感知技術(shù)突破了香農(nóng)采樣定理的局限，可用較少的數(shù)據(jù)實現(xiàn)高精度的完整信號重構(gòu)。壓縮感知技術(shù)一經(jīng)提出即引起廣大研究者的興趣，被廣泛應(yīng)用于地理、航天、通信等各個領(lǐng)域。壓縮感知基本框架如圖2所示。假定觀測值y ∈R是通過測量矩陣φ ∈R建立在原信號x ∈R上得到的，即=，壓縮感知的目的是通過稀疏的觀測值恢復(fù)原信號（欠定條件下，M＜＜N）。

圖2 壓縮感知框架Fig.2 Compressed sensing framework

信號通過稀疏矩陣可以進行稀疏表示，即：

因此，方程轉(zhuǎn)化為：

可通過下面的優(yōu)化問題近乎完美地從y 中恢復(fù)出稀疏信號s，進而恢復(fù)x：

上式即L0 范數(shù)最小化問題。若直接求解式（4），則是一個N-P hard 問題。在一定條件下，L1 范數(shù)最小化與L0 范數(shù)最小化共解，于是式（4）轉(zhuǎn)化為：

本文將POD 和CS 技術(shù)相結(jié)合，原信號x 為翼型表面完整壓力分布，觀測值y 為稀疏的壓力試驗數(shù)據(jù)，算法如下：

1）通過POD 提取CFD 計算得到的壓力場數(shù)據(jù)，得到POD 基函數(shù)，即提取計算數(shù)據(jù)的低維特征。

2）對翼型表面進行壓力點采樣，即=，為測量矩陣，y 為采樣點壓力值。

3）基于壓縮感知技術(shù)進行流場重構(gòu)（L1 范數(shù)最小化）。

2 算例驗證

2.1 亞聲速流動

采用CFD 計算了132 組NACA0012 翼型亞聲速狀態(tài)的壓力分布數(shù)據(jù)（馬赫數(shù)Ma=0.1～0.6，迎角=0°～11°），采用本征正交分解技術(shù)提取了8 組CFD 計算數(shù)據(jù)的POD 基函數(shù)，采用粒子群優(yōu)化算法得到4 個測壓點位置，并基于4 個測壓點處的數(shù)據(jù)進行重構(gòu)。CFD 計算NACA0012 翼型采用的結(jié)構(gòu)化C 型網(wǎng)格如圖3所示：

圖3 NACA0012 翼型網(wǎng)格Fig.3 NACA0012 airfoil grid

為了量化說明重構(gòu)流場與實際流場的偏差程度，定義誤差如下：

其中：RMSE 為均方根誤差；C為重構(gòu)得到的壓力系數(shù)，C為試驗測得的壓力系數(shù)；n 為翼型表面試驗壓力點的個數(shù)；C為通過重構(gòu)的壓力分布積分得到的升力系數(shù)，C為試驗數(shù)據(jù)對應(yīng)的升力系數(shù)；C表示通過重構(gòu)的壓力分布積分得到的俯仰力矩系數(shù)，C表示試驗數(shù)據(jù)對應(yīng)的俯仰力矩系數(shù)；C和C分別表示升力系數(shù)和俯仰力矩系數(shù)重構(gòu)值與試驗值之間的誤差，俯仰力矩系數(shù)是對翼型1/4 弦線處進行積分的。

圖4為NACA0012 翼型在Ma=0.3、=0.106 9°、Re=6.0×10工況與Ma=0.3、10°、Re=6.0×10工況的重構(gòu)結(jié)果及其與試驗、CFD 數(shù)據(jù)的對比，兩組工況的試驗數(shù)據(jù)與重構(gòu)數(shù)據(jù)之間的RMSE 分別為 0.042 1、 0.074 5。 圖5為 NACA0012 翼型在Ma=0.504、=4.06°、Re=3.8×10工況下的重構(gòu)結(jié)果與試驗、CFD 數(shù)據(jù)的對比，試驗數(shù)據(jù)和重構(gòu)數(shù)據(jù)之間的RMSE 為0.064 1。圖中藍(lán)色虛線為CFD數(shù)據(jù)，黑色實線為重構(gòu)曲線，紅色空心三角形為試驗數(shù)據(jù)，綠色實心倒三角形為所采用的測壓點數(shù)據(jù)。表1給出了圖5對應(yīng)試驗狀態(tài)的升力系數(shù)和俯仰力矩系數(shù)以及對應(yīng)誤差，表2給出了較多試驗工況下重構(gòu)的升力系數(shù)和俯仰力矩系數(shù)誤差上下限統(tǒng)計結(jié)果。通過對比誤差以及重構(gòu)曲線，發(fā)現(xiàn)壓力分布曲線與試驗數(shù)據(jù)匹配完整，且精度較高，升力系數(shù)誤差和俯仰力矩系數(shù)誤差在0.005 以內(nèi)。

圖4 NACA0012 亞聲速重構(gòu)結(jié)果Fig.4 Subsonic reconstruction results of NACA0012 airfoil

圖5 NACA0012 亞聲速重構(gòu)結(jié)果Fig.5 Subsonic reconstruction results of NACA0012 airfoil

表1 重構(gòu)誤差Table 1 Reconstruction error

表2 誤差區(qū)間Table 2 Error interval

圖6為基于稀疏試驗數(shù)據(jù)重構(gòu)和傳統(tǒng)插值方法重構(gòu)的結(jié)果對比，其中試驗數(shù)據(jù)與圖5相同。圖6（a）為翼型表面測壓數(shù)據(jù)均勻采樣結(jié)果，圖6（b）為升力系數(shù)和俯仰力矩系數(shù)精度隨試驗點數(shù)量的變化曲線。圖中，1–Error=1–|Intergral–Exp|，Integral表示積分獲得的升力系數(shù)或俯仰力矩系數(shù)，Exp 表示升力系數(shù)或俯仰力矩系數(shù)的試驗數(shù)據(jù)；綠線表示俯仰力矩系數(shù)精度，黑線表示升力系數(shù)精度，藍(lán)色圓圈和紅色菱形分別表示基于4 個試驗點重構(gòu)獲得的壓力分布積分得到的升力系數(shù)和俯仰力矩系數(shù)精度。如圖6（b）所示， 22 個試驗點重構(gòu)誤差的原因在于：其壓力采樣點位置沒有選到吸力峰位置，導(dǎo)致積分誤差較大，說明壓力采樣點位置的選取會影響插值重構(gòu)方法的精度?？梢钥吹剑啾扔趥鹘y(tǒng)插值重構(gòu)方法，本文重構(gòu)方法以稀疏的試驗數(shù)據(jù)獲得了高精度的升力系數(shù)和俯仰力矩系數(shù)。

圖6 稀疏重構(gòu)與插值重構(gòu)結(jié)果對比Fig.6 Comparison of sparse reconstruction and interpolation reconstruction

2.2 跨聲速流動

采用CFD 計算了100 組NACA0012 翼型跨聲速狀態(tài)的壓力分布數(shù)據(jù)（Ma=0.6～0.8，=0°～4°），通過本征正交分解技術(shù)提取了14 組CFD 計算數(shù)據(jù)的POD 基函數(shù)。為了更好地捕捉激波特性，選取了14 個試驗測壓點處的數(shù)據(jù)進行重構(gòu)。

圖7為NACA0012 翼型在 Ma=0.652、=5.86°、Re=3.0×10和Ma=0.697、=4.86°、 Re=3.0×10這2種工況下的重構(gòu)結(jié)果及其與試驗和CFD 數(shù)據(jù)的對比，表3為CFD 計算、試驗和重構(gòu)得到的升力系數(shù)及對應(yīng)誤差。圖8為NACA0012 翼型在Ma=0.756、=3.99°、Re=3.76×10，Ma=0.754、=2.98°、Re=3.80×10和Ma=0.754、=0.99°、 Re=3.96×10這3 種工況的重構(gòu)結(jié)果及其與試驗和CFD 數(shù)據(jù)的對比，表4為試驗及重構(gòu)的升力系數(shù)、俯仰力矩系數(shù)及其對應(yīng)誤差?？梢钥吹剑┭隽叵禂?shù)相對誤差較大。一方面，由于俯仰力矩系數(shù)積分點位于1/4弦線處，該點與氣動中心距離較近，所以俯仰力矩系數(shù)值較小，導(dǎo)致其相對誤差較大；另一方面，俯仰力矩系數(shù)不僅與壓力分布相關(guān)，還與摩擦力分布相關(guān)，這可能導(dǎo)致相對誤差較大。表5為較多試驗工況下重構(gòu)升力系數(shù)和俯仰力矩系數(shù)絕對誤差上下限的統(tǒng)計結(jié)果。可以看出，即使CFD 計算結(jié)果與試驗數(shù)據(jù)誤差較大，但經(jīng)過重構(gòu)得到的壓力分布曲線與試驗數(shù)據(jù)吻合較好，能夠獲得與真實試驗最完整匹配的壓力分布曲線，且對激波站位捕捉較好。通過誤差對比，能夠?qū)⑸ο禂?shù)誤差控制在0.005 以內(nèi)，俯仰力矩系數(shù)誤差控制在0.006 以內(nèi)。

表4 重構(gòu)誤差Table 4 Reconstruction error

表5 誤差區(qū)間Table 5 Error interval

圖7 NACA0012 翼型跨聲速重構(gòu)結(jié)果Fig.7 Transonic reconstruction results of NACA0012 airfoil

圖8 NACA0012 翼型跨聲速重構(gòu)結(jié)果Fig.8 Transonic reconstruction results of NACA0012 airfoil

表3 重構(gòu)誤差Table 3 Reconstruction error

2.3 變幾何變來流狀態(tài)

本文采用CST 參數(shù)化方法，選用12 個CST 參數(shù)描述翼型表面幾何形狀，采用拉丁超立方抽樣，在NACA0012 翼型CST 參數(shù)的1±30%范圍內(nèi)進行擾動，隨機抽取樣本翼型（所有樣本翼型都與NACA0012 翼型不同）。采用CFD 計算了900 組壓力分布數(shù)據(jù)，來流狀態(tài)范圍為Ma=0.6～0.75、=0°～4°。采用本征正交分解技術(shù)提取了20 組CFD 計算數(shù)據(jù)的POD 基函數(shù)，選取了18 個試驗測壓點處的數(shù)據(jù)進行重構(gòu)。

為了證明本文發(fā)展的方法對于不同翼型仍然具有很好的泛化性，對3 種翼型進行測試，如圖9所示，其中黑色點畫線表示NACA0012 翼型，紅色虛線表示MBB-A3 翼型，橙色實線表示RAE2822 翼型，淺藍(lán)色線條表示通過CST 參數(shù)化隨機采樣獲得的樣本翼型。針對不同翼型表面選取測壓點，首先將其橫坐標(biāo)歸一化在[0，1]之間，然后根據(jù)之前所得18 個測壓點位置，布點獲得對應(yīng)點的壓力數(shù)據(jù)。

圖9 翼型采樣空間Fig.9 The airfoil sampling space

圖10 為MBB-A3 翼型（Ma=0.698、=5.49°、Re=6.17×10工況）、NACA0012 翼型（Ma=0.753、=1.95°、Re=3.3×10工況）和RAE2822 翼型（Ma=0.676、=2.4°、Re=5.7×10工況）的重構(gòu)結(jié)果及其與試驗數(shù)據(jù)的對比，表6為試驗和重構(gòu)的升力系數(shù)、俯仰力矩系數(shù)以及對應(yīng)誤差。通過3 種翼型測試算例可知，與試驗壓力分布曲線相比，重構(gòu)曲線存在小幅振蕩。原因可能有兩方面：一方面，由于CFD 計算樣本量太少，提取的低維特征不足以完全捕獲任意翼型在不同來流狀態(tài)下壓力分布所構(gòu)成的特征空間；另一方面，提取的低維特征中存在高階POD 基函數(shù)使其存在微小振蕩。基于重構(gòu)的壓力分布曲線積分獲得的升力系數(shù)和俯仰力矩系數(shù)仍滿足要求，該方法針對變幾何變來流狀態(tài)算例仍然適用，但需要更多的POD 基函數(shù)和測壓點數(shù)據(jù)。表7為較多試驗工況下重構(gòu)升力系數(shù)和俯仰力矩系數(shù)絕對誤差上下限的統(tǒng)計結(jié)果，通過誤差對比，能夠?qū)⑸ο禂?shù)和俯仰力矩系數(shù)誤差控制在0.005 以內(nèi)。

圖10 變幾何變來流狀態(tài)重構(gòu)結(jié)果Fig.10 Reconstruction of variable geometry with variable flow state

表6 重構(gòu)誤差Table 6 Reconstruction error

表7 誤差區(qū)間Table 7 Error interval

3 結(jié) 論

本文提出了一種融合風(fēng)洞試驗數(shù)據(jù)和數(shù)值計算數(shù)據(jù)的新方法，基于較少的試驗數(shù)據(jù)重構(gòu)獲得與真實試驗最完整匹配的壓力分布曲線。所發(fā)展的重構(gòu)方法可在一定程度上解決空間受限稀疏觀測條件下的分布載荷精細(xì)化重構(gòu)難題。主要結(jié)論如下：

1） 通過提取壓力分布CFD 計算數(shù)據(jù)的低維特征，利用有限試驗壓力測量數(shù)據(jù)獲得的升力系數(shù)和俯仰力矩系數(shù)比直接通過試驗數(shù)據(jù)積分獲得的精度更高。

2） 該方法在固定翼型亞聲速、跨聲速算例以及變幾何變來流狀態(tài)算例中得到了驗證，對不同翼型具有良好的泛化性，重構(gòu)得到的分布載荷、升力系數(shù)和俯仰力矩系數(shù)均能精確匹配試驗結(jié)果。

3） 由于跨聲速區(qū)激波非線性的影響，需要更多的基函數(shù)和測壓點來進行壓力分布重構(gòu)。

本文雖然是對壓力分布進行重構(gòu)，但提出的基于稀疏試驗數(shù)據(jù)的融合方法可應(yīng)用于剪應(yīng)力和氣動熱測量等領(lǐng)域，對于其他學(xué)科也具有應(yīng)用價值。