基于指數(shù)插值的浸沒邊界法在可壓縮流模擬中的應(yīng)用研究

何躍龍，李 盾，劉 帥，馬漢東

（中國航天空氣動力技術(shù)研究院，北京 100074）

基于指數(shù)插值的浸沒邊界法在可壓縮流模擬中的應(yīng)用研究

何躍龍＊，李 盾，劉 帥，馬漢東

（中國航天空氣動力技術(shù)研究院，北京 100074）

由于浸沒邊界法在處理復(fù)雜外形和動態(tài)問題時(shí)具有極強(qiáng)的靈活性和便利性，因而受到廣泛的關(guān)注。在當(dāng)前浸沒邊界法的研究和應(yīng)用中，主要工作集中在關(guān)于中／低雷諾數(shù)和不可壓流動的模擬。本文引入Choi等提出的指數(shù)插值法對目標(biāo)點(diǎn)速度進(jìn)行插值，并根據(jù)多項(xiàng)式插值和Crocce－Busemann關(guān)系式等確定目標(biāo)點(diǎn)壓力和溫度，并借鑒直角網(wǎng)格法中的虛擬網(wǎng)格法計(jì)算物面內(nèi)的流動變量，建立了適用于可壓縮流動的浸沒邊界法。通過亞、跨聲速段翼型繞流和超聲速圓柱繞流算例對方法進(jìn)行驗(yàn)證，結(jié)果顯示在k值取值合適時(shí)模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)較為吻合。通過對相關(guān)算例網(wǎng)格收斂性分析，證明了本文發(fā)展的浸沒邊界法具有較高的精度。文中各算例的模擬結(jié)果表明，此方法能準(zhǔn)確捕捉亞、跨、超聲速繞流問題中的各類流動結(jié)構(gòu)，適用于可壓縮流動的模擬。

浸沒邊界法；指數(shù)插值；可壓縮流動；直接力法；虛擬網(wǎng)格法

0 引 言

由于浸沒邊界法（Immersed Boundary，IB）［1－4］、直角網(wǎng)格法［5－6］等數(shù)值模擬技術(shù)可以應(yīng)用于非貼體網(wǎng)格，在處理復(fù)雜外形或運(yùn)動邊界問題時(shí)可以大大簡化網(wǎng)格生成步驟，因而在研究和工程應(yīng)用中受到越來越廣泛的關(guān)注。

浸沒邊界法最早由Peskin［1］提出，在其發(fā)展初期，研究者借助連續(xù)力法［1］、反饋力法［3］等獲得物面附近體積力的解析表達(dá)式，再借助分布函數(shù)將體積力分布在物面附近的網(wǎ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)平均，物面邊界僅具有零階或一階精度；而通過插值方式獲取邊界附近單元的流動變量，可有效提高計(jì)算的精度。隨著離散力法在研究中受到越來越廣泛的關(guān)注，各種邊界插值法也逐漸涌現(xiàn)，如線性插值［7］、雙線性插值［8］、樣條函數(shù)插值［4］等。

目前，離散力法研究和應(yīng)用的熱點(diǎn)主要集中在低馬赫數(shù)、低雷諾數(shù)區(qū)域，這是由于常規(guī)的插值方法，如線性插值法，符合這類流動的邊界層特性，因而在這類流動的模擬中應(yīng)用較為成功［10］；而將離散力法向可壓縮流動區(qū)域推廣，需要對邊界附近插值法做相應(yīng)的改進(jìn)。Choi等［10］針對高雷諾數(shù)湍流邊界層提出了考慮切向速度二次修正的指數(shù)插值法，Keistler［11］將其發(fā)展應(yīng)用到超聲速流動中，這也是目前較為少見的浸沒邊界法在可壓縮流動中的應(yīng)用［11－14］之一。王強(qiáng)等［15］采用結(jié)合了Level－Set方法的直角網(wǎng)格法對超聲速流動問題開展了數(shù)值模擬，其中在處理物面邊界時(shí)，將速度沿法向和切向分解，分別進(jìn)行計(jì)算，以避免物面附近切向速度被抹平而形成數(shù)值邊界層，此方法與Choi［10］用不同的插值法處理切向和法向速度的思路一致。

浸沒邊界法大多采用生成簡便的直角網(wǎng)格，因而有研究者將直角網(wǎng)格法等采用直角非貼體網(wǎng)格的計(jì)算方法也歸類于浸沒邊界法。通常所說的直角網(wǎng)格法，如虛擬網(wǎng)格法［6］、切割網(wǎng)格法［5］等，往往需要對邊界處切割網(wǎng)格進(jìn)行繁復(fù)處理，因此在復(fù)雜外形問題中應(yīng)用較為困難。但直角網(wǎng)格法中某些邊界處理方式，如基于虛擬網(wǎng)格法的曲率修正對稱技術(shù)［9］（Curvature Corrected Symmetry Technique，CCST）等，對浸沒邊界法中物面內(nèi)流動變量的確定有很好的借鑒意義。

本文以離散力法為基礎(chǔ)，從Choi等［10］提出的指數(shù)插值法出發(fā)，發(fā)展了適用于可壓縮流動數(shù)值模擬的浸沒邊界法。文中采用NACA0012翼型在亞、跨聲速段繞流問題以及超聲速圓柱繞流算例對本文發(fā)展的方法進(jìn)行了驗(yàn)證，并對該方法的模擬能力和計(jì)算精度等開展了分析。

1 數(shù)值方法

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ù)雜外形動邊界問題，本文采用基于格心的有限體積法對控制方程進(jìn)行離散。從后面的分析可以知道，當(dāng)?shù)玫轿锩娓浇W(wǎng)格體積力的離散表達(dá)形式后，可在全場所有矩形網(wǎng)格上統(tǒng)一求解流動控制方程。本文計(jì)算中，采用線性MUSCL格式對網(wǎng)格單元內(nèi)流動變量進(jìn)行重構(gòu)，而單元面通量采用AUSM系列格式或Roe格式計(jì)算；時(shí)間離散可采用多步Rung－Kutta方法或LU－SGS方法。

1.3 直接力法

本文發(fā)展的方法基于Mohd－Yusof［4］提出的直接力方法。直接力法與Peskin［1］最初提出的連續(xù)力法不同，其主要特性是對離散后的動量方程中施加力源項(xiàng)，使得邊界附近網(wǎng)格流動速度與物面的速度保持一致。

若單元i中半離散的控制方程可表示為：其中RHSi包含對流通量和粘性通量。根據(jù)Yusof的直接力法，體積力fi可表示為：

這里的uBi表示根據(jù)附近流動速度和物面速度插值得到的目標(biāo)點(diǎn)流動速度。

1.4 單元劃分

本文計(jì)算中網(wǎng)格體系采用直角網(wǎng)格，根據(jù)網(wǎng)格所處位置可以將計(jì)算域內(nèi)的網(wǎng)格單元分為四種類型：流場單元、過渡單元、與物面相交的切割單元以及物面內(nèi)緊鄰切割單元的虛擬單元，且切割單元又被物面分為流場內(nèi)和流場外兩部分，圖1中所示為不同類型單元分布示意圖。流場單元中直接求解不考慮體積力的流動控制方程（ρfi＝0）；而求解過渡單元和切割單元前，需通過插值的方式獲得計(jì)算體積力所需的uBi。

圖1 單元類型示意圖Fig.1 Schematic illustrating classification of cells

由于切割單元被物面分割為流場內(nèi)和物面內(nèi)兩部分，因此插值時(shí)這兩部分將被區(qū)別對待。切割單元的流場內(nèi)區(qū)域以及過渡單元中的速度uBi直接從鄰近單元中插值（見1.6節(jié)）；而切割單元物面內(nèi)的部分及虛擬單元中的uBi，則需要借助虛擬網(wǎng)格法（見1.5節(jié)）確定。

在獲得切割單元流場內(nèi)區(qū)域和物面內(nèi)區(qū)域兩部分的流動變量后，通過體積分?jǐn)?shù)確定切割單元中各參數(shù)的平均值：

其中，q表示各流動參數(shù)，下標(biāo)I和O分別表示切割單元流場內(nèi)、外部分；為切割單元中流場內(nèi)部分占整個(gè)網(wǎng)格的體積百分?jǐn)?shù)。

1.5 虛擬網(wǎng)格法

在虛擬網(wǎng)格法中，將各物面內(nèi)節(jié)點(diǎn)向物面投影，找出此節(jié)點(diǎn)在流場中的投影點(diǎn)，原節(jié)點(diǎn)上各流動變量值根據(jù)此投影點(diǎn)確定。如圖1中A點(diǎn)位于虛擬單元的格心，B點(diǎn)處于流場中，且與A點(diǎn)關(guān)于物面對稱。在本文發(fā)展的方法中，先通過插值的方式獲得B點(diǎn)處流動速度的值，再借助Dadone等［9］提出CCST方法獲得A點(diǎn)的值。在CCST方法中，考慮曲率影響的壓力邊界條件為：

這里R是有向局部曲率半徑，曲率中心在物面內(nèi)部時(shí)為正，反之為負(fù)；Vn為切向速度。此外，還假設(shè)物面附近存在一個(gè)滿足熵與總焓局部對稱的渦流場，此時(shí)可以根據(jù)投影點(diǎn)B處的流動變量計(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)流動速度分解為沿物面切向的速度ut，B和法向速度un，B，因而目標(biāo)點(diǎn)處的流動速度可以用如下通用形式表示：

假設(shè)物面附近流動的切向速度成指數(shù)分布，而法向速度成線性分布。目標(biāo)點(diǎn)的速度在物面切向的分量可表示為節(jié)點(diǎn)距物面法向距離的函數(shù)：

其中，k為指數(shù)次數(shù)，l＝dB／dI為距離參數(shù)。中括號中第二項(xiàng)為二階精度的切向修正項(xiàng)。而目標(biāo)點(diǎn)法向速度為：

顯然，在此方法中，指數(shù)函數(shù)的次數(shù)k的取值是與問題相關(guān)的，Choi等［10］認(rèn)為，在極低雷諾數(shù)流動中Re＜1000，流動速度剖面幾乎成線性分布，可以取k＝1；而對于高雷諾數(shù)流動Re＞10000，可以取k＝1／7或k＝1／9，以避免大的流動分離。

在計(jì)算虛擬網(wǎng)格／切割單元物面內(nèi)部分的流動變量時(shí)，還需要通過插值獲得流動壓力、密度等物理量。對于壓力的插值，采用如下多項(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的位置以及流動變量是未知的。xI點(diǎn)的位置在過xB的壁面法線上選取，根據(jù)各點(diǎn)在此方向的投影距離確定；而xI處流動變量值則通過臨近單元反距離插值［10］獲得。選取臨近單元時(shí)，優(yōu)先選擇流場單元，在必須的情況下可選則過渡單元，但不選取切割單元或虛擬單元。

2 算例驗(yàn)證

為了對本文發(fā)展方法進(jìn)行驗(yàn)證，這里采用本文方法發(fā)展的程序?qū)ΧS圓柱、翼型等外形在亞、跨、超聲速階段典型繞流算例開展了數(shù)值模擬。

2.1 翼型繞流

圖3中所示為本文中NACA0012算例計(jì)算采用的網(wǎng)格，在前緣細(xì)節(jié)圖中可清晰分辨出流場單元（紅色）、過渡單元（綠色）、切割單元（藍(lán)色）以及物面（黑色實(shí)線）。本文計(jì)算采用的網(wǎng)格是經(jīng)過分層加密的直角網(wǎng)格，計(jì)算中共采用了五套網(wǎng)格，以對計(jì)算結(jié)果進(jìn)行網(wǎng)格收斂性分析。下文分析給出的流場云圖及數(shù)據(jù)都來源于最密的網(wǎng)格。

首先對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)中流動處于亞聲速階段，來流繞過翼前緣后，在背風(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ù)的對比，圖中黑色圓點(diǎn)表示實(shí)驗(yàn)結(jié)果，而各曲線表示切向速度分布指數(shù)k取不同值時(shí)，數(shù)值模擬得到的結(jié)果。從比較圖可以清晰地看出，當(dāng)k值較大（如k＝1／7）時(shí)，背風(fēng)面由于氣流膨脹形成的吸力峰值（絕對值）明顯偏?。欢S著k值的不斷減小，膨脹區(qū)的吸力峰值不斷增大，趨近于實(shí)驗(yàn)結(jié)果。從圖中不同k值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的對比看，當(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馬赫云圖。通過對不同k值模擬結(jié)果的對比容易發(fā)現(xiàn)，當(dāng)k值較大時(shí)，背風(fēng)區(qū)更容易出現(xiàn)分離區(qū)；而隨著k值的減小，流動的分離逐漸被抑制。文獻(xiàn)［10］中的分析認(rèn)為，當(dāng)k取1時(shí)，指數(shù)插值退化為線性插值，采用此指數(shù)次數(shù)模擬邊界層內(nèi)速度可近似為線性分布的小雷諾數(shù)流動更為合理。當(dāng)采用線性分布假設(shè)邊界層內(nèi)速度分布時(shí)，此時(shí)流動抗逆壓梯度的能力較低，更容易出現(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é)果最為接近。我們在早期的相關(guān)研究中也有相似的結(jié)論［17］；但從前面的分析可以發(fā)現(xiàn)，采用本文方法發(fā)展的程序在k＝1／7時(shí)的模擬結(jié)果在翼型背風(fēng)區(qū)會出現(xiàn)明顯的分離區(qū)。本文方法與我們早期的研究中采用的方法主要差別在于邊界附近溫度分布的近似方法，在研究的早期，壁面溫度分布采用等熵關(guān)系式求??；而在本文對壁面附近溫度分布進(jìn)行假設(shè)時(shí)，引入了Crocce－Busemann關(guān)系式。從最終的結(jié)果對比看，Crocce－Busemann關(guān)系式的引入（溫度分布關(guān)系的改變），對模擬結(jié)果的影響較為顯著。文獻(xiàn)［17］中采用等熵關(guān)系式時(shí)，采用較大的k值（如1／5）仍能模擬得到較為準(zhǔn)確的結(jié)果，顯然采用此溫度分布近似時(shí)對流動分離過程會有一定的抑制作用；而本文采用Crocce－Busemann關(guān)系式時(shí)，在較大k值的各模擬狀態(tài)中，背風(fēng)區(qū)流動分離較為明顯。從前文的分析可以發(fā)現(xiàn)，不同溫度分布近似對流動分離過程的影響還值得深入的研究。

Case II狀態(tài)為跨聲速狀態(tài)，雖然來流馬赫數(shù)處于亞聲速，但由于來流相對翼型具有一定的迎角，在翼型的背風(fēng)區(qū)由于流動的膨脹出現(xiàn)超聲速區(qū)，且在翼型背部有明顯的激波存在。圖7給出了采用本文方法模擬得到的Case II在插值指數(shù)次數(shù)k取1／15時(shí)的壓力云圖與馬赫云圖。

從圖8中k取不同值時(shí)，數(shù)值模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對比可以發(fā)現(xiàn)，對于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ù)對比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ū)后緣會出現(xiàn)明顯的流動分離。由于分離的出現(xiàn)，使得背風(fēng)區(qū)前部膨脹流動的膨脹角減小，膨脹區(qū)下游的主流流動速度降低，進(jìn)而使得流動膨脹引起的激波強(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翼型算例對本文方法的網(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ǔ)，對于任意的流動變量φ，其外推解φe計(jì)算方法如下：

將“精確”解通過插值再映射到各套網(wǎng)格上，可求得每套網(wǎng)格下模擬結(jié)果與外推“精確”解之間的誤差。圖10中給出了兩狀態(tài)各流動變量誤差的1范數(shù)變化曲線，其中黑色點(diǎn)線和虛線為示意斜率k＝1和k＝2的輔助線。從各曲線與斜率輔助線的對比可以看出，兩狀態(tài)中各變量在密網(wǎng)格情況下都能達(dá)到二階精度。

圖10 各變量誤差的1范數(shù)曲線Fig.10 Norm－1curves of the variables in the two cases

2.3 圓柱繞流

為檢驗(yàn)本文方法對超聲速流動的模擬能力，采用來流條件為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)中圓柱尾跡為對稱的穩(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í)的空間壓力云圖與流線圖。從圖中可以觀察到清晰的脫體激波、流動繞過圓柱時(shí)形成的膨脹區(qū)以及尾跡區(qū)尾部拖出的膨脹激波。從圖中可以看出，在超聲速條件下，來流馬赫數(shù)越高，其脫體激波的脫體距離越小，但波后壓力越高；圓柱后尾渦的長度也隨著馬赫數(shù)的增大而減小。總體來說，本文方法成功地捕捉到了超聲速條件圓柱繞流的各種典型流動特征。

圖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ù)的對比圖。從圖中可以看出，在圓柱迎風(fēng)面，選取不同k值時(shí)模擬結(jié)果幾乎相同，且數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果基本吻合；但在尾跡區(qū)，尤其是分離起始位置附近，不同k值的模擬結(jié)果存在一定差異。

通過對比可以發(fā)現(xiàn)，在尾跡區(qū)，不同k值的模擬結(jié)果之間的差異主要在于分離起始位置的子午角不同。在分離之前，不同k值模擬結(jié)果幾乎相同；如果分離時(shí)子午角增大，分離前流動速度也增大，壁面壓力越低；但在分離后壁面的壓力分布都將恢復(fù)到幾乎恒定的水平。從對比圖中可以明顯看出，當(dāng)k的取值較小時(shí)，其分離明顯靠后；而較大k值對應(yīng)的分離角明顯靠前。且從兩狀態(tài)的模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比看，當(dāng)k＝1時(shí)模擬結(jié)果與實(shí)驗(yàn)結(jié)果最為接近。從前文的分析可知，k取1時(shí)速度插值相當(dāng)于線性近似，有利于分離流動的模擬；而k取值較小時(shí)，則會抑制分離的出現(xiàn)。顯然，從本算例的模擬結(jié)果中可以得出類似的結(jié)論。在這種低速大分離的尾跡區(qū)域，采用線性近似的效果更好。

圖13 不同k值時(shí)圓柱表面壓力分布與實(shí)驗(yàn)數(shù)據(jù)對比Fig.13 Comparison of the pressure distribution

但從本算例的研究中又引出了新的問題，即當(dāng)前方法插值次數(shù)k的取值為全場統(tǒng)一，這在某些流場結(jié)構(gòu)相對單一的情況，如前述NACA 0012翼型繞流問題，可以得到很好的應(yīng)用。但當(dāng)面對的問題為混合了高／低速流動時(shí)，如何進(jìn)行動態(tài)的指數(shù)插值，是需要進(jìn)一步研究的課題。

3 結(jié) 論

本文發(fā)展了一種基于指數(shù)插值的浸沒邊界法。此方法在壁面附近根據(jù)速度、壓力和溫度等流動量自身特性，分別引入了不同的插值方式，實(shí)現(xiàn)了在非貼體網(wǎng)格上對流場進(jìn)行求解的算法。

采用亞、跨聲速段翼型繞流和超聲速圓柱繞流算例對本文方法進(jìn)行了驗(yàn)證，結(jié)果顯示在選定插值次數(shù)k的條件下，本文方法能很好的模擬亞、跨、超聲速流動問題。顯然本文方法具有很強(qiáng)的可壓縮流動模擬能力。從本文的驗(yàn)證過程可以發(fā)現(xiàn)，當(dāng)k取值較大時(shí)，插值過程接近線性近似，此時(shí)在尾流等存在流動分離或低速流動區(qū)域模擬效果更佳；而k取值較小時(shí)，對激波等強(qiáng)壓縮性區(qū)域的模擬效果更好。

從本文的分析也可以發(fā)現(xiàn)，在流動變量插值關(guān)系式的理論分析，以及此方法在強(qiáng)壓縮性和大分離流動混合的復(fù)雜流動問題的應(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)動固壁超聲速流動的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）陳懋章.粘性流體動力學(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］.空氣動力學(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

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973計(jì)劃，2014CB744100）

何躍龍＊（1985－），男，湖南湘潭人，工程師，博士，研究方向：計(jì)算流體力學(xué).E－mail：aia741＠163.com

何躍龍，李盾，劉帥，等.基于指數(shù)插值的浸沒邊界法在可壓縮流模擬中的應(yīng)用研究［J］.空氣動力學(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.