伍　新　陳志夫　尹漢鋒

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

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

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

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

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

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

0　引言

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

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

在汽車側(cè)風(fēng)、圓柱或鈍體斜向繞流等氣動(dòng)聲學(xué)問(wèn)題的數(shù)值計(jì)算中，往往存在背景斜流。斜流使流場(chǎng)內(nèi)各物理波的群速度與相位速度方向不同，這給穩(wěn)定的PML構(gòu)造帶來(lái)了挑戰(zhàn)，即如何建立正確的坐標(biāo)轉(zhuǎn)換關(guān)系[12]。針對(duì)該問(wèn)題，本文從控制方程內(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為非定常源項(xiàng)。

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

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

(2)

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

(3)

λ=ω-Maxkx-Mayky-Mazkz

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

ω-Maxkx-Mayky-Mazkz=0

(4)

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

(5)

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

(6)

(7)

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

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

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

(8)

式中，t為時(shí)間。

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

(9)

與式(9)相對(duì)應(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ù)變換方法，通過(guò)式(8)、式(13)與式(14)分別推導(dǎo)x層、y層及z層的PML吸收邊界條件。然后，根據(jù)各層PML吸收邊界條件推導(dǎo)各角層PML吸收邊界條件。對(duì)于三維計(jì)算域，整個(gè)計(jì)算域由26塊PML吸收域與一個(gè)物理計(jì)算域組成，PML吸收邊界條件由七大部分組成，分別為x層、y層、z層、xy層、yz層、xz層及xyz層，如圖2所示。

圖2　物理域與PML計(jì)算域示意圖

2.1平行層PML吸收邊界條件

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

(15)

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

(16)

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

(17)

(18)

最后，通過(guò)式(9)將新時(shí)間、空間坐標(biāo)系下的PML吸收邊界條件變回至原時(shí)間與空間坐標(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(chóng)層和y層PML吸收邊界條件即可得xy層PML吸收邊界條件，其表達(dá)式為

(25)

(26)

(27)

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

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

(28)

(29)

(30)

(31)

(32)

(33)

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

關(guān)于xyz層PML吸收邊界條件，吸收系數(shù)σx、σy與σz均不為零，因此必須結(jié)合x(chóng)層、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分別為輔助變量方程中額外增加的吸收項(xiàng)。

3　算例驗(yàn)證

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

(38)

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

3.1脈沖聲波傳播

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

(39)

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

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

3.2渦環(huán)傳播

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

(40)

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

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

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

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

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

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

3.3周期性點(diǎn)聲源傳播

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

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

(41)

歐拉計(jì)算域?yàn)閤,y,z∈[-40,40]，計(jì)算域外圍采用10層PML吸收邊界條件。

當(dāng)t=100s時(shí)，壓力數(shù)值解等值面與等值線如圖6所示?？梢钥闯?，周期性聲波能無(wú)反射地通過(guò)邊界。

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

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

4　結(jié)論

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

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

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

[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].中國(guó)機(jī)械工程, 2011, 22(16): 1938-1941.

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

[9]周正干, 魏東.時(shí)域有限差分法在超聲波聲場(chǎng)特性分析中的應(yīng)用[J].機(jī)械工程學(xué)報(bào)，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)系分析的高精度氣動(dòng)聲學(xué)計(jì)算方法研究[D].長(zhǎng)沙:湖南大學(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

國(guó)家自然科學(xué)基金資助項(xiàng)目(11302075,11002052)；湖南省教育廳高等學(xué)?？茖W(xué)研究項(xiàng)目(12C0627)

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

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