BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析

王光學，張玉倫，王運濤，＊，李 松，孟德虹，李 偉

（1.國防科技大學航天科學與工程學院，湖南長沙 410073；2.中國空氣動力研究與發(fā)展中心計算空氣動力研究所，四川綿陽 621000）

BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析

王光學1，2，張玉倫2，王運濤2，＊，李 松2，孟德虹2，李 偉2

（1.國防科技大學航天科學與工程學院，湖南長沙 410073；2.中國空氣動力研究與發(fā)展中心計算空氣動力研究所，四川綿陽 621000）

高階精度計算方法的收斂效率是目前CFD重要的研究內(nèi)容之一。應用多塊對接結(jié)構(gòu)網(wǎng)格技術(shù)，基于雷諾平均的Navier-Stokes方程（RANS）和五階空間離散精度的顯式加權(quán)緊致非線性格式（WCNS-E），針對NLR7301兩段翼型和Trap Wing梯形翼兩個低速算例，重點研究了BLU-SGS迭代方法應用于WCNS-E高階精度格式上的收斂效率問題。通過與LU-SGS迭代方法收斂效率和計算結(jié)果的比較，研究結(jié)果表明：BLU-SGS迭代方法的收斂效率明顯優(yōu)于LU-SGS迭代方法；對于收斂的流場，BLU-SGS迭代方法的計算結(jié)果與LU-SGS方法的結(jié)果基本相同。

高階精度格式；計算效率；隱式時間推進；BLU-SGS；WCNS格式；低速流動

0 引 言

隨著計算機技術(shù)和CFD方法的飛速發(fā)展，CFD已經(jīng)獲得了廣泛的應用，不僅引起了飛行器設計思想的變革，導致高性能飛行器的誕生，同時在復雜流動機理研究等方面也發(fā)揮了重要作用。然而目前廣泛應用的CFD軟件大多基于二階精度格式。對于只涉及附著流和小分離流的問題，二階精度格式能夠給出較好的模擬結(jié)果，但是對于復雜流動問題，特別是涉及到導數(shù)計算的摩阻和熱流計算，二階精度格式難以滿足科學研究和工程應用的精度要求。

高階精度方法近年來獲得了廣泛、深入的研究，發(fā)展了一系列高階精度格式，如ENO、WENO、DCS、WCNS和DG等。與二階精度格式相比較，高階精度格式具有高分辨率、低耗散等特點，能夠給出更加精細的流動結(jié)構(gòu)和更準確的計算結(jié)果，已經(jīng)被廣泛應用于流動機理、氣動聲學、大渦模擬（LES）和直接數(shù)值模擬（DNS）等科學研究領(lǐng)域。但是，針對復雜的工程問題，目前尚沒有一種實用的高階精度計算方法與數(shù)值模擬軟件。主要原因在于高階精度方法對計算網(wǎng)格的質(zhì)量要求高，計算效率低，魯棒性差等。

五階空間離散精度的WCNS［1－5］格式作為一種非線性高階精度空間離散格式，不僅具有高分辨率、低耗散、低色散和一致高階精度等特點，而且對間斷具有很好的捕捉能力。其顯式格式WCNS-E，在保持上述優(yōu)良特性的前提下，相比隱式緊致格式，進一步提高了計算效率。特別是，文獻［6］針對高階有限差分算法在滿足幾何守恒律方面的困難，發(fā)展了滿足幾何守恒律的高階精度的網(wǎng)格坐標變換導數(shù)算法，大大提高了WCNS-E格式對復雜計算網(wǎng)格的適應性和魯棒性，顯示了WCNS-E格式在復雜的工程應用方面具有潛在的優(yōu)勢。

為了將高階精度格式應用于工程實際問題，還要有可以接受的計算效率。計算效率除受空間離散格式影響外，更主要受時間推進方法影響。時間推進方法可分為顯式和隱式兩類。顯式方法受CFL條件限制，計算效率較低。與顯式方法比較，隱式方法具有更高的效率。因此，在CFD應用技術(shù)研究中得到了更多的關(guān)注，發(fā)展了許多有效的隱式方法，如Beam-Warming［7］的ADI方法，Pulliam和Chaussee［8］的DADI方法，Bardina和Lombard［9］的DDADI方法，Klopfer和Hung［10］等的D3ADI方法，Jameson和Turkel［11］的LU方法，Yoon和Jameson［12］的LUSGS方法，F(xiàn)ujii和Obayashi［13］的LU-ADI方法，Briley和Neerarambam［14］的LU-AF方法，Chen和Wang［15］的BLU＿SGS方法，Jameson和Caughey［16］的SGS-Newton迭代方法，MacCormack的GSLR［17］，Reddy等［18－19］的GSPR方法等。這些方法在二階精度的CFD計算中取得了很好的效果。但是，將這些在二階精度計算中發(fā)展的隱式時間推進方法直接應用于高階精度計算，存在一個適應性問題。比如，LU-SGS方法構(gòu)造簡單，在二階精度計算中具有良好的穩(wěn)定性和收斂性，因此得到了廣泛的應用，但是將其應用于WCNS高階精度格式上效果就不太理想。究其原因，離散方程左端項的精度低，特別當右端項采用高階精度格式離散時造成左、右兩端離散精度不匹配影響了收斂速度［10，15］。針對這個問題，文獻［20］研究了LU-SGS、點松弛、線松弛以及GMRES四種方法應用于WCNS-E高階精度格式上的適應性問題。其高超聲速數(shù)值試驗表明，與LUSGS相比，采用準確雅克比矩陣的松弛法可以顯著提高高階精度格式的收斂速度。從而證實了提高離散方程左端項的精度可以達到提高收斂速度的目的。

眾所周知，流動控制方程和CFD方法在不同的速度范圍具有不同的特性和表現(xiàn)。比如，低速問題受所謂“剛性”的影響，如果不采用特殊處理（比如，預處理技術(shù)），其計算效率往往比較低。而對于跨聲速問題，控制方程性質(zhì)處于跨界狀態(tài)，有激波與邊界層干擾，流場對很多因素都比較敏感，歷來是CFD計算的難點。低速和跨聲速問題不但非常典型，而且具有廣泛的工程背景。另一方面，BLU-SGS方法也采用了更準確的方式分裂雅克比矩陣，在二階精度計算中表現(xiàn)很優(yōu)秀。所以，本文針對低速流動問題，重點研究BLU-SGS方法應用到WCNS-E高階精度格式上的收斂效率問題，并與LU-SGS方法進行比較。

1 計算方法

1.1 時間推進方法

對非定常雷諾平均的Navier-Stokes方程（URANS）在時間方向上進行隱式離散，非定常項采用一階后向差分，其它項采用一階泰勒展開，并令δ Q＝Qn＋1－Qn，δEv＝En＋1v－Env，得：

其中Q為守恒變量，τ為時間，E、F、G為對流通量，Ev、Fv、Gv為粘性通量，A、B、C分別為E、F、G的雅克比矩陣。對式（1）左端空間離散，對流項采用一階迎風差分，擴散項采用二階中心差分，并令W＝其中p為壓力，u、v、w為x、y、z三個方向的速度分量，T為溫度，Re為雷諾數(shù)。

對于離散方程（3）可以采用不同的方法進行求解。本文采用BLU-SGS和LU-SGS兩種隱式方法進行迭代求解。采用BLU-SGS時，A±Av采用完整形式的雅克比矩陣，因此D 為 一個5×5矩陣，具體方法參見文獻［15］。采用LU-SGS時，A±＝（A±βrA），rA為A的譜半徑，Av只取對角項，因此D為對角矩陣，具體方法參見文獻［12］。

1.2 空間離散方法

對于式（1）右端項的離散，對流項采用具有五階空間離散精度的WCNS-E格式，粘性項采用基于節(jié)點／半節(jié)點交錯方法的六階中心型差分。

1.2.1 對流項的離散

對流項采用具有五階空間離散精度的WCNS-E格式離散，以E／ξ為例：

其中半節(jié)點上的通量E～j＋1／2采用Roe通量差分分裂方法計算。節(jié)點上的網(wǎng)格導數(shù)通過網(wǎng)格導數(shù)的守恒算法計算。半節(jié)點上的網(wǎng)格導數(shù)利用節(jié)點上的網(wǎng)格導數(shù)插值獲得，具體方法參見文獻［6］。

1.2.2 粘性項的離散

粘性通量采用與無粘通量相同的六階中心型格式（4）進行離散，其中粘性通量vj＋1／2包含的半節(jié)點處的原始變量通過插值得到，一階導數(shù)通過六階中心型差分格式計算。

1.2.3 湍流模型的離散

湍流模型選擇k-ωSST兩方程模型，采用與RANS方程相同的方法進行離散。即對流項采用WCNS-E五階格式離散，耗散項采用六階中心型差分格式離散；湍流模型離散方程的求解方法選擇1.1節(jié)的BLU-SGS或LU-SGS方法，流動控制方程與湍流模型方程采用松耦合方式迭代求解。

2 NLR7301兩段翼型數(shù)值模擬

2.1 NLR7301兩段翼型外形及網(wǎng)格

NLR7301兩段翼型由主翼和襟翼組成，襟翼偏角為20°，主翼和襟翼之間的間隙寬度為0.026c（c為弦長），重疊區(qū)域長度為0.053c（CFG.2）。計算網(wǎng)格如圖1所示，多塊對接網(wǎng)格，網(wǎng)格單元數(shù)為191 376，第一層壁面距離5×10－6c。

圖1 NLR7301計算網(wǎng)格Fig.1 Computational mesh of NLR7301

2.2 NLR7301兩段翼型計算結(jié)果分析

計算狀態(tài)：M∞＝0.185，α＝13.1°，Re＝2.51× 106（基于平均氣動弦長）。

2.2.1 收斂效率

圖2比較了BLU-SGS和LU-SGS兩種求解方法的殘差和升力系數(shù)和阻力系數(shù)的收斂情況。從圖2中可見，兩種求解算法的氣動力收斂結(jié)果是一致的，BLU-SGS算法的殘差收斂更佳；得到收斂的氣動力系數(shù)的迭代步數(shù)比LU-SGS節(jié)約50%，雖然單步CPU時間前者要多15%左右，但是對總CPU時間而言BLU-SGS算法仍具有明顯優(yōu)勢。

2.2.2 速度型

圖3給出了NLR7301兩段翼型四個典型站位的示意圖，主翼上表面前部、主翼上表面中部、主翼下表面后部和襟翼上表面中部分別位于層流區(qū)域、湍流充分發(fā)展區(qū)域和尾跡影響區(qū)域。圖4比較了BLUSGS和LU-SGS兩種方法在這四個站位上的速度型分布，二者給出了完全重合的結(jié)果。

圖2 NLR7301翼型的收斂性歷程曲線Fig.2 Convergence histories of NLR7301airfoil

圖3 NLR7301翼型邊界層測量站位示意圖Fig.3 Boundary layer measuring stations of NLR7301airfoil

圖4 NLR7301典型站位邊界層速度型Fig.4 Typical velocity profiles for boundary layer of NLR7301airfoil

2.2.3 表面壓力分布

圖5給出了BLU-SGS和LU-SGS兩種方法得到的表面壓力分布及局部放大圖，二者是一致的。

圖5 NLR7301翼型表面壓力分布Fig.5 Distribution of pressure coefficent for NLR7301airfoil

2.2.4 氣動特性

表1給出的采用兩種求解方法得到氣動力系數(shù)及相應的試驗結(jié)果［21］，其中CDf為摩擦阻力系數(shù)、CDp為壓差阻力系數(shù)。

表1 NLR7301的氣動力系數(shù)Table 1 Aerodynamic characteristics of NLR7301

從表1看出，BLU-SGS和LU-SGS兩種方法的結(jié)果沒有明顯的差異，不同的隱式迭代方法對升力系數(shù)和阻力系數(shù)的影響非常小。采用轉(zhuǎn)捩模型可以進一步提高計阻力系數(shù)與試驗結(jié)果的吻合程度［22］。

3 Trap Wing梯形翼數(shù)值模擬

3.1 Trap Wing梯形翼構(gòu)型及網(wǎng)格

梯形翼高升力構(gòu)型是安裝在機身上的大弦長、中等展弦比、前緣縫翼／主機翼／后緣襟翼三段構(gòu)型，機翼無扭轉(zhuǎn)、無上反角，前緣縫翼偏角為30°，后緣襟翼偏角為25°（CFG.1），半展長為2.1604m，平均氣動弦長1.0067m。計算網(wǎng)格如圖6所示，為多塊對接網(wǎng)格，網(wǎng)格單元數(shù)為14 645 760，第一層法向網(wǎng)格壁面距離4×10－6。

圖6 梯形翼構(gòu)型計算網(wǎng)格Fig.6 Computational grid for trap wing configuration

3.2 Trap Wing梯形翼計算結(jié)果及分析

計算狀態(tài)：M∞＝0.2，α＝13°，Re＝4.3×106（基于平均氣動弦長）。

3.2.1 收斂效率

圖7比較了BLU-SGS和LU-SGS兩種方法的殘差和氣動力系數(shù)的收斂情況。與二維NLR7301的表現(xiàn)一致：兩種求解算法的氣動力收斂結(jié)果是一致的，BLU-SGS算法的氣動力力系數(shù)波動更?。粡牡綌?shù)和總CPU時間來看，BLU-SGS算法都具有明顯的優(yōu)勢。殘差及氣動力系數(shù)的波動是由流場中的局部分離引起的。

圖7 梯形翼構(gòu)型氣動力系數(shù)的收斂性歷程曲線Fig.7 Convergence histories of aerodynamic characteristics of trap wing configuration

3.2.2 速度型

圖8給出了邊界層內(nèi)三個典型站位上的速度型分布，分別位于主機翼后緣、襟翼前緣和襟翼后緣，都在靠近翼梢的上表面上。從圖8可見，BLU-SGS和LU-SGS兩種方法給出了相同的結(jié)果，在這三個站位上的速度型都是重合的。

3.2.3 壓力分布

圖9給出了BLU-SGS和LU-SGS兩種求解方法得到的梯形翼構(gòu)型典型展向站位的表面壓力分布曲線。在位于機翼中部的η＝0.41和η＝0.70位置，兩種求解方法的表面壓力分布是完全一致的；在位于翼梢附近的η＝0.95和η＝0.98位置，縫翼和主機翼上的表面壓力分布也保持了一致性，只是在襟翼上表面特別是襟翼后緣存在微小差異。

3.2.4 氣動特性

表2給出了采用兩種隱式求解方法得到的梯形翼構(gòu)型的氣動力系數(shù)及相應的試驗結(jié)果［23］。

從表2給出的氣動力系數(shù)來看，兩種隱式求解方法的結(jié)果沒有明顯的差異，不同的時間推進方法對氣動力系數(shù)的影響非常小。

圖8 高升力構(gòu)型典型站位邊界層速度型Fig.8 Typical velocity profiles for boundary layer of trap wing configuration

圖9 高升力構(gòu)型表面壓力分布Fig.9 Distribution of pressure coefficent for trap-wing

表2 梯形翼構(gòu)型的氣動力系數(shù)Table 2 Aerodynamic characteristics for trap wing congfiguration

4 結(jié) 論

通過NLR-7301兩段翼型和Trap Wing梯形翼兩個低速構(gòu)型的算例分析，可以得出以下結(jié)論：

（1）BLU-SGS的收斂特性顯著優(yōu)于LU-SGS，殘差下降更快。

（2）對收斂的流場，BLU-SGS方法與LU-SGS方法的計算結(jié)果基本相同，不同隱式迭代方法對計算結(jié)果影響很小。

本項研究的下一步工作將采用高階精度計算方法，針對跨聲速構(gòu)型開展BLU-SGS方法的收斂效率研究。

［1］ Deng X G，Maekawa H.Compact high-order accurate nonlinear schemes［J］.J.Comp.Phys.，1997，130：77-91.

［2］ Deng X G，Mao M L，Weighted compact high-order nonlinear schemes for the euler equations［R］.AIAA 97-1941.

［3］Deng X G，Zhang H X.Developing high-order weighted compact nonlinear schemes［J］.J.Comp.Phys.，2000，166：22-44.

［4］Deng X G.High-order accurate weighted compact nonlinear schemes［J］.Science in China，Ser.A，2001，31（12）：1104-1117.鄧小剛，高階精度耗散加權(quán)緊致非線性格式［J］.中國科學（A輯），2001，31（12）：1104-1117.

［5］ Deng X G，Liu X，Mao M L，et al.Advances in high-order accurate weighted compact nonlinear schemes［J］.Advances in Mechanics，2007，37（3）：417-427.鄧小剛，劉欣，毛枚良，等.高精度格式加權(quán)緊致非線性格式的研究進展［J］.力學進展，2007，37（3）：417-427.

［6］ Deng X G，Mao M L，Tu G H，et al.Geometric conservation law and applications to high-order finite difference schemes with stationary grids［J］.J.Comp.Phys.，2011，230：1100-1115.

［7］ Beam R，Warming R F.An implicit factored scheme for the compressible Navier-Stokes equations［J］.AIAA Journal，1978，16：393-402.

［8］ Pulliam T H，Chaussee D S.A diagonal form of an implicit approximate factorization algorithm［J］.J.Comp.Phys.，1981，39（2）：347-363.

［9］Bardina B S，Lombard H.Three dimension hypersonic flow simulations with the CSCM implicit upwind Navier-Stokes equations［R］.AIAA 87-1114.

［10］Klopfer G H，Hung C M，Van der Wijgaart R F，et al.A diagonalized diagonal dominant alternating direction implicit（D3ADI）scheme and subiteration correction［R］.AIAA 98-2824.

［11］Jameson A，Turkel E.Implicit schemes and L-U decompositions［J］.Mathematics of Computation，1981，37（156）：385-397.

［12］Yoon S，Jameson A.Lower-upper symmetric gauss-seidel method for the Euler and Navier-Stokes equations［J］.AIAA J.，1988，26：1025-1026.

［13］Fujii K，Obayashi S.Computation of three-dimensional viscous transonic flows using the LU-ADI factored scheme［R］.NALTR-899T，1986.

［14］Briley W R，Neerarambam S S，Whitfield D L.Implicit lowerupper／approximate-factorization schemes for incompressible flows［J］.J.Comp.Phys.，1996，128：32-42.

［15］Chen R F，Wang Z J.Fast，block lower-upper symmetric gaussseidel scheme for arbitrary grids［J］.AIAA Journal，2000，38（12）：2238-2245.

［16］Jameson A，Caughey D A.How many steps are required to solve the Euler equations of steady，compressible flow：in search of a fast solution algorithm［R］.AIAA 2001-2673.

［17］MacCormack R W.Current status of numerical solutions of the Navier-Stokes equations［R］.AIAA 85-0032.

［18］Reddy K C，Jacocks J L.A locally implicit scheme for the Euler equation［R］.AIAA 87-1144.

［19］Reddy K C，Benek J A.A locally implicit scheme for 3-D compressible viscous flows［R］.AIAA 90-1525.

［20］Zhang Y F.Investigations of convergence acceleration and complex flow numerical simulation for high-order accurate scheme（WCNS）［D］.China Aerodynamics Research and Development Center，2007.張毅鋒.高階精度格式（WCNS）加速收斂和復雜流動數(shù)值模擬的應用研究［D］.中國空氣動力研究與發(fā)展中心，2007.

［21］Van Den Berg B.Boundary Layer measurements on a twodimensional wing with flap［R］.NLR-TR-79009U.Amsterdam：NLR，1979.

［22］Wang Y T，Zhang Y L，Li S，et al.Calibration of aγ-Reθtransition model and its validation in low-speed flows with high-order numerical method［J］.Chinese Journal of Aeronautics，2015，28（3）：704-711.

［23］Johnson P L，Jones K M，Madson M D.Experimental Investigation of a simplified 3Dhigh lift configuration of civil transport aircraft［R］.AIAA 2008-410.Reston：AIAA，2008.

Numerical analysis of BLU＿SGS method in WCNS high-order scheme

Wang Guangxue1，2，Zhang Yulun1，Wang Yuntao2，＊，Li Song2，Meng Dehong2，Li Wei2

（1.College of Aerospace Science and Engineering，National University of Defense Technology，Changsha 410073，China；2.Computational Aerodynamics Institute of China Aerodynamics Research and Development Center，Mianyang 621000，China）

High-order numerical method is a major hotspot of computational fluid dynamics（CFD），the convergence efficiency is one of the main subjects in the study of high-order numerical method.In this paper，the block lower-upper symmetric Gauss-Seidel（BLU-SGS）implicit algorithm is applied to the fifth-order explicit Weighted Compact Nonlinear Scheme（WCNS-E）to solve the Reynolds-averaged Navier-Stokes（RANS）equations on multi-block structured grids.Two typical test cases，subsonic flows over the NLR7301two-element airfoil and the threeelement trapezoidal wing，are adopted to investigate the convergence efficiency of the BLU-SGS approach.By comparison with the original lower-upper symmetric Gauss-Setdel（LU-SGS）implicit approach，the numerical results indicate that the efficiency of BLU-SGS is much better than that of LU-SGS in low speed flow and different implicit iterative methods have little influence on aerodynamic character，such as the typical velocity profiles for boundary layer and the distribution of pressure coefficent for two typical test cases.

high-order scheme；computational efficiency；implicit time discrete；BLU＿SGS；WCNS scheme；low speed flow

V211.3

：Adoi：10.7638／kqdlxxb-2014.0061

2014-06-30；

2014-12-02

國家重點基礎研究發(fā)展計劃（2014CB744803）

王光學（1976-），男，重慶忠縣人，副研究員，研究方向：計算空氣動力學.E-mail：gxwang＠skla.cardc.cn

王運濤＊（1967-），男，黑龍江密山人，研究員，研究方向：計算空氣動力學.E-mail：ytwang＠skla.cardc.cn

王光學，張玉倫，王運濤，等.BLU-SGS方法在WCNS高階精度格式上的數(shù)值分析［J］.空氣動力學學報，2015，33（6）：733-739.

10.7638／kqdlxxb-2014.0061 Wang G X，Zhang Y L，Wang Y T，et al.Numerical analysis of BLU＿SGS method in WCNS high-order scheme［J］.Acta Aerodynamica Sinica，2015，33（6）：733-739.

0258-1825（2015）06-0733-07