吳博文 章利特 余秋李 劉天程 馮子龍
摘 要:為研究激波透過大孔隙率結(jié)構(gòu)化球陣的氣體流動(dòng)特性,采用Fluent軟件三維雷諾平均Navier-Stokes方程和k-epsilon湍流模型,對(duì)不同顆粒(或固相)體積分?jǐn)?shù)和入射馬赫數(shù)時(shí)激波透過面心立方球陣的氣相流場(chǎng)進(jìn)行了數(shù)值模擬,重點(diǎn)分析了球陣內(nèi)部的流動(dòng)特性和反射波強(qiáng)度的參數(shù)影響。研究結(jié)果表明,不同顆粒體積分?jǐn)?shù)和入射馬赫數(shù)下FCC球陣對(duì)入射激波的反射效應(yīng)有較大不同。隨著入射馬赫數(shù)的增大,反射效應(yīng)不斷增強(qiáng),當(dāng)馬赫數(shù)增大到某一確定值后,反射效應(yīng)則趨于穩(wěn)定。在顆粒體積分?jǐn)?shù)為5%時(shí),此確定馬赫數(shù)值應(yīng)介于1.9與2.3之間。而隨著體積分?jǐn)?shù)增大,球陣對(duì)入射激波的反射效應(yīng)先增強(qiáng)后基本趨于穩(wěn)定。當(dāng)入射馬赫數(shù)固定為1.4時(shí),體積分?jǐn)?shù)為7%的FCC球陣具有最強(qiáng)的反射效應(yīng)。
關(guān)鍵詞:FCC球陣;激波;反射波;體積分?jǐn)?shù)
中圖分類號(hào):TH311 文獻(xiàn)標(biāo)識(shí)碼:A 文章編號(hào):1671-2064(2019)23-0220-04
0 引言
激波與固體顆粒群的相互作用是超聲速兩相流氣體動(dòng)力學(xué)中的一個(gè)重要研究課題[1],激波驅(qū)動(dòng)固體顆粒群理論[2]同時(shí)也被廣泛運(yùn)用于軍事,醫(yī)療,航空航天等領(lǐng)域。迄今為止,國(guó)內(nèi)外已經(jīng)有很多學(xué)者針對(duì)激波驅(qū)動(dòng)固體顆粒的流動(dòng)現(xiàn)象,采用實(shí)驗(yàn)和數(shù)值模擬的方式進(jìn)行研究。
Sun等[3]開展了對(duì)平面激波透過單球時(shí)的非穩(wěn)態(tài)阻力的實(shí)驗(yàn)測(cè)量,發(fā)現(xiàn)雷諾數(shù)和單球直徑的增大,會(huì)造成非穩(wěn)態(tài)阻力的線性遞增;Saito等[4]針對(duì)含激波攜帶氣固混合物的兩相流動(dòng)開展一系列數(shù)值模擬,研究了激波后影響非平衡流場(chǎng)結(jié)構(gòu)的因素主要為非穩(wěn)態(tài)阻力,并發(fā)現(xiàn)激波馬赫數(shù)的降低和顆粒質(zhì)量比的升高對(duì)阻力系數(shù)的時(shí)間依賴效應(yīng)的影響更深;Igra等[5]對(duì)模型球在非穩(wěn)態(tài)流場(chǎng)中的流動(dòng)進(jìn)行了數(shù)值模擬,并對(duì)模型球進(jìn)行了受力分析,發(fā)現(xiàn)球表面壓力和粘性剪切應(yīng)力決定著球的阻力;Tanno等[6]利用垂直激波管展開了對(duì)模型球的實(shí)驗(yàn)與數(shù)值模擬,發(fā)現(xiàn)激波剛接觸到模型球時(shí),所受到的非穩(wěn)態(tài)阻力最大,衍射波停滯在模型球后,會(huì)在短時(shí)間內(nèi)使阻力值下降至負(fù);Jourdan等[7]利用水平激波管開展了激波驅(qū)動(dòng)面粉顆粒群的流動(dòng)相關(guān)實(shí)驗(yàn),表明入射激波的初始條件和顆粒的物理因素決定顆粒噴射的數(shù)量;Shi等[8]通過激波管裝置對(duì)激波與球陣的相互作用進(jìn)行實(shí)驗(yàn)研究,表明球陣球與球之間存在不可忽視的干擾力,激波驅(qū)動(dòng)單球的實(shí)驗(yàn)結(jié)果不可直接運(yùn)用在多球?qū)嶒?yàn)中;Zhang等[9]對(duì)激波驅(qū)動(dòng)固體顆粒群的流場(chǎng)研究后發(fā)現(xiàn),在激波加載下的顆粒群存在一系列復(fù)雜現(xiàn)象,包括激波的透射、衍射、反射以及激波與激波之間的相互擾動(dòng)現(xiàn)象。
本文對(duì)氣固兩相不同參數(shù)下激波與結(jié)構(gòu)化球陣互相作用過程進(jìn)行了非定常數(shù)值模擬,分析了可壓縮性氣體與固體球陣相互作用的機(jī)理,主要包括入射激波馬赫數(shù),球陣體積分?jǐn)?shù)等對(duì)整體流動(dòng)的影響,解析了激波遇到不同參數(shù)結(jié)構(gòu)化球陣的反射和衍射現(xiàn)象。
1 幾何模型
由于真實(shí)的固體顆粒群形狀大小排列方式各不相同,為方便對(duì)不同顆粒群進(jìn)行研究,把固體顆粒假設(shè)成具有周期性的固定空間球陣,從物理及數(shù)學(xué)上抽象出顆粒群的簡(jiǎn)化模型,建立其表征體,即球陣模型。為了更好的描述出球陣模型的具體特性,采用自然界常見的面心立方堆積(FCC)排列方式對(duì)球陣進(jìn)行排列,其立方體單元結(jié)構(gòu)陣列如圖1所示。
本文中球體固定直徑40mm,結(jié)合單元體幾何特性,選取1%,3%,5%,7%,9%五種體積分?jǐn)?shù),計(jì)算出每種體積分?jǐn)?shù)對(duì)應(yīng)單元體的尺寸,選擇五個(gè)單元體沿流向方向并排放置,為防止前后單元體球陣的擾動(dòng),選擇中心單元體進(jìn)行監(jiān)測(cè)。采用參數(shù)化三維造型軟件Pro/E對(duì)整個(gè)FCC球陣的流道進(jìn)行建模,為了保證流動(dòng)的充分發(fā)展,對(duì)FCC球陣進(jìn)出口段進(jìn)行了延長(zhǎng)。表1展示了各體積分?jǐn)?shù)下FCC球陣的具體幾何參數(shù)。
2 數(shù)值模擬的網(wǎng)格劃分及算法
2.1 網(wǎng)格劃分
基于對(duì)稱性,采用ANSYS Workbench軟件對(duì)1/8個(gè)流道進(jìn)行非結(jié)構(gòu)化四面體網(wǎng)格劃分,為了在兼顧計(jì)算精度的同時(shí)節(jié)省計(jì)算資源,對(duì)流道整體網(wǎng)格尺寸設(shè)置為4mm,對(duì)球面網(wǎng)格尺寸加密為0.4mm,在稀疏網(wǎng)格過渡區(qū)域采用指數(shù)增長(zhǎng)方式進(jìn)行過渡。最終得到了總體網(wǎng)格質(zhì)量在0.6以上的高質(zhì)量非結(jié)構(gòu)化網(wǎng)格,整體網(wǎng)格和球面局部加密網(wǎng)格見圖2。
2.2 控制方程
采用三維非定常雷諾平均Navier-Stokes方程作為氣體(空氣)流動(dòng)的基本控制方程,湍流模型選用k-epsilon湍流模型。除連續(xù)性方程和能量方程外,動(dòng)量方程表達(dá)形式為:
(1)
其中:U為速度;ρ為液體密度;p為壓強(qiáng);SM為控制方程源項(xiàng);τ為分子應(yīng)力張量。
標(biāo)準(zhǔn)k-epsilon湍流模型的方程表達(dá)式為:
(2)
(3)
其中:μ為動(dòng)力粘度,空氣粘度采用Sutherland定律進(jìn)行計(jì)算;μt為湍流粘性系數(shù);Gk為平均速度梯度形成的湍動(dòng)能生成項(xiàng);Gb是存在浮力的影響形成的湍動(dòng)能的生成項(xiàng);ε為湍動(dòng)能耗散率;k為湍動(dòng)能;YM是可壓速湍流脈動(dòng)膨脹對(duì)總的耗散率的影響;、、為模型參數(shù)。在對(duì)空氣密度進(jìn)行求解時(shí)采用理想氣體假設(shè),其狀態(tài)方程為:
(4)
其中:R為通用氣體常數(shù),T為空氣溫度。
2.3 數(shù)值計(jì)算方法和邊界條件
為了能夠更好地捕捉精細(xì)的激波結(jié)構(gòu),選用AUSM差分格式和Roe-fds通量分裂方法,方程采用二階迎風(fēng)格式進(jìn)行離散,并采用隱式算法進(jìn)行瞬態(tài)迭代求解,固定時(shí)間步長(zhǎng)為5E-7s。在進(jìn)行數(shù)值計(jì)算時(shí),工作介質(zhì)選為空氣,采用基于密度的求解器。計(jì)算區(qū)域的進(jìn)口邊界條件設(shè)定為壓力入口,參數(shù)設(shè)置為波后總溫和總壓,具體數(shù)值根據(jù)激波管理論得到;計(jì)算域的出口邊界條件壓力出口,壓力溫度設(shè)置為環(huán)境溫度和環(huán)境壓力;球體表面設(shè)置為固壁無滑移、絕熱邊界條件。由于采用1/8個(gè)流道進(jìn)行計(jì)算,其余三個(gè)面均為對(duì)稱面,設(shè)置為對(duì)稱邊界條件。在計(jì)算開始前在球陣入口100mm處patch一個(gè)區(qū)域,將此區(qū)域溫度、壓力和速度均設(shè)置為波后靜溫、靜壓和波后氣流速度。