神经网络求解Boltzmann-BGK方程及其在微流中的应用

A Neural Network Method for the Boltzmann-BGK Equation with Applications in Microflows

摘要: 本文提出一种基于神经网络的BGK方程求解方法，特别关注在微流问题中的应用。首先，通过引入灵活辅助分布函数构造BGK方程的降维模型，从而有效降低方程维度。其次，设计全连接神经网络架构高效逼近降维分布函数，以避免时空离散化。接着针对微流问题中复杂的Maxwell边界条件，提出特殊设计的损失函数进行处理。此外，利用多尺度输入策略和Maxwellian分裂技术以提升逼近效率。最后，通过对一维Couette流和二维矩形风管流两个经典问题进行数值实验，验证了该方法的有效性。

Abstract: We consider the neural representation to solve the BGK equation, especially focusing on the application in microscopic flow problems. Firstly, a new dimension reduction model of the BGK equation with the flexible auxiliary distribution functions is deduced to reduce the problem dimension. Then, a fully connected neural network is utilized to approximate the dimension-reduced distribution with extremely high efficiency and to avoid discretization in space and time. A specially designed loss function is employed to deal with the Maxwell boundary conditions in microflow problems. Moreover, strategies such as multi-scale input and Maxwellian splitting are applied to further enhance the approximation efficiency. Finally, two classical numerical experiments, including one-dimensional Couette flow and two-dimensional duct flow, are studied to demonstrate the effectiveness of this neural representation method.

1. 引言

在航空航天和微电子机械系统领域，动理学理论的模拟已经引起广泛关注，特别是稀薄气体动理学[1]的应用。在稀薄情况下，传统的连续流体模型，如Euler方程和Navier-Stokes方程，已不再适用。此时，Boltzmann方程为研究稀薄气体动理学提供了一个基本思路。然而，Boltzmann方程的高维特征，以及复杂的二次碰撞项，为高效且准确的数值模拟带来了巨大的挑战。

对于求解Boltzmann方程的传统数值方法，通常可以划分为确定性方法和随机方法。随机方法主要包括直接模拟蒙特卡洛法[2] (DSMC)，而确定性方法则涵盖了离散速度方法[3] (DVM)、谱方法[4] [5]和矩方法[6]。近年来，随着计算能力的提升，利用神经网络求解Boltzmann方程的数值方法得到了快速发展，大致可以分为三类。第一类是基于神经网络构建代理模型，以近似碰撞算子[7]。第二类则是通过神经网络学习闭包模型，以简化Boltzmann方程的求解[8]。第三类是在物理信息神经网络[9] (PINNs)框架下，将偏微分方程及其初始和边界条件直接纳入神经网络的损失函数中，以实现对Boltzmann方程的数值求解。

本文提出了一种降维神经网络表示方法(简称为DRNR)求解Boltzmann-BGK方程。首先，基于已有研究[10] [11]，通过引入与微观速度相关的灵活辅助分布函数，构建了一种新的降维模型。当分布函数在微观速度空间的某些方向上呈现平面对称性时，该模型能够使BGK方程的微观速度空间维度与物理空间维度保持一致。在此框架下，进一步推导出了适用于微流问题的简化Maxwell边界条件[12]。尽管该降维模型增加了未知分布函数的数量，但显著减少了自变量的维数，从而大幅降低了总体计算量。

接下来，在离散速度法的框架下推导出BGK方程的半离散降维模型，并利用神经网络对降维分布函数进行建模。此外，该方法采用多尺度输入和Maxwellian分裂法[13]，利用全连接神经网络对降维后的分布函数进行参数化，并且针对微流问题，特别设计了针对简化Maxwell边界条件的损失函数，并在训练过程中引入了自适应加权策略。最后，通过一维Couette流和二维矩形风管流的数值实验，验证了该方法的有效性。

本文其余部分的结构如下：第2章主要介绍Boltzmann-BGK方程的基本性质、Maxwell边界条件，以及基于灵活辅助分布函数的降维方法。第3章详细阐述了DRNR方法，包括半离散DVM模型、网络结构和损失函数。第4章展示了相关的数值实验结果。最后，第5章给出了本文结论。

2. Boltzmann-BGK方程

2.1. 预备知识

Boltzmann-BGK方程以统计学原理为基础，用于描述粒子的运动，具体形式[14]如下

$\frac{\partial f}{\partial t}+v\cdot {\nabla }_{x}f=\frac{1}{Kn}\left(ℳ\left[f\right]-f\right),t\ge 0,x\in {ℝ}^3,v\in {ℝ}^3,$ (1)

其中$f\left(t,x,v\right)$ 表示粒子的密度函数，随时间$t$ 、空间$x$ 以及微观速度$v$ 变化，Kn是Knudsen数，描述气体的稀薄程度，$ℳ\left[f\right]$ 是局部Maxwellian分布函数，取决于局部气体密度$\rho$ 、宏观速度$u$ 和温度$T$ ，具体形式为

$ℳ\left[f\right]=\frac{\rho }{{\sqrt{2\pi T}}^3}\exp \left(-\frac{{\left|v-u\right|}^2}{2T}\right).$ (2)

而密度$\rho$ 、宏观速度$u={\left(u_1,u_2,u_3\right)}^{\text{T}}$ 和温度$T$ 等宏观变量可以由分布函数得到

$\begin{array}{l}密度:\rho \left(t,x\right)={\displaystyle {\int }_{{ℝ}^3}f}\left(t,x,v\right)\text{d}v,\\ 动量:m\left(t,x\right)\triangleq \rho \left(t,x\right)u\left(t,x\right)={\displaystyle {\int }_{{ℝ}^3}v}f\left(t,x,v\right)\text{d}v,\\ 能量:E\left(t,x\right)\triangleq \frac{3}{2}\rho T+\frac{1}{2}\rho {\left|u\right|}^2=\frac{1}{2}{\displaystyle {\int }_{{ℝ}^3}{\left|v\right|}^{2}f\left(t,x,v\right)\text{d}v}.\end{array}$ (3)

Boltzmann方程的应用涉及气体与实心壁之间的相互作用，通常采用Maxwell边界条件[12]来描述相关过程。在时间$t$ 和边界点$x$ 处，如果壁面处于静止状态，对于满足$n\cdot v>0$ 的速度$v$ ，边界分布函数$f^b\left(t,x,v\right)$ 可由固壁边界条件确定，其中$n$ 是指向气体内部的壁面法向量。当速度不满足条件时，边界分布函数则受内部分布函数的影响。假设固壁速度和温度分别是$u^W={\left(u_1^W,u_2^W,u_3^W\right)}^{\text{T}}$ 和$T^W$ ，则边界条件为

$f^b\left(t,x,v\right)=\{\begin{array}{l}{f}^W,c^W\cdot n>0,\\ f_N\left(t,x,v\right),c^W\cdot n<0,\end{array}$ (4)

其中$c^W=v-u^W$ ，$f_N\left(t,x,v\right)$ 是内部分布函数，$f^W$ 是由固壁速度$u^W$ 和温度$T^W$ 确定的Maxwellian分布函数

$f^W\left({\rho }^W,u^W,T^W\right)=\frac{{\rho }^W}{{\sqrt{2\pi T^W}}^3}\exp \left(-\frac{{\left|v-u^W\right|}^2}{2T^W}\right).$ (5)

其中，

${\rho }^W=-\frac{{\displaystyle {\int }_{c^W\cdot n<0}\left(c^W\cdot n\right)f_N\left(t,x,v\right)dv}}{{\displaystyle {\int }_{c^W\cdot n>0}\left(c^W\cdot n\right)f^W\left(1,u^W,T^W\right)dv}}.$ (6)

数值模拟中，Maxwell边界条件(4)和(6)被应用于微流问题中，且3.3节中将详细介绍对该边界条件设计的损失函数。

2.2. 降维方法

在某些数值实验中，我们会遇到物理空间的维数$D_x$ 低于微观速度维数$D_v$ 的情况，这时通常采用降维模型以降低计算成本。文献[10] [11]提出引入辅助分布函数对BGK型碰撞模型的分布函数进行降维，但是，在计算Couette流问题时，如果在物理空间中存在不同方向的宏观速度，那么微观速度中的辅助分布函数的维数要大于物理空间的$D_x$ 。为解决该问题，我们对降维方法进行改进，采用了灵活的辅助分布函数，在此框架下，微观速度空间中的辅助分布函数的维数可降低至与空间维数相同。

准确地说，考虑$D_x<D_v$ 时的情形，假设分布函数$f\left(t,x,v\right)$ 在物理空间$D_x+1$ 维到$D_v$ 维中是均匀一致的

$f\left(t,x,v\right)=f\left(t,\tilde{x},v\right),\tilde{x}={\left(x_1,\cdots ,x_{D_x}\right)}^{\text{T}}\in {ℝ}^{D_x},\frac{\partial f\left(t,x,v\right)}{\partial x_d}=0,d\in \mathcal{D},$ (7)

在微观速度空间$D_{v^{\ast }}$ 维到$D_v$ 维是平面对称的。可以推断出宏观速度$u$ 满足

$u_k=0,k\in \widehat{\mathcal{D}},$ (8)

其中，

$\mathcal{D}=\left\{D_x+1,\cdots ,D_v\right\},\widehat{\mathcal{D}}=\left\{D_{v^{\ast }}+1,\cdots ,D_v\right\},D_x\le D_{v^{\ast }}\le D_v=3.$ (9)

接下来，定义严格降维分布函数为

$g\left(t,\tilde{x},\tilde{v}\right)={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}f}\left(t,\tilde{x},v\right)d\widehat{v},h\left(t,\tilde{x},\tilde{v}\right)={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}\frac{{\left|\widehat{v}\right|}^2}{2}f\left(t,\tilde{x},v\right)d\widehat{v}},$ (10)

灵活降维分布函数为

$s_j\left(t,\tilde{x},\tilde{v}\right)={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}{v}_j}f\left(t,\tilde{x},v\right)d\widehat{v},j=D_x+1,\cdots ,D_{v^{\ast }},$ (11)

其中，

$\tilde{v}={\left(v_1,\cdots ,v_{D_x}\right)}^{\text{T}}\in {ℝ}^{D_x},\widehat{v}={\left(v_{D_x+1},\cdots ,v_{D_v}\right)}^{\text{T}}\in {ℝ}^{D_v-D_x},$ (12)

灵活降维分布函数(11)的数量是$D_{v^{\ast }}-D_x$ 。简单起见，除非另有说明，本文后续在使用$s_j$ 时均暗指 $j=D_x+1,\cdots ,D_{v^{\ast }}$ 。

在改进的降维方法中，宏观变量与降维分布函数之间的关系为

$\begin{array}{l}\rho ={\displaystyle {\int }_{{ℝ}^{D_x}}g}\left(t,\tilde{x},\tilde{v}\right)d\tilde{v},\\ \rho u_i={\displaystyle {\int }_{{ℝ}^{D_x}}{v}_i}g\left(t,\tilde{x},\tilde{v}\right)d\tilde{v},i=1,\cdots ,D_x,\\ \rho u_j={\displaystyle {\int }_{{ℝ}^{D_x}}{s}_j}\left(t,\tilde{x},\tilde{v}\right)d\tilde{v},j=D_x+1,\cdots ,D_{v^{\ast }},\\ E=\frac{3}{2}\rho T+\frac{1}{2}\rho {\left|u\right|}^2={\displaystyle {\int }_{{ℝ}^{D_x}}\left(\frac{{\left|\tilde{v}\right|}^2}{2}g\left(t,\tilde{x},\tilde{v}\right)+h\left(t,\tilde{x},\tilde{v}\right)\right)d\tilde{v}}.\end{array}$ (13)

由于维数降低，无法得到部分高阶矩，因此，我们仅计算部分应力张量${\sigma }_{ij}$ 和热通量$q_i$

$\begin{array}{l}{\sigma }_{ij}=\{\begin{array}{ll}{\displaystyle \underset{{ℝ}^{D_x}}{\int }{v}_i}{v}_{j}g\left(t,\tilde{x},\tilde{v}\right)d\tilde{v}-\rho u_i{u}_j-{\delta }_{ij}\rho T,\hfill & i,j\le D_x,\hfill \\ {\displaystyle \underset{{ℝ}^{D_x}}{\int }{v}_i}{s}_j\left(t,x,v\right)d\tilde{v}-\rho u_i{u}_j,\hfill & i\le D_x<j\le D_{v^{\ast }},\hfill \end{array}\\ q_i={\displaystyle \underset{{ℝ}^{D_x}}{\int }\left(v_i-u_i\right)}\left(\frac{1}{2}{\left|\tilde{c}\right|}^{2}g\left(t,\tilde{x},\tilde{v}\right)+h\left(t,\tilde{x},\tilde{v}\right)\right)d\tilde{v}\\ +{\displaystyle \sum _{j=D_x+1}^{D_{v^{\ast }}}\left(u_j^2{\displaystyle {\int }_{{ℝ}^{D_x}}\left(v_i-u_i\right)g\left(t,\tilde{x},\tilde{v}\right)d\tilde{v}}-2u_j{\displaystyle {\int }_{{ℝ}^{D_x}}\left(v_i-u_i\right)s_j\left(t,\tilde{x},\tilde{v}\right)d\tilde{v}}\right)},i\le D_x,\end{array}$ (14)

其中$\tilde{c}=\tilde{v}-\tilde{u}$ ，$\tilde{u}={\left(u_1,\cdots ,u_{D_x}\right)}^{\text{T}}$。

最后，通过对(1)式积分计算，可得到降维Boltzmann方程

$\begin{array}{l}\frac{\partial g}{\partial t}+\tilde{v}\cdot {\nabla }_{\tilde{x}}g=\frac{1}{Kn}\left(g^M-g\right),\\ \frac{\partial h}{\partial t}+\tilde{v}\cdot {\nabla }_{\tilde{x}}h=\frac{1}{Kn}\left(h^M-h\right),\end{array}$

$\frac{\partial s_j}{\partial t}+\tilde{v}\cdot {\nabla }_{\tilde{x}}{s}_j=\frac{1}{Kn}\left(s_j^M-s_j\right),$ (15)

其中$g^M,h^M$ 和$s_j^M$ 与Maxwellian分布函数有关，即

$g^M={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}ℳ}\left[f\right]d\widehat{v},h^M={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}\frac{{\left|\widehat{v}\right|}^2}{2}ℳ\left[f\right]d\widehat{v}},s_j^M={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}{v}_j}ℳ\left[f\right]d\widehat{v}.$ (16)

对于该降维方程，可以应用2.1节中引入的Maxwell边界条件，则由降维分布函数(10)和(11)以及边界条件(4)，可得新的边界条件为

$F^b\left(t,\tilde{x},\tilde{v}\right)=\{\begin{array}{l}{F}^W,{\tilde{c}}^W\cdot \tilde{n}>0,\\ F_N\left(t,\tilde{x},\tilde{v}\right),{\tilde{c}}^W\cdot \tilde{n}<0,\end{array}$ (17)

其中$F=g,h,s_j$ ，${\tilde{c}}^W=\tilde{v}-{\tilde{u}}^W$，${\tilde{u}}^W={\left(u_1^W,\cdots ,u_{D_x}^W\right)}^{\text{T}}$，$\tilde{n}$ 是降维后指向气体内部的壁面法向量，$F_N\left(t,\tilde{x},\tilde{v}\right)$ 是内部分布函数，由(5)式可得$F^W=g^W,h^W,s_j^W$ 的表达形式为

$g^W={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}{f}^W}d\widehat{v},h^W={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}\frac{{\left|\widehat{v}\right|}^2}{2}}{f}^{W}d\widehat{v},s_j^W={\displaystyle {\int }_{{ℝ}^{D_v-D_x}}{v}_j}{f}^{W}d\widehat{v}.$ (18)

分布函数$f^W$ 可由${\rho }^W,u^W,T^W$ 显式表示，${\rho }^W$ 可由下式计算得到

${\rho }^W=-\frac{{\displaystyle {\int }_{{\tilde{c}}^W\cdot \tilde{n}<0}\left({\tilde{c}}^W\cdot \tilde{n}\right)g_N\left(t,\tilde{x},\tilde{v}\right)d\tilde{v}}}{{\displaystyle {\int }_{{\tilde{c}}^W\cdot \tilde{n}>0}\left({\tilde{c}}^W\cdot \tilde{n}\right)g^W\left(1,u^W,T^W\right)d\tilde{v}}}.$ (19)

在接下来章节中，将使用降维方程(15)~(17)作为控制方程，用以设计神经网络表示方法。

3. 降维神经网络表示方法

3.1. 基于DVM的空间离散

根据已有研究[13]，首先对降维方程(15)在微观速度空间中进行离散，将得到的半离散系统作为控制方程，以设计神经网络结构。在此框架下，输入维度显著降低，从而简化神经网络的规模和复杂性。

首先，为了在接下来的章节中清晰呈现，将(7)和(12)中的$\tilde{x}$ 和$\tilde{v}$ 缩写为$x$ 和$v$ 。假设微观速度空间中的全离散点集为

$V=\left[v_1,v_2,\cdots ,v_{N_v}\right]\in {ℝ}^{D_x\times N_v},$ (20)

其中$N_v$ 是离散点总数，$v_k={\left(v_{1k},\cdots ,v_{D_{x}k}\right)}^{\text{T}},1\le k\le N_v$ ，因此降维BGK方程(15)的半离散形式如下

$\begin{array}{l}\frac{\partial g}{\partial t}+V\cdot {\nabla }_{x}g=\frac{1}{Kn}\left(g^M-g\right),\\ \frac{\partial h}{\partial t}+V\cdot {\nabla }_{x}h=\frac{1}{Kn}\left(h^M-h\right),\\ \frac{\partial s_j}{\partial t}+V\cdot {\nabla }_x{s}_j=\frac{1}{Kn}\left(s_j^M-s_j\right),\end{array}$ (21)

其中，$F=g,g^M,h,h^M,s_j,s_j^M$ 是离散的降维分布函数

$F=\left[F_1,F_2,\cdots ,F_{N_v}\right]\in {ℝ}^{1\times N_v},F_k=g_k,g_k^M,h_k,h_k^M,s_{jk},s_{jk}^M,k=1,\cdots ,N_v.$ (22)

在接下来章节中，将使用降维方程的半离散形式(21)作为控制方程，用以设计DRNR中的神经网络结构和损失函数。

3.2. 网络结构

本节将详细介绍用于表示降维分布函数的全连接神经网络结构。通常来说，一个全连接神经网络或前馈网络包括一个输入层、若干隐藏层和一个输出层。我们使用深度神经网络求解降维Boltzmann方程，并结合多尺度输入和Maxwellian分裂技术以提高近似效率。由于正弦函数是周期函数，具有全局性质，适合处理周期性或高频信息，能够对输入的不同区域学习相同的函数，从而更好地处理全局特征[15]，本文选择$\sigma \left(x\right)=\sin \left(x\right)$ 作为激活函数。

具体地，所有神经网络的输入变量为

${\widehat{x}}^{\text{multi}}=\left(c_1\widehat{x},c_2\widehat{x},\cdots ,c_K\widehat{x}\right)\in {ℝ}^{K\left(D_x+1\right)},$ (23)

其中，$c_i,i=1,\cdots ,K$ 是与问题相关的常数，$\widehat{x}=\left(t,x\right)\in {ℝ}^{D_x+1}$ 。为了进一步提高近似效率，我们引入了Maxwellian分裂法，通过使用该方法，可将离散的降维分布函数(22)分解为

$F\left(t,x\right)=F^{\text{eq}}\left(t,x\right)+C{F}^{\text{neq}}\left(t,x\right),F=g,h,s_j,$ (24)

其中$C$ 是可调节的参数，根据宏微观分解理论，一般选取与Kn相同的数值，$F^{\text{eq}}\left(t,x\right)\in {ℝ}^{1\times N_v}$ 是Maxwellian分布函数

$F^{\text{eq}}\left(t,x\right)=\left(F_1^{\text{eq}}\left(t,x\right),F_2^{\text{eq}}\left(t,x\right),\cdots ,F_{N_v}^{\text{eq}}\left(t,x\right)\right),$ (25)

其中

$F_k^{\text{eq}}\left(t,x\right)=\frac{{\rho }_F^{\ast }\left(t,x\right)}{{\sqrt{2\pi T_F^{\ast }\left(t,x\right)}}^{D_x}}\exp \left(-\frac{{\left|v_k-u_F^{\ast }\left(t,x\right)\right|}^2}{2T_F^{\ast }\left(t,x\right)}\right),k=1,\cdots ,N_v.$ (26)

$F^{\text{neq}}\left(t,x\right)$ 是非平衡态分布函数，假设具有以下形式

$F^{\text{neq}}\left(t,x\right)=\left(F_1^{\text{neq}}\left(t,x\right),\cdots ,F_{N_v}^{\text{neq}}\left(t,x\right)\right),F_k^{\text{neq}}\left(t,x\right)=F_k^{\text{eq}}\left(t,x\right)F_k^{\ast }\left(t,x\right),$ (27)

其中，

$F^{\ast }\left(t,x\right)=\left(F_1^{\ast }\left(t,x\right),\cdots ,F_{N_v}^{\ast }\left(t,x\right)\right)\in {ℝ}^{1\times N_v}.$ (28)

实验表明，该假设可以确保$F^{\text{eq}}\left(t,x\right)$ 和$F^{\ast }\left(t,x\right)$ 的数量级相似，从而减少训练过程中遇到的困难。因此，我们将伪宏观变量$\left({\rho }_F^{\ast },u_F^{\ast },T_F^{\ast }\right)$ 和非平衡态离散分布函数$F^{\ast }$ 作为输出变量，利用两个独立的神经网络进行训练，如(29)所示

Figure 1. Schematic of DRNR

图1. DRNR示意图

$\begin{array}{l}\left({\rho }_F^{\ast },u_F^{\ast },T_F^{\ast }\right)\triangleq {\text{NN}}_{\theta }^{\text{eq}}\left(t,x\right)\in {ℝ}^{D_x+2},\\ F^{\ast }\left(t,x\right)\triangleq {\text{NN}}_{\theta }^{\text{neq}}\left(t,x\right)\in {ℝ}^{D_v},F=g,h,s_j.\end{array}$ (29)

需要强调的是，对于每一个$F=g,h,s_j$ ，均对应两个独立的神经网络和相应的伪宏观变量，因此，网络总数为$2\left(2+D_{v^{\ast }}-D_x\right)$ 。神经网络的结构如图1所示，阴影区域中的具体结构如图1的右图所示。

在神经网络建立之后，需要设置合适的损失函数将神经网络与PDE联系起来，具体细节将在下一节中进行介绍。

3.3. 损失函数

本节介绍了一种针对Maxwell边界条件特别设计的损失函数。一般来说，损失函数包括方程的残差$L_{\text{Eq}}$ ，以及作为惩罚项的初始条件约束$L_{\text{IC}}$ 和边界条件约束$L_{\text{BC}}$ ，具体形式如下

$L_{\text{loss}}=L_{\text{Eq}}+{\lambda }_1{L}_{\text{IC}}+{\lambda }_2{L}_{\text{BC}},$ (30)

其中，${\lambda }_1$ 和${\lambda }_2$ 是与问题相关的惩罚项权重。

$L_{\text{Eq}},L_{\text{IC}},L_{\text{BC}}$ 均由以下三项构成

$L_m=L_m^g+L_m^h+{\displaystyle \sum _{j=D_x+1}^{D_{v^{*}}}{L}_m^{s_j}},$ (31)

其中$m=\text{Eq},\text{IC},\text{BC}$ 。本文的边界条件设置为Maxwell边界条件，则(18)式中的$F^b\left(t_i,x_i,v_k\right)$ 取为

$F^b\left(t_i,x_i,v_k\right)=F^W\left(t_i,x_i,v_k\right),F^W=g^W,h^W,s_j^W,$ (32)

其中壁面速度$u^W$ 和温度$T^W$ 由边界条件给定，内部分布函数可由神经网络的输出表示，根据(6)式可计算得到${\rho }^W$ ，如下，

${\rho }^W\left(t_i,x_i\right)=-\frac{{\displaystyle \sum _{\begin{array}{c}k=1\\ v_k\cdot n<0\end{array}}^{N_v}\left(v_k\cdot n\right)g_{k,\theta }^{\text{NN}}\left(t_i,x_i\right)w_k}}{{\displaystyle \sum _{\begin{array}{c}k=1\\ v_k\cdot n>0\end{array}}^{N_v}\left(v_k\cdot n\right)g_k^W\left(1,u^W\left(t_i,x_i\right),T^W\left(t_i,x_i\right)\right)w_k}}.$ (33)

注1. 稳态问题中，神经网络的输入仅仅是物理空间$x$ ，大多数稳态问题满足质量守恒，因此，我们将神经网络$g_{\theta }^{\text{NN}}$ 的输出预先调整为

$g_{\theta }^{\text{NN},\text{steady}}\left(x\right)=\frac{g_{\theta }^{\text{NN}}\left(x\right)}{{\rho }_{\text{ave}}},$ (34)

其中，

${\rho }_{\text{ave}}=\frac{{\displaystyle \sum _{i=1}^{N_{\text{Eq}}}\langle g_{\theta }^{\text{NN}}\left(x_i\right)\rangle }}{N_{\text{Eq}}},\langle g\rangle ={\displaystyle \sum _{k=1}^{N_v}{g}_k}{\omega }_k.$ (35)

$x_i,i=1,\cdots ,N_{\text{Eq}}$ 是$\Omega$ 的随机点。

由(30)式可知，损失函数是初值、边界条件以及方程残差的加权组合，然而使用DVM进行空间离散时，每个速度离散点被赋予相同的权重，实际上对于宏观量的计算，当相对速度较小时，分布函数更为重要。为此，使用下界约束不确定性加权方法[13]进行微观速度点$v_k$ 的权重计算，可得

$\begin{array}{l}{L}_{\text{Eq}}^F=\frac{1}{N_{\text{Eq}}}{\displaystyle \sum _{i=1}^{N_{\text{Eq}}}{\displaystyle \sum _{k=1}^{N_v}{\left(\frac{1}{{\left(w_{\text{Eq}}^F\right)}_k+\epsilon }{R}_k^F\left(t_i,x_i\right)+\log \left(1+{\left(w_{\text{Eq}}^F\right)}_k\right)\right)}^2}},\\ L_{\text{IC}}^F=\frac{1}{N_{\text{IC}}}{\displaystyle \sum _{i=1}^{N_{\text{IC}}}{\displaystyle \sum _{k=1}^{N_v}{\left(\frac{1}{{\left(w_{\text{IC}}^F\right)}_k+\epsilon }\left(F_{k,\theta }^{\text{NN}}\left(0,x_i\right)-F_k^0\left(x_i\right)\right)+\log \left(1+{\left(w_{\text{IC}}^F\right)}_k\right)\right)}^2}},\\ L_{\text{BC}}^F=\frac{1}{N_{\text{BC}}}{\displaystyle \sum _{i=1}^{N_{\text{BC}}}{\displaystyle \sum _{\begin{array}{c}k=1\\ v_k\cdot n>0\end{array}}^{N_v}{\left(\frac{1}{{\left(w_{\text{BC}}^F\right)}_k+\epsilon }\left(F_{k,\theta }^{\text{NN}}\left(t_i,x_i\right)-F_k^W\left(t_i,x_i\right)\right)+\log \left(1+{\left(w_{\text{BC}}^F\right)}_k\right)\right)}^2}},F=g,h,s_j,\end{array}$ (36)

其中${\left(w_m^F\right)}_k\ge 0,m=\text{Eq},\text{IC},\text{BC},k=1,\cdots ,N_v$ 是微观速度点$v_k$ 的自适应权重，可以通过神经网络进行训练，$\epsilon >0$ 是一个很小的数以防止出现零除现象。

至此，我们已经引入降维神经网络表示法。在第4章中将使用该方法求解BGK方程，并用微流问题验证方法的有效性。

4. 数值实验

本章对经典问题进行数值实验，从而验证DRNR方法的有效性和高效性，实验涉及两个稳态问题：一维Couette流和二维矩形风管流。

数值模拟中，我们使用初始学习率为${\eta }_0=0.001$ 的Adam优化器[16]，并结合热重启技术[17]，即在第$i$ 步训练中，学习率会根据余弦退火策略进行衰减

${\eta }_i=\frac{1}{2}{\eta }_0\left(1+\cos \left(\frac{i}{T_{\max }}\right)\right).$ (37)

此外，神经网络由五个隐藏层组成，每层有80个神经元，输出维数为$N_v={\displaystyle \prod _{i=1}^{D_x}{N}_i}$ ，$N_i=100$ ，$D_x$ 是物理空间的维数，(23)式中的多尺度常数取为$c_i=4^{i-1},i=1,2,3$ ，微观速度空间的计算域取为${\left[-10,10\right]}^{D_x}$ 。

4.1. 一维Couette流

本节对Couette流的变体进行研究，Couette流是一个涉及固壁边界条件的经典问题。该装置由处在

$x=\pm \frac{1}{2}$ 的两个无线平行板组成，左右板的速度方向不同，左板的速度为$u_l^W=\left(0,u_2^W,0\right)$ ，右板的速度为

$u_r^W=\left(0,0,u_3^W\right)$ ，左右板的温度为$T^W=1$ ，密度${\rho }^W$ 由(6)式确定。当使用降维BGK模型求解该问题时，由(32)式可得Couette流的边界条件(18)为

$\begin{array}{l}{g}^W\left({\rho }^W,u^W,T^W\right)=\frac{{\rho }^W}{\sqrt{2\pi T^W}}\exp \left(-\frac{{\left|v_1\right|}^2}{2T^W}\right),\\ h^W\left({\rho }^W,u^W,T^W\right)=\frac{{\left|u^W\right|}^2+2T^W}{2}{g}^W,\\ s_2^W\left({\rho }^W,u^W,T^W\right)=\left(u^W\cdot e_2\right)g^W,\\ s_3^W\left({\rho }^W,u^W,T^W\right)=\left(u^W\cdot e_3\right)g^W.\end{array}$ (38)

其中，$e_i,i=1,2,3$ 表示只有第$i$ 个元素非零的单位列向量。

数值模拟中，取$u_2^W=u_3^W=1$ ，在$x\in \left(-0.5,0.5\right)$ 中随机选择$N_{\text{PDE}}=500$ 个空间点，选择$x=\pm 0.5$ 处的固定点作为边界，那么$N_{\text{BC}}=2$ 。使用$\sigma \left(x\right)=\sin \left(x\right)$ 作为激活函数，总训练步数为10,000。Kn分别取$0.1,1,2.5$ ，参考解为DVM的数值解，则稳态时密度$\rho$ ，宏观速度$u_2,u_3$ ，温度$T$ ，应力张量${\sigma }_{13}$ 和热通量$q_1$ 的数值结果如图2所示。结果表明，对于所有变量，三个不同Kn对应的数值解和参考解都非常匹配。

Figure 2. Numerical solution of $\rho$ , $u_2$ , $u_3$ , $T$ , ${\sigma }_{13}$ and $q_1$ of the variant of Couette flow at steady state for $\text{Kn}=0.1,1,2.5$

图2. $\text{Kn}=0.1,1,2.5$ 时，稳态Couette流变体的数值解$\rho$ ，$u_2$ ，$u_3$ ，$T$ ，${\sigma }_{13}$ 和$q_1$

接下来，我们考虑微观速度降维至二维情况的Couette流，使用DRNR方法和DVM方法分别对其进行数值模拟，通过比较两者的误差在达到相同数量级时所需的时间，进行效率对比。数值模拟中，取$N_x=400$ 时DVM计算的数值解作为参考解，对于DVM方法，选择一阶迎风格式，计算时间为$t=10$ ，计算区域为$x\in \left(-0.5,0.5\right)$ ，空间网格点个数分别为$N_x=40,80$ ，对于DRNR方法，通过调整空间采样点个数和训练步数，分别计算两种方法得到的宏观量$\rho$ 与参考解的$L^2$ 误差，在DRNR方法和DVM方法得到的误差达到相同数量级时，记录两种方法的计算所需时间，相应结果如表1所示。结果表明，当空间点个数较少时，DVM的计算时间比DRNR的略少，但是随着空间点个数的增加，可以发现DRNR方法所需的计算时间要少于DVM的计算时间。在一定的精度范围内，DRNR方法随着精度的增加，逐渐比DVM方法更有优势。但是由于机器学习方法的局限性，无法达到更高的精度。

Table 1. Computational time of DRNR and DVM to achieve similar error for $\text{Kn}=1$ at $t=10$

表1. 在$\text{Kn}=1$ 、$t=10$ 的条件下，DRNR与DVM方法在达到相似误差时的计算时间

实验	实验一	实验二
$N_x$	40	80
误差	1.78e−03	8.53e−04
计算时间(DVM)	147	623
计算时间(DRNR)	267	537

4.2. 二维矩形风管流

本节研究二维矩形风管流的动力学。管道沿$z$ 轴无限延伸，计算区域为$xy$ 平面中的矩形横截面$x\in \left[0,1\right]\times \left[0,1\right]$ ，管道的四个侧面保持相同的温度$T^W$ ，管道的三个侧面静止$u_2^W=0$ ，顶面$y=1$ 处沿$z$ 轴以速度$u_1^W$ 移动，设

$u_1^W=u_{\text{max}}\left(\sin \left(\pi x\right)\right),$ (39)

其中$u_{\text{max}}=\left(0,0,0.5\right)$ 。需要注意的是，在上方两个角处速度是连续的，这样可以减轻奇异性对计算精度的影响。随着时间趋近于无穷大，流体在顶板运动的驱动下达到稳定状态。该装置目的是模拟管道内稀薄气体流的动力学，强调在可控二维环境中移动边界的影响。由于$z$ 轴不对称，该问题的边界条件为

$\begin{array}{l}{g}^W\left({\rho }^W,u^W,T^W\right)=\frac{{\rho }^W}{2\pi T^W}\exp \left(-\frac{{\left(v_1\right)}^2+{\left(v_2\right)}^2}{2T^W}\right),\\ h^W\left({\rho }^W,u^W,T^W\right)=\frac{{\left|u^W\right|}^2+T^W}{2}{g}^W,\\ s_3^W\left({\rho }^W,u^W,T^W\right)=\left(u^W\cdot e_3\right)g^W.\end{array}$ (40)

该问题的具体布局如图3所示。

Figure 3. The layout of the 2D rectangular duct flow problem

图3. 二维矩形风管问题示意图

在数值模拟中，随机选取$x\in \left(0,1\right)\times \left(0,1\right)$ 中的$N_{\text{PDE}}=2000$ 个点，每个边界使用$N_{\text{BC}}=300$ 个随机点进行离散，总训练步数为20,000。我们使用的激活函数为$\sigma \left(x\right)=\sin \left(x\right)$ ，可在特定的流体结构中表现出更好的性能。首先取Knudsen数$\text{Kn}=0.1$ ，参考解也是DVM的数值解，稳态时密度$\rho$ 、$z$ 轴的宏观速度$u_3$ 和温度$T$ 的数值结果如图4所示。结果表明，DRNR得到的数值解$\rho$ 和$u_3$ 与参考解吻合较好，然而数值解与参考解仍存在较小的差异，不过相对误差小于1%，推断这是由于Kn较小时BGK方程的非线性造成的。

接下来，Knudsen数增加到$\text{Kn}=1$ ，我们发现数值解和参考解是重叠的。将Knudsen数增加到$\text{Kn}=2.5$ ，这时数值解$u_3$ 和$T$ 与DVM得到的参考解具有很好的相关性，而密度$\rho$ 存在一些较小的差异，这表明DRNR方法具有高效性。

Figure 4. Numerical solution of $\rho$ , $u_3$ and $T$ at steady state for $\text{Kn}=0.1,1,2.5$

图4. $\text{Kn}=0.1,1,2.5$ 对应的稳态时的数值解$\rho$ ，$u_3$ 和$T$

5. 结论

本文采用神经网络的方法求解微流问题中的Boltzmann-BGK方程。首先通过分布函数的神经网络表示法构建降维模型，实现对BGK方程的高效近似。为处理微流问题中的Maxwell边界条件，设计了针对性的神经网络结构，有效地降低了网络参数的复杂性。进一步提出包含初始条件、边界条件和PDE残差的损失函数，极大地提高了神经网络的逼近效率。最后，通过一维和二维的经典问题验证了该降维神经网络表示方法的准确性和有效性。

致 谢

感谢北京大学的董彬教授、北京航空航天大学的张俊教授和北京计算科学研究中心的王艳莉研究员的宝贵建议。

基金项目

该工作得到了中国工程物理研究院院长基金项目(YZJJZQ2022017)、国家自然科学基金项目(编号：12171026、U2230402、12031013)的支持。

NOTES

^*通讯作者。