許愛國，陳 杰，宋家輝，陳大偉，陳志華

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 離散玻爾茲曼理論框架

3.1 玻爾茲曼方程簡介

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

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

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

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

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

由哈密頓正則方程：

即可得到非平衡統(tǒng)計力學(xué)的基本方程—劉維（Liouville）方程：

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

線性算法L定義為：

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

為便于描述，引入簡記：

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

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

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

劉維方程可寫為：

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

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

或

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

BBGKY方程鏈的第一、第二兩個方程為：

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

可得：

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

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

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

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

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

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

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

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

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

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

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

3.2 Chapman-Enskog多尺度分析

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

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

其中，

為麥克斯韋分布，

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

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

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

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

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

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

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

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

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

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

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

離散速度的具體取法取決于要保的不變量、基本的守恒性和必要的對稱性。可見，DBM給出的是對離散速度的物理約束，并不包含具體的離散格式。

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

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

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

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

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

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

3.4 離散玻爾茲曼建模：含外場力情形

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

引入BGK近似后成為：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

或者寫為：

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

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

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

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

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

為了更準(zhǔn)確地描述硬球間的相互作用，Carnhan和Starling在VDW狀態(tài)方程的基礎(chǔ)上，對斥力項進(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 基于分子間作用勢的DBM

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

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

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

則恩斯柯格碰撞算符可寫為：

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

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

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

則可以將簡化的恩斯柯格方程寫為：

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

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

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

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

其中，

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

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

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

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

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

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

其中，f為不考慮化學(xué)反應(yīng)時系統(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)率方程來描述：

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

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

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

對于多弛豫時間模型，Chapman-Enskog多尺度分析按如下方式展開：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 RT與RM不穩(wěn)定性演化過程中系統(tǒng)內(nèi)不同類型的非平衡強(qiáng)度與密度、溫度、流速不均勻度之間關(guān)聯(lián)程度的對比[56, 62]（更多細(xì)節(jié)參見文獻(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)定性演化過程中熱傳導(dǎo)對相關(guān)度的影響[62]（更多細(xì)節(jié)參見文獻(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)“跳變”（見圖4）。在RT不穩(wěn)定性與RM不穩(wěn)定性共存的系統(tǒng)中，重力加速度和沖擊波強(qiáng)度之間合作、競爭，共同主導(dǎo)著界面演化與物質(zhì)和能量混合過程。在RM不穩(wěn)定性主導(dǎo)的情形初始擾動界面會出現(xiàn)反轉(zhuǎn)（見圖6-圖7，更多細(xì)節(jié)參見文獻(xiàn)[62]）。

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

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

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

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

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

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

圖12 振幅 A、振幅增長率d A/dt 和形態(tài)邊界長度 L 與時間 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ù)溫度場的形態(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ī)理，并且計算結(jié)果和形態(tài)學(xué)邊界長度L是一致的。由圖15（e）～15（f）可以看出，D3,1展示了和邊界長度L相似的行為。交界面長度L越大，D3,1越大。D3,1線性增長的結(jié)束點(diǎn)，可以用來作為區(qū)分從KHI主導(dǎo)到RTI主導(dǎo)兩個階段的物理判據(jù)。具體的物理分析可參見文獻(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ī)理判斷和過渡點(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)度等對流體不穩(wěn)定性演化過程的影響，以及更加實(shí)際系統(tǒng)中流體不穩(wěn)定性的DBM建模與模擬有待進(jìn)一步研究。

4.2 相分離方面

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

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

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

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

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

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

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

圖18 黏性對熱相分離過程中SD階段持續(xù)時間和熵產(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ī)制之間的競爭與協(xié)作關(guān)系。發(fā)現(xiàn)，在熱流和黏性的作用下，兩種熵產(chǎn)生機(jī)制之間存在競爭關(guān)系；在表面張力作用下，兩種熵產(chǎn)生機(jī)制之間存在協(xié)作關(guān)系。

圖19 表面張力對熱相分離過程中SD階段持續(xù)時間和熵產(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)流方面

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

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

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

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

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

為了從理論上研究各種可能的反應(yīng)率溫度依賴關(guān)系對燃燒、爆轟等過程的影響，選取式（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]

式中

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

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

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

其中

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

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

圖23 非平衡量與熱流 jq,x 對比圖[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)生率的分布圖?？梢钥闯?，熵產(chǎn)生主要出現(xiàn)在兩個區(qū)域：在波前附近，熵產(chǎn)生主要由NOEF和NOMF構(gòu)成，后者貢獻(xiàn)大于前者；在反應(yīng)區(qū)，NOMF的貢獻(xiàn)大于NOEF，但兩者的量級均遠(yuǎn)小于化學(xué)反應(yīng)貢獻(xiàn)的熵產(chǎn)生。圖 26給出了四種情形下的全局熵產(chǎn)生率 ΔS，可以看出，化學(xué)反應(yīng)熵產(chǎn)生所占的比例遠(yuǎn)大于另外兩個方面。由于爆轟波是一個高馬赫數(shù)傳播過程，NOMF導(dǎo)致的熵產(chǎn)生大于NOEF產(chǎn)生的熵產(chǎn)生。情形3為負(fù)溫情形，負(fù)溫度系數(shù)對動力學(xué)量的作用是降低馮·諾伊曼峰的密度、壓強(qiáng)和速度，加寬反應(yīng)區(qū)，抑制化學(xué)反應(yīng)，從而使沖擊強(qiáng)度降低，爆轟更接近于等熵過程，造成 ΔSNOEF和 ΔSNOMF降低。而反應(yīng)率的降低，使爆轟過程中化學(xué)反應(yīng)的準(zhǔn)靜態(tài)程度增加，使得 ΔSCHEM降低。更多內(nèi)容可參見文獻(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)師羅開紅教授等人一起，將復(fù)雜多相流動的DBM建模研究進(jìn)一步推向深入。2017年發(fā)表了一個可用于預(yù)混、非預(yù)混或部分預(yù)混的非平衡反應(yīng)流的多組分DBM[71]，模型適用于亞聲速和超聲速流動，可用于化學(xué)反應(yīng)流或不帶化學(xué)反應(yīng)的流動，可模擬外力項存在或不存在的情況。2018年發(fā)表一個用于可壓縮放熱反應(yīng)流的多松弛時間DBM，模型具有可調(diào)的比熱比和普朗特數(shù)[72]；使用DBM研究了爆轟波面附近的動力學(xué)和熱力學(xué)非平衡（HTNE）效應(yīng)[73]，考察了化學(xué)反應(yīng)物和化學(xué)產(chǎn)物的全局和當(dāng)?shù)豀TNE特征，并分析了它們分布函數(shù)的主要特征。

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

5 小結(jié)與說明

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

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

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

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

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