左玲玉,司海青,李 耀,仇靜軒,李 根,吳曉軍,趙 煒

(1. 南京航空航天大學(xué)通用航空與飛行學(xué)院,江蘇 南京 211100;2. 中國空氣動(dòng)力研究與發(fā)展中心,四川 綿陽 621000)

近年來,隨著流體力學(xué)理論、數(shù)值方法以及高性能計(jì)算機(jī)技術(shù)的飛速發(fā)展,氣動(dòng)數(shù)值模擬技術(shù)(CFD)以其花費(fèi)低、耗時(shí)短、損耗小等優(yōu)點(diǎn)在航空航天領(lǐng)域發(fā)揮著越來越大的作用。然而,CFD技術(shù)要獲得更廣泛的應(yīng)用及效益,對(duì)其可信度的要求也越來越高,所以,對(duì)CFD計(jì)算結(jié)果開展可信度研究十分必要,而可信度研究的基本內(nèi)容和方法就是驗(yàn)證與確認(rèn)(V&V)。

驗(yàn)證主要是對(duì)數(shù)值模擬過程中的誤差進(jìn)行分析,以驗(yàn)證仿真結(jié)果的準(zhǔn)確性;確認(rèn)主要是將模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較,以評(píng)估數(shù)值模型與實(shí)際模型之間的相似性。國內(nèi)外對(duì)CFD的驗(yàn)證與確認(rèn)工作一直高度重視,迄今為止,大多數(shù)驗(yàn)證研究都是基于Richardson外推法獲得插值解從而得到數(shù)值計(jì)算的誤差估計(jì);Reed等通過對(duì)高超聲速飛行器二維和三維幾何模型的機(jī)理識(shí)別以及CFD網(wǎng)格收斂性研究,最大程度地降低高超聲速飛行器設(shè)計(jì)中各方面產(chǎn)生的不確定性;Lance等對(duì)混合對(duì)流進(jìn)行了實(shí)驗(yàn)研究,比較了實(shí)驗(yàn)數(shù)據(jù)和模型輸出,為CFD模型提供了完整的驗(yàn)證數(shù)據(jù);Islam等對(duì)四種不同的船型進(jìn)行驗(yàn)證與確認(rèn)研究,使用基于安全系數(shù)和修正系數(shù)的方法估計(jì)算例的不確定性,雖然使用了相似的網(wǎng)格拓?fù)?但所有四個(gè)仿真計(jì)算都得到了不同的結(jié)果;于沛從算例的角度對(duì)驗(yàn)證確認(rèn)的層次結(jié)構(gòu)、算例的分類、選擇等方面進(jìn)行了分析說明;陳堅(jiān)強(qiáng)等利用Richardson插值法對(duì)典型高超聲速流動(dòng)問題進(jìn)行了網(wǎng)格收斂性研究,然后,利用不確定度分析方法開展了該問題的確認(rèn);康順等總結(jié)了CFD模擬中的誤差及估計(jì)準(zhǔn)則,并以二維翼型為例介紹了如何進(jìn)行誤差估計(jì)和結(jié)果修正。通過總結(jié)分析各位專家學(xué)者們的研究發(fā)現(xiàn),在驗(yàn)證與確認(rèn)過程中針對(duì)數(shù)值模擬誤差的研究占據(jù)著相當(dāng)重要的位置,并且網(wǎng)格收斂性研究是必不可少的一部分。

1 CFD數(shù)值模擬誤差

在采用CFD軟件進(jìn)行數(shù)值模擬時(shí),模擬結(jié)果與真實(shí)結(jié)果之間存在一定誤差,主要包括以下幾個(gè)方面。1)物理建模誤差:由于模型簡化、控制方程、邊界條件等引入的誤差。2)離散誤差:主要源于對(duì)數(shù)學(xué)模型和定義域的離散,包括離散格式的截?cái)嗾`差、邊界條件設(shè)置不合理、網(wǎng)格質(zhì)量較差等情況導(dǎo)致的離散誤差增大。3)迭代誤差:CFD軟件數(shù)值模擬過程中采用松弛迭代法求解各類代數(shù)方程時(shí)隨時(shí)可能產(chǎn)生的誤差,包括不完全迭代或迭代的步數(shù)太少等。4)幾何外形:包括數(shù)模與真實(shí)模型外形之間的差異和生成網(wǎng)格時(shí)幾何修改或近似引入的誤差。5)來流條件變化:CFD仿真模擬中由于來流參數(shù)變化所導(dǎo)致的不確定性會(huì)使計(jì)算結(jié)果產(chǎn)生一定的偏差。

在各種誤差中,空間離散誤差是其中最重要且最難估計(jì)的部分,在誤差源中占較大比重,所以,離散誤差估計(jì)是我們的研究重點(diǎn),而網(wǎng)格收斂性研究是目前空間離散誤差估計(jì)最為普遍的方法。

2 Richardson外推法

本文采用Richardson外推法對(duì)離散誤差進(jìn)行分析,該方法利用離散步長將離散解與精確解的關(guān)系用泰勒展開的級(jí)數(shù)形式表示:

(1)

(2)

其中,收斂精度由下式確定:

(3)

同時(shí)可以得到網(wǎng)格上的誤差估計(jì):

(4)

而網(wǎng)格細(xì)化比由下式定義:

(5)

式中:是網(wǎng)格的網(wǎng)格單元數(shù),為空間維數(shù)。由于求解過程中會(huì)忽略高階項(xiàng),所以,誤差往往偏大,對(duì)式(2)添加一個(gè)修正因子來提高對(duì)誤差估計(jì)的精度,具體為

(6)

其中，是對(duì)高階項(xiàng)的估計(jì)。所以,修正后可以得到精確解:

(7)

Roache將網(wǎng)格數(shù)值解的網(wǎng)格收斂指標(biāo)()定義為

(8)

式中,為安全系數(shù),對(duì)于三套網(wǎng)格,取=125,最后定義收斂率為

(9)

當(dāng)0<<1時(shí),屬于單調(diào)收斂,當(dāng)<0時(shí)為振蕩收斂,而當(dāng)>1時(shí)為發(fā)散,Richardson插值法進(jìn)行離散誤差計(jì)算時(shí)只對(duì)單調(diào)收斂問題有效。

3 CFD離散誤差分析

3.1 NACA2412翼型CFD仿真計(jì)算

本文對(duì)NACA2412三維翼型進(jìn)行CFD數(shù)值仿真,其0°、4°、8°、12°攻角下升力系數(shù)和阻力系數(shù)的具體實(shí)驗(yàn)數(shù)據(jù)如表1所示。

表1 升力系數(shù)和阻力系數(shù)不同攻角下的實(shí)驗(yàn)值

對(duì)模型進(jìn)行了6種不同數(shù)量的網(wǎng)格劃分,其中,前三種采用結(jié)構(gòu)網(wǎng)格,后三種采用非結(jié)構(gòu)網(wǎng)格,圖1展示了采用O型網(wǎng)格劃分情況,圖2為非結(jié)構(gòu)網(wǎng)格的劃分情況,兩種類型網(wǎng)格均在機(jī)翼前緣和后緣進(jìn)行網(wǎng)格加密以便更好地捕捉流動(dòng)特性。

圖1 NACA2412結(jié)構(gòu)網(wǎng)格劃分細(xì)節(jié)圖

圖2 NACA2412非結(jié)構(gòu)網(wǎng)格劃分細(xì)節(jié)圖

本文采用Fluent軟件進(jìn)行定常計(jì)算,湍流方程采用Spalart-Allmaras模型,馬赫數(shù)=0.23,動(dòng)量方程及湍流耗散率方程均采用二階迎風(fēng)離散格式,計(jì)算殘差設(shè)置為10,經(jīng)過約1200步迭代計(jì)算后各個(gè)計(jì)算模型均達(dá)到收斂條件,表2和表3對(duì)應(yīng)于NACA2412翼型不同網(wǎng)格數(shù)量在0°、4°、8°、12°攻角時(shí)CFD計(jì)算得到的升力系數(shù)和阻力系數(shù)。將不同攻角條件下的CFD計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,發(fā)現(xiàn)隨著網(wǎng)格尺度的變化計(jì)算結(jié)果單調(diào)收斂且處于漸進(jìn)收斂域內(nèi),可以進(jìn)行網(wǎng)格收斂性研究。

表2 不同攻角下升力系數(shù)的CFD計(jì)算值

表3 不同攻角下阻力系數(shù)的CFD計(jì)算值

3.2 Richardson 外推法結(jié)果分析

為了展示網(wǎng)格類型以及網(wǎng)格粗密程度對(duì)誤差估計(jì)精度的影響,將6套網(wǎng)格分為兩組進(jìn)行研究,其中,網(wǎng)格U1、U2和U3為網(wǎng)格研究組1,網(wǎng)格U4、U5和U6為網(wǎng)格研究組2,表4為利用公式(9)計(jì)算得到的兩組網(wǎng)格不同攻角條件下的網(wǎng)格收斂判斷。針對(duì)升力系數(shù)的計(jì)算,在0°、4°及8°時(shí),網(wǎng)格研究組1(U1～U3)和網(wǎng)格研究組2(U4～U6)均為單調(diào)收斂,可以采用Richardson 插值法來估計(jì)其升力系數(shù)的計(jì)算誤差,但網(wǎng)格研究組2(U4～U6)在12°攻角時(shí)發(fā)散,不能進(jìn)行誤差估計(jì)。對(duì)于各個(gè)攻角的阻力系數(shù)的收斂性計(jì)算,兩個(gè)網(wǎng)格研究組合在小于最大升力攻角12°時(shí)均為單調(diào)收斂,可進(jìn)行誤差估計(jì),在12°攻角時(shí)都屬于振蕩收斂,無法進(jìn)行阻力系數(shù)的誤差估計(jì)。

表4 升阻力系數(shù)的收斂性判斷

針對(duì)0°、4°、8°攻角條件下升阻力系數(shù),計(jì)算得到網(wǎng)格研究組1、2對(duì)應(yīng)的離散誤差估計(jì)以及網(wǎng)格收斂指標(biāo)值(下標(biāo)1代表研究組1升力系數(shù)對(duì)應(yīng)的計(jì)算結(jié)果,下標(biāo)1代表研究組1的阻力系數(shù)對(duì)應(yīng)的計(jì)算結(jié)果,依此類推),具體見表5。

表5 不同攻角下網(wǎng)格組合升阻力系數(shù)的分析結(jié)果

結(jié)果表明,兩個(gè)網(wǎng)格研究組的網(wǎng)格誤差和網(wǎng)格收斂指標(biāo)計(jì)算結(jié)果接近,但結(jié)構(gòu)網(wǎng)格的網(wǎng)格誤差以及網(wǎng)格收斂指標(biāo)整體上均小于非結(jié)構(gòu)網(wǎng)格,同時(shí)網(wǎng)格類型的改變對(duì)升力系數(shù)的影響更大。

對(duì)0°、4°和8°攻角下升阻力系數(shù)計(jì)算值求解得到收斂精度如表6所示。

表6 不同攻角下網(wǎng)格組合升阻力系數(shù)的收斂精度

根據(jù)計(jì)算結(jié)果發(fā)現(xiàn)，同一攻角條件下升阻力系數(shù)收斂精度相差不大,說明CFD計(jì)算結(jié)果較為準(zhǔn)確,同時(shí)網(wǎng)格研究組1的收斂精度大于網(wǎng)格研究組2的收斂精度,可見結(jié)構(gòu)網(wǎng)格的離散解能以更快的速率接近精確解。所以,求解修正系數(shù)時(shí),對(duì)于網(wǎng)格研究組1根據(jù)經(jīng)驗(yàn)取2,而對(duì)于網(wǎng)格研究組2,假定結(jié)構(gòu)網(wǎng)格組合達(dá)到對(duì)應(yīng)的網(wǎng)格無關(guān)解,所以其各個(gè)攻角下的取其對(duì)應(yīng)的網(wǎng)格研究組1計(jì)算得到的收斂精度,由不同攻角下各網(wǎng)格組合升阻力系數(shù)收斂精度和對(duì)應(yīng)的修正系數(shù)進(jìn)一步求解得到修正后的精確解。

圖3和圖4給出了0°、4°、8°攻角下不同網(wǎng)格研究組合隨網(wǎng)格數(shù)量變化的CFD計(jì)算值、修正后的精確解及其實(shí)驗(yàn)數(shù)據(jù)的對(duì)比圖。其中，黑色實(shí)線為實(shí)驗(yàn)值,藍(lán)色點(diǎn)劃線為6套網(wǎng)格的CFD計(jì)算值,紅色雙點(diǎn)劃線和綠色劃線對(duì)應(yīng)于網(wǎng)格研究組修正后的精確解F1、F2,這兩條點(diǎn)劃線所覆蓋的網(wǎng)格長度對(duì)應(yīng)于各自進(jìn)行誤差估計(jì)時(shí)所采用的網(wǎng)格研究組。

圖3 不同攻角下升力系數(shù)對(duì)比圖

圖4 不同攻角下阻力系數(shù)對(duì)比

根據(jù)圖3中的結(jié)果可見,在0°、4°、8°這三種攻角條件下網(wǎng)格研究組1計(jì)算的升力系數(shù)的精確解與實(shí)驗(yàn)值的誤差均在2%以內(nèi),說明該種類型網(wǎng)格進(jìn)行計(jì)算時(shí)隨著網(wǎng)格數(shù)量的增加對(duì)機(jī)翼尾跡以及邊界層處的流動(dòng)情況模擬效果較為精確,結(jié)果趨于網(wǎng)格無關(guān)值,網(wǎng)格尺度對(duì)計(jì)算結(jié)果的影響逐漸消失。而網(wǎng)格研究組2由于CFD計(jì)算值隨網(wǎng)格單元數(shù)的變化有較大差異,所以計(jì)算所得的精確解與實(shí)驗(yàn)值的誤差遠(yuǎn)大于網(wǎng)格研究組1。

由圖4可知,網(wǎng)格研究組1對(duì)于阻力系數(shù)的計(jì)算值也隨著網(wǎng)格數(shù)量的增加逐漸達(dá)到網(wǎng)格無關(guān)解,但是本文主要考慮修正網(wǎng)格尺度變化所引起的誤差,忽略了可能涉及的不完全迭代誤差、時(shí)間步長誤差等情況,所以所得的精確解與實(shí)驗(yàn)值仍有差異。對(duì)于網(wǎng)格研究組2,在4°攻角時(shí)CFD計(jì)算所得的阻力系數(shù)隨著網(wǎng)格數(shù)量的增加更加趨近于網(wǎng)格無關(guān)解,而0°和8°攻角下得到的精確解與實(shí)驗(yàn)值的誤差較大,說明對(duì)于網(wǎng)格研究組2,目前的網(wǎng)格單元數(shù)仍未達(dá)到最優(yōu)效果,需要對(duì)網(wǎng)格尺寸進(jìn)行調(diào)整,通過更加精細(xì)的網(wǎng)格來獲得更加精確的計(jì)算結(jié)果。除此之外,可能還需要對(duì)模型本身或計(jì)算條件、邊界條件等進(jìn)行進(jìn)一步的改進(jìn)。

4 結(jié)束語

本文為評(píng)估網(wǎng)格類型及網(wǎng)格數(shù)量對(duì)NACA2412機(jī)翼CFD離散誤差的影響,對(duì)其三維模型分別進(jìn)行結(jié)構(gòu)、非結(jié)構(gòu)兩種類型共6套網(wǎng)格的劃分,并進(jìn)行了三維定常流場(chǎng)數(shù)值模擬,結(jié)合Richardson外推法和GCI網(wǎng)格收斂指標(biāo)法對(duì)升力系數(shù)、阻力系數(shù)進(jìn)行了離散誤差估計(jì),對(duì)比分析了CFD計(jì)算結(jié)果、修正后的精確值、實(shí)驗(yàn)數(shù)據(jù),完成了網(wǎng)格收斂性研究并得到結(jié)論:結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格模型的離散誤差有較大差異,網(wǎng)格類型的改變對(duì)升力系數(shù)的影響更大,兩種方法計(jì)算所得的非結(jié)構(gòu)網(wǎng)格升力系數(shù)的誤差為結(jié)構(gòu)網(wǎng)格誤差的兩倍左右;隨著網(wǎng)格單元數(shù)的增加,各個(gè)攻角條件下結(jié)構(gòu)網(wǎng)格的升阻力系數(shù)均比非結(jié)構(gòu)網(wǎng)格更接近于網(wǎng)格無關(guān)值。