格子玻爾茲曼方法

格子玻爾茲曼方法（Lattice Boltzmann methods，簡稱LBM）是一類用於流體動力學模擬的計算流體力學（CFD）技術。與傳統的基於納維-斯托克斯方程的求解器不同，LBM 並非直接離散化宏觀的流體控制方程，而是建立在介觀（mesoscopic）尺度的動理學理論之上。該方法通過模擬虛擬粒子群在離散格子（lattice）上的傳播（streaming）和碰撞（collision）過程，從統計物理的角度重現流體的宏觀行為。

LBM獨特的介觀視角使其在處理複雜邊界條件、多相流、多孔介質流動以及多物理場耦合問題時具有顯著優勢。該方法的演算法具有高度的局部性，極易於在現代大規模並行計算機（包括GPU）上實現高效運算。

## 歷史沿革

LBM 的發展經歷了從離散的格子氣模型到基於連續分布函數的動力學方程的理論演進。

### 格子氣自動機

LBM 的直接前身是20世紀80年代興起的格子氣自動機（Lattice Gas Automata，簡稱LGA）。LGA 採用極其簡化的粒子模型來模擬流體：虛擬粒子在規則的格子上移動，當多個粒子相遇時，它們根據預設的碰撞規則（如動量守恆）改變運動方向。1986年提出的FHP模型首次證明，在六邊形格子上運行的簡單布爾粒子系統（粒子存在或不存在）能夠在宏觀尺度上重現二維納維-斯托克斯方程描述的流體行為，這標誌著一個重要的理論突破。

然而，LGA存在一些固有缺陷，限制了其在工程中的應用：

- 統計噪聲：由於粒子狀態是離散的布爾值，需要進行大量的時空平均才能獲得平滑的宏觀物理量。
- 缺乏伽利略不變性：宏觀方程中的對流項包含非物理因子，導致模擬結果在不同的參考系下不一致。
- 各向異性：在方形格子上，模型無法保證流體性質的各向同性。

### 從LGA到LBE的演化

為了克服LGA的統計噪聲問題，1988年，研究者提出用實數值的粒子分布函數（distribution function）$f_i$來代替布爾變量，代表在特定方向$i$上粒子出現的機率或系綜平均數量。這一改進徹底消除了統計噪聲，標誌著格子玻爾茲曼方程（Lattice Boltzmann Equation，簡稱LBE）的誕生。

### LBGK模型的確立

LBE發展史上的一個關鍵里程碑是1992年LBGK模型的提出。該模型引入了BGK近似，將複雜的粒子碰撞過程簡化為一個朝向局部平衡態的單時間鬆弛過程。LBGK模型不僅極大地簡化了演算法，而且通過精心設計的平衡態分布函數，成功修正了LGA中的伽利略不變性等物理性錯誤。由於其簡潔、高效和物理完備性，LBGK至今仍是LBM中最常用和基礎的模型。

### 理論基礎的完備化

儘管LBM最初源於LGA，但其堅實的理論基礎直到1997年才被完全確立。研究證明，LBE可以不依賴於LGA，直接通過對連續玻爾茲曼方程在速度空間和時間上進行系統性離散而推導出來。這一工作將LBM從一個啟發式模型提升為一個具有堅實動理學理論基礎的數值方法，並為其後續發展（如構造高階模型和非均勻網格模型）鋪平了道路。

## 核心原理

LBM的核心思想是將流體視為大量虛擬粒子組成的系統，這些粒子在規則的格子上進行「遷移」和「碰撞」兩個基本過程。

### 粒子分布函數與格子結構

LBM不追蹤單個粒子的運動，而是追蹤粒子分布函數 $f_i(\mathbf{x}, t)$。該函數表示在時間 $t$、空間節點 $\mathbf{x}$ 處，速度方向為 $\mathbf{e}_i$ 的粒子密度。

空間被離散為均勻的格子（lattice），速度也被離散為一小組固定的向量 $\{\mathbf{e}_i\}$。這些離散速度向量必須具有足夠的對稱性，以保證宏觀流體行為的各向同性。常用的格子模型以「DnQm」命名，其中n是空間維度，m是離散速度的數量。

- D2Q9：二維空間中最常用的模型，包含1個靜止粒子和8個分別指向相鄰及次相鄰節點的粒子。
- D3Q15, D3Q19, D3Q27：用於三維空間模擬的常見模型。

### 演化方程：碰撞與遷移

LBM的演化過程通過一個離散的動力學方程——格子玻爾茲曼方程（LBE）來描述。對於LBGK模型，其方程形式為：

$$f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) + \Omega_i(\mathbf{x}, t)$$

其中，$\Omega_i$ 是碰撞算子，描述了節點上粒子因碰撞而發生的分布函數變化。對於LBGK模型，該算子表示為：

$$\Omega_i = -\frac{1}{\tau} (f_i - f_i^{\text{eq}})$$

這裡的 $\tau$ 是一個無量綱的鬆弛時間，它控制著系統趨向局部平衡態 $f_i^{\text{eq}}$ 的速率，並直接決定流體的宏觀粘度。

整個演化過程可分為兩個步驟：

- 碰撞（Collision）：在每個格點上，粒子分布函數 $f_i$ 向局部平衡態分布函數 $f_i^{\text{eq}}$ 鬆弛。這一過程是完全局部的，即每個節點的計算不依賴於其他節點。
- 遷移（Streaming）：碰撞後的粒子分布函數沿著各自的速度方向 $\mathbf{e}_i$ 移動到相鄰的格點上。這是一個簡單的資訊傳遞過程。

### 宏觀量的計算

流體的宏觀物理量（如密度 $\rho$ 和速度 $\mathbf{u}$）可以通過對分布函數求矩來獲得：

- 密度：$\rho = \sum_i f_i$
- 動量：$\rho\mathbf{u} = \sum_i f_i \mathbf{e}_i$

這些宏觀量又反過來決定了下一時刻的局部平衡態 $f_i^{\text{eq}}$，從而將微觀動力學與宏觀流體行為聯繫起來。

## 理論基礎

### 平衡態分布函數

平衡態分布函數 $f_i^{\text{eq}}$ 是LBM理論的核心。它代表了在給定宏觀密度和速度下，系統達到熱力學平衡時粒子的最可幾分布。其具體形式通常是麥克斯韋-玻爾茲曼分布在低馬赫數下的泰勒級數展開。對於D2Q9模型，其標準形式為：

$$f_i^{\text{eq}} = w_i \rho \left[ 1 + \frac{\mathbf{e}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{e}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u}^2}{2c_s^2} \right]$$

其中 $w_i$ 是與格子模型相關的權重係數，$c_s$ 是格子聲速，它定義了模型中資訊傳播的速度上限。

### 從LBE到納維-斯托克斯方程

LBM之所以能夠模擬流體，是因為通過查普曼-恩斯科格（Chapman-Enskog）展開可以證明，其在宏觀尺度和低頻極限下等價於納維-斯托克斯方程。這種多尺度分析顯示：

- 在零階近似下，LBE恢復了描述無粘流體的歐拉方程。
- 在一階近似下，LBE引入了粘性耗散項，恢復了完整的納維-斯托克斯方程。

通過該展開，可以建立起介觀參數（鬆弛時間 $\tau$）與宏觀物理量（運動粘度 $\nu$）之間的直接關係：

$$\nu = c_s^2 \left( \tau - \frac{1}{2} \right) \Delta t$$

這個關係式揭示了LBM的一個重要特性：流體的粘度可以通過調節鬆弛時間 $\tau$ 來直接控制。同時，為了保證粘度為正（物理上有意義），$\tau$ 必須大於0.5，這是LBGK模型數值穩定性的基本要求。

## 高級模型與擴展

### 多鬆弛時間（MRT）模型

標準的LBGK模型使用單一的鬆弛時間 $\tau$，這導致了一些局限性，例如固定的普朗特數和在高雷諾數下數值穩定性較差。為了克服這些問題，多鬆弛時間（MRT）模型被提出。

MRT模型的核心思想是在「矩空間」（moment space）而非速度空間進行碰撞。通過線性變換，將分布函數 $f_i$ 映射到一組具有明確物理意義的矩（如密度、動量、應力張量、熱流等）。然後對這些不同的矩使用不同的鬆弛速率。這種方法：

- 增強了數值穩定性：通過快速衰減高階非物理矩，可以有效抑制數值振盪。
- 提高了模型的靈活 性：可以獨立調節剪切粘度和體積粘度，並模擬普朗特數不為1的流體。

### 不可壓縮模型

標準LBM是弱可壓縮的，即密度會隨著壓力發生微小波動。對於嚴格不可壓縮的流動（如液體），這會引入誤差。為此，研究者開發了專門的不可壓縮LBM模型。這類模型通過引入壓力分布函數代替密度分布函數，或通過投影法來保證速度場在每一步都滿足無散條件（$\nabla \cdot \mathbf{u} = 0$）。

## 邊界條件處理

在LBM中，邊界條件的處理至關重要，它直接影響模擬的精度和物理真實性。

### 反彈格式

反彈格式（Bounce-back scheme）是處理無滑移固體邊界最簡單和最常用的方法。當一個粒子分布函數遷移到固體節點時，它會被「反彈」回流體域，即其速度方向反轉。

- 標準反彈：實現簡單，但精度僅為一階，且有效邊界位置會隨鬆弛時間變化。
- 半步長反彈：通過假設物理邊界位於流體節點與固體節點之間，可以將精度提升至二階。

由於其簡便性，反彈格式在處理具有複雜幾何形狀的系統（如多孔介質）時特別有效。

### 開放邊界條件

對於計算域的入口和出口，需要指定宏觀物理量（如速度或壓力）。Zou-He邊界條件是一種廣泛使用的方法，它基於非平衡態反彈的假設，結合質量和動量守恆，可以解析地求解出邊界節點上未知的分布函數，從而精確施加宏觀邊界條件。其他更複雜的動力學邊界條件也能實現更高的精度和普適性。

## 應用領域

LBM的介觀特性使其在許多傳統CFD方法面臨挑戰的領域中表現出色。

### 多相與多組分流

多相流是LBM最成功的應用領域之一。與需要顯式追蹤或重構相界面的傳統方法（如流體體積法或水平集方法）不同，LBM中的相界面是作為粒子間相互作用的自然結果而自動產生的。

- 偽勢模型（Shan-Chen模型）：通過引入基於密度的非局部相互作用力，該模型可以自發地實現相分離，形成液滴、氣泡等。它還能自然地處理壁面的潤濕現象。
- 自由能模型：基於嚴格的熱力學自由能泛函，該模型在描述相變動力學和複雜潤濕現象時具有更高的物理保真度。
- 顏色梯度模型：一種早期的多相流模型，通過引入「顏色」場來標記不同流體，並通過擾動碰撞算子產生表面張力。

### 熱流體動力學

為了模擬包含熱傳導和熱對流的流動，標準LBM需要擴展以求解能量方程。常見的方法是雙分布函數模型（DDF），即使用一套分布函數模擬流場，另一套獨立的分布函數模擬溫度場或內能場。這種方法可以靈活地調節普朗特數，並能包含粘性耗散和壓縮功等複雜熱物理效應。

### 多孔介質流動

LBM是模擬多孔介質流動的首選工具之一。得益於其便捷的邊界處理能力，可以直接將CT掃描得到的複雜孔隙結構圖像映射為LBM的計算網格，實現所謂的「數字岩心」分析。LBM被廣泛用於計算岩石的滲透率、研究多相驅替過程以及模擬燃料電池氣體擴散層中的物質傳輸。

### 其他應用

LBM的應用還包括：

- 微流控：模擬微通道中的滑移流、電滲流和液滴操控。
- 顆粒懸浮流：研究流體中固體顆粒的運動、沉降和流變學特性。
- 反應流與燃燒：通過與模型降階技術結合，LBM已被用於模擬複雜的燃燒過程。
- 生物流體力學：模擬血液流動、細胞變形以及植物維管系統中的水分輸運。
- 磁流體動力學（MHD）與淺水方程。

## 優缺點

### 優點

- 演算法簡單，易於實現：核心的「遷移-碰撞」步驟直觀且易於編程。
- 高度並行性：演算法的局部性使其非常適合在大規模並行計算平台（如GPU和超算集群）上高效運行。
- 複雜邊界處理：通過簡單的反彈格式，可以輕鬆處理任意複雜的幾何形狀，無需生成複雜的貼體網格。
- 多相流模擬：能夠自然地捕捉界面動力學，無需額外的人工界面追蹤演算法。
- 堅實的物理基礎：直接源於動理學理論，使其在模擬微尺度效應時具有天然優勢。

### 缺點

- 弱可壓縮性限制：標準模型是為低馬赫數（通常 $Ma < 0.3$）流動設計的，不直接適用於高速可壓縮流（如激波）。
- 均勻網格要求：標準LBM依賴於均勻的笛卡爾網格，這在需要局部網格加密以解析邊界層等問題時效率較低。（儘管非均勻網格和自適應網格的LBM已被開發）
- 記憶體消耗：需要存儲每個節點的分布函數，相比於只存儲宏觀量的傳統方法，記憶體需求通常更大。
- 熱力學模型複雜性：要精確模擬包含真實氣體狀態方程和複雜熱物理過程的系統，需要對模型進行大量擴展和修正。