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

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

0 引言

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

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

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

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

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

1.1 系統(tǒng)模型

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

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

1.2 狀態(tài)空間矩陣

系統(tǒng)模型的估計與實驗測量均存在誤差，卡爾曼濾波的數(shù)據(jù)同化手段通過狀態(tài)空間模型將誤差代入系統(tǒng)當(dāng)中：

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

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

向量yt由實驗測點對應(yīng)位置的測量值構(gòu)成。在稀疏壓力重構(gòu)中，選取實測壓力p，向量yt維數(shù)與實驗測點的數(shù)量（m 個）一致：

1.3 集合變換卡爾曼濾波方法及實現(xiàn)流程

式中：H為式（3）中的投影矩陣，在翼型網(wǎng)格和實驗?zāi)Ｐ偷膶?yīng)位置設(shè)置為1，其余位置設(shè)置為0。

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

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

2 實驗過程和結(jié)果

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

2.1 RAE 2822 翼型壓力重構(gòu)

2.1.1 網(wǎng)格和計算條件設(shè)置

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

表1 RAE 2822 網(wǎng)格質(zhì)量及節(jié)點數(shù)量Table 1 Quality and node number of RAE 2822 mesh

圖2 RAE 2822 翼型網(wǎng)格Fig.2 Mesh of RAE 2822 airfoil

2.1.2 RAE 2822 翼型ETKF 數(shù)據(jù)同化結(jié)果

2.1.2.1 迎角及馬赫數(shù)

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

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

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

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

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

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

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

2.1.2.2 壓力系數(shù)、升力系數(shù)和力矩系數(shù)

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

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

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

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

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

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

表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 翼型壓力重構(gòu)

2.2.1 網(wǎng)格和計算條件設(shè)置

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

表5 NACA0012 網(wǎng)格質(zhì)量及節(jié)點數(shù)量Table 5 Quality and node number of NACA 0012 mesh

圖9 NACA 0012 翼型網(wǎng)格Fig.9 Mesh of NACA 0012 airfoil

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

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

2.2.2 NACA 0012 翼型ETKF 數(shù)據(jù)同化結(jié)果

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

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

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

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

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

3 結(jié)論

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

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

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

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

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