關(guān)于CFD高精度保真的數(shù)值模擬研究

張涵信

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川綿陽 621000）

關(guān)于CFD高精度保真的數(shù)值模擬研究

張涵信

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川綿陽 621000）

本文由四部分組成。第一部分研究了利用高階格式、網(wǎng)格、物理模型和高性能計(jì)算機(jī)求解NS方程時(shí)，計(jì)算數(shù)據(jù)的高精度保真性的要求以及高精度保真的指標(biāo)。第二部分研究了利用大規(guī)模計(jì)算開展大渦模擬和直接數(shù)值模擬存在的問題。第三部分分析了在高雷諾數(shù)下NS方程計(jì)算方法和網(wǎng)格的關(guān)系。最后提出了建立計(jì)算數(shù)值驗(yàn)證、確認(rèn)的新方法。

保真高精度；大渦模擬（LES）；網(wǎng)格和格式的關(guān)系；驗(yàn)證確認(rèn)

0 引 言

空氣動(dòng)力學(xué)中流動(dòng)物理方程建模的研究，方程求解高階算法的發(fā)展，網(wǎng)格生成技術(shù)的進(jìn)步，超大規(guī)模計(jì)算機(jī)的迅速出現(xiàn)及并行算法的研究，再加上實(shí)驗(yàn)設(shè)備的發(fā)展和測(cè)試技術(shù)的提高，使得航空航天等工程設(shè)計(jì)者盡管在不同設(shè)計(jì)階段需要的氣動(dòng)數(shù)據(jù)精度可以不同，但為了提高飛行器的性能，擴(kuò)大飛行走廊［1］，總要求氣動(dòng)工作者所提供的CFD和實(shí)驗(yàn)數(shù)據(jù)是高精度保真的。這就提出如下問題：

1）高精度保真的指標(biāo)是什么？給出的數(shù)值結(jié)果達(dá)到什么程度就是高精度保真了？能否給出一個(gè)哪怕是大致的指標(biāo)？

2）怎樣才能做到高精度，至少以下工作要做好：①計(jì)算出發(fā)的物理方程要準(zhǔn)確，層流的方程為NS方程是準(zhǔn)確的，但對(duì)于湍流還要尋求精確的描述；②氣動(dòng)算法應(yīng)該高精度；③計(jì)算用的網(wǎng)格、算法和物理模型要匹配。

3）怎樣判定結(jié)果已達(dá)到了高精度。通過實(shí)驗(yàn)驗(yàn)證是一種方法，但存在實(shí)驗(yàn)結(jié)果必須先證明是高精度的問題，能否另找一種方法實(shí)現(xiàn)這種驗(yàn)證呢？

本文討論這些問題。

1 高精度的指標(biāo)

當(dāng)飛行器外形給定，飛行器的氣動(dòng)力、熱等任一物理量的值就確定了（我們稱它為真值，雖然它是不知道的）。設(shè)其值可表達(dá)為：當(dāng)用CFD進(jìn)行計(jì)算時(shí)，我們可給出近似值Xc，假設(shè)計(jì)算值和真值相比，它們前n位有效數(shù)值準(zhǔn)確。那么計(jì)算值的絕對(duì)誤差的最大值是：

相對(duì)誤差是：

進(jìn)一步可以給出：

例如，當(dāng)n＝2，前兩位真值準(zhǔn)確

當(dāng)n＝3，前三位真值準(zhǔn)確

現(xiàn)有的二階精度的計(jì)算結(jié)果，氣動(dòng)壓力有的可達(dá)n＝2兩位真值準(zhǔn)確。氣動(dòng)熱還遠(yuǎn)達(dá)不到n＝2真值準(zhǔn)確。我請(qǐng)教過多位設(shè)計(jì)師，問他們準(zhǔn)確到多少位就滿意了，他們?yōu)殡y地說在實(shí)際工程中有三位真值準(zhǔn)確就行了。例如運(yùn)輸機(jī)，一般升力系數(shù)CL、力矩系數(shù)Cm和阻力系數(shù)CD要求三位真值準(zhǔn)確，可以要求：

這種情況啟示我們所謂高準(zhǔn)確的CFD計(jì)算結(jié)果，可首先定為達(dá)到三位有效數(shù)據(jù)真值準(zhǔn)確作為大致的指標(biāo)。

2 湍流LES、DNS求解問題

要進(jìn)行保真的高精度計(jì)算，首先要解決流動(dòng)方程的正確描述。湍流就是一個(gè)重要問題。湍流的計(jì)算有模式理論、大渦模擬（LES）和直接數(shù)值模擬（DNS）三種。

2.1 湍流的模式理論及LES

關(guān)于湍流計(jì)算對(duì)模式理論有以下共識(shí)：

1）以雷諾應(yīng)力方程為基礎(chǔ)，采用一方程模型、kε和k－ω等RANS模式理論［2－4］，可給出不復(fù)雜外形和不大分離區(qū)（這多屬于平衡湍流）的壓力分布、摩阻、力矩及熱流分布。壓力的計(jì)算精度高些，但摩阻、力矩和熱流的計(jì)算精度尚需進(jìn)一步改進(jìn)。

2）為了能計(jì)算復(fù)雜外形和大分離區(qū)的非平衡流動(dòng)，諸多文獻(xiàn)建議采用LES方法，但應(yīng)較好地計(jì)算亞格子網(wǎng)格應(yīng)力，特別近壁面區(qū)的模型。有人建議壁面附近用RANS、外部區(qū)域用LES的混合RANS－LES的方法［5］。

但是，我們認(rèn)為，對(duì)LES，除了壁面區(qū)的模型外，還應(yīng)研究以下問題：

1）LES原是建立在空間濾波基礎(chǔ)上的，其流動(dòng)與時(shí)間有關(guān)。RANS是建立在時(shí)間平均基礎(chǔ)上的，不同體系的RANS和LES能統(tǒng)一起來嗎？?jī)烧咧g的連接有通用的方法嗎？要尋求RANS－LES統(tǒng)一的支持理論，不然就限于經(jīng)驗(yàn)的范圍。

2）對(duì)于湍流，實(shí)驗(yàn)證明在高馬赫數(shù)下，可能存在小激波，此時(shí)LES的模型還適用嗎？

3）數(shù)值求解LES方程時(shí)，大規(guī)模并行計(jì)算（HPC）出現(xiàn)問題。設(shè)特征空間為3，劃分計(jì)算網(wǎng)格，此時(shí)小網(wǎng)格尺度為：

當(dāng)Re＝105時(shí)，尺度為Δ的網(wǎng)格點(diǎn)為109，即10億。如果Re＝107，需上萬億網(wǎng)格點(diǎn)。最近日本人用他們國(guó)家最快的計(jì)算機(jī)（10PF／S，我國(guó)天河計(jì)算機(jī)為50PF／S），用LES數(shù)值求解了汽車?yán)@流，網(wǎng)格數(shù)達(dá)數(shù)萬億。如Re＝108，需數(shù)萬萬億網(wǎng)格點(diǎn)。工程上的Re高達(dá)109。這就表明，沒有HPC，求解LES是困難的，還存在結(jié)果如何平均的問題。有了HPC，甚至有更高的運(yùn)算速度，要解決工程上高雷諾數(shù)的問題，2030年前恐怕也是困難的。因此應(yīng)從物理機(jī)理的角度發(fā)展簡(jiǎn)化的LES模式理論或其他方法。

2.2 DNS的計(jì)算

一般認(rèn)為，NS方程是湍流的物理方程。對(duì)于特征空間為3，當(dāng)雷諾數(shù)為Re時(shí)，為保證 能包含大渦，Δ可能小至Kolmogorov以下尺度。計(jì)算網(wǎng)格點(diǎn)數(shù)至少應(yīng)該為Re9／4，網(wǎng)格點(diǎn)間距為Δ，也可寫成：

對(duì)于上面的分析，n＝9／4，這是Chapman的估計(jì)，Choi和Moin指出n＝37／14。這里的指數(shù)n比LES的大很多，要對(duì)NS方程求解，網(wǎng)格點(diǎn)數(shù)比LES密很多。如Re＝107，對(duì)LES為4×1012，對(duì)DNS為5.6×1015。這樣小的尺度要求，在現(xiàn)有的超大規(guī)模計(jì)算機(jī)上是可以劃分的，甚至可劃分的更小。但是按照NS方程是連續(xù)流的條件，根據(jù)連續(xù)流的要求應(yīng)有λ／Δ＜Kn＝10－2，這里λ為分子碰撞自由程，對(duì)空氣它隨空間高度增加，在地面取最小值，約為6.68×10－8m，于是要求Δ＞100λ。此時(shí)按規(guī)定劃分網(wǎng)格，就落入了滑流的范圍，流動(dòng)的計(jì)算就進(jìn)入了另一領(lǐng)域。

但是，必須指出，DNS求解能提供非常豐富的湍流特征。它能幫助我們驗(yàn)證模式理論，幫助我們驗(yàn)證LES方法和RANS－LES方法，可以啟示我們發(fā)展更新湍流結(jié)論的認(rèn)識(shí)，發(fā)展更新的湍流描述方法。

2.3 從連續(xù)流到稀薄氣體流應(yīng)發(fā)展Boltzmann算法

流體力學(xué)的發(fā)展，先是無黏性流Euler方程求解，進(jìn)一步發(fā)展到NS方程求解，下一步自然也就是Boltzmann方程求解。這個(gè)問題早已經(jīng)提出。最近，周恒和張涵信在文獻(xiàn)［8］指出，在研究近空間高超聲速滑翔機(jī)動(dòng)飛行時(shí)，流場(chǎng)內(nèi)可能同時(shí)出現(xiàn)連續(xù)流、不連續(xù)流的情況，此時(shí)也需要研究Boltzmann方程，才能正確給出所需的物理量。再加上DNS巨型網(wǎng)格算法出現(xiàn)的物理問題，證明這一趨勢(shì)是正確的。但這里Boltzmann算法的發(fā)展，應(yīng)該考慮高溫化學(xué)反應(yīng)等效應(yīng)。目前采用的諸如Shekhov類型的BGK方法，僅是對(duì)碰撞積分的簡(jiǎn)化考慮，這里仍存在碰撞積分的模型問題，應(yīng)發(fā)展全面反應(yīng)碰撞積分的方法。

3 高階計(jì)算方法及網(wǎng)格、邊界、模型的協(xié)調(diào)匹配問題

在建立保真的流動(dòng)物理模型后，采用高精度格式是必要的。現(xiàn)已有有限差分和有限體積、緊致差分以及濾波等各類［9－11］高精度的差分方法。高精度自適應(yīng)的非結(jié)構(gòu)網(wǎng)格高精度算法可參見王志堅(jiān)主編的書［9］。

現(xiàn)在討論差分計(jì)算方面的情況。在美國(guó)，Harten和Osher等提出了ENO格式（1987，1988），它是非線性插值方法，能抑制激波處的波動(dòng)。在中國(guó)，張涵信、賀國(guó)宏、陳堅(jiān)強(qiáng)、李沁等提出用熵增條件抑制激波波動(dòng)的ENN方法（1994）。傅德薰、馬延文等提出抑制激波波動(dòng)的群速度法，還可用于緊致差分上（1994）。沈孟育研究廣義緊致格式時(shí)，以三階格式為例，證明獨(dú)立發(fā)展的ENO、熵增條件和群速度三種抑制激波波動(dòng)的方法所得到的表達(dá)式是一致的（1995）。鄧小剛（1996）在Lele高階線性節(jié)點(diǎn)和半節(jié)點(diǎn)緊致格式基礎(chǔ)上，采用原始變量插值法構(gòu)造了新穎的CNS高階緊致格式。因用了類似ENO的關(guān)系，激波附近無波動(dòng)。再后，加權(quán)無波動(dòng)格式 WENO（1994）提出后，上述格式又引入了加權(quán)分析，如WENN（1998）和WCNS（1998，2000）。這些研究成果，在驗(yàn)證典型激波或其他典型簡(jiǎn)單繞流上都能很好地給出計(jì)算結(jié)果。但對(duì)于航天、航空飛行器的計(jì)算，只有個(gè)別的結(jié)果，無法給出摩阻和熱流的滿意的計(jì)算結(jié)果。特別對(duì)于高馬赫數(shù)流動(dòng)，再加上湍流方程尚有問題，要建立通用的高階格式的平臺(tái)需再加努力。

但這節(jié)我們要強(qiáng)調(diào)重視網(wǎng)格、模型和計(jì)算格式協(xié)調(diào)匹配問題。因涉及Re及邊界的影響，需要求解物理量的導(dǎo)數(shù)，這會(huì)損壞計(jì)算精度。更主要的是涉及計(jì)算格式與網(wǎng)格和模型的匹配。設(shè)三個(gè)無量綱坐標(biāo)方向計(jì)算網(wǎng)格的尺度為為計(jì)算格式的名義精度（或稱無粘性部分計(jì)算的精度），只有當(dāng)m max（α，β，γ）時(shí)，真正的粘性流的計(jì)算精度才是高的。上面關(guān)于LES的計(jì)算，無量綱的網(wǎng)格尺度為：

4 如何驗(yàn)證確認(rèn)計(jì)算精度

如果流動(dòng)的物理方程已知，計(jì)算網(wǎng)格和高階計(jì)算格式給定，那么CFD的計(jì)算結(jié)果是否達(dá)到了所要求的精度？如何驗(yàn)證、確認(rèn)？流動(dòng)方程的理論解一般給不出，用這條道路無法判定。另一方法，例如，美國(guó)組織召開的系列DPW（Drag Predictions Workshop）研討會(huì)［12］，利用和計(jì)算相同的外形、邊界條件等做實(shí)驗(yàn)，用實(shí)驗(yàn)結(jié)果和計(jì)算值比較，但這存在實(shí)驗(yàn)值是否準(zhǔn)確，因?yàn)闇y(cè)量技術(shù)不能保證非常準(zhǔn)確。我認(rèn)為可以嘗試?yán)么髷?shù)據(jù)系統(tǒng)的概率統(tǒng)計(jì)理論，即對(duì)同一問題，用不同網(wǎng)格、不同精度格式做計(jì)算，然后對(duì)關(guān)心的物理量利用大數(shù)定律求其數(shù)學(xué)期望和方差，根據(jù)其期望和方差進(jìn)一步給出修正的數(shù)學(xué)期望的估值。這一方法我們已經(jīng)做了初步應(yīng)用［13］，得到了好的結(jié)果。文獻(xiàn)［14］給出了第二次DPM會(huì)議的結(jié)果，它是眾多研究機(jī)構(gòu)利用各自的CFD軟件、不同的算法和網(wǎng)格，圍繞DLR－F6外形（不帶發(fā)動(dòng)機(jī)短艙）進(jìn)行了大量計(jì)算。以M∞＝0.75，Rec＝3×106，CL＝0.5這一狀態(tài)為例，我們對(duì)第二屆DPW的全部計(jì)算數(shù)據(jù)進(jìn)行了分析。按我們前述的方法，三位準(zhǔn)確的真值CD＝0.029406，而精細(xì)的實(shí)驗(yàn)作比較的結(jié)果為CD＝0.0295。

5 結(jié) 論

綜上所述，有以下結(jié)論：

（1）什么是保真的高精度數(shù)據(jù)？這里提出三位有效數(shù)據(jù)真值準(zhǔn)確是否可暫作一個(gè)標(biāo)準(zhǔn)。

（2）LES和RANS的結(jié)合是一個(gè)現(xiàn)今認(rèn)為最好的模式理論，但作湍流的通用模型還有進(jìn)一步發(fā)展的空間。在實(shí)際高Re情況下，即使用HPC計(jì)算，未來幾十年內(nèi)也難以用LES解決實(shí)際工程問題。NS方程可用作湍流的基本方程，可用HPC計(jì)算，但在Re很大時(shí)，要求的網(wǎng)格出現(xiàn)稀薄氣體流動(dòng)，這又出現(xiàn)了新問題。

（3）發(fā)展高階計(jì)算方法是必要的，但必須重視網(wǎng)格問題。網(wǎng)格、邊界處理和算法要協(xié)調(diào)匹配。這里僅討論了算法精度方法的問題。因流動(dòng)方程中存在Re數(shù)，用高階格式計(jì)算時(shí)，對(duì)無粘性部分為m 階精度，但粘性部分格式精度會(huì)變低。本文提出m－max（α，β，γ） 1作為粘性部分的精度，式中

（4）除用實(shí)驗(yàn)結(jié)果驗(yàn)證確認(rèn)CFD計(jì)算值外，本文提出對(duì)同一問題，用不同格式、不同網(wǎng)格做大量計(jì)算，然后利用大數(shù)定律理論求其數(shù)學(xué)期望的方法作為驗(yàn)證、確認(rèn)方法。

致謝：本文曾多次與張樹海、張來平、李沁和葉友達(dá)進(jìn)行了討論，他們提出了很好的建議，作為八十歲的我在此表示感謝。

［1］ Tinoco E N，Bogue D R，Kao T J，et al.Progress toward CFD for full flight envelope［J］.The Aeronautical Journal，2005，109：451－460.

［2］ Spalart P R，Allmaras S R.A one－equation turbulence model for aerodynamic flows［J］.La Recherche Aerospatiale，1994，1：5－21.

［3］Menter F.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32：1598－1605.

［4］ Wilcox D C.Turbulence modeling for CFD，DCW Industries［M］.3rd edition，November 2006.

［5］Direct and large－eddy simulation 9（DLES 9）［EB／OL］.ERCOFTAC workshop，April 3－9，2013，Dresden，Germany，Http：／／www.dles9.org.

［6］ Chapman D R.Computational aerodynamics development and outlook［J］.AIAA J.，1979，17：1293.

［7］ Choi H，Moin P.Grid－point requirements for large eddy simulation：chapman’s estimates revisited［J］.Phys.Fluids，2012，24：011702.

［8］ 周恒，張涵信.空氣動(dòng)力學(xué)的新問題［J］.中國(guó)科學(xué)，物理學(xué)，力學(xué)，天文學(xué).2015，45：104789.

［9］ Wang Z J.Adaptive high order methods in computational fluid dynamics［M］.Vol.2，Advances in Computational Fluid Dynamics，World Scientific.

［10］Zhang H X，et al.Some recent progress of high－order methods on structured and unstructured grids in CARDC［C］／／Invited Lecture，ICCFD8，July 14－18，2014，Chengdu，China.

［11］Deng X G.Developing high－order linear and nonlinear schemes satisfying geometric conservation law［C］／／Invited Lecture，ICCFD8，July 14－18，2014，Chengdu，China.

［12］AIAA applied aerodynamics technical committee［C］／／Drag Prediction Workshop.19th APA Conference，Anaheim，CA 2001.

［13］Kelly R.Laflin，Olaf Bradevsen，Mark Rakowitz，John C Vassberg，Edward N.Tinoco，Richard A.Wahls，Joseph H.Morrison，Jean Luc Godard.Summary of data from the second AIAA CFD drag prediction workshop（Ivited）［R］.AIAA 2004－0555，42nd AIAA Aerospace Sviences Meeting and Exhibit 5－8，January 2004，Reno，Nevada.

［14］張涵信，查俊.關(guān)于CFD驗(yàn)證確認(rèn)中的不確定度和真值估算［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28（1）：39－45.

Investigations on fidelity of high order accurate numerical simulation for computational fluid dynamics

Zhang Hanxin

（China Aerodynamics Research and Development Center，Mianyang Sichuan 621000，China）

This paper consists of four parts.First，the requirement of fidelity and its assessment indicator for numerical results with high order accurate methods，meshes，physical models and high performance computer are investigated.Second，the problems in large eddy simulation and direct numerical simulation with high performance computer are discussed.Third，the relation between numerical scheme and mesh is analyzed for solving Navier－Stokes equations with given high Reynolds number.Finally，a new method for verification and validation of CFD simulation is proposed.

fidelity of high accuracy；LES；relation between scheme and mesh；verification and validation

V211.3

Adoi：10.7638／kqdlxxb－2015.0211

0258－1825（2016）01－0001－04

2015－12－10；

2016－01－06

張涵信（1936－），男，中科院院士，長(zhǎng)期從事計(jì)算流體力學(xué)研究.

張涵信.關(guān)于CFD高精度保真的數(shù)值模擬研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2016，34（1）：1－4.

10.7638／kqdlxxb－2015.0211 Zhang H X.Investigations on fidelity of high order accurate numerical simulation for computational fluid dynamics［J］.Acta Aerodynamica Sinica，2016，34（1）：1－4.