激波加載單雙圓柱非穩態曳力的數值研究

陳婉君, 章利特, 施紅輝, 郝李娜, 黃保乾

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

激波加載單雙圓柱非穩態曳力的數值研究

陳婉君, 章利特, 施紅輝, 郝李娜, 黃保乾

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

采用Fluent軟件對激波誘導單雙圓柱模型繞流場進行二維數值計算,深入研究馬赫數為1.14的激波與直徑為40 mm的單、雙圓柱模型相互作用時繞柱流場的非穩態曳力的形成機理。結果表明:單、雙柱模型(H=1.5)的曳力系數Cd曲線都存在明顯的波峰和多個波谷結構,并以波幅不斷減小的方式逐漸趨于某個穩定正值,其中波峰峰值遠大于穩態值;相比于單柱情形,由于雙柱空間位置的限制,使得環繞每柱的激波結構可以發生相互干涉,從而影響了環繞圓柱的壓力以及剪切應力分布,導致單雙柱曳力系數曲線波動變化的差異。

激波; 單/雙柱; 非穩態曳力; 曳力系數; 相互干涉

0 引 言

激波加載固體顆粒群的問題是超音速氣固兩相流中的重要研究課題之一,掌握激波與顆粒相互作用機理對工業生產和生活至關重要,醫學界開展的無針注射[1-2]和工業上的超音速冷噴涂[3]都是基于此機理的實際應用。國內外學者針對單顆粒相互作用進行了大量的研究和數值研究,Saito等[4]對氣體—顆粒混合物中非穩定曳力對激波后流場非平衡流場結構的影響開展數值研究,發現較低的激波馬赫數和較高的顆粒質量比對阻力系數時間依賴效應的影響更大;Dan等[5]對激波與單顆粒相互作用的數值研究表明,模型球所受的阻力主要由作用于其表面的流體粘性剪切力和壓力共同決定,馬赫數和顆粒直徑的減小都導致顆粒阻力系數的增大;Tanno等[6]利用安裝傳感器直接測量激波作用在小球上的加速度,獲得與數值模擬結果相當吻合的激波加載模型球的非穩態曳力,同時應用高速攝影儀和雙曝光全系干涉儀共同揭示了在壓力測量和阻力測量中個別波的相互作用。上述研究一般都只著重于激波與單球的相互作用,未考慮在多個模型球的情況下,激波結構會發生相互干涉,以至于會出現區別于單球情形的流場結構。

由于圓柱與圓球在對稱截面上形狀一致,所受阻力(曳力)與壁面壓力、剪切應力關系類似,考慮到本文采用Fluent軟件進行流場數值模擬時,三維計算的網格數目和計算時間超出了承受范圍,而二維圓柱繞流計算并不妨礙激波結構干涉對非穩態曳力影響機制的揭示,因此本文的數值模擬基于二維進行。本文利用Fluent軟件對單雙柱繞流場進行模擬,獲得激波與單雙柱相互作用時的曳力系數曲線以及圓柱周圍流場的壓強,圓柱壁面壓力和剪切應力分布等物理特性,以圓柱所受的壓力和剪切應力為著重點深入分析繞柱流場的變化對曳力系數曲線波動的影響機制,以及單柱與雙柱曳力系數Cd曲線差異的產生原因。

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

1.1 控制方程

本文數值模擬的對象是激波管裝置內激波誘導的單、雙柱繞流場,屬于非定常可壓縮粘性流動。該流動可用笛卡爾坐標系下的N-S方程描述:

連續方程:

(1)

動量方程:

(2)

能量方程:

(3)

其中,

(4)

這里τij為粘性應力張量，ρ為密度,ui為速度分量,P為壓力,E為單位質量流體的內能,κ為熱傳導系數,T為溫度,μ為動力粘性,μT為渦旋粘性系數。

波前氣流和波后氣流均可視為理想流體,以柱直徑和波后氣流為特征參數定義的流動雷諾數為209300,為湍流流動,故采用理想氣體狀態方程和可實現k-ε湍流模型來封閉控制方程。

1.2 二維計算模型及計算參數

二維計算模型如圖1所示,利用Gambit建立單柱的2D數值計算模型,選取1.2 m×0.2 m的計算區域,模型柱直徑為0.04 m,采用結構化網格劃分求解區域,計算網格數設為238793和237278。流動通量類型選擇Roe-FDS,通過patch實現波后氣流的壓力、速度和溫度。表1顯示了依據激波管理論確定的單、雙柱繞流計算參數,其中,無量綱間距H定義為雙柱的柱心間距與柱直徑之比,本文雙柱模型取H=1.5,即雙柱的中心距離為0.06 m。聲速取22℃溫度下的對應值約為344.37 m/s。

圖1 單柱模型計算區域

激波速度VS/(m/s)激波馬赫數MS波后氣流速度v2/(m/s)波后氣流壓力P2/Pa波后氣流溫度T2/K波后氣流密度ρ2/(kg/m3)3941．1477136741321．791．48

1.3 有效阻力系數計算

文獻[5]對激波與單顆粒相互作用的數值研究表明,模型球所受的曳力主要由作用于其表面的流體粘性剪切力和壓力共同決定。壁面壓力與剪切應力跟曳力離散關系如圖2所示,通過對圓柱壁面進行受力分析,可求得非穩態力矢量FD與壁面壓力和剪切應力的積分關系,可由式(5)表示:

FD=∮(-p)cos(n,τ)dζ+∮(σ)csc(τ,i)dζ

(5)

圖2 壁面壓力與剪切應力跟曳力離散關系

本文研究的曳力指該力矢量在激波管軸向的分量FD,數值計算時,需對積分關系式(5)的右端各項離散得到關系式(6),并沿x軸投影來確定FD,再依據曳力系數公式確定Cd。

(6)

(7)

其中F為總的力矢量,由兩個環面積分的矢量和決定;P是壁面壓力;n是外法向單位矢量;S是柱面面積;σ是壁面剪切應力;τ是切向單位矢量(決定剪切應力在切面上的指向);i是x軸方向單位矢量;FD為曳力;Fx為矢量F在x軸上的投影;兩個余弦函數為方向余弦,分別是外法向單位矢量與x軸正向的余弦值,以及切向單位矢量與x軸正向的余弦值。

曳力系數公式為:

(8)

其中,ρ2和v2分別為激波后氣流密度和速度,π為圓周率,r為圓柱半徑。

1.4 流場計算準確性驗證

圖3和圖4分別顯示了不同時刻沿激波管軸線的壓力、相同時刻的速度分布和密度分布。根據圖3計算可知激波速度約為400 m/s,對應激波馬赫數為1.16。波后氣流的速度、壓力和密度分別為75.3～77.6 m/s、136 704～137 762 Pa和1.48～1.49 kg/m3。由于管壁和流體粘性的影響,波后氣流速度、壓力和密度等參數并非均勻值,而是分布在上述范圍內。與表1激波管理論所確定的參數相比,激波速度或者激波馬赫數的相對誤差僅為1.5%,且激波管理論值都包含在上述的分布范圍內。綜上所述,可以驗證本文數值模擬計算模型和方法的準確性。

圖3 不同時刻激波管軸線壓力分布

圖4 t=2.5×10-4 s時刻激波管軸線速度分布和密度分布

2 單雙柱模型數值計算結果及分析

圖5為直徑D=40 mm、Ms=1.14時的單雙圓柱曳力系數Cd與無量綱時間t′的關系。表2為單雙柱模型曳力系數曲線的波峰值和波谷值。從圖5中可以發現單、雙柱模型非穩態曳力系數曲線的趨勢基本一致,隨著無量綱時間t′(定義為圓柱直徑與激波速度之比,下文簡稱時間)的增大,曳力系數曲線先急劇上升,達到第一波峰后急劇下降,再先后出現第一負值波谷、第二波峰、第二負值波谷,繼而波幅逐漸減小,逐漸趨向相對于峰值很小的穩定正值。通過圖5和表2可以發現兩者曳力系數曲線存在一些顯著差異,a)雙柱(H=1.5)最大峰值Cd高于單柱相應值,且雙柱的第一波峰出現的時刻t′=3/4較單柱時的t′=1/2要晚;b)雙柱第二波峰的峰值比單柱相應值要高;c)雙柱曳力系數曲線波動變化持續時間更長,更難以趨于穩態正值。

圖5 單、雙圓柱曳力系數與無量綱時間的關系

類型激波馬赫數Ms第一波峰位置t′峰值Cd第一波谷位置t′谷值Cd第二波峰位置t′峰值Cd第二波谷位置t′谷值Cd單柱1．140．521．931．65-1．22．371．073．07-0．3雙柱1．140．7523．042．51-1．33．241．333．950．006

圖6為直徑D=40 mm、激波馬赫數Ms=1.14條件下,無量綱時間t′=0～2期間單柱周圍的壓力分布。圖7為相應時刻的圓柱壁面壓力和剪切應力分布。由于上、下參數分布幾乎完全重疊,因此,都只給出上半柱模擬結果。為了更好地解釋非穩態曳力系數曲線形成機理,在此作出以下四點說明:其一,t′=0是根據曳力系數上升沿起點確定的,由于這個時間很短,所以很可能會有點偏差;其二,由于激波馬赫數較小,激波達到圓柱壁面前由于柱面影響,靠近軸線的激波陣面中心區域向上游方向發生一定程度彎曲,所以激波陣面中心區到達前駐點時,其前緣已越過前駐點一段距離;其三,由于Fluent在求解粘性流動時,對激波陣面的分辨能力有限,所以激波陣面呈現圖6所示的壓力條帶狀連續變化的情形。其四,如圖7(b)所示,在t′=0～2時間范圍內,由于壁面的剪切應力數值非常小,在激波與圓柱相互作用時間這一很小的時間量級里,剪切應力只對曳力系數曲線的細節變化會有影響,不能改變基本走勢,對曳力僅產生次要影響,可忽略不計,但隨著時間的增長,進入穩態階段后,最終會占據統治地位,使曳力處于穩定的狀態。

圖6(a)對應于t′=0時刻,發現此時的激波陣面已經越過前駐點些許距離,其下游區域依然保持初始壓力,而在前駐點出現非常集中的高壓區。此時壁面的壓力分布情形為:從前駐點往后壓力持續減小,降到80°(注:0°、90°和180°分別對應于前駐點、赤道和后駐點)以后壓力幾乎沒有受到激波任何影響(見圖7(a)的t′=0時刻壁面壓力分布曲線)。由于此時赤道左半部分的壁面壓力對曳力的貢獻全部為正,導致曳力由0急劇增大。t′=1/4時刻,激波陣面分布在45～90°之間,前駐點附近的高壓區域與壁面的接觸范圍進一步擴展,且壁面壓力值相對于t′=0時刻顯著增大,即保持了曳力系數曲線大斜率上升的趨勢。

圖6 D=40 mm、Ms=1.14時,單柱周圍的壓力分布(激波從左往右傳播)

圖7 不同時刻單柱壁面壓力和剪切應力分布

圖6(c)對應于t′=1/2時刻,激波陣面分布在赤道兩側約70～110°范圍,前駐點附近的高壓區域繼續擴大,此時壓力影響范圍擴展到約120°位置(見圖7(a)),在0～90°范圍內的各點處壓力相對于t′=1/4時刻都顯著升高,且在0～30°范圍分布的壓力值相對于其他任意時刻都高。觀察圖6(d)對應于t′=3/4時刻可以發現其壓力作用范圍已超過150°,且90～150°分布的壓力比t′=1/2時刻在各點分布的壓力都高,由曳力積分關系式(1)可知,分布在0～45°范圍的壓力貢獻為正,且較45～90°范圍的壓力影響更大,而在90～180°的貢獻為負。因此,在t′=1/2時刻附近出現了曳力(系數)的最大值,對應了圖3中的曳力系數曲線波峰位置,而且此后到t′=2之間的各時刻柱面在90°～180°范圍的壓力分布對曳力呈負值貢獻,所以曳力系數曲線呈現持續下降變化。

在t′=1/2之后,由于從前駐點側反射的激波向外傳播,影響范圍不斷向外蔓延,而赤道左壁面壓力卻呈現下降趨勢(見圖7(a)),這對曳力在t′=1/2之后的下降具有重要影響。在t′=1/2～1期間,激波陣面壓力分布范圍逐漸覆蓋赤道右半區,但右半區的壓力卻沒有顯著升高(見圖7(a)),這是由于激波在跨越赤道后的沿壁面發生衍射,激波強度隨角度的增大不斷下降。在t′=1到t′=5/4之間的某時刻開始,衍射波在后駐點附近聚集,使局部壓力急劇升高,影響范圍也逐漸擴展,由于從赤道后不同位置發出的衍射波強度和速度不盡相同,因此到達后駐點形成的聚焦行為可能持續出現。又由于后駐點附近的急劇升高的聚焦壓力對曳力的貢獻為負,所以后駐點激波聚焦發生時,曳力系數曲線出現波谷,甚至是負值,而且可能多次出現(見圖5)。

圖8為D=40 mm、Ms=1.14、H=1.5時雙柱周圍的壓力分布,圖9為雙柱上、下側壁面的壓力分布。t′=0時刻,雙柱的上、下兩側激波陣面基本對稱,與單柱情形一致。但從t′=1/4時刻開始,由于雙柱相對空間位置對彼此激波結構的限制作用,雙柱的壓力分布表現出與單柱的差異,激波陣面在雙柱之間出現扭曲,每柱上、下側的壓力分布不對稱性隨時間的延長也逐漸顯現。

圖8 D=40 mm、Ms=1.144時,雙柱周圍的壓力分布(激波從左往右傳播)

由圖8(c)對應的t′=1/2和之后各時刻壓力分布可以發現,入射激波被前駐點附近壁面反射后,以各自前駐點為中心向外輻射發展,各自影響范圍逐漸擴大,這與單柱情形相同,并且在t′=1/2至t′=3/4之中某時刻起雙柱的反射激波在雙柱中心處相遇,開始發生相互干涉。在t′=1時刻,雙柱對稱中心線與雙柱前駐點連線的焦點附近出現高壓區域,稱之為干涉區。在t′=5/4時刻,各柱的反射激波影響范圍和干涉區都顯著擴展,干涉區內的壓力值進一步提升。從t′=7/4時刻開始,盡管反射激波的影響范圍繼續擴大,但壓力值呈現明顯下降,且維持高壓力值,這是由于從雙柱內側壁面不斷反射的激波向上游方向和對方柱面傳播,導致其強度和傳播方向都會有相應的改變,從而引起高壓干涉區的持續存在和發展。此外,干涉區的外層內也保持較高壓力水平,在靠近柱面部分形成尖角結構,并朝對側柱面不斷發展,相信在t′>2之后較短時間內,終將到達柱面,勢必引起赤道左半區內側柱面壓力的顯著升高。除了反射激波的干涉外,從赤道右半區內側壁面不同位置處發出的衍射波同樣也會發生干涉,使得雙柱間激波陣面前緣形狀從凹變凸,勢必改變衍射波的強度、速度和衍射波干涉區域內的壓力分布。

圖9 不同時刻雙柱壁面壓力分布

雙柱(H=1.5)第一波峰出現的原因與圖3的單柱情形基本相同,針對單柱和雙柱曳力系數曲線的差異,其根本原因是雙柱沿激波傳播垂直方向彼此之間空間位置的限制,使得隸屬于每柱的激波結構可以發生相互干涉,從而影響了環繞每柱的壓力以及柱壁面的壓力和剪切應力分布。與單柱類似,對非穩態曳力而言,柱壁面的壓力分布起關鍵作用,而壁面剪切應力只起次要影響。因此,從單、雙柱壁面壓力分布的差異和原因闡述造成上述3點曲線差異的具體原因。

首先,在t′=1/2時刻,雙柱彼此空間位置限制使其間的激波陣面扭曲,引起雙柱赤道左半區內側壁面反射的激波發生相互干涉,導致t′=1/4～2期間赤道左半區內側壁面分布的壓力值相對于單柱的各對應時刻相同位置有所提升,同時由于激波被前駐點附近壁面反射而提升局部壓力的疊加影響,這就使得t′=3/4時刻曳力(系數)在t′=1/4時刻基礎上進一步提升,導致雙柱曳力系數曲線第一波峰出現的時刻較晚,但峰值更大。

其次,雙柱第二波峰的峰值更大,可解釋為單柱的第二波峰是衍射波在后駐點的持續聚焦過后,后駐點近區壓力降低造成的,視為一種被動形式的曳力升高;而雙柱的第二波峰則是反射波經干涉后到達赤道左半區內壁面,造成該區域壁面壓力提升,從而提高了曳力,視為一種主動形式的曳力提升。

最后,雙柱曳力系數曲線波更難以趨于穩定的根本原因是源自雙柱的反射波和衍射波的相互干涉以及到達對側壁面后的再次反射和干涉,如此往復,使得繞柱流場壓力、圓柱壁面壓力和剪切應力波動變化更為復雜,這些參數分布趨于穩態需要經歷更長時間,曳力系數達到穩態值也需要更長時間。

3 結 論

基于Fluent6.3計算平臺對激波誘導單/雙圓柱的繞流場和非穩態曳力進行數值計算,分析激波的反射、衍射和聚焦對非穩態曳力形成的機理。主要結論如下:

a) 在激波馬赫數Ms=1.14條件下,單柱和雙柱模型(H=1.5)的曳力系數Cd曲線的變化規律為出現明顯的波峰和若干的負值波谷,并最終趨于穩定。

b) 在相同條件下,由于雙柱相對空間位置在垂直于激波傳播方向的彼此限制,導致反射激波、衍射波的相互干涉影響了圓柱壁面的壓力和剪切應力的分布,從而出現雙柱的第一波峰和第二波峰出現的時刻t′略晚,第一峰值和第二峰值都偏高以及雙柱曳力系數曲線波動變化持續時間更長,更難以趨于穩態正值等曳力系數曲線間的差異。

c) 對于有運動激波存在的氣固兩相流,不能簡單地采用單顆粒穩態曳力系數模型,當顆粒載荷比較大時,即無量綱間距較小時,需要考慮隸屬于鄰近的不同固體顆粒激波的波系結構的相互干涉對曳力的非穩態影響。

[1] Liu Y. Performance studies of particle acceleration for transdernal drug delivery[J]. Med Comput, 2006, 44: 551-559.

[2] Liu Y. Physical-Msther Mstical modeling of fluid and particle transportation for DNA vaccination[J]. International Journal of Engineering Science, 2006, 44: 1037-1049.

[3] 饒 瓊, 周香林, 張濟山, 等. 超音速噴涂技術及其應用[J]. 熱加工工藝, 2004(10): 49-52.

[4] Saito T, Saba M, Sun M, et al. The effect of an unsteady drag force on the structure of a non-equilibrium region behind a shock wave in a gas-particle mixture[J]. Shock Waves, 2007, 17: 255-262.

[5] Dan I, Ozer I, Lazhar H, et al. Simulation of sphere’s motion induced by shock waves[J]. Journal of Fluids Engineering, 2012, 10(134): 1-4.

[6] Tanno H, Itoh K, Saito T, et al. Shock wave interaction with a sphere in a shock tube[C]//Symposium on Interdisciplinary Shock Wave Research, 2004: 483-497.

(責任編輯： 康 鋒)

Numerical study of Unsteady Drag Force in Shock Wave’s Interaction with Single/Double Cylindrical Models

CHENWan-jun,ZHANGLi-te,SHIHong-hui,HAOLi-na,HUANGBao-qian

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

Two-dimensional numerical calculation of ambient flow field of single/double cylindrical models induced by the shock wave was conducted with Fluent software. Formation mechanism of unsteady drag force of ambient flow field was deeply studied during interactions between shock waves withMs=1.14 and single/double cylindrical models with the diameter of 40 mm. The results show that drag force coefficientCdof single/double cylindrical models (H=1.5) has significant wave crest and multiple troughs and gradually tends to a steady positive value in the form of decreasing amplitude. The wave peak value is much greater than the steady-state value. Compared with single cylinder, due to space position limit of dual cylinders, shock wave structure surrounding each cylinder may mutually interfere, which thus influences distribution of cylindrical pressure and shear stress and leads to differences of fluctuation changes of single/double cylindrical drag force coefficient curve.

shock wave; single/double cylinder; unsteady drag force; drag coefficient; interference

1673- 3851 (2015) 01- 0055- 07

2014-06-05

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

陳婉君(1990-),女,浙江東陽人,碩士研究生,主要從事可壓縮性氣固兩相流方面的研究。

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

TK121

A