何躍龍,李 盾,劉 帥,馬漢東
(中國(guó)航天空氣動(dòng)力技術(shù)研究院,北京 100074)
基于指數(shù)插值的浸沒邊界法在可壓縮流模擬中的應(yīng)用研究
何躍龍*,李 盾,劉 帥,馬漢東
(中國(guó)航天空氣動(dòng)力技術(shù)研究院,北京 100074)
由于浸沒邊界法在處理復(fù)雜外形和動(dòng)態(tài)問題時(shí)具有極強(qiáng)的靈活性和便利性,因而受到廣泛的關(guān)注。在當(dāng)前浸沒邊界法的研究和應(yīng)用中,主要工作集中在關(guān)于中/低雷諾數(shù)和不可壓流動(dòng)的模擬。本文引入Choi等提出的指數(shù)插值法對(duì)目標(biāo)點(diǎn)速度進(jìn)行插值,并根據(jù)多項(xiàng)式插值和Crocce-Busemann關(guān)系式等確定目標(biāo)點(diǎn)壓力和溫度,并借鑒直角網(wǎng)格法中的虛擬網(wǎng)格法計(jì)算物面內(nèi)的流動(dòng)變量,建立了適用于可壓縮流動(dòng)的浸沒邊界法。通過亞、跨聲速段翼型繞流和超聲速圓柱繞流算例對(duì)方法進(jìn)行驗(yàn)證,結(jié)果顯示在k值取值合適時(shí)模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)較為吻合。通過對(duì)相關(guān)算例網(wǎng)格收斂性分析,證明了本文發(fā)展的浸沒邊界法具有較高的精度。文中各算例的模擬結(jié)果表明,此方法能準(zhǔn)確捕捉亞、跨、超聲速繞流問題中的各類流動(dòng)結(jié)構(gòu),適用于可壓縮流動(dòng)的模擬。
浸沒邊界法;指數(shù)插值;可壓縮流動(dòng);直接力法;虛擬網(wǎng)格法
由于浸沒邊界法(Immersed Boundary,IB)[1-4]、直角網(wǎng)格法[5-6]等數(shù)值模擬技術(shù)可以應(yīng)用于非貼體網(wǎng)格,在處理復(fù)雜外形或運(yùn)動(dòng)邊界問題時(shí)可以大大簡(jiǎn)化網(wǎng)格生成步驟,因而在研究和工程應(yīng)用中受到越來越廣泛的關(guān)注。
浸沒邊界法最早由Peskin[1]提出,在其發(fā)展初期,研究者借助連續(xù)力法[1]、反饋力法[3]等獲得物面附近體積力的解析表達(dá)式,再借助分布函數(shù)將體積力分布在物面附近的網(wǎng)格單元上,通過這種間接的方式模擬物體對(duì)流動(dòng)的影響。但此類方法往往受到計(jì)算穩(wěn)定性的限制,且若考慮的物形為剛性物體,難以獲得有效的、解析的力源項(xiàng)表達(dá)式。Yusof等[4]提出的離散力法擴(kuò)展了浸沒邊界法的應(yīng)用范圍,提高了計(jì)算效率,大大地推進(jìn)了浸沒邊界法的發(fā)展,尤其是在固壁問題中的應(yīng)用。Fadlun等[7]研究表明,在計(jì)算浸沒單元的離散力源項(xiàng)時(shí),若不采用任何插值或僅采用體積分?jǐn)?shù)加權(quán)平均,物面邊界僅具有零階或一階精度;而通過插值方式獲取邊界附近單元的流動(dòng)變量,可有效提高計(jì)算的精度。隨著離散力法在研究中受到越來越廣泛的關(guān)注,各種邊界插值法也逐漸涌現(xiàn),如線性插值[7]、雙線性插值[8]、樣條函數(shù)插值[4]等。
目前,離散力法研究和應(yīng)用的熱點(diǎn)主要集中在低馬赫數(shù)、低雷諾數(shù)區(qū)域,這是由于常規(guī)的插值方法,如線性插值法,符合這類流動(dòng)的邊界層特性,因而在這類流動(dòng)的模擬中應(yīng)用較為成功[10];而將離散力法向可壓縮流動(dòng)區(qū)域推廣,需要對(duì)邊界附近插值法做相應(yīng)的改進(jìn)。Choi等[10]針對(duì)高雷諾數(shù)湍流邊界層提出了考慮切向速度二次修正的指數(shù)插值法,Keistler[11]將其發(fā)展應(yīng)用到超聲速流動(dòng)中,這也是目前較為少見的浸沒邊界法在可壓縮流動(dòng)中的應(yīng)用[11-14]之一。王強(qiáng)等[15]采用結(jié)合了Level-Set方法的直角網(wǎng)格法對(duì)超聲速流動(dòng)問題開展了數(shù)值模擬,其中在處理物面邊界時(shí),將速度沿法向和切向分解,分別進(jìn)行計(jì)算,以避免物面附近切向速度被抹平而形成數(shù)值邊界層,此方法與Choi[10]用不同的插值法處理切向和法向速度的思路一致。
浸沒邊界法大多采用生成簡(jiǎn)便的直角網(wǎng)格,因而有研究者將直角網(wǎng)格法等采用直角非貼體網(wǎng)格的計(jì)算方法也歸類于浸沒邊界法。通常所說的直角網(wǎng)格法,如虛擬網(wǎng)格法[6]、切割網(wǎng)格法[5]等,往往需要對(duì)邊界處切割網(wǎng)格進(jìn)行繁復(fù)處理,因此在復(fù)雜外形問題中應(yīng)用較為困難。但直角網(wǎng)格法中某些邊界處理方式,如基于虛擬網(wǎng)格法的曲率修正對(duì)稱技術(shù)[9](Curvature Corrected Symmetry Technique,CCST)等,對(duì)浸沒邊界法中物面內(nèi)流動(dòng)變量的確定有很好的借鑒意義。
本文以離散力法為基礎(chǔ),從Choi等[10]提出的指數(shù)插值法出發(fā),發(fā)展了適用于可壓縮流動(dòng)數(shù)值模擬的浸沒邊界法。文中采用NACA0012翼型在亞、跨聲速段繞流問題以及超聲速圓柱繞流算例對(duì)本文發(fā)展的方法進(jìn)行了驗(yàn)證,并對(duì)該方法的模擬能力和計(jì)算精度等開展了分析。
1.1 控制方程
考慮體積力后,積分形式的二維可壓縮N-S方程可表示為:
其中:
式中變量ρ、p、e、T和k分別表示流體密度、靜壓、內(nèi)能、溫度和熱傳導(dǎo)系數(shù);u、v分別為速度矢量V在直角坐標(biāo)系下的速度分量;而根據(jù)Stokes假設(shè),λ=-2μ/3;熱傳導(dǎo)項(xiàng)滿足qxi=-kT/xi;粘性系數(shù)μ由Sutherland公式給出。
1.2 離散方法
為了方便將此方法推廣至三維復(fù)雜外形繞流,特別是三維復(fù)雜外形動(dòng)邊界問題,本文采用基于格心的有限體積法對(duì)控制方程進(jìn)行離散。從后面的分析可以知道,當(dāng)?shù)玫轿锩娓浇W(wǎng)格體積力的離散表達(dá)形式后,可在全場(chǎng)所有矩形網(wǎng)格上統(tǒng)一求解流動(dòng)控制方程。本文計(jì)算中,采用線性MUSCL格式對(duì)網(wǎng)格單元內(nèi)流動(dòng)變量進(jìn)行重構(gòu),而單元面通量采用AUSM系列格式或Roe格式計(jì)算;時(shí)間離散可采用多步Rung-Kutta方法或LU-SGS方法。
1.3 直接力法
本文發(fā)展的方法基于Mohd-Yusof[4]提出的直接力方法。直接力法與Peskin[1]最初提出的連續(xù)力法不同,其主要特性是對(duì)離散后的動(dòng)量方程中施加力源項(xiàng),使得邊界附近網(wǎng)格流動(dòng)速度與物面的速度保持一致。
若單元i中半離散的控制方程可表示為:其中RHSi包含對(duì)流通量和粘性通量。根據(jù)Yusof的直接力法,體積力fi可表示為:
這里的uBi表示根據(jù)附近流動(dòng)速度和物面速度插值得到的目標(biāo)點(diǎn)流動(dòng)速度。
1.4 單元?jiǎng)澐?/p>
本文計(jì)算中網(wǎng)格體系采用直角網(wǎng)格,根據(jù)網(wǎng)格所處位置可以將計(jì)算域內(nèi)的網(wǎng)格單元分為四種類型:流場(chǎng)單元、過渡單元、與物面相交的切割單元以及物面內(nèi)緊鄰切割單元的虛擬單元,且切割單元又被物面分為流場(chǎng)內(nèi)和流場(chǎng)外兩部分,圖1中所示為不同類型單元分布示意圖。流場(chǎng)單元中直接求解不考慮體積力的流動(dòng)控制方程(ρfi=0);而求解過渡單元和切割單元前,需通過插值的方式獲得計(jì)算體積力所需的uBi。
圖1 單元類型示意圖Fig.1 Schematic illustrating classification of cells
由于切割單元被物面分割為流場(chǎng)內(nèi)和物面內(nèi)兩部分,因此插值時(shí)這兩部分將被區(qū)別對(duì)待。切割單元的流場(chǎng)內(nèi)區(qū)域以及過渡單元中的速度uBi直接從鄰近單元中插值(見1.6節(jié));而切割單元物面內(nèi)的部分及虛擬單元中的uBi,則需要借助虛擬網(wǎng)格法(見1.5節(jié))確定。
在獲得切割單元流場(chǎng)內(nèi)區(qū)域和物面內(nèi)區(qū)域兩部分的流動(dòng)變量后,通過體積分?jǐn)?shù)確定切割單元中各參數(shù)的平均值:
其中,q表示各流動(dòng)參數(shù),下標(biāo)I和O分別表示切割單元流場(chǎng)內(nèi)、外部分;為切割單元中流場(chǎng)內(nèi)部分占整個(gè)網(wǎng)格的體積百分?jǐn)?shù)。
1.5 虛擬網(wǎng)格法
在虛擬網(wǎng)格法中,將各物面內(nèi)節(jié)點(diǎn)向物面投影,找出此節(jié)點(diǎn)在流場(chǎng)中的投影點(diǎn),原節(jié)點(diǎn)上各流動(dòng)變量值根據(jù)此投影點(diǎn)確定。如圖1中A點(diǎn)位于虛擬單元的格心,B點(diǎn)處于流場(chǎng)中,且與A點(diǎn)關(guān)于物面對(duì)稱。在本文發(fā)展的方法中,先通過插值的方式獲得B點(diǎn)處流動(dòng)速度的值,再借助Dadone等[9]提出CCST方法獲得A點(diǎn)的值。在CCST方法中,考慮曲率影響的壓力邊界條件為:
這里R是有向局部曲率半徑,曲率中心在物面內(nèi)部時(shí)為正,反之為負(fù);Vn為切向速度。此外,還假設(shè)物面附近存在一個(gè)滿足熵與總焓局部對(duì)稱的渦流場(chǎng),此時(shí)可以根據(jù)投影點(diǎn)B處的流動(dòng)變量計(jì)算得到物面內(nèi)節(jié)點(diǎn)A處的值:
式中:Δn為線段AB的距離;Vn和Vt分別表示速度在物面法向和切向的分量。
1.6 插值方法
高效、高精度的插值方式是IB方法重點(diǎn)關(guān)注的問題。常見的插值方法有雙線性插值(三維情況為三線性插值)、最小二乘插值等;而本文研究中采用的是Choi等[10]提出的指數(shù)插值。
圖2 插值點(diǎn)模版[10]Fig.2 Interpolation stencil[10]for pointB
圖2中給出了插值的模版,xB為目標(biāo)點(diǎn),xS為目標(biāo)點(diǎn)在邊界上的投影點(diǎn),xI為插值依賴點(diǎn)。在Choi的指數(shù)插值法中,將待求點(diǎn)流動(dòng)速度分解為沿物面切向的速度ut,B和法向速度un,B,因而目標(biāo)點(diǎn)處的流動(dòng)速度可以用如下通用形式表示:
假設(shè)物面附近流動(dòng)的切向速度成指數(shù)分布,而法向速度成線性分布。目標(biāo)點(diǎn)的速度在物面切向的分量可表示為節(jié)點(diǎn)距物面法向距離的函數(shù):
其中,k為指數(shù)次數(shù),l=dB/dI為距離參數(shù)。中括號(hào)中第二項(xiàng)為二階精度的切向修正項(xiàng)。而目標(biāo)點(diǎn)法向速度為:
顯然,在此方法中,指數(shù)函數(shù)的次數(shù)k的取值是與問題相關(guān)的,Choi等[10]認(rèn)為,在極低雷諾數(shù)流動(dòng)中Re<1000,流動(dòng)速度剖面幾乎成線性分布,可以取k=1;而對(duì)于高雷諾數(shù)流動(dòng)Re>10000,可以取k=1/7或k=1/9,以避免大的流動(dòng)分離。
在計(jì)算虛擬網(wǎng)格/切割單元物面內(nèi)部分的流動(dòng)變量時(shí),還需要通過插值獲得流動(dòng)壓力、密度等物理量。對(duì)于壓力的插值,采用如下多項(xiàng)式關(guān)系:
其中b=dp/dn|n=0,是與物面加速度相關(guān)的量。而目標(biāo)點(diǎn)溫度與插值依賴點(diǎn)溫度關(guān)系根據(jù)Crocce-Busemann關(guān)系式[16]確定:
式中:Tw為壁溫,r是溫度恢復(fù)因子,在層流邊界層中可取而在湍流邊界層中取Pr1/3。
在此插值過程中,關(guān)鍵的插值依賴點(diǎn)xI的位置以及流動(dòng)變量是未知的。xI點(diǎn)的位置在過xB的壁面法線上選取,根據(jù)各點(diǎn)在此方向的投影距離確定;而xI處流動(dòng)變量值則通過臨近單元反距離插值[10]獲得。選取臨近單元時(shí),優(yōu)先選擇流場(chǎng)單元,在必須的情況下可選則過渡單元,但不選取切割單元或虛擬單元。
為了對(duì)本文發(fā)展方法進(jìn)行驗(yàn)證,這里采用本文方法發(fā)展的程序?qū)ΧS圓柱、翼型等外形在亞、跨、超聲速階段典型繞流算例開展了數(shù)值模擬。
2.1 翼型繞流
圖3中所示為本文中NACA0012算例計(jì)算采用的網(wǎng)格,在前緣細(xì)節(jié)圖中可清晰分辨出流場(chǎng)單元(紅色)、過渡單元(綠色)、切割單元(藍(lán)色)以及物面(黑色實(shí)線)。本文計(jì)算采用的網(wǎng)格是經(jīng)過分層加密的直角網(wǎng)格,計(jì)算中共采用了五套網(wǎng)格,以對(duì)計(jì)算結(jié)果進(jìn)行網(wǎng)格收斂性分析。下文分析給出的流場(chǎng)云圖及數(shù)據(jù)都來源于最密的網(wǎng)格。
首先對(duì)NACA 0012在亞、跨聲速兩個(gè)狀態(tài)進(jìn)行了模擬,分別為Case I:Ma=0.5,α=3.51°,Re=2× 105;Case II:Ma=0.703,α=3.22°,Re=3.75×106。
圖3 NACA 0012翼型計(jì)算網(wǎng)格及前緣細(xì)節(jié)圖Fig.3 Computational mesh and close-up for NACA0012
圖4中所示為Case I算例在插值指數(shù)k=1/15時(shí)模擬得到的壓力云圖與馬赫云圖。本狀態(tài)中流動(dòng)處于亞聲速階段,來流繞過翼前緣后,在背風(fēng)面膨脹,形成一個(gè)小的低壓區(qū)。
圖4 Case I壓力云圖與馬赫云圖(k=1/15)Fig.4 Pressure and Mach contour of Case I(k=1/15)
圖5中給出了該算例表面壓力系數(shù)分布數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比,圖中黑色圓點(diǎn)表示實(shí)驗(yàn)結(jié)果,而各曲線表示切向速度分布指數(shù)k取不同值時(shí),數(shù)值模擬得到的結(jié)果。從比較圖可以清晰地看出,當(dāng)k值較大(如k=1/7)時(shí),背風(fēng)面由于氣流膨脹形成的吸力峰值(絕對(duì)值)明顯偏??;而隨著k值的不斷減小,膨脹區(qū)的吸力峰值不斷增大,趨近于實(shí)驗(yàn)結(jié)果。從圖中不同k值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比看,當(dāng)k=1/15時(shí),采用本文方法模擬得到的結(jié)果與實(shí)驗(yàn)結(jié)果最為接近。
圖5 表面壓力系數(shù)數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)比較Fig.5 Comparison of Case I pressure distribution
圖6中給出的是當(dāng)k取1和1/7時(shí)模擬得到的Case I馬赫云圖。通過對(duì)不同k值模擬結(jié)果的對(duì)比容易發(fā)現(xiàn),當(dāng)k值較大時(shí),背風(fēng)區(qū)更容易出現(xiàn)分離區(qū);而隨著k值的減小,流動(dòng)的分離逐漸被抑制。文獻(xiàn)[10]中的分析認(rèn)為,當(dāng)k取1時(shí),指數(shù)插值退化為線性插值,采用此指數(shù)次數(shù)模擬邊界層內(nèi)速度可近似為線性分布的小雷諾數(shù)流動(dòng)更為合理。當(dāng)采用線性分布假設(shè)邊界層內(nèi)速度分布時(shí),此時(shí)流動(dòng)抗逆壓梯度的能力較低,更容易出現(xiàn)分離。
圖6k=1和1/7時(shí)Case I馬赫云圖Fig.6 Case I Mach contour atk=1andk=1/7
Choi等[10]提出采用指數(shù)插值的方式近似高雷諾數(shù)邊界層,文獻(xiàn)分析認(rèn)為當(dāng)k=1/7或1/9時(shí),其研究的不可壓算例數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果最為接近。我們?cè)谠缙诘南嚓P(guān)研究中也有相似的結(jié)論[17];但從前面的分析可以發(fā)現(xiàn),采用本文方法發(fā)展的程序在k=1/7時(shí)的模擬結(jié)果在翼型背風(fēng)區(qū)會(huì)出現(xiàn)明顯的分離區(qū)。本文方法與我們?cè)缙诘难芯恐胁捎玫姆椒ㄖ饕顒e在于邊界附近溫度分布的近似方法,在研究的早期,壁面溫度分布采用等熵關(guān)系式求??;而在本文對(duì)壁面附近溫度分布進(jìn)行假設(shè)時(shí),引入了Crocce-Busemann關(guān)系式。從最終的結(jié)果對(duì)比看,Crocce-Busemann關(guān)系式的引入(溫度分布關(guān)系的改變),對(duì)模擬結(jié)果的影響較為顯著。文獻(xiàn)[17]中采用等熵關(guān)系式時(shí),采用較大的k值(如1/5)仍能模擬得到較為準(zhǔn)確的結(jié)果,顯然采用此溫度分布近似時(shí)對(duì)流動(dòng)分離過程會(huì)有一定的抑制作用;而本文采用Crocce-Busemann關(guān)系式時(shí),在較大k值的各模擬狀態(tài)中,背風(fēng)區(qū)流動(dòng)分離較為明顯。從前文的分析可以發(fā)現(xiàn),不同溫度分布近似對(duì)流動(dòng)分離過程的影響還值得深入的研究。
Case II狀態(tài)為跨聲速狀態(tài),雖然來流馬赫數(shù)處于亞聲速,但由于來流相對(duì)翼型具有一定的迎角,在翼型的背風(fēng)區(qū)由于流動(dòng)的膨脹出現(xiàn)超聲速區(qū),且在翼型背部有明顯的激波存在。圖7給出了采用本文方法模擬得到的Case II在插值指數(shù)次數(shù)k取1/15時(shí)的壓力云圖與馬赫云圖。
從圖8中k取不同值時(shí),數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比可以發(fā)現(xiàn),對(duì)于Case II這種跨聲速算例,仍是當(dāng)k=1/15時(shí)數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)最為接近,此時(shí)模擬得到的機(jī)翼背風(fēng)面激波的位置和強(qiáng)度與實(shí)驗(yàn)數(shù)據(jù)都極為吻合。
圖7 Case II壓力云圖與馬赫云圖(k=1/15)Fig.7 Pressure and Mach contour of Case II(k=1/15)
圖8 Case II壁面壓力分布數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.8 Comparison of the Case II pressure distribution
從圖9中不同k值時(shí)模擬得到的Case II馬赫云圖可以發(fā)現(xiàn),當(dāng)k取值較小時(shí),與Case I中情況類似,在翼型背風(fēng)區(qū)后緣會(huì)出現(xiàn)明顯的流動(dòng)分離。由于分離的出現(xiàn),使得背風(fēng)區(qū)前部膨脹流動(dòng)的膨脹角減小,膨脹區(qū)下游的主流流動(dòng)速度降低,進(jìn)而使得流動(dòng)膨脹引起的激波強(qiáng)度減弱甚至消失。因而在圖8中可以觀察到,當(dāng)k取1/7時(shí),機(jī)翼背風(fēng)區(qū)的吸力峰值變化光滑,此時(shí)膨脹區(qū)沒有形成激波;而在k=1/9時(shí),吸力峰值后部開始出現(xiàn)跳躍;在k為1/11的狀態(tài),出現(xiàn)了明顯的激波,但此時(shí)激波強(qiáng)度比實(shí)驗(yàn)結(jié)果要弱,且位置靠前。
圖9k=1/7和1/11時(shí)Case II馬赫云圖Fig.9 Case II Mach contour atk=1/7andk=1/11
2.2 網(wǎng)格收斂性
這里選用NACA 0012翼型算例對(duì)本文方法的網(wǎng)格收斂性開展研究。如前所述,本文研究過程中共采用了五套計(jì)算網(wǎng)格。第一套(最粗)網(wǎng)格計(jì)算域?yàn)榈谝粚泳W(wǎng)格為(50×50),之后在物面附近進(jìn)行了七次分層加密。第二~五套網(wǎng)格由前一套網(wǎng)格進(jìn)行一分為四的四叉樹加密得到。圖3中所示即為第五套網(wǎng)格的細(xì)節(jié)圖。五套網(wǎng)格總的網(wǎng)格數(shù)及最小網(wǎng)格的尺度如表1所示。
表1 網(wǎng)格總數(shù)與網(wǎng)格尺度列表Table 1 List of total cell numbers and minimum mesh scales
在開展網(wǎng)格收斂性研究時(shí),根據(jù)第四和第五套網(wǎng)格的模擬結(jié)果,采用Richardson外推法[18]插值獲得“精確”解。外推過程中,以最密兩套網(wǎng)格的數(shù)值解為基礎(chǔ),對(duì)于任意的流動(dòng)變量φ,其外推解φe計(jì)算方法如下:
將“精確”解通過插值再映射到各套網(wǎng)格上,可求得每套網(wǎng)格下模擬結(jié)果與外推“精確”解之間的誤差。圖10中給出了兩狀態(tài)各流動(dòng)變量誤差的1范數(shù)變化曲線,其中黑色點(diǎn)線和虛線為示意斜率k=1和k=2的輔助線。從各曲線與斜率輔助線的對(duì)比可以看出,兩狀態(tài)中各變量在密網(wǎng)格情況下都能達(dá)到二階精度。
圖10 各變量誤差的1范數(shù)曲線Fig.10 Norm-1curves of the variables in the two cases
2.3 圓柱繞流
為檢驗(yàn)本文方法對(duì)超聲速流動(dòng)的模擬能力,采用來流條件為Re∞=2×105,M∞=1.3和1.7的二維圓柱繞流算例進(jìn)行驗(yàn)證。計(jì)算仍采用分層加密的非均勻直角網(wǎng)格。根據(jù)相關(guān)文獻(xiàn)[13,19]的實(shí)驗(yàn)和計(jì)算結(jié)果可知,這兩狀態(tài)中圓柱尾跡為對(duì)稱的穩(wěn)定渦結(jié)構(gòu),因此采用半模進(jìn)行計(jì)算。計(jì)算域大小為[100d× 50d](半模),第一層網(wǎng)格(800×400),在物面附近進(jìn)行七次加密,總的計(jì)算網(wǎng)格約為36萬。
圖11和圖12分別為采用本文方法模擬得到的Ma=1.3和1.7兩狀態(tài)在k=1時(shí)的空間壓力云圖與流線圖。從圖中可以觀察到清晰的脫體激波、流動(dòng)繞過圓柱時(shí)形成的膨脹區(qū)以及尾跡區(qū)尾部拖出的膨脹激波。從圖中可以看出,在超聲速條件下,來流馬赫數(shù)越高,其脫體激波的脫體距離越小,但波后壓力越高;圓柱后尾渦的長(zhǎng)度也隨著馬赫數(shù)的增大而減小。總體來說,本文方法成功地捕捉到了超聲速條件圓柱繞流的各種典型流動(dòng)特征。
圖11Ma=1.3狀態(tài)空間壓力云圖與流線圖(k=1)Fig.11 Pressure contour and streamline atMa=1.3(k=1)
圖12Ma=1.7狀態(tài)空間壓力云圖與流線圖(k=1)Fig.12 Pressure contour and streamline atMa=1.7(k=1)
圖13中所示是兩狀態(tài)在指數(shù)次數(shù)k取不同值時(shí)模擬得到的圓柱表面壓力分布于實(shí)驗(yàn)數(shù)據(jù)的對(duì)比圖。從圖中可以看出,在圓柱迎風(fēng)面,選取不同k值時(shí)模擬結(jié)果幾乎相同,且數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本吻合;但在尾跡區(qū),尤其是分離起始位置附近,不同k值的模擬結(jié)果存在一定差異。
通過對(duì)比可以發(fā)現(xiàn),在尾跡區(qū),不同k值的模擬結(jié)果之間的差異主要在于分離起始位置的子午角不同。在分離之前,不同k值模擬結(jié)果幾乎相同;如果分離時(shí)子午角增大,分離前流動(dòng)速度也增大,壁面壓力越低;但在分離后壁面的壓力分布都將恢復(fù)到幾乎恒定的水平。從對(duì)比圖中可以明顯看出,當(dāng)k的取值較小時(shí),其分離明顯靠后;而較大k值對(duì)應(yīng)的分離角明顯靠前。且從兩狀態(tài)的模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比看,當(dāng)k=1時(shí)模擬結(jié)果與實(shí)驗(yàn)結(jié)果最為接近。從前文的分析可知,k取1時(shí)速度插值相當(dāng)于線性近似,有利于分離流動(dòng)的模擬;而k取值較小時(shí),則會(huì)抑制分離的出現(xiàn)。顯然,從本算例的模擬結(jié)果中可以得出類似的結(jié)論。在這種低速大分離的尾跡區(qū)域,采用線性近似的效果更好。
圖13 不同k值時(shí)圓柱表面壓力分布與實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.13 Comparison of the pressure distribution
但從本算例的研究中又引出了新的問題,即當(dāng)前方法插值次數(shù)k的取值為全場(chǎng)統(tǒng)一,這在某些流場(chǎng)結(jié)構(gòu)相對(duì)單一的情況,如前述NACA 0012翼型繞流問題,可以得到很好的應(yīng)用。但當(dāng)面對(duì)的問題為混合了高/低速流動(dòng)時(shí),如何進(jìn)行動(dòng)態(tài)的指數(shù)插值,是需要進(jìn)一步研究的課題。
本文發(fā)展了一種基于指數(shù)插值的浸沒邊界法。此方法在壁面附近根據(jù)速度、壓力和溫度等流動(dòng)量自身特性,分別引入了不同的插值方式,實(shí)現(xiàn)了在非貼體網(wǎng)格上對(duì)流場(chǎng)進(jìn)行求解的算法。
采用亞、跨聲速段翼型繞流和超聲速圓柱繞流算例對(duì)本文方法進(jìn)行了驗(yàn)證,結(jié)果顯示在選定插值次數(shù)k的條件下,本文方法能很好的模擬亞、跨、超聲速流動(dòng)問題。顯然本文方法具有很強(qiáng)的可壓縮流動(dòng)模擬能力。從本文的驗(yàn)證過程可以發(fā)現(xiàn),當(dāng)k取值較大時(shí),插值過程接近線性近似,此時(shí)在尾流等存在流動(dòng)分離或低速流動(dòng)區(qū)域模擬效果更佳;而k取值較小時(shí),對(duì)激波等強(qiáng)壓縮性區(qū)域的模擬效果更好。
從本文的分析也可以發(fā)現(xiàn),在流動(dòng)變量插值關(guān)系式的理論分析,以及此方法在強(qiáng)壓縮性和大分離流動(dòng)混合的復(fù)雜流動(dòng)問題的應(yīng)用方面還需要開展相關(guān)研究。
[1]Peskin C S.Flow patterns around heart valves:A numerical method[J].Journal of Computational Physics,1972,10:252-271.
[2]Mittal R,Iaccarino G.Immersed boundary methods[J].Annu. Rev.Fluid Mech.,2005,37:239-261.
[3]Goldstein D,Handler R,Sirovich L.Modeling a no-slip flow boundary with an external force field[J].Journal of Computational Physics,1989,105:354-366.
[4]Mohd-Yusof J.Combined immersed boundaries/B-Splines methods for simulations of flows in complex geometries[R].CTR Annual Research Briefs,NASA Ames/Stanford University,1997.
[5]Udaykumar H S,Mittal R,Rampunggoon P,et al.A sharp interface Cartesian grid method for simulating flows with complex moving boundaries[J].J.Comput.Phys.,2001,174:345-380.
[6]Majumdar S,Iaccarino G,Durbin P A.RANS solver with adaptive structured boundary non-conforming grids[R].Annu.Res.Briefs,Cent.Turbul.Res.353-364.
[7]Fadlun E A,Verzicco R,Orlandi P,et al.Combined immersed-boundary finite-difference methods for 3Dcomplex flow simulations[J].J.Comput.Phys.,2000,161:25-60.
[8]Kim J,Kim D,Choi H.An immersed boundary finite volume method for simulations of flow in complex geometries[J].Journal of Computational Physics,2001,171:132-150.
[9]Dadone A,Grossman B.An immersed body methodology for inviscid flows on Cartesian grids[R].AIAA-2002-1059.
[10]Choi J I,Oberoi R C,Edwards J R,et al.An immersed boundary method for complex incompressible flows[J].Journal of Computational Physics,2007,224:757-784.
[11]Keistler P G.An immersed boundary method for supersonic flow[R].AIAA-2008-529.
[12]Cho J,Chopra J,Morris P J.Immersed boundary method for compressible high-reynolds number viscous flow around moving bodies[R].AIAA-2007-125.
[13]De Palma P,De Tullio M D,Pascazio G,et al.An immersed Boundary method for compressible viscous flows[J].Computers&Fluids,2006,35:639-702.
[14]Liu J M.The immersed boundary method in compressible flows and its applications[D].Nanjing university of aeronautics and astronautics,2010.(in Chinese)劉劍明.可壓縮流體計(jì)算中的浸入邊界方法及其應(yīng)用[D].南京航空航天大學(xué),2010.
[15]Wang Q,Hu X Y,Jiang Z L.A Cartesian mesh algorithm for supersonic flows around arbitrary moving bodies[J].Chinese Journal of Computational Physics,2009,26(4):517-526.(in Chinese)王強(qiáng),胡湘渝,姜宗林.一種含運(yùn)動(dòng)固壁超聲速流動(dòng)的Descartes網(wǎng)格算法[J].計(jì)算物理,2009,26(4):517-526.
[16]Chen M Z.Fundamentals of viscous fluid dynamics[M].Higher education press,2006.(in Chinese)陳懋章.粘性流體動(dòng)力學(xué)基礎(chǔ)[M].高等教育出版社,2006.
[17]He Y L,Li D,Liu S,et al.An immersed boundary method based on volume fraction[C]//Asia-Pacific International Symposium on Aerospace Technology,2014,Shanghai.
[18]Chen J Q,Zhang Y R.Verification and validation in CFD based on the Richardson extrapolation method[J].Acta Aerodynamica Sinica,2012,30(2):176-183.(in Chinese)陳堅(jiān)強(qiáng),張益榮.基于Richardson插值法的CFD驗(yàn)證和確認(rèn)方法的研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2012,30(2):176-183.
[19]Bashkin V A,Vaganov A V,Egorov I V,et al.Comparison of calculated and experimental data on supersonic flow past a circular cylinder[J].Fluid Dynamics,2002,37(3):473-482.
The application of power-law interpolation based immersed boundary method in compressible flow simulation
He Yuelong*,Li Dun,Liu Shuai,Ma Handong
(China Academy of Aerospace Aerodynamics,Beijing 100074,China)
The immersed boundary method(IBM)was rapid developed in recent years because of its flexibility for complex configurations and moving bodies.Up to date,most of the researches on IB method focused on the incompressible flow due to the limitation of the interpolation method.The current research applied the Choi’s power-law interpolation and the Crocce-Busemann relationship for the velocity and temperature approximation in the objective cells.The ghost cell method was also used for computing the variables at the cells in the body.The subsonic/transonic NACA0012airfoil compressible cases and 2Dsupersonic cylinder benchmark cases were chosen for method validation.The numerical result of these benchmark cases with proper parameter k compared well with the experimental data.The mesh convergence of the new method has also been proved.The analysis and discussion show that the method is suitable for compressible flow simulation.
immersed boundary;power-law interpolation;compressible flow;direct force;ghost cell method
V211.3
Adoi:10.7638/kqdlxxb-2014.0107
0258-1825(2016)04-0426-07
2014-09-17;
2014-12-16
國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃,2014CB744100)
何躍龍*(1985-),男,湖南湘潭人,工程師,博士,研究方向:計(jì)算流體力學(xué).E-mail:aia741@163.com
何躍龍,李盾,劉帥,等.基于指數(shù)插值的浸沒邊界法在可壓縮流模擬中的應(yīng)用研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2016,34(4):426-432.
10.7638/kqdlxxb-2014.0107 He Y L,Li D,Liu S,et al.The application of power-law interpolation based immersed boundary method in compressible flow simulation[J].Acta Aerodynamica Sinica,2016,34(4):426-432.