李慈應(yīng)　李姝霖

【摘? 要】在探月工程中，月球探測(cè)器的軟著陸過(guò)程中必須考慮到緩沖發(fā)動(dòng)機(jī)羽流激起月壤塵埃的問(wèn)題，因?yàn)閾P(yáng)起的月壤很可能對(duì)探測(cè)器設(shè)備產(chǎn)生多種不良的影響。本文以緩沖發(fā)動(dòng)機(jī)羽流與月壤所形成的氣固兩相流場(chǎng)為主要研究對(duì)象，介紹了如何應(yīng)用FLUENT-EDEM兩相耦合模塊對(duì)其進(jìn)行模擬仿真。以一個(gè)簡(jiǎn)化模型作為實(shí)例，闡述了該耦合方法相較于以往的兩相流仿真所具有的優(yōu)勢(shì)。

【關(guān)鍵詞】?jī)上嗔?月壤;耦合;EDEM

引言

實(shí)現(xiàn)月球探測(cè)器的軟著陸是我國(guó)探月工程第二階段的關(guān)鍵技術(shù)，而火箭發(fā)動(dòng)機(jī)反向噴氣是一種常見(jiàn)的著陸緩沖技術(shù)。但在探測(cè)器軟著陸過(guò)程中緩沖發(fā)動(dòng)機(jī)羽流激起的月壤可能會(huì)衰減導(dǎo)航敏感器件的探測(cè)信號(hào)，甚至導(dǎo)致儀器無(wú)法正常工作，激起的月塵還會(huì)造成的視覺(jué)模糊，嚴(yán)重影響探測(cè)器的安全及準(zhǔn)確著落，因此探月過(guò)程中必須考慮羽流中的月壤對(duì)探測(cè)器所能造成的影響。對(duì)氣固兩相流的研究于探月工程具有不可忽視的積極意義，本文將要討論CFD-EDEM耦合方法在氣固兩相流研究中的應(yīng)用。

1 概述

氣固兩相流的理論分析相較于單相流的分析困難得多，而描述兩相流的通用微分方程組至今尚未建立。大量理論工作主要采用以下兩類簡(jiǎn)化模型：

1）均相模型。將兩相介質(zhì)看成是一種混合得非常均勻的混合物，假定處理單相流動(dòng)的概念和方法仍然適用于兩相流，但需對(duì)它的物理性質(zhì)及傳遞性質(zhì)作合理的假定;

2）分相模型。認(rèn)為單相流的概念和方法可分別用于兩相系統(tǒng)的各個(gè)相，同時(shí)考慮兩相之間的相互作用。

目前，對(duì)如何處理氣固兩相流在數(shù)值計(jì)算上主要采用四種方法：一是將某相看成連續(xù)的，根據(jù)連續(xù)理論導(dǎo)出歐拉基本方程，稱為歐拉法;二是將某相視為不連續(xù)的離散相，對(duì)每個(gè)質(zhì)點(diǎn)進(jìn)行拉格朗日追蹤，稱為拉格朗日法。將上述的這兩種方法綜合起來(lái)，氣固兩相流的研究方法有歐拉-歐拉、歐拉-拉格朗日、拉格朗日-拉格朗日三種方法。第四種方法就是基于分子運(yùn)動(dòng)理論的直接模擬蒙特卡洛法與計(jì)算流體力學(xué)方法相結(jié)合對(duì)氣固兩相流數(shù)值計(jì)算的方法。

而本文所要探討的CFD-EDEM耦合方法相較于上述四種方法，則具有更加獨(dú)到的優(yōu)勢(shì)，我們將在下文作展開(kāi)說(shuō)明。

2 CFD-EDEM耦合方法

2.1 EDEM軟件簡(jiǎn)介

仿真軟件EDEM的主要功能是仿真、分析和觀察顆粒流的運(yùn)動(dòng)規(guī)律，它是全球第一個(gè)被專門設(shè)計(jì)為使用離散元模型技術(shù)作為基礎(chǔ)的通用CAE軟件。

EDEM的Creater（前處理模塊）可以快速、簡(jiǎn)便的為散貨顆粒進(jìn)行參數(shù)化建模。

Simulator（求解器）基于先進(jìn)的離散元方法，結(jié)合了經(jīng)典的碰撞模型，在拉格朗日框架下，離散求解系統(tǒng)的運(yùn)動(dòng)學(xué)方程、動(dòng)力學(xué)方程和本構(gòu)方程。從而獲得每個(gè)顆粒及壁面單元體的速度、受力等參數(shù)。

Analyst是EDEM的后處理模塊，含有很多工具對(duì)計(jì)算結(jié)果進(jìn)行加工處理，方便地獲得顆粒的速度、能量分布，設(shè)備的受力情況等，并可利用統(tǒng)計(jì)方法對(duì)各物理量間的相關(guān)性進(jìn)行分析。

2.2 Fluent-EDEM耦合方案

EDEM不僅可以單獨(dú)進(jìn)行顆粒力學(xué)模擬，也可以和通用CFD軟件Fluent耦合進(jìn)行流體-顆粒系統(tǒng)的模擬。其主要特點(diǎn)如下：

1）EDEM-Fluent耦合計(jì)算時(shí)，EDEM的耦合模塊內(nèi)置于Fluent軟件中，兩者無(wú)縫連接;2）EDEM能夠自動(dòng)讀取Fluent中的網(wǎng)格;3）EDEM可以使用Fluent中離散相模型和歐拉多項(xiàng)流模型進(jìn)行多相流耦合;4）采用EDEM直接對(duì)顆粒運(yùn)動(dòng)進(jìn)行模擬;5）完全雙向動(dòng)力耦合。

Fluent-EDEM可以模擬絕大多數(shù)的流體-顆粒系統(tǒng)，特別是在以下方面，更是具有不可代替的優(yōu)勢(shì)：

?顆粒尺度分布范圍廣;非球形顆粒;流體中充滿顆粒，顆粒的體積分?jǐn)?shù)大;顆粒凝聚和分離;機(jī)械設(shè)備上顆粒沉積和黏結(jié);可以直觀化地觀察計(jì)算過(guò)程中顆粒的運(yùn)動(dòng)情況。

Fluent-EDEM耦合模擬流程如下：先由Fluent將某一時(shí)間點(diǎn)的流場(chǎng)計(jì)算收斂，將流場(chǎng)信息通過(guò)曳力模型轉(zhuǎn)化為EDEM中作用在顆粒上的流體曳力，EDEM計(jì)算每個(gè)顆粒所受的外力，并由此更新顆粒的位置、速度等信息，最后這些顆粒屬性以動(dòng)量匯的形式加到CFD計(jì)算中，從而影響流場(chǎng)。

2.3 CFD-EDEM耦合理論

2.3.1 CFD-EDEM升力模型

2.3.1.1 Saffman升力模型

在邊界層，剪切層內(nèi)，或者是流體發(fā)生剪切時(shí)，會(huì)出現(xiàn)較高的速度梯度。當(dāng)顆粒處于高速流體當(dāng)中，較高的速度梯度造成了顆粒表面的壓力差，從而產(chǎn)生了升力。Saffman升力可用下式表達(dá)：

2.3.1.2 Magnus升力模型

Magnus升力是由于顆粒自旋而產(chǎn)生的升力，球體的阻力系數(shù)與雷諾數(shù)較高的情況下的Magnus升力的關(guān)聯(lián)可以表示為下式：

這個(gè)模型即使對(duì)于雷諾數(shù)高達(dá)2000時(shí)對(duì)于顆粒仍然適用。

2.3.1.3 Fluid-Induced Torque（流激誘導(dǎo)轉(zhuǎn)矩）

在以上耦合的模型中，沒(méi)有考慮由顆粒自旋產(chǎn)生的流體剪切力的影響。顯然，若在高速流體中這個(gè)力不能忽略。根據(jù)Dennis等人的研究，有以下公式：

通過(guò)這樣的處理，這一擴(kuò)展模型可用于顆粒雷諾數(shù)在32到1000的情況。

2.3.2 CFD-EDEM阻力模型

Fluent-EDEM耦合采用一種改進(jìn)的自由流阻力來(lái)計(jì)算作用在顆粒球形上的力，對(duì)于顆粒表面為非球面的情況，先用對(duì)顆粒的球型包圍盒進(jìn)行阻力的計(jì)算。這種處理方法僅適合處理兩種境況下的顆粒：（1）那些大小相同或者是體積小于流體網(wǎng)格單元的顆粒;（2）周圍的流體參數(shù)（速度、密度和粘度等）的變化基本為常量的顆粒。阻系數(shù)C，取決于雷諾數(shù)Re：

式中：r 是流體密度，h 為流體的粘度，L 是顆粒球的直徑，v 是顆粒與流體間的相對(duì)速度，a 是CFD網(wǎng)格單元的自由體積。除了阻力，顆粒的固有浮力也必須考慮。浮力的計(jì)算公式是：

一般對(duì)于Lagrangian耦合來(lái)說(shuō)，這種改進(jìn)的阻力法則可以用來(lái)考慮顆粒載荷的影響;而在Eulerian耦合中，顆粒載荷需通過(guò)結(jié)合體積分?jǐn)?shù)的方法來(lái)考慮這一影響。

3 FLUENT-EDEM耦合實(shí)例

3.1 計(jì)算模型建立

氣固兩相流數(shù)值計(jì)算過(guò)程中，顆粒在氣流中場(chǎng)中的運(yùn)動(dòng)是從靜止開(kāi)始計(jì)算的。為了簡(jiǎn)化計(jì)算，同時(shí)使流場(chǎng)設(shè)置不至于過(guò)大，分析探測(cè)器發(fā)動(dòng)機(jī)距月表較近時(shí)（發(fā)動(dòng)機(jī)噴管出口垂直月表向下）羽流與月面月塵（月壤）的相互作用。

據(jù)此，將模型的參數(shù)設(shè)置為：基本流場(chǎng)為以5mm為半徑的底面、15mm為高的圓柱，圓柱的上底面作為噴流的入口，下底面作為月面。月壤顆粒密度設(shè)為3100，泊松比為0.42，彈性模量為10MPa，環(huán)境壓強(qiáng)為0.01Pa;噴管的出口（亦即流場(chǎng)的入口）速度為1000m/s（因?yàn)閲姽艹隹趬毫υ?.02MPa左右，而環(huán)境壓強(qiáng)僅為0.01Pa，噴流在噴管內(nèi)處為極度膨脹不足的狀態(tài)，所以達(dá)到這樣大的速度也是合理的）。重力加速度為1.62。

邊界條件設(shè)置如下：流場(chǎng)入口為壓力入口邊界條件;出口則為壓力出口邊界條件;在固體壁面上，氣相流場(chǎng)采用無(wú)滑移固壁邊界條件，顆粒與壁面之間采用Hertz Mindlin（No Slip）碰撞模型。

3.2 數(shù)值計(jì)算方法

采用FLUENT與EDEM軟件耦合求解氣固兩相流場(chǎng)，氣流流場(chǎng)采用FLUENT軟件，選用不可壓縮流體、 紊流模型計(jì)算;在EDEM軟件中采用Lagrangian坐標(biāo)方法分析顆粒在彎管中運(yùn)動(dòng)和分布規(guī)律。首先由FLUENT得到原始?xì)饬鲌?chǎng)結(jié)果，并將流場(chǎng)參數(shù)傳遞給EDEM;由EDEM中DEM模塊計(jì)算在該流場(chǎng)下顆粒受力的大小，然后通過(guò)內(nèi)部迭代計(jì)算顆粒下一時(shí)間步長(zhǎng)運(yùn)動(dòng)參數(shù)，同時(shí)更新顆粒運(yùn)動(dòng)軌跡;最后將EDEM計(jì)算結(jié)果返回FLUENT進(jìn)行下一時(shí)間步長(zhǎng)迭代。

采用Lagrangian方法計(jì)算顆粒運(yùn)動(dòng)軌跡，由Rayleigh數(shù)控制迭代步長(zhǎng)。Rayleigh數(shù)由下面公式給出：

式中，R是顆粒半徑，ρ是顆粒密度，G是剪切模量，υ是泊松比。則顆粒運(yùn)動(dòng)的迭代時(shí)間步長(zhǎng)為：

FLUENT時(shí)間步長(zhǎng)應(yīng)為EDEM軟件時(shí)間步長(zhǎng)的1/10～1/100，由此來(lái)確定流場(chǎng)的網(wǎng)格尺度，故FLUENT和EDEM網(wǎng)格不同、時(shí)間尺度不同。

3.3計(jì)算結(jié)果及分析

經(jīng)過(guò)耦合計(jì)算之后，流場(chǎng)顯示收斂，保存計(jì)算結(jié)果。通過(guò)EDEM自帶的Analyst模塊，我們可以導(dǎo)出進(jìn)行加工處理后的數(shù)據(jù)，再對(duì)其進(jìn)行相關(guān)的分析討論。

計(jì)算流場(chǎng)的速度云圖如下所示：

月壤顆粒的各項(xiàng)特性隨仿真時(shí)間的推進(jìn)的變化圖如下：

為了加快仿真時(shí)間，在EDEM中的顆粒工廠只加入了1000個(gè)月壤顆粒進(jìn)行計(jì)算，由上圖可知，仿真只進(jìn)行到了0.3s左右，顆粒已經(jīng)幾乎全部被吹出了仿真區(qū)域，而且被吹起的月壤的平均速度大多能夠達(dá)到二三十米每秒，說(shuō)明激起的揚(yáng)塵足以對(duì)探測(cè)器上所攜帶的儀器產(chǎn)生干擾，會(huì)對(duì)探測(cè)器的安全及準(zhǔn)確降落造成影響，與Apollo15的錄像測(cè)量數(shù)據(jù)相吻合。

計(jì)算完成后，我們可以創(chuàng)建圖表來(lái)記錄指定的時(shí)間步長(zhǎng)或者一定時(shí)間段內(nèi)的仿真情況。以下便是不同時(shí)間點(diǎn)下一定速度范圍內(nèi)各月壤顆粒數(shù)目統(tǒng)計(jì)的直方圖，從圖中我們可以非常直觀地得到仿真過(guò)程中某一時(shí)間點(diǎn)下的不同速度的顆粒分布狀況，從而對(duì)我們的后續(xù)研究及觀察提供有利的參考依據(jù)。

以下為各個(gè)時(shí)間點(diǎn)下流場(chǎng)中不同位置的顆粒的數(shù)目分布圖，據(jù)此我們可以得出流場(chǎng)內(nèi)顆粒的分布密度，從而判定是否對(duì)探測(cè)器上的設(shè)備造成影響。

4 結(jié)論

相較于以往的兩相流計(jì)算仿真，F(xiàn)LUENT-EDEM耦合方法不需要像使用UDF添加顆粒相需要進(jìn)行編程那么麻煩，而且EDEM由于是采用表面網(wǎng)格來(lái)描述邊界表面，使其可以與FLUENT的流體網(wǎng)格的邊界表面元素的點(diǎn)對(duì)點(diǎn)耦合;加上耦合方案為雙向動(dòng)力耦合，令仿真效果更加真實(shí)可信。由于可以直接在軟件中手動(dòng)設(shè)置顆粒的相關(guān)參數(shù)，可以很方便地進(jìn)行各項(xiàng)屬性的更改，而且EDEM也提供了非球體顆粒的設(shè)置，可以讓仿真更加貼近實(shí)際情況。此外，EDEM中還加入了先進(jìn)的靜電模型來(lái)仿真帶電顆粒，這與月壤因太陽(yáng)風(fēng)的吹掃而帶電的情況相吻合。

參考文獻(xiàn)：

[1]倪晉仁，李振山.風(fēng)沙兩相流理論及其應(yīng)用[M].北京：科學(xué)出版社，2006.

[2]Oesterlé，B.and Bui Dinh，T.，Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers，Experiments in Fluids，25，pp.16-22，1998.

[3]Zhou，H.，F(xiàn)lamant，G.，Gauthier，D.，DEM-LES simulation of coal combustion in a bubbling fluidized bed.Part II：coal combustion at the particle level，Chemical Engineering Science，59 pp.4205-4215，2004.

[4]許峰，王彤，傅耀.彎管內(nèi)小尺度顆粒運(yùn)動(dòng)的數(shù)值模擬與分析.四川兵工學(xué)報(bào).2010.02