题 1

一无限流体被单位质量力 $\mu r^{-3/2}$ 的力作用，方向指向原点。若开始时流体静止且其中有一个空腔 $r=c$ ，证明：在过 $(\frac{2}{5 \mu})^{1/2} c^{5/4}$ 时间之后，空腔被流体填满。

这里假设流体不可压缩。取球坐标系 $(r, \theta, \phi)$ 进行计算，并以空腔中心为原点。

思路1：功能关系

假设当 $t=t_0$ 时，空腔边界的球壳形薄液面位于 $r=r_0$ 处。该薄液面厚度为 $\mathrm{d}r_0$ , 液面速度 $v_r (r_0, t_0) = 0$ 。

若 $t_0$ 时刻，位于 $r=r_0$ 处且厚度为 $\mathrm{d}r_0$ 的液面向内推动了 $\mathrm{d} r$ ，由于该液面的质量积为 $\mathrm{d}m = \rho (4 \pi r_0^2) \mathrm{d}r_0$ ，则这一过程中体积力（ $F_{br} = -\mu r^{-3/2}$ ）所作的功为

\mathrm{d} W = \mathrm{d}m \cdot \int_{r_0}^{+\infty} F_{br} \mathrm{d}r = \mathrm{d}m \cdot (-2 \mu r_0^{-1/2}) = -8 \mu \rho \pi r_0^{3/2} \mathrm{d}r_0 .

当 $t_1 = t_0 + \delta_t$ 时，该液面前进到 $r=r_1$ 位置（ $r_1 \le r_0$ ），此时液面速度为 $v_r(r_1, t_1) = \dfrac{\mathrm{d} r_1}{\mathrm{d} t}$ 。

所以薄液面从 $r_0$ 运动到 $r_1$ 过程中， $F_{br}$ 做的总功为：

W = \int_{r_0}^{r_1} \mathrm{d} W = \frac{16}{5} \mu \rho \pi (r_0^{5/2} - r_1^{5/2})

记速度 $\bold{v} = [v_r, 0, 0]^T$ 。代入球坐标系下的连续性方程，化简得：

\frac{\partial (r^2 v_r)}{\partial r} = 0

对 $r$ 积分，可得 $r^2 v_r = f(t)$ ，其中 $f(t)$ 为待定函数。根据这一方程，可得在 $t_1$ 时刻满足 $r_0^2 v_r(r_0, t_1) = r_1^2 v_r(r_1, t_1)$ 。

在 $t_1$ 时刻， $r=r_0$ 层液体的动能表示为

\mathrm{d} E = \frac{1}{2} \cdot \mathrm{d}m \cdot v_r^2(r_0, t_1) = 2 \rho \pi \cdot \frac{r_1^4}{r_0^2} \cdot v_r^2(r_1, t_1) \mathrm{d}r_0

这一时刻下整个系统的动能为：

E = \int_{r_1}^{+\infty} \mathrm{d} E = 2 \rho \pi r_1^3 \cdot v_r^2(r_1, t_1)

据此，我们令 $t_0$ 时刻即为初始时刻（ $t_0=0$ ），满足 $r_0 = c$ 和 $v_r(t=0)=0$ 。所以 $t_0$ 时刻的动能为 0。对 $t_0$ 到 $t_1$ 时刻的过程， $F_{br}$ 所做的功 $W$ 全部转化为流场动能 $E$ ，即： $W=E$ 。

\frac{16}{5} \mu \rho \pi (c^{5/2} - r_1^{5/2}) = 2 \rho \pi r_1^3 \cdot v_r^2(r_1, t_1) - 0

也就有：

v_r(r_1, t_1) = \frac{\mathrm{d} r_1}{\mathrm{d} t} = \sqrt{\frac{8 \mu}{5}} \cdot \sqrt{r_1^{-3/2} (c^{5/2} - r_1^{5/2})}

分离变量后解得：

\frac{16}{25} (c^{5/2} - r_1^{5/2}) = \left( \sqrt{\frac{8 \mu}{5}} t + c_1 \right)^2

其中 $c_1$ 为待定常数。当 $t=0$ 时， $r_1 = c$ ，所以 $c_1 = 0$ 。

假设当 $t=T_{m}$ 时 $r_1=0$ ，则可解得： $T_m = (\frac{2}{5 \mu})^{1/2} c^{5/4}$ 。

思路2：联立动量方程和连续性方程

记速度 $\bold{v} = [v_r, 0, 0]^T$ 。代入连续性方程，化简得：

\frac{\partial (r^2 v_r)}{\partial r} = 0

对 $r$ 积分，可得 $r^2 v_r = f(t)$ ，其中 $f(t)$ 为待定函数。

该场景下的运动方程为

\rho \frac{\mathrm{D} \bold{v}}{\mathrm{D} t} = \rho \bold{F}_b

所以，动量方程写作：

\frac{\partial v_r}{\partial t} + v_r \frac{\partial v_r}{\partial r} = F_{br}

由题目可知 $F_{br} = - \mu r^{-3/2}$ 。将 $r^2 v_r = f(t)$ 代入到运动方程中，得：

\frac{f'(t)}{r^2} + v_r \frac{\partial v_r}{\partial r} = -\mu r^{-3/2}

记 $R(t)$ 为真空球面的半径变化，则 $v_r (r=R(t)) = R'(t)$ 。对 $r$ 积分，得到：

\begin{aligned} \int_{\infty}^{R(t)} \frac{f'(t)}{r^2} + v_r \frac{\partial v_r}{\partial r} \mathrm{d}r &= -\mu \int_{\infty}^{R(t)} r^{-3/2} \mathrm{d}r \\ \left.\left[ -\frac{f'(t)}{r} + \frac{(v_r)^2}{2} \right]\right|_{\infty}^{R(t)} &= \left. \frac{2 \mu}{\sqrt{r}} \right|_{\infty}^{R(t)} \end{aligned}

代入 $v_r (r=\infty)=0$ ，整理得到：

-\frac{f'(t)}{R(t)} + \frac{(R'(t))^2}{2} = \frac{2 \mu}{\sqrt{R(t)}}

在边界上有 $f(t) = R^2(t) \cdot v_r (r=R(t)) = R^2 \cdot R'$ ，所以 $f'(t) = 2 R R'^{2} + R^2 \cdot R{''}$ 。上式化简为：

-\frac{3}{2} R'^{2} -R \cdot R{''} = \frac{2 \mu}{\sqrt{R}}

由于 $\frac{\mathrm{d} (R^{1/2} (R')^{2})}{\mathrm{d} R} =\frac{R^{-1/2}}{2} (R')^{2} + 2 R^{1/2} R{''}$ ，上式整理得：

-\frac{\mathrm{d}(R^{1/2} (R')^{2})}{8\mu + 5 (R^{1/2} (R')^{2})} = \frac{\mathrm{d}R}{2R}

分离变量后并整理出 $R'$ ，可表示为：

R' = \sqrt{\frac{R^{-1/2}}{5} (c_1 R^{-5/2} - 8 \mu) }

将 $R(t=0)=c$ 和 $R'(t=0) = 0$ 代入，解得： $c_1 = 8 \mu c^{5/2}$ 。即：

R' = \sqrt{\frac{R^{-1/2}}{5} 8 \mu (c^{5/2} R^{-5/2} - 1)} \\= \sqrt{\frac{8 \mu}{5}} \cdot \sqrt{R^{-1/2} (c^{5/2} R^{-5/2} - 1)}

综上，真空球面的总运动时长表示为：

T = \int_{0}^{c} \frac{1}{R'} \mathrm{d}R = \sqrt{\frac{5}{8 \mu}} \int_{0}^{c} \left[ R^{-1/2} \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R

因为：

\begin{aligned} & \int_{0}^{c} \left[ R^{-1/2} \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R \\ =& \int_{0}^{c} R^{1/4} \left[ \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}R \\ =& \frac{4}{5} \int_{0}^{c} \left[ \left( \left(\frac{c}{R}\right)^{5/2} - 1 \right) \right]^{-1/2} \mathrm{d}(R^{5/4}) \\ \xlongequal{t=R^{5/4}} & \frac{4}{5} \int_{0}^{c^{5/4}} \left[ \left( c^{5/2} t^{-2} - 1 \right) \right]^{-1/2} \mathrm{d}t \quad= \frac{4}{5} \int_{0}^{c^{5/4}} t \left[ \left( c^{5/2} - t^2 \right) \right]^{-1/2} \mathrm{d}t \\ =& \frac{4}{5} \cdot \left(-\frac{1}{2}\right) \int_{0}^{c^{5/4}} \left[ \left( c^{5/2} - t^2 \right) \right]^{-1/2} \mathrm{d}(c^{5/2} - t^2) \\ =& \frac{4}{5} \cdot \left(-\frac{1}{2}\right) \left.2 (c^{5/2} - t^2)^{1/2}\right|_{0}^{c^{5/4}} = \frac{4}{5} c^{5/4} \end{aligned}

所以 $T = \sqrt{\frac{5}{8 \mu}} \cdot \frac{4}{5} c^{5/4} = \sqrt{\frac{2}{5 \mu}} c^{5/4}$ .

证毕。

题 2

体积为 $\frac{4}{3} \pi a^3$ 的液体充满于两个同心球面之间，外球面有压强 $\pi$ 作用，无质量力，内球面压强为 0 。开始时液体静止，内球面半径为 $2a$ 。证明：当内球面半径变为 $a$ 时，其速度为 $\displaystyle \left(\frac{14 \pi}{3 \rho} \frac{2^{1/3}}{2^{1/3} - 1}\right)^{\frac{1}{2}}$ 。

解

这里假设流体不可压缩。取球坐标系 $(r, \theta, \phi)$ 进行计算，并以空腔中心为原点。

记初始状态下（即 $t=0$ 时刻），内球面半径为 $R_0=2 a$ ，外球面半径为 $R_1$ ；当内球面半径变为 $a$ 时（即 $t_1$ 时刻），此时的内、外球面半径为分别为 $R_{01} = a$ 和 $R_{11}$ 。由于流体体积 $V = \frac{4}{3} \pi a^3$ ，可得：

\begin{aligned} V &= \frac{4}{3} \pi (R_1^3 - R_0^3) \\ &= \frac{4}{3} \pi (R_{11}^3 - R_{01}^3) \end{aligned}

解得： $R_1 = 9^{1/3} a$ 和 $R_{11} = 2^{1/3} a$ 。

记速度 $\bold{u} = [v_r, 0, 0]^T$ 。代入连续性方程，化简得：

\frac{\partial (r^2 v_r)}{\partial r} = 0

对 $r$ 积分，可得 $r^2 v_r = f(t)$ ，其中 $f(t)$ 为待定函数。

对 $t \in [0,t_1]$ 时间段，外界压强在外表面做的功为：

W = \int_{R_{1}}^{R_{11}} (-\pi) S(r) \mathrm{d} r = \frac{28}{3} \pi^2 a^3

其中 $S(r) = 4 \pi r^2$ 为球的表面积公式。

在 $t_1$ 时刻， $r^2 v_r = f(t_1)$ 为定值。若记此时刻的内球面速度为 $\bold{u}_0 = [v_0, 0, 0]^T$ ， $r$ 处的流体速度为 $\bold{u}_r = [v_r, 0, 0]^T$ （ $r \in [R_{01}, R_{11}]$ ），则 $R_{01}^2 v_0 = r^2 v_r$ 。所以 $v_r = \left( \frac{R_{01}}{r} \right)^2 v_0$ 。

$t_1$ 时刻的流场总动能表示为：

E_{t_1} = \int_{R_{01}}^{R_{11}} \frac{1}{2} v_r^2 \cdot (\rho S(r) \mathrm{d}r) = 2 \pi \rho v_0^2 a^3 \frac{2^{1/3} - 1}{2^{1/3}}

由于初始时刻的动能 $E_0 = 0$ ，根据功能关系可知：

W = E_{t_1} - E_0

解得 $t_1$ 时刻的内球面速度为： $\displaystyle v_0= - \left(\frac{14 \pi}{3 \rho} \frac{2^{1/3}}{2^{1/3} - 1}\right)^{\frac{1}{2}}$ ，方向指向原点。

LBM笔记：Exact Difference Method 构建的源项

发表于 2024-09-07 更新于 2024-10-25 分类于格子Boltzmann方法，源项

本文同步发表于知乎和CSDN

我在阅读其他文献时，看到有文献提到 Kupershtokh 的一篇关于 LBM 外力项的文章^[1]，他提出的格式在其他文章中又被称为 Exact difference method（EDM）。

本想搜索这篇文章的PDF进行阅读，但最后只找到一个网站介绍了部分信息。这份粗略笔记的内容，都是基于其他文献^[2]^[3]^[4]中对该方法的介绍而整理的。

介绍

Boltzmann-BGK方程的近似

密度分布函数 $f$ 的 Boltzmann-BGK方程写作：

\frac{\partial f}{\partial t} + \bold{\xi} \cdot \nabla f + \bold{a} \cdot \nabla_{\bold{\xi}} f = -\frac{1}{\tau_f} \left[ f - f^{(eq)}\right] \tag{1}

其中：

f^{(eq)} = \frac{\rho}{(2 \pi R_g T)^{D/2}} \exp{\left(- \frac{(\bold{\xi} - \bold{u})^2}{2 R_g T}\right)} \tag{2}

符号含义分别为

$\bold{\xi}$ : 粒子速度；
$\rho$ : 流体密度；
$\bold{u}$ : 流场宏观速度；
$R_g$ : 气体常数；
$T$ : 温度；
$\bold{a}=\mathrm{d}\bold{u}/\mathrm{d}t$ : 流场加速度。

当Knudsen数较小时，Boltzmann方程的体力项可以基于 $f^{(eq)}$ 进行近似，即

\nabla_{\bold{\xi}} f \approx \nabla_{\bold{\xi}} f^{(eq)}

根据式(2)， $f^{(eq)}$ 是分子随机速度 $\bold{C} = \bold{\xi} - \bold{u}$ 的指数函数，则：

\frac{\partial f^{(eq)}}{\partial \bold{u}} = -\frac{\partial f^{(eq)}}{\partial \bold{\xi}} = -\nabla_{\bold{\xi}} f^{(eq)}

若假设密度 $\rho$ 为常数，则 $\frac{\mathrm{d} \rho}{\mathrm{d} t}=0$ ，那么 $f^{(eq)}$ 的全导数为：

\begin{aligned} \frac{\mathrm{d} f^{(eq)}}{\mathrm{d} t} &= \frac{\partial f^{(eq)}}{\partial \bold{u}} \cdot \frac{\mathrm{d} \bold{u}}{\mathrm{d} t} + \frac{\partial f^{(eq)}}{\partial \rho} \cdot \frac{\mathrm{d} \rho}{\mathrm{d} t} \\ &\approx \frac{\partial f^{(eq)}}{\partial \bold{u}} \cdot \frac{\mathrm{d} \bold{u}}{\mathrm{d} t} = -\bold{a} \cdot \nabla_{\bold{\xi}} f^{(eq)} \end{aligned}

也就是说，Kupershtokh 等^[1:1]^[2:1]近似的 Boltzmann-BGK方程写作：

\frac{\partial f}{\partial t} + \bold{\xi} \cdot \nabla f = -\frac{1}{\tau_f} \left[ f - f^{(eq)}\right] + \frac{\mathrm{d} f^{(eq)}}{\mathrm{d} t} \tag{3}

近似方程的离散化

将式(3)在 $[t,t+\delta_t]$ 内沿着 $\bold{e}_\alpha$ 方向积分，可得

\begin{aligned} f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t) - f_\alpha (\bold{x}, t) =& -\frac{1}{\tau} \left[ f_\alpha (\bold{x}, t) - f_\alpha^{(eq)} (\rho, \bold{u^*}) \right] \\ & + \left[ f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t)) \right] \end{aligned} \tag{4}

其中：

$\delta_t$ 为时间步长；
$\tau = \tau_f / \delta_t$ 为无量纲松弛时间；
$\bold{x}$ 为网格点位置， $\rho, \bold{u^*}$ 分别为该点在 $t$ 时刻的密度和速度（与宏观速度 $\bold{u}$ 有区别）；
$f_\alpha, f_\alpha^{(eq)}$ 分别为 $\bold{e}_\alpha$ 方向的分布函数和平衡态；

\begin{aligned} f_\alpha^{(eq)} &= w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{(\bold{e}_\alpha \cdot \bold{u})^2}{2 c_s^4} - \frac{|\bold{u}|^2}{2 c_s^2} \right] \\ &= w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{1}{2 c_s^4}(\bold{u}\bold{u}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \\ \end{aligned} \tag{4-1}

这里 $w_\alpha$ 为 $\bold{e}_\alpha$ 的权重系数， $[\bold{I}]$ 为单位矩阵；
两个同样为m*n大小的矩阵 $[\bold{A}]$ 和 $[\bold{B}]$ ，则 $\displaystyle [\bold{A}] : [\bold{B}] = \sum^{m}_{i=1}\sum^{n}_{j=1} A_{ij} B_{ij}$ ；

$F_\alpha$ 为外力项，其表达式基于外力 $\bold{F} = \rho \dfrac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} = \rho \bold{a}$ 进行构建。

F_\alpha = f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t))

并且

\begin{aligned} \bold{u^*}(t) &= \left(\sum_{\alpha} f_\alpha \bold{e}_\alpha \right) / \left(\sum_{\alpha} f_\alpha \right) \\ \bold{u^*}(t+\delta_t) &= \bold{u^*}(t) + \int_{t}^{t+\delta_t} \frac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} \mathrm{d} t \\ &= \bold{u^*}(t) + \int_{t}^{t+\delta_t} \frac{\bold{F}}{\rho} \mathrm{d} t \approx \bold{u^*}(t) + \frac{\bold{F}}{\rho} \delta_t \end{aligned} \tag{5}

这意味着加速度 $\bold{a}=\bold{F}/\rho$ 在 $[t,t+\delta_t]$ 时间段内被视为常数。

Kupershtokh等在这篇文章^[2:2]说明，在使用该源项时，流体宏观量的计算应表示为：

\rho = \sum_{\alpha} f_\alpha ,\quad \rho \bold{u} = \sum_{i} \bold{e}_\alpha f_\alpha + \frac{\delta_t \bold{F}}{2} \tag{6}

Li 等^[4:1]通过 Chapman-Enskog 展开证明：该格式向宏观方程引入的误差项为：

\mathbf{Error}_{\mathrm{EDM}} = -\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right) \tag{7}

EDM 源项的展开式

为便于分析，这里采用如下书写形式的平衡态：

f_\alpha^{(eq)} = w_\alpha \rho \cdot \left[ 1 + \frac{\bold{e}_\alpha \cdot \bold{u}}{c_s^2} + \frac{1}{2 c_s^4}(\bold{u}\bold{u}):(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right]

将其带入 EDM 源项，记 $\Delta\bold{u} = \frac{\bold{F}}{\rho} \delta_t$ ，整理得：

\begin{aligned} F_{\alpha} &= f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t)) \\&= \rho w_\alpha \left\{ \frac{\bold{e}_\alpha \cdot \Delta\bold{u}}{c_s^2} + \newline \frac{1}{2 c_s^4} \left[ (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) \right.\right. \\ &\quad \left. - \bold{u}^* \bold{u}^* \right] : \left.(\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right\} \end{aligned}

这里涉及矩阵 $(\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u})$ 和 $\bold{u}^* \bold{u}^*$ 的相减，由于：

(\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) = (u^*_{i} + \Delta u_i)(u^*_{j} + \Delta u_j), \quad \bold{u}^* \bold{u}^* = u^*_{i} u^*_{j}

其中 $u_i^*$ 和 $\Delta u_i = F_i \delta_{t} / \rho$ 分别为 $\bold{u}^*$ 和 $\Delta\bold{u}$ 在 i 轴上的分量，所以：

\begin{aligned} (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) - \bold{u}^* \bold{u}^{*} &= (u^*_{i} + \Delta u_i)(u^*_{j} + \Delta u_j) - u^*_{i} u^*_{j} \\ &= \frac{\delta_t}{\rho} (u_i^* F_j + u_j^* F_i) + \frac{\delta_t^2}{\rho^2} F_i F_j \\ &= \frac{\delta_t}{\rho} (u_i^* + \frac{F_i \delta_t}{2 \rho}) F_j + \frac{\delta_t}{\rho} (u_j^* + \frac{F_j \delta_t}{2 \rho}) F_i \end{aligned}

这里令 $\bold{v}_{\rm{EDM}} = \bold{u}^{*} + \frac{\bold{F}}{2 \rho} \delta_{t}$ ,则 EDM 源项的展开式化简为：

F_\alpha = \delta_t w_\alpha \left[ \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} + \frac{1}{2 c_s^4} \left( \bold{v}_{\rm{EDM}} \bold{F} + \bold{F} \bold{v}_{\rm{EDM}} \right) : (\bold{e}_\alpha \bold{e}_\alpha - c_s^2 [\bold{I}]) \right] \tag{8}

上式即为 EDM 源项的展开结果。

Kupershtokh, A. L. New method of incorporating a body force term into the lattice Boltzmann equation. in Proceeding of the 5th international EHD workshop 241–246 (Poitiers, France, 2004). ↩︎ ↩︎
A.L. Kupershtokh, D.A. Medvedev, & D.I. Karpov (2009). On equations of state in a lattice Boltzmann method. Computers & Mathematics with Applications, 58(5), 965-974. DOI:10.1016/j.camwa.2009.02.024. ↩︎ ↩︎ ↩︎
Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, & Q. Liu (2016). Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Progress in Energy and Combustion Science, 52, 62-105. DOI:10.1016/j.pecs.2015.10.001. ↩︎
Li, Q., Zhou, P., & Yan, H. (2016). Revised Chapman-Enskog analysis for a class of forcing schemes in the lattice Boltzmann method. Physical Review E, 94, 043313. DOI:10.1103/PhysRevE.94.043313. ↩︎ ↩︎

【代码】DEM 的颗粒流数据粗粒化

发表于 2024-08-10 更新于 2024-10-25 分类于离散单元法

基于这篇的内容，我用 Python 写了一份对单种颗粒系统进行粗粒化的代码。

其具体使用流程为：

obj = CG_mono_3D(...): 初始化。
obj.set_domain(...): 设置矩形的计算域。
obj.get_CG_data(R, CG_width): 对 R 点进行粗粒化，其粗粒化的长度为 CG_width。

DEM_CG_mono.py

"""
Function to convert DEM simulation data into coarse graining field.

NOTE:
    This function is used for graunlar system which consists of single type of particles.
"""
import typing as _T

import numpy as np
import numba as nb
import numpy.typing as npt

from tqdm import tqdm
from scipy.special import erf
from scipy.integrate import simpson


@nb.jit(forceobj=True)
def CG_Gaussian(r: npt.NDArray[np.double], w: float) -> float:
    """Gaussian Coarse Graining function.

    Parameters
    ====
    r: numpy.NDArray[double]
        Relative position vector, whose shape is `(3,)`.
    w: float
        Course grain width `w`.
                                        
    Return
    ===
    W: float
        Average weight.
    """
    c = 3 * w    # Coarse graining cut-off width
    c_2 = c**2

    # Norm of vector `r`
    r_norm: float = r[0]**2 + r[1]**2 + r[2]**2
    # define weight
    W = 0.0
    V_w = 0.0
    if r_norm <= c_2:
        V_w = np.sqrt(8) * np.pi**1.5 * w**3 * erf(c/w/np.sqrt(2)) - \
            4 * c * w**2 * np.pi * np.exp(-c**2 / w**2 / 2)
        W = np.exp(-r_norm / (2 * w**2)) / V_w
    return W


@nb.jit("f8[:, :](f8, f8, f8[:], f8[:])",
        nogil=True, nopython=True, parallel=True)
def SIGMA_q_k(mass_i: float,
              W_i: float,
              vel_i: npt.NDArray[np.double],
              U: npt.NDArray[np.double]) -> npt.NDArray[np.double]:
    """Calculate part of `sigma_qk` contributed by ball `i`

    Parameters
    ---
    mass_i: float
        Mass of ball i.
    W_i: float
        Weight coefficient of ball i.
    vel_i: numpy.NDArray[double]
        Velocity vector of ball i, whose shape is (3,).
    U: numpy.NDArray[double]
        Coarse grained velocity vector, whose shape is (3,).

    Return
    ---
    res: numpy.NDArray[double]
        3x3 Matrix of `sigma_qk`caused by ball i.
    """
    res = np.zeros((3,3), dtype=np.double)
    for i in nb.prange(3):
        for j in nb.prange(3):
            res[i, j] += mass_i * W_i * \
                (vel_i[i] - U[i]) * (vel_i[j] - U[j])
    return res


@nb.jit("f8[:, :](f8[:], f8[:], f8)",
        nogil=True, nopython=True, parallel=True)
def SIGMA_q_c(F_ij: npt.NDArray[np.double],
              l_ij: npt.NDArray[np.double],
              W_ij: float) -> npt.NDArray[np.double]:
    """Calculate part of `sigma_qk` contributed by contact from i to j.

    Parameters
    ---
    F_ij: numpy.NDArray[double]
        Contact force vector, pointing from i to j.
    l_ij: numpy.NDArray[double]
        Contact branch vector, pointing from i to j.
    W_ij: float
        Integrated weight coefficient, from i to j.
    
    Return
    ---
    res: numpy.NDArray[double]
        3x3 Matrix of `sigma_qc`caused by contact from i to j.
    """
    res = np.zeros((3,3), dtype=np.double)
    for i in nb.prange(3):
        for j in nb.prange(3):
            res[i,j] = F_ij[i] * l_ij[j] * W_ij
    return res


@nb.jit("f8[:](f8[:], f8[:], f8[:,:], b1, b1[:])",
        nopython=True)
def get_distance(R: npt.NDArray[np.double],
                 r_i: npt.NDArray[np.double],
                 lims: npt.NDArray[np.double],
                 hasPeriodic: bool,
                 isPeriodic: npt.NDArray[np.bool_]):
    dR = r_i - R
    if hasPeriodic:
        for q in range(3):
            if not isPeriodic[q]:
                continue
            # Deal with periodic boundary
            _r_2L = abs(r_i[q] - lims[q,0])
            _r_2R = abs(r_i[q] - lims[q,1])
            _R_2L = abs(R[q] - lims[q,0])
            _R_2R = abs(R[q] - lims[q,1])
            r_near_left = _r_2L < _r_2R
            R_near_left = _R_2L < _R_2R
            r2R = abs(dR[q])
            if R_near_left and (not r_near_left):
                tmp = _r_2R + _R_2L
                if tmp < r2R: dR[q] = -tmp
            elif (not R_near_left) and r_near_left:
                tmp = _r_2L + _R_2R
                if tmp < r2R: dR[q] = tmp
    return dR   


class CG_mono_3D:
    def __init__(self,
                 rho: float, 
                 pos: npt.NDArray[np.double],
                 vel: npt.NDArray[np.double],
                 radius: npt.NDArray[np.double],
                 ball_id: npt.NDArray[np.int_],
                 contact_force: npt.NDArray[np.double],
                 contact_branch: npt.NDArray[np.double],
                 contact_end_ball: npt.NDArray[np.double],
                 in_flow: _T.Optional[npt.NDArray[np.bool_]] = None,
                 dp: _T.Optional[float] = None, *,
                 progress_bar: bool = False
                ) -> None:
        """
        Parameters
        ----
        rho: float
            Density of all particles.
        pos: numpy.NDArray[double]
            Position vector, whose shape is `[particle_num, 3]`.
        vel: numpy.NDArray[double]
            Velocity vector, whose shape is `[particle_num, 3]`.
        radius: numpy.NDArray[double]
            Radius list, whose shape is `(particle_num,)`.
        ball_id: numpy.NDArray[int]
            ID number of all balls.
        contact_force: numpy.NDArray[double]
            Contact force list, whose shape is `[contact_num, 3]`.
        contact_branch: numpy.NDArray[double]
            List of contact branch, whose shape is `[contact_num, 3]`.
        contact_end_ball: numpy.NDArray[double]
            List of contacted balls tuple, whose shape is `[contact_num, 2]`.
        in_flow: numpy.NDArray[np.bool_], optional
            Mask array of ball-in-flow, whose shape is `(particle_num,)`.
        dp: float, optional
            Mean diameter of all particles.
        progress_bar: bool, optional, keyword argument only.
            Switch to control whether show the process. Default is False.

        Note
        ----
        1. Take the `i`-th ball for example, its ball ID is `ball_id[i]`.
           And its radius and velocity are `radius[i]` and `vel[i,:]`, 
           respectively. Note that `ball_id[i]` is not equal to  `i`.
        2. For the `j`-th contact, The ID of two contacted balls are 
           `b1 = contact_end_ball[j, 0]`, and `b2 = contact_end_ball[j ,1]`.
        """
        self.particle_num = pos.shape[0]
        self.contact_num  = contact_force.shape[0]

        self._id    = ball_id  # Ball id
        self.pos    = pos      # Position vector
        self.vel    = vel      # Velocity vector
        self.radius = radius   # Ball radius

        # Ball density
        # |=> NOTE: Since this function is used for mono-disperse system,
        #    The density of particles are constant.
        self.rho = rho

        # => Contact properties
        self.c_force  = contact_force     # Contact force
        self.c_branch = contact_branch    # Contact branch
        self.c_ends   = contact_end_ball  # Contact end ball tuple

        # => Optional parameters
        # Ball-in-flow mask array.
        self.inflow = in_flow
        # Mean particle diameter
        if dp is None:
            dp = np.mean(radius) * 2
        else:
            self.dp = dp

        # => Mass array
        self.mass = 4/3 * np.pi * np.power(radius, 3) * rho

        # => Ball-ID position array.
        #    For example, `self._c_end_id[j,0]` corresponds to the first end
        #    ball of `j`-th contact (whose ID is `self.c_ends[j,0]`).
        #    And we have `self.c_ends[j,0] = self._id[self._c_end_id[j,0]]`.
        self._c_end_id = [
            [np.where(ball_id == contact_end_ball[i,0])[0][0],
             np.where(ball_id == contact_end_ball[i,1])[0][0]]
            for i in range(self.contact_num)
        ]
        self._c_end_id = np.array(self._c_end_id, dtype=np.intc)

        # => STATUS of domain config
        self._domain_init = False

        # => Max search width
        self._MAX_width = 7

        # => Option to show progress bar
        self._use_tqdm = progress_bar
        pass

    def set_domain(self,
                   xlim: _T.Tuple[float, float],
                   ylim: _T.Tuple[float, float],
                   zlim: _T.Tuple[float, float],
                   periodic_flag: _T.Tuple[bool, bool, bool]):
        """Function to set cuboid domain of coarse graining.
        
        Parameters
        ----
        xlim: (float, float)
            Range of x-coordinate, contains with `(xMin, xMax)`.
        ylim: (float, float)
            Range of y-coordinate, contains with `(yMin, yMax)`.
        zlim: (float, float)
            Range of z-coordinate, contains with `(zMin, zMax)`.
        periodic_flag: (bool, bool, bool)
            Flag tuple to determinate which axis is periodic.
        """
        xlim = (min(*xlim), max(*xlim))
        ylim = (min(*ylim), max(*ylim))
        zlim = (min(*zlim), max(*zlim))
        self._lims = np.array([xlim, ylim, zlim])

        self._isPeriodic  = np.array(periodic_flag, dtype=np.bool_)
        self._hasPeriodic = np.any(periodic_flag)
    
        # Update state
        self._domain_init = True

    def get_CG_data(self,
                    R: npt.NDArray[np.double],
                    w: float):
        """
        Get coarse grainning data at center node `R`.

        Parameter
        ---
        R: numpy.NDArray[double]
            Center of coarse graining node.
        w: float
            Width of coarse graining.
        
        Return
        ---
        rho: float
            Coarse grained density at node `R`.
        U: numpy.NDArray[double]
            Coarse grained velocity at node `R`.
        sigma: numpy.NDArray[double]
            Coarse grained stress tensor at node `R`.
        """
        assert R.shape == (3,)
        if self._domain_init is False:
            raise Exception("Exception: `self.set_domain` haven't run to set domain.")

        rho, U   = self.get_CG_vel(R, w)
        sigma_qk = self.get_CG_sigma_qk(R, w, U)
        sigma_qc = self.get_CG_sigma_qc(R, w)
        sigma    = sigma_qk + sigma_qc
        return rho, U, sigma

    # ============================ #
    # === Sub-process function === #
    # ============================ #

    def get_CG_vel(self,
                   R: npt.NDArray[np.double],
                   CG_width: float) -> _T.Tuple[float, npt.NDArray[np.double]]:
        """Get velocity and density.

        Parameter
        ---
        R: numpy.NDArray[double]
            Center of coarse graining node.
        CG_width: float
            Width of coarse graining.

        Return
        ---
        rho: float
            Coarse grained density at node `R`.
        U: numpy.NDArray[double]
            Coarse grained velocity at node `R`.
        """
        MAX_width = self._MAX_width * CG_width

        # Variable to store result
        rho, Jx, Jy, Jz = 0.0, 0.0, 0.0, 0.0

        # Start loop
        _ball_iter = range(self.particle_num)
        if self._use_tqdm:
            print("|> [Get coarse graining velocity]")
            _ball_iter = tqdm(_ball_iter, unit="Ball")

        for i in _ball_iter:
            # Skip nodes not in flow group
            if (self.inflow is not None) and (not self.inflow[i]): continue

            # Get distance vector
            dR = self._get_distance_vec(R, self.pos[i,:])

            # Skip node which too far from center node `R`.
            if np.linalg.norm(dR) >= MAX_width: continue

            # Get weight
            W_i = CG_Gaussian(-dR, CG_width)

            # Assign value
            rho += self.mass[i] * W_i
            Jx  += self.mass[i] * W_i * self.vel[i, 0]
            Jy  += self.mass[i] * W_i * self.vel[i, 1]
            Jz  += self.mass[i] * W_i * self.vel[i, 2]
            pass

        # Get velocity
        if rho != 0:
            U = np.array([Jx / rho, Jy / rho, Jz / rho], dtype=np.double)
        else:
            U = np.zeros((3,))
        return rho, U

    def get_CG_sigma_qk(self,
                        R: npt.NDArray[np.double],
                        CG_width: float,
                        U: npt.NDArray[np.double]) -> npt.NDArray[np.double]:
        """
        Get stress tensor contributed by kinetic energy.
        
        Parameter
        ---
        R: numpy.NDArray[double]
            Center of coarse graining node.
        CG_width: float
            Width of coarse graining.
        U: numpy.NDArray[double]
            Coarse grained velocity at node `R`.
        
        Return
        ---
        sigma_qk: numpy.NDArray[double]
            Stress tensor contributed by kinetic energy.
        
        Note
        ---
        $\sigma_{\\alpha \\beta}^{q(\mathrm{k})} = \sum_{i \in q} m_i v_{i \\alpha}^{'} v_{i \\beta}^{'} W(\\boldsymbol{r} - \\boldsymbol{r}_i (t))$

        where:
        * $v_{i \\alpha}^{'} = u_{\\alpha}^{q}(\\boldsymbol{r}, t) - u_{i,\\alpha}$ is the $\\alpha$ component of the fluctuation velocity of particle i.
        """
        MAX_width = self._MAX_width * CG_width

        # Result array
        sigma_qk = np.zeros((3,3), dtype=np.double)

        # Start loop
        _ball_iter = range(self.particle_num)
        if self._use_tqdm:
            print("|> [Get stress tensor contributed by kinetic energy]")
            _ball_iter = tqdm(_ball_iter, unit="Ball")

        for i in _ball_iter:
            # Skip nodes not in flow group
            if (self.inflow is not None) and (not self.inflow[i]): continue

            # Get distance vector
            dR = self._get_distance_vec(R, self.pos[i,:])

            # Skip node which too far from center node `R`.
            if np.linalg.norm(dR) >= MAX_width: continue

            W_i = CG_Gaussian(-dR, CG_width)
            sigma_qk += SIGMA_q_k(self.mass[i], W_i, self.vel[i, :], U)

        return sigma_qk

    def get_CG_sigma_qc(self,
                        R: npt.NDArray[np.double],
                        CG_width: float) -> npt.NDArray[np.double]:
        """
        Get stress tensor contributed by contact.
        
        Parameter
        ---
        R: numpy.NDArray[double]
            Center of coarse graining node.
        CG_width: float
            Width of coarse graining.
        
        Return
        ---
        sigma_qc: numpy.NDArray[double]
            Stress tensor contributed by contact.

        Note
        ---
        $\sigma_{\\alpha \\beta}^{q(\mathrm{ck})} = \sum_{i \in q} \sum\limits_{j \in \bar{Q}}^{j \neq i} F_{ij,\\alpha} l_{ij,\\beta} \int_{0}^{1} W(\\boldsymbol{r} - \\boldsymbol{r}_i + s \\boldsymbol{l}_{ij}) \mathrm{d}s$

        where:
        * $v_{i \\alpha}^{'} = u_{\\alpha}^{q}(\\boldsymbol{r}, t) - u_{i,\\alpha}$ is the $\\alpha$ component of the fluctuation velocity of particle i.
        * $\\bar{Q}$ is the union of all particles classes, $q$.
        * $\\boldsymbol{l}_{ij}$ is as the branch vector between two particles, i and j. (from i to j)
        * $\\boldsymbol{F}_{ij}$ is contact force vector between two particles, i and j. (from i to j)
        """
        MAX_width = self._MAX_width * CG_width

        # Result array
        sigma_qc = np.zeros((3,3), dtype=np.double)

        # Start loop
        _c_iter = range(self.contact_num)
        if self._use_tqdm:
            print("|> [Get stress tensor contributed by contact]")
            _c_iter = tqdm(_c_iter, unit="Contact")

        for i in _c_iter:
            # ID of two balls on the i-th contact
            ball_i1 = self._c_end_id[i, 0]
            ball_i2 = self._c_end_id[i, 1]

            # contact force (from i1 to i2)
            F_i12 = self.c_force[i, :]

            # contact branch (from i1 to i2)
            l_i12 = self.c_branch[i, :]

            # Deal with ball_i1
            dR1 = self._get_distance_vec(R, self.pos[ball_i1, :])
            if np.linalg.norm(dR1) < MAX_width:
                W_ij = self._W_ij_4_sigma_qc(dR1, l_i12, CG_width)
                if W_ij != 0:
                    sigma_qc += SIGMA_q_c(F_i12, l_i12, W_ij)

            # Deal with ball_i2
            dR2 = self._get_distance_vec(R, self.pos[ball_i2, :])
            if np.linalg.norm(dR2) < MAX_width:
                W_ij = self._W_ij_4_sigma_qc(dR2, -l_i12, CG_width)
                if W_ij != 0:
                    sigma_qc += SIGMA_q_c(-F_i12, -l_i12, W_ij)

        return sigma_qc

    # ========================== #
    # === Auxiliary function === #
    # ========================== #

    def _W_ij_4_sigma_qc(self, dR, l_ij, CG_width) -> float:
        """
        Parameters
        ---
        dR: numpy.NDArray[double]
            Relative position vector, whose shape is `(3,)`.
        l_ij: numpy.NDArray[double]
            Contact branch vector, whose shape is `(3,)`.
        CG_width: float
            Width of coarse graining.
        
        Return
        ----
        W_ij: float
            Integrated weight coefficient used in the calculation of
            `sigma_qc`.
        """
        _func = lambda s: CG_Gaussian(-dR + s * l_ij, CG_width)
        _x = np.linspace(0, 1, 101)
        _y = np.array([_func(s) for s in _x])
        W_ij = simpson(_y, x=_x)
        return W_ij

    def _get_distance_vec(self,
                          R: npt.NDArray[np.double],
                          r_i: npt.NDArray[np.double]) -> npt.NDArray[np.double]:
        """Get distance between center of coarse graining node `R` and 
        loaction of particle's center `r_i`.

        Parameter
        ---
        R: numpy.NDArray[double]
            Center of coarse graining node.
        r_i: numpy.NDArray[double]
            Loaction of particle's center.
        
        Return
        ---
        dR: numpy.NDArray[double]
            Distance between `R` and `r_i`, with periodic boundary considered.
            Pointing from `R` to `r_i`.
        """
        dR = get_distance(R, r_i, self._lims, self._hasPeriodic, self._isPeriodic)
        return dR

    @property
    def progress_bar(self):
        """Bool flag to control whether show the process by tqdm."""
        return self._use_tqdm

    @progress_bar.setter
    def progress_bar(self, value: bool):
        self._use_tqdm = bool(value)
    
    @property
    def search_width_coef(self):
        """Coefficient of search radius, and the search radius is
        `self.search_width_coef * CG_width`.
        """
        return self._MAX_width
    
    @search_width_coef.setter
    def search_width_coef(self, value: float):
        if value <= 0:
            raise ValueError("`search_width_coef` should be positive")
        elif value <= 3:
            raise ValueError("`search_width_coef` should be larger than cut off width coefficient (3.0)")
        else:
            self._MAX_width = value

离散单元法的颗粒数据粗粒化

发表于 2024-07-30 更新于 2024-10-25 分类于离散单元法

离散单元法（Discrete Element Method，简称DEM）是颗粒流数值模拟的重要数值方法，被用于各种各样的颗粒流定量研究中。颗粒流的粗粒化（Coarse Graining），便是将颗粒流流场的大量数据转化为连续场数据的数值方法。

在开始介绍之前，先简要列出主要的参考资料：

Breard, E. C. P., Dufek, J., Fullard, L., & Carrara, A. (2020). The Basal Friction Coefficient of Granular Flows With and Without Excess Pore Pressure: Implications for Pyroclastic Density Currents, Water‐Rich Debris Flows, and Rock and Submarine Avalanches. Journal of Geophysical Research: Solid Earth, 125(12), e2020JB020203. DOI:10.1029/2020JB020203.
Lacaze, L., & Kerswell, R. R. (2009). Axisymmetric Granular Collapse: A Transient 3D Flow Test of Viscoplasticity. Physical Review Letters, 102(10), 108305. DOI:10.1103/PhysRevLett.102.108305.
Weinhart, T., Labra, C., Luding, S., & Ooi, J. Y. (2016). Influence of coarse‐graining parameters on the analysis of DEM simulations of silo flow. Powder Technology, 293, 138–148. DOI:10.1016/j.powtec.2015.11.052

单种颗粒的粗粒化

为简单起见，先从仅存在单种球形颗粒的颗粒流开始说起，并且只考虑获取空间内一点 $\boldsymbol{r}$ 在 $t$ 时刻的粗粒化数据。

粗粒化函数

在这里，引入高斯函数

W(\boldsymbol{r}) = \begin{cases} \frac{1}{V_w} \exp{\left(- \frac{\vert\boldsymbol{r}\vert^2}{2w} \right)} ,& {\vert\boldsymbol{r}\vert} < c \\ 0 ,& \text{otherwise} \end{cases}

其中： $w$ 为粗粒化宽度； $c = 3 w$ 为粗粒化截断长度； $V_w$ 是确保其密度积分等于总质量的系数

V_w = 2 \sqrt{2} \cdot w^3 \pi^{\frac{3}{2}} \mathrm{erf}\left( \frac{c \sqrt{2}}{2 w} \right) - 4 c w^2 \pi \exp{\left( \frac{- c^2}{2 w^2} \right)}

Weinhart 等建议 $w$ 的取值范围应在 $[0.75 d_{\mathrm{mean}}, 1.25 d_{\mathrm{mean}}]$ 之间，其中 $d_{\mathrm{mean}}$ 为平均颗粒直径。

密度和动量的计算

在 $t$ 时刻，单种颗粒的集合为 $q$ 。对于其中的第 $i$ 个球，其位置矢量为 $\boldsymbol{r}_i$ ，质量为 $m_i$ ，速度为 $\boldsymbol{u}_i$ 。

密度的粗粒化表达式为：

\rho(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

动量的粗粒化表达式为：

\boldsymbol{j}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i \boldsymbol{u}_i W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

所以，速度的计算需要通过 $\boldsymbol{j}$ 和 $\rho$ 完成：

\boldsymbol{u}(\boldsymbol{r}, t) = \boldsymbol{j}(\boldsymbol{r}, t) / \rho(\boldsymbol{r}, t)

根据 $W(\boldsymbol{r})$ 的定义，当 $\vert \boldsymbol{r} - \boldsymbol{r}_{i}(t) \vert \ge c$ 时， $W(\boldsymbol{r} - \boldsymbol{r}_{i}(t)) = 0$ 。
因此在实际计算中可根据这一点跳过某些球的计算。

应力张量的计算

应力张量 $[\boldsymbol{\sigma}] = \sigma_{\alpha\beta}$ 被表示为由接触作用产生的部分 $\sigma_{\alpha\beta}^{(c)}$ ，以及由动能作用产生的部分 $\sigma_{\alpha\beta}^{(k)}$ 。即：

\sigma_{\alpha\beta} = \sigma_{\alpha\beta}^{(k)} + \sigma_{\alpha\beta}^{(c)}

对于 $\sigma_{\alpha\beta}^{(k)}$ ，有：

\sigma_{\alpha\beta}^{(k)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} m_i v_{i,\alpha}^{'} v_{i,\beta}^{'} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t))

其中 $v_{i,\alpha}^{'} = u_{i,\alpha} - u_{\alpha}(\boldsymbol{r}, t)$ 为第 $i$ 个球的速度波动值。

对于 $\sigma_{\alpha\beta}^{(c)}$ ，考虑第 $i$ 个球的所有接触：

\sigma_{\alpha\beta}^{(c)}(\boldsymbol{r}, t) = \sum\limits_{i \in q} \sum\limits_{j \in \bar{Q}}^{j \neq i} F_{ij,\alpha} l_{ij,\beta} \int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s

式中

$j$ 是整个颗粒系统 $\bar{Q}$ 中的一个球（ $j \neq i$ ），并且球 $i$ 与球 $j$ 发生接触；
$\boldsymbol{F}_{ij}$ 是球 $i$ 所受到的来自球 $j$ 接触的力， $F_{ij,\alpha}$ 为其 $\alpha$ 方向分量；
$\boldsymbol{l}_{ij}$ 是球 $i$ 与球 $j$ 的接触矢量（从 $i$ 向 $j$ ）， $l_{ij,\alpha}$ 为其 $\alpha$ 方向分量。

在实际操作中，对于

\int_{0}^{1} W(\boldsymbol{r} - \boldsymbol{r}_{i}(t) + s \boldsymbol{l}_{ij}) \mathrm{d}s

我选择使用复化 Simpson 积分计算其数值，以简化操作。

多种颗粒的粗粒化

假设一个颗粒流系统 $\bar{Q}$ 中有 $Q$ 类颗粒，对于第 $q$ 类颗粒 $q \in Q$ ，可以使用单类颗粒粗粒化的公式计算出：第 $q$ 类颗粒的粗粒化密度 $\rho^{q}(\boldsymbol{r}, t)$ 、速度 $\boldsymbol{u}^{q}(\boldsymbol{r}, t)$ 和应力张量 $\sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t)$ 。

所以整个颗粒流系统 $\bar{Q}$ 的粗粒化结果为：

\begin{aligned} \rho(\boldsymbol{r}, t) &= \sum_{q \in Q} \rho^{q}(\boldsymbol{r}, t) \\ \boldsymbol{u}(\boldsymbol{r}, t) &= \sum_{q \in Q} \boldsymbol{u}^{q}(\boldsymbol{r}, t) \\ \sigma_{\alpha\beta}(\boldsymbol{r}, t) &= \sum_{q \in Q} \sigma_{\alpha\beta}^{q}(\boldsymbol{r}, t) \\ \end{aligned}

使用粗粒化数据提取颗粒流流变参数

这里以 $\mu(I)$ 流变关系为例子，以上文的文献作为参照。

记速度剪切率张量 $\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}$ ，其中 $u_i$ 源自粗粒化速度场。令：

\gamma_{ij}^{d} = \gamma_{ij} - \frac{\delta_{ij}}{3} \mathrm{tr}([\gamma])

其中 $\mathrm{tr}([\gamma]) = \gamma_{xx} + \gamma_{yy} + \gamma_{zz}$ 为 $[\gamma]$ 的迹。

所以剪切应变率定义为：

|\dot{\gamma}| = \sqrt{\frac{1}{2} {\gamma}_{ij}^{d} {\gamma}_{ij}^{d} }

定义压强为 $p = \frac{1}{3} \mathrm{tr}([\sigma]) = (\sigma_{xx} + \sigma_{yy} + \sigma_{zz}) / 3$ ，偏应力张量为： $[\tau] = \tau_{ij} = \sigma_{ij} - p \cdot \delta_{ij}$ ，所以切应力为：

|\tau| = \sqrt{\frac{1}{2} \tau_{ij} \tau_{ij} }

简单的粗粒化代码示意

我用 Python 写了一份对单种颗粒系统进行粗粒化的代码，详情见此网页。

我的 Github Page 建站记录

发表于 2024-07-28 更新于 2024-10-25 分类于杂谈

我的这个 Github Page 大约花了我一天时间从完全不懂到凑活着用。
这里简单分享一下我使用 hexo + hexo-next-theme + giscus 实现这个 Github Page 的参考资料。

建站：
从零开始用Hexo+GithubPage搭建个人网站（保姆级） - 哔哩哔哩
使用 Github Action：
如何优雅的使用Github Action服务来将Hexo部署到Github Pages - CSDN博客
 GitHub Pages | Hexo
Katex 渲染公式：
Hexo博客渲染KaTeX数学公式 | Feng’s Blog
Giscus 评论系统：
Hexo 使用 giscus 评论系统 | 404频道
 Hexo 评论系统（Giscus）

Literature Review: 简化的格子 Boltzmann 方法

发表于 2024-05-02 更新于 2024-10-25 分类于格子Boltzmann方法

注：该文章同步发在知乎专栏和CSDN
这里所说的 “简化的格子 Boltzmann 方法” 其实是 Simplified and Highly Stable Lattice Boltzmann Method (SHSLBM)(^[1], ^[2], ^[3], ^[4])。该方法的特色是仅使用平衡态分布函数、宏观密度、宏观速度进行 LBM 计算，并保证无条件稳定。
下面的内容只是我对相关文献进行简单阅读后的部分整理

单松弛格子 Boltzmann 方法

基本方程

针对流场的 LBGK 演化方程为

f_{\alpha} (\boldsymbol{r} + \boldsymbol{c}_\alpha + \Delta t, t + \Delta t) - f_{\alpha} (\boldsymbol{r}, t) = -\frac{1}{\tau} \left[ f_{\alpha} (\boldsymbol{r}, t) - f_{\alpha}^{(eq)} \right] \tag{1}

其中 $f_{\alpha} (\boldsymbol{r}, t)$ 为 $t$ 时刻在 $\boldsymbol{r}$ 处沿着 $\boldsymbol{c}_\alpha$ 方向的密度分布函数， $\tau$ 为 BGK 碰撞项的松弛时间。 $f_{\alpha}^{(eq)}$ 为 $\boldsymbol{c}_\alpha$ 方向的平衡态分布

f_{\alpha}^{(eq)} = \rho w_\alpha \left[ 1 + \frac{\boldsymbol{c}_\alpha \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_\alpha \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{\boldsymbol{u} \cdot \boldsymbol{u}}{2 c_s^2} \right] \tag{2}

式 (2) 中 $\rho$ 为流体密度， $\boldsymbol{u} = (u_x, u_y, u_z)^T$ 为流体宏观速度， $w_\alpha$ 为 $\boldsymbol{c}_\alpha$ 的权重系数， $c_s$ 为格子声速。 $w_\alpha$ 、 $\boldsymbol{c}_\alpha$ 和 $c_s$ 的数值由选取的离散速度模型决定。

对于 LBM 而言，宏观量表示为

\rho = \sum_{\alpha} f_{\alpha} = \sum_{\alpha} f_{\alpha}^{(eq)}, \quad \rho\boldsymbol{u} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha} = \sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \tag{3}

Chapman-Enskog 展开

Boltzmann 方程和 LBGK 方程都可以通过多尺度展开还原为 Navier-Stokes 方程。这里以 Chapman-Enskog 为例：

f_{\alpha} = f_{\alpha}^{(eq)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ... , \quad \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_0} + \epsilon^2 \frac{\partial}{\partial t_1}, \quad \frac{\partial}{\partial \boldsymbol{r}} = \epsilon \frac{\partial}{\partial \boldsymbol{r}_0}

上述展开保留至 $\epsilon^{1}$ ，即可还原为^[5]

\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0 \tag{4}

和

\frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \Pi = \boldsymbol{0} \tag{5}

式 (5) 中张量 $\Pi$ 表示为

\Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) \epsilon f_{\alpha}^{(1)} \right]

由于只保留到 $\epsilon^{1}$ ，所以 $\epsilon f_{\alpha}^{(1)} \approx f_{\alpha} - f_{\alpha}^{(eq)} = f_{\alpha}^{(neq)}$ 近似为非平衡态分布函数。此时

\Pi_{\beta \gamma} = \sum_{\alpha} (\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} \left[ f_{\alpha}^{(eq)} + \left(1 - \frac{1}{2 \tau}\right) f_{\alpha}^{(neq)} \right] \tag{6}

并且存在如下关系

f_{\alpha}^{(neq)} \approx \epsilon f_{\alpha}^{(1)} = -\tau \left[ \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla \right] f_{\alpha}^{(eq)} \Delta t \tag{7}

令 $\mathrm{D} = \frac{\partial}{\partial t} + \boldsymbol{c}_{\alpha} \cdot \nabla$ ，则 $f_{\alpha}^{(neq)} \approx -\tau \cdot (\Delta t) \cdot \mathrm{D}f_{\alpha}^{(eq)}$ 。

流体运动粘度 $\nu$ 和松弛时间 $\tau$ 的关系为

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

SHSLBM

式 (7) 说明：在保留至 $\epsilon^{1}$ 的情况下， $f_{\alpha}^{(neq)}$ 和 $f_{\alpha}^{(eq)}$ 存在联系。考虑到式 (1) 中由 $\frac{-1}{\tau}$ 松弛的部分即为 $f_{\alpha}^{(neq)}$ ，因此式 (1) 理应可以被写作完全由 $f_{\alpha}^{(eq)}$ 进行表示的形式。这便是 SHSLBM 的核心思想。

基本算法

SHSLBM 的单个时间步迭代分为预报和校正两部分，具体写作：

预报步：计算 $(\boldsymbol{r},t)$ 的预估密度 $\rho^{*}$ 和预估速度 $\boldsymbol{u}^{*}$

\rho^{*} = \sum_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t), \quad \rho^{*} \boldsymbol{u}^{*} = \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) \tag{8}

校正步：根据 $\rho^{*}$ 和 $\boldsymbol{u}^{*}$ 得出校正值

\begin{cases} \rho = \rho^{*} \\ \rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})] \end{cases} \tag{9-1}

式 (9-1) 中非平衡态的计算通过对式 (7) 的差分实现，即：

f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t) \right] + O((\Delta t)^3)

同理：

f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) = -\tau \left[ f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) - f_{\alpha}^{(eq)}(\boldsymbol{r}, t - \Delta t) \right] + O((\Delta t)^3)

在这两条式子中， $f_{\alpha}^{(eq)}(\boldsymbol{r}, t)$ 和 $f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t)$ 都是未知的。参照Shu等^[6]的做法，对上面两条公式中的这两个量使用预报步数据进行近似，即

f_{\alpha}^{(eq)}(\boldsymbol{r}, t) = f_{\alpha}^{(eq)}(\rho^{*}(\boldsymbol{r}, t), \boldsymbol{u}^{*}(\boldsymbol{r}, t))

式 (9-1) 中 $\rho \boldsymbol{u}$ 的计算还可以进一步化简。注意到，对于 $f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2})$ ，有

\begin{aligned} \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r} - \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= -\tau \left[ \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \right.\\ &\quad\left. -\sum\limits_{\alpha} \boldsymbol{c}_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t - \Delta t)\right] + O((\Delta t)^3) \\ &= -\tau (\rho^{*} \boldsymbol{u}^{*} - \rho^{*} \boldsymbol{u}^{*}) + O((\Delta t)^3) \\ &= O((\Delta t)^3) \end{aligned}

所以式 (9-1) 可进一步化简为式 (9-2)

\begin{cases} \rho = \rho^{*} \\ \rho \boldsymbol{u} = \rho^{*} \boldsymbol{u}^{*} - \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) \end{cases} \tag{9-2}

对应的宏观流程

事实上，式 (8) 和式 (9-1), (9-2) 的预报和校正存在着对应的宏观偏微分方程。下面简要说明一下。

预报步

式 (8) 对应的宏观方程为

\begin{cases} \frac{\partial \rho}{\partial t} + \nabla \cdot \left(\sum\limits_{\alpha} \boldsymbol{c}_\alpha f_{\alpha}^{(eq)} \right) = 0 \\ \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(eq)} + \frac{1}{2 \tau} \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0} \end{cases} \tag{10}

首先，对 $f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t)$ 在 $(\boldsymbol{r},t)$ 进行展开，并代入式 (7) ，得到：

\begin{aligned} f_{\alpha}^{(eq)} (\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t, t -\Delta t) &= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\ &\quad+ \frac{(\Delta t)^2}{2} \mathrm{D}^{2} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\ &= f_{\alpha}^{(eq)} (\boldsymbol{r}, t) - (\Delta t) \mathrm{D} f_{\alpha}^{(eq)} (\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \mathrm{D} f_{\alpha}^{(neq)} (\boldsymbol{r}, t) + O((\Delta t)^3) \\ \end{aligned} \tag{11}

将式 (11) 代入式 (8) 的第一条公式，有：

\begin{aligned} \rho^{*} &= \sum_{\alpha}f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \end{aligned} \tag{12}

式 (12) 右侧第一项与 $\rho^{*}$ 相等，右侧第三项根据非平衡态的守恒律可知其为 0 。所以在消去式 (12) 中的 $\Delta t$ 后， $\rho^{*}$ 的计算与式 (10) 中第一条公式对应，并具有 $O((\Delta t)^2)$ 精度。

将式 (11) 代入式 (8) 的第二条公式，整理得：

\begin{aligned} \rho^{*} \boldsymbol{u}^{*} &= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \boldsymbol{c}_{\alpha} \mathrm{D} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \\ &= \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) - (\Delta t) \left[ \frac{\partial}{\partial t} \sum_{\alpha} \boldsymbol{c}_{\alpha} f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \right. \\ &\quad \left. \nabla \cdot \sum_{\alpha} (\boldsymbol{c}_{\alpha})_\beta (\boldsymbol{c}_{\alpha})_\gamma \left( f_{\alpha}^{(eq)}(\boldsymbol{r}, t) + \frac{1}{2 \tau} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) \right) \right] \\ &\quad - \frac{\Delta t}{2 \tau} \sum_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t) + O((\Delta t)^3) \end{aligned} \tag{13}

式 (13) 右侧第一项与 $\rho^{*} \boldsymbol{u}^{*}$ 相等，右侧 $\sum\limits_{\alpha} \frac{\partial}{\partial t} \boldsymbol{c}_{\alpha} f_{\alpha}^{(neq)}(\boldsymbol{r}, t)$ 项根据非平衡态的守恒律可知其为 0 。所以在消去式 (13) 中的 $\Delta t$ 后， $\rho^{*} \boldsymbol{u}^{*}$ 的计算与式 (10) 中第二条公式对应，并具有 $O((\Delta t)^2)$ 精度。

校正步

式 (9) 对应的宏观方程为

\begin{cases} \frac{\partial \rho}{\partial t} = 0 \\ \frac{\partial \rho \boldsymbol{u}}{\partial t} + \nabla \cdot \left[ \left( 1 - \frac{1}{2 \tau} \right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)} \right] = \boldsymbol{0} \end{cases} \tag{14}

这里主要介绍一下动量方程的还原。对 $f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2})$ 进行 Taylor 展开，得到

\begin{aligned} f_{\alpha}^{(neq)}(\boldsymbol{r} \pm \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2}) &= f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\ &\quad \pm \frac{\Delta t}{2} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \\ &\quad + \frac{(\Delta t)^2}{8} (\boldsymbol{c}_{\alpha} \cdot \nabla)^2 f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^3) \end{aligned}

这里我们以未化简的式 (9-1) 为例。在式 (9-1) 的动量校正步中，有：

\begin{aligned} &\quad \left(1 - \frac{1}{\tau}\right) \frac{1}{\Delta t} \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} [f_{\alpha}^{(neq)}(\boldsymbol{r}+\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2}) - f_{\alpha}^{(neq)}(\boldsymbol{r}-\frac{\Delta t}{2}\boldsymbol{c}_{\alpha},t-\frac{\Delta t}{2})] \\ &= \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha} \boldsymbol{c}_{\alpha} (\boldsymbol{c}_{\alpha} \cdot \nabla) f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) + O((\Delta t)^2) \\ &= \nabla \cdot \left[ \left(1 - \frac{1}{\tau}\right) \sum\limits_{\alpha}(\boldsymbol{c}_\alpha)_{\beta} (\boldsymbol{c}_\alpha)_{\gamma} f_{\alpha}^{(neq)}(\boldsymbol{r}, t - \frac{\Delta t}{2}) \right] + O((\Delta t)^2) \\ \end{aligned}

此时联立式 (13) ，即可还原式 (14) 中的动量方程。

SHSLBM 的稳定性分析

这里令 $\mathbf{Y} = (y_1, y_2, y_3, y_4)^T$ 。其中

y_1 = \rho,\quad y_2 = \rho u_x ,\quad y_3 = \rho u_y ,\quad y_4 = \rho u_z

将这些符号代换到式 (2) 的 $f_{\alpha}^{(eq)}$ 中，得：
$\begin{aligned} f_{\alpha}^{(eq)} &= y_{1 \alpha} w_{\alpha} + \frac{w_\alpha}{c_s^2} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} ) \\ &\quad + \frac{w_\alpha}{2 y_{1 \alpha} c_s^4} ( c_{\alpha x} y_{2 \alpha}+ c_{\alpha y} y_{3 \alpha} + c_{\alpha z} y_{4 \alpha} )^2 \\ &\quad - \frac{w_\alpha}{2 y_{1 \alpha} c_s^2} (y_{2 \alpha}^2 + y_{3 \alpha}^2 + y_{4 \alpha}^2) \end{aligned} \tag{15}$
where the subscript $\alpha$ denotes the quantities at a space level of $(\boldsymbol{r} - \boldsymbol{c}_{\alpha} \Delta t)$ .^[2:1]

预报步

对于第 $n$ 步的结果 $\mathbf{Y}^{n}$ ，式 (8) 的预报步写作 $\mathbf{Y}^{*} = \mathbf{G}(\mathbf{Y}^{n})$ 。

为便于进行诺依曼稳定性分析，对上式进行求导，有

\mathrm{d \mathbf{Y}^{*}} = \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_1} \right]^n \mathrm{d}y_1^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_2} \right]^n \mathrm{d}y_2^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_3} \right]^n \mathrm{d}y_3^n + \left[ \frac{\partial \mathbf{Y}^{*}}{\partial y_4} \right]^n \mathrm{d}y_4^n

上标 $n$ 代表第 $n$ 个时间步。由于 $\mathrm{d \mathbf{Y}^{*}} = (\mathrm{d} y_1^{*}, \mathrm{d} y_2^{*}, \mathrm{d} y_3^{*}, \mathrm{d} y_4^{*})^T$ ，所以，对于预报步结果的每个分量而言有： $\mathrm{d} y_{\lambda}^{*} = \sum_{\kappa = 1}^{4}\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n \mathrm{d}y_{\kappa}^n$ （ $\lambda, \kappa \in [1,2,3,4]$ ）。

对于本文列出的方程组， $\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n$ 多达 16 个。因此下文仅拿 $\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n$ 举例进行计算。

\begin{aligned} \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n &= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\ &= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^4} ( c_{\alpha x} y_{2 \alpha}^n + c_{\alpha y} y_{3 \alpha}^n + c_{\alpha z} y_{4 \alpha}^n )^2 \right. \\ &\quad\left. + \frac{w_\alpha}{2 (y_{1\alpha}^n)^2 c_s^2} ((y_{2 \alpha}^n)^2 + (y_{3 \alpha}^n)^2 + (y_{4 \alpha}^n)^2) \right] \end{aligned} \tag{16}

假设式 (16) 中的 $\mathrm{d} y^{*}_{i}$ , $\mathrm{d} y^{n}_{i}$ , $y^{n}_{i \alpha}$ 均可被展开为如下形式（ $i = [1,2,3,4]$ ， $j=\sqrt{-1}$ ）

\mathrm{d} y^{*}_{i} = \overline{\mathrm{d} y_{i}}^{*} \exp(j \boldsymbol{k \cdot r}) ,\quad \mathrm{d} y^{n}_{i} = \overline{\mathrm{d} y_{i}}^{n} \exp(j \boldsymbol{k \cdot r}) ,\quad y^{n}_{i \alpha} = y_{i}^{n} \exp(j \boldsymbol{k \cdot r})

代入到式 (16) 和 $\mathrm{d \mathbf{Y}^{*}}$ 表达式，得

\begin{aligned} \left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n &= \frac{\partial}{\partial y_{1}} \sum_{\alpha} f_{\alpha}^{(eq)} \\ &= \sum_{\alpha} \left[ w_{\alpha} - \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^4} ( c_{\alpha x} y_{2}^n + c_{\alpha y} y_{3}^n + c_{\alpha z} y_{4}^n )^2 \right. \\ &\quad\left. + \frac{w_\alpha}{2 (y_{1}^n)^2 c_s^2} ((y_{2}^n)^2 + (y_{3}^n)^2 + (y_{4}^n)^2) \right] \end{aligned} \tag{17}

格子张量的一些性质^[2:2]
$\begin{aligned} & \sum_{\alpha} w_{\alpha} = 1 .&\qquad \alpha \in [0,1, ..., q-1]\\ & \sum_{\alpha} w_{\alpha} c_{\alpha \beta} = 0 .&\qquad \beta \in [x,y,z] \\ & \sum_{\alpha} w_{\alpha} c_{\alpha \beta} c_{\alpha \gamma} = c_s^2 \delta_{\beta \gamma} .&\qquad \beta, \gamma \in [x,y,z] \end{aligned}$
注: $q$ 为离散速度数目

基于格子张量的性质，式 (17) 右侧仅剩下 $\sum\limits_{\alpha} w_{\alpha}$ ，即：

\left[ \frac{\partial y_{1}^{*}}{\partial y_{1}} \right]^n = 1

同理，其他项也可进行类似分析，并得出

\left[ \frac{\partial y_{\lambda}^{*}}{\partial y_{\kappa}} \right]^n = \delta_{\lambda \kappa}

其中 $\delta_{\lambda \kappa}$ 为狄拉克函数，且 $\lambda,\kappa \in [1,2,3,4]$ 。所以我们可以得到 $\overline{\mathrm{d} \mathbf{Y}}^{*} = \overline{\mathrm{d} \mathbf{Y}}^{n}$ 。

校正步

令 ${\mathrm{d} y_i}^{n+1} = \overline{\mathrm{d} y_i}^{n+1} \exp(j \boldsymbol{k \cdot r})$ （ $i \in [1,2,3,4]$ ）。由于校正步（式 (9-2) ）中并未对 $\rho$ 做修改，易得：

\overline{\mathrm{d} y_1}^{n+1} = \overline{\mathrm{d} y_1}^{*} = \overline{\mathrm{d} y_1}^{n}

根据式 (9-2) 的校正步，以及 $f_{\alpha}^{(neq)}(\boldsymbol{r} + \frac{\Delta t}{2} \boldsymbol{c}_{\alpha}, t - \frac{\Delta t}{2})$ 的差分表示，有

\begin{aligned} y_2^{n+1} = y_2^{*} - (\tau - 1) y_2^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ y_3^{n+1} = y_3^{*} - (\tau - 1) y_3^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha y} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ y_4^{n+1} = y_4^{*} - (\tau - 1) y_4^{n} + (\tau - 1) \sum\limits_{\alpha} c_{\alpha z} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \end{aligned} \tag{18}

这里同样以 $y_2^{n+1}$ 的计算为例，记

\begin{aligned} S &= \sum\limits_{\alpha} c_{\alpha x} f_{\alpha}^{(eq)}(\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t, t) \\ &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} y_{1(-\alpha)}^{*} \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{c_s^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*}) \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 y_{1(-\alpha)}^{*}} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 y_{1(-\alpha)}^{*}} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2 \\ \end{aligned}

其中下标 $(-\alpha)$ 代表来自 $\boldsymbol{r} + \boldsymbol{c}_{\alpha} \Delta t$ 的值。则

y_2^{n+1} = y_2^{*} + (\tau - 1) (S - y_2^{n}) \tag{19}

所以其差分表示为

\mathrm{d} y_2^{n+1} = \mathrm{d} y_2^{*} + (\tau - 1) (\mathrm{d} S - \mathrm{d} y_2^{n}) \tag{20}

对于 $\mathrm{d} S$ ，有

\mathrm{d} S = \left( \frac{\partial S}{\partial y_1} \right)^{*} \mathrm{d}y_1^{*} + \left( \frac{\partial S}{\partial y_2} \right)^{*} \mathrm{d}y_2^{*} + \left( \frac{\partial S}{\partial y_3} \right)^{*} \mathrm{d}y_3^{*} + \left( \frac{\partial S}{\partial y_4} \right)^{*} \mathrm{d}y_4^{*}

以 $\left( \frac{\partial S}{\partial y_1} \right)^{*}$ 为例，有

\begin{aligned} \left( \frac{\partial S}{\partial y_1} \right)^{*} &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1(-\alpha)}^{*})^2} (c_{\alpha x} y_{2(-\alpha)}^{*} + c_{\alpha y} y_{3(-\alpha)}^{*} + c_{\alpha z} y_{4(-\alpha)}^{*})^2 \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1(-\alpha)}^{*})^2} ( (y_{2(-\alpha)}^{*})^2 + (y_{3(-\alpha)}^{*})^2 + (y_{4(-\alpha)}^{*})^2 )^2 \end{aligned}

与预报步的做法类似，将 $y_{i(-\alpha)}^{*}$ 近似为 $y_{i}^{*}$ ，得到：

\begin{aligned} \left( \frac{\partial S}{\partial y_1} \right)^{*} &= \sum\limits_{\alpha} c_{\alpha x} w_{\alpha} \\ &\quad - \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^4 (y_{1}^{*})^2} (c_{\alpha x} y_{2}^{*} + c_{\alpha y} y_{3}^{*} + c_{\alpha z} y_{4}^{*})^2 \\ &\quad + \sum\limits_{\alpha} \frac{c_{\alpha x} w_{\alpha}}{2 c_s^2 (y_{1}^{*})^2} ( (y_{2}^{*})^2 + (y_{3}^{*})^2 + (y_{4}^{*})^2 )^2 \end{aligned}

此时所有 $y_{i}^{*}$ 均可被视为常数提取出来。根据格子张量的性质，可得： $\left( \frac{\partial S}{\partial y_1} \right)^{*} = 0$ 。

同理

\left( \frac{\partial S}{\partial y_2} \right)^{*} = 1,\quad \left( \frac{\partial S}{\partial y_3} \right)^{*} = \left( \frac{\partial S}{\partial y_4} \right)^{*} = 0

此时 $\overline{\mathrm{d} S } = \overline{\mathrm{d} y_2}^{*}$ ，则式 (20) 可化简为：

\overline{\mathrm{d} y_2}^{n+1} = \overline{\mathrm{d} y_2}^{*} = \overline{\mathrm{d} y_2}^{n}

同理，对 $y_3^{n+1}$ 和 $y_4^{n+1}$ 计算可得： $\overline{\mathrm{d} y_3}^{n+1} = \overline{\mathrm{d} y_3}^{n}$ 和 $\overline{\mathrm{d} y_4}^{n+1} = \overline{\mathrm{d} y_4}^{n}$ 。

记 $\mathrm{d \mathbf{Y}^{n}} = (\mathrm{d} y_1^{n}, \mathrm{d} y_2^{n}, \mathrm{d} y_3^{n}, \mathrm{d} y_4^{n})^T$ ， $\mathrm{d \mathbf{Y}^{n+1}} = (\mathrm{d} y_1^{n+1}, \mathrm{d} y_2^{n+1}, \mathrm{d} y_3^{n+1}, \mathrm{d} y_4^{n+1})^T$ ，那么：

\overline{\mathrm{d \mathbf{Y}}}^{n+1} = \mathbf{I} \, \overline{\mathrm{d \mathbf{Y}}}^{n}

其中 $\mathbf{I}$ 为单位矩阵。显然，特征矩阵的特征值全是 1 ，代表着该格式无条件稳定。

Chen, Z., Shu, C., Wang, Y., Yang, L. M., & Tan, D. (2017). A Simplified Lattice Boltzmann Method without Evolution of Distribution Function. Advances in Applied Mathematics and Mechanics, 9(1), 1–22. DOI:10.4208/aamm.OA-2016-0029 ↩︎
Z. Chen, C. Shu, D. Tan; Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids. 1 May 2017; 29 (5): 053601. DOI:10.1063/1.4983339 ↩︎ ↩︎ ↩︎
Chen Z, Shu C, Tan D, Wu C. On improvements of simplified and highly stable lattice Boltzmann method: Formulations, boundary treatment, and stability analysis. Int J Numer Meth Fluids. 2018; 87: 161–179. DOI:10.1002/fld.4485 ↩︎
Chen, Z., & Shu, C. (2020). Simplified and Highly Stable Lattice Boltzmann Method. WORLD SCIENTIFIC. DOI:10.1142/12047 ↩︎
本文中所有 $\nabla$ 均表示空间偏导 $\frac{\partial}{\partial \boldsymbol{r}}$ 。 ↩︎
C. Shu, Y. Wang, C. Teo, and J. Wu. Development of lattice Boltzmann flux solver for simulation of incompressible flows. Advances In Applied Mathematics And Mechanics. 6, 436 (2014). ↩︎

从格子Boltzmann方程到宏观流体问题的常用还原方法

发表于 2023-04-08 更新于 2024-10-25 分类于格子Boltzmann方法

这是我之前的一次期末大报告作业，里面的内容也算是从其他文章里搬运的，在此重新誊抄出来。内容如有错误，欢迎各位指正。

[注] 1. 下文提到的 $x$ 、 $u$ 和 $\xi$ 无论有没有加粗均为向量。
2. 本文已经同步发表在本人的知乎和CSDN上，可在其他网站获得同样的阅读效果。

一、单松弛格子Boltzmann方程

1.1 Boltzmann方程

格子Boltzmann方法基于Boltzmann方程（即(1)式），该方程引入速度分布函数，描述了非热力学平衡状态的热力学系统统计行为：

\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \boldsymbol{x}} + \frac{F}{m} \cdot \frac{\partial f}{\partial \mathbf{\xi}}=\Omega(f) \tag{1}

(1)式中 $\Omega(f)$ 为碰撞项， $\mathbf{\xi}=\frac{\partial x}{\partial t}$ 。速度分布函数的定义是：在时刻 $t$ ，以 $x$ 为中心的空间微元 $\mathrm{d}\mathit{x}$ 内，速度在 $\xi$ 和 $(\xi+\mathrm{d}\it{\xi})$ 之间的分子的数量为 $f \mathrm{d}\it{x} \mathrm{d}\it{\xi}$ 。记分子质量为 $m$ ，流体宏观速度为 $u$ ，单位体积内能为 $e$ ，分子数密度为 $n$ ，则：

\rho = mn = m \int{f} \mathrm{d}\xi

\rho u = m \int{\xi f} \mathrm{d}\xi

\rho E = \rho e + \frac{1}{2} \rho u^2 = m \int{ f \frac{\xi^2}{2} } \mathrm{d} \xi

碰撞项描述了分子碰撞对整个系统运动产生的影响，它具有两个特点：①满足质量、动量和能量守恒；②能够反映系统趋于平衡态的趋势。Bhatnagur、Gross和Krook三人在此基础上提出了最为简单且著名的BGK模型，即(2)式。

\Omega(f) = -\frac{1}{\tau_c} \left ( f - f^{(eq)} \right ) \tag{2}

其中 $\tau_c$ 为松弛时间， $f^{(eq)}$ 为平衡态分布函数。可根据Boltzmann H定理计算得：

f^{(eq)} = \frac{n}{(2 \pi RT)^{3/2}} \mathrm{exp}\left( -\frac{(\mathbf{\xi-u})^2}{2RT} \right) \tag{3}

1.2 格子Boltzmann方程

为了将Boltzmann方程应用于实际的数值模拟，需要将流体视为大量的离散粒子。离散粒子于流体分子不同，它们驻留在网格上，并按一定的规则在网格上发生碰撞和迁移。因此，1.2节将 $mf$ 记为 $f$ ，此时 $f$ 为密度分布函数。对密度分布函数 $f$ ，不考虑 $F$ 的Boltzmann-BGK方程为式。

\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \mathit{x}} = - \frac{1}{\tau_c} \left[ f - f^{(eq)} \right] \tag{4}

其宏观量按如下方式确定

\rho = \int { f } \mathrm{d} \xi ，\rho u = \int { \xi f } \mathrm{d} \xi ，\rho E = \int { f \frac{\xi^2}{2} } \mathrm{d} \xi

对 $f$ 和其对应的 $f^{(eq)}$ 可进行如下离散： $f_i = w_i f$ ， $f_i^{(eq)} = w_i f^{(eq)}$ 。其中下标 $i$ 表示特征方向 $\boldsymbol{c}_i$ ， $f_i$ 和 $f_i^{(eq)}$ 均沿着 $\boldsymbol{c}_i$ 方向运动。据此，沿着特征方向离散(4)式得到的差分格式为

f_i (x+ξ_i Δt,t+Δt)-f_i (x,t)=-\frac{1}{τ} \left[f_i (x,t)-f_i^{(eq)} (x,t) \right] \tag{5}

其中 $\tau = \frac{\tau_c}{\Delta t}$ 为无量纲松弛时间。由于对应的 $f^{(eq)}$ 为：

f^{(eq)} = \rho \cdot \mathrm{exp} ( - \frac{(\xi - \it{u})^2}{2\it{RT}} ) \tag{6-1}

所以将 $f^{(eq)}$ 展开至2阶并离散，可得：

f^{(eq)}_i = w_i \rho \left[ 1 + \frac{\xi_i \cdot u}{c_s^2} + \frac{(\xi_i \cdot u)^2}{2 c_s^4} - \frac{u \cdot u}{2 c_s^2} \right] \tag{6-2}

其中 $c_s = \sqrt{RT}$ 在LBM模拟中一般被称作格子声速。

二、Chapman-Enskog展开分析

将基本碰撞不变量 $φ_i=1,ξ,\frac{ξ^2}{2}$ 乘以(1)式两端，并对 $\xi$ 积分有：

\frac{∂ρ}{\partial t}+∇⋅(ρ \mathbf{u})=0 \tag{7-1}

\frac{∂(ρu)}{\partial t}+∇⋅(ρ \mathbf{uu}+P) = ρ \mathit{\bf{a}} \tag{7-2}

\frac{∂(\rho E)}{\partial t} + ∇⋅(\rho \mathbf{u}E + P⋅\mathbf{u} + \mathbf{Q})=\rho \mathbf{a⋅u} \tag{7-3}

Chapman-Enskog展开（下文简称C-E展开）是处理Boltzmann方程的一种分析方法，该方法引入小量 $\varepsilon$ 作为展开因子（一般与Knudsen数同阶），将BGK方程在不同的尺度上展开，即：

f = \sum_{n=0}^{\infty} \varepsilon^n f^{(n)} \tag{8}

并将分布函数、导数、物理量等都按照 $\varepsilon$ 的不同阶次展开。该方法对 $f^{(n)}$ 进行了限制，使其满足： $\int f^{(n)} φ_i \mathrm{d}ξ=0$ ， $φ_i (i=1,..,5)$ 为基本碰撞不变量。将 $f$ 视为 $\rho$ 、 $\mathbf{u}$ 和 $T$ 及其各阶导数的泛函，而时间偏导项则展开为

\frac{\partial }{\partial t} = \sum_{n=0}^{\infty } \frac{\partial _n}{\partial t} \tag{9}

将(8)式和(9)式代入到(7)式中。得 $O({\varepsilon}^0)$ 方程组为：

\frac{∂_0 ρ}{\partial t}+∇⋅(ρ \mathbf{u})=0 \tag{10-1}

\frac{∂_0 (ρu)}{\partial t}+∇⋅(ρ \mathbf{uu}+P^{(0)}) = ρ \mathbf{a} \tag{10-2}

\frac{∂_0(\rho E)}{\partial t} + ∇⋅(\rho \mathbf{u}E + P^{(0)} ⋅ \mathbf{u} + \mathbf{Q}^{(0)})=\rho \mathbf{a⋅u} \tag{10-3}

以及 $O({\varepsilon}^n)$ 方程组为 $(n>0)$ ：

\frac{∂_n ρ}{\partial t}=0 \tag{11-1}

\frac{∂_n (ρu)}{\partial t}+∇ ⋅ P^{(n)} = 0 \tag{11-2}

\frac{∂_n (\rho E)}{\partial t} + ∇⋅( P^{(n)} ⋅ \mathbf{u} + \mathbf{Q}^{(n)} ) = 0 \tag{11-3}

其中

C=\boldsymbol{ξ-u}, P^{(n)}=m∫CCf^{(n)} \mathrm{d}ξ , Q^{(n)}=m∫C \frac{C^2}{2} f^{(n)} \mathrm{d}ξ

只需取精确到 $O({\varepsilon}^1)$ 的宏观方程组进行分析，即可还原出Naiver-Stokes方程。且具有如下关系：

P=P^{(0)} + \varepsilon P^{(1)} = p \mathbf{I} - 2 \varepsilon \mu (\mathbf{S} - \frac{\mathrm{Tr}(S)}{3} \mathbf{I})

Q = Q^{(0)} + \varepsilon Q^{(1)} = - \varepsilon \kappa \nabla{T}

其中： $\kappa$ 为热传导系数， $\mu$ 为粘性系数。

离散的LBGK方程也可采用类似的方法展开，还原出宏观流体运动方程（见附录A）。展开方式如下：

f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)} + ...

\frac{\partial}{\partial x} = \varepsilon \frac{\partial}{\partial x_1} , x_1 = \varepsilon x

略去 $O({\varepsilon}^3)$ 项，并对 $O({\varepsilon}^0)$ 项、 $O({\varepsilon}^1)$ 项和 $O({\varepsilon}^2)$ 项求0阶和1阶速度矩，即可得到不可压缩流体的Navier-Stokes方程。由于在分析过程中略去了 $O({\varepsilon}^3)$ 项和速度的3阶项，导致结果缺失了部分物理信息。

三、 Grad-13矩方法及其发展

3.1 Grad的矩分析方法

Grad -13矩是一种由数学家Harold Grad提出的分析方法，其基于Hermite多项式进行分析。对于一个长度为 $N$ 的向量 $\boldsymbol{x}$ ，若定义权函数

\omega( \boldsymbol{x}) = \frac{1}{(2 \pi)^{N/2}} \mathrm{exp} \left(- \frac{ \boldsymbol{x} \cdot \boldsymbol{x}}{2} \right) \tag{12}

则其对应的 $n$ 阶Hermite多项式写作

\mathcal{H} ^{(n)} ( \boldsymbol{x}) = \frac{(-1)^n}{\omega ( \boldsymbol{x})} \nabla^{n} \omega ( \boldsymbol{x}) \tag{13}

由Hermite多项式的正交性，可证明得到：

∫_{-∞}^{+∞} \omega ( \boldsymbol{x}) H_i^{(m)}( \boldsymbol{x}) H_j^{(n)}( \boldsymbol{x}) \mathrm{d} \boldsymbol{x} = \delta_{mn} \delta_{ij} \tag{14}

据此，可将任意一个函数 $f(x)$ 使用Hermite多项式进行展开^[1]。即：

f( \boldsymbol{x})=ω^{\frac{1}{2}} \sum_{n=0}^∞ \sum_{i=1}^N a_i^{(n)} H_i^{(n)}( \boldsymbol{x}) \tag{15}

所以

\int_{-\infty}^{-\infty} \omega^{1/2} f(x) \mathcal{H}_j^{(m)} (x) \mathrm{d} \it{x} = \sum_{i=\mathrm{1}}^{N} a_i^{\mathrm{(}m\mathrm{)}} \delta_{ij}^m = m\mathrm{!} a_i^{m} \mathrm{\tag{16}}

其中 $a_i^{(m)}$ 可借由Hermite多项式的正交性进行求解。 $\delta_{ij}^{m}$ 是一个记号，表示 $m$ 个 $\delta_{ij}$ 的乘积之和， $i$ 和 $j$ 分别来自下标集 $(i_1 , i_2 , ... , i_m)$ 和 $(j_1 , j_2 , ... , j_m)$ 。
在此数学基础之上，Harold Grad将分布函数 $f$ 展开为^[2]：

f(x, \xi, t) = \omega \sum_{n=0}^{\infty} \frac{1}{n!} a^{(n)} H^{(n)}(\xi) \tag{17}

根据式(16)， $a^{(n)}$ 的表达式为：

a^{(n)} = \int_{-\infty}^{+\infty} f(x, \xi, t) H^{(n)}(\xi) d\xi \tag{18}

可得 $a^{(n)}$ 的前4项为：

a^{(0)} = \int_{-\infty}^{+\infty} f(x, \xi, t) d\xi= \rho \tag{19-1}

a^{(1)} = \int_{-\infty}^{+\infty} f(x, \xi, t) \xi d\xi = \rho u \tag{19-2}

a^{(2)} = \int_{-\infty}^{+\infty} f(x, \xi, t) (\xi^2 - \delta) d\xi = \rho (\mathbf{u}^2 - \mathbf{\delta}) + \mathbf{P} \tag{19-3}

a^{(3)} = \int_{-\infty}^{+\infty} f(x, \xi, t) (\xi^3 - \xi \delta) d\xi = (D-1) \rho \bf{u}^3 + \bf{u} \it{a}^{\mathrm{(2)}} \mathrm{+} \mathbf{Q} \tag{19-4}

注意到从 $a^{(0)}$ 到 $a^{(3)}$ 都描述了流场的宏观量，所以通过上述表达式的值可以解出流场宏观量。与C-E分析方法不同的是分析过程中没有舍去高阶量。因此这种展开方法得到的式子合理地引入了更高阶的矩，物理意义更加明确。

3.2 Gauss-Hermite积分的引入和离散化BGK方程

Grad的思路开创了描述的新方式，而单肖文、袁学锋、陈沪东三位学者在此基础上进一步发展，于2006年在 Journal of Fluid Mechanics 上提出了将原思路拓展到离散的Boltzmann-BGK方程中的方法^[3]，简化了求解过程。
首先，基于式(16)将 $f$ 截断至前 $N$ 阶，得：

f \approx f^N = \omega \sum_{n=0}^{N} \frac{1}{n!} a^{(n)} H^{(n)}(\xi) \tag{20}

$f^N$ 和 $f$ 有着完全相同的前 $N$ 阶矩。通过截断操作将BGK方程对宏观流动方程的近似保持在动力学层面而不是热力学层面。此时 $f^N$ 可表示为：

f^{N}(x, \xi , t) H^{(n)}(\xi) = \omega(\xi) p(x, \xi, t) \tag{21}

其中 $p(x, \xi, t)$ 为阶数不超过 $2N$ 的多项式。使用Gauss-Hermite积分可将 $a^{(n)}$ 精确地展开为 $p(x, \xi, t)$ 的加权和。

a^{(n)} =\int \omega(\xi) p(x, \xi, t) d \xi = \sum_{a=1}^{d} w_a p(x, \xi_a, t) = \sum_{a=1}^{d} \frac{w_a}{\omega(\xi_a)} f^{N}(x, \xi_a, t) H^{(n)}(\xi_a) \tag{22}

式(22)中 $w_a$ 和 $ \boldsymbol{\xi}_a$ 是Gauss-Hermite积分的权重和特征向量（ $a=1,2,...,d$ ）。所以使用 $f^N (x, \boldsymbol{\xi}_a, t)$ 表示 $f^N$ 的前 $N$ 阶矩在数学上可行。
据此可定义离散化的分布函数 $f_a$ （ $a=1,2,...,d$ ），其对应的权重和特征方向为 $w_a$ 和 $\xi_a$ ：

f_a (x,t) = \frac{w_a}{\omega( \boldsymbol{\xi}_a)} f(x, \boldsymbol{\xi}_a, t) \tag{23}

则宏观量可表达为：

\rho = \sum_{a=1}^{d} f_a , \rho \mathbf{u} = \sum_{a=1}^{d} f_a \mathbf{\xi}_a , \rho \mathbf{uu} + \mathbf{P} = \sum_{a=1}^{d} f_a \mathbf{\xi}_a \mathbf{\xi}_a \tag{24}

下面将(1)式的碰撞项换为(2)式，并把 $\frac{\mathbf{F}}{m} \frac{\partial f}{\partial \xi}$ 直接记作外力项 $-\mathbf{F} (\xi)$ ，得

\frac{\partial f}{\partial t} + \mathbf{\xi} \cdot \frac{\partial f}{\partial \mathit{x}} = \frac{-1}{\tau} ( f - f^{(eq)} ) + \mathbf{F} (\xi) \tag{25}

对式(25)，令 $\boldsymbol{\xi} = \boldsymbol{\xi}_a$ ，并在两端乘上 $\frac{w_a}{\omega( \boldsymbol{\xi}_a)}$ ，即可得到使用 $f_a$ 离散的方程（ $a=1,2,...,d$ ）：

\frac{\partial f_a}{\partial t} + \xi_a \cdot \nabla{f_a} = -\frac{1}{\tau} ( f_a - f_a^{(eq)} ) + F_a \tag{26-1}

f_a^{(eq)} = \frac{w_a}{\omega (\xi_a)} f^{(eq)} (\xi_a) \tag{26-2}

F_a = \frac{w_a}{\omega (\xi_a)} \mathbf{F} (\xi_a) \tag{26-3}

此时只需将 $f^{(eq)}$ 和 $F$ 进行类似的展开，即可完成对整个系统的离散。

附录A．使用C-E展开方法将离散的LBGK方程还原为Navier-Stokes方程

LBE方程为(A-1.1)式。源项 $\boldsymbol{F}$ 和 $f^{(\mathrm{eq})}$ 为由郭照立、郑楚光和施保昌提出的格式，即(A-1.2)式和(A-1.3)式。[注意，这里的 $f$ 是密度分布函数]

f_i (x+ξ_i Δt,t+Δt)-f_i (x,t)=-\frac{1}{τ} \left[f_i (x,t)-f_i^{(eq)} (x,t) \right]+Δt⋅F_i (x,t) \tag{A-1.1}

F_i=(1 - \frac{Δt}{2τ}) w_i \left[ \frac{\mathbf{\xi_i-u} }{c_s^2 } + \frac{\mathbf{(\xi_i \cdot u) \xi_i} }{c_s^4} \right] ⋅ F \tag{A-1.2}

f_i^{(eq)} (x,t) = \rho w_i \left[ 1 + \frac{\mathbf{ξ_i⋅u}}{c_s^2} +\frac{ \mathbf{uu} : ( \mathbf{ξ_i ξ_i} - c_s^2 \mathbf{I}) }{2c_s^4} \right] \tag{A-1.3}

并且对于 $f_i$ 和 $f_i^{(eq)}$ 有：

\begin{cases} \rho = \sum_i {f_i} \\ \rho \mathbf{u} = \sum_i {\mathbf{\xi_i} f_i} + \frac{\Delta t}{2} F \end{cases}

和

\begin{cases} \rho = \sum_i {f_i^{(eq)}} \\ \rho \mathbf{u} = \sum_i {\mathbf{\xi_i} f_i^{(eq)}} \\ \rho (\mathbf{uu} + c_s^2 \mathbf{I}) = \sum_i {\mathbf{\xi_i \xi_i} f_i^{(eq)}} \end{cases}

将(A-1.1)式在 $(\mathbf{x} , t)$ 处展开，取至 2 阶项得：

\begin{aligned} \Delta t \frac{\partial f_i}{\partial t} + ξ_i \Delta t ⋅ \frac{\partial f_i}{∂x} + &\\ \frac{1}{2} \left[\Delta t^2 \frac{\partial^2 f_i}{\partial t^2} + 2\Delta t(\xi_i \Delta t) ⋅ \frac{\partial^2 f_i}{\partial t∂x} + \Delta t^2 (\xi_i \xi_i): \frac{\partial^2 f_i}{∂x^2} \right] + &\\ \frac{1}{\tau } \left[f_i - f_i^{(eq)} \right] - \Delta t F_i &= 0 \tag{A-2} \end{aligned}

代入多尺度展开等式，并舍去 $\varepsilon^3$ 项和更高阶的项，可得：

\begin{aligned} (\varepsilon \frac{\partial f_i^{(0)}}{\partial t_1} + \varepsilon^2 \frac{\partial f_i^{(1)}}{\partial t_1} + \varepsilon^2 \frac{\partial f_i^{(0)}}{\partial t_2} ) + \xi_i ⋅ (\varepsilon \frac{\partial f_i^{(0)}}{\partial x_1} + \varepsilon^2 \frac{\partial f_i^{(0)}}{\partial x_1}) + &\\ \frac{Δt}{2} \left[ \varepsilon^2 \frac{\partial^2 f_i^{(0)}}{\partial t_1^2} + 2 \xi_i \varepsilon^2 ⋅ \frac{\partial^2 f_i^{(0)}}{\partial t_1 \partial x_1} + (\xi_i \xi_i) : (\varepsilon^2 \frac{\partial^2 f_i^{(0)}}{\partial x^2} ) \right] &\\ + \frac{1}{τΔt} \left[ f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)} -f_i^{(eq)} \right] - \varepsilon F_i^{(1)} + O(\varepsilon^3 ) &=0 \tag{A-3} \end{aligned}

从(A-3)式中提取 $O({\varepsilon}^0)$ 项、 $O({\varepsilon}^1)$ 项和 $O({\varepsilon}^2)$ 项，并令 $D_k = \frac{\partial}{\partial t_k} + \xi_k \cdot \frac{\partial}{\partial \mathbf{x}_k}$ 。(A-3)式对任意 $\varepsilon$ 都成立，则关于各阶的系数均为0，即：

O({\varepsilon}^0) = \frac{1}{\tau \Delta t} (f_i^{(0)} - f_i^{(eq)}) = 0 \tag{A-4.1}

O(\varepsilon^1) = D_1 f_i^{(0)} + \frac{ f_i^{(1)} }{\tau \Delta t} - F_i^{(1)} = 0 \tag{A-4.2}

O(\varepsilon^2) = \frac{\partial f_i^{(0)}}{\partial t_2} - (1 - \frac{1}{2\tau}) D_1 f_i^{(1)} + \frac{ f_i^{(2)} }{\tau \Delta t} + \frac{\Delta t}{2} D_1 F_i^{(1)} = 0 \tag{A-4.3}

根据多尺度展开的变量关系，可得：
$\sum_i {f_i^{(1)}} = \sum_i {f_i^{(2)}} = 0$ ， $\sum_i { \xi_i f_i^{(2)}} = 0$ ， $\sum_i { \xi_i f_i^{(1)}} = - \frac{\Delta t}{2} F_i^{(1)}$

\left\{\begin{matrix} \sum_{i} {F_i} = 0 \\ \sum_{i} {\xi_i F_i} = (1 - \frac{1}{2 \tau}) F \\ \sum_{i} {\xi_i \xi_i F_i} = (1 - \frac{1}{2 \tau}) 2 u F \end{matrix}\right.

，

\left\{\begin{matrix} \sum_{i} {F_i^{(1)}} = 0 \\ \sum_{i} {\xi_i F_i^{(1)}} = (1 - \frac{1}{2 \tau}) F^{(1)} \\ \sum_{i} {\xi_i \xi_i F_i^{(1)}} = (1 - \frac{1}{2 \tau}) 2 u F^{(1)} \end{matrix}\right.

为方便下文分析，记： $\Pi^{(k)} = \sum_i { \xi_i \xi_i f_i^{(k)} }$ ， $\Gamma^{(k)} = \sum_i { \xi_i \xi_i \xi_i f_i^{(k)} }$ 。
对 $O({\varepsilon}^1)$ 项求1阶和2阶速度矩，得：

\frac{\partial \rho}{\partial t_1} + \frac{\partial (\rho u)}{\partial x_1} = 0 \tag{A-5.1}

\frac{\partial (\rho u)}{\partial t_1} + \frac{\partial \Pi^{(0)}}{\partial x_1} = 0 \tag{A-5.2}

\Pi^{(0)} = \rho (\it{\bf{uu}} \mathrm{+} c_s^2 \bf{I}) \tag{A-5.3}

对 $O({\varepsilon}^2)$ 项求1阶和2阶速度矩，得：

\frac{\partial {\rho}}{\partial t_{2}} = 0 \tag{A-6.1}

\frac{\partial (\rho \boldsymbol{u})}{\partial t_{2}} + (1 - \frac{1}{2 \tau}) \frac{\partial \Pi^{(1)}}{\partial \boldsymbol{x}_{1}} + \Delta t (1 - \frac{1}{2 \tau}) \frac{\partial \boldsymbol{F}^{(1)}}{\partial \boldsymbol{x}_{1}} = 0 \tag{A-6.2}

由(A-4.2)式，得： $f_i^{(1)} = \tau \Delta t \{ D_1 f_i^{(0)} - F_i^{(1)}\}$ ，所以

\Pi^{(1)} = \sum_{i} \boldsymbol{\xi}_{i} \boldsymbol{\xi}_{i} f_{i}^{(1)} = -\tau \Delta t [\frac{\partial (\rho \boldsymbol{u} \boldsymbol{u})}{\partial t_{1}} + c_s^2 \frac{\partial (\rho \boldsymbol{I})}{\partial t_{1}} + \frac{\partial \Gamma^{(0)}}{\partial \boldsymbol{x}_{1}} - (1 - \frac{1}{2 \tau}) 2\boldsymbol{uF^{(1)}} ] \tag{A-7}

$\boldsymbol{uF^{(1)}}$ 是Grad的论文中的记号， $(\bf{\it{uF}}^{(1)} )_{\it{lm}} = \mathrm{\frac{1}{2}} (\it{ u_l F_m^{(1)} + u_m F_l^{(1)} } \mathrm{)}$ 。下文为方便分析，取 $\Pi_{lm}^{(1)}$ 分析。

对第1项，有： $\frac{\partial (\rho u_l u_m)}{\partial t_1} = \rho u_l \frac{\partial u_m}{\partial t_1} + u_m \frac{\partial (\rho u_l)}{\partial t_1}$ 。由(A-5.2)式，可得^[4]：

\rho \frac{\partial u_m}{\partial t_1} + u_m \frac{\partial \rho}{\partial t_1} + \frac{\partial}{\partial x_{1n}} \left( \rho u_m u_n + \rho c_s^2 \delta_{mn} \right) = F_m^{(1)} \tag{A-8.1}

\frac{\partial (\rho u_l)}{\partial t_1}+ \frac{\partial}{\partial x_{1n}} \left( \rho u_l u_n + \rho c_s^2 \delta_{ln} \right) = F_l^{(1)} \tag{A-8.2}

并且

\frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} = - u_l u_m \frac{\partial (\rho u_n)}{\partial x_{1n}} + u_l \frac{\partial (\rho u_m u_n)}{\partial x_{1n}} + u_m \frac{\partial (\rho u_l u_n)}{\partial x_{1n}} \tag{A-8.3}

所以

\frac{\partial (\rho u_l u_m)}{\partial x_{1n}} = (u_l F_m^{(1)} + u_m F_l^{(1)}) - \frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} - c_s^2 \left( u_l \frac{\partial \rho}{\partial x_{1m}} + u_m \frac{\partial \rho}{\partial x_{1l}} \right) \tag{A-8.4}

所以第2项的分量为： $c_s^2 \delta_{lm} \frac{\partial \rho}{\partial t_1}$ 。

对第3项的分量，由于

\Gamma_{lmn}^{(0)} = \rho c_s^2 ( u_l \delta_{mn} + u_n \delta_{lm} + u_m \delta_{ln} )

，所以：

\frac{\partial \Gamma_{lmn}^{(0)}}{\partial x_{1n}} = c_s^2 ( \delta_{mn} \frac{\partial \rho u_l}{\partial x_{1n}} + \delta_{lm} \frac{\partial \rho u_n}{\partial x_{1n}} + \delta_{ln} \frac{\partial \rho u_m}{\partial x_{1n}} ) \tag{A-8.5}

所以

\Pi_{lm}^{(1)} = - \tau \Delta t \left[ - \frac{\partial (\rho u_l u_m u_n)}{\partial x_{1n}} + \frac{1}{2 \tau} (u_l F_m^{(1)} + u_m F_l^{(1)}) + \rho c_s^2 ( \frac{\partial u_l}{\partial x_{1m}} + \frac{\partial u_m}{\partial x_{1l}} ) \right] \tag{A-8.6}

略去 $\Pi_{lm}^{(1)}$ 中的速度的3阶项（即 $\frac{\partial}{\partial x_{1n}} \{ \rho u_l u_m u_n \}$ ），有

\Pi^{(1)} = - \Delta t \cdot \boldsymbol{u F}^{(1)} - \tau \rho c_s^2 \Delta t \left[ \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} + \left( \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} \right)^{\mathrm{T}} \right] \tag{A-9}

回代（A-6.2）式，得等价形式为

\frac{\partial (\rho \boldsymbol{u})}{\partial t_2} + c_s^2 \Delta t (\tau - \frac{1}{2}) \frac{\partial}{\partial \boldsymbol{x_1}} \left\{ \rho \left[ \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} + \left( \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{x_1}} \right)^{\mathrm{T}} \right] \right\} = \boldsymbol{0} \tag{A-10}

将(A-10)式和(A-5.2)式联立，得

\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u u)}{\partial x} = - \frac{\partial (\rho c_s^2)}{\partial x} + c_s^2 \Delta t (\tau - \frac{1}{2}) \frac{\partial}{\partial x}\{\rho [\frac{\partial u}{\partial x} + (\frac{\partial u}{\partial x})^T ]\} + F \tag{A-11}

此时令 $p = c_s^2 \rho$ ， $\nu = c_s^2 \Delta t (\tau - \frac{1}{2})$ ，即可还原回Navier-Stokes方程

附录B．Gauss-Hermite积分

对于任意一维函数 $f(\xi)$ ，高斯积分通过寻找积分点 $\xi_k$ 和其对应权重 $w_k$ 实现对 $\int_{a}^{b} \omega(\xi) f(\xi) d \xi$ 的近似，即

\int_{a}^{b} \omega(\xi) f(\xi) d \xi \approx \sum_{k=1}^{n} w_k \xi_k \tag{B-1}

其中 $\omega (\xi)$ 是任意的权重函数。n 点高斯求积的积分点正是对应的第n个正交多项式的根，此时 $\xi_k$ 对应的权重 $w_k$ 为

w_k = \frac{ <P_{n-1}, P_{n-1}> }{ P_{n-1}(\xi_k) P_{n}^{'}(\xi_k) } \tag{B-2}

其中 $P_{n}^{'} = \frac{d P_n}{d x}$ 。(B-2)式的精度为 $2n-1$ 。
所以在一维Gauss-Hermite积分中，取 $P_n = \mathcal{H}^{(n)}$ ， $\mathcal{H}^{(n)}$ 的权函数和表达式见(12)式和(13)式（下文同理）。 n点高斯求积的积分点即为 $\mathcal{H}^{(n)}$ 的根。因为 $\frac{d \mathcal{H}^{(n)}}{d \xi} = \xi \mathcal{H}^{(n)} - \mathcal{H}^{(n+1)} = n \mathcal{H}^{(n-1)}$ ，所以有

w_k = n! / [n \mathcal{H}^{(n-1)} (\xi_k)]^2 \tag{B-3}

对于更高维的情况，可做类似构造^[3:1]。分析如下积分

\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi \tag{B-4}

其中 $\xi$ 是长度为D的向量。 $p(\xi)$ 是D维n次多项式，可写为(B-5)式

p(\xi) = \sum\limits_{n_1 + n_2 + ... +n_D \le n} {a_{n_1 n_2 ... n_D}} \prod\limits_{j=1}\limits^{D}{\xi_{j}^{n_j}} \tag{B-5}

其中 $a_{n_1 n_2 ... n_D}$ 为常数。

考虑(B-5)式中的任意一项 $\prod\limits_{j=1}\limits^{D}{\xi_{j}^{n_j}}$ （由于 $a_{n_1 n_2 ... n_D}$ 为常数所以可暂不处理，在最后一步被归并到常数项中），回代到式(B-4)中，得

\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \prod\limits_{j=1}\limits^{D} \left( \frac{1}{\sqrt{2 \pi}} \int \mathrm{exp}(- \frac{\xi_{j}^2}{2}) \xi_j^{n_j} \mathrm{d}\xi_j \right) = \prod\limits_{j=1}\limits^{D} ( \sum\limits_{k=1}\limits^{n} w_k \xi_k ) \tag{B-6.1}

其中 $w_k$ 和 $\xi_k$ 是一维n次多项式高斯积分的权重和积分点。积分结果被转化为多个一维Gauss-Hermite积分的连乘。注意到(B-6.1)式等价于

\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \sum\limits_{k_1=1}\limits^{n} ... \sum\limits_{k_D=1}\limits^{n} [(w_{k_1} w_{k_2} ... w_{k_D}) (\xi_{k_1}^{n_1} \xi_{k_2}^{n_2} ... \xi_{k_D}^{n_D})] \tag{B-6.2}

如果定义 $\xi_{k_1 ... k_D} = (\xi_{k_1}, \xi_{k_2}, ... , \xi_{k_D})$ 和 $w_{k_1 ... k_D} = w_{k_1} w_{k_2} ... w_{k_D}$ ，那么(B-4)式等价为

\frac{1}{(2 \pi)^{D/2}} \int \mathrm{exp}(- \frac{\xi \cdot \xi}{2}) p(\xi) \mathrm{d} \xi = \sum w_{k_1 ... k_D} p(\xi_{k_1 ... k_D}) \tag{B-7}

此时则将(B-4)式转换为与一维Gauss-Hermite积分类似的形式，即(B-7)式。

参考

GRAD H. Note on N-Dimensional Hermite Polynomials[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 325–330. http://doi.org/10.1002/cpa.3160020402 ↩︎
GRAD H. On the Kinetic Theory of Rarefied Gases[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 331–407. http://doi.org/10.1002/cpa.3160020403 ↩︎
SHAN X, YUAN X-F, CHEN H. Kinetic Theory Representation of Hydrodynamics: A Way beyond the Navier–Stokes Equation[J]. Journal of Fluid Mechanics, 2006, 550(1): 413. http://doi.org/10.1017/S0022112005008153 ↩︎ ↩︎
CAO W. Investigation of the Applicability of the Lattice Boltzmann Method to Free-Surface Hydrodynamic Problems in Marine Engineering[D/OL]. Laboratoire de recherche en Hydrodynamique, Énergétique et Environnement Atmosphérique (LHEEA): École centrale de Nantes, 2019. https://tel.archives-ouvertes.fr/tel-02383174 ↩︎

题 1