一種基于動模態分解的翼型流動轉捩預測新方法

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

一種基于動模態分解的翼型流動轉捩預測新方法

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

西北工業大學 航空學院 翼型／葉柵空氣動力學國家級重點實驗室，西安 710072

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

轉捩預測；動模態分解；eN方法；翼型；計算流體力學

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

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

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

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

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

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

1 翼型流場數據的動模態分解方法

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

首先，提取流場信息。將從實驗或數值計算中獲得的流場數據表示成如下快照序列的形式：

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

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

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

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

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

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

同時，可以得到

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

接下來的工作就是從利用A的伴隨矩陣提取基于數據序列Vn1的動態特征。可以參照文獻［22］通過對Vn1的QR分解來解出S的具體形式，其中Q為正交矩陣；R為上三角矩陣。Schmid［22］指出這種解法雖然理論上正確，但實際應用中卻是一個無法提取出多于第一階或前兩階主導模態的病態方法。尤其當快照數據被噪聲或其他因素污染時，這種現象更加嚴重。為了解決上述問題，選擇魯棒性更好的滿秩矩陣珘S代替S。

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

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

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

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

最后，對DMD分析得到的A的特征值進行對數映射，得到空間模式下的特征值為

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

2 基于DMD的翼型轉捩預測方法

2．1 DMD／eN方法

本文運用DMD取代線性穩定性分析，并與eN轉捩預測方法結合，提出了一種適用于翼型定常與非定常流動的轉捩預測新方法。

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

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

本文提出如下將動模態分解方法和eN轉捩方法結合的過程：

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

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

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

5）逐段向后推進求解，找出擾動放大因子超過轉捩閥值時對應的段標號l。通過Nl與轉捩閥值的比較，線性插值得到閥值對應的轉捩點位置xtran。

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

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

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

理論分析和證明DMD／eN方法的理論依據，以及與基于線性穩定性分析的eN方法（LST／eN）的聯系。

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

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

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

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

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

結合式（13）和式（14）可得

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

參照LST方法，可以構造擾動的表達式為

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

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

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

這就證明了DMD分析得到的特征值為某頻率的擾動在該段上的平均放大率。回顧證明過程，式（13）為DMD理論合理假設，而式（16）為類似于LST理論的合理假設，其余均為嚴格的數學推導。

在完成了DMD進行穩定性分析的理論推導后，下面將對DMD／eN轉捩預測方法的原理做一理論解釋。

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

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

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

3 算 例

3．1 算例1：NLF0416翼型定常流動轉捩預測

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

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

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

表1 DMD計算的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計算的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方法的轉捩點計算結果與實驗值進行對比。上表面和下表面的結果及誤差分析分別如表3和表4所示。誤差計算方法為其中，實驗測得的是翼型表面層流區和湍流區位置，轉捩點應位于兩者之間。由表3和表4可以看出，由100快照得到的上下表面轉捩點位置均在實驗值對應的層流區處，與線性穩定性分析的結果相近，是合理的結果。將得到的轉捩點位置回帶至RANS求解器，迭代至流場收斂。得到的阻力系數CD與實驗值及全湍流的比較如表5所示，可以看出自由轉捩計算所得的阻力系數相比全湍流更加接近實驗值。同時將翼型表面的壓力分布，與全湍流計算的結果及實驗值進行比較，如圖11所示。圖中：方塊為實驗測得的壓力系數Cp分布；虛線為全湍流計算的結果；實線為DMD分析得到的轉捩點帶回RANS求解器得到的壓力系數分布。可以看出，考慮轉捩和全湍流計算的壓力系數分布與實驗值均吻合良好，但是在上表面轉捩位置附近和上表面后緣處，自由轉捩的結果與實驗值更加吻合。圖11（b）給出了壓力系數分布的局部放大，可以看出在壓力系數分布拐折處自由轉捩結果與實驗值更接近。

表3 NLF0416翼型上表面轉捩點位置計算值與實驗比較Table 3 Comparison of predicted transition location with experimental data for upper surface of NLF0416 airfoil

表4 NLF0416翼型下表面轉捩點位置計算值與實驗比較Table 4 Comparison of predicted transition location with experimental data for lower surface of NLF 0416 airfoil

表5 NLF0416翼型計算阻力系數與實驗值的比較Table 5 Comparison of calculated drag coefficient with experimental value for NLF0416airfoil

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

3．2 算例2：S809翼型定常流動轉捩預測

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

表6 DMD計算的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計算的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翼型上、下表面轉捩點位置計算值與實驗值比較Table 8 Comparison of predicted transition location with experimental data for upper and lower surface of S809airfoil

圖14和圖15分別為自由轉捩和全湍流計算得到的S809翼型的整體及分離泡附近的局部放大流線圖。可以看出，自由轉捩計算很好地模擬出了分離泡的存在，而全湍流則沒能模擬出流動中的分離泡。圖16是全湍流及考慮轉捩計算的壓力系數分布與實驗值的比較，圖中方塊表示實驗測得的壓力系數分布，點畫線是全湍流計算的結果，實線是流場快照為100時得到的轉捩位置回帶至RANS求解器的結果。可以看出，自由轉捩的計算的壓力系數分布更加接近實驗值，尤其在分離泡處的模擬結果明顯優于全湍流。其中，翼型上表面半弦長處有明顯的分離泡，下表面在半弦長靠前的位置也存在小的分離泡，這與文獻［27］中實驗結果的描述也完全一致。總之，通過對NLF0416及S809翼型的轉捩點位置的計算及與實驗值的比較，驗證了DMD／eN方法用于翼型定常流動轉捩預測的可行性和正確性，同時也表明DMD／eN方法能夠處理含層流分離泡的轉捩預測問題。

3．3 SD7003翼型非定常流動轉捩預測算例

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

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

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

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

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

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

4 結 論

1）DMD方法能直接從流場數據中分析翼型邊界層內的擾動放大，無需求解穩定性方程。多個驗證算例表明該方法預測的轉捩位置均與實驗值吻合良好，且與線性穩定性分析的結果相當，說明本文提出的方法是正確的。

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

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

4）對于非定常效應不顯著的翼型繞流，采用本文發展的轉捩預測方法，對每個時刻進行空間DMD分析和轉捩判定，仍然可以得到較合理的結果。

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

致 謝

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

［1］ 符松，王亮．湍流轉捩模式研究進展［J］．力學進展，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］ 周恒．關于轉捩和湍流的研究［C］／2003空氣動力學前沿研究論文集．北京：中國空氣動力學會，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轉捩模型在繞翼型流動問題中 的 應 用 ［J］．空 氣 動 力學學報，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］ 符松，王亮．基于雷諾平均方法的高超音速邊界層轉捩模擬［J］．中國科學：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］．空氣動力學學報，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］ 劉方良．低雷諾數翼型流動轉捩判斷與優化設計方法［D］．西安：西北工業大學，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方程的機翼邊界層轉捩判斷及應用研究［D］．西安：西北工業大學，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；網絡出版時間：2016－08－08 08：41

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

國家自然科學基金 （11302177）；民機專項 （MJ－2015－F－016）

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

韓忠華，王紹楠，韓莉，等．一種基于動模態分解的翼型流動轉捩預測新方法［J］．航空學報，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．

（責任編輯：鮑亞平）

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