王 雅，蔣靈杰，鄧小龍,2

（1北京計(jì)算科學(xué)研究中心，北京 100193；2弗吉尼亞大學(xué)，美國(guó)弗吉尼亞州夏洛茨維爾 22904）

高速氣體和固體顆粒群的相互作用是一類典型的可壓縮多相流問(wèn)題，廣泛存在于天文、自然災(zāi)害、工業(yè)安全、醫(yī)療工業(yè)和國(guó)防等領(lǐng)域，如超新星爆炸[1]、火山爆發(fā)[2]、粉塵爆炸[3]、無(wú)針注射[4]和炸彈爆炸[5]等。在高速顆粒流中，顆粒體積分?jǐn)?shù)αd是一個(gè)重要參數(shù)[6]。當(dāng)αd? 1 時(shí)，顆粒之間彼此遠(yuǎn)離，顆粒間的碰撞效應(yīng)可以忽略不計(jì)[7]；當(dāng)αd≥ 0.5時(shí)，顆粒之間彼此靠近，顆粒間的碰撞是其運(yùn)動(dòng)的主要機(jī)制，流體對(duì)固體顆粒的作用可以忽略不計(jì)[8]。本研究的固體體積分?jǐn)?shù)為15%，屬于0.001 ＜αd＜ 0.5范圍[8]，此時(shí)流體與顆粒的相互作用以及顆粒之間的相互作用變得尤為重要[9-10]。

在物理實(shí)驗(yàn)方面，Rogue等[11]在垂直激波管中進(jìn)行了激波與水平顆粒床相互作用的實(shí)驗(yàn)。Wagner等[12]使用多相流激波管（Multiphase shock tube,MST）對(duì)激波與顆粒簾的相互作用進(jìn)行實(shí)驗(yàn)研究，觀察到反射激波和透射激波的產(chǎn)生和傳播，以及顆粒簾寬度擴(kuò)張和傳播過(guò)程。Wagner等[13]利用X射線測(cè)量技術(shù)改進(jìn)了激波管，以觀察在整個(gè)實(shí)驗(yàn)過(guò)程中顆粒簾內(nèi)部體積分?jǐn)?shù)的分布變化。Theofanous等[14-15]使用ASOS激波管研究激波與顆粒簾的相互作用，發(fā)現(xiàn)了此類實(shí)驗(yàn)中存在時(shí)間標(biāo)度（Time scaling）現(xiàn)象，進(jìn)而總結(jié)出顆粒簾寬度在實(shí)驗(yàn)前期加速擴(kuò)張、實(shí)驗(yàn)后期勻速擴(kuò)張的規(guī)律。

在數(shù)值模擬方面，研究顆粒流的常用方法有點(diǎn)顆粒模型法[16-17]和直接數(shù)值模擬（Direct numerical simulation,DNS）方法，DNS方法又包括基于網(wǎng)格重構(gòu)的拉格朗日-歐拉移動(dòng)網(wǎng)格法（Arbitray Lagrangian-Eulerian,ALE）[18]和浸入邊界法（Immersed boundary method,IBM）[19-21]等。Zhu等[22]使用界面解析的DNS方法研究了不同形狀（橢球和圓球）的顆粒在湍流通道中的流動(dòng)。Zastawny等[23]使用改進(jìn)的鏡像浸入邊界法進(jìn)行DNS，推導(dǎo)了流動(dòng)中4種非球形顆粒的阻力和升力系數(shù)以及扭矩系數(shù)。鄒立勇等[24]實(shí)驗(yàn)研究了在馬赫數(shù)為1.18的平面激波沖擊作用下，雙橢圓界面R-M不穩(wěn)定性演化的動(dòng)力學(xué)過(guò)程。2018年Jiang等[25]基于分層流模型[26-27]在笛卡爾背景網(wǎng)格下，進(jìn)行了高速氣體與二維圓柱云相互作用前期的系統(tǒng)性數(shù)值模擬，并對(duì)高速氣體與三維圓球云的數(shù)值模擬進(jìn)行了初步探究。Deng等[28]進(jìn)一步探討了高速氣體與三維圓球云相互作用后期的相關(guān)規(guī)律，發(fā)現(xiàn)了與Theofanous等[14-15]實(shí)驗(yàn)中獲得的時(shí)間標(biāo)度類似的現(xiàn)象。

本研究基于分層流模型[26-27]，在前人[25,28]的基礎(chǔ)上拓展研究，對(duì)平面激波與橢圓柱云的相互作用進(jìn)行DNS，重點(diǎn)關(guān)注橢圓柱橫截面不同長(zhǎng)短軸之比（λ）以及橢圓柱橫截面長(zhǎng)軸與來(lái)流方向成不同角度（θ）時(shí)對(duì)流場(chǎng)的影響程度。圖1為橢圓柱橫截面幾何示意圖，其中a為橢圓柱橫截面長(zhǎng)軸長(zhǎng)，b為橢圓柱橫截面短軸長(zhǎng)，λ=a/b表示長(zhǎng)短軸之比，θ表示x軸（流場(chǎng)來(lái)流方向）與長(zhǎng)軸之間的夾角。

圖1 橢圓柱橫截面幾何示意圖Fig.1 Illustration of the geometry for the cross-section of the elliptical cylinder

1 數(shù)值方法

本研究采用的數(shù)值方法是基于Chang和Liou[27]提出的分層流模型，其控制方程如下

式中：下標(biāo)“g”表示氣相；ρ為密度；v為速度矢量；n為控制體邊界的單位外法向量；p為壓力；E和H分別為每個(gè)控制體的總能量和總焓；Sg為控制體中氣體相的表面積；Vg=αgVi，其中Vg為控制體中的氣相體積，αg為氣體的體積分?jǐn)?shù)，Vi為控制體的體積。應(yīng)用理想氣體狀態(tài)方程來(lái)封閉式（1）。采用有限體積法（Finite volume method,FVM）離散控制方程，空間重構(gòu)使用三階TVD格式，由于不同控制體的體積分?jǐn)?shù)αg不都相同，每個(gè)控制體界面可以重構(gòu)為氣-氣、固-固和氣-固3個(gè)部分，如圖2所示。其中氣-氣界面之間的通量使用AUSM+-up近似黎曼解法器[27,29-30]計(jì)算，時(shí)間推進(jìn)采用三階龍格庫(kù)塔（Runge-Kutta）方法，計(jì)算網(wǎng)格使用笛卡爾網(wǎng)格。本研究主要針對(duì)激波與橢圓柱云相互作用的前期階段，可認(rèn)為此時(shí)橢圓柱固定不動(dòng)。由于流場(chǎng)的流速較高，因此橢圓柱所受的合外力僅通過(guò)圓柱表面對(duì)壓力積分獲得，而忽略流體黏性對(duì)橢圓柱受力的影響。氣-固界面采用滑移邊界條件。使用MPI實(shí)現(xiàn)并行計(jì)算，以便進(jìn)行大規(guī)模數(shù)值模擬。

2 平面激波與橢圓柱云相互作用研究

2.1 方法驗(yàn)證

通過(guò)測(cè)試網(wǎng)格的收斂性檢驗(yàn)本方法的正確性。采用4種不同分辨率的網(wǎng)格進(jìn)行數(shù)值模擬。表1列出了不同分辨率下x和y方向上的網(wǎng)格數(shù)Nx和Ny，以及解析橢圓短軸所使用的網(wǎng)格數(shù)nb，計(jì)算區(qū)域設(shè)置見(jiàn)圖3(a)。橢圓圓心位于原點(diǎn) (0 ,0)處，入射激波位于x=-0.04處，計(jì)算域?yàn)閤∈[-0.140，0.140]，y∈[-0.084，0.084]，z∈[-0.020，0.020]。為了節(jié)省計(jì)算資源，當(dāng)x∈[-0.056，0.056]且y∈[-0.028，0.028]時(shí)，網(wǎng)格間距是均勻的，且在該區(qū)域Δx=Δy；除此之外，x和y方向采用不等間距的拉伸網(wǎng)格，z方向上的均勻網(wǎng)格個(gè)數(shù)設(shè)置為2，且該方向上兩側(cè)為周期性邊界條件。圖3(a)的上下邊界條件為周期性邊界條件，左邊界為入口邊界條件，右邊界為出口邊界條件，入口邊界和出口邊界都按照氣體的初始條件設(shè)置為固定值。氣體初始條件為：波前 (p,T,u)R=(8.234 9×104Pa,294.9 K,0.0 m/s)，波后 (p,T,u)L=(2.548 9×105Pa,425.2 K,309.1 m/s)。

表1 網(wǎng)格收斂性分析實(shí)驗(yàn)中使用的4種網(wǎng)格Table 1 Four meshes used in the convergence analysis experiment

圖3 網(wǎng)格收斂性分析實(shí)驗(yàn)示意圖Fig.3 Illustration of the convergence analysis experiment

在此初始條件下可以產(chǎn)生馬赫數(shù)Ma=1.67的運(yùn)動(dòng)激波。圖3(b)顯示了網(wǎng)格收斂性的數(shù)值模擬結(jié)果，主要考察了x方向上橢圓柱所受的外力Fx隨時(shí)間的演化規(guī)律。從圖3(b)中可以看出：隨著nb的增加，數(shù)值模擬結(jié)果顯示出很好的收斂性，且當(dāng)nb=32時(shí)，數(shù)值模擬結(jié)果與nb=64時(shí)的結(jié)果吻合很好。因此，以下均使用nb=32的網(wǎng)格進(jìn)行DNS。

2.2 直接數(shù)值模擬

設(shè)入射激波馬赫數(shù)Ma=1.67，氣體初始條件與2.1節(jié)中的初始條件相同，網(wǎng)格分辨率nb取為32。初始流場(chǎng)設(shè)置見(jiàn)圖4，入射激波位于x=-2.5，計(jì)算域?yàn)閤∈[-3.0，4.0]，y∈[-0.5，0.5]，z∈[-0.1，0.0]。為了節(jié)省計(jì)算資源，當(dāng)x∈[-0.65，0.65]時(shí)，采用等間距網(wǎng)格，除此之外，x方向采用拉伸網(wǎng)格，y方向采用均勻網(wǎng)格，z方向上的網(wǎng)格個(gè)數(shù)設(shè)置為1。橢圓柱云位于x∈[-0.5，0.5]區(qū)域，其寬度L=1。橢圓柱個(gè)數(shù)Np=440，橢圓柱的橫截面積均相同，其排布參照J(rèn)iang等[25]的排布方案。表2列出了λ為2、3和4時(shí)網(wǎng)格的設(shè)置，a為橢圓柱橫截面長(zhǎng)軸長(zhǎng)，b為橢圓柱橫截面短軸長(zhǎng)，Δx為均勻網(wǎng)格區(qū)域網(wǎng)格的寬度，Nx和Ny分別表示x和y方向上的網(wǎng)格總數(shù)，N為整個(gè)計(jì)算區(qū)域內(nèi)的網(wǎng)格總數(shù)。

圖4 x-y平面計(jì)算區(qū)域設(shè)置示意圖（右圖為初始橢圓柱云分布圖，藍(lán)色表示低壓區(qū)域，紅色表示高壓區(qū)域）Fig.4 Illustration of the computational domain setting in the x-y plane (The right plot shows the initial distribution of the elliptical cylinder cloudThe red and blue regions represent the high-pressure and low-pressure regions,respectively.)

表2 平面激波與橢圓柱云相互作用數(shù)值模擬使用的網(wǎng)格設(shè)置Table 2 Mesh settings in numerical simulation of the interaction between plane shock and elliptical column cloud

在數(shù)值模擬過(guò)程中，可以觀察到類似于高速氣體與圓柱云相互作用的物理現(xiàn)象[25]。圖5顯示了當(dāng)λ=2、θ=0°時(shí)，無(wú)量綱時(shí)間t為1.3、1.5、1.8、2.4和3.5時(shí)流場(chǎng)中的壓強(qiáng)分布，其中壓強(qiáng)采用初始時(shí)刻波前壓強(qiáng)（8.234 9×104Pa）進(jìn)行了無(wú)量綱化處理。從圖5中可以看出：在激波與橢圓柱云相互作用的過(guò)程中產(chǎn)生了一道向流場(chǎng)上游運(yùn)動(dòng)的反射激波和一道向流場(chǎng)下游運(yùn)動(dòng)的透射激波；當(dāng)無(wú)量綱時(shí)間為2.4和3.5時(shí)，在橢圓柱云內(nèi)部和流場(chǎng)下游部分區(qū)域流場(chǎng)擾動(dòng)較大，原因是此區(qū)域存在激波的反射以及激波之間的相互作用。

圖5 當(dāng)λ=2、θ=0°時(shí)，不同無(wú)量綱時(shí)間下流場(chǎng)的無(wú)量綱壓強(qiáng)分布Fig.5 Distributions of the dimensionless pressure at different dimensionless time when λ=2 and θ=0°

為了定量分析流場(chǎng)，先給出本研究所涉及的部分變量的定義[8]，變量 φ的體積平均可以定義為

式中：V為連續(xù)相和離散相的總體積。連續(xù)相的相平均（或雷諾平均）定義為

式中：Vc為V中連續(xù)相的體積。質(zhì)量平均定義為

式中：下標(biāo)i表示速度的3個(gè)方向表示速度的質(zhì)量平均表示的平方。連續(xù)相的體積平均后的內(nèi)能、動(dòng)能和湍動(dòng)能的定義[8]分別為

為了便于描述，將計(jì)算域分為3部分，即上游區(qū)域、橢圓柱云區(qū)域和下游區(qū)域，如圖6(e)所示，其中橢圓柱云的邊界分為橢圓柱云上游邊界（UFC）和橢圓柱云下游邊界（DFC）。圖6展示了在t=3.5，λ為2、3、4時(shí)，橢圓柱橫截面長(zhǎng)軸與來(lái)流方向呈不同角度時(shí)的流場(chǎng)速度（u）、流場(chǎng)內(nèi)能和流場(chǎng)動(dòng)能在計(jì)算域中沿著x方向的分布。從圖6中可以看出：隨著θ從0°增大到135°，入射激波與橢圓柱云正面沖擊的有效面積先增大后減小，橢圓柱云對(duì)入射激波的反射效果也先增強(qiáng)后減弱；當(dāng)θ達(dá)到90°時(shí)，入射激波與橢圓柱云正面沖擊的有效面積最大，橢圓柱云對(duì)入射激波的反射效果最強(qiáng)。從圖6中來(lái)流速度u的分布可以看出，當(dāng)θ=90°時(shí)，反射激波位置離UFC最遠(yuǎn)，而透射激波位置離DFC最近。從內(nèi)能分布圖中也可以看出，與其他角度相比，當(dāng)θ=90°時(shí)，反射激波與UFC之間的區(qū)域中的內(nèi)能相對(duì)于初始時(shí)刻的內(nèi)能增加更大，而下游區(qū)域中內(nèi)能增加最小。從動(dòng)能分布中可以看出，當(dāng)θ=90°時(shí)，入射激波對(duì)下游的影響相比其他角度來(lái)說(shuō)更小，說(shuō)明流場(chǎng)中的能量更多以內(nèi)能形式保留在從反射激波面到UFC之間的區(qū)域中，而較少輸入流場(chǎng)的下游。由于橢圓柱云的分布具有沿中心軸上下近似對(duì)稱的性質(zhì)，如圖4所示，因此當(dāng)θ=45°和135°時(shí)，流場(chǎng)中主流速度(u)、內(nèi)能和動(dòng)能的分布也具有相似性。

圖6 t=3.5時(shí)不同λ下θ分別為0°、45°、90°、135°時(shí)的流場(chǎng)速度、流場(chǎng)內(nèi)能和流場(chǎng)動(dòng)能分布（灰色矩形區(qū)域表示橢圓柱云，RS、TS、UFC、DFC分別表示反射激波、透射激波、橢圓柱云上游邊界、橢圓柱云下游邊界）Fig.6 Distributions of the fluid velocity,fluid internal energy and fluid kinetic energy with different λ when θ equals to 0°,45°,90°,135° at dimensionless time t=3.5 (The gray rectangular regions stand for the elliptical cylinder cloud.Hereafter,RS,TS,UFC and DFC mean reflected shock,transmitted shock,the upstream front of elliptical column cloud,and the downstream front of elliptical column cloud,respectively.)

圖7顯示了在無(wú)量綱時(shí)間t=3.5且θ不變，λ分別為2、3、4時(shí)的主流速度、流場(chǎng)內(nèi)能和流場(chǎng)動(dòng)能在x方向上的分布。當(dāng)θ=0°時(shí)，隨著λ從2增加到4，相同面積的橢圓逐漸變得細(xì)長(zhǎng)，入射激波與橢圓柱云正面沖擊的有效面積減小，橢圓柱云對(duì)入射激波的反射效果減弱。從圖7中來(lái)流速度u的分布來(lái)看，當(dāng)θ=0°、λ=2時(shí)，反射激波間斷面離UFC最遠(yuǎn)，透射激波的間斷面離DFC最近。從內(nèi)能分布圖中也可以看出，隨著λ的增大，反射激波與UFC之間的區(qū)域中的內(nèi)能相對(duì)于初始時(shí)刻內(nèi)能上升的幅度減小，傳遞到下游的內(nèi)能相對(duì)于初始時(shí)刻內(nèi)能逐漸增加，流場(chǎng)下游區(qū)域的內(nèi)能占流場(chǎng)總內(nèi)能的比例逐漸增加。從動(dòng)能分布來(lái)看，隨著λ的增加，動(dòng)能整體增加。而當(dāng)θ=90°時(shí)，隨著λ的增大，反射激波與UFC之間的區(qū)域中的內(nèi)能相對(duì)于初始時(shí)刻內(nèi)能上升的幅度增加，傳遞到下游的內(nèi)能相對(duì)于初始時(shí)刻內(nèi)能逐漸減少，流場(chǎng)上游區(qū)域中的內(nèi)能占流場(chǎng)總內(nèi)能的比例逐漸增加，且此時(shí)流場(chǎng)中的動(dòng)能隨著λ的增加而減小，但差別不大。當(dāng)θ=45°時(shí)，不同的λ下，入射激波與橢圓柱云正面沖擊的有效面積差異較小，因此λ的不同對(duì)主流速度、流場(chǎng)內(nèi)能和流場(chǎng)動(dòng)能的影響較小。

圖7 t=3.5時(shí)不同θ下λ分別為2、3、4時(shí)，流場(chǎng)速度、流場(chǎng)內(nèi)能和流場(chǎng)動(dòng)能的分布（灰色矩形區(qū)域表示橢圓柱云）Fig.7 Distributions of the fluid velocity,fluid internal energy and fluid kinetic energy with different θ,when λ equals to 2,3,4,at dimensionless time t=3.5,where the gray rectangular regions stand for the elliptical cylinder cloud

圖8顯示了在無(wú)量綱時(shí)間t=3.5時(shí)，相同θ下，當(dāng)λ從2增加到4，沿x和y方向上的RMS速度（u′′，v′′）和湍動(dòng)能的分布，虛線之間的區(qū)域表示橢圓柱云?？梢杂^察到，當(dāng)θ為0°、45°和135°時(shí)，λ的變化對(duì)u′′、v′′和湍動(dòng)能分布的影響較小，此時(shí)u′′、v′′和湍動(dòng)能數(shù)值較大的區(qū)域分布在DFC附近。當(dāng)θ=90°時(shí)，隨著λ的增加，u′′、v′′和湍動(dòng)能的分布區(qū)域逐漸收窄，此時(shí)湍動(dòng)能數(shù)值較大的區(qū)域分布在UFC附近。

圖8 t=3.5，θ=0°,45°,90°,135°時(shí)流場(chǎng)RMS速度u′′、v′′以及湍動(dòng)能k在不同 λ下沿x方向的分布Fig.8 Distributions of the fluid RMS velocity u′′,v′′ and turbulent kinetic energy k in the x-direction at different λ,when θ is equal to 0°,45°,90°,135° at dimensionless time t=3.5

為了探究更詳細(xì)的規(guī)律，需要對(duì)流場(chǎng)內(nèi)能、流場(chǎng)動(dòng)能和流場(chǎng)湍動(dòng)能進(jìn)行定量分析，將計(jì)算域x∈[-3.0,4.0]分為3部分，分別是計(jì)算域上游區(qū)域（x∈[-3.0,-0.5]）、橢圓柱云區(qū)域（x∈[-0.5,0.5]）和計(jì)算域下游區(qū)域（x∈[0.5,4.0]），在每個(gè)計(jì)算區(qū)域內(nèi)對(duì)內(nèi)能、動(dòng)能和湍動(dòng)能進(jìn)行積分，如圖9所示，圖中標(biāo)出了每個(gè)區(qū)域具體積分值以便分析。

從內(nèi)能分布來(lái)看，對(duì)于不同的λ和θ，在上游區(qū)域中的內(nèi)能占整個(gè)流場(chǎng)中內(nèi)能的比例最大。隨著λ的增大，主要差異體現(xiàn)在：當(dāng)θ=0°時(shí)，上游區(qū)域內(nèi)能從23.8減小到20.6；而當(dāng)θ=90°時(shí)，變化趨勢(shì)則相反，內(nèi)能從29.4增加到32.9，同時(shí)內(nèi)能在此區(qū)域中的值比其他角度在此區(qū)域中的值都大，在橢圓柱云區(qū)域，內(nèi)能從4.2減小到1.9，在下游區(qū)域，內(nèi)能從7.8減小到6.8；當(dāng)θ=45°，135°時(shí)，在不同λ下，每個(gè)區(qū)域的內(nèi)能對(duì)λ的變化不敏感。

從動(dòng)能分布來(lái)看，當(dāng)θ=0°時(shí)，隨著λ的增大，動(dòng)能在上游區(qū)域從1.18增加到1.54，在橢圓柱云區(qū)域從0.10增加到0.23，在下游區(qū)域從0.83增加到1.53；當(dāng)θ=90°時(shí)，趨勢(shì)則相反，隨著λ的增大，動(dòng)能在上游區(qū)域、橢圓柱云區(qū)域和下游區(qū)域都逐漸減??；當(dāng)θ=45°，135°時(shí)，隨著λ的增大，動(dòng)能在3個(gè)區(qū)域的值緩慢增加。

從湍動(dòng)能分布來(lái)看，當(dāng)θ處于不同角度下，湍動(dòng)能在上游的值均很小，主要分布在橢圓柱云區(qū)域和下游區(qū)域。當(dāng)θ=0°,45°,135°時(shí)，隨著λ的增加，湍動(dòng)能在橢圓柱云區(qū)域和下游區(qū)域的值都逐漸減小；當(dāng)θ=90°時(shí)，湍動(dòng)能主要分布在橢圓柱云區(qū)域。

圖9 t=3.5時(shí)不同θ和λ下流場(chǎng)內(nèi)能、流場(chǎng)動(dòng)能和流場(chǎng)湍動(dòng)能在計(jì)算域上游區(qū)域（x∈[-3.0,-0.5]）、橢圓柱云區(qū)域（x∈[-0.5,0.5]）和計(jì)算域下游區(qū)域（x ∈[0.5,4.0]）分布Fig.9 Distributions of the fluid internal energy,fluid kinetic energy and fluid turbulent kinetic energy at different θ and λ in three different regions,that is the upstream area of the domain x ∈[-3.0,-0.5],elliptical column cloud area x∈[-0.5,0.5],the downstream area of the domain x ∈[0.5,4.0]at dimensionless time t=3.5

圖10 不同角度θ和不同λ下流場(chǎng)內(nèi)能、流場(chǎng)動(dòng)能和流場(chǎng)湍動(dòng)能隨無(wú)量綱時(shí)間t的變化Fig.10 Variations of the fluid internal energy,fluid kinetic energy and fluid turbulent kinetic energy with dimensionless time t at different θ and λ

圖10展示了流場(chǎng)總內(nèi)能、流場(chǎng)總動(dòng)能和流場(chǎng)湍動(dòng)能在t∈[1.2,3.5]的演化過(guò)程。因?yàn)閠∈[0,1.2]時(shí)，流場(chǎng)中的激波未與橢圓柱云相互作用，因此該過(guò)程中的能量變化規(guī)律并不重要，可以忽略。這里主要討論當(dāng)x∈[1.2,3.5]時(shí)3個(gè)量的變化。3個(gè)量中流場(chǎng)總內(nèi)能總是隨著時(shí)間的增加而逐漸增大。當(dāng)θ=90°時(shí)，流場(chǎng)總內(nèi)能增長(zhǎng)得最快；λ越大，θ越小，流場(chǎng)中總內(nèi)能增長(zhǎng)得越慢；到θ=0°時(shí)，流場(chǎng)總內(nèi)能增長(zhǎng)得最慢。總動(dòng)能方面，當(dāng)θ=0°、λ=3或者θ=0°、λ=4時(shí)，計(jì)算域中的總動(dòng)能隨時(shí)間的增加而逐漸增大，而其他情況下則隨時(shí)間的增加而逐漸減小。當(dāng)θ=90°時(shí)，計(jì)算域中的總動(dòng)能減小得較快。就湍動(dòng)能而言，總體上都隨著時(shí)間的增加而逐漸增大，而且在相同θ下，λ越小，湍動(dòng)能越大。綜合來(lái)看，在相同λ的條件下，當(dāng)θ=90°時(shí)，入射激波與橢圓柱云正面沖擊的有效截面積最大，整個(gè)流場(chǎng)中總內(nèi)能增長(zhǎng)得最快，湍動(dòng)能增長(zhǎng)得最慢，而總動(dòng)能減少得最快。

2.3 一維體積平均模型

本節(jié)討論改進(jìn)Jiang等[25]的一維體積平均模型與DNS的結(jié)果對(duì)比。改進(jìn)后的一維體積平均模型的具體形式如下

式中：Cd為人工有效阻力系數(shù)，αd為固體體積分?jǐn)?shù)。表3給出了在t=3.5時(shí)不同λ和θ下最優(yōu)Cd的取值，其中θ=90°時(shí)以整體擬合趨勢(shì)為準(zhǔn)，其他角度則以一維體積平均模型的反射激波和透射激波位置與DNS結(jié)果的激波位置擬合最優(yōu)為準(zhǔn)。圖11展示了λ=2，3，4，且θ=0°，45°，90°和135°時(shí)的速度分布。從圖11中可以看出，一維體積平均模型與當(dāng)前DNS結(jié)果的擬合效果較好。

圖11 t=3.5時(shí)不同的λ和θ下一維體積平均模型與DNS擬合結(jié)果Fig.11 Fitting results of the 1-D volume-averaged model and DNS at different λ and θ when the dimensionless time is equal to 3.5

表3 人工有效阻力系數(shù)Cd的最優(yōu)取值Table 3 Optimal values of artificial effective drag coefficient Cd

利用表3中的數(shù)據(jù)得到如圖12所示的Cd分布。由圖12可以看出：θ∈[0°,45°]時(shí)，λ越大，Cd越?。划?dāng)θ∈[60°,120°]時(shí)，λ越大，Cd越大；當(dāng)θ=90°時(shí)，Cd比相同λ、其他角度下的Cd都要大。在目前擬合的θ和λ范圍內(nèi)，Cd最小值出現(xiàn)在圖12中A處，最大值出現(xiàn)在圖12中B處。

圖12 人工有效阻力系數(shù)Cd的最優(yōu)取值分布Fig.12 Distribution for the optimal value of artificial effective drag coefficient Cd

3 結(jié) 論

研究了不同橢圓柱橫截面長(zhǎng)短軸之比λ和不同橢圓柱橫截面長(zhǎng)軸與來(lái)流方向所成角度θ對(duì)流場(chǎng)的影響，得到如下結(jié)論。

（1）在t=3.5時(shí)，保持λ不變，θ從0°增大到90°，激波與橢圓柱云正面沖擊的有效橫截面積逐漸增大，反射激波的位置逐漸遠(yuǎn)離橢圓柱云，透射激波的位置逐漸靠近橢圓柱云，此時(shí)流場(chǎng)中的能量更多以內(nèi)能形式保留在從反射激波面到橢圓柱云上邊界之間的區(qū)域中。

（2）在t=3.5時(shí)，保持θ不變，當(dāng)θ=0°時(shí)，隨著λ從2增大到4，反射激波的位置逐漸靠近橢圓柱云，透射激波的位置逐漸遠(yuǎn)離橢圓柱云，流場(chǎng)下游區(qū)域的內(nèi)能占流場(chǎng)總內(nèi)能的比例逐漸增加。θ=90°時(shí)，λ從2增加到4，此時(shí)流場(chǎng)上游區(qū)域中的內(nèi)能占流場(chǎng)總內(nèi)能的比例逐漸增加。θ=45°的內(nèi)能分布和動(dòng)能分布對(duì)λ的變化不敏感。

（3）當(dāng)t=3.5、θ=90°時(shí)，隨著λ的增加，u′′、v′′和湍動(dòng)能的分布區(qū)域逐漸收窄，此時(shí)湍動(dòng)能數(shù)值較大的區(qū)域分布在橢圓柱云上邊界附近。而其他角度下，u′′、v′′和湍動(dòng)能的分布對(duì)λ的變化不敏感，且u′′、v′′和湍動(dòng)能數(shù)值較大區(qū)域分布在橢圓柱云下邊界附近。

（4）當(dāng)t∈[1.2,3.5]時(shí)，流場(chǎng)總內(nèi)能隨時(shí)間逐漸增大，相同λ下，當(dāng)θ=90°時(shí)，流場(chǎng)總內(nèi)能增長(zhǎng)最快，θ=0°時(shí)，流場(chǎng)總內(nèi)能增長(zhǎng)最慢。當(dāng)θ=0°、λ=3或者θ=0°、λ=4時(shí)，計(jì)算域中的總動(dòng)能隨時(shí)間的增加而逐漸增大，而其他情況則隨時(shí)間的增加而逐漸減小。當(dāng)θ=90°時(shí)，計(jì)算域中的總動(dòng)能減小較快。湍動(dòng)能隨時(shí)間的增加而逐漸增大，相同θ下，λ越小，湍動(dòng)能越大。在相同λ的條件下，當(dāng)θ=90°時(shí)，動(dòng)能轉(zhuǎn)換成內(nèi)能的效率最高。

（5）一維體積平均模型在采用合適的有效人工阻力系數(shù)Cd時(shí)，可以較好地?cái)M合當(dāng)前DNS結(jié)果。