高速氣體與橢圓柱云相互作用的數(shù)值研究

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

（1北京計算科學研究中心，北京 100193；2弗吉尼亞大學，美國弗吉尼亞州夏洛茨維爾 22904）

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

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

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

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

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

1 數(shù)值方法

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

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

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

2.1 方法驗證

通過測試網(wǎng)格的收斂性檢驗本方法的正確性。采用4種不同分辨率的網(wǎng)格進行數(shù)值模擬。表1列出了不同分辨率下x和y方向上的網(wǎng)格數(shù)Nx和Ny，以及解析橢圓短軸所使用的網(wǎng)格數(shù)nb，計算區(qū)域設置見圖3(a)。橢圓圓心位于原點 (0 ,0)處，入射激波位于x=-0.04處，計算域為x∈[-0.140，0.140]，y∈[-0.084，0.084]，z∈[-0.020，0.020]。為了節(jié)省計算資源，當x∈[-0.056，0.056]且y∈[-0.028，0.028]時，網(wǎng)格間距是均勻的，且在該區(qū)域Δx=Δy；除此之外，x和y方向采用不等間距的拉伸網(wǎng)格，z方向上的均勻網(wǎng)格個數(shù)設置為2，且該方向上兩側為周期性邊界條件。圖3(a)的上下邊界條件為周期性邊界條件，左邊界為入口邊界條件，右邊界為出口邊界條件，入口邊界和出口邊界都按照氣體的初始條件設置為固定值。氣體初始條件為：波前 (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)格收斂性分析實驗中使用的4種網(wǎng)格Table 1 Four meshes used in the convergence analysis experiment

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

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

2.2 直接數(shù)值模擬

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

圖4 x-y平面計算區(qū)域設置示意圖（右圖為初始橢圓柱云分布圖，藍色表示低壓區(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)格設置Table 2 Mesh settings in numerical simulation of the interaction between plane shock and elliptical column cloud

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

圖5 當λ=2、θ=0°時，不同無量綱時間下流場的無量綱壓強分布Fig.5 Distributions of the dimensionless pressure at different dimensionless time when λ=2 and θ=0°

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

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

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

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

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

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

圖7 t=3.5時不同θ下λ分別為2、3、4時，流場速度、流場內(nèi)能和流場動能的分布（灰色矩形區(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顯示了在無量綱時間t=3.5時，相同θ下，當λ從2增加到4，沿x和y方向上的RMS速度（u′′，v′′）和湍動能的分布，虛線之間的區(qū)域表示橢圓柱云。可以觀察到，當θ為0°、45°和135°時，λ的變化對u′′、v′′和湍動能分布的影響較小，此時u′′、v′′和湍動能數(shù)值較大的區(qū)域分布在DFC附近。當θ=90°時，隨著λ的增加，u′′、v′′和湍動能的分布區(qū)域逐漸收窄，此時湍動能數(shù)值較大的區(qū)域分布在UFC附近。

圖8 t=3.5，θ=0°,45°,90°,135°時流場RMS速度u′′、v′′以及湍動能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

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

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

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

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

圖9 t=3.5時不同θ和λ下流場內(nèi)能、流場動能和流場湍動能在計算域上游區(qū)域（x∈[-3.0,-0.5]）、橢圓柱云區(qū)域（x∈[-0.5,0.5]）和計算域下游區(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 不同角度θ和不同λ下流場內(nèi)能、流場動能和流場湍動能隨無量綱時間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展示了流場總內(nèi)能、流場總動能和流場湍動能在t∈[1.2,3.5]的演化過程。因為t∈[0,1.2]時，流場中的激波未與橢圓柱云相互作用，因此該過程中的能量變化規(guī)律并不重要，可以忽略。這里主要討論當x∈[1.2,3.5]時3個量的變化。3個量中流場總內(nèi)能總是隨著時間的增加而逐漸增大。當θ=90°時，流場總內(nèi)能增長得最快；λ越大，θ越小，流場中總內(nèi)能增長得越慢；到θ=0°時，流場總內(nèi)能增長得最慢。總動能方面，當θ=0°、λ=3或者θ=0°、λ=4時，計算域中的總動能隨時間的增加而逐漸增大，而其他情況下則隨時間的增加而逐漸減小。當θ=90°時，計算域中的總動能減小得較快。就湍動能而言，總體上都隨著時間的增加而逐漸增大，而且在相同θ下，λ越小，湍動能越大。綜合來看，在相同λ的條件下，當θ=90°時，入射激波與橢圓柱云正面沖擊的有效截面積最大，整個流場中總內(nèi)能增長得最快，湍動能增長得最慢，而總動能減少得最快。

2.3 一維體積平均模型

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

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

圖11 t=3.5時不同的λ和θ下一維體積平均模型與DNS擬合結果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°]時，λ越大，Cd越??；當θ∈[60°,120°]時，λ越大，Cd越大；當θ=90°時，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 結 論

研究了不同橢圓柱橫截面長短軸之比λ和不同橢圓柱橫截面長軸與來流方向所成角度θ對流場的影響，得到如下結論。

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

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

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

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

（5）一維體積平均模型在采用合適的有效人工阻力系數(shù)Cd時，可以較好地擬合當前DNS結果。