基于本征正交分解方法的多段翼型流動分析

范晨麟， 李孝偉

(上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上海200072)

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

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

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

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

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

1 控制方程

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

式中，

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

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

2 POD方法

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

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

式中，Ω為計算流場域，

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

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

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

關(guān)于POD理論的更多介紹請參見文獻［1］.

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

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

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

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

4 算例與分析

4.1 方法驗證

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

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

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

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

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

本研究對迎角變化的多段翼型流場進行了數(shù)值模擬，選用的多段翼型為帶30%襟翼的GAW-2翼型結(jié)構(gòu).采用的嵌套網(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對應(yīng)的65階本征模態(tài)各本征值占全部動能的比例.在snapshots計算結(jié)果中，占支配地位的本征模態(tài)通常代表流場中具有較高動能的大尺度湍流結(jié)構(gòu).從圖中可以看出，在翼型周圍的流場中，一階模態(tài)本征值幾乎占總動能的87%-u和77%-v，而二階模態(tài)占總動能的大約8%-u和19%-v，三階模態(tài)占總動能的大約2%-u和2%-v.由此可見，高階模態(tài)對總動能的貢獻相對較低，因此，可以認(rèn)為低階模態(tài)是整個流動結(jié)構(gòu)的載體，而較高階模態(tài)則包含了小尺度的復(fù)雜流動結(jié)構(gòu)的信息.

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

圖4 原流場Fig.4 Original flow field

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

圖6 三階重構(gòu)流場Fig.6 Rank 3 approximation of flow field

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

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

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

5 結(jié)論

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

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

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

(3)利用POD方法能進一步獲得隱藏在平均流場背景中低階模態(tài)的復(fù)雜相干結(jié)構(gòu).同一位置在二、三階模態(tài)中分別出現(xiàn)了順時針脫體渦結(jié)構(gòu)和逆向流動結(jié)構(gòu)，但由于脫體渦出現(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è)大學(xué)，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］.中國科學(xué)：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］.上海：上海大學(xué)，2007：22-29.

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