基于多源數(shù)據(jù)融合的翼型表面壓強精細化重構方法

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

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

0 引 言

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

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

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

1 研究模型及方法

1.1 研究模型

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

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

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

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

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

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

1.2 研究方法

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

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

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

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

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

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

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

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

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

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

因此，方程轉化為：

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

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

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

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

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

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

2 算例驗證

2.1 亞聲速流動

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

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

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

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

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

圖4 NACA0012 亞聲速重構結果Fig.4 Subsonic reconstruction results of NACA0012 airfoil

圖5 NACA0012 亞聲速重構結果Fig.5 Subsonic reconstruction results of NACA0012 airfoil

表1 重構誤差Table 1 Reconstruction error

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

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

圖6 稀疏重構與插值重構結果對比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°），通過本征正交分解技術提取了14 組CFD 計算數(shù)據(jù)的POD 基函數(shù)。為了更好地捕捉激波特性，選取了14 個試驗測壓點處的數(shù)據(jù)進行重構。

圖7為NACA0012 翼型在 Ma=0.652、=5.86°、Re=3.0×10和Ma=0.697、=4.86°、 Re=3.0×10這2種工況下的重構結果及其與試驗和CFD 數(shù)據(jù)的對比，表3為CFD 計算、試驗和重構得到的升力系數(shù)及對應誤差。圖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 種工況的重構結果及其與試驗和CFD 數(shù)據(jù)的對比，表4為試驗及重構的升力系數(shù)、俯仰力矩系數(shù)及其對應誤差?？梢钥吹?，俯仰力矩系數(shù)相對誤差較大。一方面，由于俯仰力矩系數(shù)積分點位于1/4弦線處，該點與氣動中心距離較近，所以俯仰力矩系數(shù)值較小，導致其相對誤差較大；另一方面，俯仰力矩系數(shù)不僅與壓力分布相關，還與摩擦力分布相關，這可能導致相對誤差較大。表5為較多試驗工況下重構升力系數(shù)和俯仰力矩系數(shù)絕對誤差上下限的統(tǒng)計結果?？梢钥闯?，即使CFD 計算結果與試驗數(shù)據(jù)誤差較大，但經(jīng)過重構得到的壓力分布曲線與試驗數(shù)據(jù)吻合較好，能夠獲得與真實試驗最完整匹配的壓力分布曲線，且對激波站位捕捉較好。通過誤差對比，能夠將升力系數(shù)誤差控制在0.005 以內(nèi)，俯仰力矩系數(shù)誤差控制在0.006 以內(nèi)。

表4 重構誤差Table 4 Reconstruction error

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

圖7 NACA0012 翼型跨聲速重構結果Fig.7 Transonic reconstruction results of NACA0012 airfoil

圖8 NACA0012 翼型跨聲速重構結果Fig.8 Transonic reconstruction results of NACA0012 airfoil

表3 重構誤差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°。采用本征正交分解技術提取了20 組CFD 計算數(shù)據(jù)的POD 基函數(shù)，選取了18 個試驗測壓點處的數(shù)據(jù)進行重構。

為了證明本文發(fā)展的方法對于不同翼型仍然具有很好的泛化性，對3 種翼型進行測試，如圖9所示，其中黑色點畫線表示NACA0012 翼型，紅色虛線表示MBB-A3 翼型，橙色實線表示RAE2822 翼型，淺藍色線條表示通過CST 參數(shù)化隨機采樣獲得的樣本翼型。針對不同翼型表面選取測壓點，首先將其橫坐標歸一化在[0，1]之間，然后根據(jù)之前所得18 個測壓點位置，布點獲得對應點的壓力數(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工況）的重構結果及其與試驗數(shù)據(jù)的對比，表6為試驗和重構的升力系數(shù)、俯仰力矩系數(shù)以及對應誤差。通過3 種翼型測試算例可知，與試驗壓力分布曲線相比，重構曲線存在小幅振蕩。原因可能有兩方面：一方面，由于CFD 計算樣本量太少，提取的低維特征不足以完全捕獲任意翼型在不同來流狀態(tài)下壓力分布所構成的特征空間；另一方面，提取的低維特征中存在高階POD 基函數(shù)使其存在微小振蕩?；谥貥嫷膲毫Ψ植记€積分獲得的升力系數(shù)和俯仰力矩系數(shù)仍滿足要求，該方法針對變幾何變來流狀態(tài)算例仍然適用，但需要更多的POD 基函數(shù)和測壓點數(shù)據(jù)。表7為較多試驗工況下重構升力系數(shù)和俯仰力矩系數(shù)絕對誤差上下限的統(tǒng)計結果，通過誤差對比，能夠將升力系數(shù)和俯仰力矩系數(shù)誤差控制在0.005 以內(nèi)。

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

表6 重構誤差Table 6 Reconstruction error

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

3 結 論

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

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

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

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

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