范晨麟， 李孝偉

(上海大學上海市應用數(shù)學和力學研究所，上海200072)

飛機起飛和著陸時的穩(wěn)定性和安全性問題已越來越受到人們的重視，所以，有關增升裝置的研究一直都是飛機設計中的重要環(huán)節(jié).多段翼型是增升裝置的二維構型，研究其流動結構及氣動性能對于提高增升裝置的設計水平有著十分重要的意義.

流動結構的復雜特征使得對多段翼型進行實驗研究和理論研究都有很大困難，所以，近年來數(shù)值模擬方法已成為研究多段翼型流動的首要手段.但是，在流動狀態(tài)參量如迎角、馬赫數(shù)等發(fā)生變化時，通過數(shù)值模擬將會得到海量的流場數(shù)據(jù).如何從這些海量的數(shù)據(jù)中獲得流動的主要信息，甚至從已知的信息出發(fā)去預測一些新的信息，已成為數(shù)值處理研究中的一個重要課題.

本征正交分解方法(POD)，又稱為Karhunen-Loeve展開，是新近發(fā)展起來的一種處理海量數(shù)據(jù)的有效手段［1］.它提供了一種獨特的描述隨機場的方法，即將隨機場寫為基函數(shù)——本征模態(tài)的級數(shù)展開式，而這些本征模態(tài)取決于隨機場本身.由于這些本征模態(tài)在均方意義上是統(tǒng)計最優(yōu)的，因而，在展開式中只需要少量的項數(shù)即可較準確地描述隨機場，因此它又被作為數(shù)值模型的一種有力的降解和簡化工具.POD方法中的本征模態(tài)是按照對應本征值所含流場能量的大小進行排序的，因此，通過捕捉高能量本征模態(tài)能很好地描述和分析隨機場的特性;同時，由于能夠得到在狀態(tài)參量變化的有效范圍內(nèi)的所有流場信息，因此，POD方法兼具一定的預測功能.基于此，POD方法被廣泛地應用于流體力學的研究中［2］.

在復雜外形的流場計算中，為了使計算網(wǎng)格的生成更簡單，同時能保持高質(zhì)量，一般都會采用多塊分區(qū)網(wǎng)格技術，如本工作中的多段翼型流動計算就是采用多塊嵌套網(wǎng)格技術.運用嵌套網(wǎng)格技術，意味著會存在“挖洞”建立內(nèi)邊界的過程，此時，會留下一些沒有實際意義且不參與流場計算的“洞中點”.如果用POD方法來處理采用嵌套網(wǎng)格技術得到的數(shù)值結果，首先要面對的問題就是如何處理“洞中點”上攜帶的信息.如果對這些點都采用POD方法進行處理分析，則會在一定程度上影響流場的重構和模態(tài)分析的準確性，所以需要對這些點進行特殊處理.

本工作首先采用數(shù)值模擬技術對迎角變換過程中的多段翼型的流場進行計算.在網(wǎng)格的布局方面采用嵌套網(wǎng)格技術，在得到一系列的數(shù)值解以后，運用POD方法對計算數(shù)據(jù)進行后處理，以幫助理解流場的主要結構和信息.在處理“洞中點”時，主要運用預處理方法對無意義點進行篩選和剔除，避免這些點對整個POD后處理分析的干擾和影響.

1 控制方程

流動控制方程［3］為Navier-Stokes(N-S)方程：

式中，

式中，ρ，(u，v)，E，H，p，T分別為流體的密度、速度q在絕對坐標系下的2個分量、總能、總焓、壓強、溫度，q，qb分別為流體質(zhì)點的絕對速度、控制體表面的絕對速度，Ix，Iy分別為慣性系的沿坐標軸方向的單位矢量.

本工作采用中心有限體積格式進行空間離散［4］，湍流模型采用Baldwin-Lomax(B-L)模型.

2 POD方法

POD方法的一個重要特點就是可以將相干結構及其包含的能量聯(lián)系起來，即從能量的角度對速度場進行分析，并識別流場中含有最大能量的模態(tài).POD方法的核心思想［5］就是通過計算得到函數(shù)空間{vn(x)∈Ω}的一組“最優(yōu)”正交基{φ1，φ2，…，φ∞}，這里的“最優(yōu)”意味著使函數(shù)投影到正交基上產(chǎn)生的誤差達到最小，這一目的等同于以下最優(yōu)問題：

式中，{Ω}為函數(shù)空間的定義域，該空間的內(nèi)積為(·，·)，〈·〉表示對集合中n個元素進行算術平均.可以認為集合{vn(x)∈Ω}是不同次實驗采集的矢量數(shù)據(jù)集，也可以認為它是在不同時刻采集的數(shù)據(jù)集合.用變分方法求解這個最優(yōu)問題，可以得到如下特征值問題：

式中，Ω為計算流場域，

其中R(x，x')為半正定的關系矩陣.同時，由于φi(x)為正交基，所以

而集合{Vn∈Ω}中的元素都可以由下式表示：

由式(6)，可得式(7)中的系數(shù)ai(t)為

關于POD理論的更多介紹請參見文獻［1］.

由于本工作中處理的數(shù)據(jù)集合點數(shù)較為龐大，如果直接求解特征問題(4)，其求解過程會非常困難.為了解決這個問題，本工作采用Sirovich等［6］提出的快照技術(snapshots)［1］，因為該方法能夠在很大程度上減少計算對于內(nèi)存的消耗，即采用原函數(shù)空間元素vn的線性組合來表示特征模態(tài)，即

同時，對R(x，x')進行離散，可得

將式(9)和(10)代入式(4)，可得

其中當N足夠大時，式(11)可化簡為

繼續(xù)簡化式(12)，可得

最后，可以得到

即

3 嵌套網(wǎng)格與數(shù)據(jù)篩選和還原

嵌套網(wǎng)格［7］主要包括2個部分：①將計算域分成多個互相有重疊部分的子域，人為地給定重疊區(qū)域的內(nèi)、外邊界;②建立各子域間流場信息的傳遞.

在本工作中，襟翼網(wǎng)格的物面邊界條件包括物面無滑移條件、物面絕熱壁條件和物面的法向梯度為0.嵌合埋入邊界上的流動參數(shù)由子域之間的相互干擾來確定，遠場邊界上每個網(wǎng)格的中心值都由背景網(wǎng)格的流場值插值得到.為實現(xiàn)兩段翼型流場間的信息交流，需要在主段翼型流場中給出一個包圍襟翼的內(nèi)邊界即嵌合埋入邊界，該內(nèi)邊界將襟翼流場的信息傳遞給主段翼型流場.為此，須定義一個“洞”，該洞的邊界圍繞襟翼且存在于襟翼網(wǎng)格中.如果主段翼型流場中的點處于該洞中，則稱其為“洞中點”;如果主段翼型流場中的點與洞相毗鄰且不在洞中，則稱其為內(nèi)邊界點或插值點，這樣的點就形成了主段翼型流場的一個內(nèi)邊界.洞中點不參與流場的計算，而插值點上的流場參數(shù)在計算中不變，且在求解主翼流場之前由襟翼流場插值計算得到.另外，在計算襟翼流場之前，還需通過主翼流場插值計算確定其外邊界上的流場參數(shù).

在數(shù)據(jù)篩選和還原方面，主翼和襟翼均采用靜態(tài)的嵌套網(wǎng)格［8］系統(tǒng)，即主翼網(wǎng)格和襟翼網(wǎng)格相對靜止且對坐標系絕對靜止.由于嵌套網(wǎng)格系統(tǒng)的特點，計算時會產(chǎn)生一部分無意義點——洞邊界內(nèi)部和物面內(nèi)部的點，而這些無意義點若參與snapshots方法后處理分析，則會對正確描述流動結構產(chǎn)生干擾作用.因此，在網(wǎng)格生成過程中，首先，將洞邊界內(nèi)部和物面內(nèi)部的無意義點進行標記;然后，在應用snapshots方法處理之前，對每個變化時段里的無意義的標記點進行預處理篩選并剔除在計算之外;最后，在通過snapshots方法處理完成有效數(shù)據(jù)之后，再將洞邊界內(nèi)部和物面內(nèi)部的無意義點的數(shù)據(jù)還原至原坐標位置.這樣既保證了參與后處理數(shù)據(jù)的有效性，又保證了分析和對比流場時對應數(shù)據(jù)的完整性.

4 算例與分析

4.1 方法驗證

為了驗證本工作snapshots方法的適用性，對曲面方程進行數(shù)值重構模擬，設

定義空間步長x為均分的25個點，時間步長t為均分的50個點，函數(shù)z的原型和重構如圖1所示.首先，構造相關矩陣Z=z×z'，通過POD方法分解關系矩陣Z，分別對Z進行前三階的重構擬合.圖2為本征模態(tài)對應本征值的排列.綜上可得，前三階的模態(tài)已經(jīng)包含了全部能量的99%左右.由圖1(d)可以很明顯地看出，前三階模態(tài)的重構已經(jīng)能很好地擬合曲面函數(shù)，同時，和文獻［9］中的方法相比較可證明本工作snapshots方法的有效性.

圖1 原函數(shù)與重構擬合Fig.1 Actual surface and approximation

4.2 迎角變化的多段翼型襟翼的流動數(shù)值模擬和POD分析

圖2 模態(tài)能量值Fig.2 Model value

本研究對迎角變化的多段翼型流場進行了數(shù)值模擬，選用的多段翼型為帶30%襟翼的GAW-2翼型結構.采用的嵌套網(wǎng)格布局如圖3所示，對主翼和襟翼分別產(chǎn)生一個C形網(wǎng)格.同時為了保證計算精度，對流動比較劇烈的襟翼倉附近的網(wǎng)格進行橫向加密，在2個網(wǎng)格中都建立了人工洞邊界.在縫道流動區(qū)內(nèi)都進行縱向局部加密，這樣既可以滿足縫道流動的粘性特點，又能保證網(wǎng)格嵌套時重疊區(qū)有足夠的網(wǎng)格層數(shù).算例的來流馬赫數(shù)Ma=0.3，雷諾數(shù)Re=2.2×106，迎角起始角為α=－2°.迎角的變換范圍為－2°≤α≤10.8°，變換間隔為Δα=0.2°，共取65個瞬態(tài)流場，并且對其進行snapshots分析.圖4為迎角變換至10.8°時的瞬態(tài)流場，而圖5是u，v對應的65階本征模態(tài)各本征值占全部動能的比例.在snapshots計算結果中，占支配地位的本征模態(tài)通常代表流場中具有較高動能的大尺度湍流結構.從圖中可以看出，在翼型周圍的流場中，一階模態(tài)本征值幾乎占總動能的87%-u和77%-v，而二階模態(tài)占總動能的大約8%-u和19%-v，三階模態(tài)占總動能的大約2%-u和2%-v.由此可見，高階模態(tài)對總動能的貢獻相對較低，因此，可以認為低階模態(tài)是整個流動結構的載體，而較高階模態(tài)則包含了小尺度的復雜流動結構的信息.

圖3 嵌套C形網(wǎng)格Fig.3 Overlapped C-grid

圖4 原流場Fig.4 Original flow field

圖5 模態(tài)能量百分比Fig.5 Modal value percentage

圖6 三階重構流場Fig.6 Rank 3 approximation of flow field

由于u，v的前三階模態(tài)分別包含了總能量的97%和98%，因此，由前三階模態(tài)重構的流場已經(jīng)涵蓋了整個流場的絕大多數(shù)的能量，與原流場擬合得很好(見圖6).圖7為前三階本征模態(tài)的速度場，每個模態(tài)代表包含在流動中的不同尺度的流動結構：①一階模態(tài)代表流動中包含了絕大部分能量的結構，它的流動結構就是平均流，同時在每個模態(tài)中，主翼和襟翼縫隙間都出現(xiàn)了渦旋結構，因此，同平均流一樣，在主翼和襟翼縫隙間出現(xiàn)的渦旋在整個變化過程中對流場的影響起主要作用;②二階模態(tài)中主翼的頭部、尾部和襟翼的下方有3條明顯的分離線，而在主翼上方出現(xiàn)了明顯的順時針脫體渦，同時，由于二階模態(tài)中v方向上占總能量比重相對要大，所以速度場受v的影響要大，因此，由圖7(b)可以發(fā)現(xiàn)，整個流動結構呈向上環(huán)流形態(tài);③三階模態(tài)同樣在主翼的前上方和尾部上方有2條較為明顯的分離線，在分離線之間形成逆向的回流，且回流結構位置與二階模態(tài)中出現(xiàn)的脫體渦位置幾乎一致，但流動方向相反.而渦旋結構則主要出現(xiàn)在主翼和襟翼的尾部，由圖7(c)可以看到2個明顯的尾渦.圖8為u，v方向的前三階模態(tài)權重在整個變化過程的分布規(guī)律.可以發(fā)現(xiàn)，低階模態(tài)所含能量與迎角變化有弱相關性，而隨著模態(tài)階數(shù)的上升，高階模態(tài)所占能量對于迎角變化十分敏感，即不同的小尺度復雜脈動結構對于整個流場的影響會隨著迎角的變化而變化，而大尺度的流動結構的影響隨迎角的變化相對不明顯.

圖7 前三階模態(tài)Fig.7 Modals of the first three rank

圖8 模態(tài)權重Fig.8 Modal weight

5 結論

本工作通過嵌套網(wǎng)格技術對迎角變化情況下的多段翼型主翼和襟翼周圍的流場進行了數(shù)值模擬，并對連續(xù)采集的瞬態(tài)速度場進行篩選處理;最后，將具備研究意義的數(shù)據(jù)點進行POD后處理統(tǒng)計分析.通過研究，得到如下結論.

(1)POD方法不同，內(nèi)積方式的選擇對結果分析不會產(chǎn)生本質(zhì)的影響，但對結果的收斂速度和模態(tài)分析會產(chǎn)生一定的影響.本工作采用了將u，v分離的能量內(nèi)積方法對多段翼型主翼和襟翼周圍的速度場進行POD分析，相對減緩能量譜的收斂速度，以確保低階模態(tài)具有較為簡單的結構，這樣更有利于對流場的分析.對于POD方法而言，內(nèi)積方式的選取不存在單一的準則，可以根據(jù)研究需要進行選取.

(2)迎角變化情況下的多段翼型主翼和襟翼周圍的流場平均流占據(jù)總能量的絕大部分，因而通過POD處理后，第一階模態(tài)捕捉到的結構就是平均流，其他模態(tài)分析的是脈動結構.隨著模態(tài)階數(shù)的升高，POD所捕捉的流場結構更為復雜和不規(guī)則，分布在模態(tài)間的能量收斂開始變慢，即高階模態(tài)所捕捉到的復雜結構對于總體的影響開始趨向于平均，各自占有的能量比重十分接近.

(3)利用POD方法能進一步獲得隱藏在平均流場背景中低階模態(tài)的復雜相干結構.同一位置在二、三階模態(tài)中分別出現(xiàn)了順時針脫體渦結構和逆向流動結構，但由于脫體渦出現(xiàn)在更低階的模態(tài)中，因此，其對整個流場的影響更大.而同樣對流場影響較大的還有主翼和襟翼的尾部出現(xiàn)的同向尾渦.

［1］ WANGA X，MAY C，YANW J.The proper orthogonal decomposition method for Navier-Stokes equations［J］.ACAD J XJTU，2008，20(3)：141-148.

［2］ VOLKWEINS.Proper orthogonal decomposition(POD) for nonlinear dynamical systems ［R］. Austria：University of Graz Disc Summerschool，2005：1-26.

［3］ 李孝偉.基于嵌套網(wǎng)格的全機帶襟、副翼繞流的數(shù)值模擬［D］.西安：西北工業(yè)大學，1999：8-10.

［4］ JAMESONA，SCHMIDTW，TURKELE.Numerical solution of the Euler equations by finite volume methods with Runge-Kuuta time stepping schemes［C］∥ AIAA Paper.1981：81-1259.

［5］ 楊琴，符松.超音速平面混合層的POD分析［J］.中國科學：G輯，2008，38(3)：319-336.

［6］ EVERSONR，SIROVICHL.Karhunen-Loeve procedure for gappy data［J］.J Opt Soc Am，1995，12(8)：1657-1664.

［7］ BENEKJ A，BUNINGP G，STEGERJ L.A 3-D chimera grid embedding technique［C］∥ AIAA Paper.1985：322-331.

［8］ 杜超.多段翼型非定常粘性繞流數(shù)值研究［D］.上海：上海大學，2007：22-29.

［9］ CHATTERJEE A. An introduction to the proper orthogonal decomposition［J］.Current Science，2000，78(7)：808-817.