激波透過大孔隙率結(jié)構(gòu)化球陣的數(shù)值模擬

吳博文 章利特 余秋李 劉天程 馮子龍

摘 要：為研究激波透過大孔隙率結(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)域。迄今為止，國內(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è)置為波后靜溫、靜壓和波后氣流速度。