汪東旭 張尊華 梁俊杰 衡怡君 李格升
(武漢理工大學(xué)能源與動力工程學(xué)院 武漢 430063)
激波管是一種利用激波實現(xiàn)對工質(zhì)進行加熱獲得目標(biāo)范圍內(nèi)溫度和壓力的實驗裝置[1-2],其包括兩端封閉的圓管和膜片分隔開的驅(qū)動段和被驅(qū)動段.激波管有結(jié)構(gòu)簡單,可在較寬范圍調(diào)節(jié)試驗區(qū)的溫度和壓力的優(yōu)點,廣泛應(yīng)用于研究可壓縮流體現(xiàn)象、高溫氣體、氣相燃燒反應(yīng)和超音速燃燒等領(lǐng)域.
隨著國內(nèi)外流體數(shù)值模擬軟件的發(fā)展,計算流體力學(xué)(CFD)數(shù)值計算技術(shù)在激波與介質(zhì)相互作用機理研究中占據(jù)較為重要的地位,可應(yīng)用計算流體動力學(xué)流場可視化方式研究和解釋激波管內(nèi)的流場結(jié)構(gòu)[3].激波管在設(shè)計優(yōu)化前,需預(yù)先了解激波管中的流動和激波運行過程,由于數(shù)值模擬可以在較廣的參數(shù)范圍內(nèi)較快給出定量效果,所以可在設(shè)計初期利用數(shù)值模擬來優(yōu)化激波管性能參數(shù),選擇最佳結(jié)構(gòu)和驅(qū)動方案;使用CFD數(shù)值模擬方法可設(shè)置多測點監(jiān)測激波的傳播及參數(shù)變化曲線,能夠直觀預(yù)測激波運行過程中內(nèi)部氣流之間、流體與管壁之間和流動隨時間變化過程中氣體間的相互作用而產(chǎn)生的熱力學(xué)特性,獲得空氣動力參數(shù),如溫度、壓力、密度等參數(shù)分布;利用數(shù)值方法對激波管內(nèi)流場流動進行計算對實際設(shè)備的傳感器測量過程也有參考意義.
近年來,國內(nèi)外學(xué)者針對激波管內(nèi)激波和流體運動特性進行了一系列試驗與模擬研究,無黏模型和層流黏性模型在激波管的流動數(shù)值計算中得到廣泛應(yīng)用.Lamnaouer等[4]針對高壓激波管,模擬了激波在非反應(yīng)流和反應(yīng)流中的傳播和反射,對激波管側(cè)壁傳熱和反射激波/邊界層相互作用等非理想狀態(tài)進行了量化研究,模型能夠準(zhǔn)確地模擬激波和膨脹波的傳播和反射以及激波反射后的流動不均勻性;Amrollah等[5]在兩種不同驅(qū)動氣體壓力和膜片厚度的情況下,基于激波管二維軸對稱模型研究了正激波在層流粘性流體中的運行過程,仿真激波通過的時間間隔數(shù)據(jù)與試驗數(shù)據(jù)一致性良好;Davidson等[6]利用激光診斷技術(shù)對激波管邊界層效應(yīng)進行了試驗研究,結(jié)果表明:邊界層與激波分岔發(fā)生在湍流區(qū)域,它們受激波管結(jié)構(gòu)和激波前后條件影響;何中偉等[7]通過試驗研究發(fā)現(xiàn)實際激波管中管壁壁面邊界層基本上是湍流邊界層占主導(dǎo),激波與邊界層干擾中,大多是激波與湍流邊界層的相互干擾.為更準(zhǔn)確反映激波管中的物理現(xiàn)象,有必要對激波管進行湍流流動模擬研究并驗證其結(jié)果的可信度,其結(jié)果可為下一步開展激波管流體力學(xué)與化學(xué)動力學(xué)耦合的參數(shù)化研究提供重要理論依據(jù).
利用數(shù)值模擬的方法,模擬理想氣體空氣在總長為2.08 m、直徑為67.00 mm的激波管內(nèi)湍流粘性流場的流動過程,并與相關(guān)試驗數(shù)據(jù)進行對比,驗證模擬方法的可靠性.依據(jù)不同時刻下激波傳播位移和入射激波后溫度、壓力等參數(shù)隨時間的變化,總結(jié)激波運行規(guī)律,獲得激波管最大試驗時間限度.
研究對象為軸對稱結(jié)構(gòu)激波管,采用二維模型對計算區(qū)域進行建模,見圖1.該模型參考了阿米爾卡比爾理工大學(xué)的小尺寸激波管全幾何形狀,激波管長度為2.08 m,直徑為67.00 mm,主要由高壓段、中間膜片、低壓段三個部分組成,破膜方式為壓差破膜,當(dāng)高壓段和低壓段達到一定壓差時膜片破裂產(chǎn)生激波.
圖1 激波管計算區(qū)域
X=0處定義為膜片位置,高壓驅(qū)動段長度為0.64 m,低壓被驅(qū)動段長度為1.44 m,為記錄不同位置的流體壓力隨時間的變化,在高壓段靠近膜片處至低壓段尾端依次分別設(shè)置了9個壓力傳感器監(jiān)測點Sn,監(jiān)測點位置見表1.
表1 監(jiān)測點位置
為描述激波管內(nèi)的流體流動,采用非定常雷諾平均納維-斯托克斯方程進行數(shù)學(xué)建模.其流體力學(xué)控制方程包括連續(xù)性方程、動量方程和能量方程,為
連續(xù)方程:
(1)
動量方程:
(2)
能量方程:
(3)
式中:keff為有效熱傳導(dǎo)系數(shù);Jj為組分j的擴散流量;Sh為其他體積熱源項,不考慮化學(xué)反應(yīng)的情況下,Sh=0,以上方程右邊的前三項分別表示熱傳導(dǎo)、組分?jǐn)U散和黏性耗散引起的能量運輸.
控制方程通過有限體積法轉(zhuǎn)換為可以用數(shù)值方法求解的代數(shù)方程,在每個控制體內(nèi)積分控制方程,進而產(chǎn)生基于控制體的每一個變量都守恒的離散方程.對與任意控制體積V的積分形式的方程有:
(4)
式中:ρ為密度;v為速度矢量;A為曲面面積矢量;Γφφ為φ的擴散系數(shù);φ為φ的梯度,對二維問題即為(?φ/?x)i+(?φ/?y)j;Sφ為單位體積φ的源項.對于每個網(wǎng)格區(qū)域:
(5)
式中:Nfaces為封閉單元面的數(shù)量;φf為通過表面f的對流量;ρfvf·Af為通過表面的質(zhì)量流量;Af為表面f的面積;φf為表面f上φ的梯度;V為單元體積.
湍流是一種高度復(fù)雜的三維非穩(wěn)態(tài)、帶旋轉(zhuǎn)的不規(guī)則流體運動[8].超聲速流動模擬中常用的有Baldwin-Lomax代數(shù)渦黏性模型,Spalart-Allmaras湍流模型,k-ε雙方程模型,為仿真激波管內(nèi)多尺度復(fù)雜流動的激波及非定常流場流動,選用模擬復(fù)雜流動和高壓力梯度流動均有較好性能的雙方程湍流模型,該模型用于計算有旋的均勻剪切流,平面混合流,平面射流,圓形射流,管道內(nèi)充分發(fā)展流動,都取得了與實驗數(shù)據(jù)比較一致的結(jié)果[9].壓力基求解器主要適用于低速不可壓縮流體,故采用適用于高速可壓縮流體的密度基隱式求解器.求解器的時間類型設(shè)定為瞬態(tài),F(xiàn)lux-Type選擇Roe-FDS差分格式,Roe格式可在流動不連續(xù)性附近產(chǎn)生振蕩解,標(biāo)量梯度采用基于格林高斯網(wǎng)格節(jié)點方法,空間離散化采用二階迎風(fēng)格式,為提高計算精度,選擇雙精度并行模式進行數(shù)值模擬.在求解初始化過程中,利用patch命令對高低壓區(qū)分別設(shè)置相應(yīng)的初始壓力和溫度,仿真迭代時間步長為1×10-7s,最大迭代次數(shù)為20.
運用前處理軟件ICEM劃分網(wǎng)格,為了加速計算收斂,采用良好的結(jié)構(gòu)化四邊形網(wǎng)格對計算區(qū)域進行離散化,在壁面附近創(chuàng)建精細(xì)的網(wǎng)格層以捕捉邊界層效應(yīng),局部計算網(wǎng)格見圖2.
圖2 局部計算網(wǎng)格
FLUENT軟件數(shù)值模擬結(jié)果與網(wǎng)格劃分質(zhì)量有很強的關(guān)聯(lián)性,為減少數(shù)值誤差對模擬結(jié)果的影響,進行了網(wǎng)格敏感性研究.在長為2.08 m、直徑為67 mm的模型中,繪制網(wǎng)格數(shù)目分別為700×15,1 000×20,2 000×25,2 500×40和2 999×80的五種計算網(wǎng)格,X軸方向基礎(chǔ)網(wǎng)格尺寸分別約為2.971,2.080,1.040,0.832和0.694 mm,Y方向全局加密比例因子設(shè)置為1.111,以上網(wǎng)格正交質(zhì)量接近1,網(wǎng)格最大長寬比均低于邊界值80,處于理想范圍內(nèi)[10].選用相同的邊界條件(air-air,p4=500 kPa,p1=50 kPa,T4,1=300 K)和數(shù)值方法進行數(shù)值模擬,記錄t=1.0 ms時刻激波管徑向Y=0.03 m位置所對應(yīng)的3區(qū)和2區(qū)接觸面處溫度邊界層內(nèi)等值線軸向坐標(biāo)值,見圖3.由圖3可知:發(fā)現(xiàn)隨著網(wǎng)格數(shù)增大,徑向Y=0.03 m位置所對應(yīng)的等溫線軸向坐標(biāo)值差異越小,當(dāng)網(wǎng)格數(shù)增大到2 000×25時,激波管溫度邊界層內(nèi)溫度分布幾乎不受網(wǎng)格變化影響,綜合考慮計算成本和求解精確性等方面,選擇網(wǎng)格數(shù)為2 500×40(X軸方向基礎(chǔ)網(wǎng)格尺寸為0.832 mm)的網(wǎng)格劃分方案開展二維激波管內(nèi)流體流動數(shù)值模擬研究.
圖3 徑向Y=0.03 m時,激波管內(nèi)溫度等值線軸向坐標(biāo)分布
實際激波管中可觀測到湍流和層流兩種邊界層,層流邊界層向湍流邊界層轉(zhuǎn)變是通過臨界雷諾數(shù)定義的,層流邊界層只有在被驅(qū)動段初始壓力很低、激波較弱等條件下可能存在,邊界層大部分區(qū)域是湍流的.y+值為管壁有界流動的無量綱壁面距離,是評判各種模型壁面邊界層網(wǎng)格分辨率的一個重要參數(shù)[11].采用標(biāo)準(zhǔn)壁面函數(shù)劃分湍流邊界層網(wǎng)格,合理的y+值應(yīng)在(30,300)范圍內(nèi),以保證流場結(jié)構(gòu)解析到粘性子層,壁面邊界層網(wǎng)格y+值計算公式為
(6)
式中:Δy為第一層網(wǎng)格與壁面之間的距離;ρ為流體密度;UT為流體速度;μ為黏性系數(shù),隨著網(wǎng)格分辨率的增加,y+值會減小,利用這種規(guī)律使網(wǎng)格分辨率達到理想值.通過調(diào)整管壁第一層網(wǎng)格高度后經(jīng)計算獲得.由于本文流場是瞬態(tài)流動,需利用不同時刻壁面邊界層y+值判斷邊界層網(wǎng)格設(shè)置質(zhì)量,值設(shè)定為0.05 mm時,記錄Case1激波運行全過程壁面邊界層y+值,不同時刻壁面邊界層y+值均在(30,36)范圍內(nèi),故認(rèn)為本文湍流邊界層網(wǎng)格設(shè)置是合理的,第一層網(wǎng)格與壁面之間的距離值設(shè)置為5×10-5.
由圖2可知,模型共有壁面(wall)、軸對稱(axis-symmetry)、內(nèi)部面(interior)三種邊界類型,忽略浮力和重力影響:①激波管左右兩端端壁和兩側(cè)管壁均設(shè)定為絕熱、無滑移wall邊界;②管體設(shè)定為axis-symmetry邊界;③膜片部分設(shè)定為interior邊界.
數(shù)值模擬初始條件見表2.本文研究背景為無反應(yīng)流場流動,為避免非理想介質(zhì)對模擬結(jié)果精確性的影響,高低壓段氣體工質(zhì)均設(shè)定為理想空氣,黏性常數(shù)為1.789 4×10-5kg/(m·s).基于以下邊界條件,研究激波的傳播規(guī)律,分析激波運行和激波后溫度邊界層的分布規(guī)律.
表2 數(shù)值模擬的邊界條件
為驗證湍流粘性流體數(shù)值模擬方法的適用性,需將數(shù)值計算得到的結(jié)果與實驗數(shù)據(jù)進行對照分析,文獻[5]中激波管幾何模型與本文相同,實驗設(shè)定的邊界條件見表3.本文激波管湍流模型采用相同的實驗條件進行模擬計算,湍流模擬算例記為Case2和Case3,激波位移計算結(jié)果情況見圖4.
表3 文獻實驗條件
圖4 不同工況下入射激波所抵達的位置
由圖4可知,Case2條件下,t=0.002 2 s時,正激波抵達位置的橫軸坐標(biāo)為X1=0.966 m,入射激波速度為vs=452.73 m/s;Case3條件下,t=0.002 0 s時,正激波抵達位置的橫軸坐標(biāo)為X2=0.991 m,入射激波速度為vs=495.5 m/s.
將湍流模擬不同時刻入射激波速度與文獻數(shù)據(jù)進行對比,見表4.由表4可知,本文湍流數(shù)值模擬算例2在0.002 2 s時刻得到的激波速度相對文獻試驗數(shù)據(jù)的誤差為0.291%;算例3在0.002 0 s時刻激波速度相對試驗數(shù)據(jù)的誤差為2.641%.數(shù)值模擬結(jié)果與試驗結(jié)果存在一定的誤差,考慮到試驗過程中激波管管壁為非絕熱壁面,同時實際流體的黏性效應(yīng)和激波管壁面的摩擦效應(yīng)會使激波強度有一定的損失,而模擬時假定激波管管壁為絕熱恒溫?zé)o滑移壁面,因此模擬結(jié)果與試驗結(jié)果必然會產(chǎn)生一定的誤差.湍流模擬結(jié)果相對于層流模擬更接近試驗數(shù)據(jù),且兩個算例模擬結(jié)果產(chǎn)生的誤差數(shù)值小于5%,均在數(shù)值模擬誤差可接受范圍內(nèi),綜上所述,可認(rèn)為該湍流粘性流體數(shù)值模擬方法適合用于預(yù)測激波管內(nèi)激波的運動特性.
表4 湍流模擬結(jié)果與文獻試驗、模擬值對比
在仿真初始狀態(tài),高低壓段充入不同壓力的空氣,此時管內(nèi)包含兩個區(qū)域:初始高壓氣體區(qū)域(4區(qū))、初始低壓氣體區(qū)域(1區(qū)).模擬過程從膜片破裂時刻(t=0 s)開始,假設(shè)X=0 m處的膜片瞬時破裂消失,見圖5.由圖5可知,初始時刻,4區(qū)與1區(qū)之間的壓力間斷面稱為接觸面.
圖5 t=0 s時刻,膜片破裂時激波管內(nèi)瞬時壓力(Pa)云圖(Case4:p4=500 kPa,p1=20 kPa,T4,1=300 K)
高低壓氣體之間的壓差會使氣體進行非定常運動,同時形成正激波向低壓段傳播,入射激波和反射激波對工質(zhì)進行壓縮加熱加壓,由于激波對介質(zhì)的加熱時間極短,可認(rèn)為整個壓縮加熱過程是絕熱、非等熵、均勻加熱的.圖6為t=1.50 ms時刻沿軸向位置管內(nèi)入射激波后管體溫度局部分布云圖.由圖6可知:入射激波壓縮區(qū)域溫度高于初始溫度,膨脹波壓縮區(qū)域溫度低于初始溫度,其中具有溫度梯度的薄層稱為溫度邊界層[12].壁面邊界層溫度高于中心主流區(qū)溫度,這是由于粘性影響,氣體在邊界層內(nèi)減速,動能轉(zhuǎn)化成熱能,邊界層內(nèi)的流動為非等熵過程,同時存在動能交換和熱交換過程,因此溫度邊界層處的溫度要高于中心主氣流.
圖6 t=1.50 ms時,沿軸向位置管內(nèi)溫度(K)局部分布(Case4:p4=500 kPa,p1=20 kPa,T4,1=300 K)
此時,激波管內(nèi)形成四個區(qū)域,沿軸向位置溫度分布見圖7,膜片破裂后,產(chǎn)生的正激波由膜片位置向低壓段傳播,正激波與高低壓段氣體接觸面之間的區(qū)域稱為2區(qū),產(chǎn)生的膨脹波向高壓段傳播,膨脹波與接觸面之間的區(qū)域稱為3區(qū);正激波到達右側(cè)尾端端壁后會發(fā)生反射產(chǎn)生反射激波,反射激波經(jīng)過的區(qū)域稱為5區(qū),即試驗區(qū).粘性流體流經(jīng)管體壁面時,入射激波后壁面附近的流體區(qū)域會形成一定厚度的邊界層,邊界層內(nèi)激波速度會迅速降低,在管壁處的速度為零,激波管內(nèi)2區(qū)、3區(qū)及5區(qū)壁面因流體粘性效應(yīng)附近均會出現(xiàn)邊界層,1區(qū)和4區(qū)的氣體處于靜止?fàn)顟B(tài)不存在邊界層.
圖7 t=1.50 ms時,沿軸向位置管內(nèi)溫度(K)分布范圍(Case4:p4=500 kPa,p1=20 kPa,T4,1=300 K)
采用不同時刻的壓力和溫度云圖來反映激波移動的過程,見圖8~9,給出了激波產(chǎn)生、傳播及反射過程中流體的壓力和溫度變化云圖.t=0.1 ms時,膜片消失后,入射激波向低壓段傳播,經(jīng)過的2區(qū)壓力和溫度迅速升高,膨脹波向高壓段傳播,膨脹波后的區(qū)域壓力和溫度會降低;t=1.0 ms時,由壓力云圖可以清晰觀察到膨脹波和入射激波面的運行位置,根據(jù)溫度云圖可觀察到入射激波與接觸面之間的2區(qū),隨著時間推移,2區(qū)范圍會逐漸增大,這是因為激波運行速度大于接觸面前進速度;當(dāng)t=2.5 ms時,入射激波抵達右尾端端壁并產(chǎn)生反射激波向左傳播,經(jīng)過的5區(qū)壓力和溫度再一次升高達到穩(wěn)定狀態(tài)并持續(xù)一段時間,膨脹波到達左端端壁后發(fā)生反射產(chǎn)生反射膨脹波向右傳播;t=3.0 ms時,反射激波向左傳播與接觸面相遇并產(chǎn)生二次反射激波,二次反射激波向右傳播;t=3.2 ms時刻,二次反射激波會對5區(qū)進行壓縮作用,使5區(qū)壓力和溫度不均衡上升,試驗區(qū)的氣體狀態(tài)發(fā)生改變,試驗區(qū)環(huán)境被破壞,距尾端最遠的測點壓力穩(wěn)定狀態(tài)最先被破壞,當(dāng)二次反射激波行至尾端端壁,5區(qū)有效試驗時間結(jié)束;t=3.5 ms時,二次反射激波已到達尾端端壁,再次反射產(chǎn)生三次反射激波,這個過程會持續(xù)多次,直至激波管內(nèi)整體壓力處于平衡狀態(tài).
圖8 Case4激波管內(nèi)流體瞬時壓力(Pa)云圖
圖9 Case4激波管內(nèi)流體溫度(K)云圖
以空氣(理想狀態(tài))為例,計算移除高壓段與低壓段之間膜片后可以獲得的激波管5區(qū)有效試驗時間,將反射激波后的流動特性穩(wěn)定的時間定義為試驗時間,試驗時間是評估激波管設(shè)備性能的重要參數(shù),可為進行著火延時測量提供重要理論參考依據(jù).
激波管中流體運動位置-時間平面圖是一種重要的工具,它可以很好地估計在移除膜片后可獲得的有效試驗測試時間.圖10~11為Case4(air-air,p4=500 kPa,p1=20 kPa,T4,1=300 K)運行時壓力和溫度分布的預(yù)測x-t圖.圖10壓力分布x-t圖中接觸面沒有出現(xiàn),因為它連續(xù)的跟隨入射激波后,圖11為整個激波管設(shè)備的溫度變化,激波、膨脹波和接觸面都呈現(xiàn)在圖中,激波之后是接觸面,直到反射激波與接觸面相互作用,然后被進一步反射.膨脹波到達左端壁面會產(chǎn)生反射膨脹波,由圖11可知,膨脹波1曲線斜率小于和接觸面的曲線斜率,即接觸面?zhèn)鞑ニ俣却笥诜瓷渑蛎洸ǎ虼?,高壓段長度足夠條件下,反射膨脹波到達低壓段末端之前,激波就已與接觸面相交.二次反射激波到達右端壁面時,標(biāo)志激波管測試時間結(jié)束.
圖10 壓力云圖中,激波的位置隨時間變化曲線
圖11 溫度云圖中,激波和高低壓段間的接觸面的位置隨時間變化曲線
在激波管壁面設(shè)置了8個監(jiān)測點,另外在低壓段端面中心設(shè)置了1個監(jiān)測點,圖12為Case4工況下激波管壁面監(jiān)測點壓力隨時間變化情況.開始時,高壓段測點S1保持初始壓力不變,膨脹波經(jīng)過后壓力持續(xù)降低,待反射膨脹波達到時再次降低;低壓段壁面測點S2~S8入射激波到達時壓力均有階躍升高,S6~S8處反射激波經(jīng)過,S6處距尾端較遠,反射激波已遇到接觸面產(chǎn)生二次反射激波,反射強度下降,因此壓力二次升高值低于測點S7、S8.測點S9不同于側(cè)壁測點S8,只有一次壓力階躍過程,入射激波與尾端壁面相遇時刻就是5區(qū)的開始.
圖12 激波管壁面監(jiān)測點壓力變化
實際試驗過程中,通常有兩種方式表示測量時間,分別為側(cè)壁傳感器壓力穩(wěn)定時間t側(cè)和端壁傳感器壓力穩(wěn)定時間t端,壓力測點相當(dāng)于實際設(shè)備上的壓力傳感器,以上方式也是5區(qū)試驗時間限度的重要參考參數(shù).表5為激波管試驗時間限度計算結(jié)果比較.
表5 激波管試驗時間限度計算結(jié)果比較
由表5可知,該工況下幾種試驗時間限度估計方式得出的結(jié)果差異較小,側(cè)壁測點S8得到的壓力穩(wěn)定時間為1.039 ms,比端壁測點時間S9壓力穩(wěn)定時間1.073 ms要短,這是因為邊界層效應(yīng)的影響,反射激波會發(fā)生激波分岔和激波衰減,離反射端面處的距離越近,受邊界層影響越小,側(cè)壁測點壓力開始上升至穩(wěn)定的時間稍長于尾端端壁測點.然而,激波管在該算例下x-t圖得到的測試時間為1.169 ms,大于端壁測點的壓力穩(wěn)定時間,綜上所述,該算例理論有效試驗時間為1.169 ms,實際測量有效試驗時間最大為1.073 ms.
1) 采用FLUENT中的湍流模型進行激波管流體流動數(shù)值模擬,計算結(jié)果與層流模型模擬結(jié)果以及實驗結(jié)果進行比較,與實驗結(jié)果吻合較好,表明湍流模型對激波管流動具有較好的模擬效果,本文CFD數(shù)值方法是高效可靠的.
2) 詳細(xì)模擬了激波的產(chǎn)生和傳播過程,給出了激波形成、傳播、反射和激波與接觸面相互作用的整個過程,接觸面與反射激波相遇時刻和二次反射激波行至端壁的時間是有效試驗時間限度的重要影響因素.
3) 本數(shù)值方法能夠提供x-t關(guān)系圖,從該圖中可以確定激波管實驗測量有效持續(xù)時間,激波管5區(qū)內(nèi)側(cè)壁、端壁壓力傳感器測得的試驗時間長度均小于x-t圖獲得的最大的測試時間,且差值較小,符合預(yù)期結(jié)果.