基于數據同化的氣動壓力稀疏重構方法

黃俊，郭雨欣，冀晶晶，黃永安

華中科技大學 機械科學與工程學院，武漢 430074

0 引言

風洞作為空氣動力學研究的重要地面實驗設備，在先進飛行器研制和基礎空氣動力學問題研究中發揮了不可替代的作用。實驗流體力學（EFD）著眼于風洞實驗數據的采集、分析和處理，是風洞系統工程中不可分割的部分。作為風洞技術測量目標之一，表面壓力的測量直接關系到飛行器的升力系數與姿態控制，同時還可以通過壓力分布判斷飛行器表面的分離、轉捩等信息?，F有的壓敏材料（PSP）利用非接觸方式獲得連續、大范圍的模型表面壓力分布[1]，是風洞實驗中有效的測壓方法，然而，基于光學原理的PSP 存在觀測死角（如飛機短艙內表面）問題。傳統的在模型表面布設測壓孔的方法技術成熟度高、測量準確性好、可按需布置[2]，最新發展的柔性智能蒙皮可以粘貼于模型內外表面，實現多物理量同步測量[3-4]。但測壓孔和柔性傳感陣列往往數量有限，需發展基于稀疏測點的場重構技術，實現壁面繞流及空間流場的感知。

流場數據具有多維度且復雜的特點，而測量值來自壁面稀疏測點，如何構建有限數據與全空間流場的映射關系是流場重構的難點所在。國內外學者從流場物理模型與數據驅動的角度進行了理論嘗試，Vamsi Krishna 等[5]基于快速畸變理論和Taylor假設建立雙向加權模型，基于先驗速度場對湍流演化進行重構。Callaham 等[6]探索了基于二維流場稀疏表示和數據庫的流場重構法，嘗試了數據驅動的流場重構技術，并討論了需要的隨機測點最小數量。Sun 等[7]結合Bayesian 神經網絡與先驗物理模型重構了二維鐘形與T 形流場，并對結果不確定性進行估計。李靜等[8]采用本征正交分解方法，以較少階模態高精度再現非定常圓柱繞流的完整流場。這些重構方法大多局限在簡單幾何和穩態流場重構，與實際工程應用之間仍存在一定距離。得益于計算流體力學（CFD）的發展，通過數值方法求解Navier–Stokes 方程[9]，能夠得到豐富的流場先驗知識，為全流場重構提供良好的條件。然而，CFD 的計算過程通常不能充分考慮風洞實際流動過程中的不確定性[10]，導致確定性的計算機仿真難以復現風洞實驗對物理參數的影響，從而無法實現與風洞實測數據的匹配和融合。因此，結合CFD 與實測稀疏數據對風洞的真實流場進行重構極富挑戰。

數據同化（Data Assimilation，DA）方法能夠結合估計值和測量值進行預測[11]，可以作為聯系CFD理論先驗信息和EFD 風洞實測數據的橋梁，其在流場重構中已有發展。Chandramouli 等[12]提出一種變分數據同化方法，通過圓柱繞流在過渡狀態下的三維湍流尾流案例進行方案驗證，與PIV 實驗數據對比，驗證了變分數據同化在不可壓縮湍流重構中的有效性。Belligoli 等[13]通過數據同化方法實現了二維翼型的迎角及馬赫數的修正，降低了實驗測量值與RANS 模擬計算值之間的誤差。集合變換卡爾曼濾波（Ensemble Transform Kalman Filter，ETKF）是卡爾曼濾波（Kalman Filter，KF）針對非線性系統模型的擴展方法[14]，在復雜系統中有良好的數據同化效果。ETKF 以蒙特卡洛方法為依托，以貝葉斯原理為核心，利用稀疏測量數據所包含的信息對物理場先驗預報數據進行濾波，借此給出對不確定參數的連續后驗估計[15]。該方法源自地球物理領域，是基于海面離散探測數據對海洋與大氣的真實演化狀態進行連續修正的重要工具[16]。近年來，國內外學者將該方法引入空氣動力學領域，美國加州大學通過無限長薄板周圍隨機擾動無黏渦流模擬，印證了ETKF 在流場數據同化中的有效性[17]。日本宇航研究開發機構（JAXA）應用ETKF 對風洞中存在不確定性的馬赫數、迎角、湍流黏度進行估計[18]，成功重構剛性機翼繞流場，證實了該方法對風洞中復雜流場的估計能力。國內上海交通大學在預測直升機轉子三維流場特征時，采用了ETKF 優化剪切應力傳輸模型常數，為修正逆壓梯度下邊界層流動分離提供了參考[19]。

本文面向風洞實驗的壓力測量應用，使用集合變換卡爾曼濾波方法，以二維翼型RAE 2822 和二維對稱翼型NACA 0012 為研究對象，進行機翼表面氣動壓力重構，以達到風洞實驗對全域數據重構精度的要求。先通過CFD 計算得到全域先驗分布，結合機翼表面有限數量的壓力測量值，對機翼的迎角及馬赫數進行修正，重構得到高精度的機翼周圍全域壓力場。再采用ETKF 進行流場數據同化：一方面，可以充分利用稀疏測量數據，將這些高精度數據泛化到整個流域，使得最終展示的流場更加接近真實的流場，為空間全域感知提供可能；另一方面，ETKF 可以作為一種流場風洞干擾的修正方法，基于真實測量數據，直接進行流場分布的修正，獲取高精度的流場分布。

1 集合變換卡爾曼濾波方法

1.1 系統模型

在流場計算中，從邊界條件到物理量的分布不是簡單的線性系統，而是離散的非線性系統，其中狀態變量估計可通過求解Navier–Stokes 方程得到。本文通過軟件FLUENT 對流場進行數值計算，選用Spalart–Allmaras（S–A）湍流模型。采用Navier–Stokes 方程作為流動控制方程，并將此方程作為卡爾曼濾波中的系統模型：

式中：W 為守恒狀態矢量，包含了密度、速度和能量；V 為控制體體積；Fc為對流通量；Fv為黏性通量；S 為控制體表面積；n 為控制體外法線方向的單位矢量；τ為時間。

1.2 狀態空間矩陣

系統模型的估計與實驗測量均存在誤差，卡爾曼濾波的數據同化手段通過狀態空間模型將誤差代入系統當中：

式中：下標t 表示迭代次數，xt、yt向量分別表示系統模型和實驗測量數據的狀態向量，vt、ωt向量分別表示系統模型的噪聲和觀測噪聲，假定其符合高斯分布。非線性算子F 是從第t-1次迭代到第t 次迭代的映射，在數據同化中，由系統模型計算得到。H矩陣是將系統估計的數據矩陣投影到實驗觀測數據的投影矩陣。

向量xt包含了迎角α、馬赫數Ma 和所有（n 個）計算網格節點上的密度ρ、笛卡爾速度分量u、v 和壓力p，其維數l=4n+2，向量xt的表達式為：

向量yt由實驗測點對應位置的測量值構成。在稀疏壓力重構中，選取實測壓力p，向量yt維數與實驗測點的數量（m 個）一致：

1.3 集合變換卡爾曼濾波方法及實現流程

式中：H為式（3）中的投影矩陣，在翼型網格和實驗模型的對應位置設置為1，其余位置設置為0。

圖1 ETKF 流程圖Fig.1 The flowchart of ETKF

在各物理量符合收斂條件判定準則之后，將最后一次迭代的集合元素均值作為整個數據同化過程的來流條件修正值，代入FLUENT 進行重新計算，得到修正后的流場分布。

2 實驗過程和結果

RAE 2822 是一個二維跨聲速湍流流動的經典翼型，被許多國外的項目合作組（如EUROVAL）選作經典算例[21]；NACA 0012 作為經典的對稱翼型，有實驗數據可用來驗證其數值模擬的準確性，故本文實驗使用二維翼型RAE 2822 和NACA 0012 進行驗證。RAE 2822 翼型實驗主要用來展現重構的收斂過程，ETKF 對迎角、馬赫數的影響，以及與其他方法在精度上的對比；NACA 0012 翼型實驗則面向風洞測量中基于有限測點重構的需求，用以探究影響ETKF 重構精度的因素。

2.1 RAE 2822 翼型壓力重構

2.1.1 網格和計算條件設置

使用ICEM 軟件對RAE 2822 進行網格劃分（圖2），采用四邊形結構化網格，機翼模型弦長c=1 m，計算域在弦長方向（即x 方向）為40 m，弦長垂直方向（即y 方向）為40 m，網格質量如表1所示。本文笛卡爾坐標系的原點位于機翼最前緣點。將FLUENT 軟件的計算模型設置為S–A 湍流模型，此模型對跨聲速流動的求解有較好的效果，且計算成本相對較低。自由流空氣設置為可壓縮空氣，計算方式為穩態計算。實驗數據來源為Cook 等[22]的案例6 中RAE 2822 跨聲速實驗壓力數據，其實驗工況為：Ma=0.725，α=2.92°，Re=6.2×106。共有103 個測量點，其中一部分測點與機翼網格點難以匹配，會給數據同化帶來較大誤差，將其去除后剩余的壓力測點個數為75（即m=75）。

表1 RAE 2822 網格質量及節點數量Table 1 Quality and node number of RAE 2822 mesh

圖2 RAE 2822 翼型網格Fig.2 Mesh of RAE 2822 airfoil

2.1.2 RAE 2822 翼型ETKF 數據同化結果

2.1.2.1 迎角及馬赫數

圖3 集合成員迎角在ETKF 同化前后的分布Fig.3 Distribution of ensemble angel of attack before and after ETKF assimilation

圖4 集合成員馬赫數在ETKF 前后的分布Fig.4 Distribution of ensemble Mach number before and after ETKF assimilation

圖5 集合的迎角及馬赫數均值隨迭代過程的變化Fig.5 Mean of angle of attack and Mach number changes with iterations

圖6 為迭代過程中所有網格點的密度（ρ）、速度分量（u、v）和壓力系數（Cp）的均方誤差MSE 變化曲線。從圖中可以看出，在進行了第3 次迭代之后，所有物理量的均方誤差已經小于第一次迭代的1%，這表明在第3 次迭代時ETKF 就已經滿足收斂要求，繼續迭代發現均方誤差可以進一步降低，于是迭代繼續進行。

圖6 ρ、u、v、Cp 的MSE 隨迭代過程的變化Fig.6 MSE of ρ,u,v and Cp changes with iterations

隨著迭代的進行，集合中所有成員的邊界條件都在產生變化并且更加集中（圖3 和4）。計算集合成員的方差如表2所示，可以看到在數據同化過程中，隨機抽樣產生的所有集合成員在向更為準確的修正值變化，且馬赫數的集中效果要比迎角的集中效果更好。最終，集合成員迎角及馬赫數均值在第10 次迭代之后分別收斂為2.434°和0.732 8（圖5）。

表2 ETKF 前后集合成員迎角、馬赫數的均值及方差Table 2 Mean and variance of ensemble angle of attack and Mach number before and after ETKF

2.1.2.2 壓力系數、升力系數和力矩系數

基于圖2 網格，將收斂后集合成員的迎角及馬赫數均值（表3）作為邊界條件代入FLUENT 進行重新計算，得到ETKF 同化后的壓力系數Cp曲線分布，此結果即為ETKF 數據同化的最終結果。將ETKF 同化結果與線性理論修正值進行對比，結果如圖7 和8所示。

表3 RAE 2822 case 6 邊界條件Table 3 Boundary condition of RAE 2822 case 6

圖7 ETKF 同化后的壓力系數曲線與線性理論修正曲線對比Fig.7 Comparison of the pressure coefficient curve after assimilation of ETKF with the linear theory correction

圖8 激波位置的壓力系數對比Fig.8 Comparison of pressure coefficient at shock position

從圖7 和8 中可以看出：與線性理論修正后的計算結果相比，ETKF 同化后的壓力系數曲線更加靠近實驗數據；線性理論的迎角修正值過大，導致其計算誤差與ETKF 相比更大，尤其是在激波位置。從數值上分析（式（25））可得，ETKF 同化的壓力系數平均相對誤差e 比線性理論修正后降低了約3%。e 的計算方法如下：

式中：j 為實驗測點的數量，Cp,exp為實驗測點的壓力系數。

表4 ETKF 同化和線性理論修正后CL 和Cm、激波位置Cp 平均相對誤差與實驗值的對比Table 4 Comparison of CL,Cm and Cp average relative error near the shock wave position between ETKF,linear theory and experiment

2.2 NACA 0012 翼型壓力重構

2.2.1 網格和計算條件設置

如圖9 和表5所示，NACA 0012 翼型采用的計算網格為四邊形結構化網格，機翼模型弦長c=0.152 4 m，計算域在弦長方向為4.5 m，在垂直弦長方向為5 m，網格質量高，FLUENT 計算時選用S–A湍流模型。實驗數據來源于NASA 進行的NACA 0012 翼型風洞實驗[23]，選取TEST 119 的壓力數據進行數據同化的測試，其實驗工況為：Ma=0.402 2，α=4.020 8°，Re=6.070 6×106。共有46 個壓力測量點，由于靠近機翼前緣位置的數據在同化時難以與網格節點進行匹配，容易帶來較大誤差，故在同化時將其去除，僅使用其中44 個數據（即m≤44）。

表5 NACA0012 網格質量及節點數量Table 5 Quality and node number of NACA 0012 mesh

圖9 NACA 0012 翼型網格Fig.9 Mesh of NACA 0012 airfoil

NACA 0012 翼型實驗主要面向風洞測量中基于有限測點重構的應用需求，對比分析基于少數測點進行壓力稀疏重構的可行性。圖10 給出了實驗壓力測點的位置及序號。實驗共進行4 組，其中3 組實驗分別選用上翼面（序號1～22）的不同位置6 個實驗測點進行數據同化，第4 組實驗作為同化結果的對標，使用所有數據測點進行同化處理。

圖10 NACA 0012 實驗測點位置及序號Fig.10 Test point location and serial number of NACA 0012

2.2.2 NACA 0012 翼型ETKF 數據同化結果

在數據同化之前進行初始集合的抽樣，因為此實驗無具體參考修正值，故基于實驗原始工況進行初始值的抽樣，并設定較大的抽樣范圍。設置集合成員數量k=40，在3.52～4.52 范圍進行迎角抽樣，在0.372 2～0.432 2 范圍進行馬赫數抽樣，得到初始集合的迎角均值為4.020 8°，馬赫數均值為0.402 1。采用第1 節的理論與流程進行數據同化，同化實驗1～4 均從此初始值開始分別進行迭代，以對比使用不同實驗測量數據得到的迭代效果。4 組實驗所用的壓力測點及迭代至收斂所需的迭代次數見表6。

表6 4 組實驗的壓力測點序號及迭代次數Table 6 The number of pressure measuring points and iteration times of 4 groups of experiments

實驗結果如表7 和圖11所示，表7 中壓力系數平均相對誤差由式（25）計算得到，4 組同化實驗均使用上翼面實驗測量值進行誤差的計算（即j=22）。如圖11所示，相比于深綠色點劃線表示的未同化實驗數據，經過ETKF 同化的壓力系數均有一定的優化效果。同化實驗4 使用了所有實驗測量點進行同化，其結果的誤差最小，只有1.05%，說明ETKF 方法的同化精度與使用的數據量有關，測點數據越多，得到的同化效果越好；只采用6 個測點進行同化的實驗1～3 也達到了小于6%的平均相對誤差，證實了基于有限測點進行壓力稀疏重構、實現全域感知的可行性。且對比同化實驗1～3 的結果可以發現，重構精度與初始測點位置的選取密切相關，初始測點位置不同，相對誤差從2.42%到5.74%不盡相同，這說明在風洞實驗時可以基于重構算法進行傳感測點位置的優化布置。

表7 4 組實驗ETKF 后的迎角、馬赫數及壓力系數平均相對誤差Table 7 Angle of attack,Mach number and pressure coefficient average relative error after ETKF in 4 groups of experiments

圖11 NACA 0012 的ETKF 同化結果與實驗值的對比Fig.11 Comparison of ETKF assimilation results with experimental values in NACA 0012

3 結論

本文針對風洞測量中有限測點進行氣動壓力重構開展研究，使用集合變換卡爾曼濾波方法對二維翼型RAE 2822 和NACA 0012 的湍流流場進行數據同化，結合稀疏的實驗測量壓力數據，對兩種翼型的實驗迎角及馬赫數進行修正，通過修正后的邊界條件進行重新計算得到了全域流場。將ETKF 同化后的機翼表面壓力系數及機翼的升力系數、力矩系數與修正前數據和線性理論修正值作對比，得到以下結論：

1）ETKF 可以通過有限的測量信息進行流場的全域重構，其重構的精度與測量點的選取有直接關系，在實驗測量位置進行壓力系數的修正，其精度比經典線性理論更加靠近實驗測量值。

2）使用ETKF 修正后的迎角及馬赫數計算得到的機翼升力系數及力矩系數與實驗測量值之間的誤差很小，這表明了基于ETKF 數據同化的有效性。

3）由于ETKF 進行同化的數據基礎是實驗測量數據，其對于流場的估計精度取決于此測量值的準確性。

本文研究表明，作為一種用于湍流流場的數據同化方法，集合變換卡爾曼濾波能夠有效利用CFD數值計算數據和EFD 高精度稀疏測量數據，重構得到更高精度的全域流場，適用于湍流場穩態計算的分布預測。本文在給定了實驗迎角及馬赫數的條件下進行了壓力場稀疏重構，后續可以進一步進行其他分布條件下（如剪應力分布）的數據同化實驗，進行多數據融合下的流場重構，推進風洞實驗測量技術向高效化發展。