基于開源計算流體力學軟件SU2的鈍體翼尖渦研究

崔 青,錢曉輝,劉劍明

(1.中國航空工業(yè)空氣動力研究院，沈陽 110034;2.南京航空航天大學 航空學院，南京 210016;3.江蘇師范大學 數(shù)學與統(tǒng)計學院，江蘇 徐州 221116)

眾所周知，升力面后面的流動如機翼，會導致持續(xù)且強烈對稱的有組織尾渦渦旋結構或翼尖渦旋,而近翼尖渦旋，對于流體力學和空氣動力學研究具有重要的應用。因此，人們對翼尖渦旋進行了廣泛的實驗和數(shù)值研究。Devenport與Rife等[1]通過風洞實驗使用熱線給出了矩形NACA0012翼下游遠場方向尾渦。Giuni 與Green[2]使用煙霧可視化技術研究了NACA0012翼方形翼尖與圓形翼尖的近場渦形成，揭示了低雷諾數(shù)3000下的近場主渦、次渦的相互作用。除風洞實驗外，計算流體力學CFD工具在工程應用中也越來越流行，計算流體力學已經(jīng)成為除風洞試驗外，最經(jīng)濟有效的方法。比如張明華與劉薇等[3-4]開展不同湍流模型翼尖渦模擬研究。最近，Liang與Xue等[5]通過分離渦模擬(DES)方法研究了雷諾數(shù) 1.8×105攻角 α=8°、10° 時，矩形NACA0015翼方形翼尖渦，計算得到的翼尖流向渦量與實驗結果進行了比較，但因渦量包含剪切或拉伸[6,7]，計算結果并不能客觀反映渦旋運動。Lee與Han等[8]針對NACA 0012矩形翼，通過數(shù)值方法研究了雷諾平均(RANS)下，當雷諾數(shù)Re=1×105攻角 13° 與 30° 時的方形翼尖渦生成，同樣使用渦量也混合了剪切與拉伸，增加了分辨次渦與主渦的相互作用的難度。因此，為了進一步研究翼尖渦的形成與相互作用，有必要考察新的渦表示方式。

渦是世界上普遍存在的自然現(xiàn)象。渦旋運動在流體中起關鍵作用，但如何解決這個問題仍然是有爭議的，而缺乏統(tǒng)一的渦旋定義，往往導致渦旋的可視化和理解湍流中的渦旋結構的困惑[6-7，9]。實際上，湍流的顯著特征之一是流體充滿各種大小和強度的無數(shù)旋渦。旋渦的形成與發(fā)展和湍流的生成與發(fā)展密切相關。但是，缺乏普遍接受的渦旋定義嚴重阻礙人們深入理解湍流的機理[6-7，9]。最近，Liu C等[9-10]指出渦量不能代表當?shù)亓黧w的剛性旋轉，而應該將渦量進一步分解為旋轉部分和非旋轉部分?；谶@一概念Liu C等[4,6-7]給出了一種剛性旋轉渦向量Rortex/Liutex，詳細討論也可參見作者及他人近期的論文[11-13]。

本文通過開源計算平臺SU2(Stanford University Unstructured)，通過Rortex/Liutex渦定義，研究了矩形NACA0012翼方形翼尖純剛體渦的生成發(fā)展，揭示了次渦與主渦生成與相互作用。

1 SU2計算平臺與軟件結構

目前CFD的開源軟件中，openfoam是最為熱門的軟件，吸引了大量的科研人員的使用與繼續(xù)開發(fā)，但openfoam比較繁雜，具有一定的學習難度。最近，對于航空航天相關問題，斯坦福大學的開源軟件SU2因其程序架構清晰，獲得了廣泛的關注，引起了很多航空設計課題組的興趣[14-16]。作為越來越流行的SU2軟件，它提供了在非結構網(wǎng)格下求解基于偏微分方程(PDE)問題的軟件一般框架，目前其主要的兩個功能是CFD和氣動優(yōu)化[8]。SU2軟件框架由C++程序及Python腳本集成而成，通過C++類面向?qū)ο筇卣鳎褂美^承、多態(tài)等實現(xiàn)軟件不同功能及多物理場模塊，支持任意變形、運動網(wǎng)格與旋轉坐標系，結合外部自適應網(wǎng)格模塊實現(xiàn)高效的網(wǎng)格自適應。

SU2的高層次C++類結構中，作為最具代表性的SU2軟件的功能模塊是SU2_CFD。在SU2_CFD中，最頂層定義了一個驅(qū)動程序類 CDriver,如圖1所示[14]，它主要控制多物理場模擬的解過程。在此類中，實例化幾何、物理以及求解特定問題的數(shù)值方法。CDriver類中包括下面類的實例：CConfig、COutput、CIntegration、CIteration，其中最核心的類就是CIntegration和CIteration。在SU2_CFD中使用CIntegration的指針通過多態(tài)調(diào)用子類CMultiGridIntegration或CSingleGridIntegration來集成特定的控制方程，并開展多重網(wǎng)格加速計算，通過CGeometry、CSolver和CNumerics子類的實例進行時間和空間上的積分，而CIteration 類及其子類完成不同物理模塊的單步迭代。在SU2軟件中，計算工具的核心功能就是分別嵌入管理幾何的CGeometry類，求解器功能的CSolver類和數(shù)值方法CNumerics類中。

CGeometry類讀入與處理網(wǎng)格，輸入的網(wǎng)格結構是一種特有的以SU2為后綴的網(wǎng)格，這種網(wǎng)格已經(jīng)可以通過商用軟件比如Pointwise直接生成。CGeometry類的子類CPhysicalGeometry 建立了SU2軟件中格點有限體積方法的對偶網(wǎng)格。CSolver類定義了求解的過程，它的每一個子類用來求解一個特定的控制方程，比如CEulerSolver 類給出求解可壓縮無粘性Euler方程的求解器，而CTurbSolver類用于求解湍流模型。CNumerics類使用輸入文件中設置的方案參數(shù)，通過不同指定方案離散化控制方程組，再經(jīng)由有幾個子類為對流通量、粘性通量以及給定的偏微分方程中可能存在的任何源項提供廣泛的離散化技術。

2 數(shù)值方法與Rortex/Liutex渦定義

SU2是關于偏微分方程模型的集成解決平臺，可以求解各種通過偏微分方程描述的物理工程問題，作為例子，我們僅介紹CFD相關部分。

考慮微分形式的Navier-Stokes方程

(1)

其積分形式是

(2)

對公式(2)，使用對偶網(wǎng)格格點有限體積方法，采用二階近似可以得到半離散形式

(3)

其中對流通量積分可以使用HLLC、Roe等數(shù)值通量，通過使用MUSCL方法得到二階精度，梯度計算可以采用Green-Gauss方法或者最小二乘法，而粘性項則采用中心格式，時間推進采用Runge-Kutta顯格式或者LU-SGS隱格式。為了封閉RANS方程(1)或者(2)的求解，需要求解湍流模型，在SU2中主要提供SA湍流模型與SST湍流模型，作為推薦，一般對于湍流模型的求解，只使用一階精度迎風格式。如果考慮大分離流動，可以使用脫體渦模擬方法DES/DDES[17]。

因為渦量混合了剪切與拉伸，不能區(qū)分層流邊界層與旋轉。為了顯示計算產(chǎn)生的渦結構，需要使用特殊的渦識別方法[9,10]。目前工程中比較流行的渦識別方法是Q標準[18]，表示為速度梯度張量▽V的第2個伽利略不變量

(4)

其中B、A分別是速度梯度張量的反對稱張量與對稱張量。使用Q標準，在Q>0的區(qū)域，認為有旋轉，并通過其等值面表示渦的形態(tài)，但其不能區(qū)分剪切，從而會造成剪切污染，且渦的形態(tài)和閾值相關，不同的閾值可能得到完全不同的渦形態(tài)，并產(chǎn)生錯誤的結論[9-10]。渦應該既有大小也有方向，也就是旋轉軸，但目前流行的渦表示方法Q標準、λ2與Δ等都是標量，且并沒有明確的物理意義[7,9,10]。

最近LiuC等[7,9,11]提出了一種渦向量定義，Rortex/Liutex向量，不但給出了旋轉方向，還有剛體旋轉大小，并能夠區(qū)分剪切層與旋轉。Rortex/Liutex向量的大小可以用于渦顯示，其渦形態(tài)不會被剪切污染，而且具有強的渦形態(tài)保持能力[13]。根據(jù)Rortex/Liutex渦理論，當速度梯度張量有一個實特征值與兩個共軛復特征值時，代表瞬時流線呈現(xiàn)圓形或螺旋型，即存在旋渦結構，此與Δ方法和λci方法在判斷是否存在渦結構是一致的，但Δ方法和λci方法僅使用到速度梯度張量的不變量和特征值，并沒有給出明確的旋轉方向的概念，而且用于表示渦的量沒有明確的物理意義，且缺乏特征向量方向所提供的信息。當速度梯度張量有3個實特征值時，代表流體微團在3個特征方向只有拉伸或壓縮，沒有旋轉運動，也就是沒有渦，此時Rortex/Liutex向量為0；當速度梯度張量有一個實特征值λr和兩個共軛的復特征值λcr±iλci時，實特征值對應的特征向量方向只有拉伸或壓縮，而渦旋轉運動只能發(fā)生在垂直于實特征向量r的平面內(nèi)，此時的實特征向量方向即為當?shù)亓黧w微團的旋轉軸，而大小定義為垂直于實特征向量r的平面內(nèi)兩倍的最小角速度。用簡單的顯公式，可以將Rortex向量表示為[13]

(5)

其中ω表示渦量，r表示實特征值對應單位特征向量方向(取和渦量內(nèi)積大于0的方向)，λci表示復特征值的虛部。

3 數(shù)值實驗與分析

3.1 代碼驗證

為驗證SU2軟件的有效性，首先考察二維NACA0012攻角與升力系數(shù)的關系圖[19]。此標準模型問題的流動狀態(tài)中，雷諾數(shù)Re=6×106，馬赫數(shù)為0.2，計算非結構網(wǎng)格如圖2a所示，計算采用二階Roe格式離散對流項，并使用SA湍流模型封閉RANS方程的計算。圖2b給出了升力系數(shù)CL與攻角α的計算結果以及實驗比較。目前的計算除了在高攻角時因分離流動湍流模型的緣故與實驗有偏差外，攻角一直到 15° 升力系數(shù)結果都和實驗能夠很好吻合。此外，升力系數(shù)CL與攻角α也較好顯示了升力系數(shù)CL與攻角α的2π斜率規(guī)律。

圖2 數(shù)值驗證

3.2 方形翼尖渦識別

對于三維復雜流體，本文考慮三維NACA0012翼型鈍體翼尖渦流，雷諾數(shù)是Re=2×106，馬赫數(shù)為0.2，攻角12°，半翼展長為2.5倍弦長。計算采用混合六面體與三棱柱網(wǎng)格，共大約200萬個非結構網(wǎng)格，最小網(wǎng)格壁面距離約為2×10-6。圖3a為翼端壁面網(wǎng)格，圖3b為翼端2.5%弦長截面位置(如圖3a所示)的壓強系數(shù)與實驗結果的比較，結果除了因上翼面主旋轉渦的存在，壓強系數(shù)與實驗大約在x/c=0.35和x/c=0.75 附近有偏差外，其他位置和實驗相符。

圖3 NAC0012翼型鈍體翼尖網(wǎng)格與壓強系數(shù)

圖4a為使用渦識別Q標準給出的翼端渦等值圖，清晰顯示渦在翼梢卷起旋轉拖出去的長尾渦，且流線圍繞渦等值面。圖4b為渦向量Rortex大小R的等值面圖，此等值面清晰反映了鈍體翼尖渦。比較圖4兩圖可以發(fā)現(xiàn)，使用Q標準，翼面上很多地方Q>0，而 Rortex量R只在翼尖與前緣顯示，這與Q標準會被剪切污染，而Rortex是純剛體旋轉，沒有剪切污染的理論相符[7，9]。

圖4 鈍體翼尖渦Q標準與渦向量Rortex量等值面圖

圖5 方形翼尖不同截面位置渦量云圖

圖6 方形翼尖不同截面位置Rortex量云圖

圖7 Rortex剛體旋轉量R在渦量量中的占比

4 結論

本文采用開源軟件SU2，使用新的渦定義Rortex/Liutex向量，針對矩形翼NACA0012的方形翼尖，在雷諾數(shù)Re=2×106，馬赫數(shù)為0.2，攻角12°，半翼展長為2.5倍弦長條件下，研究了方形翼尖渦的不同位置渦的形態(tài)與相互作用，并對比了和渦量結果的關系，得到如下結論：

(1)通過典型二維與三維物體的數(shù)值實驗，分析和驗證SU2軟件，說明SU2軟件能夠較好解決相關問題，結果和實驗能夠匹配;

(2)對方形翼尖渦的研究發(fā)現(xiàn)Q標準會出現(xiàn)剪切污染，而Rortex向量的量沒有剪切污染，能找到清晰的旋轉渦結構;

(3)方形翼尖渦的側面奇異邊產(chǎn)生渦量與Rortex/Liutex向量，但渦量在側邊主要貢獻剪切，其中一部分才是旋轉渦，而且Rortex/Liutex次渦與主渦開始相互作用的位置要晚于渦量。