郝李娜，章利特，王天航，施紅輝

(浙江理工大學機械與自動控制學院，杭州 310018)

?

激波誘導模型球陣非穩(wěn)態(tài)阻力的數(shù)值模擬研究

郝李娜，章利特，王天航，施紅輝

(浙江理工大學機械與自動控制學院，杭州 310018)

采用CFD軟件，對激波誘導單排四球模型球陣繞流場進行三維數(shù)值模擬，在計算模型和方法準確性驗證基礎上，分析模型球非穩(wěn)態(tài)阻力的形成機理，并揭示激波馬赫數(shù)Ms和無量綱間距H對非穩(wěn)態(tài)阻力和激波結構的影響。結果表明：模型球阻力系數(shù)Cd曲線存在一個尖銳的最大值波峰和一個波谷的結構，且最終趨于某一穩(wěn)態(tài)正值；Ms越小，Cd峰值越大，谷值越小，Cd曲線波動幅度越大；激波反射、衍射、聚焦和干涉行為共同影響非穩(wěn)態(tài)阻力系數(shù)；無量綱間距H越小，相鄰球的激波干涉越顯著，Cd峰值越大。

激波；模型球陣；非穩(wěn)態(tài)阻力；阻力系數(shù)；激波干涉

0　引　言

激波與固體顆粒群相互作用的現(xiàn)象普遍存在于自然界和工業(yè)生產(chǎn)中，對該現(xiàn)象的研究一直是超聲速氣-固兩相流領域的重要課題之一。目前，國內外研究學者已針對激波管內可壓縮性氣-固兩相流的現(xiàn)象，采用實驗和數(shù)值模擬等手段進行了一系列研究。Sun等[1]對平面激波加載單模型球時的非穩(wěn)態(tài)阻力進行了實驗測量，發(fā)現(xiàn)在模型球直徑或雷諾數(shù)較大的情況下，非穩(wěn)態(tài)阻力會較大。Saito等[2]通過對含激波的氣體-顆?；旌衔飪上嗔鲃舆M行數(shù)值模擬，研究了不穩(wěn)定阻力對激波后非平衡流場結構的影響，發(fā)現(xiàn)在激波與顆粒群相互作用早期非穩(wěn)態(tài)阻力對流動結構的影響較大。Igra等[3]對激波誘導的單顆粒繞流場進行了數(shù)值模擬后得出，作用在小球上的阻力主要取決于球表面壓力和粘性剪切應力分布以及尾跡湍流結構。施紅輝等對激波與顆粒群相互作用時的反射、透射和衍射行為進行了實驗研究[4-6]，發(fā)現(xiàn)激波的這些動力學行為導致了顆粒群阻力系數(shù)顯著的非穩(wěn)態(tài)性，表現(xiàn)為其具有明顯的波峰、波谷結構。Shi等[7]對激波與模型球陣相互作用時的反射和聚集機理進行了實驗研究，發(fā)現(xiàn)相鄰球之間存在明顯的激波干涉行為，認為顆粒群阻力系數(shù)建模時必須將其加以考慮。陳婉君等[8]對激波誘導的雙圓柱繞流場和阻力系數(shù)進行了二維數(shù)值計算，發(fā)現(xiàn)了明顯的雙圓柱激波結構的干涉行為，這對圓柱面壓力、剪切應力分布并最終對圓柱阻力產(chǎn)生顯著的影響。Jourdan等[9]利用光學方法捕捉了單顆粒球的運動軌跡，并進行了顆粒運動參數(shù)的推導，據(jù)此研究激波驅動下顆粒的運動規(guī)律。Parmar等[10]利用插值方法得到了雷諾數(shù)、馬赫數(shù)適用范圍都很廣且精度相當高的單球阻力系數(shù)經(jīng)驗公式，該公式可直接用于固體顆粒相為稀相時的氣-固兩相流建模。

上述研究主要關注激波與單球的相互作用，一些學者(如Parmar等[10])所建立的阻力系數(shù)模型在稀顆粒相氣-固兩相流中的適用性已經(jīng)得到充分驗證。從先前學者對激波與多球(圓柱)相互作用的研究中不難得出，鄰近球(圓柱)的激波結構存在干涉行為，它勢必會對顆粒阻力產(chǎn)生影響，但這方面的物理機理和規(guī)律研究仍不夠全面深入。鑒于顆粒阻力在兩相流相間耦合建模方面的重要性，有必要開展激波誘導球陣的繞流場及阻力的數(shù)值計算，探索顆粒群非穩(wěn)態(tài)阻力形成機理和影響規(guī)律。為此，本文將介紹激波誘導球陣(單排四球)繞流場和阻力的計算模型與方法，并探討激波馬赫數(shù)Ms、無量綱間距H對非穩(wěn)態(tài)阻力的影響。

1　三維計算模型與計算方法

1.1控制方程

激波誘導模型球陣繞流場屬于三維非定常可壓縮粘性流動，該流動可用Navier-Stokes方程描述。

連續(xù)性方程：

(1)

動量守恒方程：

(2)

能量方程：

=-

(3)

其中E表示單位質量總能，

(4)

(5)

1.2三維計算模型及計算參數(shù)

三維計算模型如圖1所示，計算區(qū)域總尺寸為0.7m×0.2m×0.2m，其中波后氣流區(qū)域尺寸為0.2m×0.2m×0.2m，模型球(陣)中心到計算區(qū)域入口端面距離為0.3m，模型球直徑均為0.04m。計算區(qū)域的截面尺寸與實際測試段一致，但為了節(jié)約計算時間，同時又不影響模型球陣非穩(wěn)態(tài)阻力的分析，球陣前后均只截取的部分長度，計算結果在入口或出口反射波到達球陣前有效。采用壓力進口、壓力出口和無滑移壁面(包括管和球壁面)條件。計算區(qū)域的網(wǎng)格劃分采用了結構化網(wǎng)格，球的近壁面區(qū)域作了網(wǎng)格加密處理，網(wǎng)格總數(shù)約為90萬。表1顯示了根據(jù)激波管理論確定的流場計算參數(shù)，其中：p2、V2、T2、p0和T0分別表示波后氣流的壓力、速度、溫度、滯止壓力和滯止溫度； p1和T1分別表示波前氣流的壓力和溫度，分別取為實驗條件下的平均值0.1013MPa和28 ℃。激波前、后氣體均為空氣。激波馬赫數(shù)Ms定義為入射激波速度與被驅動初始聲速(約為347.85m/s)之比，無量綱間距H定義為兩相鄰球球心間距L(如圖7所示)與球直徑D之比，Ms和H取值范圍均來自于實驗工況。

圖1　單排四球模型球陣繞流場的計算區(qū)域

MsHV2/(m·s-1)ρ2/(kg·m-3)p2/PaT2/Kp0/PaT0/K1.06361.18091.26981.41061.6/1.8/2.035.81.298116840313.7117706314.3596.81.534147963336.1156147341.30139.81.715173717352.9194165364.30203.42.002218331379.9275205405.89

1.3非穩(wěn)態(tài)阻力系數(shù)的計算

Igra等[3]通過對激波與單模型球相互作用的數(shù)值模擬發(fā)現(xiàn)，模型球所受的阻力主要由作用于其表面的流體粘性剪切力和壓力共同決定，故本文進行了模型球表面壓力和剪切應力的分析。圖2中，Δσ為球表面任一局部小面積；γο為外法向單位矢量，其指向與球面壓力Pγ方向相反；φo為球面緯線切向單位矢量，與緯線剪切力τφ同向；θο為球面經(jīng)線切向單位矢量，與經(jīng)線剪切力τθ同向；i是x軸方向單位矢量；j是y軸方向單位矢量；k是z軸方向單位矢量；θ為球心角，即為外法向單位矢量γο與z軸正向之間的夾角；φ為球面緯線切向單位矢量φo與x軸負向之間的夾角。

由于球阻力等于x方向的合外力，即球面壓力Pγ、緯線剪切力τφ和經(jīng)線剪切力τθ對應的力沿x軸投影的疊加，可建立非穩(wěn)態(tài)阻力FD與球面壓力和剪切應力的積分關系：

FD=?

(6)

設球面小面積總數(shù)為N，可得到阻力FD的離散形式：

FD=

(7)

圖2　球面壓力與剪切應力分析圖

(8)

阻力系數(shù)定義為：

(9)

其中：ρ2和v2分別為波后氣流密度和速度，π為圓周率，R為球半徑。

1.4計算準確性驗證

圖3(a)和圖(3)(b)分別顯示了Ms=1.2698時，t=0和50 μs兩個時刻激波管軸線上的壓力和速度分布的計算結果，其中t=0對應了激波恰好到達球前駐點的時刻，激波陣面此時位于x=0.2 m處。根據(jù)圖3(a)和圖3(b)中的壓力和速度間斷位置可知，激波在50 μs內的傳播距離為0.02215 m，故激波速度和馬赫數(shù)的數(shù)值計算結果分別為443 m/s和1.2735。其中，數(shù)值計算的激波馬赫數(shù)與理論值Ms=1.2698相比，相對誤差僅為0.29%。

圖3　沿激波管軸線的壓力分布和速度分布

圖4為Ms=1.2698時的Cd實驗曲線和計算結果的比較，可以發(fā)現(xiàn)兩條波形曲線具有極相似的變化規(guī)律，即Cd隨著時間的增大，先出現(xiàn)一個尖銳的最大正值波峰，繼而快速減小到一個負值波谷，之后上升回復到正值，最終趨于一個穩(wěn)定正值。數(shù)值計算的Cd峰值為5.033，實驗測量值為4.967，相對誤差僅為1.03%。文中提及的實驗Cd曲線和實驗紋影圖均由本文實驗獲得。由于實驗采用的激波管裝置內徑較大，為200 mm裝置，膜片的不規(guī)則破裂引起了尾隨入射激波的小波結構，又由于此時激波馬赫數(shù)較小，小波衍射和聚集長時間存在，這導致了Cd實驗曲線在負值波谷之后出現(xiàn)較明顯的小幅(相對于峰值)振蕩階段。然而在數(shù)值模擬中，采用了正激波間斷關系，相當于采用了理想的規(guī)則破膜假設，即不存在上述小波結構，因此計算的Cd曲線在上述相應階段未出現(xiàn)明顯振蕩。

圖4　Ms=1.2698時，單排四球模型Cd實驗曲線和計算結果的比較

圖5為Ms=1.2698時，不同時刻(圖5 (a)t=20 μs、圖5 (b)t=45 μs、圖5 (c)t=90 μs數(shù)值計算的密度梯度圖與實驗紋影圖的對比。由于實驗紋影圖視角垂直于激波管軸線，四球中距離觀察者最近和最遠的兩球影像重疊地位于上、下兩球的正中間處，因此實驗圖像顯示為并列縱排的三球影像。而數(shù)值計算密度梯度圖為激波管軸截面圖，故其顯示為并列縱排兩球的影像。通過對比可以發(fā)現(xiàn)，數(shù)值計算跟紋影拍攝的透射和反射激波結構在不同時刻都很好地吻合。為定量比較起見，定義數(shù)值計算的密度梯度圖(左半)和實驗紋影圖(右半)中球心到球面反射激波面前緣的距離分別為Xc和Xe(見圖5 (b))，通過測量可以發(fā)現(xiàn)，在t=45 μs時刻，前后兩者之間的相對誤差僅為2.6%。

圖5　Ms=1.2698時單排四球模型密度梯度數(shù)值模擬圖(左半)與實驗紋影圖(右半)的對比

以上分別從流場參數(shù)、模型球阻力系數(shù)和激波結構三方面，進行了數(shù)值模擬結果與實驗數(shù)據(jù)(或激波管理論值)的對比，充分驗證了本文數(shù)值計算模型和方法的準確性。

2　數(shù)值模擬結果及分析

2.1阻力系數(shù)Cd與時間t的關系

圖6顯示了H=1.6時，不同Ms下模型球阻力系數(shù)Cd與時間t的關系。從圖6中可以發(fā)現(xiàn)，不同Ms下阻力系數(shù)曲線的變化趨勢基本一致，即：隨著t的增大，Cd曲線先急劇上升，出現(xiàn)一個最大正值波峰后急速減小，出現(xiàn)一個波谷，繼而漸進趨于一個穩(wěn)定值，嚴格來講，最終趨于一個穩(wěn)態(tài)階段的正值。且Ms越小，Cd峰值越大，谷值越小，波動幅度越大，Cd越難趨于穩(wěn)定。由圖6中最小激波馬赫數(shù)Ms=1.0636時的Cd曲線可以發(fā)現(xiàn)，在1000μs計算時間內，Cd仍未回復正值，但通過更長時間(約2000μs)的計算可以發(fā)現(xiàn)，Cd最終可以回復正值，見圖9。Ms對Cd產(chǎn)生上述影響的原因是，強度較弱的激波在球赤道后的衍射和球后駐點附近區(qū)域的聚集行為會持續(xù)更長時間，從而使Cd的非穩(wěn)態(tài)階段維持的時間更長。

圖6　不同激波馬赫數(shù)下單排四球模型阻力系數(shù)與時間的關系曲線

2.2阻力系數(shù)Cd與壓力云圖的同步比對

為了深入研究不同球之間激波結構相互干涉的影響，取間距較小的相鄰兩球為觀測對象，并選取如圖7所示的球陣斜截面為壓力云圖提取截面。圖8為激波馬赫數(shù)Ms=1.1809時，不同時刻球軸截面外半圓周上的壓力和剪切應力分布，此處外半圓周指相對兩球(見圖7中的上、下兩球)軸截面外側半圓周線，可以看出，在本文關注的激波與球陣相互作用期間，壓力比剪切應力值大4～5個量級，這說明該階段阻力FD取決于球周圍壓力分布，而受球面剪切應力影響極小，完全可以忽略不計，因此本文只進行壓力云圖與阻力系數(shù)Cd曲線的同步對比。

圖7　單排四球模型球陣數(shù)據(jù)提取截面示意

圖8　不同時刻球軸截面外半圓周上的壓力和剪切應力分布

圖9和圖10分別為Ms=1.1809時模型球Cd曲線和對應于圖9中各標注時刻球陣斜截面上的壓力云圖。由圖10可以發(fā)現(xiàn)，激波陣面后的壓力呈條帶狀分布，而并非預期的大片區(qū)域均勻分布，這可能是以下兩方面原因造成的：一是數(shù)值模擬時考慮了流體粘性的影響，二是所采用的二階迎風差分格式對激波陣面的分辨能力有限。通過更細地時間同步對比可發(fā)現(xiàn)：在a時刻，入射激波還未到達模型球前駐點，球面壓力依然保持初始狀態(tài)(見圖10(a))，阻力系數(shù)為0；在b時刻，入射激波開始與模型球相互作用，在小球前駐點附近開始出現(xiàn)高壓區(qū)(見圖10(b))，阻力系數(shù)由0增大到2左右；在c時刻，入射激波到達模型球赤道附近，球面反射激波發(fā)生干涉，兩個球的高壓區(qū)域相交(見圖10(b))，阻力值急劇增大到峰值；在d、e時刻，激波繼續(xù)向下游移動并發(fā)生衍射，赤道下游球面激波投射區(qū)域擴大，同時赤道上游半球的高壓區(qū)域向外擴大、壓力下降(見圖10(d)和(e))，阻力下降；在f時刻，入射激波波前已遠離后駐點，處于圖10壓力云圖視野范圍之外，球面反射激波已遠離赤道上游球面，之前由于球面反射造成的激波后高壓區(qū)已消失，另一方面，由于衍射激波在后駐點附近的持續(xù)聚集，該處變?yōu)楦邏簠^(qū)(見圖10(f))，從而導致阻力最小值的出現(xiàn)；在g時刻，在前駐點附近發(fā)生激波的二次反射，使得該處又出現(xiàn)高壓區(qū)域，此時衍射波在后駐點近區(qū)的聚集行為已經(jīng)結束，該處的高壓區(qū)也已完全消失(見圖10(g))，導致了阻力第二峰值的出現(xiàn)。

圖9　H=1.6，Ms=1.1809時的Cd曲線

圖10　對應于圖9中各標注時刻球陣斜截面上的壓力云圖

2.3球面壓力分布曲線

圖11顯示了不同激波馬赫數(shù)Ms下Cd曲線代表時刻球軸截面外半圓周上的壓力分布，其中0°、90°和180°分別對應了模型球的前駐點、赤道線和后駐點位置。圖11(a)為不同Ms下Cd波峰時刻球軸截面外半圓周上的壓力分布，從圖11中可以看出，在Cd波峰時刻，前駐點處壓力值最大，壓力值隨位置的后移持續(xù)減小，在球后駐點前某處達到最小值，即初始壓力101325 Pa，且Ms越大，球面上分布的壓力值越大。由圖11(b)可以看出，Cd波谷時刻，球面分布的壓力值隨位置的后移，先減小后增大，壓力最小值出現(xiàn)在約100°位置，且Ms越大，壓力變化越劇烈，前、后駐點處的壓力差值越大。由圖11(c)可以看出，Cd第二波峰時刻，前駐點處壓力值最大，但該值顯著小于Cd波峰時刻的相應值，球面壓力隨角度增大，先減后增，在后駐點附近可能再次減小，壓力最小值出現(xiàn)在約80°位置處，且Ms越大，壓力變化越劇烈。

圖11　不同Ms下，Cd曲線代表時刻球軸截面外半圓周上的壓力分布

2.4H對激波干涉和Cd峰值的影響

圖12為激波馬赫數(shù)Ms= 1.2698時，不同無量綱間距H下Cd峰值時刻球軸截面外半圓周上的壓力分布，可以發(fā)現(xiàn)，此時入射激波前緣均到達赤道面附近，無法分辨位置差異，無量綱間距H=1.6和H=1.8時反射激波在雙球內側發(fā)生明顯的干涉行為，表現(xiàn)為雙球反射激波陣面的交匯和波后高壓區(qū)的合并，且H=1.6時上述現(xiàn)象相對更為顯著，其結果是對雙球內側壓力分布乃至球的阻力影響更大，而當H=2.0時，未見明顯的雙球反射激波干涉行為。綜上可得出，相鄰球無量綱間距越小，反射激波干涉越顯著，進而對球阻力的影響越大。

圖12　Ms= 1.2698時，不同H下Cd峰值時刻球軸截面外半圓周上的壓力分布

圖13為Ms=1.2698時，不同無量綱間距下的曳力系數(shù)Cd曲線。發(fā)現(xiàn)隨著H的減小，由于激波間的干涉行為逐漸變得更加劇烈，Cd峰值增大，在出現(xiàn)第二峰值之后，Cd曲線漸漸平穩(wěn)，并趨近于一個穩(wěn)定的正值。

圖13　Ms= 1.2698時不同H下的Cd曲線

3　結　論

本文進行了單排四球球陣在激波作用下的繞流場和非穩(wěn)態(tài)阻力的數(shù)值計算和分析，主要結論如下：

a)模型球阻力系數(shù)Cd曲線先后出現(xiàn)一個尖銳的最大值波峰和最小值波谷，激波馬赫數(shù)Ms較小時可能出現(xiàn)負的Cd谷值，隨后逐漸回復并最終趨于某一穩(wěn)態(tài)正值。

b)Ms越小，Cd峰值越大，谷值越小，Cd曲線波動幅度越大，需要更長時間趨于穩(wěn)態(tài)。

c)激波與球陣相互作用的過程中，出現(xiàn)的激波反射、衍射、聚焦和干涉行為共同影響球面的壓力分布，進而影響非穩(wěn)態(tài)阻力(系數(shù))。

d)無量綱間距H會影響球面反射激波的干涉行為，從而改變球面壓力分布，進而導致Cd曲線的差異；H越小，激波干涉越顯著，Cd峰值越大。

[1]SUNM,SAITOT,TAKAYAMAK,etal.Unsteadydragonaspherebyshockwaveloading[J].ShockWaves,2005,14(1/2):3-9.

[2]SAITOT,SABAM,SUNM,etal.Theeffectofanunsteadydragforceonthestructureofanon-equilibriumregionbehindashockwaveinagas-particlemixture[J].ShockWaves,2007,17(4):255-262.

[3]IGRAD,IGRAO,HOUASL,etal.Simulationofsphere’smotioninducedbyshockwaves[J].JournalofFluidsEngineering,2012,134(10):21-24.

[4]ZHANGLT,SHIHH,WANGC,etal.Aerodynamiccharacteristicsofsolidparticles’accelerationbyshockwaves[J].ShockWaves,2011,21(3):243-252.

[5] 章利特,施紅輝,王超,等.激波與可運動顆粒群相互作用反射與透射機理實驗研究[J].應用力學學報,2010,27(2):280-285.

[6] 張曉娜,岳樹元,章利特,等.激波驅動的氣固兩相流力學特性研究[J].水動力學研究與進展:A輯,2008,23(5):538-545.

[7]SHIHH,YAMAMURAK.Theinteractionbetweenshockwavesandsolidspheresarraysinashocktube[J].ActaMechanicaSinica,2004,20(3):219-227.

[8] 陳婉君,章利特,施紅輝,等.激波加載單雙圓柱非穩(wěn)態(tài)曳力的數(shù)值研究[J].浙江理工大學學報,2015,30(2):208-212.

[9]JOURDANG,HOUASL,IGRAO,etal.Shocktubestudyofthedragcoefficientofasphere[J].ShockWaves,2009,7:521-526.

[10]PARMARM,HASELBACHERA,BALACHANDARS.Improveddragcorrelationforspheresandapplicationtoshock-tubeexperiments[J].AIAAJournal,2010,48(6):1273-1276.

(責任編輯： 康鋒)

Numerical Simulation Study on Unsteady Drag Force Acting on Model-Sphere Array Induced by Shock Wave

HAOLina,ZHANGLite,WANGTianhang,SHIHonghui

(Faculty of Mechanical Engineering & Automation, Zhejiang Sci-Tech University, Hangzhou 310018, China)

In this study, the 3D numerical simulation of a gas flow around a four-spheres-in-one-row array induced by shock wave was performed with CFD software. Based on verifying the accuracy of calculation model and method, the formation mechanism of unsteady drag force exerted on a model sphere and the effect of the Mach numberMsof shock wave and the non-dimensional interval distanceHon unsteady drag force and shock wave structure were analyzed. The result shows that each curve of the drag coefficientCdhas a sharp maximum value peak and a valley, and finallyCdtends to be a steady-state positive value. Furthermore, the smaller the value ofMsis, the larger the peak value and the fluctuating amplitude of the curve, and the smaller the valley value ofCd. Besides, the combination of shock wave reflection, diffraction, focusing and interference behaviors affects the unsteady drag coefficient. And the smaller the value ofHis, the more remarkable the shock interference between adjacent spheres, and the larger the peak value ofCd.

shock wave; model-sphere array; unsteady drag force; drag coefficient; shock wave interference

10.3969/j.issn.1673-3851.2016.09.0016

2015-11-16

國家自然科學基金項目(51006091)；浙江省自然科學基金項目(LY13E060011)；流體機械及工程省重點學科及流體工程技術創(chuàng)新團隊項目(11130031201301)；流動腐蝕與防控技術創(chuàng)新團隊(浙理工科〔2013〕13號)

郝李娜(1990-)，女，山西太原人，碩士研究生，主要從事湍流與復雜流動以及可壓縮性與瞬態(tài)流動方面的研究。

章利特，E-mail: langzichsh@zstu.edu.cn

TK121

A

1673- 3851 (2016) 05- 0726- 08 引用頁碼： 090403