韓忠華，王紹楠，韓莉，劉方良，許建華，宋文萍

一種基于動(dòng)模態(tài)分解的翼型流動(dòng)轉(zhuǎn)捩預(yù)測(cè)新方法

韓忠華＊，王紹楠，韓莉，劉方良，許建華，宋文萍

西北工業(yè)大學(xué) 航空學(xué)院 翼型／葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，西安 710072

考慮自由轉(zhuǎn)捩的定常／非定常流動(dòng)Navier－Stokes方程數(shù)值求解，對(duì)于翼型流動(dòng)細(xì)節(jié)的精確模擬和氣動(dòng)力的精確預(yù)測(cè)均具有十分重要的意義。采用動(dòng)模態(tài)分解（DMD）方法進(jìn)行流動(dòng)穩(wěn)定性分析，再結(jié)合eN方法，提出了一套適用于翼型繞流的轉(zhuǎn)捩預(yù)測(cè)新方法，稱為DMD／eN方法。相比于傳統(tǒng)的線性穩(wěn)定性分析方法，DMD方法不需要求解附面層方程和線性穩(wěn)定性方程，也沒有引入平行流假設(shè)，具有更好的理論適用性和算法魯棒性。開展了NLF0416、S809和SD7003等翼型的轉(zhuǎn)捩預(yù)測(cè)數(shù)值驗(yàn)證研究，通過(guò)與實(shí)驗(yàn)結(jié)果以及與傳統(tǒng)的基于線性穩(wěn)定性分析的eN方法的比較，驗(yàn)證了本文所發(fā)展的轉(zhuǎn)捩預(yù)測(cè)新方法在預(yù)測(cè)翼型的定常流動(dòng)和非定常流動(dòng)轉(zhuǎn)捩方面的正確性，也表明了該方法具有解決含層流分離泡的翼型繞流轉(zhuǎn)捩預(yù)測(cè)的能力。

轉(zhuǎn)捩預(yù)測(cè)；動(dòng)模態(tài)分解；eN方法；翼型；計(jì)算流體力學(xué)

流動(dòng)由層流轉(zhuǎn)變成湍流的過(guò)程稱為轉(zhuǎn)捩，轉(zhuǎn)捩現(xiàn)象普遍存在于流體流動(dòng)中［1］。高雷諾數(shù)下的機(jī)翼、大前緣半徑的風(fēng)力機(jī)葉片、高空螺旋槳和微型飛行器等的表面流動(dòng)都存在轉(zhuǎn)捩現(xiàn)象。由于影響因素眾多、物理機(jī)理復(fù)雜，轉(zhuǎn)捩問題一直是流體力學(xué)中尚未完全解決的前沿問題之一［2］。層流和湍流兩種流態(tài)下物體表面的摩擦阻力、熱傳導(dǎo)速率和流動(dòng)分離位置等特性大為不同，故轉(zhuǎn)捩位置直接影響飛行器氣動(dòng)力計(jì)算的精確性，其中對(duì)阻力預(yù)測(cè)的影響尤為顯著。因此，發(fā)展定常／非定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)方法，對(duì)于飛行器流動(dòng)細(xì)節(jié)的精確模擬和氣動(dòng)力的準(zhǔn)確預(yù)測(cè)均具有十分重要的意義。

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，大渦模擬（LES）和直接數(shù)值模擬（DNS）被應(yīng)用于轉(zhuǎn)捩研究［3－5］。但由于其計(jì)算量過(guò)于龐大，只能用于低雷諾數(shù)平板和翼型的研究，仍不適合工程應(yīng)用。目前工程上常用的轉(zhuǎn)捩預(yù)測(cè)方法主要有基于雷諾平均的湍流／轉(zhuǎn)捩模式和基于線性穩(wěn)定性理論的eN方法。前者在近年來(lái)有了較快發(fā)展，具有代表性的是基于當(dāng)?shù)刈兞康腉amma－Theta轉(zhuǎn)捩模型［6－7］，且在其基礎(chǔ)上考慮橫流的轉(zhuǎn)捩預(yù)測(cè)方法［8］也有了一定進(jìn)展。符松和王亮［9］發(fā)展的適用于超聲速邊界層轉(zhuǎn)捩預(yù)測(cè)的k－ω－γ三方程模型，也成為超聲速流動(dòng)轉(zhuǎn)捩研究的熱點(diǎn)。盡管如此，較早提出的基于線性穩(wěn)定性理論的eN方法仍是目前工程中最常使用的轉(zhuǎn)捩預(yù)測(cè)方法。其廣泛用于翼型／機(jī)翼的邊界層轉(zhuǎn)捩預(yù)測(cè)［10－11］。

相比于定常流動(dòng)，非定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)更加復(fù)雜。德國(guó)國(guó)家航空航天研究院 （DLR）的Krumbein等［12］將拓展的eN轉(zhuǎn)捩預(yù)測(cè)方法與非定常流動(dòng)的雷諾平均Navier－Stokes（RANS）方程相結(jié)合，并以做俯仰振蕩翼型進(jìn)行了算例驗(yàn)證。德國(guó)布倫瑞克工業(yè)大學(xué)的 Radespiel等［13－14］應(yīng)用該方法，實(shí)現(xiàn)了低雷諾數(shù)下翼型非定常層流分離流動(dòng)的轉(zhuǎn)捩位置自動(dòng)預(yù)測(cè)。在國(guó)內(nèi)，將eN轉(zhuǎn)捩預(yù)測(cè)方法應(yīng)用于翼型、機(jī)翼繞流的非定常數(shù)值模擬研究并不多見。劉方良［15］通過(guò)直接求解Navier－Stokes方程為預(yù)測(cè)邊界層轉(zhuǎn)捩提供高精度解，將eN方法與非定常RANS求解器耦合，研究了低雷諾數(shù)螺旋槳翼型的流動(dòng)特性。然而，傳統(tǒng)的基于線性穩(wěn)定性分析的eN轉(zhuǎn)捩預(yù)測(cè)方法，不能預(yù)測(cè)轉(zhuǎn)捩過(guò)程中非線性效應(yīng)的影響；同時(shí)由于其建立在平行流假設(shè)的基礎(chǔ)上，忽略了流動(dòng)橫向速度和流向參數(shù)的弱增長(zhǎng)。

為了克服線性穩(wěn)定性分析的缺點(diǎn)，PSE（Parabolized Stability Equations）方法［16－17］和全局穩(wěn)定性分析方法［18］被提出。同時(shí)模態(tài)分解方法也逐漸成為穩(wěn)定性分析和流場(chǎng)相干結(jié)構(gòu)分析的有力工具。其中，本征正交分解［19－21］（Proper Orthogonal Decomposition，POD）方法和動(dòng)模態(tài)分解［22］（Dynamic Mode Decomposition，DMD）方法是目前應(yīng)用比較廣泛的兩種模態(tài)分解方法。

2010年，Schmid［22］在 Koopman分析［23］的基礎(chǔ)上提出了動(dòng)模態(tài)分解方法，隨后便成為眾多研究者研究的熱點(diǎn)［24－25］。該方法的基本思想是直接從實(shí)驗(yàn)或數(shù)值模擬得到的流場(chǎng)中提取出流動(dòng)的動(dòng)態(tài)信息；它既可用來(lái)表征流動(dòng)的全局穩(wěn)定性和動(dòng)力學(xué)特征，又可用來(lái)解釋流場(chǎng)的物理流動(dòng)機(jī)理。不同于POD方法的每個(gè)模態(tài)對(duì)應(yīng)著多個(gè)頻率，DMD方法可以提取出隨著單一頻率變化的空間模態(tài)，因而可以表征不同頻率的擾動(dòng)在流動(dòng)中沿著空間方向的增長(zhǎng)。同時(shí)，DMD適用于包含線性增長(zhǎng)或衰減的復(fù)雜頻率，這使得該方法可以分析和解決流動(dòng)不穩(wěn)定性問題。目前DMD已被應(yīng)用于邊界層流動(dòng)的穩(wěn)定性分析研究。但是，國(guó)內(nèi)外還沒有被應(yīng)用于邊界層轉(zhuǎn)捩預(yù)測(cè)研究的報(bào)道。

本文在基于線性穩(wěn)定性分析的eN方法基礎(chǔ)上，提出了一種用DMD替代線性穩(wěn)定性分析的轉(zhuǎn)捩預(yù)測(cè)新方法，稱為DMD／eN方法。運(yùn)用該方法，開展了繞二維翼型的定常及非定常效應(yīng)不顯著的準(zhǔn)定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)方法研究，驗(yàn)證了該方法的可行性和正確性，以及在處理含層流分離泡轉(zhuǎn)捩問題方面的能力。

1 翼型流場(chǎng)數(shù)據(jù)的動(dòng)模態(tài)分解方法

這里簡(jiǎn)要介紹DMD方法的數(shù)值求解過(guò)程，詳細(xì)理論推導(dǎo)及算法介紹參見文獻(xiàn)［22］。DMD方法可分為時(shí)間DMD和空間DMD，兩者主要是快照劃分方向的不同。由于本文關(guān)心的是流動(dòng)在空間方向的發(fā)展規(guī)律，且通過(guò)分析，空間DMD更適用于本文研究的轉(zhuǎn)捩預(yù)測(cè)方法，因而下面主要介紹空間DMD方法。

首先，提取流場(chǎng)信息。將從實(shí)驗(yàn)或數(shù)值計(jì)算中獲得的流場(chǎng)數(shù)據(jù)表示成如下快照序列的形式：

式中：n為快照數(shù)；vi為第i個(gè)抽樣數(shù)據(jù)的流場(chǎng)信息，稱為快照（Snapshots）；Vn1為流場(chǎng)快照形成的矩陣，下標(biāo)和上標(biāo)分別為起始和終止的流場(chǎng)快照編號(hào)。

空間DMD方法假設(shè)相鄰兩個(gè)快照之間的空間間隔（Δx）是相同的，故需要將原有計(jì)算網(wǎng)格的流場(chǎng)信息插值到新的等間距分布的DMD網(wǎng)格線上，從而獲得所需的流場(chǎng)快照信息。本文用于穩(wěn)定性分析的流場(chǎng)狀態(tài)量是速度信息，即（u，v），其中u和v分別為直角坐標(biāo)系下x方向和y方向的速度。具體表示為

其次，基于流場(chǎng)快照，形成包含系統(tǒng)流動(dòng)特征結(jié)構(gòu)及時(shí)空演化信息的矩陣A。如果采樣的樣本Vn1大到可以完全描述原流場(chǎng)的時(shí)空演化特征，可假設(shè)兩個(gè)相鄰的流場(chǎng)快照vi和vi＋1可以由線性變化A聯(lián)系起來(lái)，即

同時(shí)，假設(shè)這種變換關(guān)系在Δx足夠小時(shí)對(duì)于整個(gè)抽樣空間［0，（n－1）Δx］均適用。由于分布間隔相同，該流動(dòng)序列可由Krylov序數(shù)表示，即

DMD理論假設(shè)當(dāng)快照的個(gè)數(shù)超過(guò)某個(gè)值后，式（4）的向量為線性相關(guān)，即超過(guò)某一臨界個(gè)數(shù)后，整體向量空間將保持不變，此時(shí)即使再增加更多的快照也不會(huì)對(duì)其產(chǎn)生影響。達(dá)到臨界值后，第n個(gè)向量vn可由前n－1個(gè)向量線性表示，

式中：aT＝［a1a2… an－1］；r為殘值向量。

同時(shí)，可以得到

式中：en－1∈Rn－1為n－1階單位向量；Vn2和Vn－11均從CFD數(shù)值模擬結(jié)果中提取。由于矩陣A是一個(gè)反映流場(chǎng)內(nèi)在動(dòng)力學(xué)機(jī)制，且具體形式未知的高維矩陣，所以DMD方法采用A的伴隨矩陣S來(lái)作為A的低維近似。S的表達(dá)式為

接下來(lái)的工作就是從利用A的伴隨矩陣提取基于數(shù)據(jù)序列Vn1的動(dòng)態(tài)特征?？梢詤⒄瘴墨I(xiàn)［22］通過(guò)對(duì)Vn1的QR分解來(lái)解出S的具體形式，其中Q為正交矩陣；R為上三角矩陣。Schmid［22］指出這種解法雖然理論上正確，但實(shí)際應(yīng)用中卻是一個(gè)無(wú)法提取出多于第一階或前兩階主導(dǎo)模態(tài)的病態(tài)方法。尤其當(dāng)快照數(shù)據(jù)被噪聲或其他因素污染時(shí)，這種現(xiàn)象更加嚴(yán)重。為了解決上述問題，選擇魯棒性更好的滿秩矩陣珘S代替S。

下面簡(jiǎn)述滿秩矩陣珘S的求解過(guò)程：

對(duì)由流場(chǎng)快照組成的矩陣Vn－11進(jìn)行奇異值分解，得到

式中：U為Vn－11的左奇異矩陣；Σ為半正定對(duì)角矩陣，對(duì)角線上的元素是奇異值，Σ＝diag（σ1，σ2，…，σn－1），其中σi按從大到小的順序進(jìn)行排列；WT為的右奇異矩陣。代入式（6），兩端左乘UT和右乘以WΣ－1得到

這樣就得到了滿秩矩陣珘S的具體形式。隨后采用QR分解求珘S的特征值和特征向量。珘S＝Y(jié)ΛY－1代入式（9）得到AUY＝UYΛ，令＝UY。Y為珘S的特征向量矩陣；為A的特征向量矩陣。也就是說(shuō)，矩陣A和矩陣珘S的特征值相同，均為Λ＝diag（μ1，μ2，…，μn－1）。這樣，通過(guò)滿秩矩陣珘S就得到了A的特征值和特征模態(tài)。

最后，對(duì)DMD分析得到的A的特征值進(jìn)行對(duì)數(shù)映射，得到空間模式下的特征值為

式中：μ＝μre＋μimi為DMD方法得到的流動(dòng)特征值；λre為指數(shù)放大率或衰減率；λim為空間波數(shù)，i為虛數(shù)單位（i2＝－1）。

2 基于DMD的翼型轉(zhuǎn)捩預(yù)測(cè)方法

2．1 DMD／eN方法

本文運(yùn)用DMD取代線性穩(wěn)定性分析，并與eN轉(zhuǎn)捩預(yù)測(cè)方法結(jié)合，提出了一種適用于翼型定常與非定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)新方法。

基于線性穩(wěn)定性分析的eN方法，由于平行流假設(shè)的引入，使得流動(dòng)的橫向速度和流向參數(shù)的弱增長(zhǎng)被忽略。該方法不僅對(duì)層流解的精度要求較高，穩(wěn)定性方程求解的魯棒性欠佳，并且不能直接用于非定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)。而本文采用的DMD方法，直接從實(shí)驗(yàn)或者CFD數(shù)值模擬結(jié)果中提取數(shù)據(jù)，不對(duì)數(shù)據(jù)做任何假設(shè)，克服了線性穩(wěn)定性分析的以上缺點(diǎn)。文獻(xiàn)研究表明，DMD的分析結(jié)果等價(jià)于全局穩(wěn)定性分析的結(jié)果［22］，可以很好捕捉系統(tǒng)的動(dòng)力學(xué)特征，得到不穩(wěn)定的主導(dǎo)模態(tài)及其對(duì)應(yīng)的特征值。DMD方法的優(yōu)勢(shì)在于無(wú)須解龐大的全局穩(wěn)定性矩陣即可獲得全局動(dòng)態(tài)特征。

2010年Schmid提出了基于瞬時(shí)流場(chǎng)的DMD方法，并且比較了平面Poiseuille流線性小擾動(dòng)的Orr－Sommerfeld方程的精確特征值與DMD分析所得的特征值，驗(yàn)證了DMD方法的收斂性和有效性［22］。

本文提出如下將動(dòng)模態(tài)分解方法和eN轉(zhuǎn)捩方法結(jié)合的過(guò)程：

1）將翼型上下表面沿翼型弦線方向等分L段。理論上分成的段數(shù)越多得到的轉(zhuǎn)捩點(diǎn)位置越精確，但是需要的流場(chǎng)快照數(shù)也越多；而為保證流場(chǎng)快照的精度則要求更大的網(wǎng)格量，從而帶來(lái)更大的計(jì)算量。若是分段過(guò)少，則轉(zhuǎn)捩位置精度很難保證。本文作為方法的初步探索，將翼型上、下表面各分為10段。

2）對(duì)翼型表面每一段的流場(chǎng)快照數(shù)據(jù)進(jìn)行DMD計(jì)算，得到每一段表征穩(wěn)定性的特征值。值得注意的是，在DMD實(shí)際應(yīng)用中，常常要去除流場(chǎng)中的噪聲以及不重要的流動(dòng)結(jié)構(gòu)對(duì)DMD計(jì)算的干擾。而對(duì)流場(chǎng)快照矩陣進(jìn)行奇異值分解后，得到的奇異值大小代表所對(duì)應(yīng)的流場(chǎng)結(jié)構(gòu)對(duì)流場(chǎng)的貢獻(xiàn)，奇異值σi越大表示對(duì)應(yīng)的模態(tài)在流場(chǎng)中所占比重越大。其中，奇異值過(guò)小的模態(tài)對(duì)流場(chǎng)的貢獻(xiàn)很小，有時(shí)可能是由數(shù)值誤差引起的。因此奇異值的大小可以作為預(yù)處理的閥值。根據(jù)流動(dòng)特點(diǎn)確定合適的閥值σ，選擇σi／σ1≥σ的奇異值對(duì)應(yīng)的模態(tài)進(jìn)行分析，目的是降低數(shù)值誤差對(duì)結(jié)果的影響。本文選擇σ＝10－8進(jìn)行分析計(jì)算。

3）確定流場(chǎng)演化矩陣A的有效特征值。本文參照文獻(xiàn)［15］中通過(guò)對(duì)每一段內(nèi)取不同快照進(jìn)行DMD分析（例如取80個(gè)和100個(gè)流場(chǎng)快照），選擇不穩(wěn)定特征值（λim，λre）分布重合度最好的特征值作為本段的有效特征值，對(duì)應(yīng)的模態(tài)即是本段的特征模態(tài)。研究中發(fā)現(xiàn)數(shù)值噪聲對(duì)DMD分

5）逐段向后推進(jìn)求解，找出擾動(dòng)放大因子超過(guò)轉(zhuǎn)捩閥值時(shí)對(duì)應(yīng)的段標(biāo)號(hào)l。通過(guò)Nl與轉(zhuǎn)捩閥值的比較，線性插值得到閥值對(duì)應(yīng)的轉(zhuǎn)捩點(diǎn)位置xtran。

本文發(fā)展的繞翼型定常流動(dòng)的轉(zhuǎn)捩預(yù)測(cè)具體過(guò)程如圖2所示。首先采用固定轉(zhuǎn)捩方法對(duì)繞翼型的流動(dòng)進(jìn)行數(shù)值模擬，在保證較大層流范圍的前提下使流場(chǎng)充分收斂；其次，提取邊界層內(nèi)的速度信息，插值形成流場(chǎng)快照；而后，采用DMD／eN方法進(jìn)行轉(zhuǎn)捩預(yù)測(cè)；最后，將得到轉(zhuǎn)捩位置回帶至RANS求解器，直至得到收斂的解。析的結(jié)果影響較大，在這里選取兩種快照數(shù)重疊性最好的特征值的目的是為了避免數(shù)值噪聲干擾。由于80快照和100快照的特征值中重疊最好的兩個(gè)也存在些微小差異（并不是完全重合），故以流場(chǎng)快照數(shù)多的100快照的結(jié)果為準(zhǔn)。

4）利用式（11）對(duì)DMD分析得到的空間放大率進(jìn)行積分，得到從翼型前緣開始到第l段末尾的擾動(dòng)累計(jì)放大因子Nl（下標(biāo)l為段標(biāo)號(hào)）：對(duì)于非定常效應(yīng)不是十分顯著的翼型繞流，本文采用如圖3所示的轉(zhuǎn)捩預(yù)測(cè)過(guò)程。首先，進(jìn)行固定轉(zhuǎn)捩的定常RANS求解，得到一個(gè)初步的流場(chǎng)；以定常計(jì)算的初始流場(chǎng)為非定常計(jì)算0時(shí)刻的初值，進(jìn)行非定常RANS的迭代，直至得到周期性變化的力系數(shù)，則認(rèn)為流場(chǎng)收斂，保存不同物理時(shí)刻的流場(chǎng)數(shù)值。其次，提取非定常計(jì)算中一個(gè)周期內(nèi)不同物理時(shí)刻的速度流場(chǎng)數(shù)據(jù)，分別采用DMD／eN方法進(jìn)行轉(zhuǎn)捩預(yù)測(cè)。最后，將不同時(shí)刻的轉(zhuǎn)捩位置帶回非定常RANS進(jìn)行計(jì)算，直至得到收斂的解。

2．2 DMD／eN方法原理分析

理論分析和證明DMD／eN方法的理論依據(jù)，以及與基于線性穩(wěn)定性分析的eN方法（LST／eN）的聯(lián)系。

LST／eN方法是在定常層流邊界層中引入一個(gè)小擾動(dòng)，通過(guò)解線性穩(wěn)定性方程來(lái)確定中性曲線上不同頻率的擾動(dòng)在流向不同站位處特征值（空間波數(shù)和放大率）；然后，對(duì)每個(gè)頻率下的放大率沿流向進(jìn)行積分，得到若干N值增長(zhǎng)曲線，并通過(guò)與閥值的比較得到轉(zhuǎn)捩點(diǎn)。

而本文提出的DMD／eN方法采用類似的思路，其理論的合理性取決于如下一個(gè)結(jié)論：DMD分析得到的每段流場(chǎng)演化矩陣A的特征值（對(duì)數(shù)映射后的實(shí)部），是該段“特征擾動(dòng)”的平均放大率。所謂特征擾動(dòng)，是指由DMD分析的特征值所對(duì)應(yīng)的擾動(dòng)。下面對(duì)這一結(jié)論進(jìn)行理論證明：

在翼型表面劃分的某一段流場(chǎng)內(nèi)，DMD分析的流場(chǎng)快照矩陣為

式中：珋vi為層流邊界層的速度。根據(jù)DMD理論，假設(shè)流場(chǎng)快照數(shù)足夠多，包含了所有空間演化信息，則兩個(gè)相鄰快照間的關(guān)系可以寫為

假設(shè)在第i站位處平均速度項(xiàng)中存在擾動(dòng)δi，沿流向進(jìn)行演化得到i＋1站位的δi＋1。假設(shè)擾動(dòng)足夠小，使得加上擾動(dòng)后的流場(chǎng)仍然滿足關(guān)系式（13），則有

結(jié)合式（13）和式（14）可得

式（15）表明，第i＋1站位處的擾動(dòng)可由矩陣A與第i站位處的擾動(dòng)聯(lián)系起來(lái)。A是由流場(chǎng)的內(nèi)在動(dòng)力學(xué)機(jī)制（主控方程）決定。而A的特征值μ＝e（λre＋iλim）Δx和 特 征 向 量 反 映 了 擾 動(dòng) 在 空 間 的演化。

參照LST方法，可以構(gòu)造擾動(dòng)的表達(dá)式為

式中：k為由矩陣A得到的圓頻率為ωk的擾動(dòng)的特征向量，下標(biāo)k為自定義的不同時(shí)間頻率標(biāo)號(hào)；λre和λim與式（10）中含義一致，分別為 DMD分析得到的特征值μ對(duì)數(shù)映射后的放大率和空間波數(shù)。

由式（15）和特征值的定義可知，在第i站位處的擾動(dòng)，到i＋1站位處演化為

照數(shù)從1到n），平均的擾動(dòng)放大率為

這就證明了DMD分析得到的特征值為某頻率的擾動(dòng)在該段上的平均放大率。回顧證明過(guò)程，式（13）為DMD理論合理假設(shè)，而式（16）為類似于LST理論的合理假設(shè)，其余均為嚴(yán)格的數(shù)學(xué)推導(dǎo)。

在完成了DMD進(jìn)行穩(wěn)定性分析的理論推導(dǎo)后，下面將對(duì)DMD／eN轉(zhuǎn)捩預(yù)測(cè)方法的原理做一理論解釋。

DMD／eN方法的原理可以概述為：從定常層流解中提取邊界層速度信息，分析得到流向每段上“特征擾動(dòng)”模態(tài)和對(duì)應(yīng)的特征值；由于所提取的特征值實(shí)部對(duì)應(yīng)了擾動(dòng)在該段上的放大率，因此對(duì)其放大率進(jìn)行積分可以得到N值增長(zhǎng)曲線，并通過(guò)與閥值的對(duì)比來(lái)預(yù)測(cè)轉(zhuǎn)捩點(diǎn)位置。與LST／eN方法不同的是，DMD是數(shù)據(jù)驅(qū)動(dòng)的穩(wěn)定性分析方法，它不需要通過(guò)求解線性穩(wěn)定性方程得到一系列圓頻率下的擾動(dòng)增長(zhǎng)曲線。雖然在構(gòu)造擾動(dòng)時(shí)也考慮了時(shí)間頻率項(xiàng)，但計(jì)算過(guò)程中并不需要計(jì)算該頻率具體是多少，而只關(guān)心“特征擾動(dòng)”的放大率信息。

在運(yùn)用DMD進(jìn)行穩(wěn)定性分析的基礎(chǔ)上，本文提出的轉(zhuǎn)捩預(yù)測(cè)方法實(shí)際上類似于eN包絡(luò)線方法。所計(jì)算得到的每段上重疊最好的特征值，一般來(lái)講對(duì)應(yīng)的LST／eN方法中包絡(luò)線上平均放大率最大的擾動(dòng)。需要說(shuō)明的是，每段上的有效特征值對(duì)應(yīng)的并不是同一個(gè)時(shí)間頻率，這與eN包絡(luò)線法中包絡(luò)線上對(duì)應(yīng)的有時(shí)不是同一時(shí)間頻率是類似的。

總之，所提出的DMD／eN轉(zhuǎn)捩預(yù)測(cè)方法，不需要算出不同時(shí)間頻率下擾動(dòng)的增長(zhǎng)情況，而是采用了類似于eN包絡(luò)線方法的原理（算例部分將會(huì)與基于線性穩(wěn)定性分析的eN包絡(luò)線方法作比較）。也就是先計(jì)算出每段“特征擾動(dòng)”的放大情況，進(jìn)行疊加得到一條類似包含若干擾動(dòng)的包絡(luò)線，然后將閥值與包絡(luò)線上的N值增長(zhǎng)做對(duì)比，預(yù)測(cè)出轉(zhuǎn)捩點(diǎn)位置。

3 算 例

3．1 算例1：NLF0416翼型定常流動(dòng)轉(zhuǎn)捩預(yù)測(cè)

為 驗(yàn) 證 DMD／eN方 法 的 有 效 性，選 擇NLF0416翼型［26］進(jìn)行數(shù)值模擬。計(jì)算狀態(tài)為：馬赫數(shù)Ma＝0．1，雷諾數(shù)Re＝4×106，迎角α＝0°。采用RANS求解器，湍流模型為Spalart－Allmaras模型。圖4給出了NLF0416的網(wǎng)格圖，網(wǎng)格單元數(shù)為640×240，其中翼型上下表面共有512個(gè)網(wǎng)格單元，第1層網(wǎng)格距物面高度為2×10－6c。其中，c為翼型弦長(zhǎng)。采用課題組自主研發(fā)的二維流動(dòng)求解程序PMNS2D對(duì)翼型進(jìn)行固定轉(zhuǎn)捩數(shù)值模擬，上下表面的初始固定轉(zhuǎn)捩點(diǎn)位置設(shè)在約60%弦長(zhǎng)處（必須保證足夠的層流區(qū)用于DMD分析）。當(dāng)流場(chǎng)充分收斂后，提取流場(chǎng)信息，通過(guò)插值形成流場(chǎng)快照。圖5給出了翼型上表面0．1c～0．2c區(qū)間，流場(chǎng)快照取80，σ＝10－8的模態(tài)特征值（μre和μim分別為矩陣A的特征值μ的實(shí)部和虛部）。位于單位圓之外的特征值是不穩(wěn)定模態(tài)對(duì)應(yīng)的特征值；位于單位圓上的特征值對(duì)應(yīng)的模態(tài)是平均流動(dòng)；位于單位圓之內(nèi)的特征值是穩(wěn)定特征值，其對(duì)應(yīng)的模態(tài)是穩(wěn)定的模態(tài)。圖6給出了分別采用50、80和100快照數(shù)分析得到的特征值分布圖。上三角、方框與倒三角分別表示流場(chǎng)快照數(shù)n＝100，80，50的特征值分布?？梢钥闯?種流場(chǎng)快照數(shù)分析得到的特征值呈收斂的趨勢(shì)。由于橫坐標(biāo)和縱坐標(biāo)的尺度相差較大，故圖中看似重疊最好的點(diǎn)并不是距離最小的點(diǎn)，實(shí)線框中的3個(gè)為距離最小的不穩(wěn)定特征值，代表著該段的特征模態(tài)。無(wú)論采用50快照或80快照標(biāo)定都得到是同樣的有效特征值。這里采用了80快照來(lái)選擇重疊性最好的100快照的特征值作為有效特征值。

圖7 給出了有效特征值的選取示意圖，實(shí)線方框選中的點(diǎn)是不穩(wěn)定特征值中重疊性最好的點(diǎn)。根據(jù)2．1節(jié)的特征值選取標(biāo)準(zhǔn)，這里選擇100快照對(duì)應(yīng)的特征值為有效特征值，80快照相當(dāng)于標(biāo)定作用。按照2．1節(jié)敘述的方法，分別得到翼型上下表面每一段80快照和100快照對(duì)應(yīng)的擾動(dòng)放大率及N值發(fā)展過(guò)程（上表面的結(jié)果見表1，下表面的結(jié)果見表2）。表中x／c為翼型表面的坐標(biāo)；λre為對(duì)應(yīng)站位的擾動(dòng)放大率；N為對(duì)擾動(dòng)放大率沿弦向積分得到的N 值。這里只展示了初始固定轉(zhuǎn)捩層流區(qū)內(nèi)（60%弦長(zhǎng)范圍）的N值增長(zhǎng)，此后為湍流區(qū)，按照原理不必進(jìn)行DMD分析。轉(zhuǎn)捩閥值取6，以快照數(shù)多的（100快照）N值發(fā)展曲線作為轉(zhuǎn)捩位置判斷依據(jù)，分別將超過(guò)轉(zhuǎn)捩閥值的位置及對(duì)應(yīng)的轉(zhuǎn)捩閥值用黑體字表示。從表中右半部分100快照的對(duì)應(yīng)列可看出上表面N值在40%弦長(zhǎng)處超過(guò)閥值，下表面在50%弦長(zhǎng)處超過(guò)閥值。同時(shí)，給出了線性穩(wěn)定性分析的結(jié)果，轉(zhuǎn)捩閥值N也取為6。翼型上表面和下表面的擾動(dòng)增長(zhǎng)曲線分別如圖8和圖9所示，圖中：xtr＿up和xtr＿low分別為上表面和下表面轉(zhuǎn)捩閥值的位置。

圖10給出了NLF0416算例中LST／eN方法與DMD／eN方法的擾動(dòng)增長(zhǎng)曲線對(duì)比。黑色粗虛線是DMD／eN方法的N 值增長(zhǎng)曲線；其他為L(zhǎng)ST／eN方法不同特征頻率f擾動(dòng)波的N 值增長(zhǎng)曲線。從圖中看出，DMD／eN方法的N 值增長(zhǎng)曲線很接近LST／eN方法的包絡(luò)線。雖然個(gè)別站位處有些差異，但大體的增長(zhǎng)趨勢(shì)是一致的，轉(zhuǎn)捩點(diǎn)位置比LST／eN方法計(jì)算的略微靠前，但誤差很小，在合理范圍內(nèi)。這也在一定程度上驗(yàn)證了DMD／eN方法的正確性。

表1 DMD計(jì)算的NLF0416翼型上表面N值（Ma＝0．1，Re＝4×106，α＝0°，σ＝10－8，N＝6）Table 1 Nfactor calculated by DMD for upper surface of NLF0416airfoil（Ma＝0．1，Re＝4×106，α＝0°，σ＝10－8，N＝6）

表2 DMD計(jì)算的NLF0416翼型下表面N值（Ma＝0．1，Re＝4×106，α＝0°，σ＝10－8，N＝6）Table 2 Nfactor calculated by DMD for lower surface of NLF0416airfoil（Ma＝0．1，Re＝4×106，α＝0°，σ＝10－8，N＝6）

將 DMD／eN方法與 LST／eN方法的轉(zhuǎn)捩點(diǎn)計(jì)算結(jié)果與實(shí)驗(yàn)值進(jìn)行對(duì)比。上表面和下表面的結(jié)果及誤差分析分別如表3和表4所示。誤差計(jì)算方法為其中，實(shí)驗(yàn)測(cè)得的是翼型表面層流區(qū)和湍流區(qū)位置，轉(zhuǎn)捩點(diǎn)應(yīng)位于兩者之間。由表3和表4可以看出，由100快照得到的上下表面轉(zhuǎn)捩點(diǎn)位置均在實(shí)驗(yàn)值對(duì)應(yīng)的層流區(qū)處，與線性穩(wěn)定性分析的結(jié)果相近，是合理的結(jié)果。將得到的轉(zhuǎn)捩點(diǎn)位置回帶至RANS求解器，迭代至流場(chǎng)收斂。得到的阻力系數(shù)CD與實(shí)驗(yàn)值及全湍流的比較如表5所示，可以看出自由轉(zhuǎn)捩計(jì)算所得的阻力系數(shù)相比全湍流更加接近實(shí)驗(yàn)值。同時(shí)將翼型表面的壓力分布，與全湍流計(jì)算的結(jié)果及實(shí)驗(yàn)值進(jìn)行比較，如圖11所示。圖中：方塊為實(shí)驗(yàn)測(cè)得的壓力系數(shù)Cp分布；虛線為全湍流計(jì)算的結(jié)果；實(shí)線為DMD分析得到的轉(zhuǎn)捩點(diǎn)帶回RANS求解器得到的壓力系數(shù)分布?？梢钥闯觯紤]轉(zhuǎn)捩和全湍流計(jì)算的壓力系數(shù)分布與實(shí)驗(yàn)值均吻合良好，但是在上表面轉(zhuǎn)捩位置附近和上表面后緣處，自由轉(zhuǎn)捩的結(jié)果與實(shí)驗(yàn)值更加吻合。圖11（b）給出了壓力系數(shù)分布的局部放大，可以看出在壓力系數(shù)分布拐折處自由轉(zhuǎn)捩結(jié)果與實(shí)驗(yàn)值更接近。

表3 NLF0416翼型上表面轉(zhuǎn)捩點(diǎn)位置計(jì)算值與實(shí)驗(yàn)比較Table 3 Comparison of predicted transition location with experimental data for upper surface of NLF0416 airfoil

表4 NLF0416翼型下表面轉(zhuǎn)捩點(diǎn)位置計(jì)算值與實(shí)驗(yàn)比較Table 4 Comparison of predicted transition location with experimental data for lower surface of NLF 0416 airfoil

表5 NLF0416翼型計(jì)算阻力系數(shù)與實(shí)驗(yàn)值的比較Table 5 Comparison of calculated drag coefficient with experimental value for NLF0416airfoil

因湍流黏性系數(shù)遠(yuǎn)大于層流黏性系數(shù)，故層流與湍流兩種流態(tài)下物面的摩擦阻力差異很大。在圖12中給出了NLF0416翼型全湍流和自由轉(zhuǎn)捩狀態(tài)下的物面摩擦力系數(shù)Cf分布。其中虛線是全湍流摩阻系數(shù)，點(diǎn)畫線是LST／eN方法得到的結(jié)果，實(shí)線是DMD／eN方法得出的結(jié)果。由圖可看出，LST／eN和 DMD／eN兩種方法得到的摩阻系數(shù)分布吻合良好。翼型上下表面的摩阻系數(shù)分別在各自轉(zhuǎn)捩點(diǎn)附近有一個(gè)突增，對(duì)應(yīng)位置在0．35和0．50左右。轉(zhuǎn)捩點(diǎn)位置之前是層流區(qū)，全湍流的模擬結(jié)果與自由轉(zhuǎn)捩差異很大，Cf遠(yuǎn)大于自由轉(zhuǎn)捩的結(jié)果，而轉(zhuǎn)捩點(diǎn)之后進(jìn)入湍流，Cf的值開始逐步接近。

3．2 算例2：S809翼型定常流動(dòng)轉(zhuǎn)捩預(yù)測(cè)

為了進(jìn)一步驗(yàn)證本文方法的正確性，選擇了轉(zhuǎn)捩點(diǎn)位置實(shí)驗(yàn)值豐富的風(fēng)力機(jī)翼型S809［27］進(jìn)行DMD分析與轉(zhuǎn)捩預(yù)測(cè)。計(jì)算狀態(tài)：Ma＝0．1、Re＝2×106、α＝0°。轉(zhuǎn)捩后的湍流區(qū)域采用Spalart－Allmaras湍流模型。S809翼型計(jì)算網(wǎng)格如圖13所示，計(jì)算網(wǎng)格量為640×256，第一層網(wǎng)格高度為2×10－6c。根據(jù)以往的經(jīng)驗(yàn)，在利用eN方法進(jìn)行風(fēng)力機(jī)翼型繞流的轉(zhuǎn)捩預(yù)測(cè)時(shí)，閥值N取為6較為合適。在進(jìn)行DMD分析時(shí)為了排除數(shù)值噪聲的影響，σ取10－8。分別取50和100快照對(duì)翼型每段進(jìn)行DMD分析，得到的上表面和下表面擾動(dòng)放大率及N值的發(fā)展過(guò)程分別如表6和表7所示?？梢钥闯?00個(gè)流場(chǎng)快照對(duì)應(yīng)的上表面和下表面超過(guò)轉(zhuǎn)捩閥值的點(diǎn)，分別在60%弦長(zhǎng)和50%弦長(zhǎng)處。對(duì)N值插值后的上、下表面的轉(zhuǎn)捩點(diǎn)位置與實(shí)驗(yàn)值［27］的比較如表8所示。可以看出，本文得到的上、下表面轉(zhuǎn)捩點(diǎn)位置均與實(shí)驗(yàn)值吻合良好。

表6 DMD計(jì)算的S809翼型上表面N值（Ma＝0．1，Re＝2×106，α＝0°，σ＝10－8，N＝6）Table 6 Nfactor calculated by DMD for upper surface of S809airfoil（Ma＝0．1，Re＝2×106，α＝0°，σ＝10－8，N＝6）

表7 DMD計(jì)算的S809翼型下表面N值（Ma＝0．1，Re＝2×106，α＝0°，σ＝10－8，N＝6）Table 7 Nfactor calculated by DMD for lower surface of S809airfoil（Ma＝0．1，Re＝2×106，α＝0°，σ＝10－8，N＝6）

表8 S809翼型上、下表面轉(zhuǎn)捩點(diǎn)位置計(jì)算值與實(shí)驗(yàn)值比較Table 8 Comparison of predicted transition location with experimental data for upper and lower surface of S809airfoil

圖14和圖15分別為自由轉(zhuǎn)捩和全湍流計(jì)算得到的S809翼型的整體及分離泡附近的局部放大流線圖?？梢钥闯觯杂赊D(zhuǎn)捩計(jì)算很好地模擬出了分離泡的存在，而全湍流則沒能模擬出流動(dòng)中的分離泡。圖16是全湍流及考慮轉(zhuǎn)捩計(jì)算的壓力系數(shù)分布與實(shí)驗(yàn)值的比較，圖中方塊表示實(shí)驗(yàn)測(cè)得的壓力系數(shù)分布，點(diǎn)畫線是全湍流計(jì)算的結(jié)果，實(shí)線是流場(chǎng)快照為100時(shí)得到的轉(zhuǎn)捩位置回帶至RANS求解器的結(jié)果?？梢钥闯?，自由轉(zhuǎn)捩的計(jì)算的壓力系數(shù)分布更加接近實(shí)驗(yàn)值，尤其在分離泡處的模擬結(jié)果明顯優(yōu)于全湍流。其中，翼型上表面半弦長(zhǎng)處有明顯的分離泡，下表面在半弦長(zhǎng)靠前的位置也存在小的分離泡，這與文獻(xiàn)［27］中實(shí)驗(yàn)結(jié)果的描述也完全一致?？傊?，通過(guò)對(duì)NLF0416及S809翼型的轉(zhuǎn)捩點(diǎn)位置的計(jì)算及與實(shí)驗(yàn)值的比較，驗(yàn)證了DMD／eN方法用于翼型定常流動(dòng)轉(zhuǎn)捩預(yù)測(cè)的可行性和正確性，同時(shí)也表明DMD／eN方法能夠處理含層流分離泡的轉(zhuǎn)捩預(yù)測(cè)問題。

3．3 SD7003翼型非定常流動(dòng)轉(zhuǎn)捩預(yù)測(cè)算例

本文以低雷諾數(shù)下SD7003翼型做俯仰運(yùn)動(dòng)的非定常流動(dòng)為例，進(jìn)行了基于DMD的非定常轉(zhuǎn)捩預(yù)測(cè)。為了保證邊界層信息的精度，選擇網(wǎng)格量為640×256，第1層網(wǎng)格高度為1×10－6c，網(wǎng)格示意圖如圖17所示。計(jì)算狀態(tài)：Ma＝0．05，Re＝6×104，運(yùn)動(dòng)規(guī)律：α＝5．5°±3°sin（ωt），其中減縮頻率為k＝ωc／2U∞＝0．52。

首先進(jìn)行了定常固定轉(zhuǎn)捩計(jì)算，為非定常流動(dòng)提供初始流場(chǎng)。然后進(jìn)入非定常固定轉(zhuǎn)捩計(jì)算，初始固定轉(zhuǎn)捩點(diǎn)位置上、下表面均取為0．6c。其中每個(gè)運(yùn)動(dòng)周期劃分為300個(gè)物理時(shí)間步，每個(gè)物理時(shí)間步內(nèi)子迭代20步，共計(jì)算6個(gè)周期。得到具有周期性變化的力系數(shù)和流場(chǎng)，則認(rèn)為流場(chǎng)收斂。

對(duì)于這種頻率較小，非定常效應(yīng)不是特別明顯的翼型俯仰運(yùn)動(dòng)算例，轉(zhuǎn)捩預(yù)測(cè)過(guò)程中采用準(zhǔn)定常的處理方法，與文獻(xiàn)［15］中采用的處理方法一致。

在非定常流場(chǎng)充分收斂的情況下，分別提取（k／8）T（k＝0，1，2，…，8）時(shí)刻的流場(chǎng)數(shù)據(jù)，進(jìn)行坐標(biāo)轉(zhuǎn)化、流場(chǎng)插值和DMD計(jì)算，提取有效特征值并結(jié)合eN方法得到各自時(shí)刻的轉(zhuǎn)捩點(diǎn)位置，進(jìn)而得到一個(gè)運(yùn)動(dòng)周期內(nèi)的轉(zhuǎn)捩曲線，如圖18所示。這里選取有效特征值的方法與定常算例所采用方法相同，是采用50和100快照重疊性最好的作為有效特征值。

圖18中圓圈是實(shí)驗(yàn)測(cè)得的非定常轉(zhuǎn)捩點(diǎn)位置，虛線是文獻(xiàn)［15］利用線性穩(wěn)定性分析程序結(jié)合eN方法得到的轉(zhuǎn)捩點(diǎn)位置。點(diǎn)畫線是文獻(xiàn)［13］，即Radespiel等通過(guò)在非定常RANS方程求解器中耦合基于線性穩(wěn)定性理論拓展的eN轉(zhuǎn)捩預(yù)測(cè)方法，對(duì)該算例的模擬結(jié)果。其在湍流模型的選擇，穩(wěn)定性分析方法和轉(zhuǎn)捩預(yù)測(cè)方法上與本文有所不同。實(shí)線是本文方法所得的非定常流動(dòng)1個(gè)周期內(nèi)的轉(zhuǎn)捩點(diǎn)位置。可以看出，本文得到的結(jié)果在總體趨勢(shì)上和文獻(xiàn)［13］是一致的，并且在大部分時(shí)刻都與其結(jié)果相近，但在t／T＝2／8，3／8時(shí)刻的轉(zhuǎn)捩點(diǎn)位置與實(shí)驗(yàn)值及文獻(xiàn)值［13］相比稍微靠后。因?yàn)槲墨I(xiàn)［13］采用的是拓展的eN方法，考慮到了流動(dòng)的完全非定常效應(yīng)，而本文與文獻(xiàn)［15］采用的都是標(biāo)準(zhǔn)的eN方法，用準(zhǔn)定常的方式處理，故在一定程度上降低了預(yù)測(cè)精度，這也是在個(gè)別時(shí)刻轉(zhuǎn)捩點(diǎn)位置沒有文獻(xiàn)［13］與實(shí)驗(yàn)值吻合好的原因。

將本文得到的轉(zhuǎn)捩點(diǎn)位置回帶至流場(chǎng)求解程序中，迭代求解至流場(chǎng)收斂。得到典型時(shí)刻的壓力系數(shù)分布及流場(chǎng)的壓力云圖和流線圖，如圖19所示?？梢钥闯鲆硇驮诔跏嘉恢茫é粒?．5°）有一個(gè)小的層流分離泡，隨著翼型的抬頭運(yùn)動(dòng)，分離泡逐漸變大。當(dāng)翼型在達(dá)到最大迎角（α＝8．5°）后做下俯運(yùn)動(dòng)，分離泡逐漸變長(zhǎng)并向翼型后緣移動(dòng)。達(dá)到最低位置（α＝2．5°）后，翼型又進(jìn)行抬頭運(yùn)動(dòng)回到初始位置，分離泡也逐步脫落，最終回到初始狀態(tài)。從圖中能較為清楚地看出1個(gè)周期內(nèi)翼型做俯仰運(yùn)動(dòng)的分離泡的運(yùn)動(dòng)發(fā)展過(guò)程，以及對(duì)應(yīng)的壓力系數(shù)分布的變化情況。通過(guò)與文獻(xiàn)［13］及文獻(xiàn)［15］的對(duì)比可知，本文的方法也能較好地模擬出翼型俯仰運(yùn)動(dòng)中的非定常特性，以及轉(zhuǎn)捩位置和分離泡的發(fā)展變化。其對(duì)于研究和分析低雷諾數(shù)下非定常流動(dòng)特性也具有一定價(jià)值。更加精確的非定常轉(zhuǎn)捩預(yù)測(cè)需要在后續(xù)工作中考慮完全非定常效應(yīng)對(duì)轉(zhuǎn)捩預(yù)測(cè)的影響。

4 結(jié) 論

1）DMD方法能直接從流場(chǎng)數(shù)據(jù)中分析翼型邊界層內(nèi)的擾動(dòng)放大，無(wú)需求解穩(wěn)定性方程。多個(gè)驗(yàn)證算例表明該方法預(yù)測(cè)的轉(zhuǎn)捩位置均與實(shí)驗(yàn)值吻合良好，且與線性穩(wěn)定性分析的結(jié)果相當(dāng)，說(shuō)明本文提出的方法是正確的。

2）在進(jìn)行定常流動(dòng)轉(zhuǎn)捩預(yù)測(cè)時(shí)，通過(guò)選取不同流場(chǎng)快照數(shù)（例如快照數(shù)為50、80和100），對(duì)DMD分析得到的特征值進(jìn)行比較，選取重疊度最好的不穩(wěn)定特征值作為有效特征值，是一種合理的有效特征值選取方法。

3）對(duì)DMD分析得到的翼型表面各段的擾動(dòng)放大N因子進(jìn)行線性插值，可提高轉(zhuǎn)捩位置預(yù)測(cè)的精度。

4）對(duì)于非定常效應(yīng)不顯著的翼型繞流，采用本文發(fā)展的轉(zhuǎn)捩預(yù)測(cè)方法，對(duì)每個(gè)時(shí)刻進(jìn)行空間DMD分析和轉(zhuǎn)捩判定，仍然可以得到較合理的結(jié)果。

本文工作展望如下：針對(duì)三維流動(dòng)機(jī)翼TS波誘導(dǎo)的轉(zhuǎn)捩，本文方法應(yīng)該是能夠直接推廣應(yīng)用的。因?yàn)镈MD在計(jì)算中不依賴邊界條件和控制方程，且能較靈活地應(yīng)用于復(fù)雜流動(dòng)和子域分析。初步設(shè)想為，將機(jī)翼沿展向劃分為若干截面，在每個(gè)截面處應(yīng)用類似二維的轉(zhuǎn)捩預(yù)測(cè)方法得到對(duì)應(yīng)展向位置處的轉(zhuǎn)捩點(diǎn)，通過(guò)插值得到機(jī)翼表面的轉(zhuǎn)捩線，類似于eN方法對(duì)機(jī)翼轉(zhuǎn)捩預(yù)測(cè)的處理［28］。若進(jìn)一步考慮CF波的情況，由于轉(zhuǎn)捩機(jī)理更為復(fù)雜，還有待于進(jìn)一步深入研究。

致 謝

本文所提出的轉(zhuǎn)捩預(yù)測(cè)新方法，是在與西班牙馬德里理工大學(xué)（UPM）的Esteban Ferrer博士的合作過(guò)程中得到啟發(fā)，主要由我們課題組教師和研究生歷時(shí)3年完成的。在此感謝Esteban Ferrer博士的大力支持，感謝他提供了DMD分析程序，并給予我們很多有益的建議。

［1］ 符松，王亮．湍流轉(zhuǎn)捩模式研究進(jìn)展［J］．力學(xué)進(jìn)展，2007，37（3）：409－416．FU S，WANG L．Progress in turbulence／transition modelling［J］．Advances in Mechanics，2007，37（3）：409－416（in Chinese）．

［2］ 周恒．關(guān)于轉(zhuǎn)捩和湍流的研究［C］／2003空氣動(dòng)力學(xué)前沿研究論文集．北京：中國(guó)空氣動(dòng)力學(xué)會(huì)，2003：87－93．ZHOU H．Studies on transition and turbulence［C］／2003 Advanced Research Papers on Aerodynamics．Beijing：Aerodynamic Society of China，2003：87－93（in Chinese）．

［3］ FASEL H F，MEITZ H L，BACHAN C R．DNS and LES for investigating transition and transition control：AIAA－1997－1820［R］．Reston：AIAA，1997．

［4］ SCHLATTER S．Assessment of direct numerical simulation data of turbulent boundary layers［J］．Journal of Fluid Mechanics，2010，659：116－126．

［5］ BRAZELL M J，KIRBY A，STOELLINGER M，et al．Using LES in a Discontinuous Galerkin method with constant and dynamic SGS models：AIAA－2015－0060［R］．Reston：AIAA，2015．

［6］ 陳奕，高正紅．Gamma－Theta轉(zhuǎn)捩模型在繞翼型流動(dòng)問題中 的 應(yīng) 用 ［J］．空 氣 動(dòng) 力學(xué)學(xué)報(bào)，2009，27（4）：411－419．CHEN Y，GAO Z H．Application of Gamma－Theta transition model to flows around airfoils［J］．Acta Aerodynamic Sinica，2009，27（4）：411－419（in Chinese）．

［7］ MENTER F R，LANGTRY R B，LIKKI S R，et al．A correlation－based transition model using local variables—Part I：model formulation［J］．Journal of Turbomachinery，2006，128（3）：413－422．

［8］ GRABE C，KRUMBEIN A．Extension of theγ－Reθtmodel for prediction of crossflow transition：AIAA－2014－1269［R］．Reston：AIAA，2014．

［9］ 符松，王亮．基于雷諾平均方法的高超音速邊界層轉(zhuǎn)捩模擬［J］．中國(guó)科學(xué)：G輯，2009，39（4）：617－626．FU S，WANG L．Modelling flow transition in a hypersonic boundary layer with Reynolds－averaged Navier－Stokes approach［J］．Science in China：Series G，2009，39（4）：617－626（in Chinese）．

［10］ STOCK H W，HAASE W．Navier－Stokes airfoil computations with eNtransition prediction including transitional flow regions［J］．Journal of Aircraft，2000，38（11）：2059－2066．

［11］ 張坤，宋文萍．NS方程計(jì)算中耦合轉(zhuǎn)捩自動(dòng)判斷的阻力精確計(jì)算方法初探［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（4）：400－404．ZHANG K，SONG W P．Accurate drag calculation by coupling automatic prediction of transition point to the Navier－Stokes method［J］．Acta Aerodynamic Sinica，2009，27（4）：400－404（in Chinese）．

［12］ KRUMBEIN A，KRIMMELBEIN N，SEYFERT C．Automatic transition prediction in unsteady airfoil flows using an unstructured CFD code：AIAA－2011－3365［R］．Reston：AIAA，2011．

［13］ RADESPIEL R，WINDTE J，SCHOLZ U．Numerical and experimental flow analysis of moving airfoils with laminar separation bubbles［J］．AIAA Journal，2007，45（6）：1346－1356．

［14］ WINDTE J，RADESPIEL R．Propulsive efficiency of a moving airfoil at transitional low Reynolds numbers［J］．AIAA Journal，2008，46（9）：2165－2177．

［15］ 劉方良．低雷諾數(shù)翼型流動(dòng)轉(zhuǎn)捩判斷與優(yōu)化設(shè)計(jì)方法［D］．西安：西北工業(yè)大學(xué)，2014：33－67．LIU F L．Transition prediction and optimization design method for low－Reynolds－number airfoil［D］．Xi’an：Northwestern Polytechnical University，2014：33－67 （in Chinese）．

［16］ HERBERT T．Parabolized stability equations［J］．Annual Review of Fluid Mechanics，1997，291（1）：245－283．

［17］ BERTOLOTTI F P，HERBERT T．Analysis of the linear stability of compressible boundary layers using the PSE［J］．Theoretical and Computational Fluid Dynamics，1991，3（2）：117－124．

［18］ NOACK B R，MORZYNSKI M，TADMOR G．Reducedorder modelling for flow control［M］．Udine：CISM，2011：77－110．

［19］ CHATTERJEE A．An introduction to the proper orthogonal decomposition［J］．Computational Science，2000，78（7）：808－807．

［20］ SIROVICH L．Turbulence and the dynamics of coherent structures［J］．Quarterly of Applied Mathematics，1987，45（3）：561－590．

［21］ BERKOOZ G，HOLMES P，LUMLEY J L．The proper orthogonal decomposition in the analysis of turbulent flows［J］．Annual Review of Fluid Mechanics，1993，25（1）：539－575．

［22］ SCHMID P J．Dynamic mode decomposition of numerical and experimental data［J］．Journal of Fluid Mechanics，2010，656：5－28．

［23］ ROWLEY C，MEZIC I，BAGHERI S，et al．Spectral analysis of nonlinear flows［J］．Journal of Fluid Mechanics，2009，641：115－127．

［24］ SCHMID P J，LI L，JUNIPER M，et al．Applications of the dynamic mode decomposition［J］．Theoretical and Computational Fluid Dynamics，2011，25（1－4）：249－259．

［25］ SCHMID P J．Application of the dynamic mode decomposition to experimental data［J］．Experiments in Fluids，2011，50（4）：1123－1130．

［26］ SOMERS D M．Design and experimental results for a natural－laminar－flow airfoil for general aviation applications：NASA－TP－1861［R］．Washington，D．C．：NASA，1981．

［27］ SOMERS D M．Design and experimental results for the S809airfoil：NREL／SR－440－6918［R］．Colorado：NREL，1997．

［28］ 張坤．基于NS方程的機(jī)翼邊界層轉(zhuǎn)捩判斷及應(yīng)用研究［D］．西安：西北工業(yè)大學(xué)，2011：97－157．ZHANG K．Automatic transition prediction and application to 3Dswept wings［D］．Xi’an：Northwestern Polytechnical University，2011：97－157（in Chinese）．

A novel method for automatic transition prediction of flows over airfoils based on dynamic mode decomposition

HAN Zhonghua＊，WANG Shaonan，HAN Li，LIU Fangliang，XU Jianhua，SONG Wenping
National Key Laboratory of Science and Technology on Aerodynamic Design and Research，School of Aeronautics，Northwestern Polytechnical University，Xi’an 710072，China

Transition prediction is crucial for the simulation of steady and unsteady flows，since it can improve the accuracy of predicting the aerodynamic forces as well as capturing the flow phenomena．By combining dynamic mode decomposition（DMD）and eNmethod，a novel transition prediction method for flows over airfoils is proposed．Compared with conventional linear stability－analysis－based eNmethod，DMD requires neither the solution of boundary layer and linear stability equations，nor the assumption of parallel flows，and has better applicability in theory and is more algorithmically robust．Transition prediction of steady flows around NLF0416and S809airfoils and unsteady flow around SD7003airfoil are carried out．The predicted transition locations are in reasonably good agreement with the experimental data and the results of eNmethod based on linear stability analysis．It is shown that the proposed DMD／eNmethod is feasible for transition prediction for steady and unsteady flows over airfoils，including the flows with laminar separation bubbles．

transition prediction；dynamic mode decomposition；eNmethod；airfoil；computational fluid dynamics

2016－01－11；Revised：2016－01－31；Accepted：2016－07－21；Published online：2016－08－08 08：41

URL：www．cnki．net／kcms／detail／11．1929．V．20160808．0841．002．html

s：National Natural Science Foundation of China（11302177）；Civil Aircraft Project（MJ－2015－F－016）

V211．3

A

1000－6893（2017）01－120034－17

http：／hkxb．buaa．edu．cn hkxb＠buaa．edu．cn

10．7527／S1000－6893．2016．0225

2016－01－11；退修日期：2016－01－31；錄用日期：2016－07－21；網(wǎng)絡(luò)出版時(shí)間：2016－08－08 08：41

www．cnki．net／kcms／detail／11．1929．V．20160808．0841．002．html

國(guó)家自然科學(xué)基金 （11302177）；民機(jī)專項(xiàng) （MJ－2015－F－016）

＊通訊作者 ．E－mail：hanzh＠nwpu．edu．cn

韓忠華，王紹楠，韓莉，等．一種基于動(dòng)模態(tài)分解的翼型流動(dòng)轉(zhuǎn)捩預(yù)測(cè)新方法［J］．航空學(xué)報(bào)，2017，38（1）：120034．HAN Z H，WANG S N，HAN L，et al．A novel method for automatic transition prediction of flows over airfoils based on dynamic mode decomposition［J］．Acta Aeronautica et Astronautica Sinica，2017，38（1）：120034．

（責(zé)任編輯：鮑亞平）

＊Corresponding author．E－mail：hanzh＠nwpu．edu．cn