多相流系統(tǒng)的離散玻爾茲曼研究進(jìn)展

許愛(ài)國(guó)，陳 杰，宋家輝，陳大偉，陳志華

(1. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100088；2. 北京大學(xué) 應(yīng)用物理與技術(shù)研究中心和高能量密度物理數(shù)值模擬教育部重點(diǎn)實(shí)驗(yàn)室，北京 100871；3. 北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081；4. 南京理工大學(xué) 瞬態(tài)物理國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210094；5. 北京理工大學(xué) 宇航學(xué)院，北京 100081)

0 引 言

在自然界和工程技術(shù)領(lǐng)域存在著大量形式各異的多相復(fù)雜流動(dòng)，例如超新星爆發(fā)、巖土礦石開(kāi)采中的多孔介質(zhì)流動(dòng)、石油開(kāi)采中原油在運(yùn)輸管道中的運(yùn)輸、航空發(fā)動(dòng)機(jī)中的氣液燃燒、武器物理、慣性約束核聚變，等等。多相復(fù)雜流動(dòng)研究具有重要的基礎(chǔ)與現(xiàn)實(shí)意義。

多相復(fù)雜流動(dòng)系統(tǒng)研究，需要不同層次、不同視角的方法和認(rèn)識(shí)。目前，在微觀、介觀和宏觀層面都有相應(yīng)的描述方法。在微觀層面，主要建模和模擬方法是分子動(dòng)力學(xué)（Molecular Dynamics, MD）[1-8]，其中分子間作用勢(shì)的構(gòu)建是模型構(gòu)建的關(guān)鍵，分子間作用勢(shì)有效半徑的選取是模擬過(guò)程中的關(guān)鍵。有效半徑的選取，以滿足需求且最小為最佳原則，一般通過(guò)模擬結(jié)果擬合已知的關(guān)鍵物性參數(shù)，根據(jù)模擬結(jié)果的合理性來(lái)判斷。MD的優(yōu)點(diǎn)是提供的信息完備，但缺點(diǎn)是能夠模擬的時(shí)間和空間尺度有限，目前主要用于一些分子、原子層次的機(jī)理研究[9]。宏觀模型主要是指基于納維-斯托克斯（Navier-Stokes，NS）方程的各類(lèi)多相流模型，模擬方法以傳統(tǒng)計(jì)算流體力學(xué)方法為主[10-11]，近年來(lái)，光滑粒子方法、格子玻爾茲曼方法（Lattice Boltzmann Method, LBM）等得到了快速發(fā)展和廣泛應(yīng)用[12-13]。多相流系統(tǒng)的介觀模型，一般是指基于非平衡統(tǒng)計(jì)物理學(xué)中動(dòng)理學(xué)理論構(gòu)建的動(dòng)理學(xué)模型。

經(jīng)過(guò)30年左右的迅猛發(fā)展，LBM已經(jīng)快速成長(zhǎng)為計(jì)算流體力學(xué)領(lǐng)域中的重要組成部分[14-26]。LBM方法與應(yīng)用方面的研究論文，其作者的學(xué)習(xí)、研究背景涵蓋領(lǐng)域極廣，因而同一個(gè)詞匯在不同的文獻(xiàn)里承載的內(nèi)涵也許并不相同。文獻(xiàn)中的很多LBM是以恢復(fù)或者求解流體力學(xué)方程為目的和功能的。但作為一種全新的數(shù)值解法，一大批非流體方程例如Benjamin方程、Korteweg-de Vries (KdV)-Burgers方程、Ginzburg-Landau方程、波動(dòng)方程、Poisson方程、Laplace方程、Lorenz方程、Fisher方程、Kuramoto-Sivashinsky方程、Richards方程、Schr?dinger方程等的LBM解法也引起了廣泛的興趣，獲得了廣泛研究。因?yàn)榛謴?fù)或者求解的方程并不是一般意義下的流體力學(xué)方程（組），所以這部分LBM所用的“玻爾茲曼方程”、“矩關(guān)系”也并不是非平衡統(tǒng)計(jì)物理學(xué)意義下的玻爾茲曼方程、矩關(guān)系。這部分LBM是純計(jì)算數(shù)學(xué)意義下的偏微分方程（Partial Differential Equation, PDE）數(shù)值解法。即便是恢復(fù)或者求解流體力學(xué)方程（組）的LBM中，有些LBM使用的“玻爾茲曼方程”和“矩關(guān)系”是與非平衡統(tǒng)計(jì)物理學(xué)理論一致的，有些則不一致或者部分不一致。這些具體處理方法的不同，展現(xiàn)的是LBM方法內(nèi)涵的多樣性；只要處理得當(dāng)，LBM便可滿足相應(yīng)的需求。

在復(fù)雜流體數(shù)值研究中，物理模型構(gòu)建和方程解法設(shè)計(jì)是不可或缺的兩個(gè)重要環(huán)節(jié)。在傳統(tǒng)計(jì)算流體力學(xué)中，物理建模和算法設(shè)計(jì)的界限是清晰的，后者為數(shù)值求解前者所獲方程（方程組）提供離散格式和步驟。近幾十年來(lái)，元胞自動(dòng)機(jī)-格子氣-格子玻爾茲曼系列方法的出現(xiàn)，讓二者在某些情形的界限變得模糊。因?yàn)閺牟煌囊暯强?，這些方法便具有不同的功能?；蛘哒f(shuō)，這些方法本身就可以朝著不同的方向發(fā)展。一方面，從物理學(xué)視角，長(zhǎng)期以來(lái)，元胞自動(dòng)機(jī)-格子氣模型就是一大類(lèi)復(fù)雜系統(tǒng)研究的理想化模型，統(tǒng)計(jì)物理與復(fù)雜系統(tǒng)的很多研究就是基于這類(lèi)理想化模型的。另一方面，近30年來(lái)，格子氣-格子玻爾茲曼又被作為一種全新的方程解法，在計(jì)算數(shù)學(xué)的規(guī)范下獲得了廣泛的研究。因?yàn)槟繕?biāo)不同，決定了構(gòu)建規(guī)則不可能完全相同，所以朝著不同目標(biāo)發(fā)展的這兩個(gè)方向成為這類(lèi)方法研究的重要內(nèi)容。在目前的絕大多數(shù)文獻(xiàn)中，LBM的功能是恢復(fù)或者求解相應(yīng)的偏微分方程，因而LBM在很大程度上成為“偏微分方程數(shù)值解法LBM”的簡(jiǎn)稱(chēng)。

在物理建模方面，又可粗略地分為兩種情形：A.傳統(tǒng)模型存在、合理、物理功能足夠且方便有效的情形；B.傳統(tǒng)模型不存在（以前未涉獵的新領(lǐng)域）、不再合理或物理功能不足的情形。元胞自動(dòng)機(jī)-格子氣-格子玻爾茲曼系列方法在情形A和情形B均適用。在情形A，這些方法又可視為求解傳統(tǒng)模型方程或方程組的一種全新方法（其求解思路與傳統(tǒng)解法完全不同）。因?yàn)閿?shù)值解法研究和物理建模研究目標(biāo)不同，所以需要遵循的規(guī)則不會(huì)完全相同；即使是從同一起點(diǎn)出發(fā)，即使有重疊，二者的發(fā)展軌跡也不會(huì)完全相同；二者時(shí)而近時(shí)而遠(yuǎn)，在相互啟發(fā)中發(fā)展。

在本文中，基于或借鑒動(dòng)理學(xué)理論的方法粗略地劃分為兩類(lèi)：方程解法設(shè)計(jì)和物理模型構(gòu)建方法，重點(diǎn)討論流體系統(tǒng)的物理建模，重點(diǎn)關(guān)注傳統(tǒng)模型描述不了或描述不好而MD又由于適用尺度受限無(wú)能為力的“介尺度”、兩難情形。這類(lèi)“介尺度”物理問(wèn)題的研究需求是離散玻爾茲曼方法（Discrete Boltzmann Method, DBM）產(chǎn)生的背景和驅(qū)動(dòng)力。在不產(chǎn)生歧義的情形，DBM又可視為離散玻爾茲曼模型（Discrete Boltzmann Model）、離散玻爾茲曼建模方法（Discrete Boltzmann Modeling method）的簡(jiǎn)稱(chēng)。

原則上，廣義地講，將玻爾茲曼方程在某些方面做些離散而做進(jìn)一步研究的方法，甚至基于玻爾茲曼方程發(fā)展起來(lái)的離散方法都可以顧名思義地稱(chēng)為離散玻爾茲曼方法（DBM）。這些DBM方法可以是理論物理中的建模方法，也可以是計(jì)算數(shù)學(xué)中的方程解法。由于可以根據(jù)理論或模擬研究需要，分別將時(shí)間、空間和粒子速度這三個(gè)自變量之一、之二或之三進(jìn)行離散，所以DBM的含義可以很廣。有些DBM是包含離散速度圖像的（例如LBM），也有些并不包含（例如，計(jì)算流體力學(xué)中的有些動(dòng)理學(xué)方法充分借助麥克斯韋分布函數(shù)動(dòng)理學(xué)矩的解析解，并不使用離散速度圖像，其中的離散指的可以是時(shí)間和空間的離散）。為方便描述，在本文中，如果沒(méi)有特殊說(shuō)明，則DBM特指針對(duì)上述“介尺度”、“兩難”情形而構(gòu)建的一類(lèi)相對(duì)具體的理論模型或方法，其中的離散指的是粒子速度空間的離散。

具體來(lái)說(shuō)，本文要重點(diǎn)介紹的離散玻爾茲曼建模，是非平衡統(tǒng)計(jì)物理學(xué)粗?；７椒ㄔ诹黧w力學(xué)領(lǐng)域的具體應(yīng)用，是理論物理范疇下的模型構(gòu)建方法。它根據(jù)研究需求，抓主要矛盾，選取一個(gè)視角，研究系統(tǒng)的一組動(dòng)理學(xué)性質(zhì)，因而它要求描述這組性質(zhì)的動(dòng)理學(xué)矩其計(jì)算結(jié)果不能因?yàn)槟Ｐ秃?jiǎn)化而改變。除了分布函數(shù)的守恒矩（質(zhì)量、動(dòng)量和能量），DBM同時(shí)關(guān)注部分關(guān)系最密切的非守恒矩的時(shí)空演化，從一個(gè)更寬的視角來(lái)考察和描述系統(tǒng)。除了為實(shí)現(xiàn)模擬而進(jìn)行的抓主要矛盾和粗?；＃珼BM同時(shí)關(guān)注模擬之后的復(fù)雜物理場(chǎng)分析（復(fù)雜物理場(chǎng)分析，也需要通過(guò)建模來(lái)提取更多有價(jià)值的信息），它本身自帶一套復(fù)雜物理場(chǎng)分析方法或技術(shù)[27-28]。

廣義的DBM包含LBM，LBM是其中使用離散速度的一類(lèi)。本文要重點(diǎn)介紹的這類(lèi)具體的DBM也可以視為廣義的LBM系列中物理建模這一類(lèi)的進(jìn)一步發(fā)展[29]。文獻(xiàn)[29]指出，在模型構(gòu)建與非平衡統(tǒng)計(jì)物理學(xué)理論一致的情形，可以借助 (f-feq)的非守恒矩來(lái)描述系統(tǒng)偏離平衡的方式和幅度，提取非平衡信息和效應(yīng)。隨后，在 (f-feq)非守恒矩張開(kāi)的相空間及其子空間中，借助到原點(diǎn)的距離提出相應(yīng)的“非平衡強(qiáng)度”概念；借助兩點(diǎn)間距離d來(lái)描述兩個(gè)非平衡狀態(tài)的差異，借助其倒數(shù)提出“非平衡狀態(tài)相似度”的概念（S=1/d）；借助一段時(shí)間內(nèi)兩點(diǎn)間距離的平均值dˉ來(lái)描述兩個(gè)動(dòng)理學(xué)過(guò)程之間的差異，借助其倒數(shù)提出“動(dòng)理學(xué)過(guò)程相似度”的概念（Sp=1/dˉ）。從而，使得一些以前無(wú)法獲取的非平衡信息得以分層次、定量化研究[27-28,30-31]。鑒于目前LBM在很大程度上已經(jīng)成為“偏微分方程解法LBM”的簡(jiǎn)稱(chēng)，所以在本文中，如果沒(méi)有特殊說(shuō)明，在不引起誤解的情形，我們將沿用這一簡(jiǎn)稱(chēng)。在這些約定下，DBM和LBM各自的內(nèi)涵是具體的、清晰的。

本文結(jié)構(gòu)如下：第1節(jié)，介紹流體系統(tǒng)的微觀、介觀和宏觀三種建模方法；第2節(jié)，介紹從格子氣模型到離散玻爾茲曼方法的發(fā)展歷程；第3節(jié)，給出離散玻爾茲曼方法的理論框架，然后分別介紹含外場(chǎng)力情形、含分子間作用勢(shì)情形、含化學(xué)反應(yīng)情形和等離子體情形的DBM建模；第4節(jié)，介紹基于DBM模擬的多相流機(jī)理研究，包括在流體不穩(wěn)定性研究、相分離研究、化學(xué)反應(yīng)流研究等方面的進(jìn)展。第5節(jié)，給出本文的小結(jié)與說(shuō)明。

1 流體系統(tǒng)的微觀、介觀與宏觀建模

物質(zhì)世界是無(wú)限可分的，微觀、介觀、宏觀的界定是相對(duì)的。對(duì)于流體系統(tǒng)來(lái)說(shuō)，微觀描述一般是指基于MD的描述。在MD中，分子間相互作用勢(shì)的建立是模型構(gòu)建的關(guān)鍵，有效作用半徑的截取是模擬計(jì)算的關(guān)鍵。宏觀描述一般是指基于歐拉（Euler）方程、NS方程等連續(xù)介質(zhì)力學(xué)建模的描述。在這個(gè)層面上，人們關(guān)注的系統(tǒng)行為通過(guò)密度、流速、溫度、壓強(qiáng)、應(yīng)力、熱流等物理量表征，其控制方程是代表質(zhì)量、動(dòng)量和能量守恒的流體演化方程組。其本構(gòu)關(guān)系往往是根據(jù)經(jīng)驗(yàn)或唯象給出的。因?yàn)榉瞧胶饨y(tǒng)計(jì)物理學(xué)可以聯(lián)系微觀和宏觀連續(xù)描述，所以經(jīng)常稱(chēng)為介觀描述。其中，使用比較多的是基于玻爾茲曼方程的描述。在這類(lèi)描述中，描述的起點(diǎn)是理想氣體模型，恩斯柯格方程可以看作是玻爾茲曼方程的推廣，可以根據(jù)具體系統(tǒng)引入合適的分子間作用力或勢(shì)。

對(duì)于多相流等復(fù)雜流體系統(tǒng)來(lái)說(shuō)，系統(tǒng)本身可能是宏觀尺度的，但其內(nèi)部存在大量的中間尺度的空間結(jié)構(gòu)和動(dòng)理學(xué)模式，這些結(jié)構(gòu)和模式的存在和演化極大地影響著系統(tǒng)的物理性能和功能。多相復(fù)雜流體系統(tǒng)內(nèi)部往往具有大量的界面，包括物質(zhì)界面和力學(xué)界面，系統(tǒng)內(nèi)部的受力和響應(yīng)過(guò)程非常復(fù)雜。對(duì)于這些復(fù)雜系統(tǒng)的研究，NS方程等宏觀流體模型只能對(duì)系統(tǒng)中大尺度的和緩變的行為進(jìn)行較好的描述，而對(duì)于一些銳利的界面（例如沖擊波和爆轟波等）和快變模式的描述，則不能滿足要求，這至少體現(xiàn)在兩個(gè)方面：（1）只考慮線性響應(yīng)（一階非平衡效應(yīng)）的NS方程給出的應(yīng)力和熱流幅度可能偏大或偏小[32-33]；甚至在某些情形，NS方程給出的應(yīng)力、熱流、速度連方向都是錯(cuò)誤的[34]；（2）從動(dòng)理學(xué)理論視角，銳利界面的精細(xì)物理結(jié)構(gòu)描述不僅需要（分布函數(shù)的）守恒矩，還需要部分關(guān)系最密切的非守恒矩，而NS方程描述的只是守恒矩及其演化[27]。同時(shí)，發(fā)展相對(duì)成熟、大家相對(duì)熟悉的MD和直接模擬蒙特卡羅（Direct Simulation Monte Carlo, DSMC）方法，理論上相對(duì)完備，但由于運(yùn)算量極大，對(duì)計(jì)算機(jī)內(nèi)存的要求極高，所以在絕大多數(shù)情況下，它們能夠模擬的系統(tǒng)大小和時(shí)間尺度都遠(yuǎn)遠(yuǎn)不能滿足需求。

到這里，我們面對(duì)的狀況是，當(dāng)系統(tǒng)的克努森（Knudsen,Kn）數(shù)①分子的平均自由程λ與平均分子間距l(xiāng)成正相關(guān)，所以克努森數(shù)Kn可視為重新標(biāo)度的平均分子間距，描述的是系統(tǒng)的離散程度，其倒數(shù)描述的是系統(tǒng)的連續(xù)程度。對(duì)于非平衡流動(dòng)，Kn又可視為兩個(gè)分子發(fā)生碰撞的平均時(shí)間間隔與以分子平均速度通過(guò)關(guān)注的宏觀特征尺度所需時(shí)間之比，因而可視為重新標(biāo)度的分子碰撞平均時(shí)間間隔。因其與熱力學(xué)弛豫時(shí)間成正相關(guān)，所以可進(jìn)一步視為重新標(biāo)度的熱力學(xué)弛豫時(shí)間。從這個(gè)意義上，Kn描述的是系統(tǒng)偏離熱力學(xué)平衡的程度。很小時(shí)，連續(xù)介質(zhì)假設(shè)合理程度很高，傳統(tǒng)流體力學(xué)模型可以很好地描述；當(dāng)Kn很大時(shí)，系統(tǒng)的離散性很高，在關(guān)注的尺度內(nèi)粒子數(shù)較少；少到一定程度，就可以使用MD或DSMC來(lái)進(jìn)行模擬。而當(dāng)Kn介于這兩種情形之間（介尺度情形）時(shí)，傳統(tǒng)連續(xù)模型不再合理，而MD和DSMC又由于適用的時(shí)空尺度受限而無(wú)能為力。同時(shí)，隨著進(jìn)入低壓稀疏、微介觀或者快變模式情形，Kn逐漸升高，系統(tǒng)偏離平衡程度增加，系統(tǒng)行為的復(fù)雜度急劇上升，因而使用分布函數(shù)非守恒矩對(duì)系統(tǒng)行為進(jìn)行描述的必要性在逐漸增加。DBM就是在這個(gè)背景下，應(yīng)這些物理問(wèn)題的研究需求而產(chǎn)生的理論模型[28]。因?yàn)檫@些以前知之甚少的非平衡行為蘊(yùn)含著大量尚待開(kāi)發(fā)的物理功能，所以隨著時(shí)間的推進(jìn)，使用分布函數(shù)非守恒矩描述系統(tǒng)行為的收益正在增加，這在后面介紹DBM應(yīng)用時(shí)將會(huì)看到。

2 粗?；＃簭母褡託饽Ｐ偷诫x散玻爾茲曼

元胞自動(dòng)機(jī)又稱(chēng)為格子氣（元胞）自動(dòng)機(jī)（Lattice Gas（Cellular）Automata, LG（C）A）或 格 子 氣 模 型（Lattice Gas Model，LGM）。LGM的使用在統(tǒng)計(jì)物理學(xué)的研究中有著悠久的歷史。1952年，李政道和楊振寧發(fā)表在《物理評(píng)論》的《Statistical Theory of Equations of State and Phase Transitions. II. Lattice Gas and Ising Model》[35]就是一個(gè)關(guān)于LGM的工作。

從20世紀(jì)60年代開(kāi)始，人們?cè)O(shè)想使用一種稱(chēng)之為“格子氣”或“元胞自動(dòng)機(jī)”的簡(jiǎn)單離散模型來(lái)作為連續(xù)流體系統(tǒng)的理想化模型。這類(lèi)模型通常由一個(gè)背景網(wǎng)格、分布在網(wǎng)格結(jié)點(diǎn)上的一群“虛擬粒子”和一套簡(jiǎn)單作用規(guī)則這三個(gè)部分構(gòu)成。其中，背景網(wǎng)格提供空間粗?；鴺?biāo)。分布在網(wǎng)格結(jié)點(diǎn)上的“虛擬粒子”具有相同的質(zhì)量，并且只具有少數(shù)和網(wǎng)格點(diǎn)的連接方式綁定的運(yùn)動(dòng)方向和速度大小。常用的簡(jiǎn)單規(guī)則為—在一個(gè)時(shí)間步長(zhǎng)內(nèi)，粒子只能從當(dāng)前格點(diǎn)沿格點(diǎn)的連線方向運(yùn)動(dòng)到相鄰格點(diǎn)上，這個(gè)過(guò)程稱(chēng)為“傳播”。當(dāng)來(lái)自不同方向的粒子在某個(gè)網(wǎng)格點(diǎn)相遇時(shí)，它們則發(fā)生“碰撞”。碰撞的設(shè)計(jì)要使得碰撞前后系統(tǒng)局域的質(zhì)量守恒、動(dòng)量守恒以及能量守恒。為了能夠模擬連續(xù)的流體系統(tǒng)，特別是讓在微觀層次上破缺的對(duì)稱(chēng)性在宏觀上得到恢復(fù)，格子氣模型的構(gòu)建必須滿足對(duì)稱(chēng)性約束條件。到20世紀(jì)80年代后期，LGM開(kāi)始得到快速發(fā)展[36-37]。其中，以Frisch-Hasslacher-Pomeau（FHP）模型為典型代表[38]。

物理圖像簡(jiǎn)單、直觀、易于編程等優(yōu)點(diǎn)，使得LGM在粗略模擬流體行為（特別是復(fù)雜流體行為）方面發(fā)揮著重要作用。NS等偏微分方程的數(shù)值解法和非平衡復(fù)雜流動(dòng)的粗粒化物理建模是LGM的兩個(gè)發(fā)展方向。這兩類(lèi)LGM目標(biāo)不同，所以構(gòu)建規(guī)則不同，使用的判據(jù)也不同。

LBM是在LGM的基礎(chǔ)上經(jīng)過(guò)幾個(gè)主要的演變發(fā)展而來(lái)的。LGM的簡(jiǎn)單離散規(guī)則使得它和以Euler、NS方程為代表的連續(xù)模型描述的流體行為之間出現(xiàn)差異。由于連續(xù)模型在其適用的范圍內(nèi)的可靠性是經(jīng)過(guò)大量的實(shí)踐檢驗(yàn)的，為了減小模擬的誤差，人們對(duì)LGM的合理構(gòu)造展開(kāi)了廣泛的研究。經(jīng)過(guò)引入局域平衡態(tài)分布函數(shù)、線性碰撞算符、單弛豫時(shí)間模型等，LGM逐漸發(fā)展成了現(xiàn)在文獻(xiàn)中普遍使用的LBM。后來(lái)人們發(fā)現(xiàn)，LBM可以看作是單弛豫時(shí)間線性化玻爾茲曼方程在時(shí)間、空間和粒子速度空間的一種巧妙的離散化形式。隨后，多弛豫時(shí)間模型也得到了廣泛的研究。根據(jù)對(duì)玻爾茲曼方程輸運(yùn)項(xiàng)所采用的不同的計(jì)算方式，LBM可以分為有限差分LBM、有限體積LBM、有限元LBM等等。相應(yīng)的，由LGM發(fā)展而來(lái)的、繼承了“傳播+碰撞”這一簡(jiǎn)單演化規(guī)則的LBM通常稱(chēng)為標(biāo)準(zhǔn)LBM[29]。不同研究背景和知識(shí)結(jié)構(gòu)的研究人員都結(jié)合自己的研究興趣和需求，對(duì)LBM進(jìn)行了廣泛的推廣和發(fā)展?？傮w來(lái)講，LBM也有著兩個(gè)不同的發(fā)展方向：一是NS等偏微分方程的全新的數(shù)值解法，二是非平衡復(fù)雜流動(dòng)的動(dòng)理學(xué)建模方法。

在現(xiàn)有文獻(xiàn)中，多數(shù)LBM的功能是從一個(gè)全新的角度求解流體方程或其他偏微分方程。“LBM”已經(jīng)被普遍接受為一種全新的偏微分方程的數(shù)值解法，并且往往是標(biāo)準(zhǔn)LBM的代稱(chēng)。作為流體力學(xué)偏微分方程數(shù)值解法的LBM，模擬結(jié)果在數(shù)值誤差范圍內(nèi)要與需要求解的流體方程相符。LBM數(shù)值解法并不局限于流體力學(xué)方程，只要能找到一個(gè)模擬中的可控量與Kn相對(duì)應(yīng)，再將分布函數(shù)和動(dòng)理學(xué)矩的概念根據(jù)需要求解的方程作合理的延伸（例如，使所有離散“分布函數(shù)”之和等于某個(gè)需要求解的量），就可以通過(guò)在格子“玻爾茲曼方程”中添加合適的修正項(xiàng)來(lái)求解一些其他的偏微分方程。但是，概念延伸后的“分布函數(shù)”和“動(dòng)理學(xué)矩”可能已經(jīng)不再具有統(tǒng)計(jì)力學(xué)中分布函數(shù)和動(dòng)理學(xué)矩的物理含義，概念延伸后的“玻爾茲曼方程”也可能不再是非平衡統(tǒng)計(jì)力學(xué)中用于描述非平衡行為的玻爾茲曼方程。它們只是形式上的“分布函數(shù)”、“動(dòng)理學(xué)矩”、“玻爾茲曼方程”，因而只是一種叫法上的延續(xù)。相對(duì)于作為非平衡復(fù)雜流動(dòng)的動(dòng)理學(xué)建模方法而言，作為偏微分方程數(shù)值解法的格子玻爾茲曼，在算法的設(shè)計(jì)上可以不拘泥于嚴(yán)格的物理對(duì)應(yīng)，具有更大的靈活性。

相對(duì)于數(shù)值解法這一方向來(lái)說(shuō)，LBM在微介觀與非平衡效應(yīng)和建模方面并沒(méi)有得到同樣快速的發(fā)展。在作為非平衡復(fù)雜流動(dòng)動(dòng)理學(xué)建模方法時(shí)，根據(jù)關(guān)注點(diǎn)不同，LBM又可以分為兩類(lèi)。一類(lèi)主要關(guān)注質(zhì)量、動(dòng)量和能量三個(gè)守恒動(dòng)理學(xué)矩及其演化；另一類(lèi)則兼顧三個(gè)守恒矩和部分與之密切相關(guān)的非守恒矩，其中非平衡狀態(tài)的描述和非平衡行為的研究是核心（例如，通過(guò)非守恒矩與其相應(yīng)的平衡態(tài)值的差異研究熵增等不可逆機(jī)制，研究當(dāng)前非平衡狀態(tài)對(duì)下一步流動(dòng)行為的影響等）。

為了與作為偏微分方程數(shù)值解法的LBM相區(qū)別，同時(shí)由于第二類(lèi)作為非平衡流動(dòng)動(dòng)理學(xué)建模方法的模型一般不再使用來(lái)自格子氣方法的“虛擬粒子”、“傳播+碰撞”圖像（在一個(gè)時(shí)間步內(nèi)，虛擬粒子由當(dāng)前格點(diǎn)，沿著網(wǎng)格線方向移動(dòng)到相鄰的格點(diǎn)上去，來(lái)自相鄰格點(diǎn)的虛擬粒子在當(dāng)前格點(diǎn)發(fā)生碰撞），“格子”稱(chēng)謂已經(jīng)失去了原有的與格子氣方法對(duì)應(yīng)的含義，但因仍然使用離散的玻爾茲曼方程，所以在近期的文獻(xiàn)中逐漸改稱(chēng)為離散玻爾茲曼模型或建模方法（Discrete Boltzmann Model/Modeling method），簡(jiǎn)稱(chēng)DBM。這里，對(duì)空間導(dǎo)數(shù)的離散格式和控制方程的時(shí)間積分方法不做特別約束，可以根據(jù)具體情況合理選取，“離散玻爾茲曼”中的“離散”指的是分子速度空間的離散[27-28]。

3 離散玻爾茲曼理論框架

3.1 玻爾茲曼方程簡(jiǎn)介

在統(tǒng)計(jì)物理中，對(duì)于N個(gè)力學(xué)性質(zhì)完全相同的粒子組成的系統(tǒng)，假設(shè)每個(gè)粒子的力學(xué)狀態(tài)由r個(gè)廣義坐標(biāo)qm(m=1,2,···,r)和r個(gè)共軛的廣義動(dòng)量pm來(lái)確 定，那 么 我 們 可 以 用 2Nr維 ?？?間 中 的 一 個(gè) 點(diǎn)(q1,q2,···,qNr;p1,p2,···,pNr)來(lái)描述系統(tǒng)的狀態(tài)。

假設(shè)系統(tǒng)宏觀量為B(x,t)，微觀力學(xué)量為b(q,p;x,t)。則宏觀系統(tǒng)力學(xué)函數(shù)的觀測(cè)值等于相應(yīng)微觀函數(shù)的系綜平均值：

其中，F(xiàn)(q,p)為相空間的分布函數(shù)。

對(duì)于保守力學(xué)體系，哈密頓量等于系統(tǒng)總能量，等于常數(shù)，即：

若使用統(tǒng)計(jì)力學(xué)的“薛定諤圖像”，即系統(tǒng)的力學(xué)函數(shù)（b(q,p;x) ）給定，則系統(tǒng)的演化B(x,t)由分布函數(shù)F(q,p;t)隨時(shí)間的變化來(lái)描述。因?yàn)橄嗫臻g中代表點(diǎn)的數(shù)目在運(yùn)動(dòng)過(guò)程中保持不變，所以分布函數(shù)F(q,p;t)的 物質(zhì)導(dǎo)數(shù)為0，可以得到：

由哈密頓正則方程：

即可得到非平衡統(tǒng)計(jì)力學(xué)的基本方程—?jiǎng)⒕S（Liouville）方程：

上式又經(jīng)常寫(xiě)為：

線性算法L定義為：

其 中 [···]P為 泊 松 括 號(hào)。對(duì) 于 兩 個(gè) 力 學(xué) 函 數(shù)b(q,p)和c(q,p)，泊松括號(hào)的定義為：

為便于描述，引入簡(jiǎn)記：

引入s約化分布函數(shù)fs(x1,x2,···,xs)(s≤N)：

其中s是 0到N之間的整數(shù)。s粒子約化分布函數(shù)可以理解為在某一時(shí)刻在相空間中s個(gè)點(diǎn)同時(shí)發(fā)現(xiàn)s個(gè)粒子的聯(lián)合概率。對(duì)于N個(gè)質(zhì)量為m的全同粒子組成的經(jīng)典系統(tǒng)，假設(shè)系統(tǒng)的體積是有限的。系統(tǒng)的哈密頓量可做如下分解：

劉維量和哈密頓量成線性關(guān)系，所以劉維算法可作類(lèi)似的分解：

劉維方程可寫(xiě)為：

根據(jù)系統(tǒng)內(nèi)粒子數(shù)守恒和粒子在分布函數(shù)中的對(duì)稱(chēng)性，得：

上式可進(jìn)一步寫(xiě)為：

或

這是一個(gè)確定約化分布函數(shù)的方程鏈（或譜系），因其作者的名字（Bogoliubov-Born-Green-Kirkwood-Yvon）而稱(chēng)為BBGKY譜系。由于在引入過(guò)程中沒(méi)有作任何的近似，所以BBGKY譜系與劉維方程完全等價(jià)。

BBGKY方程鏈的第一、第二兩個(gè)方程為：

其中

如果我們引入五個(gè)假設(shè)對(duì)研究系統(tǒng)進(jìn)行約束：

假設(shè)一：當(dāng)s≥3時(shí) ，fs≡0（即認(rèn)為，相對(duì)于2個(gè)粒子的聯(lián)合概率，3個(gè)及以上粒子的聯(lián)合概率小到可以忽略不計(jì)）。于是，方程（19）簡(jiǎn)化為：

假設(shè)二：三個(gè)及以上粒子間的相互作用可以忽略，兩個(gè)粒子間的相互作用只與它們之間的距離有關(guān)，即：

假設(shè)三：外場(chǎng)保守假設(shè)，即粒子的加速度：

假設(shè)四：分子混沌假設(shè)，即在同一時(shí)刻在x1處找到一個(gè)粒子而在x2處找到另外一個(gè)粒子的聯(lián)合概率f2， 簡(jiǎn)單地等于在x1處 找到一個(gè)粒子的概率與在x2處找到另外一個(gè)粒子的概率的乘積：

假設(shè)五：剛球模型和彈性碰撞假設(shè)，即分子之間的碰撞用剛球之間的彈性碰撞近似。

由式（18）、式（19），再根據(jù)角動(dòng)量守恒，碰撞前后粒子“瞄準(zhǔn)距離”相等，可推導(dǎo)得到：

其中，b為瞄準(zhǔn)距離，χ為散射角，即兩粒子碰撞前后相對(duì)動(dòng)量之間的夾角，p、p1和p′、分別表示兩粒子碰撞前后的動(dòng)量。式（26）即為著名的玻爾茲曼方程。

從BBGKY方程鏈出發(fā)，經(jīng)過(guò)一系列假設(shè)，我們得到了玻爾茲曼方程。下面，我們從物理圖像的角度對(duì)玻爾茲曼方程進(jìn)行比較直觀的解釋。

分子動(dòng)理學(xué)提出，物質(zhì)是由大量粒子組成的，物質(zhì)在宏觀上的性質(zhì)是由粒子的種類(lèi)和粒子運(yùn)動(dòng)決定的。對(duì)于氣體來(lái)說(shuō)，氣體的宏觀性質(zhì)主要受到分子與分子之間以及分子和壁面之間的碰撞的影響。分子間發(fā)生三體碰撞的概率遠(yuǎn)低于分子間發(fā)生二體碰撞的概率，所以我們暫且只考慮三體碰撞概率或效應(yīng)小到可以忽略的稀薄情形，這時(shí)只需考慮分子間的二體碰撞。

分子二體碰撞中，最簡(jiǎn)單的碰撞為彈性碰撞，即分子間的碰撞沒(méi)有發(fā)生平動(dòng)能和內(nèi)能之間的能量交換，碰撞過(guò)程中機(jī)械能守恒。分子對(duì)碰撞前后的速度的大小可由動(dòng)量守恒和能量守恒來(lái)確定，而要確定分子對(duì)碰撞后速度的方向，則需要根據(jù)引入的不同的分子間作用勢(shì)模型來(lái)確定分子對(duì)的瞄準(zhǔn)距離b和碰撞偏轉(zhuǎn)角χ。

瞄準(zhǔn)距離b定義為尚未受到分子間的作用力時(shí)兩分子相對(duì)運(yùn)動(dòng)軌道之間的距離。b越小，兩分子間的碰撞效果越明顯，b=0即對(duì)應(yīng)正碰撞。

圖1所示為質(zhì)心坐標(biāo)中二體碰撞的分子運(yùn)動(dòng)軌跡示意圖。其中，vr和分別為分子對(duì)碰撞前的相對(duì)速度和碰撞后的相對(duì)速度。由動(dòng)量守恒和能量守恒可以得到vr=。分子碰撞軌道偏轉(zhuǎn)角χ由分子間作用勢(shì)決定。vr和為分子對(duì)碰撞前后相對(duì)速度的大小。假設(shè)分子碰撞前的瞄準(zhǔn)距離為b，碰撞后的瞄準(zhǔn)距離為b*，則由角動(dòng)量守恒可得b=b*。

圖1 質(zhì)心坐標(biāo)系中二體碰撞的分子運(yùn)動(dòng)軌跡Fig. 1 The molecular trajectory of two-body collision in the centroid coordinate system.

通過(guò)引入不同的分子間作用勢(shì)，可以得到不同的分子碰撞的散射特征。圖2為三維情形下分子碰撞和散射示意圖。

圖2 三維分子碰撞截面和散射示意圖Fig. 2 Schematic of 3D molecular collision cross section and scattering

假設(shè)兩個(gè)分子的相對(duì)速度為vr，通過(guò)微分截面bdbdω的 分子在碰撞之后散射到附近的立體角dΩ內(nèi)，其中：

設(shè) σ為單位立體角所對(duì)應(yīng)的截面積，則有：

可得：

由式（28）和式（29）可得總的碰撞截面 σT為：

其中，偏轉(zhuǎn)角χ以及瞄準(zhǔn)距離b的積分值由分子間作用勢(shì)模型所決定。選取不同的分子間作用勢(shì)模型，可以得到不同的分子模型的碰撞截面。其中，ω為碰撞平面與某一參考平面的夾角。

對(duì)于三維空間中的氣體系統(tǒng)，約化的分布函數(shù)f1(x1,t) 常 取 為f(r,v,t)。 對(duì)f(r,v,t)在 速 度空間 積 分 即可得到t時(shí)刻在空間位置r處的分子數(shù)密度n，即：

分布函數(shù)在速度空間的一階矩和二階縮并矩則分別表示t時(shí)刻位置空間r處的動(dòng)量和能量：

其中u和T為系統(tǒng)的宏觀速度和溫度，m為氣體分子的質(zhì)量，k為玻爾茲曼常數(shù)。對(duì)于t時(shí)刻相空間 (r,v)處的分子，假設(shè)其經(jīng)過(guò) dt時(shí) 間后移動(dòng)到了(r+dr,v+adt)處，其中a為分子受到外力所產(chǎn)生的加速度。則如果不考慮分子之間的相互碰撞，t時(shí)刻在相空間 (r,v)附近相體積 drdv內(nèi) 的分子將全部遷移到 (r+dr,v+adt)附近的相體積 drdv內(nèi)：

對(duì)等式(34)的左端進(jìn)行泰勒展開(kāi)可以得到：

其中f為t時(shí) 刻在 (r,v)處的分子的速度分布函數(shù)，等式（35）即為無(wú)碰撞時(shí)分布函數(shù)f(r,v,t) 隨時(shí)間的演化。若考慮分子間碰撞的影響，則演化方程變?yōu)椋?/p>

等式右端的碰撞項(xiàng)表示由于分子間的碰撞導(dǎo)致的分布函數(shù)的變化。對(duì)于碰撞項(xiàng)的具體形式，我們做如下分析。假設(shè)在空間體積 dr中，速度為v的分子和速度為v1的分子發(fā)生碰撞，碰撞后兩分子的速度分別變?yōu)関*和速度分布函數(shù)分別由f和f1變 為f*和。設(shè)碰撞前兩分子之間的相對(duì)速度為vr，則相對(duì)于速度為v1的 分子，速度為v的 分子在 dt時(shí)間內(nèi)運(yùn)動(dòng)所掃過(guò)的體積為vrσdΩdt。 空間體積 dr中，速度為v1的分子數(shù)密度為f1dv1， 則在時(shí)間 dt內(nèi) ，一個(gè)速度為v的分子與速度為v1的 分子之間發(fā)生碰撞的次數(shù)為vrσf1dΩdv1dt。由相空間體積 drdv內(nèi) 速度為v的 分子數(shù)為fdrdv，可以得到在時(shí)間 dt內(nèi) ，相空間體積 drdv內(nèi) 速度為v的分子和速度為v1的 分子之間發(fā)生碰撞的次數(shù)為vrσff1dΩdvdv1drdt。這個(gè)碰撞過(guò)程使得速度為v的分子數(shù)減少，稱(chēng)為正過(guò)程。

類(lèi)似地，如果碰撞前分子速度分別為v*和，分子速度分布函數(shù)分別為f*和，碰撞后分子速度分別為v和v1， 分子速度分布函數(shù)分別為f和f1，則可以得到在時(shí)間 dt內(nèi)，兩種分子在相空間體積 drdv*內(nèi)發(fā)生碰撞的次數(shù)為這個(gè)碰撞過(guò)程使得速度為v的 分子數(shù)增加，稱(chēng)為逆過(guò)程。根據(jù)分子碰撞的 對(duì) 稱(chēng) 性 有又 有所以發(fā)生逆碰撞的次數(shù)可寫(xiě)為

所以，在時(shí)間dt內(nèi)相空間體積 drdv中，由于和速度為v1的 分子發(fā)生碰撞而導(dǎo)致的速度為v的 分子數(shù)目的變化為對(duì)速度空間v1積分，則可以得到在時(shí)間 dt內(nèi) 相空間體積 drdv中，由于分子間的碰撞而導(dǎo)致的速度為v的分子數(shù)目的變化：

將碰撞項(xiàng)代入式(36)，可以得到完整形式的玻爾茲曼方程：

3.2 Chapman-Enskog多尺度分析

Chapman-Enskog多尺度分析是Chapman和Enskog在1910年至1920年提出的一種多尺度分析技術(shù)[24]。從數(shù)學(xué)角度看，Chapman-Enskog多尺度分析實(shí)際上就是將分布函數(shù)f、分布函數(shù)的時(shí)間導(dǎo)數(shù)和空間導(dǎo)數(shù)都看作是Kn的函數(shù)，將它們?cè)贙n=0處做泰勒展開(kāi)，代入玻爾茲曼方程，然后分別求密度、動(dòng)量和能量三個(gè)動(dòng)理學(xué)矩，以此來(lái)獲得流體動(dòng)力學(xué)方程組。

我們以一維情形為例來(lái)說(shuō)明。在Kn較小時(shí)，將分布函數(shù)及其時(shí)間導(dǎo)數(shù)和空間導(dǎo)數(shù)在Kn=0處做泰勒展開(kāi)：

其中，

為麥克斯韋分布，

將展開(kāi)式（39）～（41）代入分布函數(shù)f的演化方程，分別在方程兩側(cè)求同樣的動(dòng)理學(xué)矩（密度矩、動(dòng)量矩和能量矩等）。令方程兩側(cè)Kn同階項(xiàng)的系數(shù)相等，便可獲得該近似下f(n)用f(n-1),f(n-2),···,f(0)表示，最終用f(0)表示的關(guān)系式，而f(0)就是麥克斯韋分布，為已知量。代入相應(yīng)的密度矩、動(dòng)量矩和能量矩演化方程，即可得到該近似下的流體力學(xué)方程組（對(duì)應(yīng)質(zhì)量守恒、動(dòng)量守恒和能量守恒）。

下面，對(duì)時(shí)間、空間變化率的多尺度展開(kāi)方法蘊(yùn)含的測(cè)量逐步細(xì)化的思路做些理解性分析。首先，由式（46）可以看到， ?/?tn描述的是在將觀測(cè)所用的時(shí)間單元尺度由tn-1減 小到tn后，在之前的基礎(chǔ)上進(jìn)一步多觀測(cè)到的更高頻運(yùn)動(dòng)信息。因?yàn)榛跁r(shí)間單元t0（即n=0時(shí) ）觀測(cè)到的時(shí)間變化率為0，所以可以想到，t0對(duì)應(yīng)的應(yīng)該是系統(tǒng)行為演化跨越的總時(shí)間。對(duì)空間變化率的多尺度展開(kāi)方法可以做類(lèi)似的理解，?/?xm描 述的是將觀測(cè)所用的空間單元尺度由xm-1減小到xm后，在之前基礎(chǔ)上進(jìn)一步多觀測(cè)到的更細(xì)微結(jié)構(gòu)信息。因?yàn)榛诳臻g單元x0（即m=0時(shí)）觀測(cè)到的空間變化率為0，所以可以想到，x0對(duì)應(yīng)的應(yīng)該是系統(tǒng)在空間x方向跨越的總尺度，即系統(tǒng)在x方向上的大小。在目前文獻(xiàn)中的Chapman-Enskog多尺度分析，空間變化率一般只取一個(gè)有效觀測(cè)單元尺度x1。

從擾動(dòng)論角度看，在Chapman-Enskog多尺度分析中，Kn對(duì)應(yīng)的是施加給系統(tǒng)的擾動(dòng)。Chapman-Enskog多尺度展開(kāi)收斂的情形對(duì)應(yīng)系統(tǒng)可以回到平衡態(tài)的情形；如果擾動(dòng)過(guò)強(qiáng)，導(dǎo)致泰勒展式發(fā)散，則意味著該擾動(dòng)有可能引發(fā)新的結(jié)構(gòu)或模式。

3.3 離散玻爾茲曼建模：理想氣體情形

非平衡統(tǒng)計(jì)力學(xué)是聯(lián)系微觀與宏觀的橋梁。玻爾茲曼方程使用單體分布函數(shù)描述系統(tǒng)動(dòng)理學(xué)行為隨時(shí)間的演化。盡管相對(duì)于劉維方程來(lái)說(shuō)，玻爾茲曼方程已經(jīng)是個(gè)高度粗?；哪Ｐ?，但其碰撞項(xiàng)仍然非常復(fù)雜，以至于在絕大多數(shù)情形下仍然不方便求解。為了進(jìn)行有效使用，不得不根據(jù)具體情形，繼續(xù)對(duì)其進(jìn)行簡(jiǎn)化（抓住主要矛盾，有所保，有所丟）。

與動(dòng)理學(xué)宏觀建模（從動(dòng)理學(xué)理論出發(fā)推導(dǎo)宏觀模型方程組）不同，DBM是一種根據(jù)系統(tǒng)性質(zhì)研究需求直接建模的方法。在這種建模方式中，只是根據(jù)系統(tǒng)性質(zhì)研究需求，借助某種方式（例如Chapman-Enskog多尺度分析的物理圖像、經(jīng)過(guò)驗(yàn)證的唯象理論等）快速確定建模過(guò)程中要保持不變的底層動(dòng)理學(xué)矩關(guān)系，原則上無(wú)需經(jīng)過(guò)繁瑣的理論推導(dǎo)獲得最終的宏觀流體動(dòng)力學(xué)方程組。DBM模型對(duì)應(yīng)的宏觀流體動(dòng)力學(xué)方程組可以是后知的結(jié)果，而不必是DBM建模與模擬的前提。同時(shí)，DBM自動(dòng)提供部分關(guān)系最密切的非守恒動(dòng)理學(xué)矩及其演化，自動(dòng)彌補(bǔ)宏觀流體動(dòng)力學(xué)方程組在非平衡行為和效應(yīng)描述方面的功能不足。在借助Chapman-Enskog多尺度分析物理圖像的情形，Chapman-Enskog多尺度分析理論是這類(lèi)建模思路合理的理論依據(jù)。

從玻爾茲曼方程到目前版本DBM的建立，包括三個(gè)步驟：（1）碰撞算符的簡(jiǎn)約化；（2）分子速度空間的離散化；（3）非平衡狀態(tài)描述與非平衡信息的提取。其中，前兩步是針對(duì)待研系統(tǒng)的粗?；锢斫?；第三步是針對(duì)復(fù)雜物理場(chǎng)的描述與粗?；畔⑻崛?，是DBM建模方法的目的和核心。

由于玻爾茲曼方程的碰撞項(xiàng)非常復(fù)雜，在絕大多數(shù)情況下，不得不對(duì)其進(jìn)行簡(jiǎn)化。簡(jiǎn)化的方法之一就是引入一個(gè)形式上的局域目標(biāo)分布函數(shù)f*，將原來(lái)的碰撞項(xiàng)寫(xiě)成一個(gè)線性化碰撞算符的形式 (f*-f)/τ。它的物理含義是—分子碰撞的效果是使得分布函數(shù)f趨于f*， 其快慢由弛豫時(shí)間 τ來(lái)控制。碰撞項(xiàng)的線性化是一個(gè)粗?；５倪^(guò)程，會(huì)使得一部分物理信息丟失，在這個(gè)過(guò)程中，需要遵循的原則是：所關(guān)心的物理特征量使用簡(jiǎn)化前和簡(jiǎn)化后的模型計(jì)算，所得的結(jié)果一致。根據(jù)f*所選取的形式不同，該線性化模型習(xí)慣上分別稱(chēng)為BGK模型[39]、橢圓統(tǒng)計(jì)BGK模型[40]、Shakhov模型[41]、Rykov模型[42]、Liu模型[43]等。

DBM借助離散速度，將原本連續(xù)、積分形式的動(dòng)理學(xué)矩轉(zhuǎn)化為求和進(jìn)行計(jì)算。由于任何一個(gè)分子都可以朝向任何一個(gè)方向運(yùn)動(dòng)，且其范圍都是 (-∞,+∞)，所以（時(shí)間和位置空間常用的）常規(guī)離散方式對(duì)分子的速度空間并不適用。也就是說(shuō)，這里，fi并不具有確定的物理含義，即并不代表速度為vi的概率。所以在分析系統(tǒng)時(shí)所使用的并不是fi的具體數(shù)值，而是分布函數(shù)f的動(dòng)理學(xué)矩。DBM建模要求—關(guān)心的動(dòng)理學(xué)矩轉(zhuǎn)換為求和進(jìn)行計(jì)算后，得到的結(jié)果必須相同。即：

其中m為分子質(zhì)量，n為分子數(shù)密度，u為流體宏觀速度，D為空間維數(shù)，k為玻爾茲曼常數(shù)。

由Chapman-Enskog多尺度分析，f動(dòng)理學(xué)矩的計(jì)算可以歸結(jié)為feq動(dòng)理學(xué)矩的計(jì)算。所以，在分子速度空間離散化時(shí)所要遵循的約束為：

離散速度的具體取法取決于要保的不變量、基本的守恒性和必要的對(duì)稱(chēng)性?？梢?jiàn)，DBM給出的是對(duì)離散速度的物理約束，并不包含具體的離散格式。

非平衡信息的提取和分析技術(shù)是DBM的目的和核心。在宏觀流體模型中，克努森數(shù)、黏性、熱傳導(dǎo)、密度、溫度、壓強(qiáng)等宏觀量的梯度等都是常用的非平衡程度表征量，都從各自的角度來(lái)描述系統(tǒng)的非平衡程度。但它們都是將相關(guān)信息高度濃縮的、平均化、粗粒化描述，很多物理上感興趣的具體的非平衡信息（例如：不同自由度上的內(nèi)能，黏性應(yīng)力，熱流或更高階動(dòng)理學(xué)矩）通常它們是看不到且無(wú)法直接研究的。DBM借助 (f-feq)的非守恒動(dòng)理學(xué)矩對(duì)復(fù)雜流動(dòng)系統(tǒng)非平衡行為進(jìn)行更加細(xì)致的描述。

用Mn=Mn(f)來(lái) 表示分布函數(shù)f關(guān)于分子速度v的n階 動(dòng) 理學(xué) 矩，令使 用Mm,n表 示由m階張量縮并成的n階張量。在系統(tǒng)趨于或離開(kāi)熱力學(xué)平衡態(tài)的過(guò)程中，只有密度、動(dòng)量和能量三個(gè)動(dòng)理學(xué)矩是守恒的，其余的動(dòng)理學(xué)矩都是非守恒的。所以，

描述的就是相應(yīng)視角下流體系統(tǒng)偏離熱力學(xué)平衡態(tài)的具體信息。動(dòng)理學(xué)矩Mn和非平衡特征量 Δn中同時(shí)包含了流體分子的平均行為（即流體動(dòng)力學(xué)行為）和純粹的熱漲落行為（即熱力學(xué)行為）。所以，稱(chēng)非平衡特征量 Δn為熱動(dòng)非平衡（Thermo-Hydrodynamic Non-Equilibrium, THNE）特征量。進(jìn)而，使用描述分布函數(shù)f關(guān)于分子漲落速度 (v-u) 的n階動(dòng)理學(xué)矩，即中心動(dòng)理學(xué)矩。令則由中心矩定義的非平衡特征量為：

這些非守恒矩（或非平衡特征量）張量皆由若干分量構(gòu)成，其中只有部分分量是獨(dú)立的。這些張量構(gòu)成一個(gè)集合，其中的獨(dú)立分量也構(gòu)成一個(gè)集合。這里可以使用非平衡特征量集合{Δn}或{}的獨(dú)立分量張開(kāi)一個(gè)高維相空間。在這個(gè)相空間中，坐標(biāo)原點(diǎn)對(duì)應(yīng)熱動(dòng)（或熱力學(xué)）平衡態(tài)，其余任何一點(diǎn)都對(duì)應(yīng)一個(gè)具體的熱動(dòng)（熱力學(xué)）非平衡態(tài)。圖3（a）展示的是由熱力學(xué)非平衡特征量的獨(dú)立分量張開(kāi)的熱力學(xué)非平衡相空間示意圖（以具有三個(gè)獨(dú)立分量為例）。通過(guò)圖3（a），可以清楚地觀測(cè)（研究）系統(tǒng)從一個(gè)熱力學(xué)非平衡態(tài)（Thermodynamic non-equilibrium state 1，狀態(tài)1）到另一個(gè)熱力學(xué)非平衡態(tài)（Thermodynamic non-equilibrium state 2，狀態(tài)2）的演化過(guò)程。

通過(guò)這些非平衡特征量，我們可以研究系統(tǒng)的熵增，進(jìn)而通過(guò)熵增研究物質(zhì)混合，研究系統(tǒng)內(nèi)不同非平衡行為特征之間的空間關(guān)聯(lián)、時(shí)間關(guān)聯(lián)、時(shí)空關(guān)聯(lián)、競(jìng)爭(zhēng)與協(xié)作[44]。在非平衡特征量張開(kāi)的相空間中，某點(diǎn)到原點(diǎn)的距離可以定義為該狀態(tài)的非平衡程度或強(qiáng)度。在圖3（b）中，線段D*的長(zhǎng)度描述的是狀態(tài)1的非平衡程度或強(qiáng)度。在這個(gè)描述下，只要在一個(gè)球面上，非平衡強(qiáng)度就相同；所以非平衡強(qiáng)度相同的非平衡狀態(tài)有無(wú)窮多。進(jìn)一步，如圖3（c）所示，該相空間中兩狀態(tài)點(diǎn)之間的距離d描述的是兩個(gè)狀態(tài)偏離熱力學(xué)平衡態(tài)的差異，其倒數(shù)（S=1/d）描述的是這兩個(gè)狀態(tài)偏離平衡態(tài)的相似度。如果兩點(diǎn)間距離d為零，即兩點(diǎn)重合，即代表在當(dāng)前粗?；枋鱿逻@兩個(gè)狀態(tài)偏離平衡的方式完全相同，因而相似度無(wú)窮高。如果這兩個(gè)狀態(tài)都隨時(shí)間演化，如圖3（d）所示，那么某段時(shí)間內(nèi)兩點(diǎn)間距離的平均值可用于表述這兩個(gè)動(dòng)理學(xué)過(guò)程的差異，其倒數(shù)（Sp=1/dˉ）則可定義為這兩個(gè)動(dòng)理學(xué)過(guò)程之間的相似度，等等[27-28,30-31]。

圖3 非平衡特征量張開(kāi)的非平衡相空間Fig. 3 Phase space opened by non-equilibrium characteristic quantities

3.4 離散玻爾茲曼建模：含外場(chǎng)力情形

當(dāng)流場(chǎng)中的外場(chǎng)力不可忽略時(shí)（例如重力場(chǎng)存在下的瑞利-泰勒不穩(wěn)定性（Rayleigh-Taylor instability,RTI）問(wèn)題，電場(chǎng)力、磁場(chǎng)力存在下的等離子體輸運(yùn)問(wèn)題等），需要考慮外場(chǎng)力對(duì)流體介質(zhì)的作用。包含外場(chǎng)力的玻爾茲曼方程具有如下形式：

引入BGK近似后成為：

如果分布函數(shù)偏離平衡部分即 (f-feq)在外力項(xiàng)中的效應(yīng)不顯著，則可以將玻爾茲曼方程外力項(xiàng)中的f用feq近似代替，完成對(duì)粒子速度的求導(dǎo)運(yùn)算后，再將粒子速度空間離散，便得到如下的DBM演化方程：

其中，fi(r,vi,t)是 離散分布函數(shù)，t為時(shí)間變量，vi為離散速度，下標(biāo)i=1,2,···,N為離散速度的序號(hào)，r為空間變量，a為外場(chǎng)力產(chǎn)生的加速度，u為流體宏觀流速，T為流體溫度，τ為弛豫時(shí)間，為對(duì)應(yīng)平衡態(tài)分布函數(shù)的離散形式。

3.5 離散玻爾茲曼建模：含分子間作用勢(shì)情形

3.5.1 基于范德瓦爾斯?fàn)顟B(tài)方程的DBM

作為一個(gè)粗粒化模型和諸多問(wèn)題的研究起點(diǎn)，玻爾茲曼方程描述的是理想氣體系統(tǒng)。對(duì)于理想氣體系統(tǒng)，在確定的溫度和壓強(qiáng)下，只能有一個(gè)密度，系統(tǒng)基態(tài)永遠(yuǎn)是單一態(tài)，是不可能發(fā)生相變的。因此，在涉及相變的多相流問(wèn)題研究中，一個(gè)關(guān)鍵的問(wèn)題就是分子間作用力的引入。

從微觀分子相互作用勢(shì)角度來(lái)看，分子間的相互作用可以通過(guò)成對(duì)的勢(shì) φ(r)來(lái)描述，作用勢(shì)的大小只依賴(lài)于兩分子之間的距離。當(dāng)分子之間的距離大于某一值 σ時(shí)，作用勢(shì)表現(xiàn)為吸引作用；當(dāng)分子間距離小于 σ時(shí)，作用勢(shì)表現(xiàn)為排斥作用[45]。這兩種形式可以統(tǒng)一包含在Lenard-Jone勢(shì)函數(shù)中：

ε和 σ是兩個(gè)待定參數(shù)，對(duì)于不同氣體分子具有不同的值。在經(jīng)典力學(xué)中，N個(gè)質(zhì)量為m的全同粒子組成的系統(tǒng)的哈密頓量為：

其 中pi表 示 第i個(gè) 分 子 的 動(dòng) 量 矢 量，rij表 示 第i個(gè) 分 子和第j個(gè)分子之間的距離， 〈i,j〉表示對(duì)所有分子對(duì)求和。對(duì)于正則系綜，N粒子的配分函數(shù)ZN可以寫(xiě)成[45]：

其中D為空間維數(shù)， β=1/(kT)，k為玻爾茲曼常數(shù)（在下面的討論中一般取為1），hˉ為約化的普朗克常數(shù)，λth為德布羅意波長(zhǎng)：

由配分函數(shù)，可得自由能F的表達(dá)式：

以及壓強(qiáng)P、內(nèi)能密度e和單位粒子熵s的表達(dá)式：

其中下標(biāo)表示求導(dǎo)時(shí)的不變量。

范德瓦爾斯（Van Der Waals，VDW）理論提出對(duì)式（55）中的Lenard-Jone勢(shì)進(jìn)行簡(jiǎn)化，簡(jiǎn)化后的N粒子系統(tǒng)配分函數(shù)可寫(xiě)成[45]：

其中n=N/V表示粒子數(shù)密度，式中用到了近似公式：

再由式（60）可得VDW氣體狀態(tài)方程：

或者寫(xiě)為：

其中b=v0,a=εv0， 比體積v=1/n。

將VDW氣體狀態(tài)方程和表面張力引入DBM，可以得到如下形式的離散分布函數(shù)演化方程[46]：

其 中fi為離 散速 度vi對(duì) 應(yīng)的 分布 函數(shù)，為 相 應(yīng) 的 局域平衡態(tài)分布函數(shù)。i為離散速度的編號(hào)。Ii則用來(lái)表征分子間的相互作用力的影響[47]：

A、B、C和Cq為由表面張力和狀態(tài)方程決定的四個(gè)待定參數(shù)，ci=vi-u為 第i個(gè)離散速度與系統(tǒng)宏觀速度的矢量差。

3.5.2 基于C-S狀態(tài)方程的DBM

為了更準(zhǔn)確地描述硬球間的相互作用，Carnhan和Starling在VDW狀態(tài)方程的基礎(chǔ)上，對(duì)斥力項(xiàng)進(jìn)行了修正，提出Carnhan-Starling（C-S）狀態(tài)方程[48]：

其中 η=bρ/4。a和b分別為表示分子間引力和斥力的參數(shù)。

基于C-S狀態(tài)方程的DBM同樣具有式（68）和式（69）所示的形式，但是與基于VDW狀態(tài)方程的DBM不同，基于C-S狀態(tài)方程的DBM中Ii的系數(shù)分別為：

ρ、u、T分別是局域密度、流速和溫度。

Λ為壓強(qiáng)張量，I是單位張量，K是表面張力系數(shù)，ζ是體黏性系數(shù)。

系統(tǒng)總的能量密度為：

3.5.3 基于分子間作用勢(shì)的DBM

在對(duì)玻爾茲曼方程進(jìn)行介紹時(shí)提到，玻爾茲曼方程的碰撞項(xiàng)是基于理想氣體假設(shè)的，是忽略分子體積效應(yīng)的。在處理稠密氣體或液體時(shí)，這個(gè)假設(shè)并不合理。隨著分子數(shù)密度的增加，相對(duì)于平均分子間距，分子大小不再可以忽略。考慮分子的體積效應(yīng)，對(duì)玻爾茲曼碰撞算符進(jìn)行修正，可以得到恩斯柯格(Enskog)碰撞算符：

其中d0為 硬球分子直徑，χ為考慮分子體積效應(yīng)的碰撞概率修正，為分子中心相對(duì)位置單位矢量。對(duì)碰撞項(xiàng)在r處進(jìn)行泰勒展開(kāi)，并保持一階導(dǎo)數(shù)，可以得到：

其中χ、和f1均 為在位置r處的值。如果將(78)中后兩項(xiàng)中的f近似為局域平衡分布函數(shù)feq：

則恩斯柯格碰撞算符可寫(xiě)為：

使用BGK模型對(duì)上式右端第一項(xiàng)進(jìn)行簡(jiǎn)化，可得：

由式（82）則可以得到近似不可壓、等溫系統(tǒng)的簡(jiǎn)化恩斯柯格方程：

如果將外力項(xiàng)中的分布函數(shù)f用feq近似：

則可以將簡(jiǎn)化的恩斯柯格方程寫(xiě)為：

右端最后一項(xiàng)代表分子間斥力的作用，即為分子的體積效應(yīng)。

恩斯柯格方程引入了分子間斥力作用。如果在恩斯柯格方程中引入分子間的引力，就可以得到基于恩斯柯格方程的多相流模型。由平均場(chǎng)理論，分子間的相互吸引可用平均力勢(shì)[17]，

來(lái) 描 述，其 中r12=|r1-r2|是r1和r2兩 點(diǎn) 間 的 距 離，φattr(r12) 代 表 分 子 間 的 吸 引 勢(shì)。將 ρ(r2)在 點(diǎn)r1展 開(kāi)，忽略二階以上高階項(xiàng)， Φ(r1)可近似為：

將F表達(dá)式代入式（85），可得基于恩斯柯格方程的多相流模型：

其中，

F′表達(dá)式中第一項(xiàng)與狀態(tài)方程有關(guān)，第二項(xiàng)與表面張力有關(guān)。

由式（89）則可以得到相應(yīng)的DBM演化方程為：

3.6 離散玻爾茲曼建模：含化學(xué)反應(yīng)情形

在模擬含化學(xué)反應(yīng)的系統(tǒng)時(shí)，通過(guò)引入化學(xué)反應(yīng)項(xiàng)考慮化學(xué)反應(yīng)對(duì)系統(tǒng)的影響。以BGK模型為例，考慮化學(xué)反應(yīng)的玻爾茲曼-BGK方程具有如下形式：

其中：f為分布函數(shù)；feq為對(duì)應(yīng)的平衡態(tài)分布函數(shù)；τ為熱力學(xué)弛豫時(shí)間；v為 分子速度；C為化學(xué)反應(yīng)項(xiàng)，描述由于化學(xué)反應(yīng)而引起的分布函數(shù)f的變化率。

通常情況下，系統(tǒng)化學(xué)反應(yīng)的時(shí)間尺度tC是遠(yuǎn)大于分子熱力學(xué)弛豫的時(shí)間尺度 τ的。例如，常溫常壓條件下氫氣系統(tǒng)的熱力學(xué)弛豫時(shí)間為1 0-10s量級(jí)，而氫氣爆燃或爆轟的時(shí)間尺度為1 0-5s量級(jí)。因此，可以近似認(rèn)為，在化學(xué)反應(yīng)過(guò)程中，系統(tǒng)是始終處于熱力學(xué)平衡態(tài)的。這樣，可以得到：

其中，f為不考慮化學(xué)反應(yīng)時(shí)系統(tǒng)的平衡態(tài)分布函數(shù)，f*eq=f*eq(ρ,u,T*)為考慮化學(xué)反應(yīng)貢獻(xiàn)后系統(tǒng)的平衡態(tài)分布函數(shù)。

若化學(xué)反應(yīng)不可逆，則反應(yīng)進(jìn)程可由下面的反應(yīng)率方程來(lái)描述：

其中Q為單位質(zhì)量的反應(yīng)物發(fā)生反應(yīng)可以釋放的熱量。由式（94）～（96）可以看出，化學(xué)反應(yīng)項(xiàng)C的強(qiáng)弱不僅取決于反應(yīng)進(jìn)程，還受到反應(yīng)放熱Q的影響。即使化學(xué)反應(yīng)速率很快，當(dāng)Q很小時(shí)，化學(xué)反應(yīng)項(xiàng)C的貢獻(xiàn)也可能很小。

由于單弛豫時(shí)間模型是多弛豫時(shí)間模型的特例，所以模擬燃燒系統(tǒng)的DBM演化方程可統(tǒng)一由下式描述[27]：

其中，i=1,2,···,N為離散速度編號(hào)，N為離散速度的數(shù)目是fi和的動(dòng)理學(xué)矩；RN) 為碰撞參數(shù)，描述趨于平衡態(tài)值的快慢，其倒數(shù)給出相應(yīng)模式的弛豫時(shí)間；為保證離散玻爾茲曼模型能夠描述正確的流動(dòng)行為而做的必要修正。這是 因 為，盡 管 從 數(shù) 學(xué) 角 度 來(lái) 說(shuō)的 松 弛 因 子可以獨(dú)立調(diào)節(jié)，但從物理角度來(lái)說(shuō)，不同動(dòng)理學(xué)模式之間可能存在耦合，需要通過(guò)Chapman-Enskog多尺度分析或其他方法（例如實(shí)驗(yàn)標(biāo)定、經(jīng)過(guò)驗(yàn)證的唯象理論或模型等），來(lái)找回丟失的關(guān)聯(lián)[49]。修正項(xiàng)Ai必須是Kn的一階項(xiàng)。對(duì)于單弛豫時(shí)間模型 ，

對(duì)于多弛豫時(shí)間模型，Chapman-Enskog多尺度分析按如下方式展開(kāi)：

3.7 離散玻爾茲曼建模：多介質(zhì)情形

與宏觀二流體、多流體模型相對(duì)應(yīng)，DBM也有二流體、多流體模型。N流體模型使用N個(gè)分布函數(shù)描述系統(tǒng)的狀態(tài)，每個(gè)分布函數(shù)描述一種介質(zhì)。根據(jù)趨于平衡的次序，有單步（弛豫）碰撞模型和多步（弛豫）碰撞模型（一般是二步碰撞模型）；根據(jù)弛豫時(shí)間的個(gè)數(shù)，有單弛豫時(shí)間模型和多弛豫時(shí)間模型。

對(duì)于多介質(zhì)流體系統(tǒng)，引入速度空間和動(dòng)理學(xué)矩相空間的離散分布函數(shù)和。這里的下標(biāo)i(=1,2,···,m)對(duì)應(yīng)離散速度，α表示坐標(biāo)分量，如在三維空間直角坐標(biāo)系中代表xyz。上標(biāo) σ為流體系統(tǒng)中介質(zhì)的編號(hào)。分別是介質(zhì) σ的局域質(zhì)量密度、（摩爾或）粒子數(shù)密度、動(dòng)量和流速。

其中mσ是（摩爾或）粒子質(zhì)量?；旌衔锞钟虻目傎|(zhì)量密度 ρ、（摩爾或）粒子數(shù)密度n、 總動(dòng)量Jα和平均流速u(mài)α分別由下面的關(guān)系式得到：

介質(zhì) σ的局域能量和混合物總局域能量分別為：

比單介質(zhì)情形復(fù)雜的是，內(nèi)能（溫度）的定義依賴(lài)于作為參考的流速，而在多介質(zhì)情形，既有介質(zhì) σ的流速，又有混合物的（平均）流速u(mài)α。首先，借助混合物的平均流速u(mài)α定 義介質(zhì) σ的溫度和混合物（系統(tǒng)）溫度：

D為空間維數(shù)，Iσ是介質(zhì) σ的額外自由度數(shù)目。同時(shí)，定義介質(zhì) σ相對(duì)于自己流速的溫度：

至此，引入了三種溫度。溫度和壓強(qiáng)之間由狀態(tài)方程相聯(lián)系。有幾種溫度，就有幾種壓強(qiáng)。鑒于問(wèn)題的復(fù)雜性，在本節(jié)討論中暫且忽略分子間作用力，使用理想氣體狀態(tài)方程，給出一種建模思路。對(duì)于理想氣體情形，首先定義介質(zhì) σ基于混合物（平均）流速的壓強(qiáng)：

則混合物的壓強(qiáng)：

同時(shí)，定義介質(zhì) σ 基于自己流速的壓強(qiáng)：

可見(jiàn)，如果將混合物（平均）流速作為參考，則混合物的總壓強(qiáng)等于各介質(zhì)的分壓強(qiáng)之和。

平衡態(tài)分布函數(shù)依賴(lài)于粒子數(shù)密度、流速和溫度，溫度的定義依賴(lài)于流速，而我們有兩種流速—介質(zhì) σ 的流速和混合物的（平均）流速u(mài)α。所以，針對(duì)介質(zhì) σ ，可以引入三種平衡態(tài)分布函數(shù)：

與此對(duì)應(yīng)，動(dòng)理學(xué)矩空間的分布函數(shù)也有三種：

二步（弛豫）碰撞模型的思路是：介質(zhì)內(nèi)先平衡，介質(zhì)間再平衡。演化方程可寫(xiě)為：

3.8 示蹤粒子與DBM的耦合建模

在單流體模型中，只有一套流體力學(xué)量 (ρ、u、T)，描述的是系統(tǒng)局域的總密度、平均流速和平均溫度，無(wú)法識(shí)別混合過(guò)程中物質(zhì)粒子來(lái)源于哪一組分。受到顆粒示蹤實(shí)驗(yàn)的啟發(fā)，張戈等發(fā)展了示蹤粒子與DBM的耦合建模，實(shí)現(xiàn)了在單流體模型框架下，識(shí)別物質(zhì)粒子的來(lái)源[51]。更重要的是，示蹤粒子的引入，為流場(chǎng)動(dòng)理學(xué)研究提供了一個(gè)嶄新的視角。

在含示蹤粒子的系統(tǒng)中，我們可以使用斯托克斯數(shù)（Stokes number,St）來(lái)描述該粒子的動(dòng)力學(xué)狀態(tài)：

其中 τP是粒子的特征弛豫時(shí)間，u0是當(dāng)?shù)氐牧鲃?dòng)速度，l0是特征長(zhǎng)度（通常選取顆粒的直徑）。當(dāng)St＞1時(shí)，粒子將由于慣性脫離當(dāng)?shù)亓鲃?dòng)而運(yùn)動(dòng)，特別是流速變化劇烈的情況下；當(dāng)St?1 時(shí)，粒子緊緊貼著流線運(yùn)動(dòng)。通過(guò)調(diào)整St到足夠小的數(shù)量級(jí)，則該粒子能夠作為流場(chǎng)的示蹤粒子使用。假設(shè)粒子的弛豫時(shí)間τP接近0，則其慣性完全可以忽略，其運(yùn)動(dòng)速度瞬間達(dá)到當(dāng)?shù)亓魉?，因而完全跟隨流體運(yùn)動(dòng)。

我們往往需要引入一批示蹤粒子。沿著每一個(gè)示蹤粒子的軌跡，都可以給出拉格朗日視角的基于示蹤粒子的描述。對(duì)于第k個(gè)粒子來(lái)說(shuō)，其運(yùn)動(dòng)方程如下：

其中，rk是 第k個(gè)示蹤粒子的空間位置，uP為示蹤粒子的運(yùn)動(dòng)速度。

示蹤粒子往往需要盡可能的小。假設(shè)其體積與質(zhì)量極小，在流場(chǎng)中用一個(gè)點(diǎn)來(lái)表示，其與流體之間的動(dòng)量交換在瞬間完成，那么示蹤粒子的速度就直接由局域的流動(dòng)速度決定：

其中a為示蹤粒子在流場(chǎng)中所處位置的微元，D表示對(duì)示蹤粒子速度具有影響的流場(chǎng)積分區(qū)域，?為Dirac函數(shù)，在模擬中通常使用其離散形式 ψ代替。第k個(gè)示蹤粒子的速度將根據(jù)它所處的位置以及當(dāng)?shù)氐牧鲃?dòng)速度決定，例如經(jīng)過(guò)時(shí)間間隔 Δt，示蹤粒子速度將從變化為為了更新點(diǎn) 顆 粒 的 位 置，可使用四階龍格-庫(kù)塔格式求解離散顆粒的運(yùn)動(dòng)方程。

因?yàn)槭聚櫫Ｗ釉谶\(yùn)動(dòng)過(guò)程中，很難剛好落在網(wǎng)格點(diǎn)上，所以，其速度需要根據(jù)它附近的流體網(wǎng)格點(diǎn)的速度確定。具體而言，就是將附近的網(wǎng)格點(diǎn)上的速度根據(jù)網(wǎng)格點(diǎn)位置與示蹤粒子位置的距離遠(yuǎn)近而作加權(quán)平均，在數(shù)學(xué)上通過(guò)使用離散的Dirac函數(shù)（ ψ）來(lái)實(shí)現(xiàn)。在二維情況下，

ψ 函數(shù)可以被分解為兩個(gè)部分：

其中，i,j為網(wǎng)格點(diǎn)編號(hào)。張戈等在文獻(xiàn)[51]中所應(yīng)用的權(quán)重函數(shù) φ為：

據(jù)此，示蹤粒子從流場(chǎng)中獲得了自己的運(yùn)動(dòng)速度。

3.9 離散玻爾茲曼建模：等離子體情形

等離子體是由大量處于非束縛態(tài)的帶電粒子組成的具有宏觀時(shí)間和空間尺度的體系。在等離子體中，時(shí)空尺度以及密度、溫度等宏觀狀態(tài)量可以跨越10個(gè)數(shù)量級(jí)及以上的范圍，這使得：（1）相關(guān)的特征時(shí)間和空間尺度極其豐富，例如研究中常用到的空間尺度有德拜半徑、拉莫爾半徑等，常用到的時(shí)間尺度有朗繆爾振蕩周期等；（2）研究中使用的物理模型和方法從宏觀到微觀種類(lèi)繁多，如磁流體力學(xué)、弗拉索夫動(dòng)理學(xué)方程、粒子模擬方法（如質(zhì)點(diǎn)網(wǎng)格法，Particle in Cell，PIC）等。

受控?zé)岷司圩兪墙鉀Q人類(lèi)能源問(wèn)題的關(guān)鍵，其中等離子體的運(yùn)動(dòng)處于高溫高密高壓的極端環(huán)境，同時(shí)涉及到多時(shí)空尺度的多場(chǎng)耦合及帶電粒子間的碰撞輸運(yùn)過(guò)程，是人們關(guān)注的重點(diǎn)。

在等離子體介觀動(dòng)理學(xué)模擬中，一般采用德拜半徑作為特征空間尺度將等離子體間的相互作用劃分為自洽電磁場(chǎng)（長(zhǎng)程外力項(xiàng)部分）和碰撞（短程輸運(yùn)部分），描述方程如式(127)所示：

其中α、β代表等離子體中不同種類(lèi)的帶電粒子，左側(cè)第三項(xiàng)為帶電粒子所受的除洛倫茲力外的合力，左側(cè)第四項(xiàng)為自洽電磁場(chǎng)部分（需要耦合求解Maxwell方程組），右側(cè)為帶電粒子間的碰撞效應(yīng)，其中i、e分別表示離子和電子。

由于作用的時(shí)空尺度不同，人們往往在長(zhǎng)程作用及碰撞輸運(yùn)作用中選其一進(jìn)行研究，這種取舍一方面是基于部分物理現(xiàn)實(shí)問(wèn)題的尺度分離（具有合理性），另一方面也是由于不同尺度耦合問(wèn)題的復(fù)雜性和處理方法的不足甚至是缺乏（具有無(wú)奈性）。但是，對(duì)于有些問(wèn)題如靜電激波陣面結(jié)構(gòu)、慣性約束聚變中的流體不穩(wěn)定性問(wèn)題等，兩種作用的時(shí)空尺度接近，無(wú)法將其中之一作為微擾，因而兩種作用需要同時(shí)加以考慮。

等離子體中自洽電磁場(chǎng)同等離子體運(yùn)動(dòng)相耦合，采用有限差分的Maxwell方程組解法主要有時(shí)域有限差分（Finite difference time domain, FDTD或Yee格式）、交替方向隱式迭代時(shí)域有限差分（Alternating direction implicit-FDTD, ADI-FDTD）及分裂算子時(shí)域有限差分（Splitting operator-FDTD, S-FDTD）等。

等離子體中粒子間多體碰撞可以看作發(fā)生在德拜半徑以?xún)?nèi)，根據(jù)適用的碰撞算子的復(fù)雜性可以分為Boltzmann算子、Fokker-Planck算子、Landau算子以及BGK算子，其中BGK算子由于較為簡(jiǎn)單適用性最為廣泛。和多介質(zhì)DBM類(lèi)似，等離子體DBM也可以分為單步（弛豫）建模及多步（弛豫）建模。

忽略中間動(dòng)理學(xué)過(guò)程，采用由Andries[52]等發(fā)展的AAP單步弛豫模型的等離子體動(dòng)理學(xué)模型為：

若考慮不同種粒子間的碰撞弛豫（中間動(dòng)理學(xué)過(guò)程），采用由Haack[53]等發(fā)展的等離子體動(dòng)理學(xué)模型為：

其 中 vαβ(c) 為 不 同 類(lèi) 型 粒 子 間 碰 撞 的 頻 率，Tαβ和uαβ分別為不同類(lèi)型粒子間碰撞的中間弛豫溫度和速度。對(duì)于同種粒子間的碰撞，Tαβ和uαβ直接取該種粒子的宏觀溫度和速度。對(duì)于不同種粒子間的碰撞，需要耦合相應(yīng)的溫度及速度模型，如：

式（131）和式（132）分別為滿足動(dòng)量弛豫約束的BGKEM模型和滿足能量弛豫約束的BGK-ET模型。目前，基于式（130）及有限差分的Maxwell方程組解法正用于等離子體靜電激波及相關(guān)問(wèn)題的研究。

如前面所述，將離子和電子分別作為不同的流體介質(zhì)，可以使用雙流體玻爾茲曼模型。當(dāng)然，還存在其他不同粗?；潭鹊奈锢斫Ｐ问剑绮捎梅植己瘮?shù)描述離子的行為特征，而假設(shè)電子始終處于局域熱力學(xué)平衡態(tài)或?qū)﹄娮硬捎昧黧w描述方法。

在等離子體系統(tǒng)的動(dòng)理學(xué)描述中，系統(tǒng)行為由相應(yīng)的分布函數(shù)fα來(lái)描述。而要確切地掌握分布函數(shù)fα，等價(jià)于確切地掌握分布函數(shù)fα的所有可能的動(dòng)理學(xué)矩。這對(duì)于很多實(shí)際情形，是既不可能，又不必要的。對(duì)于這些情形，只需要根據(jù)研究需求，抓主要矛盾，確切掌握分布函數(shù)fα的部分動(dòng)理學(xué)性質(zhì)。因此需要進(jìn)一步對(duì)模型進(jìn)行簡(jiǎn)化。簡(jiǎn)化過(guò)程中要遵循的原則是—描述這部分動(dòng)理學(xué)性質(zhì)的動(dòng)理學(xué)矩其結(jié)果不因模型的簡(jiǎn)化而改變。這正是等離子體系統(tǒng)DBM建模的初衷和任務(wù)。

DBM主要針對(duì)的是非平衡效應(yīng)較強(qiáng)，以至于傳統(tǒng)流體建模失效，同時(shí)粒子之間碰撞效應(yīng)又不能忽略的情形。目前，等離子體系統(tǒng)的DBM建模與模擬研究尚處于起步階段，后期我們將開(kāi)展進(jìn)一步研究。

4 多相流機(jī)理研究：基于DBM

作為系統(tǒng)行為粗?；枋龅囊环N物理模型構(gòu)建方法，DBM本身并不包含具體的數(shù)值離散格式。獲得了DBM模型，跟獲得了NS模型一樣，在數(shù)值模擬中還需根據(jù)問(wèn)題特點(diǎn)選取合適的數(shù)值方法。從物理問(wèn)題研究的角度，可以選用滿足當(dāng)次物理問(wèn)題研究需求的任何一種離散格式。在下面物理結(jié)果的原始文獻(xiàn)中使用的離散格式只是滿足當(dāng)次研究需求的多種離散格式之一。

4.1 流體不穩(wěn)定性方面

流體界面不穩(wěn)定性（經(jīng)常簡(jiǎn)稱(chēng)流體不穩(wěn)定性）與物質(zhì)混合現(xiàn)象廣泛存在于自然界、慣性約束核聚變和武器物理等領(lǐng)域。慣性約束核聚變和武器物理等領(lǐng)域重點(diǎn)關(guān)注的流體不穩(wěn)定性主要有三種：RM不穩(wěn)定性（Richtmyer-Meshkov instability, RMI）、RT不穩(wěn)定性（Rayleigh-Taylor instability, RTI）和KH不穩(wěn)定性（Kelvin-Helmholtz instability, KHI）。

流體不穩(wěn)定性系統(tǒng)通常具有以下特點(diǎn)：系統(tǒng)的本身可能是宏觀尺度的，但其內(nèi)部存在大量的中間尺度的空間結(jié)構(gòu)和動(dòng)理學(xué)模式；這些結(jié)構(gòu)和模式的存在與演化極大地影響著系統(tǒng)的物理性能和功能。系統(tǒng)內(nèi)部往往具有大量的界面，包括物質(zhì)界面和力學(xué)界面，系統(tǒng)內(nèi)部的受力和響應(yīng)過(guò)程非常復(fù)雜。對(duì)于這些系統(tǒng)的研究，NS方程只能描述大尺度緩變行為，而對(duì)于一些銳利的邊界，流體的平均分子間距相對(duì)于界面尺度已經(jīng)不再是可以忽略的小量，基于連續(xù)介質(zhì)假設(shè)的NS方程受到挑戰(zhàn)。在一些快變流動(dòng)模式描述方面，流體系統(tǒng)趨于熱力學(xué)平衡的弛豫時(shí)間相對(duì)于該流動(dòng)的特征時(shí)間來(lái)講，不再是可以忽略的小量，熱力學(xué)非平衡效應(yīng)較強(qiáng)，只考慮一階熱力學(xué)非平衡效應(yīng)的NS描述不再能滿足需求。

對(duì)于流體不穩(wěn)定性問(wèn)題，DBM可以根據(jù)需要研究的非平衡程度和選定研究的具體非平衡行為特征構(gòu)建滿足需求的物理模型。在只考慮一階熱力學(xué)非平衡效應(yīng)的情況，除了NS可以描述的流體行為，借助所定義的各種非平衡特征量，DBM可以給出NS所遺漏的一系列非平衡行為的描述。

2016年，利用單弛豫時(shí)間DBM，賴(lài)惠林等重點(diǎn)研究了流體界面不穩(wěn)定性演化過(guò)程中的熱力學(xué)非平衡效應(yīng)[54-55]。研究發(fā)現(xiàn)，可壓性對(duì)RT界面演化呈現(xiàn)出“先抑制、減速，后促進(jìn)、加速”的階段性，這一階段性可從能量轉(zhuǎn)換角度獲得很好的解釋。在每個(gè)時(shí)刻，所有非平衡動(dòng)理學(xué)模式均隨可壓性增強(qiáng)而增強(qiáng)；隨可壓性增強(qiáng)，可觀測(cè)的非平衡效應(yīng)越來(lái)越豐富，后期高階非平衡動(dòng)理學(xué)模式慢慢凸顯出來(lái)；在不同階段，非平衡模式之間相對(duì)強(qiáng)弱會(huì)發(fā)生改變；某些非平衡動(dòng)理學(xué)模式的強(qiáng)度始終較小。

陳鋒等使用多弛豫時(shí)間（Multiple Relaxation Time,MRT）DBM從宏觀行為和非平衡特征兩個(gè)角度研究Rayleigh-Taylor不穩(wěn)定性問(wèn)題，尤其探討了以下兩方面問(wèn)題：（1）系統(tǒng)內(nèi)密度、流速、溫度、壓強(qiáng)等宏觀量的不均勻度與各種不同形式的非平衡行為之間的關(guān)聯(lián)度；（2）黏性、熱傳導(dǎo)、Prandtl數(shù)對(duì)界面擾動(dòng)增長(zhǎng)過(guò)程、對(duì)上述各類(lèi)關(guān)聯(lián)的影響[56]。

2019年，借助DBM，甘延標(biāo)等研究了KH不穩(wěn)定性系統(tǒng)的流變和形態(tài)特征[57]。重點(diǎn)研究了傳統(tǒng)流體建模所忽略、而MD模擬因適用時(shí)空尺度受限而無(wú)法直接研究的熱力學(xué)非平衡行為和效應(yīng)。同時(shí)，為解決KH不穩(wěn)定性演化過(guò)程中各類(lèi)復(fù)雜物理場(chǎng)的分析問(wèn)題，他們提出了通過(guò)追蹤非平衡行為特征和形態(tài)分析技術(shù)相結(jié)合[58-60]，進(jìn)行特征結(jié)構(gòu)或模式的物理甄別與追蹤技術(shù)設(shè)計(jì)，定量表征KH混合層寬度和發(fā)展速率的新途徑。借助雙介質(zhì)DBM，林傳棟等人更加細(xì)致地研究了RT不穩(wěn)定性系統(tǒng)的非平衡行為特性[61]。

陳鋒等人的研究[62]表明，在RM不穩(wěn)定性或RT不穩(wěn)定性演化過(guò)程中，系統(tǒng)內(nèi)的溫度不均勻度和無(wú)組織能量流（Non-Organized Energy Flux, NOEF）之間、密度不均勻度和TNE強(qiáng)度之間的相關(guān)度較高，接近1；速度不均勻度和無(wú)組織動(dòng)量流（Non-Organized Momentum Flux，NOMF）強(qiáng)度之間相關(guān)度也相對(duì)較高（見(jiàn)圖4）；熱傳導(dǎo)是影響相關(guān)度的重要因素（見(jiàn)圖5）。

圖4 RT與RM不穩(wěn)定性演化過(guò)程中系統(tǒng)內(nèi)不同類(lèi)型的非平衡強(qiáng)度與密度、溫度、流速不均勻度之間關(guān)聯(lián)程度的對(duì)比[56, 62]（更多細(xì)節(jié)參見(jiàn)文獻(xiàn)[62]）Fig. 4 Comparison of the correlation degree between different types of non-equilibrium strength and density, temperature and flow rate unevenness in the evolution of RT and RM instability[56, 62] (see Ref. [62] for more details)

圖5 RM不穩(wěn)定性演化過(guò)程中熱傳導(dǎo)對(duì)相關(guān)度的影響[62]（更多細(xì)節(jié)參見(jiàn)文獻(xiàn)[62]）Fig. 5 Effect of heat conduction on the degree of correlation in the evolution of RM instability[62] (see Ref. [62] for more details)

相比于RT不穩(wěn)定性系統(tǒng)，在RM不穩(wěn)定性系統(tǒng)中，沖擊波的透射和反射使得相關(guān)度演化曲線出現(xiàn)“跳變”（見(jiàn)圖4）。在RT不穩(wěn)定性與RM不穩(wěn)定性共存的系統(tǒng)中，重力加速度和沖擊波強(qiáng)度之間合作、競(jìng)爭(zhēng)，共同主導(dǎo)著界面演化與物質(zhì)和能量混合過(guò)程。在RM不穩(wěn)定性主導(dǎo)的情形初始擾動(dòng)界面會(huì)出現(xiàn)反轉(zhuǎn)（見(jiàn)圖6-圖7，更多細(xì)節(jié)參見(jiàn)文獻(xiàn)[62]）。

RT不穩(wěn)定性和RM不穩(wěn)定性分別主導(dǎo)的參數(shù)區(qū)間如圖8所示。圖9展示了全局平均TNE強(qiáng)度、全局平均密度梯度、全局平均NOMF強(qiáng)度以及全局平均NOEF強(qiáng)度隨時(shí)間的演化；重力加速度對(duì)非平衡行為特征的影響表現(xiàn)出階段性（見(jiàn)圖9）；隨著馬赫數(shù)的增加，系統(tǒng)的整體非平衡強(qiáng)度指數(shù)增加，而溫度不均勻度與NOEF的相關(guān)度指數(shù)衰減（見(jiàn)圖10）。圖8～圖10為數(shù)值模擬結(jié)果，更多細(xì)節(jié)參見(jiàn)文獻(xiàn)[62]。

圖6 RTI與RMI共存系統(tǒng)界面反轉(zhuǎn)現(xiàn)象出現(xiàn)機(jī)制的示意圖[62]Fig. 6 The schematic diagram of the collaboration and competition relations between RM and RT instability[62]

圖7 RTI與RMI共存系統(tǒng)界面反轉(zhuǎn)現(xiàn)象出現(xiàn)與否的數(shù)值模擬密度云圖[62]Fig. 7 Numerical simulation density nephogram of interface inversion in RTI and RMI coexistence system[62]

圖8 RTI與RMI共存系統(tǒng)宏觀特征[62]Fig. 8 Macrocharacteristics of the RTI and RMI coexistence system[62]

圖9 不同重力場(chǎng)下RTI與RMI共存系統(tǒng)的非平衡行為特征[62]Fig. 9 Non-equilibrium characteristics of RTI and RMI coexistence system under different gravity fields[62]

2020年，基于多弛豫時(shí)間DBM，結(jié)合形態(tài)學(xué)分析、時(shí)空關(guān)聯(lián)等方法，陳鋒等進(jìn)一步對(duì)耦合瑞利-泰勒-開(kāi)爾文-亥姆霍茲不穩(wěn)定性（RTKHI）系統(tǒng)進(jìn)行了研究[63]。研究發(fā)現(xiàn)：溫度場(chǎng)圖靈斑圖的總邊界長(zhǎng)度L和平均熱通量強(qiáng)度D3,1均可用于測(cè)量浮力與剪切強(qiáng)度之比，并定量評(píng)估RTKHI系統(tǒng)早期階段的主要機(jī)理。針對(duì)早期、KHI主導(dǎo)，后期RTI主導(dǎo)的耦合RTKHI系統(tǒng)，形態(tài)學(xué)邊界長(zhǎng)度L可以很好地捕獲從KHI主導(dǎo)到RTI主導(dǎo)的過(guò)渡點(diǎn)。L線性增加的終點(diǎn)可以作為區(qū)分這兩個(gè)階段的幾何準(zhǔn)則。TNE量之一，熱通量強(qiáng)度D3,1與邊界長(zhǎng)度L表現(xiàn)出相似的行為，并且在早期呈現(xiàn)很強(qiáng)的正相關(guān)。因此，D3,1線性增加的終點(diǎn)可以作為區(qū)分這兩個(gè)階段的物理準(zhǔn)則。形態(tài)邊界長(zhǎng)度L是高溫和低溫（輕質(zhì)和重質(zhì)）流體之間的界面長(zhǎng)度。它反映了不穩(wěn)定發(fā)展和材料混合的程度。TNE量D3,1反映了系統(tǒng)不同區(qū)域偏離平衡的程度，并且可以更清楚準(zhǔn)確地定位界面的位置。L和D3,1這兩個(gè)標(biāo)準(zhǔn)從不同角度出發(fā)，但是是一致的，具有各自的優(yōu)勢(shì)，可以互相補(bǔ)充。采用這兩個(gè)標(biāo)準(zhǔn)有助于識(shí)別耦合RTKHI系統(tǒng)的主要機(jī)制和關(guān)鍵時(shí)間點(diǎn)。

圖11為重力加速度g=0.005、 切向速度差u0=0.05的耦合RTKHI系統(tǒng)的溫度圖像、 圖靈斑圖以及u0=0.05的純KHI的圖靈斑圖?？梢钥闯?，和單純的KHI系統(tǒng)相比，耦合RTKHI系統(tǒng)具有傾斜的、非對(duì)稱(chēng)的氣泡、尖釘以及蘑菇狀結(jié)構(gòu)，單純的KHI系統(tǒng)的演化大大落后于耦合RTKHI系統(tǒng)的演化。在此情況下，RTI扮演了一個(gè)主要角色，切向速度的存在主要破壞了RTI結(jié)構(gòu)的對(duì)稱(chēng)性。更多情況的算例及物理分析可參見(jiàn)文獻(xiàn)[63]。

圖10 馬赫數(shù)對(duì)RTI與RMI共存系統(tǒng)TNE強(qiáng)度、宏觀量梯度與非平衡特征相關(guān)度的影響[62]Fig. 10 Influence of Mach number on thermodynamic non-equilibrium strength, macroscopic gradient and non-equilibrium characteristic correlation of RTI and RMI coexistence system[62]

圖11 耦合RTKHI系統(tǒng)（g =0.005, u 0=0.05）的溫度圖像、圖靈斑圖，以及u 0=0.05 的純KHI的圖靈斑圖 [63]Fig. 11 A coupled RTKHI with g =0.005 and u 0=0.05. (a) and(b) are the temperature and the corresponding Turing pattern(T th=1.0) of the RTKHI, respectively. (c) is the temperature Turing pattern of pure KHI with u 0=0.05[63]

圖12展示了擾動(dòng)幅值A(chǔ)、擾動(dòng)幅值的增長(zhǎng)率dA/dt、形態(tài)學(xué)邊界長(zhǎng)度L隨時(shí)間的演化。對(duì)于KHI在早期階段占主導(dǎo)、RTI在后期階段占主導(dǎo)的情形，演化過(guò)程可以粗略地劃分為兩個(gè)階段，分別為剪切主導(dǎo)階段和浮力主導(dǎo)階段。由剪切主導(dǎo)階段到浮力主導(dǎo)階段的過(guò)渡狀態(tài)稱(chēng)為過(guò)渡點(diǎn)。從圖12中可以看到，在開(kāi)始階段，LRTKHI首先呈指數(shù)增長(zhǎng)，然后呈線性增長(zhǎng)（虛線所示）。LRTKHI線性增長(zhǎng)結(jié)束的點(diǎn)近似等于RTKHI系統(tǒng)幅值增長(zhǎng)率最小值點(diǎn)以及相關(guān)KHI系統(tǒng)的幅值最大值點(diǎn)。從這一時(shí)刻開(kāi)始，RTI開(kāi)始扮演一個(gè)主要角色。這一時(shí)刻稱(chēng)為過(guò)渡點(diǎn)，在圖中使用豎直虛線和圓標(biāo)記。因此，LRTKHI線性增長(zhǎng)結(jié)束的點(diǎn)可以作為區(qū)別兩個(gè)階段的幾何判據(jù)。剪切率越大，過(guò)渡點(diǎn)越早出現(xiàn)。

圖13為耦合RTKHI系統(tǒng)早期主要機(jī)制的形態(tài)分析曲線。對(duì)不穩(wěn)定性發(fā)展程度和材料混合程度，形態(tài)學(xué)總邊界長(zhǎng)度的值是一個(gè)很有用和有效的指標(biāo)。更進(jìn)一步，溫度場(chǎng)圖靈斑圖的邊界長(zhǎng)度L可以用來(lái)測(cè)量浮力和剪切強(qiáng)度的比值，因此，它可以用來(lái)定量地判斷RTKHI系統(tǒng)早期階段的主要作用機(jī)理。特別的，當(dāng)KHI（RTI）主導(dǎo)時(shí)，LKHI＞LRTI（LKHI＜LRTI）；當(dāng)KHI和RTI相互平衡時(shí)，LRTI=LKHI。

圖14展示了t=150 時(shí)刻和t=250時(shí)刻的非平衡量以及相關(guān)的NOEF強(qiáng)度d3,1。在圖14中，由d3,1可以看到一個(gè)非常清楚的雙渦結(jié)構(gòu)。這些信息和結(jié)構(gòu)不能從（或者不容易從）溫度云圖（圖11（a））看出來(lái)。與個(gè)別分量相比，NOEF強(qiáng)度d3,1提供了一個(gè)高分辨率的交界面，可以用來(lái)描述RTKHI模擬中的完整界面。

圖12 振幅 A、振幅增長(zhǎng)率d A/dt 和形態(tài)邊界長(zhǎng)度 L 與時(shí)間 t 的關(guān)系[63]Fig. 12 Perturbation amplitude A, growth rate d A/dt, and morphological boundary length L vs time t [63]

圖13 根據(jù)溫度場(chǎng)的形態(tài)分析早期主要機(jī)制[63]Fig. 13 Morphological analysis of the main mechanism in the early stage[63]

圖15（a）～15（d）展示了非平衡特征在早期主要機(jī)理判斷中的應(yīng)用。和全局平均TNE強(qiáng)度DTNE相比，全局平均NOEF強(qiáng)度D3,1可以更精確地判斷早期階段的主要機(jī)理，并且計(jì)算結(jié)果和形態(tài)學(xué)邊界長(zhǎng)度L是一致的。由圖15（e）～15（f）可以看出，D3,1展示了和邊界長(zhǎng)度L相似的行為。交界面長(zhǎng)度L越大，D3,1越大。D3,1線性增長(zhǎng)的結(jié)束點(diǎn)，可以用來(lái)作為區(qū)分從KHI主導(dǎo)到RTI主導(dǎo)兩個(gè)階段的物理判據(jù)。具體的物理分析可參見(jiàn)文獻(xiàn)[63]。

圖14 RTKHI系統(tǒng)(g =0.005, u 0=0.1)，不同視角的非平衡特征[63]Fig. 14 Non-equilibrium characteristics of the RTKHI system(g =0.005, u 0=0.1) [63]

圖15 耦合RTKHI系統(tǒng)的非平衡特征在早期主要機(jī)理判斷和過(guò)渡點(diǎn)捕獲中的應(yīng)用[63]Fig. 15 The applications of non-equilibrium characteristics in the early main mechanism judgment and transition point capture[63]

由于流體界面不穩(wěn)定性系統(tǒng)的DBM建模與模擬尚處于起步階段，所以以上研究均是基于理想氣體模型的。表面張力、材料強(qiáng)度等對(duì)流體不穩(wěn)定性演化過(guò)程的影響，以及更加實(shí)際系統(tǒng)中流體不穩(wěn)定性的DBM建模與模擬有待進(jìn)一步研究。

4.2 相分離方面

當(dāng)介質(zhì)由高溫均勻態(tài)突然降溫到兩相共存區(qū)域時(shí)，原均勻態(tài)會(huì)因自由能太高而失穩(wěn)，從而發(fā)生相變和兩相的分離，這個(gè)過(guò)程通常稱(chēng)為淬火。在這個(gè)過(guò)程中，流體系統(tǒng)經(jīng)歷兩個(gè)動(dòng)力學(xué)階段—前期的亞穩(wěn)相分解或旋節(jié)線分解（spinodal decomposition, SD）階段和后期的相疇增長(zhǎng)（domain growth, DG）階段。兩個(gè)階段分別對(duì)應(yīng)于單相區(qū)域的形成和相區(qū)的融合增長(zhǎng)。

在早期相分離研究中，怎樣準(zhǔn)確地劃分亞穩(wěn)相分解和相疇增長(zhǎng)兩個(gè)階段一直是個(gè)開(kāi)放的問(wèn)題。一種方法是，兩個(gè)階段的臨界時(shí)間可以通過(guò)特征尺度的增長(zhǎng)特征來(lái)大致給出，即將標(biāo)度率開(kāi)始的時(shí)刻作為兩個(gè)階段的臨界點(diǎn)，但這樣區(qū)分兩個(gè)階段是非常粗略的。

2012年，甘延標(biāo)等發(fā)現(xiàn)：密度場(chǎng)圖靈斑圖的邊界總長(zhǎng)度L是描述相分離進(jìn)程的有效特征量—邊界總長(zhǎng)度L首先隨著時(shí)間逐漸增長(zhǎng)，后期則隨著時(shí)間逐漸減??；結(jié)合相分離兩個(gè)階段的其他物理特征，提出使用邊界總長(zhǎng)度L極大值點(diǎn)作為劃分亞穩(wěn)相分解和相疇增長(zhǎng)兩個(gè)階段的幾何判據(jù)；基于該判據(jù)，對(duì)亞穩(wěn)相分解的統(tǒng)計(jì)特征給出一些新的認(rèn)識(shí)[64]。但總體來(lái)說(shuō)，亞穩(wěn)相分解階段的熱力學(xué)非平衡行為特征還遠(yuǎn)遠(yuǎn)沒(méi)有獲得充分的理解。近期的研究[44,65]表明，熱力學(xué)非平衡強(qiáng)度和熵增率均在亞穩(wěn)相分解階段隨時(shí)間逐漸增大，在相疇增長(zhǎng)階段隨時(shí)間減小，所以熱力學(xué)非平衡強(qiáng)度和熵增率極大值點(diǎn)均可作為劃分這兩個(gè)階段的物理判據(jù)。根據(jù)時(shí)間先后，分別稱(chēng)之為第一和第二物理判據(jù)。2016年之前的相關(guān)研究綜述可參見(jiàn)文獻(xiàn)[66]。

受篇幅限制，這里只簡(jiǎn)單介紹一下2019年張玉東等發(fā)表在Soft Matter上的一項(xiàng)工作[44]。研究表明，部分非平衡特征和熵產(chǎn)生速率之間存在較為緊密的關(guān)系。相分離過(guò)程中的熵產(chǎn)生主要有兩種機(jī)制，分別來(lái)源于NOMF和NOEF。圖16所示為等溫與非等溫相分離過(guò)程熵產(chǎn)生速率演化特征的對(duì)比。由圖16可以看出，熵產(chǎn)生速率在相分離的第一階段（亞穩(wěn)相分解階段）隨時(shí)間增加，而在第二階段（相疇增長(zhǎng)階段）隨時(shí)間降低，因此熵產(chǎn)生速率的極大值點(diǎn)可以作為劃分相分離兩個(gè)階段的一個(gè)物理判據(jù)。

圖17至圖19所示為相分離過(guò)程中熱流、黏性和表面張力對(duì)相分離SD階段的持續(xù)時(shí)間和熵產(chǎn)生特性的影響。如圖17至圖19所示，熱流的作用是加快相分離的SD階段，黏性和表面張力的作用是延緩SD階段，并且熱流和黏性對(duì)SD持續(xù)時(shí)間的影響是（近似）指數(shù)形式的，而表面張力的影響是線性的。熱流對(duì)NOEF部分的熵產(chǎn)生速率起抑制作用，對(duì)NOMF部分的熵產(chǎn)生速率起促進(jìn)作用；黏性對(duì)NOEF部分的熵產(chǎn)生速率起促進(jìn)作用，對(duì)NOMF部分的熵產(chǎn)生速率起抑制作用；表面張力對(duì)NOEF和NOMF的熵產(chǎn)生速率都起抑制作用。其物理原因可以從流場(chǎng)中溫度梯度和速度梯度的變化來(lái)獲得解釋?zhuān)嗉?xì)節(jié)參見(jiàn)文獻(xiàn)[44]。

圖16 等溫與非等溫相分離過(guò)程熵產(chǎn)生速率演化特征的對(duì)比[44]Fig. 16 Comparison of entropy generation rate evolution characteristics between isothermal and non-isothermal phase separation processes[44]

圖17 熱流對(duì)熱相分離過(guò)程中SD階段持續(xù)時(shí)間和熵產(chǎn)生速率的影響[44]Fig. 17 Effect of heat flow on duration of SD stage and entropy generation rate in thermal phase separation[44]

圖18 黏性對(duì)熱相分離過(guò)程中SD階段持續(xù)時(shí)間和熵產(chǎn)生速率的影響[44]Fig. 18 Effect of viscosity on duration of SD stage and entropy generation rate in thermal phase separation[44]

圖20給出了不同熱流、黏性和表面張力情形下，兩種熵產(chǎn)生機(jī)制之間的競(jìng)爭(zhēng)與協(xié)作關(guān)系。發(fā)現(xiàn)，在熱流和黏性的作用下，兩種熵產(chǎn)生機(jī)制之間存在競(jìng)爭(zhēng)關(guān)系；在表面張力作用下，兩種熵產(chǎn)生機(jī)制之間存在協(xié)作關(guān)系。

圖19 表面張力對(duì)熱相分離過(guò)程中SD階段持續(xù)時(shí)間和熵產(chǎn)生速率的影響[44]Fig. 19 Effect of surface tension on duration of SD stage and entropy generation rate in thermal phase separation process[44]

圖20 兩種熵產(chǎn)生速率幅值之間的關(guān)系曲線，圖中箭頭指向熱傳導(dǎo)系數(shù)（1/Pr）、黏性系數(shù)(τ)和表面張力系數(shù)（K）增大的方向[44]Fig. 20 The relationship curve between the amplitude of two kinds of entropy production rate, the arrow in the figure points to the increasing direction of heat conduction coefficient (1/Pr), viscosity coefficient (τ) and surface tension coefficient (K) [44]

4.3 化學(xué)反應(yīng)流方面

近年來(lái)，DBM在燃燒領(lǐng)域的應(yīng)用取得了一系列的進(jìn)展。

2013年，閆鉑等提出了一個(gè)模擬燃燒和爆轟現(xiàn)象的二維DBM模型，研究了爆轟波的精細(xì)物理結(jié)構(gòu)和附近的非平衡效應(yīng)[67]。這是2012年許愛(ài)國(guó)等提出“借助 (f-feq)的非守恒矩來(lái)描述非平衡狀態(tài)、提取非平衡信息”這一思路[29]以來(lái)的第一個(gè)具體應(yīng)用。該模型中的化學(xué)反應(yīng)使用Lee-Tarver模型描述,反應(yīng)熱和流體行為自然耦合。使用算子分裂的方法對(duì)爆轟系統(tǒng)進(jìn)行了成功的模擬。根據(jù)所定義的非平衡特征量，對(duì)系統(tǒng)偏離熱動(dòng)平衡態(tài)所引發(fā)的效應(yīng)獲得一些新的理解。2014年，林傳棟等提出了一個(gè)用于爆轟現(xiàn)象的極坐標(biāo)系下的DBM[68]，研究了一些典型的內(nèi)爆和外爆過(guò)程。為了探索燃燒過(guò)程中的流體動(dòng)力學(xué)非平衡和熱力學(xué)非平衡，許愛(ài)國(guó)等于2015年提出了一個(gè)二維多弛豫時(shí)間DBM；受統(tǒng)計(jì)物理學(xué)“使用粒子的位置x和速度v分量張開(kāi)相空間”描述方法的啟發(fā)，提出使用 (f-feq)非守恒矩的獨(dú)立分量張開(kāi)相空間，使用該相空間來(lái)描述系統(tǒng)偏離熱力學(xué)平衡引發(fā)的各種行為特征，進(jìn)一步在該相空間及其子空間中借助到原點(diǎn)的距離提出相應(yīng)非平衡強(qiáng)度的概念[30]。模型具有可調(diào)的比熱比和普朗特?cái)?shù)，化學(xué)反應(yīng)釋放的熱量可以自動(dòng)耦合到流體系統(tǒng)中。模型適用于亞聲速流動(dòng)燃燒和超聲速流動(dòng)爆轟模擬。使用提出的模型，初步對(duì)穩(wěn)態(tài)和非穩(wěn)態(tài)一維爆轟過(guò)程中爆轟波陣面附近的各種非平衡行為（包括各流體動(dòng)力學(xué)非平衡行為之間、各熱力學(xué)非平衡之間以及動(dòng)力學(xué)非平衡和熱力學(xué)非平衡之間的各種復(fù)雜影響）進(jìn)行了探索研究。

上述燃燒DBM均是單流體模型。2016年，林傳棟等提出一個(gè)二維雙流體DBM[69]。使用兩個(gè)分布函數(shù)，分別描述反應(yīng)物和產(chǎn)物，其演化方程是兩個(gè)相耦合的離散玻爾茲曼方程。該模型可用于亞聲速和超聲速燃燒現(xiàn)象的模擬，比熱比可調(diào)。原文中證明，宏觀流體模型（例如帶有化學(xué)反應(yīng)的NS方程、Fick第一和第二定律、Stefan-Maxwell擴(kuò)散定律等）所描述的行為均包含在DBM描述范圍之內(nèi)。除此之外，通過(guò)DBM還可以很方便地觀測(cè)、研究各種熱力學(xué)非平衡效應(yīng)。

下面介紹2016年張玉東等在帶有爆轟的反應(yīng)流動(dòng)理學(xué)建模與模擬方面的一項(xiàng)工作[70]。首先從理論上推導(dǎo)了一套新的流體力學(xué)方程組，如式（133）所示：

對(duì)比（133）和NS方程組可以看出，NS方程中的黏性應(yīng)力和熱流分別用文中定義的兩個(gè)非平衡量（NOMF）和（ NOEF）進(jìn)行了代替，其中對(duì)應(yīng)動(dòng)量方程中黏性應(yīng)力張量項(xiàng)Π；對(duì)應(yīng)能量方程中的熱流項(xiàng)。

為了從理論上研究各種可能的反應(yīng)率溫度依賴(lài)關(guān)系對(duì)燃燒、爆轟等過(guò)程的影響，選取式（133）所示的反應(yīng)率模型。反應(yīng)速率系數(shù)按式（134）所示的函數(shù)形式構(gòu)造。

圖21 四種反應(yīng)速率系數(shù)隨溫度變化特征圖[70]Fig. 21 Profiles of k in function of temperature for four different cases[70]

式中

通過(guò)選取如圖21所示的四種反應(yīng)速率系數(shù)，研究了四種情況下起爆過(guò)程中的非平衡效應(yīng)（包括熵產(chǎn)生特性），具體研究結(jié)果如下。

圖22和圖23分別反映了四種反應(yīng)率情形下爆轟波過(guò)程中NS黏性應(yīng)力和、 NS熱流和間的對(duì)比關(guān)系?？梢钥闯觯涸谶h(yuǎn)離爆轟波陣面處（系統(tǒng)處于熱力學(xué)平衡或近平衡態(tài)），兩組參量符合得很好；在爆轟波陣面附近（系統(tǒng)顯著偏離熱力學(xué)平衡態(tài)），兩組參量出現(xiàn)了可觀測(cè)的偏差。從建模角度來(lái)看，在由直接從玻爾茲曼方程求密度矩、動(dòng)量矩和能量矩（不做任何近似）得到的廣義NS方程組中，和Δ中包含了分布函數(shù)中偏離平衡態(tài)的高階項(xiàng)。而一般的NS方程組中Πxx和jq,x是保留了分布函數(shù)偏離平衡態(tài)的一階項(xiàng)但忽略了二階以上的高階項(xiàng)得到的。在系統(tǒng)處于平衡態(tài)附近時(shí)，高階項(xiàng)的作用很微弱，但當(dāng)系統(tǒng)顯著偏離平衡態(tài)時(shí)，這些高階項(xiàng)的作用就表現(xiàn)得比較顯著。因而和比Πxx和jq,x包含了更多信息，越是在系統(tǒng)偏離平衡態(tài)較大的區(qū)域，兩者的差別越大。

兩個(gè)非平衡特征量和熵產(chǎn)生率之間的關(guān)系如式（138 ）所示：

其中

式（139）中JS和 σ分別稱(chēng)作熵流和熵產(chǎn)生率?？梢钥闯?，（當(dāng)?shù)?、局域）熵產(chǎn)生有三個(gè)來(lái)源：NOEF（），NOMF（）和化學(xué)反應(yīng)。

圖22 非平衡量 與黏性應(yīng)力張量Π xx 對(duì)比圖[70]Fig. 22 Comparisons of non-equilibrium quantity and viscous stressΠ xx[70]

圖23 非平衡量與熱流 jq,x 對(duì)比圖[70]Fig. 23 Comparisons of non-equilibrium quantity and heat flux jq,x[70]

圖24 情形1與情形2下的三種局域熵產(chǎn)生率分布[70]Fig. 24 Three kinds of profiles of entropy productions for case1 and case2[70]

圖24和圖 25為四種情形下三種局域熵產(chǎn)生率的分布圖?？梢钥闯觯禺a(chǎn)生主要出現(xiàn)在兩個(gè)區(qū)域：在波前附近，熵產(chǎn)生主要由NOEF和NOMF構(gòu)成，后者貢獻(xiàn)大于前者；在反應(yīng)區(qū)，NOMF的貢獻(xiàn)大于NOEF，但兩者的量級(jí)均遠(yuǎn)小于化學(xué)反應(yīng)貢獻(xiàn)的熵產(chǎn)生。圖 26給出了四種情形下的全局熵產(chǎn)生率 ΔS，可以看出，化學(xué)反應(yīng)熵產(chǎn)生所占的比例遠(yuǎn)大于另外兩個(gè)方面。由于爆轟波是一個(gè)高馬赫數(shù)傳播過(guò)程，NOMF導(dǎo)致的熵產(chǎn)生大于NOEF產(chǎn)生的熵產(chǎn)生。情形3為負(fù)溫情形，負(fù)溫度系數(shù)對(duì)動(dòng)力學(xué)量的作用是降低馮·諾伊曼峰的密度、壓強(qiáng)和速度，加寬反應(yīng)區(qū)，抑制化學(xué)反應(yīng)，從而使沖擊強(qiáng)度降低，爆轟更接近于等熵過(guò)程，造成 ΔSNOEF和 ΔSNOMF降低。而反應(yīng)率的降低，使爆轟過(guò)程中化學(xué)反應(yīng)的準(zhǔn)靜態(tài)程度增加，使得 ΔSCHEM降低。更多內(nèi)容可參見(jiàn)文獻(xiàn)[70]。

圖25 情形3與情形4下的三種局域熵產(chǎn)生率分布[70]Fig. 25 Three kinds of profiles of entropy productions for case3 and case4[70]

圖26 四種反應(yīng)速率情形下的全局熵產(chǎn)生率[70]Fig. 26 Three kinds of profiles of global entropy productions for four cases[70]

林傳棟后期與清華大學(xué)合作導(dǎo)師羅開(kāi)紅教授等人一起，將復(fù)雜多相流動(dòng)的DBM建模研究進(jìn)一步推向深入。2017年發(fā)表了一個(gè)可用于預(yù)混、非預(yù)混或部分預(yù)混的非平衡反應(yīng)流的多組分DBM[71]，模型適用于亞聲速和超聲速流動(dòng)，可用于化學(xué)反應(yīng)流或不帶化學(xué)反應(yīng)的流動(dòng)，可模擬外力項(xiàng)存在或不存在的情況。2018年發(fā)表一個(gè)用于可壓縮放熱反應(yīng)流的多松弛時(shí)間DBM，模型具有可調(diào)的比熱比和普朗特?cái)?shù)[72]；使用DBM研究了爆轟波面附近的動(dòng)力學(xué)和熱力學(xué)非平衡（HTNE）效應(yīng)[73]，考察了化學(xué)反應(yīng)物和化學(xué)產(chǎn)物的全局和當(dāng)?shù)豀TNE特征，并分析了它們分布函數(shù)的主要特征。

2019年，張玉東等提出了一個(gè)用于模擬爆轟的一維DBM[74]?；谛履Ｐ?，對(duì)反應(yīng)速率負(fù)溫度系數(shù)對(duì)爆轟的影響進(jìn)行了進(jìn)一步研究，發(fā)現(xiàn)了在負(fù)溫度系數(shù)條件下，會(huì)發(fā)生周期性出現(xiàn)兩個(gè)波面的反常爆轟現(xiàn)象。2020年，林傳棟等研究了初始擾動(dòng)幅值和波長(zhǎng)以及化學(xué)反應(yīng)放熱對(duì)具有非平衡效應(yīng)的非穩(wěn)定爆轟演化的影響[75]。

5 小結(jié)與說(shuō)明

對(duì)于多相復(fù)雜流體系統(tǒng)的研究，宏觀、微觀和介觀都有相應(yīng)的模型和描述方法。微觀方法由于時(shí)間和空間尺度的限制，目前主要用于從原子（或分子）層次進(jìn)行一些機(jī)理性的研究[76-78]。而只考慮一階離散效應(yīng)和非平衡效應(yīng)的NS方程組只能對(duì)系統(tǒng)中大尺度和緩變行為進(jìn)行較好的描述，對(duì)于流體系統(tǒng)中存在的一些銳利界面和快變模式的描述則不能滿足要求。從物理描述能力角度，DBM描述的動(dòng)理學(xué)性質(zhì)比玻爾茲曼方程要少（研究視角相對(duì)較窄），但比NS等宏觀流體力學(xué)方程組要多（研究視角相對(duì)較寬）；要保的動(dòng)理學(xué)性質(zhì)可以根據(jù)離散或非平衡程度的研究需求來(lái)選取。所以，DBM為多相復(fù)雜流體系統(tǒng)中各種空間尺度和時(shí)間尺度上的各類(lèi)非平衡行為的建模與模擬提供了一條簡(jiǎn)潔有效的思路和方法。

在復(fù)雜物理場(chǎng)分析方面，DBM在(f-feq)非守恒矩張開(kāi)的相空間及其子空間中描述系統(tǒng)偏離平衡的方式和幅度，描述相關(guān)方面的動(dòng)理學(xué)輸運(yùn)性質(zhì)；進(jìn)一步，借助兩點(diǎn)間距離描述兩個(gè)非平衡狀態(tài)的差別，借助其倒數(shù)描述這兩個(gè)非平衡狀態(tài)的相似度；借助一段時(shí)間內(nèi)兩點(diǎn)間距離的平均值描述兩個(gè)動(dòng)理學(xué)過(guò)程之間的差別，借助其倒數(shù)描述這兩個(gè)動(dòng)理學(xué)過(guò)程之間的相似度；借助非守恒矩及其演化來(lái)描述復(fù)雜流動(dòng)過(guò)程中引發(fā)熵增的主要因素及其相對(duì)的重要性；等等。示蹤粒子的引入，使得可以在單流體理論框架下對(duì)混合過(guò)程中的粒子來(lái)源進(jìn)行識(shí)別與追蹤；部分或全部示蹤粒子在其速度狀態(tài)空間的分布與演化為復(fù)雜流場(chǎng)的研究提供一個(gè)嶄新的視角。

系統(tǒng)內(nèi)不同非平衡行為特征之間的空間關(guān)聯(lián)、時(shí)間關(guān)聯(lián)、時(shí)空關(guān)聯(lián)、競(jìng)爭(zhēng)與協(xié)作等是多相復(fù)雜流動(dòng)過(guò)程中特征、機(jī)制與規(guī)律研究的重要內(nèi)容。這些以前知之甚少的“介尺度”非平衡行為特征，蘊(yùn)含著大量有待開(kāi)發(fā)的物理功能。隨著研究逐步深入，DBM將在連續(xù)介質(zhì)建模失效或物理功能不足、而微觀MD方法因適用尺度受限而無(wú)能為力的各類(lèi)復(fù)雜流動(dòng)研究中發(fā)揮更加重要的作用。

最后需要說(shuō)明的是，與NS模型一樣，DBM模型是粗粒化描述系統(tǒng)行為的一種理論模型，它對(duì)系統(tǒng)動(dòng)理學(xué)性質(zhì)描述的廣度與深度根據(jù)具體問(wèn)題研究需求而調(diào)整。獲得了DBM模型，跟獲得了NS模型一樣，還需要選擇合適的數(shù)值計(jì)算方法，才能進(jìn)行數(shù)值實(shí)驗(yàn)研究。從物理問(wèn)題研究需求的角度，可以選用滿足物理問(wèn)題研究需求的且當(dāng)前硬軟件環(huán)境允許的任何一種數(shù)值計(jì)算方法。對(duì)數(shù)值計(jì)算方法感興趣的讀者可參閱更專(zhuān)業(yè)的文獻(xiàn)著作。希望這篇理論物理視角的多相流研究綜述能與探討多相流實(shí)驗(yàn)、多相流計(jì)算方法的研究論文和綜述論文形成互補(bǔ)與借鑒。

致謝：感謝甘延標(biāo)、陳鋒、林傳棟、張玉東、賴(lài)惠林、李華、應(yīng)陽(yáng)君、謝侃、陶建軍、馬天宇、劉枝朋、張戈、張德佳、單奕銘、孫光蘭、陳鋮、李晗蔚等在研究過(guò)程中不同形式的有益討論。