三維斜流線性完全耦合層吸收邊界條件

伍　新　陳志夫　尹漢鋒

1.湖南大學(xué)汽車車身先進(jìn)設(shè)計制造國家重點實驗室，長沙，4100822.湖南工程學(xué)院，湘潭，4111043.廣州汽車集團(tuán)股份有限公司汽車工程研究院，廣州，511434

三維斜流線性完全耦合層吸收邊界條件

伍新1,2陳志夫3尹漢鋒1

1.湖南大學(xué)汽車車身先進(jìn)設(shè)計制造國家重點實驗室，長沙，4100822.湖南工程學(xué)院，湘潭，4111043.廣州汽車集團(tuán)股份有限公司汽車工程研究院，廣州，511434

采用傅里葉與拉普拉斯變換方法分析了三維斜流背景下聲波、渦波與熵波的色散關(guān)系；根據(jù)各物理波的色散軌跡特征，結(jié)合頻率變化的時空坐標(biāo)變換方法，給出了一組時間與空間坐標(biāo)變換關(guān)系式，并將三維斜流線性歐拉方程變換至新坐標(biāo)系；采用復(fù)數(shù)變換方法，引入阻尼，分別構(gòu)建了x層、y層、z層及角層的完全耦合層(PML)吸收邊界條件,給出了吸收項的施加原則；最后通過三維脈沖聲波、對稱渦環(huán)與周期性點聲源在斜時均流中的傳播問題驗證了該吸收邊界條件的正確性。研究結(jié)果表明：所提出的坐標(biāo)變換關(guān)系能夠有效解決各物理波相位速度與群速度不一致的問題；在斜背景流下，該P(yáng)ML吸收邊界條件能較好地吸收物理波，有效抑制邊界反射，可用于氣動聲學(xué)計算。

完全耦合層；邊界條件；計算氣動聲學(xué)；歐拉方程；色散

0　引言

在開放區(qū)域氣動聲學(xué)問題的數(shù)值計算中，無限計算域需要人工截斷，形成一種特殊的邊界條件，即無反射邊界條件，它既要使計算域內(nèi)的各種物理波能無反射或者較小反射地通過邊界，又能讓計算域外的物理波能順利通過邊界進(jìn)入計算域，同時，還能阻止計算域外的非物理波傳入計算域。完全耦合層(perfectly matched layer,PML)吸收邊界條件作為最優(yōu)秀的無反射邊界條件之一，在近年取得了極大的發(fā)展。

文獻(xiàn)[1-5]首次將Berenger提出的PML技術(shù)引入計算氣動聲學(xué)領(lǐng)域，提出了一系列穩(wěn)定的PML吸收邊界條件。Lin等[6]建立了一個適用于平行流計算的非線性和黏性PML吸收邊界條件。結(jié)合譜差分方法，Zhou等[7]采用該黏性PML吸收邊界條件求解了圓柱繞流等典型氣動聲學(xué)問題。柳占新等[8]從聲學(xué)角度推導(dǎo)了笛卡兒坐標(biāo)系和柱坐標(biāo)下全歐拉方程的PML吸收邊界條件，并將其應(yīng)用于渦扇發(fā)動機(jī)進(jìn)氣道流場模擬。周正干等[9]應(yīng)用PML技術(shù)分析了超聲波聲場特性。Parrish等[10-11]通過兩次坐標(biāo)變換，提出了穩(wěn)定的二維線性和非線性斜流PML吸收邊界條件，并將其拓展至圓柱坐標(biāo)系。雖然上述PML吸收邊界條件具有較高的數(shù)值精度，但是，它們無法求解三維任意方向入流問題。

在汽車側(cè)風(fēng)、圓柱或鈍體斜向繞流等氣動聲學(xué)問題的數(shù)值計算中，往往存在背景斜流。斜流使流場內(nèi)各物理波的群速度與相位速度方向不同，這給穩(wěn)定的PML構(gòu)造帶來了挑戰(zhàn)，即如何建立正確的坐標(biāo)轉(zhuǎn)換關(guān)系[12]。針對該問題，本文從控制方程內(nèi)各物理波的色散關(guān)系軌跡出發(fā)，構(gòu)建恰當(dāng)?shù)淖鴺?biāo)轉(zhuǎn)換關(guān)系，修正各物理波的群速度與相位速度方向，發(fā)展三維斜流線性PML吸收邊界條件。

1　控制方程色散關(guān)系分析與坐標(biāo)變換

1.1控制方程

笛卡兒坐標(biāo)系下，三維量綱一斜流線性歐拉方程為

(1)

其中，ρ為密度，u、v、w分別為x、y、z方向的速度，Max、May、Maz分別為x、y、z方向的馬赫數(shù)，p為壓力，H為非定常源項。

1.2色散關(guān)系分析

根據(jù)波數(shù)分析理論，對式(1)等號兩邊進(jìn)行傅里葉與拉普拉斯變換：

(2)

其中，ω為頻率，kx、ky、kz分別為x、y、z方向的波數(shù)。經(jīng)該變換后即可獲得三維斜流線性歐拉方程的色散關(guān)系矩陣Ψ：

(3)

λ=ω-Maxkx-Mayky-Mazkz

通過求解色散關(guān)系矩陣Ψ的零點λi(i=1,2,…,5)即可獲得各物理波的色散關(guān)系。熵波與渦波的色散關(guān)系為

ω-Maxkx-Mayky-Mazkz=0

(4)

聲波的色散關(guān)系為

(5)

為了分析各物理波的群速度與相位速度方向，將式(4)與式(5)等號兩邊同時除以頻率ω，得

(6)

(7)

渦波、熵波及聲波在坐標(biāo)系(kx/ω,ky/ω,kz/ω)內(nèi)的色散關(guān)系軌跡如圖1a所示。圖1a中，渦波與熵波的色散關(guān)系軌跡為一個斜置的平面，聲波的色散關(guān)系軌跡為一個中心不在坐標(biāo)系原點的斜置橢球。根據(jù)色散關(guān)系軌跡穩(wěn)定性分析方法[13]，三維斜流線性歐拉方程中，不論是渦波、熵波還是聲波，其色散關(guān)系軌跡均存在不穩(wěn)定區(qū)域，因此，在進(jìn)行PML變換之前，必須采用適當(dāng)?shù)淖鴺?biāo)變換來修正域內(nèi)各物理波的群速度與相位速度方向。

1.3坐標(biāo)變換關(guān)系的構(gòu)建

如果將z方向的斜流消去，則該問題可以轉(zhuǎn)化為二維斜流線性歐拉方程的坐標(biāo)變換關(guān)系式的推導(dǎo)。參考文獻(xiàn)[13]中的結(jié)論，x層坐標(biāo)變換的關(guān)系式可以表示為

(8)

式中，t為時間。

由式(8)可獲得各物理量關(guān)于新坐標(biāo)與原坐標(biāo)偏導(dǎo)數(shù)之間的關(guān)系:

(9)

與式(9)相對應(yīng)的新波數(shù)、新頻率與原波數(shù)、原頻率之間的關(guān)系為

(10)

(11)

聲波的色散關(guān)系為

(12)

同理，可以獲得用于y層與z層PML吸收邊界條件推導(dǎo)的坐標(biāo)變換關(guān)系式，其表達(dá)式分別為

(13)

(14)

(a)變換前(b)x層

(c)y層(d)z層圖1　物理波在坐標(biāo)系(kx/ω,ky/ω,kz/ω)與內(nèi)的色散關(guān)系軌跡

2　PML吸收邊界條件推導(dǎo)

采用復(fù)數(shù)變換方法，通過式(8)、式(13)與式(14)分別推導(dǎo)x層、y層及z層的PML吸收邊界條件。然后，根據(jù)各層PML吸收邊界條件推導(dǎo)各角層PML吸收邊界條件。對于三維計算域，整個計算域由26塊PML吸收域與一個物理計算域組成，PML吸收邊界條件由七大部分組成，分別為x層、y層、z層、xy層、yz層、xz層及xyz層，如圖2所示。

圖2　物理域與PML計算域示意圖

2.1平行層PML吸收邊界條件

首先，將式(9)代入式(1)，將控制方程變換至新坐標(biāo)系：

(15)

其次，引入阻尼，將新時間、空間坐標(biāo)系下的控制方程轉(zhuǎn)換至頻域：

(16)

再次，引入輔助變量q1，將頻域PML吸收邊界條件變換至新的時間與空間坐標(biāo)系：

(17)

(18)

最后，通過式(9)將新時間、空間坐標(biāo)系下的PML吸收邊界條件變回至原時間與空間坐標(biāo)系，獲得x層PML吸收邊界條件：

(19)

(20)

同理，也可獲得y層與z層PML吸收邊界條件：

(21)

(22)

(23)

(24)

2.2角層PML吸收邊界條件

關(guān)于xy層PML吸收邊界條件，吸收系數(shù)σx與σy均不為零，因此結(jié)合x層和y層PML吸收邊界條件即可得xy層PML吸收邊界條件，其表達(dá)式為

(25)

(26)

(27)

為了使xy層PML吸收邊界條件穩(wěn)定，在輔助變量方程(式(26)與式(27))等號左邊，分別增加額外的吸收項σyq1與σxq2，該吸收項對穩(wěn)定性的影響詳見文獻(xiàn)[13]。

同理，也可獲得yz層和xz層PML吸收邊界條件，其表達(dá)式分別為

(28)

(29)

(30)

(31)

(32)

(33)

為使各角層PML吸收邊界條件穩(wěn)定，在式(29)、式(30)、式(32)與式(33)等號的左邊依次增加吸收項σzq2、σyq3、σzq1與σxq3。

關(guān)于xyz層PML吸收邊界條件，吸收系數(shù)σx、σy與σz均不為零，因此必須結(jié)合x層、y層和z層PML吸收邊界條件即可得xyz層PML吸收邊界條件，其表達(dá)式為

[σx(I+βxA)+σy(I+βyB)+σz(I+βzC)]u=0

(34)

(35)

(36)

(37)

式(35)～式(37)中，(σy+σz)q1、(σx+σz)q2與(σx+σy)q3分別為輔助變量方程中額外增加的吸收項。

3　算例驗證

為驗證本文構(gòu)建的三維斜流線性PML吸收邊界條件的正確性，選用三維脈沖聲波、對稱渦環(huán)與周期性點聲源在斜時均流中的傳播問題作為測試算例。計算域四周均采用PML吸收邊界條件，空間離散格式采用改進(jìn)的7點色散保持有限差分格式[13]，時間推進(jìn)格式采用6級4階低耗散低色散RK顯式格式[14]。為了消除短波對數(shù)值計算的影響，采用人工黏性耗散[15]進(jìn)行計算。計算時空間網(wǎng)格尺寸Δx=Δy=Δz=1，時間步長取0.1 s，吸收系數(shù)σx、σy、σz分別為

(38)

D為PML吸收寬度，xb、yb、zb分別為PML計算域與歐拉計算域的交界位置，根據(jù)文獻(xiàn)[10]，σm、α分別取值2、3。

3.1脈沖聲波傳播

假定該脈沖聲波初始位置在坐標(biāo)原點，并且處于速度場為(0.5,0.5,0.5)的時均流場中，該問題的初始條件可描述為

(39)

設(shè)該問題的歐拉計算域為x,y,z∈[-20,20]，吸收寬度取值為D=10Δx=10Δy=10Δz，在t=40 s時的壓力等值面及截面等值線如圖3所示?？梢钥闯?，當(dāng)三維脈沖聲波達(dá)到PML吸收邊界時，聲波迅速衰減，在邊界處未見明顯反射。

(a)壓力等值面圖(b)平面z=20上的壓力等值線圖圖3　t=40 s時壓力圖

3.2渦環(huán)傳播

渦環(huán)動力學(xué)及其發(fā)聲機(jī)理是目前一個非?；钴S的研究方向，也是理論流體力學(xué)與氣動聲學(xué)領(lǐng)域中的難點。以對稱渦環(huán)在斜時均流中的傳播問題作為測試算例進(jìn)行測試。設(shè)渦環(huán)的半徑r0=10，渦核的半徑b=3，初始時刻渦環(huán)對稱中心位于坐標(biāo)系原點，將該渦環(huán)置于(0.5,0.5,0.5)的時均流場中，如圖4所示，該問題的初始條件可描述為

(40)

ux(r)=εr0r-1(r-r0)e-α[x2+(r-r0)2]

ur(r)=-εr0r-1(r-r0)xe-α[x2+(r-r0)2]

圖4　對稱渦環(huán)示意圖

采用與3.1節(jié)中相同的計算方案，獲得渦量隨時間的變化情況，如圖5所示。該渦環(huán)隨時均流對流，可以看出，當(dāng)渦環(huán)到達(dá)PML計算域內(nèi)時，渦量迅速衰減，且未見明顯反射。由此可知，當(dāng)任意方向斜流存在時，渦環(huán)能順利穿過計算邊界，該P(yáng)ML吸收邊界條件能較好地吸收渦環(huán)中的物理波。

(a)t=0時渦量等值面圖(b)t=0時平面x=0上的渦量等值線圖

(c)t=40 s時渦量等值面圖(d)t=40 s時平面x=20上的渦量等值線圖圖5　渦量隨時間變化圖

3.3周期性點聲源傳播

為了測試PML吸收邊界條件對周期性聲波的吸收性能，在式(1)中的壓力方程右端施加一個周期性點聲源，該點聲源為

s(x,y,z,t)=sin(0.2πt)e-(ln2)[x2+y2+z2]/32

(41)

歐拉計算域為x,y,z∈[-40,40]，計算域外圍采用10層PML吸收邊界條件。

當(dāng)t=100s時，壓力數(shù)值解等值面與等值線如圖6所示。可以看出，周期性聲波能無反射地通過邊界。

(a)聲壓等值面圖(b)平面x=0上的聲壓等值線圖

(c)平面y=0上的聲壓等值線圖(d)平面z=0上的聲壓等值線圖圖6　t=100 s時聲壓等值面與等值線圖

4　結(jié)論

(1)三維斜流線性歐拉方程中，渦波、熵波與聲波的色散關(guān)系軌跡均存在不穩(wěn)定區(qū)域。

(2)提出的坐標(biāo)變換關(guān)系式能有效修正各物理波的群速度與相位速度方向，使其色散關(guān)系軌跡均穩(wěn)定。

(3)提出的PML吸收邊界條件能較好地吸收熵波、渦波與聲波，未見明顯反射。

[1]BerengerJP.APerfectlyMatchedLayerfortheAbsorptionofElectromagneticWaves[J].JournalofComputationalPhysics, 1994, 114(2): 185-200.

[2]HuFQ.OnAbsorbingBoundaryConditionsforLinearizedEulerEquationsbyaPerfectlyMatchedLayer[J].JournalofComputationalPhysics, 1996, 129(1): 201-219.

[3]HuFQ.AStablePerfectlyMatchedLayerforLinearizedEulerEquationsinUnsplitPhysicalVariables[J].JournalofComputationalPhysics, 2001, 173(2): 455-480.

[4]HuFQ.APerfectlyMatchedLayerAbsorbingBoundaryConditionforLinearizedEulerEquationswithaNon-uniformMeanFlow[J].JournalofComputationalPhysics, 2005, 208(2): 469-492.

[5]HuFQ,LiXD,LinDK.AbsorbingBoundaryConditionforNonlinearEulerandNavier-stokesEquationsBasedonthePerfectlyMatchedLayerTechnique[J].JournalofComputationalPhysics, 2008, 227(9): 4398-4424.

[6]LinDK,LiXD,HuFQ.AbsorbingBoundaryConditionforNonlinearEulerEquationsinPrimitiveVariablesBasedonthePerfectlyMatchedLayerTechnique[J].Computers&Fluids, 2011, 40(1): 333-337.

[7]ZhouY,WangZJ.AbsorbingBoundaryConditionsfortheEulerandNavier-stokesEquationswiththeSpectralDifferenceMethod[J].JournalofComputationalPhysics, 2010, 229(23): 8733-8749.

[8]柳占新, 高頻, 仝志勇.全歐拉方程的理想匹配層邊界條件[J].中國機(jī)械工程, 2011, 22(16): 1938-1941.

LiuZhanxin,GaoPin,TongZhiyong.PMLBoundaryConditionsforFullEulerEquations[J].ChinaMechanicalEngineering, 2011, 22(16): 1938-1941.

[9]周正干, 魏東.時域有限差分法在超聲波聲場特性分析中的應(yīng)用[J].機(jī)械工程學(xué)報，2010, 46(2): 9-13.

ZhouZhenggan,WeiDong.AnalysisofUltrasonicSoundFieldCharacteristicwithFDTD[J].ChineseJournalofMechanicalEngineering，2010, 46(2): 9-13.

[10]ParrishSA,HuFQ.PMLAbsorbingBoundaryConditionsfortheLinearizedandNonlinearEulerEquationsintheCaseofObliqueMeanFlow[J].InternationalJournalforNumericalMethodsinFluids, 2009, 60(5): 565-589.

[11]ParrishSA.AnalysisandApplicationofPerfectlyMatchedLayerAbsorbingBoundaryConditionsforComputationalAeroacoustics[D].Norfolk:Univ.ofOldDominion, 2008.

[12]HuFQ.DevelopmentofPMLAbsorbingBoundaryConditionsforComputationalAeroacoustics:aProgressReview[J].Computers&Fluids, 2008, 37(4): 336-348.

[13]陳志夫.基于色散關(guān)系分析的高精度氣動聲學(xué)計算方法研究[D].長沙:湖南大學(xué), 2013.

[14]HuFQ,HussainiMY,MantheyJL.Low-dissipationandLow-dispersionRunge-KuttaSchemesforComputationalAcoustics[J].JournalofComputationalPhysics, 1996, 124(1): 177-191.

[15]TamCKW,JayC,ZhongD.AStudyoftheShortWaveComponentsinComputationalAcoustics[J].JournalofComputationalAcoustics, 1993, 1(1): 1-30.

(編輯陳勇)

Three Dimensional Linear PML Absorbing Boundary Conditions with an Oblique Mean Flow

Wu Xin1,2Chen Zhifu3Yin Hanfeng1

1.State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University,Changsha,410082 2.Hunan Institute of Engineering,Xiangtan,Hunan,411104 3.Guangzhou Automobile Group Co.,Ltd. Automotive Engineering Institute,Guangzhou,511434

For three dimensional linear Euler equations in the case of oblique mean flow, the dispersion relations of acoustic, vortex and entropy wave were first analyzed by using Fourier and Laplace transform method. Then the hypothesis for changed frequency was employed, a proper space-time transformation was presented for deriving three dimensional linear Euler equations in transformed coordinates. A complex change was applied to the new equations and a damping parameter was introduced. A three linear PML absorbing boundary conditions in the case of oblique mean flow forxlayer,ylayer,zlayer and corner layer were derived. In addition, the importance of added absorption term was emphasized. Finally, the effectiveness of linear PML absorbing boundary conditions was validated by computing the computational aeroacoustics benchmark problems. The results prove that: the presented space-time transformation can solve the problem of direction inconsistence in group and phase velocity of physical wave; in the case of oblique mean flow, the proposed PML absorbing boundary conditions can absorb the physical wave with little or no reflection.Therefore,it also can be applied to aeroacoustic computation.

perfectly matched layer(PML); boundary condition; computational aeroacoustics; Euler equation; dispersion

2014-05-19

國家自然科學(xué)基金資助項目(11302075,11002052)；湖南省教育廳高等學(xué)校科學(xué)研究項目(12C0627)

V211.3DOI：10.3969/j.issn.1004-132X.2015.01.001

伍新，男，1976年生。湖南大學(xué)汽車車身先進(jìn)設(shè)計制造國家重點實驗室博士研究生，湖南工程學(xué)院機(jī)械工程學(xué)院講師。主要研究方向為計算聲學(xué)、振動與噪聲控制。陳志夫，男，1986年生。廣州汽車集團(tuán)股份有限公司汽車工程研究院工程師。尹漢鋒，男，1982年生。湖南大學(xué)汽車車身先進(jìn)設(shè)計制造國家重點實驗室講師、博士。