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

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

0 引 言

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

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

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

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

1 計算模型及網(wǎng)格

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

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

AeCW-1有兩個必算算例：

Case1: 網(wǎng)格收斂性研究

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

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

Case2: 抖振特性研究

Ma=0.78;

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

采用中等規(guī)模網(wǎng)格;

Case2a：常規(guī)外形Config1;

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

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

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

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

氣動系數(shù)的參考量見表1。

表1 氣動系數(shù)參考量Table 1 Reference quantity of aerodynamic coefficient

AeCW-1的網(wǎng)格生成準則如下：

(1) 網(wǎng)格收斂性算例:

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

(2) 網(wǎng)格生成指南:

a. 物面法向網(wǎng)格尺度:

◆ 粗網(wǎng)格y+<1.0,

◆ 中等網(wǎng)格y+<2/3,

◆ 細網(wǎng)格y+<4/9;

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

c. 網(wǎng)格收斂算例的網(wǎng)格量增長率為3倍;

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

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

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

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

h. 對于中等網(wǎng)格：

◆ 機翼前后緣的弦向網(wǎng)格尺寸約為0.1%的當?shù)叵议L,

◆ 機翼翼根展向網(wǎng)格尺度約為0.1%的半展長,

◆ 機翼翼尖展向網(wǎng)格尺度約為0.1%的半展長,

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

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

j. 機翼后緣網(wǎng)格：

◆ 粗網(wǎng)格至少有8個單元,

◆中等網(wǎng)格至少有12個,

◆細網(wǎng)格至少有16個;

k.中等網(wǎng)格的網(wǎng)格量應滿足工業(yè)應用阻力預測的要求;

l. 滿足多重網(wǎng)格計算要求;

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

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

圖2 CHN-T1網(wǎng)格Fig.2 CHN-T1 grid

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

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

圖3 η=0.17處的網(wǎng)格剖面Fig.3 Grids at η=0.17

圖4 各向異性三角形網(wǎng)格
Fig.4 Anisotropic triangle grid

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

2 數(shù)值方法

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

表2 網(wǎng)格信息Table 2 Information of grids

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

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

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

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

3 網(wǎng)格收斂性研究

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

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

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

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

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

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

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

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

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

圖6 網(wǎng)格收斂特性Fig.6 Grid-convergence properties

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

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

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

4 抖振特性研究

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

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

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

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

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

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

圖14和圖15是Config1外形上下翼面不同迎角下的分離情況。在迎角3.5°時，機翼上翼面在激波后出現(xiàn)明顯的沿展向分布的小分離氣泡，在4°時，變?yōu)榇蠓蛛x，從而引起升力系數(shù)曲線發(fā)生拐折，升力線斜率變小。在迎角為-2°時，下翼面產(chǎn)生明顯的大分離。

圖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°時平尾的壓力分布云圖。由于尾支撐的存在，平尾下翼面靠近機身的低壓區(qū)范圍明顯變小，壓力值增大；而上翼面的低壓區(qū)范圍變大，壓力值減小。同時尾支桿對垂尾也有一定的影響，使垂尾的低壓區(qū)范圍變大，壓力值減小。

圖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外形的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)的差量。在小迎角時，兩者計算結果很接近，當迎角變大(大于3.5°)時，兩者差異的差量明顯增大。

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

6 結 論

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

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

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

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

(4) 考慮靜氣動彈性之后，機翼的當?shù)赜亲冃。瑝毫︻A測捕捉到了這一變化。

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