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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

同時，可以得到

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

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

下面簡述滿秩矩陣珘S的求解過程：

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

式中：U為Vn－11的左奇異矩陣；Σ為半正定對角矩陣，對角線上的元素是奇異值，Σ＝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的特征向量矩陣。也就是說，矩陣A和矩陣珘S的特征值相同，均為Λ＝diag（μ1，μ2，…，μn－1）。這樣，通過滿秩矩陣珘S就得到了A的特征值和特征模態(tài)。

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

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

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

2．1 DMD／eN方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算 例

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

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

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

圖10給出了NLF0416算例中LST／eN方法與DMD／eN方法的擾動增長曲線對比。黑色粗虛線是DMD／eN方法的N 值增長曲線；其他為LST／eN方法不同特征頻率f擾動波的N 值增長曲線。從圖中看出，DMD／eN方法的N 值增長曲線很接近LST／eN方法的包絡(luò)線。雖然個別站位處有些差異，但大體的增長趨勢是一致的，轉(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)行對比。上表面和下表面的結(jié)果及誤差分析分別如表3和表4所示。誤差計(jì)算方法為其中，實(shí)驗(yàn)測得的是翼型表面層流區(qū)和湍流區(qū)位置，轉(zhuǎn)捩點(diǎn)應(yīng)位于兩者之間。由表3和表4可以看出，由100快照得到的上下表面轉(zhuǎn)捩點(diǎn)位置均在實(shí)驗(yàn)值對應(yīng)的層流區(qū)處，與線性穩(wěn)定性分析的結(jié)果相近，是合理的結(jié)果。將得到的轉(zhuǎn)捩點(diǎn)位置回帶至RANS求解器，迭代至流場收斂。得到的阻力系數(shù)CD與實(shí)驗(yàn)值及全湍流的比較如表5所示，可以看出自由轉(zhuǎn)捩計(jì)算所得的阻力系數(shù)相比全湍流更加接近實(shí)驗(yàn)值。同時將翼型表面的壓力分布，與全湍流計(jì)算的結(jié)果及實(shí)驗(yàn)值進(jìn)行比較，如圖11所示。圖中：方塊為實(shí)驗(yàn)測得的壓力系數(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)附近有一個突增，對應(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翼型定常流動轉(zhuǎn)捩預(yù)測

為了進(jìn)一步驗(yàn)證本文方法的正確性，選擇了轉(zhuǎn)捩點(diǎn)位置實(shí)驗(yàn)值豐富的風(fēng)力機(jī)翼型S809［27］進(jìn)行DMD分析與轉(zhuǎn)捩預(yù)測。計(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ù)測時，閥值N取為6較為合適。在進(jìn)行DMD分析時為了排除數(shù)值噪聲的影響，σ取10－8。分別取50和100快照對翼型每段進(jìn)行DMD分析，得到的上表面和下表面擾動放大率及N值的發(fā)展過程分別如表6和表7所示?？梢钥闯?00個流場快照對應(yīng)的上表面和下表面超過轉(zhuǎn)捩閥值的點(diǎn)，分別在60%弦長和50%弦長處。對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翼型的整體及分離泡附近的局部放大流線圖?？梢钥闯?，自由轉(zhuǎn)捩計(jì)算很好地模擬出了分離泡的存在，而全湍流則沒能模擬出流動中的分離泡。圖16是全湍流及考慮轉(zhuǎn)捩計(jì)算的壓力系數(shù)分布與實(shí)驗(yàn)值的比較，圖中方塊表示實(shí)驗(yàn)測得的壓力系數(shù)分布，點(diǎn)畫線是全湍流計(jì)算的結(jié)果，實(shí)線是流場快照為100時得到的轉(zhuǎn)捩位置回帶至RANS求解器的結(jié)果?？梢钥闯觯杂赊D(zhuǎn)捩的計(jì)算的壓力系數(shù)分布更加接近實(shí)驗(yàn)值，尤其在分離泡處的模擬結(jié)果明顯優(yōu)于全湍流。其中，翼型上表面半弦長處有明顯的分離泡，下表面在半弦長靠前的位置也存在小的分離泡，這與文獻(xiàn)［27］中實(shí)驗(yàn)結(jié)果的描述也完全一致?？傊?，通過對NLF0416及S809翼型的轉(zhuǎn)捩點(diǎn)位置的計(jì)算及與實(shí)驗(yàn)值的比較，驗(yàn)證了DMD／eN方法用于翼型定常流動轉(zhuǎn)捩預(yù)測的可行性和正確性，同時也表明DMD／eN方法能夠處理含層流分離泡的轉(zhuǎn)捩預(yù)測問題。

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

致 謝

本文所提出的轉(zhuǎn)捩預(yù)測新方法，是在與西班牙馬德里理工大學(xué)（UPM）的Esteban Ferrer博士的合作過程中得到啟發(fā)，主要由我們課題組教師和研究生歷時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空氣動力學(xué)前沿研究論文集．北京：中國空氣動力學(xué)會，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)捩模型在繞翼型流動問題中 的 應(yīng) 用 ［J］．空 氣 動 力學(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］．中國科學(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)捩自動判斷的阻力精確計(jì)算方法初探［J］．空氣動力學(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ù)翼型流動轉(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ò)出版時間：2016－08－08 08：41

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

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

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

韓忠華，王紹楠，韓莉，等．一種基于動模態(tài)分解的翼型流動轉(zhuǎn)捩預(yù)測新方法［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