基于格子Boltzmann方法的多相態(tài)流動模型研究

范 俊

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

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

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

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

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

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

1 數(shù)值方法

1.1 格子Boltzmann方法

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

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

其中:τ是松弛時間參數(shù)，v為運動黏性系數(shù)，δt為時間步長，cs是與格子模型有關的格子聲速.

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

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

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

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

LBM模型三要素中的演化方程可分解為兩個步驟進行計算:

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

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

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

郭照立等人對LBM方程中的外力項進行研究并提出了一種二階矩處理模型，應用到了MRT模型上［2］:

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

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

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

1.2 偽勢模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)語

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

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

［2］郭照立，鄭楚光.格子Boltzmann方法的原理及應用［M］.北京:科學出版社，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.