基于本征正交分解方法的多段翼型流動(dòng)分析

范晨麟， 李孝偉

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

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

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

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

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

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

1 控制方程

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

式中，

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

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

2 POD方法

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

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

式中，Ω為計(jì)算流場(chǎng)域，

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

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

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

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

由于本工作中處理的數(shù)據(jù)集合點(diǎn)數(shù)較為龐大，如果直接求解特征問(wèn)題(4)，其求解過(guò)程會(huì)非常困難.為了解決這個(gè)問(wèn)題，本工作采用Sirovich等［6］提出的快照技術(shù)(snapshots)［1］，因?yàn)樵摲椒軌蛟诤艽蟪潭壬蠝p少計(jì)算對(duì)于內(nèi)存的消耗，即采用原函數(shù)空間元素vn的線性組合來(lái)表示特征模態(tài)，即

同時(shí)，對(duì)R(x，x')進(jìn)行離散，可得

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

其中當(dāng)N足夠大時(shí)，式(11)可化簡(jiǎn)為

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

最后，可以得到

即

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

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

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

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

4 算例與分析

4.1 方法驗(yàn)證

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

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

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

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

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

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

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

圖4 原流場(chǎng)Fig.4 Original flow field

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

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

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

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

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

5 結(jié)論

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

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

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

(3)利用POD方法能進(jìn)一步獲得隱藏在平均流場(chǎng)背景中低階模態(tài)的復(fù)雜相干結(jié)構(gòu).同一位置在二、三階模態(tài)中分別出現(xiàn)了順時(shí)針脫體渦結(jié)構(gòu)和逆向流動(dòng)結(jié)構(gòu)，但由于脫體渦出現(xiàn)在更低階的模態(tài)中，因此，其對(duì)整個(gè)流場(chǎng)的影響更大.而同樣對(duì)流場(chǎng)影響較大的還有主翼和襟翼的尾部出現(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)格的全機(jī)帶襟、副翼繞流的數(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］.中國(guó)科學(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.