工位上换了新电脑，于是需要重新安装 RDPWrap 让机子上的 Windows 11 家庭版能够打开 RDP 远程连接。

1. 安装 RDPWrap

在 RDPWrap 的 Github 网址上下载它的最新 Release 【其实已经很久没更新了】。 解压后运行 install.bat 进行安装。

2. 出现 Listener state: Not-listening 的解决方案。

1）首先在 sebaxakerhtc/rdpwrap.ini 这个仓库里下一份最新的 rdpwrap.ini，一般来说大多数Windows OS 版本号它都有配置。
2）然后将自行下载的 rdpwrap.ini 替换掉原始文件^[1]，原始文件位于 C:\Program Files\RDP Wrapper 文件夹。
3）替换完成后，回到解压 RDPWrap 的文件夹并在此处打开终端，运行 .\RDPWInst.exe -r ^[2]。这之后我的电脑就设置完成了。

工位上换了新电脑，于是需要给新电脑安装 Linux 子系统（WSL 2）用来跑代码。
这里简单记录一下我的安装过程。

1. 安装 WSL 2

这里主要是按照微软官方的教程^[1][2]来做。

先在 powershell 里启用 WSL 和虚拟机功能

1
2

dism.exe /online /enable-feature /featurename:Microsoft-Windows-Subsystem-Linux /all /norestart  
dism.exe /online /enable-feature /featurename:VirtualMachinePlatform /all /norestart  

然后安装适用于 x64 计算机的 WSL2 Linux 内核更新包，并在 powershell 里将 WSL 2 设置为默认。

1

wsl --set-default-version 2

WSL 2 的自定义安装位置参考的是CSDN上给的做法：我下载了Ubuntu 24.04 的 AppxBundle，从包里提取 Ubuntu_2404.0.5.0_x64.Appx ，然后在我想要的位置将它当成 zip 解压。最后运行一下 ubuntu2404.exe 即可。

2. 导入Windows系统的字体

这里主要参照知乎上给的方法^[3]，做了个软链接。

1
2
3

sudo mkdir /usr/share/fonts/win11 
sudo ln -s /mnt/c/Windows/Fonts/* /usr/share/fonts/win11
fc-cache -fv

3. 安装CUDA

Nvidia 的安装教程 和它给的下载链接主要是针对 WSL-Ubuntu 的。

我用的网络下载，在shell里就应该是：

1
2
3
4
5

wget https://developer.download.nvidia.com/compute/cuda/repos/wsl-ubuntu/x86_64/cuda-keyring_1.1-1_all.deb
sudo dpkg -i cuda-keyring_1.1-1_all.deb
sudo apt-get update
sudo apt-get -y install cuda-toolkit-12-6
sudo apt-get -y install cuda

即使是这么安装完，还是会有非常幽默的 nvcc 找不到、 libcuda.so 找不到等问题。我就参照网上的做法^[4][5]在 ~/.bashrc 里补充了这么一段：

1
2
3
4
5
6

# >>> Nvidia CUDA initialize >>>
export CUDA_HOME=/usr/local/cuda
export PATH=${CUDA_HOME}/bin:${PATH}
export LD_LIBRARY_PATH=${CUDA_HOME}/lib64:${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}
export LD_LIBRARY_PATH=/usr/lib/wsl/lib/:${LD_LIBRARY_PATH}
# <<< Nvidia CUDA initialize <<<

搞定。

题 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}}$ ，方向指向原点。

本文同步发表于知乎和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 + \boldsymbol{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$ : 温度；
• $\boldsymbol{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} = -\boldsymbol{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 源项的展开结果。

源项误差的 Chapman-Enskog 分析

这里的推导主要参考 Li 等^[4:2]的文章，并补充了一些自己的笔记。

需要提醒的是，下文的速度记号和前文不同。

源项的各阶速度矩

根据 LBM 格子张量对称性，有前 5 阶表达式：

$$\sum_{\alpha} w_{\alpha} = 1 ,\quad \sum_{\alpha} w_{\alpha} e_{\alpha i} = 0 ,\\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} = c_s^2 \delta_{ij} ,\quad \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} = 0 \\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} e_{\alpha l} = c_s^4 (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{kj}) \\ \sum_{\alpha} w_{\alpha} e_{\alpha i} e_{\alpha j} e_{\alpha k} e_{\alpha l} e_{\alpha m} = 0$$

其中 $e_{\alpha i}$ 表示 D 维向量 $\bold{e}_{\alpha}$ 的第 $i$ 个分量， $w_\alpha$ 为 $\bold{e}_{\alpha}$ 的权重系数。 $\delta_{ij}$ 为克罗内克符号。

所以，对于 EDM 格式的源项 $F_{\alpha,\mathrm{EDM}}$ ，其各阶速度矩为：

$$\sum_{\alpha} F_{\alpha,\mathrm{EDM}} = 0 ,\quad \sum_{\alpha} \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} = \bold{F} ,\quad\\ \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}} = \delta_t \left( (\bold{u}^* \bold{F} + \bold{F}\bold{u}^*) + \delta_t \frac{\bold{FF}}{\rho} \right) \tag{9}$$

技术细节可见附录（A）。

LBM 的常规 Chapman-Enskog 展开

令展开参数 $\epsilon$ 是与克努森数同阶的小量，则将分布函数展开为

$$f_{\alpha} = f_{\alpha}^{(0)} + \epsilon f_{\alpha}^{(1)} + \epsilon^2 f_{\alpha}^{(2)} + ...,$$

源项展开为 $F_{\alpha} = \epsilon F_{1 \alpha}$，时间和空间偏导展开为

$$\frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_1} + \epsilon^2 \frac{\partial}{\partial t_2} ,\quad \frac{\partial}{\partial \bold{x}} = \epsilon \frac{\partial}{\partial \bold{x}_1}.$$

将式(4)中的 $f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t)$ 在 $f_\alpha (\bold{x}, t)$ 处进行 Taylor 展开（至2阶项），并把上述展开式代入，可得 $\epsilon$ 各阶项的方程为：

$$\begin{aligned} O(\epsilon^0):\quad & f_{\alpha}^{(0)} = f_{\alpha}^{(eq)} , &(10-1)\\ O(\epsilon^1):\quad & \mathrm{D}_{1 \alpha} f_{\alpha}^{(0)} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(1)} + \frac{F_{1 \alpha}}{\delta_t} , &(10-2)\\ O(\epsilon^2):\quad & \frac{\partial f_{\alpha}^{(0)}}{\partial t_2} + \mathrm{D}_{1 \alpha} f_{\alpha}^{(1)} + \frac{\delta_t}{2} \mathrm{D}_{1 \alpha}^2 f_{\alpha}^{(0)} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(2)} . &(10-3)\\ \end{aligned}$$

其中 $\displaystyle \mathrm{D}_{1 \alpha} = \frac{\partial}{\partial t_1} + \bold{e}_{\alpha} \cdot \frac{\partial}{\partial \bold{x}_1}$。将式(10-2)代入式(10-3)中，式(10-3)可被化简为式(10-4)：

$$\begin{aligned} O(\epsilon^2):\quad & \frac{\partial f_{\alpha}^{(0)}}{\partial t_2} + \mathrm{D}_{1 \alpha} \left(1 - \frac{1}{2 \tau}\right) f_{\alpha}^{(1)} + \frac{1}{2} \mathrm{D}_{1 \alpha} F_{1 \alpha} = -\frac{1}{\tau \delta_t} f_{\alpha}^{(2)} . &(10-4) \end{aligned}$$

（1）考虑对式(10-3)和式(10-4)计算零阶矩，得：

$$\frac{\partial \rho}{\partial t_1} + \nabla_1 \cdot (\rho \bold{u}^*) = 0 \\ \frac{\partial \rho}{\partial t_2} + \frac{1}{2} \nabla_1 \cdot (\delta_t \bold{F}_1) = 0 \\$$

此处 $\bold{u}^* = (\sum_{i} \bold{e}_\alpha f_\alpha) / (\sum_{i} f_\alpha)$ 并非真实流场速度，而 $\hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho)$ 为真实流场速度（也就是式(6)给出的形式）。

基于前面假设的偏微分关系式，我们可以得到 $\displaystyle \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}}) = 0$ 。

（2）同理，对式(10-3)和式(10-4)计算一阶矩，得：

$$\frac{\partial \rho \bold{u}^*}{\partial t_1} + \nabla_1 \cdot (\rho \bold{u}^* \bold{u}^*) = -\nabla_1 p + \bold{F}_1 \\ \frac{\partial \rho \bold{u}^*}{\partial t_2} + \nabla_1 \cdot \left[\left(1 - \frac{1}{2 \tau}\right) \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(1)} \right] + \frac{\delta_t}{2} \frac{\partial \bold{F}_1}{\partial t_1} + \frac{1}{2} \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right) = 0$$

其中 $p = \rho c_s^2$ 。

根据式(10-2)，有：

$$\begin{aligned} \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(1)} &= -\tau \left[ \delta_t \frac{\partial}{\partial t_1}\left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) \right. +\\ &\quad \left. \delta_t \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) - \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right] \end{aligned}$$

记

$$\begin{aligned} \Pi^{(1)} &= -\delta_t \left( \tau - \frac{1}{2} \right) \left[ \frac{\partial}{\partial t_1}\left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) + \nabla_1 \cdot \left( \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \bold{e}_\alpha f_{\alpha}^{(0)} \right) \right] \\ &= -\delta_t \left( \tau - \frac{1}{2} \right) \{\rho c_s^2 [\nabla_1 \bold{u}^* + (\nabla_1 \bold{u}^*)^T] + \bold{u}^*\bold{F}_1 + \bold{F}_1 \bold{u}^* \} \end{aligned}$$

则式(10-4)的一阶矩表示为： $\displaystyle \frac{\partial \rho \bold{u}^*}{\partial t_2} + \nabla_1 \Pi^{(1)} + \frac{\delta_t}{2} \frac{\partial \bold{F}_1}{\partial t_1} + \nabla_1 \cdot \left( \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{1 \alpha} \right) = 0$ 。

那么同样根据假设的微分关系，即可还原出一阶矩对应的方程。这里直接写出使用真实流场速度 $\hat{\bold{u}} = \bold{u}^* + \bold{F} \delta_t / (2 \rho)$ 表示的形式：

$$\begin{aligned} \frac{\partial (\rho \hat{\bold{u}})}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}} \hat{\bold{u}}) &= -\nabla p + \bold{F} + \nabla \cdot \{ \rho \nu [\nabla_1 \bold{u}^* + (\nabla_1 \bold{u}^*)^T \} \\ &\quad - \nabla \cdot \{ \rho \nu [\nabla_1 \bold{a} + (\nabla_1 \bold{a})^T] \} + \frac{\delta_t}{2} \epsilon^2 \frac{\partial \bold{F}}{\partial t_2} \\ &\quad+ \nabla \cdot \left[ \tau \delta_t (\bold{u^* F} + \bold{F u^*}) + \delta_t^2 \frac{\bold{FF}}{4 \rho} - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha} \right] \end{aligned}$$

其中 $\bold{a} = \delta_t \bold{F} / 2 \rho$ ， $\nu = (\tau - \frac{1}{2}) c_s^2 \delta_t$ 为流体运动粘度。 通过对比N-S方程，可见其误差项为：

$$\text{Err} = - \nabla \cdot \{ \rho \nu [\nabla_1 \bold{a} + (\nabla_1 \bold{a})^T] \} + \frac{\delta_t}{2} \epsilon^2 \frac{\partial \bold{F}}{\partial t_2} + \nabla \cdot \left[ \tau \delta_t (\bold{u^* F} + \bold{F u^*}) + \delta_t^2 \frac{\bold{FF}}{4 \rho} - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha} \right]$$

但上述的展开方式和Wagner等^[5]使用Taylor展开的结果并不一致。 Li 等^[4:3] 指出这种错误的来源是：上文将 $f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \bold{u}^*)$，从而得到了 $\sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = 0$。 然而实际上 $f_{\alpha}^{(0)} = f_{\alpha}^{(eq)}(\rho, \hat{\bold{u}})$ ，从而有 $\sum_{\alpha} \bold{e}_{\alpha} f_{\alpha}^{(1)} = -\delta_t \bold{F}_1 / 2$ 。

修正后的 Chapman-Enskog 展开

因此，需要将原本的 LBE 等效地修改为：

$$f_\alpha (\bold{x} + \bold{e}_\alpha \delta_t, t+\delta_t) - f_\alpha (\bold{x}, t) = \hat{F}_{\alpha} -\frac{1}{\tau} \left[ f_\alpha (\bold{x}, t) - f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) \right] \tag{11}$$

其中 $\hat{F}_{\alpha} = F_{\alpha} - \frac{1}{\tau} (f_\alpha^{(eq)} (\rho, \hat{\bold{u}}) - f_\alpha^{(eq)} (\rho, \bold{u}^*))$ 为基于原始源项 $F_{\alpha}$ 修改而来的等效源项。并且需要强调的是，式(11)中的平衡态采用真实流场速度 $\hat{\bold{u}}$ 计算。 在式(11)这一等效公式上，执行与上一节一致的 Chapman-Enskog 展开后，可得到还原的N-S方程为：

$$\begin{aligned} \frac{\partial (\rho \hat{\bold{u}})}{\partial t} + \nabla \cdot (\rho \hat{\bold{u}} \hat{\bold{u}}) &= -\nabla p + \bold{F} + \nabla \cdot \{ \rho \nu [\nabla_1 \hat{\bold{u}} + (\nabla_1 \hat{\bold{u}})^T \} \\ &\quad + \nabla \cdot \left[ \left(\tau - \frac{1}{2}\right) \delta_t (\hat{\bold{u}}\bold{F} + \bold{F}\hat{\bold{u}}) - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \hat{F}_{\alpha} \right] \end{aligned} \tag{12}$$

式(12)的误差项为：

$$\text{Err} = \nabla \cdot \left[ \left(\tau - \frac{1}{2}\right) \delta_t (\hat{\bold{u}}\bold{F} + \bold{F}\hat{\bold{u}}) - \tau \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \hat{F}_{\alpha} \right] \tag{13}$$

将 EDM 的源项代入式(13)中，得到： $\text{Err}=-\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right)$

附录

（A）源项各阶速度矩的推导

这里仅以 $\displaystyle \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha F_{\alpha,\mathrm{EDM}}$ 的计算举例来简单说明。

$$\begin{aligned} \sum_{\alpha} \bold{e}_{\alpha} \bold{e}_{\alpha} F_{\alpha,\mathrm{EDM}} &= \sum_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot \delta_t w_\alpha \left[ \frac{\bold{e}_\alpha \cdot \bold{F}}{c_s^2} + \right.\\ &\quad\left. \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] \\ &= \frac{\delta_t}{c_s^2} \sum_{\alpha} \left[ \bold{e}_\alpha \bold{e}_\alpha (w_\alpha \bold{e}_\alpha \cdot \bold{F}) \right] + \\ &\quad \frac{\delta_t}{2 c_s^4} \sum_{\alpha} w_{\alpha} \bold{e}_\alpha \bold{e}_\alpha \cdot [\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}])] \end{aligned}$$