基于非結構混合網格的CHN-T1標模氣動特性預測

張耀冰, 唐 靜, 陳江濤, 鄧有奇

(中國空氣動力研究與發展中心 計算空氣動力研究所， 綿陽 621000)

0 引 言

近些年來CFD取得了很大的進步，包括網格生成技術、流場求解和高性能計算機等，評估三維真實運輸機外形成為可能。CFD預測工業相關外形的氣動性能取得了長足的進步[1-4],隨之而來的CFD可信度、數值計算結果驗證和確認以及數值結果的不確定度分析在國內外CFD業界受到廣泛關注和高度重視。AIAA為此專門召開了六次阻力預測會議[5-7]和三次高升力預測會議[8]，邀請世界范圍內的大學、研究所和工業部門共同評估當前CFD方法預測運輸類飛行器氣動力和力矩系數的發展水平，明確未來的發展方向。

國內也多次召開了相關的研討會，中國空氣動力研究與發展中心(CARDC)與中國航空工業空氣動力研究院聯合組織了DLR-F4翼身組合體、NLR7301兩段翼型數值模擬研討會和CT-1標模的大迎角數值模擬技術研討會[9]，2010年和2013年分別召開了第一、二屆航空CFD可信度開放式專題研討會，極大地促進了國內CFD的發展[10-13]。

非結構網格因為其比較容易處理復雜外形而得到了廣泛應用。對于復雜外形來說，生成非結構混合網格的時間比生成多塊對接結構網格少很多。因此計算準備時間大大減少，而且只需要較少的人為干預。非結構另一個吸引人的方面是可以方便地使用基于流場的網格自適應技術[14]。超過一半的AIAA阻力預測會議參與者使用了非結構網格技術[6]。

為了評估國內航空CFD軟件的技術狀態，明確CFD方法與軟件的下一步發展方向，推進國內CFD驗證與確認工作穩步發展，為大飛機研制提供技術參考，CARDC組織召開了第一屆CHN-T1標模CFD可信度研討會(AeCW-1)，從基本氣動力預測、氣動彈性影響、網格技術等方面，對各種CFD求解器進行統一的確認研究，促進并發展CFD對復雜流動現象的計算模擬能力。為了評估自行研制的基于非結構混合網格的流場解算器程序MFlow[15-19]對飛行器力和力矩的預測能力，課題組參加了這次會議。本文是對課題組完成的工作，即CHN-T1標模的非結構混合網格生成和采用MFlow程序完成的計算結果的分析和總結。

1 計算模型及網格

計算外形為CARDC研制的用于風洞試驗和CFD可信度驗證的，具有窄機身超臨界機翼特征的運輸機標模CHN-T1。該模型包含機身、機翼、平尾、垂尾、短艙、翼吊、起落架整流包等部件，詳細參數見文獻[20]。在研討會的標模計算時，不考慮短艙和翼吊，如圖1所示。該模型的機身采用單通道窄體尺度，代表目前中短航程商用運輸機的機身特征。機翼采用亞聲速高氣動效率的超臨界翼型下單翼。設計巡航狀態是Ma=0.78，CL=0.5。

圖1 CHN-T1外形Fig.1 Configuration of CHN-T1

AeCW-1有兩個必算算例：

Case1: 網格收斂性研究

Ma=0.78，CL=0.500(±0.001);

采用粗、中、細三套網格，網格量從約600萬到5 000萬。

Case2: 抖振特性研究

Ma=0.78;

α=-2°,-1°,0°,1°,2°,3°,3.5°,3.75°,4°,4.25°,4.5°;

采用中等規模網格;

Case2a：常規外形Config1;

Case2b：尾支撐裝置影響外形Config2;

Case2c：尾支撐加靜氣動彈性影響外形Config3;

Config1、Config2和Config3的外形比較如圖1所示。

所有計算狀態的雷諾數均為Re=3.3×106(基于平均氣動弦長)，參考溫度為300 K。

氣動系數的參考量見表1。

表1 氣動系數參考量Table 1 Reference quantity of aerodynamic coefficient

AeCW-1的網格生成準則如下：

(1) 網格收斂性算例:

CHNT-1 Wing-Body-Tail外形要求3套不同層次的網格。

(2) 網格生成指南:

a. 物面法向網格尺度:

◆ 粗網格y+<1.0,

◆ 中等網格y+<2/3,

◆ 細網格y+<4/9;

b. 建議物面法向前兩層網格尺度保持不變;

c. 網格收斂算例的網格量增長率為3倍;

d. 對于結構網格，在每個坐標方向網格增長率保持1.5倍;

e. 網格收斂算例的網格必須保持相同的網格分布，即保持相同的拉伸因子、相同的拓撲結構等等;

f. 附面層法向增長比<1.25;

g. 遠場取100倍的平均氣動弦長;

h. 對于中等網格：

◆ 機翼前后緣的弦向網格尺寸約為0.1%的當地弦長,

◆ 機翼翼根展向網格尺度約為0.1%的半展長,

◆ 機翼翼尖展向網格尺度約為0.1%的半展長,

◆ 機身頭部和后體約為2%的平均氣動弦長;

i. 對于粗網格和細網格，以上的值應該做相應的縮放;

j. 機翼后緣網格：

◆ 粗網格至少有8個單元,

◆中等網格至少有12個,

◆細網格至少有16個;

k.中等網格的網格量應滿足工業應用阻力預測的要求;

l. 滿足多重網格計算要求;

m. 對于結果網格，粗網格約為六百萬個點，中等網格約為一千五百萬個點，細網格約為五千萬個點。對于基于節點型解算器的非結構網格，網格量要求與結構網格相同，網格空間尺度的要求為內部節點間的距離。對于基于格心型解算器的非結構網格，網格空間尺度的要求為格心或面心的距離，網格點的數量約為上面的要求的1/3。

計算網格由作者使用網格生成軟件Pointwise生成，為三棱柱/四面體/金字塔組成的非結構混合網格，這套網格適合用于格心型流場求解器使用。圖2為表面網格和空間網格示意圖，在物面附近的附面層區域使用三棱柱網格，空間采用四面體，兩者交界的區域使用金字塔作為過渡。

圖2 CHN-T1網格Fig.2 CHN-T1 grid

網格收斂性研究需要的基礎網格有3套，基準網格為中等網格，然后在此基礎上分別進行細化和粗化，得到細網格和粗網格。圖3是粗、中、細三套網格在機翼剖面η=0.17處空間網格的比較。由圖可見，從物面到遠場，網格過渡平緩，沒有突變。從粗網格到細網格，網格整體加密。在附面層區域有足夠多的三棱柱單元，附面層第一層網格的高度變小，法向增長率變小，而三棱柱區域的高度基本保持不變。在機翼的前后緣，網格加密，達到網格生成指南要求的當地弦長的0.1%的尺寸。

根據會議提供的網格生成指南，網格量的比例因子為3，粗網格的y+<1.0。網格的一維單元數比例因子為31/3≈1.44。遠場取為100倍的平均氣動弦長，約為20 m。如圖4所示，在翼前后緣及翼根翼梢采用各向異性的三角形面網格單元，在滿足模擬精度的同時盡量減少網格量，提高計算效率。

圖3 η=0.17處的網格剖面Fig.3 Grids at η=0.17

圖4 各向異性三角形網格
Fig.4 Anisotropic triangle grid

由于網格生成軟件的技術限制，未能滿足網格指南中(2)b條的要求，即附面層前兩層法向取相同的尺寸。網格信息見表2。

2 數值方法

本文計算采用課題組自行發展的雷諾平均NS方程流場求解器MFlow。MFlow是基于非結構混合網格的亞跨超聲速流場解算器，可以處理任意形狀的網格單元，具有較強的靈活性。采用有限體積法對控制方程進行空間離散，未知變量存儲于網格單元的體心。

表2 網格信息Table 2 Information of grids

在本文計算中，對流通量采用Roe通量差分裂方法進行離散，具有很高的間斷和黏性分辨率，熵修正采用課題組自有的改進型Harten熵修正[15],修正參數取0.025。單元內使用線性重構使得空間離散具有二階精度，使用Venkatakrishnan限制器來抑制間斷附近的振蕩，限制參數取50。變量梯度求解使用節點型Gauss方法[15,18]。黏性通量采用中心差分離散[15]。使用當地時間步長的一階歐拉后差來達到定常狀態。Jacobian通量使用一階迎風格式得到，采用LU-SGS方法求解離散方程，CFL數取200。使用基于FMG(Full Multigrid)方法的4重“V”型多重網格方法加速收斂，粗網格的步數均為2000步。計算采用基于MPI的大規模并行計算技術加速收斂,CPU核數為128核。

計算假定流場為全湍流流場，湍流模型主要采用原始的SA一方程模型，同時采用QCR修正的SA模型計算了Case2狀態作為對比。

如圖5所示，在本文所有的計算中，收斂情況都比較好。密度的殘差都下降了十個量級以上，氣動力系數前六位有效數字保持不變。

圖5 密度殘差和力系數收斂歷程Fig.5 Convergence history of density residual and force coefficient

3 網格收斂性研究

網格收斂性研究的主要目標是估算力和力矩的網格收斂解，也就是估計當網格量趨近于無窮大時力和力矩的值。

表3是本文三套網格的計算結果、CARDC FL-26風洞試驗結果[21]和文獻[22]中采用網格量為104億的結構網格計算結果的比較。計算的阻力系數與試驗結果比較接近，細網格的阻力系數與試驗結果相差約10.2 counts(1 counts=0.0001)。俯仰力矩系數兩者差異較大，根據Case2的計算結果分析，可能是尾支桿的影響造成的。迎角相比，細網格比試驗結果小0.1568°。

“∞”表示采用中等網格和細網格線性插值得到的網格尺度無窮小，網格量無窮大的結果。與“∞”相比，細網格的阻力系數相差約-3.7 counts。“∞”與文獻結果相比，總阻力相差-3.038 counts，壓差阻力相差3.556 counts，摩擦阻力相差-0.517 counts，俯仰力矩相差-0.000 997 7，迎角相差0.114 37°。

圖6是總阻力、壓差阻力、摩擦阻力、俯仰力矩和迎角的網格收斂圖。橫坐標是N-2/3，N代表網格單元數。采用N-2/3是基于MFlow程序的數值方法為二階精度，對于由粗到細的一套自相似網格，表示其網格尺寸的二次方，因此直線就表示空間網格收斂為二階精度，N-2/3為0時表示網格量為無窮大。

由圖可見，本文計算結果當網格量趨于無窮時力系數和迎角的收斂都是單調的，而且收斂精度接近二階，證實流場解位于漸進解區域。壓差阻力CDP隨著網格加密減小，摩擦阻力CDV略有增大，因此總阻力減小。摩擦阻力隨著網格變化較小，表明當附面層網格加密到一定程度(y+<1.0)后可以準確地捕捉邊界層特性。

由于計算網格不同、解算器的計算參數等原因，文獻[22]的氣動特性計算結果與本文計算結果的網格收斂曲線不符，兩者有一定的差距。

圖7是網格對CHN-T1模型機翼的壓力分布的影響，包括三個典型展向站位，分別為0.17、0.50和0.85，包含了從翼根到翼尖各個不同的區域。圖8是平尾的三個展向站位，即0.28、0.50和0.95。

該狀態(CL=0.5)下，機翼上的流動基本沒有分離。隨著網格加密，壓力分布變化很小。網格密度對剖面壓力分布的影響主要體現在吸力峰值和對激波的分辨能力上，在機翼的其他位置相差很小。在激波附近，網格越密，激波越陡。平尾的上翼面壓力大于下翼面壓力，從而產生負升力，對飛機起抬頭作用。網格加密對平尾壓力分布的影響很小。

表3 case1狀態的CFD計算結果與試驗結果Table 3 Results of CFD and experiment for case1

圖6 網格收斂特性Fig.6 Grid-convergence properties

圖7 不同網格機翼剖面壓力分布的比較Fig.7 Comparison of Cp prediction with grid refinement of wing section

圖8 不同網格平尾剖面壓力分布的比較Fig.8 Comparison of Cp prediction on horizontal tail with grid refinement

圖9是翼身結合處物面附近的分離情況比較。CFX表示表面摩擦系數在流向的分量，經常被用于分析近壁區的流動分離情況。圖中藍色區域表示CFX<0的位置，即存在分離的區域。隨著網格的加密，模擬得到的翼身結合處分離氣泡尺度越大，總的來說，分離泡的尺寸很小，流向長度不超過當地弦長的2.5%。

4 抖振特性研究

圖10顯示了外形變化對升力系數、阻力系數和俯仰力矩系數的影響。

升力和阻力系數與試驗結果[21]符合得較好，Config1外形的俯仰力矩系數與試驗相比有一個較大的平移量，約為0.04左右，Config2和Config3的俯仰力矩系數與試驗比較接近。

考慮尾支撐(Config2)時，升力系數略有增大；阻力系數變化很小，在10 counts以內，其中壓阻系數增大，摩阻系數減小；俯仰力矩系數有一個低頭力矩增量，且隨迎角增大增量值減小。

考慮靜氣動彈性變形(Config3)的情況下，機翼的實際迎角減小，造成升力系數減小，且隨迎角增大，減小的量值增大，在3.5°以后，量值變化趨勢改變；阻力系數減小，其中壓阻系數減小明顯，摩阻系數變化很小；俯仰力矩系數有一個很小的抬頭力矩增量。

圖11和圖12是Config1和Config2外形的機翼和平尾壓力分布剖面的比較，迎角為2°。由圖可見，尾支撐對機翼的影響較小，對平尾的影響很明顯。由于尾支撐的存在，平尾產生的負升力明顯減小，從而飛機的總升力增大，由于平尾位于飛機質心之后，負升力的減小產生了一個低頭俯仰力矩增量。

圖13是Config2和Config3外形的機翼壓力分布剖面比較，迎角為3°。由于靜氣動彈性變形，機翼上翼面的激波位置前移，吸力峰值變小，從而造成機翼產生的升力減小，俯仰力矩產生一個小的抬頭增量。

圖14和圖15是Config1外形上下翼面不同迎角下的分離情況。在迎角3.5°時，機翼上翼面在激波后出現明顯的沿展向分布的小分離氣泡，在4°時，變為大分離，從而引起升力系數曲線發生拐折，升力線斜率變小。在迎角為-2°時，下翼面產生明顯的大分離。

圖9 翼身結合處的分離氣泡比較
Fig.9 Comparison of separation bubble on wing body junction

圖10 不同構型氣動特性比較Fig.10 Comparison of aerodynamics characteristic among different configurations

圖11 機翼壓力分布比較
Fig.11 Comparison ofCpprediction of wing section

圖12 平尾壓力分布比較
Fig.12 Comparison ofCpprediction of horizontal tail

圖13 機翼壓力分布比較
Fig.13 Comparison ofCpprediction of wing section

圖14 上翼面分離情況
Fig.14 Separation on wing upper surface

圖16是Config1和Config2在迎角3°時平尾的壓力分布云圖。由于尾支撐的存在，平尾下翼面靠近機身的低壓區范圍明顯變小，壓力值增大；而上翼面的低壓區范圍變大，壓力值減小。同時尾支桿對垂尾也有一定的影響，使垂尾的低壓區范圍變大，壓力值減小。

圖15 下翼面分離情況
Fig.15 Separation on wing lower surface

圖16 尾支撐對平尾壓力分布的影響Fig.16 Influence of horizontal with sting on Cp prediction

5 QCR修正的影響

表4是QCR修正SA模型和原始SA模型計算的Config1外形的升力系數、阻力系數和俯仰力矩系數的差量。在小迎角時，兩者計算結果很接近，當迎角變大(大于3.5°)時，兩者差異的差量明顯增大。

表4 QCR修正差量Table 4 Increment of QCR on the prediction of aerodynamic characteristic

6 結 論

本文使用自行發展的非結構混合網格流場解算器MFlow對AeCW-1提供的CHN-T1客機標模進行了數值模擬，研究了網格收斂性、氣動特性曲線、壓力分布和表面流動分離情況，分析結果顯示，MFlow解算器對CHN-T1標模具有較好的模擬能力。具體結論如下：

(1) 本文計算得到了近似線性的網格收斂特性。

(2) 網格加密對激波和分離氣泡的模擬精度有比較明顯的影響。

(3) 支撐系統對平尾壓力分布和俯仰力矩系數預測影響較大，對機翼影響較小。

(4) 考慮靜氣動彈性之后，機翼的當地迎角變小，壓力預測捕捉到了這一變化。

(5) QCR修正的SA模型與原始SA模型計算結果在小迎角時相差較小，大迎角下有明顯差別。