基于S-A與SSG/LRR-ω兩種湍流模型的CHN-T1標模計算與分析

聶勝陽, 王 垠, 劉志強, 金 朋, 焦 瑾

(1. 西安理工大學 土木建筑工程學院, 西安 710058; 2. 西安航空學院 飛行器學院, 西安 710077)

0 引 言

計算流體力學(CFD)和實驗流體力學是航空領域特別倚重的兩種研究航空器氣動特性的手段。實驗流體力學方法經典但昂貴，計算流體力學需要實驗流體力學技術對計算流體力學程序進行充分的驗證與確認。國際上對CFD的驗證與確認工作已經進行了很多，如已經完成的六次AIAA 阻力預測會議(https://aiaa-dpw.larc.nasa.gov/)，三次高升力計算會議(https://hiliftpw.larc.nasa.gov/)，用可重復的高水平實驗，來驗證和確認不同研究機構開發的CFD程序，提供了大量的寶貴經驗并為CFD程序的改進提供了非常有價值的參考意見。相關工作在中國并未大規模地、有系統地、多方參與地組織起來，在這個背景下，中國空氣動力學會主辦的第一屆航空CFD可信度研討會(AeCW-1)，以CHN-T1單通道客機標模[1]為研究對象，開展以CFD程序驗證和確認為目標的系統研究工作，邀請全國CFD開發者參與進來，標模試驗見參考文獻[2]。

工程計算上最常用的計算方法是求解雷諾平均Navier-Stokes方程(RANS)，湍流行為用湍流模型模擬，如常用的 Spalart-Alamars模型(S-A)[3]和Menter[4]開發的SST k-ω兩方程模型。兩種模型都是基于渦黏性假設構造而成，對湍流非線性特征明顯的流動，如旋渦流動、流線彎曲、二次渦等，因超出渦黏性假設成立的范圍，模擬能力不足，進而出現了很多針對這些不同流動類型的修正，如Shur等[5]提出的旋渦修正，Mani等[6]提出的流線彎曲修正，Spalart等[7]和Mani等[8]提出的二次本構關系(Quadratic Constitutive Relation，QCR)對S-A模型加入湍流各向異性效應的修正等。以上的修正方法提升了兩種常用湍流模型的應用范圍，尤其是QCR修正，在第四次AIAA會議[9-11]以后，有QCR修正的湍流模型因對翼身聯結處的角渦具有很好的預測能力而大受歡迎。但是，不添加任何修正，天然考慮流線彎曲、湍流各項異性等特點的雷諾應力模型(RSM)，由于控制方程數目多、模型構造復雜、計算耗費資源較多、且存在著計算魯棒性的問題，未得到工業界足夠的重視，工程應用較少。對于復雜的外形的工程應用計算，僅見Eisfeld開發的SSG/LRR-ω雷諾應力模型[12-14](RSM)和Wilcox的Stress-ω模型[15]的有限應用，如NASA的Lee-Rausch等[16]用前述兩種RSM模型計算了第二屆AIAA高升力會議的構型。目前文獻[12]中提出的SSG/LRR-ω雷諾應力模型，對激波誘導分離流動有較好的預測能力，正取得了工業界的重視。一方面，Rumsey等[17]發現，在ONERA M6機翼上，SSG/LRR-ω雷諾應力模型在激波位置和分離區的大小的預測上表現穩健；DLR的Keye等[11]比較了SSG/LRR-ω模型和添加QCR修正的渦黏性模型對翼身結合處的角渦的預測，發現SSG/LRR計算效果和引入QCR的渦黏性模型相當，顯示出該湍流模型是一個很有潛力進行工程計算的高級湍流模型；同時，RSM在模型的構造上，主要利用速度一階導數信息，因此模型對網格的敏感性比k-ω類方法要低，達到網格收斂所需網格數目較少。另一方面，NASA的Rumsey[17]和Dudek[18]等在一些其他的常見流動類型的研究表明，相對于使用對應修正的渦黏性模型，RSM模型并沒有明顯的壓倒性的優勢。綜上，對RSM的研究和未來發展仍有不少爭議，因此NASA的CFD2030遠景[19]也提到，需要對RSM模型需進行更充分的研究論證，以進一步開發該模型的潛力。近年來，對雷諾應力模型的開發進展緩慢，但是RSM的工程應用仍在不停發展。Togiti等[20]開發的基于g變量的SSG/LRR雷諾應力模型的變種，具有更好的數值魯棒性。Nie等[21-22]將RSM和轉捩模型耦合，開發出具有自動轉捩能力的雷諾應力模型。國內對雷諾應力模型的開發和復雜外形上的計算應用較少。聶勝陽等[23]研究過RSM在ONERA M6機翼上的激波分離流中的預測能力以及聶勝陽等[24]研究RSM在旋渦流動中的預測能力，董義道等[25]對SSG/LRR-ω雷諾應力模型在典型算例和DLR-F6上進行了計算分析，都論證了RSM具有較好的復雜湍流的計算能力。本次CHN-T1模型上的角渦流動很小，因此本文分別使用了標準Spalart-Alamars模型(S-A)和SSG/LRR-ω雷諾應力模型對CHN-T1的所有工況進行計算,用來考察渦黏性模型和不使用渦黏性假設的雷諾應力模型在CHN-T1運輸機標模上的模擬能力，對所得結果進行詳細的分析，為湍流模型的開發和應用評估提供參考。

1 計算狀態和計算網格

AeCW-1研討會組委會提供的基準非結構網格[26-27]，采用Pointwise軟件生成。網格生成策略是用block-to-block方法先生成表面網格，然后外推三棱柱。在三棱柱外的空間用陣面推進法填充剩余空間。粗中細網格由block上的點數增加1.414倍得到。粗中細網格的三棱柱層的高度一致。最終的粗網格有230萬個網格點，650萬網格單元，其中三棱柱有350萬；中等網格有595萬網格點，1680萬網格單元，其中三棱柱有840萬；密網格有1725萬個網格點，4950萬個網格單元，其中三棱柱有2648萬。Config1的中等網格和細網格在飛機上的分布如圖1 所示。Config2和Config3上的計算網格與Config1的中等網格一致。計算構型的信息和計算狀態分別參見表1和表2。

表1 計算構型匯總Table 1 Summary of the configurations

表2 計算狀態匯總Table 2 Summary of the computation settings

圖1 case1上中等網格和細網格在機翼機身結合處附近的表面網格分布Fig.1 Illustration of the grid distribution on medium mesh(left) and fine mesh (right)on the wing-body junction for case1

2 計算方法與湍流模型

2.1 計算方法

本文使用自研的基于非結構網格的流場求解器——湍流模型開發平臺(UTMDP)，對AeCW-1研討會組委會提供的基準非結構網格進行了計算。流場控制方程為RANS。湍流模型對流項的離散格式為二階Roe格式；無黏通量的離散格式為中心格式，添加適量的矩陣人工黏性，二階人工黏性系數為0.5，四階黏性系數為1/128，不采用墑修正；用多重網格技術加速收斂，LU-SGS隱式時間推進，定常計算，CFL數為5.0，計算分區為24塊用于并行計算。

2.2 湍流模型

湍流模型選擇計算很魯棒的標準S-A模型和SSG/LRR-ω模型。S-A模型采用去掉轉捩控制項的SA-noft2模型，在本文中用S-A簡稱。SSG/LRR-ω模型若采用通用擴散模型(Generalized Gradient Diffusion Model，GGDH)無法得到穩定的計算結果，所以擴散項采用擴散項模型，即SSG/ LRR-RSM-w2012-SD，在本文中用RSM簡稱。遠場邊界條的自由來流湍流度設為0.2%，湍流黏性比設為0.001。計算發現，對于RSM，需要對輸運ω的最小值進行一定的限制，以獲得穩定的計算結果。SA-noft2模型和SSG/ LRR-RSM-w2012-SD模型的具體信息可以參看NASA Turbulence Model Resource(TMR)[28]。

2.3 計算收斂性

采用S-A模型計算，在Config1，Config2和Config3構型上都獲得了收斂的解。 RSM在Config1上可以得到穩定收斂解。如圖2所示，在Config1構型上，即使在密網格計算，使用S-A模型計算5000步后，密度的平均殘值的對數可以下降到-6，力系數接近收斂。使用RSM計算則需要更多的迭代步數才能使密度平均殘值下降到S-A模型計算5000步后的水平。

而RSM在有尾支撐的構型上(Config2和Config3)，密度平均殘值的對數下降到-2之后，就難以進一步下降，力系數震蕩。在較小迎角的工況，升阻力系數的變化在1 count以內。當迎角大于3.5°，升阻力系數變化范圍擴大，但是升力系數的變化不超過10 counts，阻力系數不超過5 counts，典型的RSM計算的力系數隨迭代步數的變化見圖3。此時迎角較大，機翼上出現了激波誘導的較大的分離，飛機處在抖振邊界附近，力系數會出現震蕩。由于在Config1上RSM可以得到收斂的結果，因此引起Config2和Config3上小迎角狀工況的計算仍難以收斂的原因應來源于尾支撐的存在。

圖4展示了使用S-A和RSM模型計算得到的尾撐上的表面極限流線分布。組委會提供的非結構網格，基于的幾何模型是對尾撐進行了截斷處理的模型，在尾撐下游的流動出現大分離。S-A模型計算的分離渦是一個穩定的大渦。而RSM解析出更復雜的旋渦結構，進行非定常地脫落，導致計算難以收斂，力系數出現小幅震蕩。

圖2 在Config1上使用密網格計算2°狀態下的收斂歷程Fig.2 Convergence history of the computation using fine mesh at 2° on Config1

圖3 4°工況下RSM在Config2上的力系數隨迭代步數變化Fig.3 Force coefficient varies with iteration steps computed by RSM on Config2 at angle of attack 4°

圖4 Config2上0°工況下S-A和RSM計算的尾撐上的表面極限流線圖Fig.4 Surface skin-friction lines computed by S-A and RSM on Config2 at angle of attack 0°

3 結果分析

3.1 網格收斂性研究

表3是采用RSM模型和S-A模型對Config1做的網格收斂性研究，參考了文獻[2]中的實驗結果和文獻[29]中的基于104億網格的RANS計算結果。計算表明，在密網格上仍未達到網格收斂，收斂趨勢比較明確,通過外插到網格尺度為0的值也在表中列出。S-A模型和RSM模型具有相似的網格收斂過程，但RSM模型具有更快的網格收斂趨勢，如圖5所示。RSM計算出的定升力迎角和S-A模型計算的定升力迎角有0.1°的差距。和104億網格計算結果對比發現，RSM模型更接近高密度網格的模擬結果，只是阻力偏大。

表3 Case1上使用S-A模型和RSM的計算結果匯總Table 3 Summary of results by S-A and RSM

圖5 定升力條件下的阻力系數隨網格單元數目的變化圖Fig.5 Drag coefficient with different grid types and turbulence models for case1

圖6展示的是機翼上部分典型站位的壓力分布，取展向站位17%、70%和95%處的機翼以及平尾上3個站位進行比較，左側為S-A模型在粗中細網格上的計算結果，右側為RSM計算結果。隨著網格密度的增加，激波的分辨率提高，寬度越來越窄,吸力峰值稍微升高，總體變化幅度很小。在翼梢附近的弱激波，隨著網格密度的增加，越來越明顯。水平尾翼上的壓力吸力峰值隨著網格密度增加而小幅增加，其他位置上的壓力分布受網格分布的影響很小。

圖6 Case1上由S-A和RSM計算的不同展位的壓力分布Fig.6 Pressure distribution at different spanwise locations with S-A and RSM models for Case1

圖7展示的是S-A模型和RSM在粗網格和密網格上計算的壓力分布的比較。主要的區別是：在上表面，RSM計算的激波位置較靠近上游，在激波位置的上游，RSM計算的負壓大于S-A，在激波下游，RSM計算的負壓小于S-A模型的計算結果；在下表面的機翼后加載區域，RSM計算的負壓小于S-A模型的計算結果。這些細微差別不隨網格密度增加而改變，表明這些差別來源于湍流模型本身對湍流分布和強度的解析不同所致。

圖7 S-A和RSM在粗網格和密網格上的計算結果對比Fig.7 Pressure coefficients calculated by coarse and fine meshes with S-A and RSM models for Case1

3.2 抖振特性研究

3.2.1 計算結果綜合分析

抖振特性的計算研究包含對Config1、Config2和Config3三種構型的計算研究。由于在網格收斂性研究中發現，中等網格上的計算結果和密網格上的計算結果有些微細小差別，但是在抖振邊界附近，流動中出現較多復雜分離湍流結構，因此對Case2a的計算，除了按要求在中等網格上進行了計算，還增加了在密網格上的計算分析，用來研究在抖振邊界附近工況下網格的影響。圖8使用S-A模型和RSM計算的力系數隨迎角的變化圖。對升力系數而言，S-A模型在Config1、Config2和Config3上計算的升力系數較實驗數據都稍微偏大一些。Config2和Config3考慮了尾撐效應，計算的升力系數在線性區更大一些，在非線性區，Config3上計算的升力系數更接近實驗數據。對Config1構型使用中等網格的計算，升力系數在線性區和非線性區上沒有明顯區別。在非線性區，中等網格上的計算結果更接近實驗數據，密網格上計算的升力系數明顯偏大。對RSM而言，在不同構型上計算的升力系數都很接近實驗數據。主要差別在非線性區，Config1上計算的升力系數最大，Config3上計算的升力系數最小，都明顯小于實驗數據。但是Config1上使用密網格計算的升力系數更接近實驗數據。對理想阻力系數來說，所有構型上計算的阻力系數都小于實驗觀測。Config1上計算的阻力系數明顯高于Config2和Config3；RSM計算結果和實驗值的差距大于S-A模型計算的差距；在Config1上的S-A模型和RSM計算的阻力系數更接近實驗數據。而俯仰力矩系數的計算，S-A模型和RSM的計算結果都和實驗數據差別較大。在Config1構型上，當迎角小于3.75°時，網格密度影響不大，只有高于3.75°時，密網格上計算的俯仰力矩系數略小于中等網格上的計算結果；Config2上的計算結果和實驗數據的偏差最小；在較大迎角的工況，考慮了尾撐效應和機翼靜氣彈變形的Config3構型，其上計算的俯仰力矩系數和實驗數據的差別大于沒有考慮機翼靜氣彈變形的Config2，且RSM在Config3上計算的俯仰力矩系數在較大迎角時，產生的偏差甚至大于在Config1上的結果；在Config2和Config3上，RSM計算的結果都較S-A模型計算的結果偏離實驗數據。

圖8 S-A 模型和RSM計算的升力，理想阻力(CD_ideal=CD-CL/(2πAR))和俯仰力矩系數隨迎角的變化Fig.8 CL, CD_ideal and Cm vs. angle of attack computed by S-A and RSM models for Case2

S-A和RSM模型在同一構型上，計算出的流動現象有一些差別。當迎角小于3.5°時，機翼上沒有激波誘導的分離；迎角大于3.5°時，出現激波誘導分離，且分離區隨迎角增加而增大。當迎角增加到4.5°時，機翼中段出現了大范圍的激波分離流。RSM計算的激波誘導分離小于S-A模型計算結果，分離區也小于S-A模型計算結果。圖9展示了S-A模型和RSM模型在Config1上計算的摩阻系數分布，計算工況從左到右分別是迎角為3°、3.75°和4.25°。在3.75°工況上，RSM計算的激波誘導分離明顯小于S-A模型的計算結果。隨著迎角的增大，分離區擴大，RSM計算的分離區始終小于S-A計算的分離區。

圖9 S-A和RSM表面摩阻分布和表面極限流線Fig.9 Surface skin-friction coefficient and stream lines computed by S-A and RSM models

CHN-T1標模使用的機身限制了翼根處的角渦的發展，角渦大小有限。隨著迎角的增加，翼根處的角渦增大。圖10顯示的是S-A模型和RSM計算的翼根處角渦在4.25°時細節的比較,RSM計算的角渦較大一些，網格加密后，角渦增大，但是相對機翼尺度而言，角渦一直很小，對氣動特性的影響可忽略。

圖11是Config1構型上機翼展向50%和95%站位處計算的壓力分布比較，分別使用中等網格和密網格，計算工況是3.75°。可以看出，S-A模型計算的激波位置靠后，RSM計算的激波靠前。在翼梢附近，RSM計算的激波位置仍稍微在上游。網格加密對激波位置無影響，激波隨著網格加密寬度變窄。

圖10 Case2a在中等網格是由S-A和 RSM計算的角渦Fig.10 Wing root trailing edge separation bubble computed by S-A and RSM models for Case2a

(a) η=50%

(b) η=95%

3.2.2 尾撐效應和機翼靜氣彈效應

由于S-A和RSM在CHN-T1標模上計算的現象，除了誘導分離區的大小有明顯差別外，其他的現象都和接近。因此研究尾撐效應影響的時候，只展示S-A模型的的計算結果。圖12展示了機翼和平尾部分站位上的壓力分布的比較。可以看出，在機翼上，壓力分布幾乎完全一致，尾撐不影響機翼上的流動。但在平尾上，壓力分布出現較大的差別，表現為Config2上平尾的上下壓力差變小，吸力峰值下降，越靠近平尾翼根尾撐的影響越大。CHN-T1構型的平尾產生負的升力，貢獻給抬頭力矩。由于平尾的上下壓力差變小，平尾對抬頭力矩的貢獻變小，因而計算的總的俯仰力矩系數降低。計算的結果和Config1相比，更接近實驗數據。

(a) 機翼

(b) 平尾

將有尾撐的Config2和Config3比較，發現考慮機翼靜氣彈變形的Config3構型上，機翼壓力分布出現較大的變化，如圖13(a)所示。靜氣彈變形帶來了激波前吸力下降，激波后移，越靠近翼梢越明顯。而平尾上的壓力分布受靜氣彈影響較小，吸力峰值微微增加，上下壓力差變大。而由于吸力下降，激波變弱，激波誘導的分離也變小，如圖14所示，在3.75°工況下，Config3上激波誘導的分離區明顯變小。

機翼前緣附近吸力下降，會引起俯仰力矩系數下降。而平尾上的上下壓力差產生的負升力稍微增加，帶來總的俯仰力矩系數的增加。然而，計算的力矩系數，如圖8(c)所示，俯仰力矩系數因靜氣彈效應而增加。導致這一現象的原因是尾撐和靜氣彈效應對機身尾部作用所引起的俯仰力矩的變化更復雜。圖15展示在三種構型上的機身尾部，由S-A模型計算的表面流線分布和表面摩擦系數在-2°工況的分布。可以看出，因為尾撐的存在，機身尾部上的流動出現顯著的變化，形成一個低速旋渦。低速旋渦在考慮靜氣彈效應的Config3上變小，位置向上游移動。由于靜氣彈效應，流過機翼的氣流被機翼上表面加速的程度稍微下降，被加速的氣流從水平尾翼下流過尾撐，引起尾撐和機身間的旋渦因注入能量下降而旋渦尺度和位置都發生了變化。該旋渦位置和能量的變化，因距俯仰力矩參考點較遠，帶來俯仰力矩的較大變化。由RSM計算的結果也有同樣的現象，如圖16所示。只是RSM計算的旋渦要小很多。隨著迎角的增加，如圖17和18，低速旋渦的尺寸縮小。對于考慮靜氣彈效應，計算的力矩系數反而比Config2上更大，迎角越大，和實驗數據的差距越大。由于CHN-T1標模的尾撐相對尺寸較大，安裝位置離平尾很近，帶來平尾上壓力分布顯著變化，引起俯仰力矩的變化；同時，尾撐對機尾上流動的影響更復雜，誘導出一個分離旋渦，耦合靜氣彈效應，引起上游來流能量下降，從而導致旋渦位置變化和機尾上的壓力分布變化，進而帶來俯仰力矩系數的改變。綜合考慮S-A和RSM在Config2和Config3上計算的力矩系數和流動現象，可以得出，S-A模型計算出的平尾和機身和尾撐結合處的旋渦更大，可能更接近實驗現象。

(a) 機翼

(b) 平尾

(a) Config2 (b)Config3

圖15 S-A計算的Case2a 、Case2b和 Case2c在 -2°工況下表面摩阻和極限流線分布Fig.15 Surface skin-friction coefficient and streamlines caluated by S-A model at angle of attack -2° for Case2a, Case2b and Case2c

圖16 RSM計算的Case2a 、Case2b和Case2c在 -2°工況下表面摩阻和極限流線分布Fig.16 Surface skin-friction coefficient and streamlines caluated by RSM model at angle of attack -2° for Case2a, Case2b and Case2c

圖17 S-A計算的Case2a、Case2b和Case2c在 4.25°工況下表面摩阻和極限流線分布Fig.17 Surface skin-friction coefficient and streamlines caluated by S-A model at angle of attack 4.25° for Case2a, Case2b and Case2c

圖18 RSM計算的Case2a、Case2b和Case2c在4.25°工況下表面摩阻和極限流線分布Fig.18 Surface skin-friction coefficient and streamlines caluated by RSM model at angle of attack 4.25° for Case2a, Case2b and Case2c

3.3 雷諾數影響

表4展示的是計算的雷諾數對運輸機標模的影響。在低雷諾數時，達到CL=0.5的迎角比高雷諾數工況所需迎角高約0.25°。高雷諾數下的阻力系數明顯小于低雷諾數，主要是黏性阻力大幅下降。使用S-A模型和RSM計算的定升力系數下的阻力系數和所需迎角在高雷諾數時更接近一致，表明黏性作用在高雷諾數工況時影響下降。總體而言，高雷諾數時S-A計算結果更接近實驗測量，低雷諾數時RSM計算結果更接近。圖19展示的是機翼上典型站位的壓力分布的比較：高雷諾數時吸力峰值下降，但是機翼后加載翼段上的負壓增加。因此高雷諾數時帶來的力矩系數小于低雷諾數。

表4 Case3上S-A計算結果匯總Table 4 Summary of the computations with S-A model for Case3

圖19 Case3上機翼部分站位的壓力分布Fig.19 Pressure distributions on the wing for Case3

4 結論與展望

對本文重要的結論梳理如下：

(1) 對計算收斂性研究發現：S-A模型在所有構型上都具有更好的收斂性，而RSM僅在無尾撐的構型上，收斂性好，在其他構型上計算的力系數出現小幅震蕩。引起力系數震蕩的主要原因是組委會提供的非結構網格基于的幾何模型，其尾撐沿中間某一橫截面截斷, RSM在橫斷面后計算出復雜的非定常的旋渦系統，引起力系數的震蕩。

(2) 基于Config1構型的定升力系數網格收斂性研究表明：S-A和RSM模型都能趨向于網格收斂。RSM計算的定升力阻力系數更接近實驗數據，所需迎角也接近實測數據，比S-A更快趨向于網格收斂。但是計算的力矩系數和實測數據差距明顯。使用S-A模型和RSM計算的壓力分布很接近；兩種湍流模型計算的吸力峰值隨網格密度增加而微微增加；在最密的網格上，都可以分辨出翼梢附近的弱激波。

(3) 對Config1構型的計算表明：RSM計算的激波位置比S-A模型計算的位置更靠近上游，激波引起的分離區較S-A模型計算的小，且不隨網格密度增加而改變。

(4) 對抖振特性的計算表明：在Config1，Config2和Config3上 S-A和RSM計算出的升阻力系數和主要流動現象都很接近。較小的差別體現在：所有構型上，S-A計算的升力系數較實驗數據稍微大，RSM計算的升力系數更接近實驗數據，在考慮尾撐效應和靜氣彈效應后，計算阻力系數偏小。計算力矩系數和實驗數據的偏差，在考慮尾撐效應后改善；但是考慮靜氣彈效應后，力矩系數的偏差再一次增加。主要原因是尾撐和機尾的相互作用復雜，是復雜湍流問題，并且計算用的網格在此區域可能不夠密，準確捕捉機尾上的復雜湍流現象是一大難點.

(5) 對雷諾數效應研究表明：當雷諾數從Re=3.3×106增加Re=15×106， S-A模型和RSM計算的結果趨于一致，阻力下降，俯仰力矩系數也下降。

對本次CHN-T1標模進行實驗研究發現，尾撐所引起的效應帶來的干擾較大，流動現象復雜，使得數值模擬很難再現實驗結果，因此在未來類似構型的實驗，應適當改進尾撐與機身的連接方式，減少尾撐對機尾區域流動的影響。