多相流系統的離散玻爾茲曼研究進展

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

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

0 引 言

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

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

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

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

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

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

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

具體來說，本文要重點介紹的離散玻爾茲曼建模，是非平衡統計物理學粗粒化建模方法在流體力學領域的具體應用，是理論物理范疇下的模型構建方法。它根據研究需求，抓主要矛盾，選取一個視角，研究系統的一組動理學性質，因而它要求描述這組性質的動理學矩其計算結果不能因為模型簡化而改變。除了分布函數的守恒矩（質量、動量和能量），DBM同時關注部分關系最密切的非守恒矩的時空演化，從一個更寬的視角來考察和描述系統。除了為實現模擬而進行的抓主要矛盾和粗粒化建模，DBM同時關注模擬之后的復雜物理場分析（復雜物理場分析，也需要通過建模來提取更多有價值的信息），它本身自帶一套復雜物理場分析方法或技術[27-28]。

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

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

1 流體系統的微觀、介觀與宏觀建模

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

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

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

2 粗粒化建模：從格子氣模型到離散玻爾茲曼

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

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

物理圖像簡單、直觀、易于編程等優點，使得LGM在粗略模擬流體行為（特別是復雜流體行為）方面發揮著重要作用。NS等偏微分方程的數值解法和非平衡復雜流動的粗粒化物理建模是LGM的兩個發展方向。這兩類LGM目標不同，所以構建規則不同，使用的判據也不同。

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

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

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

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

3 離散玻爾茲曼理論框架

3.1 玻爾茲曼方程簡介

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

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

其中，F(q,p)為相空間的分布函數。

對于保守力學體系，哈密頓量等于系統總能量，等于常數，即：

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

由哈密頓正則方程：

即可得到非平衡統計力學的基本方程—劉維（Liouville）方程：

上式又經常寫為：

線性算法L定義為：

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

為便于描述，引入簡記：

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

其中s是 0到N之間的整數。s粒子約化分布函數可以理解為在某一時刻在相空間中s個點同時發現s個粒子的聯合概率。對于N個質量為m的全同粒子組成的經典系統，假設系統的體積是有限的。系統的哈密頓量可做如下分解：

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

劉維方程可寫為：

根據系統內粒子數守恒和粒子在分布函數中的對稱性，得：

上式可進一步寫為：

或

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

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

其中

如果我們引入五個假設對研究系統進行約束：

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

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

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

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

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

由式（18）、式（19），再根據角動量守恒，碰撞前后粒子“瞄準距離”相等，可推導得到：

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

從BBGKY方程鏈出發，經過一系列假設，我們得到了玻爾茲曼方程。下面，我們從物理圖像的角度對玻爾茲曼方程進行比較直觀的解釋。

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

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

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

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

圖1 質心坐標系中二體碰撞的分子運動軌跡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

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

設 σ為單位立體角所對應的截面積，則有：

可得：

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

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

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

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

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

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

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

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

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

所以，在時間dt內相空間體積 drdv中，由于和速度為v1的 分子發生碰撞而導致的速度為v的 分子數目的變化為對速度空間v1積分，則可以得到在時間 dt內 相空間體積 drdv中，由于分子間的碰撞而導致的速度為v的分子數目的變化：

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

3.2 Chapman-Enskog多尺度分析

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

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

其中，

為麥克斯韋分布，

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

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

從擾動論角度看，在Chapman-Enskog多尺度分析中，Kn對應的是施加給系統的擾動。Chapman-Enskog多尺度展開收斂的情形對應系統可以回到平衡態的情形；如果擾動過強，導致泰勒展式發散，則意味著該擾動有可能引發新的結構或模式。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

引入BGK近似后成為：

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

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

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

3.5.1 基于范德瓦爾斯狀態方程的DBM

作為一個粗粒化模型和諸多問題的研究起點，玻爾茲曼方程描述的是理想氣體系統。對于理想氣體系統，在確定的溫度和壓強下，只能有一個密度，系統基態永遠是單一態，是不可能發生相變的。因此，在涉及相變的多相流問題研究中，一個關鍵的問題就是分子間作用力的引入。

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

ε和 σ是兩個待定參數，對于不同氣體分子具有不同的值。在經典力學中，N個質量為m的全同粒子組成的系統的哈密頓量為：

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

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

由配分函數，可得自由能F的表達式：

以及壓強P、內能密度e和單位粒子熵s的表達式：

其中下標表示求導時的不變量。

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

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

再由式（60）可得VDW氣體狀態方程：

或者寫為：

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

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

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

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

3.5.2 基于C-S狀態方程的DBM

為了更準確地描述硬球間的相互作用，Carnhan和Starling在VDW狀態方程的基礎上，對斥力項進行了修正，提出Carnhan-Starling（C-S）狀態方程[48]：

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

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

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

Λ為壓強張量，I是單位張量，K是表面張力系數，ζ是體黏性系數。

系統總的能量密度為：

3.5.3 基于分子間作用勢的DBM

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

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

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

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

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

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

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

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

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

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

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

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

其中，

F′表達式中第一項與狀態方程有關，第二項與表面張力有關。

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

3.6 離散玻爾茲曼建模：含化學反應情形

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

其中：f為分布函數；feq為對應的平衡態分布函數；τ為熱力學弛豫時間；v為 分子速度；C為化學反應項，描述由于化學反應而引起的分布函數f的變化率。

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

其中，f為不考慮化學反應時系統的平衡態分布函數，f*eq=f*eq(ρ,u,T*)為考慮化學反應貢獻后系統的平衡態分布函數。

若化學反應不可逆，則反應進程可由下面的反應率方程來描述：

其中Q為單位質量的反應物發生反應可以釋放的熱量。由式（94）～（96）可以看出，化學反應項C的強弱不僅取決于反應進程，還受到反應放熱Q的影響。即使化學反應速率很快，當Q很小時，化學反應項C的貢獻也可能很小。

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

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

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

3.7 離散玻爾茲曼建模：多介質情形

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

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

其中mσ是（摩爾或）粒子質量。混合物局域的總質量密度 ρ、（摩爾或）粒子數密度n、 總動量Jα和平均流速uα分別由下面的關系式得到：

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

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

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

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

則混合物的壓強：

同時，定義介質 σ 基于自己流速的壓強：

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

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

與此對應，動理學矩空間的分布函數也有三種：

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

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

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

在含示蹤粒子的系統中，我們可以使用斯托克斯數（Stokes number,St）來描述該粒子的動力學狀態：

其中 τP是粒子的特征弛豫時間，u0是當地的流動速度，l0是特征長度（通常選取顆粒的直徑）。當St＞1時，粒子將由于慣性脫離當地流動而運動，特別是流速變化劇烈的情況下；當St?1 時，粒子緊緊貼著流線運動。通過調整St到足夠小的數量級，則該粒子能夠作為流場的示蹤粒子使用。假設粒子的弛豫時間τP接近0，則其慣性完全可以忽略，其運動速度瞬間達到當地流速，因而完全跟隨流體運動。

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

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

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

其中a為示蹤粒子在流場中所處位置的微元，D表示對示蹤粒子速度具有影響的流場積分區域，?為Dirac函數，在模擬中通常使用其離散形式 ψ代替。第k個示蹤粒子的速度將根據它所處的位置以及當地的流動速度決定，例如經過時間間隔 Δt，示蹤粒子速度將從變化為為了更新點 顆 粒 的 位 置，可使用四階龍格-庫塔格式求解離散顆粒的運動方程。

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

ψ 函數可以被分解為兩個部分：

其中，i,j為網格點編號。張戈等在文獻[51]中所應用的權重函數 φ為：

據此，示蹤粒子從流場中獲得了自己的運動速度。

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

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

受控熱核聚變是解決人類能源問題的關鍵，其中等離子體的運動處于高溫高密高壓的極端環境，同時涉及到多時空尺度的多場耦合及帶電粒子間的碰撞輸運過程，是人們關注的重點。

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

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

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

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

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

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

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

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

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

如前面所述，將離子和電子分別作為不同的流體介質，可以使用雙流體玻爾茲曼模型。當然，還存在其他不同粗粒化程度的物理建模形式，如采用分布函數描述離子的行為特征，而假設電子始終處于局域熱力學平衡態或對電子采用流體描述方法。

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

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

4 多相流機理研究：基于DBM

作為系統行為粗粒化描述的一種物理模型構建方法，DBM本身并不包含具體的數值離散格式。獲得了DBM模型，跟獲得了NS模型一樣，在數值模擬中還需根據問題特點選取合適的數值方法。從物理問題研究的角度，可以選用滿足當次物理問題研究需求的任何一種離散格式。在下面物理結果的原始文獻中使用的離散格式只是滿足當次研究需求的多種離散格式之一。

4.1 流體不穩定性方面

流體界面不穩定性（經常簡稱流體不穩定性）與物質混合現象廣泛存在于自然界、慣性約束核聚變和武器物理等領域。慣性約束核聚變和武器物理等領域重點關注的流體不穩定性主要有三種：RM不穩定性（Richtmyer-Meshkov instability, RMI）、RT不穩定性（Rayleigh-Taylor instability, RTI）和KH不穩定性（Kelvin-Helmholtz instability, KHI）。

流體不穩定性系統通常具有以下特點：系統的本身可能是宏觀尺度的，但其內部存在大量的中間尺度的空間結構和動理學模式；這些結構和模式的存在與演化極大地影響著系統的物理性能和功能。系統內部往往具有大量的界面，包括物質界面和力學界面，系統內部的受力和響應過程非常復雜。對于這些系統的研究，NS方程只能描述大尺度緩變行為，而對于一些銳利的邊界，流體的平均分子間距相對于界面尺度已經不再是可以忽略的小量，基于連續介質假設的NS方程受到挑戰。在一些快變流動模式描述方面，流體系統趨于熱力學平衡的弛豫時間相對于該流動的特征時間來講，不再是可以忽略的小量，熱力學非平衡效應較強，只考慮一階熱力學非平衡效應的NS描述不再能滿足需求。

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

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

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

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

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

圖4 RT與RM不穩定性演化過程中系統內不同類型的非平衡強度與密度、溫度、流速不均勻度之間關聯程度的對比[56, 62]（更多細節參見文獻[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不穩定性演化過程中熱傳導對相關度的影響[62]（更多細節參見文獻[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不穩定性系統，在RM不穩定性系統中，沖擊波的透射和反射使得相關度演化曲線出現“跳變”（見圖4）。在RT不穩定性與RM不穩定性共存的系統中，重力加速度和沖擊波強度之間合作、競爭，共同主導著界面演化與物質和能量混合過程。在RM不穩定性主導的情形初始擾動界面會出現反轉（見圖6-圖7，更多細節參見文獻[62]）。

RT不穩定性和RM不穩定性分別主導的參數區間如圖8所示。圖9展示了全局平均TNE強度、全局平均密度梯度、全局平均NOMF強度以及全局平均NOEF強度隨時間的演化；重力加速度對非平衡行為特征的影響表現出階段性（見圖9）；隨著馬赫數的增加，系統的整體非平衡強度指數增加，而溫度不均勻度與NOEF的相關度指數衰減（見圖10）。圖8～圖10為數值模擬結果，更多細節參見文獻[62]。

圖6 RTI與RMI共存系統界面反轉現象出現機制的示意圖[62]Fig. 6 The schematic diagram of the collaboration and competition relations between RM and RT instability[62]

圖7 RTI與RMI共存系統界面反轉現象出現與否的數值模擬密度云圖[62]Fig. 7 Numerical simulation density nephogram of interface inversion in RTI and RMI coexistence system[62]

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

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

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

圖11為重力加速度g=0.005、 切向速度差u0=0.05的耦合RTKHI系統的溫度圖像、 圖靈斑圖以及u0=0.05的純KHI的圖靈斑圖。可以看出，和單純的KHI系統相比，耦合RTKHI系統具有傾斜的、非對稱的氣泡、尖釘以及蘑菇狀結構，單純的KHI系統的演化大大落后于耦合RTKHI系統的演化。在此情況下，RTI扮演了一個主要角色，切向速度的存在主要破壞了RTI結構的對稱性。更多情況的算例及物理分析可參見文獻[63]。

圖10 馬赫數對RTI與RMI共存系統TNE強度、宏觀量梯度與非平衡特征相關度的影響[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系統（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、擾動幅值的增長率dA/dt、形態學邊界長度L隨時間的演化。對于KHI在早期階段占主導、RTI在后期階段占主導的情形，演化過程可以粗略地劃分為兩個階段，分別為剪切主導階段和浮力主導階段。由剪切主導階段到浮力主導階段的過渡狀態稱為過渡點。從圖12中可以看到，在開始階段，LRTKHI首先呈指數增長，然后呈線性增長（虛線所示）。LRTKHI線性增長結束的點近似等于RTKHI系統幅值增長率最小值點以及相關KHI系統的幅值最大值點。從這一時刻開始，RTI開始扮演一個主要角色。這一時刻稱為過渡點，在圖中使用豎直虛線和圓標記。因此，LRTKHI線性增長結束的點可以作為區別兩個階段的幾何判據。剪切率越大，過渡點越早出現。

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

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

圖12 振幅 A、振幅增長率d A/dt 和形態邊界長度 L 與時間 t 的關系[63]Fig. 12 Perturbation amplitude A, growth rate d A/dt, and morphological boundary length L vs time t [63]

圖13 根據溫度場的形態分析早期主要機制[63]Fig. 13 Morphological analysis of the main mechanism in the early stage[63]

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

圖14 RTKHI系統(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系統的非平衡特征在早期主要機理判斷和過渡點捕獲中的應用[63]Fig. 15 The applications of non-equilibrium characteristics in the early main mechanism judgment and transition point capture[63]

由于流體界面不穩定性系統的DBM建模與模擬尚處于起步階段，所以以上研究均是基于理想氣體模型的。表面張力、材料強度等對流體不穩定性演化過程的影響，以及更加實際系統中流體不穩定性的DBM建模與模擬有待進一步研究。

4.2 相分離方面

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

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

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

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

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

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

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

圖18 黏性對熱相分離過程中SD階段持續時間和熵產生速率的影響[44]Fig. 18 Effect of viscosity on duration of SD stage and entropy generation rate in thermal phase separation[44]

圖20給出了不同熱流、黏性和表面張力情形下，兩種熵產生機制之間的競爭與協作關系。發現，在熱流和黏性的作用下，兩種熵產生機制之間存在競爭關系；在表面張力作用下，兩種熵產生機制之間存在協作關系。

圖19 表面張力對熱相分離過程中SD階段持續時間和熵產生速率的影響[44]Fig. 19 Effect of surface tension on duration of SD stage and entropy generation rate in thermal phase separation process[44]

圖20 兩種熵產生速率幅值之間的關系曲線，圖中箭頭指向熱傳導系數（1/Pr）、黏性系數(τ)和表面張力系數（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 化學反應流方面

近年來，DBM在燃燒領域的應用取得了一系列的進展。

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

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

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

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

為了從理論上研究各種可能的反應率溫度依賴關系對燃燒、爆轟等過程的影響，選取式（133）所示的反應率模型。反應速率系數按式（134）所示的函數形式構造。

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

式中

通過選取如圖21所示的四種反應速率系數，研究了四種情況下起爆過程中的非平衡效應（包括熵產生特性），具體研究結果如下。

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

兩個非平衡特征量和熵產生率之間的關系如式（138 ）所示：

其中

式（139）中JS和 σ分別稱作熵流和熵產生率。可以看出，（當地、局域）熵產生有三個來源：NOEF（），NOMF（）和化學反應。

圖22 非平衡量 與黏性應力張量Π 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下的三種局域熵產生率分布[70]Fig. 24 Three kinds of profiles of entropy productions for case1 and case2[70]

圖24和圖 25為四種情形下三種局域熵產生率的分布圖。可以看出，熵產生主要出現在兩個區域：在波前附近，熵產生主要由NOEF和NOMF構成，后者貢獻大于前者；在反應區，NOMF的貢獻大于NOEF，但兩者的量級均遠小于化學反應貢獻的熵產生。圖 26給出了四種情形下的全局熵產生率 ΔS，可以看出，化學反應熵產生所占的比例遠大于另外兩個方面。由于爆轟波是一個高馬赫數傳播過程，NOMF導致的熵產生大于NOEF產生的熵產生。情形3為負溫情形，負溫度系數對動力學量的作用是降低馮·諾伊曼峰的密度、壓強和速度，加寬反應區，抑制化學反應，從而使沖擊強度降低，爆轟更接近于等熵過程，造成 ΔSNOEF和 ΔSNOMF降低。而反應率的降低，使爆轟過程中化學反應的準靜態程度增加，使得 ΔSCHEM降低。更多內容可參見文獻[70]。

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

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

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

2019年，張玉東等提出了一個用于模擬爆轟的一維DBM[74]。基于新模型，對反應速率負溫度系數對爆轟的影響進行了進一步研究，發現了在負溫度系數條件下，會發生周期性出現兩個波面的反常爆轟現象。2020年，林傳棟等研究了初始擾動幅值和波長以及化學反應放熱對具有非平衡效應的非穩定爆轟演化的影響[75]。

5 小結與說明

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

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

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

最后需要說明的是，與NS模型一樣，DBM模型是粗粒化描述系統行為的一種理論模型，它對系統動理學性質描述的廣度與深度根據具體問題研究需求而調整。獲得了DBM模型，跟獲得了NS模型一樣，還需要選擇合適的數值計算方法，才能進行數值實驗研究。從物理問題研究需求的角度，可以選用滿足物理問題研究需求的且當前硬軟件環境允許的任何一種數值計算方法。對數值計算方法感興趣的讀者可參閱更專業的文獻著作。希望這篇理論物理視角的多相流研究綜述能與探討多相流實驗、多相流計算方法的研究論文和綜述論文形成互補與借鑒。

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