范 俊

(西北工業(yè)大學(xué)航空學(xué)院流體力學(xué)系，西安710072)

近20年里，格子Boltzmann方法(lattice Boltzmann method，LBM)在理論與應(yīng)用研究上取得了顯著的進(jìn)展，并逐漸成為相關(guān)領(lǐng)域的研究熱點(diǎn).與基于宏觀Navier－Stokes(N－S)方程組的傳統(tǒng)計(jì)算流體力學(xué)方法不同，格子Boltzmann方法是一種介觀方法［1－2］.

以動(dòng)理學(xué)模型為基礎(chǔ)的介觀方法，既具有微觀方法假設(shè)條件較少的特點(diǎn)，又具有宏觀方法不關(guān)心分子運(yùn)動(dòng)細(xì)節(jié)的優(yōu)勢(shì).因此在處理多尺度、多物理的復(fù)雜流動(dòng)問題中，介觀方法具有較大優(yōu)勢(shì)和潛力.格子Boltzmann方法是由格子氣自動(dòng)機(jī)(lattice gas automata，LGA)發(fā)展而來的，是一種不同于傳統(tǒng)數(shù)值方法的流體計(jì)算和建模方法，其主要優(yōu)點(diǎn)有:演化過程物理清晰，邊界條件容易處理，計(jì)算簡(jiǎn)單，容易編程實(shí)現(xiàn)，并且計(jì)算是局部的，具有良好的并行性和擴(kuò)展性，有利于大規(guī)模流動(dòng)問題的計(jì)算.該方法在宏觀上是離散方法，在微觀上是連續(xù)方法.在許多傳統(tǒng)數(shù)值方法難以勝任的領(lǐng)域，如多相多質(zhì)流、微/納米尺度流、多孔介質(zhì)流、非牛頓流體、磁流體、生物流體粒子懸浮流、湍流、化學(xué)反應(yīng)流、燃燒問題、晶體生長(zhǎng)等，格子Boltzmann方法都取得了成功的應(yīng)用，并揭示了多種復(fù)雜現(xiàn)象的機(jī)理，推動(dòng)了相關(guān)科學(xué)的發(fā)展［3］.

1988年，McNamara和Zanetti提出了第一篇關(guān)于格子Boltzmann方法的論文，出現(xiàn)了最早的格子Boltzmann模型，但仍沿用格子氣自動(dòng)機(jī)的碰撞方式，具有計(jì)算指數(shù)復(fù)雜性.1991到1992年，幾個(gè)研究小組分別獨(dú)立提出了一種更簡(jiǎn)單的模型，即流行的單松弛(Single Relaxation Time，SRT)模型，又被稱為 LBGK(lattice Bhatnagar－Gross－Krook)模型［4－5］.LBGK 模型極大地簡(jiǎn)化了計(jì)算量，并且在一定條件下可以還原出正確的Navier－Stokes方程組，有效地克服了格子氣自動(dòng)機(jī)的不足.同時(shí)期，法國(guó)學(xué)者D·d’Humières提出了更一般性的廣義多松弛時(shí)間(multiple－relaxation－time，MRT)模型［6］.2000 年，P·Lallemand和 Luo對(duì) MRT 模型做了細(xì)致的理論分析，研究表明其在物理原理、參數(shù)選取和數(shù)值穩(wěn)定性方面都有很大的優(yōu)勢(shì).

多相態(tài)流體是包含明顯分界面的流體系統(tǒng)，自然界和現(xiàn)實(shí)工程中的流體系統(tǒng)或多或少地包含了其它物質(zhì)，完全單相態(tài)單組份的流體是幾乎不存在的.多相態(tài)流動(dòng)是化工、能源、航空、航天、航海、生命科學(xué)等領(lǐng)域廣泛存在的流動(dòng)現(xiàn)象.在LBM發(fā)展的過程中，很多研究人員提出了各種基于LBM的應(yīng)用方案，有幾種多相態(tài)流體模型，分別是顏色模型［7］、偽勢(shì)模型［8－9］、自由能模型［10］、動(dòng)理學(xué)模型.其中的偽勢(shì)模型是Shan和Chen于1993年提出的刻畫粒子間相互作用力的LBM模型，利用一種偽勢(shì)函數(shù)來反映粒子間相互作用.當(dāng)選擇了正確的偽勢(shì)函數(shù)，粒子間相互作用力保證了多相態(tài)分離流動(dòng)的自動(dòng)發(fā)生.

對(duì)于傳統(tǒng)的LBM模型，Boltzmann方程中的外力項(xiàng)離散與引入是一個(gè)關(guān)注熱點(diǎn)與難題.對(duì)于多相態(tài)流動(dòng)，正確地引入與處理外力項(xiàng)是非常重要.因?yàn)檫@對(duì)數(shù)值精度和計(jì)算穩(wěn)定性有顯著的影響.本文采用了LBGK模型和MRT模型，結(jié)合偽勢(shì)模型和三種不同的外力模型，模擬多相態(tài)流體的分離流動(dòng)現(xiàn)象，并對(duì)比不同數(shù)值模型的粒子間相互作用強(qiáng)度范圍和偽速度情況.

1 數(shù)值方法

1.1 格子Boltzmann方法

完整的LBM模型通常由格子(即離散速度模型)、平衡態(tài)粒子分布函數(shù)和演化方程三部分組成.首先，LBM輸運(yùn)方程的基本演化形式為

其中:fα(x，t)是t時(shí)刻在x位置處以離散速度eα運(yùn)動(dòng)的粒子分布函數(shù)，表示的是該狀態(tài)下的粒子的統(tǒng)計(jì)學(xué)概率;a為外力引起的加速度，Ωα(f)為表示粒子間碰撞引起的變化量的碰撞算子.求解Boltzmann方程最主要的難點(diǎn)在于碰撞項(xiàng)Ωα(f)的處理，它是關(guān)于粒子分布函數(shù)的非線性積分.因此，通常的想法是考慮能否用簡(jiǎn)化近似形式代替原本的碰撞項(xiàng)，這是LBM模型分類的依據(jù).LBGK模型和MRT模型都采用線性碰撞算子處理，LBGK模型只有一個(gè)自由參數(shù)，與運(yùn)動(dòng)黏性系數(shù)相關(guān)，體黏性與運(yùn)動(dòng)黏性固定相關(guān);而MRT模型將LBM方程發(fā)展至矩空間中處理，形式更復(fù)雜，但可調(diào)節(jié)的自由參數(shù)更豐富，其中包括了體黏性系數(shù).兩種模型碰撞算子的細(xì)節(jié)可以在很多已發(fā)表的文獻(xiàn)中得到，松弛時(shí)間參數(shù)與運(yùn)動(dòng)黏性系數(shù)的關(guān)系式為

其中:τ是松弛時(shí)間參數(shù)，v為運(yùn)動(dòng)黏性系數(shù)，δt為時(shí)間步長(zhǎng)，cs是與格子模型有關(guān)的格子聲速.

圖1 D2Q9不可壓格子模型離散速度圖

1992年，Qian等人提出了目前流行的DdQm系列格子模型［11］，d代表的是空間維數(shù)，m代表的是離散速度個(gè)數(shù)，對(duì)應(yīng)的平衡態(tài)粒子分布函數(shù)形式如下

其中:ωα是與格子模型有關(guān)的權(quán)系數(shù).本文采用的D2Q9是最常用的二維不可壓格子模型，如圖1所示，定義式為

其中:c=δx/δt，δx為格子步長(zhǎng).對(duì)于 D2Q9 格子模型，權(quán)系數(shù) ωα的設(shè)置為 w0=4/9，w1，2，3，4=1/9，=1/36，格子聲速 cs為格點(diǎn)處的宏觀物理量的計(jì)算，可通過將粒子分布函數(shù)求離散速度矩得到，不考慮外力作用的情況下，

LBM模型三要素中的演化方程可分解為兩個(gè)步驟進(jìn)行計(jì)算:

其中:式(8)是碰撞步，式(9)是遷移步.從上式可以看出，碰撞步驟的計(jì)算是在當(dāng)前格點(diǎn)處局部進(jìn)行的，遷移步驟只與相鄰的格點(diǎn)有關(guān)，這體現(xiàn)出LBM良好的計(jì)算局部性和優(yōu)秀的并行計(jì)算適用性.

在LBM模型中，外力項(xiàng) a·▽vfα(x，t)無法直接應(yīng)用求解.在Shan和Chen的原始模型中，通過改變平衡態(tài)粒子分布函數(shù)的宏觀速度實(shí)現(xiàn)外力作用的影響，平衡態(tài)宏觀速度為

從式(10)可以發(fā)現(xiàn)，相同的粒子間相互作用力下，隨著松弛時(shí)間參數(shù)的改變，平衡態(tài)宏觀速度是不一樣的，這樣的影響將會(huì)在后面的算例結(jié)果與分析中得到驗(yàn)證.

郭照立等人對(duì)LBM方程中的外力項(xiàng)進(jìn)行研究并提出了一種二階矩處理模型，應(yīng)用到了MRT模型上［2］:

其中:M－1是MRT模型的逆變換矩陣，S是包含了多松弛時(shí)間參數(shù)的對(duì)角矩陣，^F是作用力的各個(gè)速度矩，定義為

Kupershtokh等人［12］提出了一種精度較高、穩(wěn)定性和通用性較優(yōu)秀的外力模型，也是通過平衡態(tài)粒子分布函數(shù)實(shí)現(xiàn)外力作用的影響.然而，他們的外力模型和郭照立等人的類似，在碰撞算子以外處理，由于是基于平衡態(tài)粒子分布函數(shù)，該外力模型可適用于不同的碰撞算子:

上文介紹到的三種外力模型，分別采取三條完全不同的思路以實(shí)現(xiàn)外力項(xiàng)在離散Boltzmann方程中的引入，然而在外力作用的影響下，三種外力模型的流體宏觀速度同時(shí)都是取為碰撞前速度與碰撞后速度的平均值:

1.2 偽勢(shì)模型

Shan和Chen提出的偽勢(shì)模型，又被簡(jiǎn)稱為SC模型.他們提出了一種偽勢(shì)函數(shù)Ψ(x，t)來表示粒子間相互作用，偽勢(shì)函數(shù)與當(dāng)?shù)孛芏扔嘘P(guān)，而粒子間相互作用力是基于偽勢(shì)函數(shù)求得:

其中:G是一個(gè)表示粒子間相互作用強(qiáng)度的參數(shù).從公式(16)可以發(fā)現(xiàn)，粒子間相互作用力是對(duì)離格點(diǎn)處最近和較近位置上的粒子相互作用勢(shì)求和.本文采用的偽勢(shì)函數(shù)為Shan和Chen所提出的最原始方程:

其中:ρ0是一個(gè)自由參數(shù).

代入式(16)和(17)，流體系統(tǒng)會(huì)得到一個(gè)非理想狀態(tài)方程:

2 數(shù)值結(jié)果與分析

偽勢(shì)模型的一個(gè)優(yōu)點(diǎn)在于多相態(tài)流體的分離流動(dòng)自動(dòng)發(fā)生，不需要其他傳統(tǒng) CFD方法(如VOF、Level Set)所用到的動(dòng)態(tài)重構(gòu)、界面追蹤等技術(shù)，而且它們還難以捕捉跟蹤大量細(xì)小、分散的相界面.為了測(cè)試多相態(tài)流體的平衡特性，本文對(duì)多相態(tài)流體的分離流動(dòng)進(jìn)行數(shù)值模擬，分別選取了不同的相互作用強(qiáng)度G和運(yùn)動(dòng)黏性系數(shù)v.

整個(gè)數(shù)值模擬是在周期邊界和一套密度為128×128的笛卡爾網(wǎng)格下完成，流場(chǎng)平衡態(tài)的數(shù)據(jù)結(jié)果統(tǒng)一在500 000迭代時(shí)間步后獲取.初始化流場(chǎng)速度為0，流場(chǎng)各格點(diǎn)密度為一個(gè)固定值加上一個(gè)格點(diǎn)處隨機(jī)波動(dòng)，該隨機(jī)波動(dòng)是流體由于周圍相互作用力的不平衡而自動(dòng)發(fā)生分離流動(dòng)的原因.對(duì)于MRT模型，除了與運(yùn)動(dòng)黏性系數(shù)相關(guān)的松弛時(shí)間參數(shù)，還有多個(gè)自由松弛時(shí)間參數(shù)需要確定，本文對(duì)它們的定義如下［13］.

sα=(1，0.7，1.5，1，0.8，1，0.8，1/τ，1/τ).

圖2所展示的是粒子間相互作用強(qiáng)度為－5時(shí)，兩相態(tài)流體的分離流動(dòng)過程的密度分布截圖，時(shí)間步分別等于 200、800、1 600、1 800、4 000、20 000.在圖2(A)可以清晰地看到前期大量細(xì)小、分散相界面的存在，隨著時(shí)間步的推進(jìn)，液滴不停地融合，最后形成一個(gè)穩(wěn)定平衡的大液滴.

圖2 兩相態(tài)流體分離過程，粒子間相互作用強(qiáng)度

如上文所提到，LBGK模型的體黏性與運(yùn)動(dòng)黏性固定相關(guān)，兩者相等，由于穩(wěn)定性問題，這限制了LBGK模型在運(yùn)動(dòng)黏性較小(或雷諾數(shù)較高)的應(yīng)用.從圖3中可以看到，即使固定了運(yùn)動(dòng)黏性系數(shù)，MRT模型對(duì)粒子間相互作用強(qiáng)度的適用范圍比LBGK模型的更廣，這意味著更大的密度比，在圖4中得到了驗(yàn)證.我們?cè)鴩L試在LBGK模型上采用郭照立和A·L·Kupershtokh的外力模型，效果仍不理想，MRT模型的上限是5.75，而LBGK模型的上限是 5.25.

偽速度，是數(shù)值模擬中多相態(tài)流動(dòng)平衡狀態(tài)下的最大流場(chǎng)速度，一般都表現(xiàn)為相界面附近的旋渦流動(dòng)，這是一種不健康的非物理現(xiàn)象.伴隨著偽密度振蕩的發(fā)生，過大的偽速度會(huì)降低數(shù)值精度和影響計(jì)算穩(wěn)定性，極大地限制了模型的應(yīng)用范圍.因此，偽速度經(jīng)常作為多相態(tài)流動(dòng)數(shù)值模擬的評(píng)判標(biāo)準(zhǔn)和對(duì)比依據(jù).圖5展示的是兩種外力模型的偽速度隨粒子間相互作用強(qiáng)度變化的曲線，兩種外力模型的結(jié)果差異不大，隨著密度比的增大而單調(diào)遞增.

圖3 不同模型下的液相態(tài)和氣相態(tài)密度隨粒子間相互作用強(qiáng)度變化曲線

圖4 不同模型下的液相態(tài)和氣相態(tài)之間密度比隨粒子間相互作用強(qiáng)度變化曲線

圖5 不同外力模型的偽速度隨粒子間相互作用強(qiáng)度變化曲線

SC模型的外力模型對(duì)粒子間相互作用力的引入受松弛時(shí)間參數(shù)的影響，即不同的運(yùn)動(dòng)黏性系數(shù)下，多相態(tài)流體平衡特性可能會(huì)不相同.因此，本文進(jìn)行了一系列的測(cè)試.三種數(shù)值架構(gòu)的粒子間相互作用力都是按照式(16)計(jì)算，在圖6中，LBGK模型的平衡態(tài)密度特性隨著運(yùn)動(dòng)黏性系數(shù)的變化而改變，這一種現(xiàn)象在一些已發(fā)表的文獻(xiàn)中曾被提到［13－14］，在圖7中的平衡態(tài)密度比曲線的變化表現(xiàn)得更明顯.由于外力的引入方式不合理而導(dǎo)致的這一種現(xiàn)象是不受歡迎的，這不利于應(yīng)用過程中對(duì)平衡態(tài)特性的設(shè)定.另一方面，從理論上，通過多尺度展開(Chapman－Enskog分析)，加入了外力模型的LBM方程可還原為傳統(tǒng)的宏觀N－S方程組:

圖6 不同模型下的液相態(tài)和氣相態(tài)密度隨運(yùn)動(dòng)黏性系數(shù)變化曲線

其中:Rρ和Rv分別是質(zhì)量方程和動(dòng)量方程偏差量.SC外力模型存在因引入外力而導(dǎo)致的偏差量，這一方法主要適用于外力為常數(shù)的情況，而郭照立和A·L·Kupershtokh的偏差量為0，完全滿足宏觀方程的還原.

圖7 不同模型下的液相態(tài)和氣相態(tài)之間密度比隨運(yùn)動(dòng)黏性系數(shù)變化曲線

3 結(jié)語(yǔ)

本文采用LBGK模型和MRT模型，結(jié)合三種外力模型對(duì)多相態(tài)流體的分離流動(dòng)進(jìn)行數(shù)值模擬.不同的粒子間相互作用強(qiáng)度和運(yùn)動(dòng)黏性系數(shù)下的測(cè)試，證明了LBGK模型對(duì)多相態(tài)流動(dòng)模擬在數(shù)值精度與精算穩(wěn)定性上的局限;MRT模型能夠模擬更大密度比的流動(dòng);同時(shí)更合理的外力模型能夠改善MRT模型在多相態(tài)流動(dòng)的平衡態(tài)特性，避免受到運(yùn)動(dòng)黏性系數(shù)變化的影響，可以更方便地設(shè)定多相態(tài)流體的平衡特性.最后，這兩種不同的外力模型的數(shù)值結(jié)果沒有表現(xiàn)出明顯的優(yōu)劣差異，也都能成功地還原到對(duì)應(yīng)的宏觀N－S方程.

［1］何雅玲，王 勇，李 慶.格子Boltzmann方法的理論及應(yīng)用［M］.北京:科學(xué)出版社，2009.

［2］郭照立，鄭楚光.格子Boltzmann方法的原理及應(yīng)用［M］.北京:科學(xué)出版社，2009.

［3］AIDUN C K，CLAUSEN J R.Lattice－Boltzmann method for complex flows［J］.Annual Review of Fluid Mechanics，2010，42:439–472.

［4］CHEN S，CHEN H，MARTINEZ D，etal.Lattice Boltzmann model for simulation of magneto － hydrodynamics［J］.Physical Review Letters，1991，67(27):3776–3779.

［5］QIAN Y H，D’HUMI RES D，LALLEMAND P.Lattice BGK models for Navier － Stokes equation［J］.Europhysics Letters，1992，17(6):479–484.

［6］D’HUMI RESD.Generalized lattice Boltzmann equations［J］.Rarefied Gas Dynamics:Theory and Simulations，1994，159:450–458.

［7］GUNSTENSEN A K，ROTHMAN D H，ZALESKI S，etal.Lat-tice Boltzmann model of immiscible fluids［J］.Physical Review A，1991，43(8):4320.

［8］SHAN X W，CHEN H D.Lattice Boltzmann model for simulating flows with multiple phases and components［J］.Physical Review E，1993，47(3):1815.

［9］SHAN X W，CHEN H D.Simulation of nonideal gases and liquid － gas phase transitions by the lattice Boltzmann equation［J］.Physical Review E，1994，49(4):2941.

［10］SWIFT M R，ORLANDINI E，OSBORN W R，etal.Lattice Boltzmann simulations of liquid－gas and binary fluid systems［J］.Physical Review E，1996，54(5):5041.

［11］QIAN Y H，D’HUMI RESD，LALLEMAND P.Lattice BGK models for Navier－ Stokes equation［J］.Europhysics Letters，1992，17(6):479–484.

［12］KUPERSHTOKH A L，MEDVEDEV D A，KARPOV D I.On equations of state in a lattice Boltzmann method［J］.Computers＆ Mathematics with Applications，2009，58(5):965－974.

［13］YU Z，F(xiàn)AN L S.Multirelaxation－time interaction－potential－based lattice Boltzmann model for two－ phase flow［J］.Physical Review E，2010，82(4):046708.

［14］HE X，SHAN X，DOOLEN G D.Discrete Boltzmann equation model for nonideal gases［J］.Physical Review E，1998，57(1):13－16.