一種基于CFD-DEM的箔條散射運(yùn)動數(shù)值模擬方法

王棟歡，肖洪，吳丁毅

（西北工業(yè)大學(xué)動力與能源學(xué)院，陜西 西安 710129）

0 引言

在雷達(dá)無源對抗中，箔條無源干擾是自二戰(zhàn)以來使用最為廣泛的無源干擾技術(shù)之一［1］，它通過大量噴撒鋁質(zhì)或者表面涂覆金屬的條狀纖維構(gòu)成的空中干擾物，能夠在一定的空間范圍內(nèi)產(chǎn)生干擾回波，以掩護(hù)己方目標(biāo)不被攻擊［2］。

目前針對箔條無源干擾技術(shù)的研究主要分為兩個方面，一方面是箔條云團(tuán)運(yùn)動擴(kuò)散的研究，另一方面則是關(guān)于箔條云團(tuán)雷達(dá)散射截面（RCS）以及回波多普勒特性的研究［3］。散射箔條云團(tuán)的運(yùn)動是一種特殊的氣固兩相流動。眾多學(xué)者致力于發(fā)展數(shù)值方法來研究氣固兩相流動，這些方法大體可以分為以下3種：歐拉-歐拉法、歐拉-拉格朗日法、拉格朗日-拉格朗日法［4-7］。歐拉-歐拉法提出將離散粒子群視為一種偽流體的假設(shè)，由此得到了描述準(zhǔn)粒子流體及其相互作用的3個準(zhǔn)粒子流體守恒方程［8-10］。歐拉-拉格朗日法采用離散元方法（DEM）或離散相等模型（DPM）處理離散顆粒的運(yùn)動，其大致方法是將牛頓運(yùn)動方程應(yīng)用于單個顆?；蛱囟w粒群，通過對單個或多個顆粒建立不同的阻力模型，實(shí)現(xiàn)氣固兩相耦合［11-14］。拉格朗日-拉格朗日法在處理粒子運(yùn)動問題上與歐拉-拉格朗日法基本類似，不同的是對氣體運(yùn)動的描述采用基于納維-斯托克斯方程或格子-玻爾茲曼模擬進(jìn)行求解［15-17］。與歐拉-歐拉方法相比，后兩種方法對顆粒運(yùn)動的描述是基于單個或多個真實(shí)粒子的，因而具有更高的精度和真實(shí)性。然而，受限于計算能力和計算資源，拉格朗日-拉格朗日法在工業(yè)生產(chǎn)及工程應(yīng)用中使用較少，相比較而言，基于歐拉-拉格朗日方法的CFD-DEM耦合模型在氣固兩相流數(shù)值研究中顯示出越來越多的優(yōu)越性。

在歐拉-拉格朗日方法中，基于CFD-DEM的耦合模型被廣泛應(yīng)用［18-20］，DEM在處理單個顆粒物運(yùn)動上優(yōu)勢大，然而在針對類似于箔條等不規(guī)則條狀或大顆粒物占據(jù)多個CFD網(wǎng)格的氣固耦合中仍舊面臨巨大挑戰(zhàn)，如圖1所示。本文基于CFD常用軟件Fluent的用戶自定義函數(shù)（UDF），結(jié)合DEM，建立了適用于氣體-條狀物混合流動的CFD-DEM耦合模型，探究了CFD-DEM耦合算法的詳細(xì)計算流程，同時應(yīng)用該模型對箔條散射結(jié)構(gòu)進(jìn)行了預(yù)測，給出了箔條散射運(yùn)動過程的合理化解釋。

圖1 條狀/固體顆粒物占據(jù)多個CFD網(wǎng)格示意圖

1 數(shù)學(xué)模型

1.1 氣體流動控制方程

應(yīng)用CFD-DEM研究氣體-顆粒的混合流動，氣體流動遵循質(zhì)量守恒和動量守恒，其矢量形式的守恒方程為：

式中，氣體粘性應(yīng)力張量tg為：

式中，α為空隙率，即氣相體積分?jǐn)?shù)；μg表示剪切粘性的動力粘性系數(shù)；Sp是動量交換源項，代表顆粒相施加于流體相單位體積的阻力，該阻力與顆粒所受曳力為相互作用力；?為哈密頓算符，其展開式為：

Sp體現(xiàn)了兩相耦合關(guān)系，其影響因素眾多，可由下式計算得出：

由式（5）可知，動量交換源項Sp等于作用在網(wǎng)格單元內(nèi)流體阻力FD的總和。式中V為CFD網(wǎng)格單元體積。由牛頓第三定律可知FD為流體對顆粒曳力的反作用力，在CFD-DEM耦合中，典型的曳力模型為3類，將在后文具體給出。

氣相體積分?jǐn)?shù)α可用式α=1-αs簡單求得，其中αs是顆粒的體積分?jǐn)?shù)，通過抽取樣本點(diǎn)的方法計算顆粒的體積分?jǐn)?shù)簡單有效。首先對每個顆粒單元抽取均勻分布的樣本點(diǎn)Nsample，然后對每個樣本點(diǎn)進(jìn)行檢驗，確定出顆粒i位于某一CFD網(wǎng)格單元中的樣本點(diǎn)數(shù)Si，因此可得到顆粒在某一CFD網(wǎng)格單元中的體積分?jǐn)?shù)為：

式中，n表示某一CFD網(wǎng)格單元中顆粒的數(shù)目。顯然，樣本點(diǎn)個數(shù)Nsample越大，αs的值越精確。對于顆粒體積小于網(wǎng)格體積60%的情況，可以直接利用公式求解(Nsample=1)即能達(dá)到滿足要求的求解精度。

1.2 湍流模型

考慮到湍流的影響，基于雷諾平均法（RANS）對氣體流動控制方程進(jìn)行時均化處理，并采用Reynolds應(yīng)力方程模型（RSM）。為了簡化計算，忽略浮力、系統(tǒng)旋轉(zhuǎn)及固體壁面反射的影響。由于RSM模型適用于充分發(fā)展的湍流，對于近壁面的流動，Re很低，RSM模型不再適用，因此，采用避面函數(shù)法求解近壁區(qū)域的流體速度。

1.3 不規(guī)則大顆粒物運(yùn)動的數(shù)學(xué)模型

在DEM中，應(yīng)用牛頓第二定律對顆粒的運(yùn)動進(jìn)行描述。質(zhì)量為mi、體積為Vpi的顆粒i平移及旋轉(zhuǎn)運(yùn)動的控制方程可表示為：

運(yùn)用運(yùn)動學(xué)方程更新每個顆粒位置即：

式 中，F(xiàn)cn，ij和Fct，ij分 別 表 示 顆 粒i與 顆 粒j的 法 向 碰 撞力和切向碰撞力；Ii、wi、Ri和nij分別表示顆粒i的轉(zhuǎn)動慣量、角速度矢量、顆粒球心到碰撞平面的距離以及顆粒i與顆粒j之間的法向單位矢量。

這里，(∑Fi)表示顆粒i受到的合力，包含曳力、重力、浮力、壓力梯度力、附加質(zhì)量力、Basset力、Saffman升 力、Magnus升 力、熱 泳 力 等。在 本 文 的CFD-DEM數(shù)值模型中只考慮重力和曳力。此外本文采用簡單的線性模型即Hertz-Mindlin非滑動接觸模型描述不規(guī)則大顆粒物之間的相互作用。其相互作用力通過滑塊、減振器和彈簧進(jìn)行模擬［21-23］。

考慮到不規(guī)則大顆粒物占據(jù)多個CFD網(wǎng)格的情況，在自由流阻力模型的基礎(chǔ)上引入當(dāng)?shù)鼐W(wǎng)格投影面積Alocal，因此，作用在單個顆粒上的曳力描述如下：

式中，m為顆粒占據(jù)的網(wǎng)格數(shù)；ug為氣體速度；us為顆粒運(yùn)動速度；CD為單顆粒曳力系數(shù)，其大小與顆粒雷諾 數(shù)Rep有 關(guān)［24-25］；Alocal是 顆粒在當(dāng)?shù)鼐W(wǎng)格 下的投 影面積，如圖2所示。

圖2 條狀/不規(guī)則大顆粒物占據(jù)多個CFD網(wǎng)格投影面積示意圖

對顆粒接觸力的求解采用DEM，利用牛頓第二定律進(jìn)行各個顆粒的運(yùn)動求解時采用了整體坐標(biāo)系，然而顆粒接觸力的計算是在局部坐標(biāo)系下完成的，因此在顆粒接觸力計算部分需要進(jìn)行坐標(biāo)變換。將整體坐標(biāo)中的相對位移轉(zhuǎn)換為局部坐標(biāo)中的相對位移，之后才能直接應(yīng)用相關(guān)接觸模型進(jìn)行接觸力的計算，并且在求解得到接觸力的值后，還需進(jìn)行第二次坐標(biāo)轉(zhuǎn)換，將法向、切向方向上的接觸力分量投影變換到整體坐標(biāo)系當(dāng)中，進(jìn)行力的求和計算［26］。顆粒接觸力的計算大體可分為：

第一步：載入模型的特征參數(shù)；

第二步：讀取接觸配對單元在整體坐標(biāo)系下的坐標(biāo)值，計算相對位移并進(jìn)行坐標(biāo)轉(zhuǎn)換；

第三步：利用接觸模型中的接觸力-位移關(guān)系計算彈性力和阻尼力；

第四步：計算切向力和顆粒單元的接觸力在絕對坐標(biāo)系上的投影分量，以及旋轉(zhuǎn)力矩；

第五步：判斷接觸對象是顆粒還是邊界，計算接觸對象的受力。

1.4 CFD-DEM耦合模型

氣體與顆粒之間的耦合可分為3種：單相耦合、雙向耦合和四相耦合。單相耦合只考慮氣體對顆粒的作用力；雙向耦合同時考慮氣體對顆粒、顆粒對氣體的相互作用力；四相耦合則是在雙相耦合的基礎(chǔ)上還考慮了顆粒與顆粒間的相互作用。對于顆粒體積分?jǐn)?shù)小于10%的稀相流，考慮氣體-顆粒的雙相耦合，應(yīng)用Fluent自帶的DPM模型可以得到較滿意的結(jié)果；而對于顆粒體積分?jǐn)?shù)較大的氣體-顆?；旌狭鲃樱杩紤]氣體-顆粒間的四相耦合。CFD-DEM耦合模型的最大優(yōu)點(diǎn)是同時考慮了氣固間的動量交換、顆粒體積分?jǐn)?shù)對流場的影響以及顆粒間的相互作用，其適用性更廣，求解精度更高。

在CFD-DEM耦合計算過程中，耦合模塊基于Fluent的UDF函數(shù)及DEM模塊的輸出接口，實(shí)現(xiàn)了Fluent與DEM模塊之間的數(shù)據(jù)交換。該耦合計算過程中，F(xiàn)luent首先在一個時間步長內(nèi)求解氣相流動，同時調(diào)用CFD-DEM耦合模塊計算顆粒的曳力并傳遞給DEM模塊，接著DEM模塊計算每個顆粒的運(yùn)動并通過耦合接口將空隙率及動量源項傳遞給Fluent，最后Fluent加載源項及空隙率并進(jìn)行下一時間步內(nèi)的迭代求解，其耦合計算流程如圖3所示。一般情況下，在求解流場的一個時間步內(nèi)，DEM模塊迭代計算10到100次。CFD計算時間步長（TF）與DEM模塊時間步長（TDEM）的大小關(guān)系為：

圖3 CFD-DEM耦合計算流程

2 數(shù)值模擬

采用發(fā)展的CFD-DEM耦合方法，對箔條散射過程中的箔條云運(yùn)動進(jìn)行了數(shù)值模擬。選擇拉瓦爾噴管附加直徑280 mm、長度390 mm的圓柱腔作為箔條散射的研究模型。拉瓦爾噴管長110 mm，入口直徑28 mm，出口直徑50 mm，喉部直徑20 mm。幾何結(jié)構(gòu)如圖4所示。采用Block塊網(wǎng)格劃分方法，整個計算域被劃分為結(jié)構(gòu)化的六面體網(wǎng)格，圖5為計算模型結(jié)構(gòu)及網(wǎng)格。

圖5 計算模型結(jié)構(gòu)化網(wǎng)格圖

選取箔條數(shù)量為1 000個，初始任意堆積于噴管近出口內(nèi)部，如圖4所示。箔條形狀簡化為圓柱體結(jié)構(gòu)，密度為500 kg/m3，直徑為0.2 mm，長度為15 mm。具體參數(shù)及計算的邊界條件為：拉法爾噴管進(jìn)口直徑為28 mm，長度L為110 mm，出口直徑為50 mm，網(wǎng)格數(shù)量為745 324個，喉部直徑為20 mm，箔條直徑為0.2 mm，長度為15 mm，數(shù)量為1 000個；邊界條件進(jìn)氣總壓為2 000 Pa，Ma數(shù)為0.1。

圖4 噴管結(jié)構(gòu)圖

圖6 給出了不同時刻下箔條的散射狀態(tài)，由圖可得初始靜止?fàn)顟B(tài)的箔條在氣體曳力的作用下瞬間從噴管噴出，在百分之一秒的時間尺度上即可呈現(xiàn)出散狀分布。初始靜止的一堆箔條經(jīng)0.02 s后平均速度達(dá)到8 m/s，箔條云形狀似一半球。隨著時間的推進(jìn)，中心區(qū)域箔條被吹向最前端并逐漸散開，箔條云中心區(qū)域速度最大，邊緣處箔條速度較小。在0.02 s時刻，箔條云變?yōu)殄F形，最大速度出現(xiàn)在錐頂處，速度大約為13 m/s。0.04 s時刻下錐狀箔條云被拉長，最前端箔條的速度可達(dá)到16 m/s。在氣流的作用下，中心箔條的運(yùn)動處于一個加速的過程，當(dāng)中心箔條被完全吹開后速度逐漸降低，受圓柱形空腔長度的限制，當(dāng)前構(gòu)型尺寸下未能顯現(xiàn)出箔條云完全被吹開時箔條的最大位移。

圖6 不同時刻箔條散射結(jié)構(gòu)圖

3 結(jié)束語

本文提出了一種適用于求解箔條散射運(yùn)動的CFD-DEM耦合方法，利用該方法對1 000數(shù)量箔條散射結(jié)構(gòu)進(jìn)行了仿真模擬。研究結(jié)果表明，箔條散射結(jié)構(gòu)初始為半球體，進(jìn)而發(fā)展為椎體，箔條云中心區(qū)域速度最大，邊緣處箔條速度較小。模擬結(jié)果初步表明本文提出的CFD-DEM耦合方法可以模擬箔條散射運(yùn)動，且能得到符合預(yù)期的箔條散射結(jié)構(gòu)。

本文對于箔條散射的數(shù)值模擬只停留在數(shù)值仿真層面，并未涉及與之對應(yīng)的實(shí)驗驗證。此外，現(xiàn)有的箔條散射數(shù)值模擬算例工況較為單一，只給出了單一速度工況下箔條云團(tuán)隨時間的初步散射結(jié)構(gòu)特性。因此，今后還需通過實(shí)驗和具體工程問題，進(jìn)行多種工況的仿真計算。