CFD-EDEM耦合方法在氣固兩相流研究中的應(yīng)用

李慈應(yīng)　李姝霖

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

【關(guān)鍵詞】兩相流;月壤;耦合;EDEM

引言

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

1 概述

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

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

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

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

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

2 CFD-EDEM耦合方法

2.1 EDEM軟件簡介

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

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

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

Analyst是EDEM的后處理模塊，含有很多工具對計(jì)算結(jié)果進(jìn)行加工處理，方便地獲得顆粒的速度、能量分布，設(shè)備的受力情況等，并可利用統(tǒng)計(jì)方法對各物理量間的相關(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ì)算時，EDEM的耦合模塊內(nèi)置于Fluent軟件中，兩者無縫連接;2）EDEM能夠自動讀取Fluent中的網(wǎng)格;3）EDEM可以使用Fluent中離散相模型和歐拉多項(xiàng)流模型進(jìn)行多相流耦合;4）采用EDEM直接對顆粒運(yùn)動進(jìn)行模擬;5）完全雙向動力耦合。

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

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

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

2.3 CFD-EDEM耦合理論

2.3.1 CFD-EDEM升力模型

2.3.1.1 Saffman升力模型

在邊界層，剪切層內(nèi)，或者是流體發(fā)生剪切時，會出現(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)可以表示為下式：

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

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

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

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

2.3.2 CFD-EDEM阻力模型

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

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

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

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

3.1 計(jì)算模型建立

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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