DDES方法在復雜旋翼流場計算中的應用

董軍，葉靚

中國航空工業空氣動力研究院， 沈陽 110034

旋翼流場是一個典型的復雜非定常流場。一方面，槳葉的運動規律復雜，在有來流速度時，需要調整槳葉的周期變距操縱，實現旋翼氣動力和力矩平衡；同時在周期性氣動力和慣性力作用下，槳葉還要附加揮舞和擺振運動，槳葉前行側氣流壓縮、后行時的動態失速等都是旋翼特有的氣動現象。另一方面旋翼的工作狀態較多，除去常規的懸停和平飛巡航狀態，還可能進行近地、小速度下降和機動飛行等，在這些工作狀態下，流場的復雜度進一步增加，預測的難度進一步加大。為準確計算旋翼的氣動力，旋翼流場中的流動細節捕捉(激波、槳尖渦和槳葉表面附近空間的氣流分離)則成為關鍵因素。

在先前的計算研究中，基于雷諾平均 Navier-Stokes(RANS)方程的流場預測占據著主導地位[1]。隨著計算流體力學的發展，一些計算量更大、計算精度更高的計算方法越來越得到重視[2-3]。

近年來，兼顧RANS和大渦模擬(LES)方法的優勢，在不同的空間區域，選擇不同特征長度，能夠更好地預測分離流動的RANS/LES混合方法開始廣泛應用于固定翼飛行器流場和氣動力的預測[4-7]。在近壁面區域，該方法等同于RANS方法，在遠離物面的區域，計算模型為LES方法。該方法基于特征尺度判斷實現兩種計算模式的切換。與RANS方法相比，混合方法能夠較好地預測近壁面附近的分離流動，在較大范圍氣流分離發生時，預測效果更好；與LES方法相比，其網格數量要求大大降低，因而成為一種可接受的計算手段。目前RANS/LES混合方法發展出了脫體渦模擬(DES)、延遲脫體渦模擬(DDES)、IDDES(Improved DDES)等多個版本[8-11]，其中DDES方法設置了延遲函數，克服了網格加密時可能形成的附面層內LES方法過早啟動的問題。

隨著計算數值手段的發展進步，精確預測流場中的流動分離和尾跡發展成為進一步提高旋翼流場和氣動力計算精度的關鍵制約因素。為更好地模擬在較大的空間范圍內螺旋發展的旋翼尾跡，渦核半徑內需要布置多個網格點，這使得旋翼流場計算總體網格數量膨脹速度較快。在計算機硬件水平及計算程序并行能力迅速提升的大背景下，RANS/LES混合方法也開始應用在旋翼流場計算方面。特別是在2014—2017年的AIAA會議上，陸續有研究者針對旋翼的懸停流場預測問題，采用RANS/LES混合方法進行了數值計算研究[12-16]。采用混合方法時，與傳統的基于RANS方法的旋翼流場預測比較，新方法本身在分離流動的預測能力上有所加強，且預測工作一般采用了更大的網格數量，預測出了更多的流場細致結構。在國內，旋翼流場的預測一般基于RANS方法[17-20]，較少查閱到基于RANS/LES混合方法的旋翼流場計算文獻。

鑒于此，本文采用DDES方法進行了旋翼流場計算研究。計算網格系統構建方面，為有效控制計算規模，用于尾跡捕捉的背景網格使用了自適應直角網格，并根據特定的策略控制網格分布。為提高計算速度，采用了分布式并行的計算方法，無論是在流場求解還是在網格組裝等方面，均實現了并行化處理。基于這些方法，對旋翼多個狀態下的流場和氣動力進行了計算，這些計算工況對應了槳葉表面基本為附著流動，槳葉表面附近空間發生小范圍和大面積的氣流分離情況。討論了對于旋翼流場和氣動力預測問題，不同工況下計算結果與傳統RANS方法的計算差異，比較的內容既包括流場結構，也包含宏觀的氣動力差異。計算結果體現出RANS/LES混合方法在槳葉表面附近存在氣流分離情況下的特定優勢，采用此方法在一定程度上可以提升旋翼流場和氣動力的計算精度。

1 計算網格系統

為有效捕捉旋翼槳葉表面附近流動和遠場尾跡發展，在本文的嵌套網格系統中，貼體網格部分使用了純六面體網格，以提高網格質量，背景網格使用了可自適應的直角網格。

背景直角網格按照以下步驟生成：

步驟1根據預定的流動區域外場大小，各方向均勻布置網格點，生成一個較為稀疏的初始網格。

步驟2根據槳葉運動掃掠過的空間區域范圍，進行一次網格加密。

步驟3對預估的尾跡區域空間網格，進行適度加密調整。

步驟4根據流場計算結果，使用渦量探針，加密流場中渦量集中區域的網格。

采用這樣的方法，貼體與背景網格插值單元的尺度比較接近，有利于維持數值格式的精度；網格自適應的第4步可以多次逐級進行，用于有效捕捉遠場尾跡發展，且使得計算網格的數量可控。

在直角網格每次自適應加密后，還要進行“洞”和“層差”的判斷，消除網格中的不光滑區域。實際旋翼流場計算時，可以根據槳葉的弦長預估可能的渦核的尺寸，從而確定初始和最終的網格加密后尺度。

在數據格式構造方面，均以非結構形式對貼體和背景網格進行數據存儲。由于貼體和背景網格采用同一數據結構[21]，在空間格式構建等方面實現了統一處理，計算程序得到了極大的簡化。

2 流場數值計算方法

2.1 控制方程形式

采用格心格式的有限體積法，地面坐標系下積分形式的雷諾平均Navier-Stokes方程可以寫為

(1)

式中：W為守恒變量；S為積分面面積；Ω為控制體體積；Fc和Fv分別為對流和黏性通量；p為壓強；ρ為密度；V為對流速度；u、v、w為3個方向速度；nx、ny、nz為網格交界面單位外法矢；H為單位質量總焓；τ為應力；T為溫度；K為傳熱系數。

不考慮與轉捩有關項，Spalart-Allmaras湍流模型[22]可以寫為

(2)

fw定義為

g=r+cw2(r6-r)

渦黏性υt定義為

模型常數定義為：cb1=0.135 5；σ=2/3；cb2=0.622；Von Karman常數κ=0.41；cw1=cb1/κ2+(1+cb2)/σ；cw2=0.3；cw3=2.0；cυ1=7.1。 采用RANS/LES混合方法時，長度尺度定義為

(3)

式中：d為空間網格到最近壁面的距離；Δ=max(Δx,Δy,Δz)為當前單元和它所有鄰居單元中心間最大距離；CDES=0.65。

對于典型的RANS網格，在近壁面區域，RANS/LES混合模型中的RANS模式被激活。因為在該區域，壁面法向距離通常小于各方向網格間距的最大值。如果流向和展向的網格尺度在近壁面區域和法向尺度近似，LES模型則在邊界層內被激活，導致邊界層過早分離。為克服這一問題，Spalart等又提出了DDES模型，長度尺度改寫為

(4)

fd=1-tanh(8rd)3

2.2 求解方法

非定常流場計算方面，主控和湍流方程以松耦合的方式，采用雙時間[23]結合LU-SGS[24]方法進行時間推進，物理時間步長選擇為旋翼旋轉周期的1/1 440。在最終的計算網格上，完成11個旋轉周期數值模擬后，再進行3個周期的計算取結果進行平均，獲得最后的計算結果。

空間格式方面，采用二階迎風格式[25]計算無黏通量。由于網格采用無結構化處理，物理量的梯度需要采用重構方法實現(見文獻[26])，黏性通量采用中心格式進行計算。

3 數據并行方法

3.1 網格分區方法

分布式并行計算時，首先要對網格進行分區，該項工作在網格每次自適應調整后進行一次，使得各處理器上的負載盡量均衡。

旋翼流場計算采用多塊嵌套網格系統，這里按照處理器線程數對每塊網格進行分區，網格按照非結構的方式進行存儲，因而可以使用metis[27]工具來進行網格的切分。

設網格塊的標號為B1,B2,…;網格B1分區后的小塊網格標號分別為B1(1),B1(2)，…;網格B2分區后的小塊網格標號分別為B2(1),B2(2)，…。在線程1上分布的網格塊為B1(1)，B2(1)，…，在線程N上分布的網格塊為B1(N)，B2(N)，…。網格分區完成后，為構建數值格式方便，每塊網格的分區邊界單元向外擴展2層虛擬單元，形成處理器之間網格信息傳遞邊界。

3.2 網格組裝分布式并行方法

網格組裝工作分為洞切割、貢獻單元搜索和數據傳遞區域構建等幾步進行。為增加整個流場解算的并行段比重，提高計算效率，這里對這些步驟進行并行化處理。

洞切割方面，首先在所有網格塊中依據距離物面的遠近，根據已有網格面形成一個洞構造表面，并把該構造表面的網格信息廣播到所有的處理器上。這樣在各個處理器上，均可按照該構造面進行洞切割。

貢獻單元搜索方面，對每個處理器上洞切割后形成的插值邊界，均發送給其他處理器，進行貢獻單元搜索，然后返回相應結果。這樣就實現了貢獻單元搜索的并行化。

對每個處理器，在其他多個處理器上進行貢獻單元搜索的結果回收后，在自身處理器上進行結果匯總，并最終確定正確的數據交換關系。

3.3 流場計算相關分布式并行方法

在分布式并行情況下，采用隱式時間推進時，由于在每個處理器上無法獲得全場的流動信息更新情況，因而原始的LU-SGS方法需要改進。這種改進主要體現在各處理器分區邊界的流動變量更新計算方面。實際計算時，可以先用退化的方式(雅克比迭代)計算出傳值邊界元的流動變量更新量，然后各個處理器進行數據交換，再進行局部的LU-SGS迭代求解，具體的執行方法見文獻[28]。

對于一些流場計算中需要的幾何信息，也可以進行并行化處理。如壁面距離確定方面，可以通過壁面坐標的廣播實現在每個處理器上的并行計算。

4 計算方法驗證

4.1 翼型振蕩氣動力計算

為驗證本文方法在非定常氣動力計算方面的能力，選取了NACA0015翼型強迫振蕩狀態進行數值模擬。翼型弦長為0.304 8 m，來流馬赫數為0.29。翼型迎角α隨時間t的變化規律為

α=4.0+4.2sin(20πt)

(5)

圖1給出了計算得到的升力系數CL、阻力系數CD和力矩系數Cm隨迎角α變化的遲滯曲線與試驗值[29]的對比。計算結果表明，采用RANS和DDES方法計算出的升力與試驗結果符合得更好，阻力的結果相對差些，這可能與湍流模型的選取等多種因素有關；兩者計算的結果均能夠有效模擬出氣動力隨迎角變化的趨勢，其相互間差異隨著翼型迎角的增大而增大。

圖1 NACA0015翼型振蕩狀態非定常氣動特性Fig.1 Unsteady aerodynamic characteristics of oscillation of NACA0015 airfoil

4.2 槳尖附近區域激波誘導氣流分離

為驗證本文方法在旋翼氣動力計算方面的能力，選擇了有試驗值可供對比的“Caradonna & Tung旋翼”。該旋翼槳葉數為2片，展長為1.143 m，展弦比為6，無扭轉尖削，槳葉截面段翼型為NACA0012。選取的計算狀態為槳尖馬赫數0.794，12°總距角。計算發現，在該狀態下，在靠近槳尖的槳葉上表面附近空間區域出現較強激波，在激波后方形成一定范圍的氣流分離。本文對比了采用非定常RANS和DDES方法的計算結果。圖2為計算得到的不同槳葉方位角ψ時槳葉上表面壓力的對比。可以看出，采用RANS方法時，計算的表面壓力幾乎不隨時間(槳葉方位角)變化，而采用DDES方法時，則預測到了槳葉上表面壓力分布時變特性(靠近槳尖區域)。這體現出DDES方法在分離流動方面的預測能力。

圖2 槳葉上表面壓力對比(左：RANS，右：DDES)Fig.2 Comparison of blade upper surface pressure(left: RANS; right: DDES)

圖3 槳葉表面壓力系數計算值與試驗值的對比Fig.3 Comparison of calculated and experimental surface pressure coefficients of blade

圖3為計算得到的靠近槳尖的部分截面時均表面壓力系數Cp計算值和試驗值[30]的對比，r為槳葉截面段所在位置半徑，R為旋翼半徑,c為弦長。可以發現，DDES方法預測了激波后方的非定常分離流動，激波后上表面的壓力值更小；相對RANS方法，DDES方法預測的激波位置也更靠近槳葉前緣。對于槳葉內側截面，由于其流動為附著流動，兩種方法預測結果一致，這里不再給出。

5 算例分析

5.1 雙旋翼干擾流場

本算例中，計算模型為南京航空航天大學的一副試驗共軸雙旋翼[31]，該旋翼的上下旋翼相距0.175 m，各由兩片槳葉組成，槳葉半徑為0.945 m，無扭轉，翼型為NACA0012，弦長為0.076 m，槳葉r=0.95R處向外有尖削，尖削比為1/3，試驗狀態槳尖速度為118 m/s，上下旋翼槳葉槳距角分別為9°和10.27°，預錐角為5°。

計算網格如圖4所示，旋翼旋轉軸方向為Y軸，左圖是平行Y軸的網格截面圖，右圖為旋翼下方垂直Y軸的網格截面圖。網格尺度方面，旋翼下方尾跡區的網格尺度L設定為0.1C(C為槳葉根弦長)，按照渦量自適應的網格尺度為0.05C，空間網格尺度隨著與旋翼距離增加逐漸放寬。

貼體網格單元為純六面體形式，每片槳葉貼體網格單元數約為430萬，背景直角網格單元數約為6 900萬。計算采用700個CPU并行，每次流場內迭代時間約為2.8 s。在該計算狀態下，由于槳葉截面相對來流速度較低，迎角不高，槳葉表面觀察不到明顯的分離流動。

圖4 共軸旋翼流場計算網格示意圖Fig.4 Schematic of grid for coaxial rotor flowfield calculation

圖5 時均軸向誘導速度Fig.5 Time-averaged axial induced velocity

圖5給出了計算得到的槳盤下方不同截面處(圖中y為槳盤下方截面的位置)軸向誘導速度時均值Vy與試驗值的對比。在有試驗結果的區域，計算與試驗值較吻合(兩種計算方法得到的軸向誘導速度重合)。采用DDES方法與RANS方法的計算差別主要體現在靠近槳根和槳尖的區域，在靠近槳根區域，DDES方法預測出稍大一些的下洗速度。這可能是由于在槳葉靠近兩端的區域，氣流繞過槳葉，展向流動效應加強，流動復雜度增加。

圖6為計算得到的下旋翼槳葉整體氣動力(拉力系數CT和扭矩系數CQ)的比較。在槳葉方位角ψ=100°～130°(定義上下旋翼俯視重合時的槳葉方位角為0°)，下旋翼受到上旋翼尾跡影響較為嚴重，無論是槳葉拉力和扭矩都出現波動。兩種方法計算得到的拉力差別體現在峰值大小上，而扭矩差別僅體現在出現氣動干擾的方位上。總體來說，兩種方法計算出的宏觀氣動力差別較小。這可能與槳葉表面空間附近無明顯分離有關。對于上旋翼槳葉，兩種計算方法得到的結果基本重合，這里不再給出。

圖6 下旋翼的拉力和扭矩系數Fig.6 Trust and torque coefficients of lower rotor

圖7 截面渦量等值圖Fig.7 Contour of vorticity in cutting plane

圖7為采用DDES方法計算得到流場空間特征截面上的渦量等值圖，上下旋翼的尾跡在經過下旋翼的下方后，開始融合并一起向下發展。采用本文的網格系統和數值方法，計算得到了較為精細的流動結構。

5.2 下降狀態旋翼流場

下降狀態是直升機的一個重要工作狀態，該狀態下流場結構更為復雜，容易發生氣動干擾現象，引發氣動和噪聲問題。

在該算例中，仍然選擇了“Caradonna & Tung旋翼”構型作為研究對象。計算狀態為槳尖馬赫數0.439， 8°總距角。軸向誘導速度Vy分別選擇為20 m/s和40 m/s。

貼體網格單元為純六面體形式，每片槳葉貼體網格單元數約為300萬，背景直角網格單元數約為3 400萬。計算采用700個CPU并行，每次流場內迭代時間約為1.2 s。

圖8為下降速度20 m/s時，采用DDES方法計算得到的不同槳葉方位角時流場特征截面上的渦量等值圖，可以發現計算的渦量具備明顯的時變特征。需要注意的是，由于來流方向為槳葉下方，在槳葉的向下作用力和向上來流共同作用下，槳葉脫出的尾跡在經過180°方位角時(到達另一片槳葉空間位置附近)并不能迅速遠離槳盤平面。

圖9為下降速度20 m/s時，采用DDES方法和RANS方法計算得到的特征截面(x=0 m)上的流線圖(0°槳葉方位角)。可以發現，采用RANS方法計算的流場具備良好的對稱性，而采用DDES方法結果則存在明顯的左右不對稱。

圖8 不同槳葉方位角時的流場截面渦量等值圖(x=0 m)Fig.8 Contour of vorticity in flowfields cutting plane with different blade azimuthal angles (x=0 m)

圖9 截面流線(x=0 m)Fig.9 Streamlines in cutting plane (x=0 m)

由于DDES方法相比RANS方法在分離流動與旋渦結構捕捉能力方面更強些，因而其計算結果在流場中存在更多復雜渦結構的情況下更接近于真實物理情況。

圖10給出了該狀態下，特征截面上的渦捕捉情況對比。可見DDES方法解析了更多的流場渦結構，且其集中渦量的計算強度普遍高于RANS方法結果。

圖11和圖12分別為Vy=20，40 m/s，槳葉方位角為0°時，采用DDES方法和RANS方法計算得到的不同槳葉展向截面上的流線圖(已經處理成相對槳葉截面速度)。可以發現在不同槳葉截面段上，Vy=20 m/s時，兩種方法計算流動情況類似，即均無明顯分離流動現象。在Vy=40 m/s時，兩種方法均計及槳葉表面大范圍的分離流動，但采用DDES方法和RANS方法計算出不同的流動狀態，尤其是靠近槳尖附近的截面，兩種方法計算得到的結果差異更大。

圖13給出了采用不同方法計算得到的槳葉升力展向變化情況。在Vy=20 m/s狀態下，RANS方法計算出時均的槳葉升力大于DDES方法的結果。可以發現，即使槳葉截面段上觀察不到明顯分離流動(圖11)，RANS和DDES方法計算的結果差異也很明顯。分析其原因可能是，在該狀態下，旋翼整體尾跡和槳葉干擾作用較強，一方面旋翼仍提供正向升力，另一方面，旋翼尾跡向上方發展，不能迅速遠離槳盤平面，由于兩種預測方法計算出的尾跡發展情況不同，因而尾跡與槳葉的作用情況也必然不同。還可以看出，在該狀態下，兩種方法預測的氣動力都隨槳葉方位角變化明顯，這表明在這種情況下，即便是軸流入流狀態，旋翼的氣動力也具備較強的時變特性。

圖10 截面渦量 (x=0 m)Fig.10 Vorticity in cutting plane (x=0 m)

圖11 Vy=20 m/s時槳葉截面流線(左：DDES，右：RANS)Fig.11 Streamlines in blade cross-sections for Vy=20 m/s (left: DDES; right: RANS)

圖12 Vy=40 m/s時槳葉截面流線(左：DDES，右：RANS)Fig.12 Streamlines in blade cross-sections for Vy=40 m/s (left: DDES; right: RANS)

圖13 槳葉展向升力比較Fig.13 Comparison of spanwise lift of blade

在Vy=40 m/s狀態下，RANS計算出時均的槳葉升力小于DDES方法的結果。這應該與兩種方法預測出的槳葉表面附近的分離形式差異有關。從圖13還可以看出，采用RANS方法計算出的氣動力相對穩定，槳葉不同方位角時的氣動力比較接近。采用DDES方法計算的結果則體現出較強的時間特性，這體現出DDES方法在物體表面附近存在強分離時的計算特性。

6 結 論

在可自適應嵌套網格系統下，基于DDES方法，進行了不同狀態下復雜旋翼的流場和氣動力計算，并與相應的RANS結果進行比較。

1) 在槳葉表面無明顯分離，且不存在較強的旋翼/尾跡間相互作用情況下，DDES和RANS方法計算出的槳葉宏觀氣動力方面差異較小。

2) 旋翼/尾跡間存在較強的相互作用情況下，即便槳葉截面流場觀察不出明顯分離流動，兩種方法計算結果也存在明顯差異，這可能與尾跡預測的強度等因素有關。在該情況下，即便是在軸流入流狀態，無論是RANS方法還是DDES方法，計算出的旋翼氣動力都體現出較強的時變特性。

3) 在槳葉表面附近空間存在大范圍分離流動的情況下，采用DDES方法計算出更多的流場結構，氣動力的時變差異也更為顯著，體現出經典的混合RANS/LES方法的計算優勢。

4) 本文采用的網格系統和計算方法，能夠較好地預測流場中的流動細節和旋翼氣動力。

值得注意的是，對存在復雜分離流動的旋翼流場預測來說，除去本文考慮的DDES方法研究外，提高計算格式的數值精度、網格數量和質量，進行湍流模型的選擇及轉捩研究等都是進一步提升旋翼流場計算精度的關鍵因素，也是未來工作的發展方向。