【注：如果被CSDN莫名其妙设置为VIP文章，就直接去我主页或知乎看 】
这是自己的一份笔记，主要是和格子Boltzmann方法的DnQb离散速度模型相关的。还请读者确定自己已经知道常见的DnQb模型长什么样再来读这篇文章，这样会好理解一些。
Gauss-Hermite积分与DnQb模型
一般来说，密度分布函数的Maxwell-Boltzmann平衡态分布为：
f
e
q
=
ρ
(
2
π
R
T
)
D
/
2
exp
[
−
(
ξ
−
u
)
2
2
R
T
]
,
(1-1)
f^{eq} = \\frac{\\rho}{ (2 \\pi R T)^{D/2} } \\exp \\left[ – \\frac{(\\boldsymbol{\\xi} – \\boldsymbol{u})^2}{2 R T} \\right], \\tag{1-1}
feq=(2πRT)D/2ρexp[−2RT(ξ−u)2],(1-1)
其中
ξ
\\boldsymbol{\\xi}
ξ 为分子速度；
R
R
R 为气体常数；
D
D
D为问题的维度数；
ρ
\\rho
ρ 为流场密度；
u
\\boldsymbol{u}
u 为流场速度。而 LBM 里常用的离散平衡态分布一般写作（沿着
c
α
\\boldsymbol{c}_\\alpha
cα 方向，权重为
W
α
W_\\alpha
Wα，
c
s
2
=
R
T
c_s^2 = RT
cs2=RT）：
f
α
e
q
=
ρ
W
α
×
[
1
+
ξ
α
⋅
u
c
s
2
+
(
ξ
α
⋅
u
)
2
2
c
s
4
−
∣
u
∣
2
2
c
s
2
]
.
(1-2)
f^{eq}_{\\alpha} = \\rho W_{\\alpha} \\times \\left[ 1 + \\frac{\\boldsymbol{\\xi}_{\\alpha} \\cdot \\boldsymbol{u}}{c_s^2} + \\frac{(\\boldsymbol{\\xi}_{\\alpha} \\cdot \\boldsymbol{u})^2}{2 c_s^4} – \\frac{|\\boldsymbol{u}|^2}{2 c_s^2} \\right]. \\tag{1-2}
fαeq=ρWα×[1+cs2ξα⋅u+2cs4(ξα⋅u)2−2cs2∣u∣2].(1-2)
那么这两者是怎么联系起来的呢？参考文献 [1] 讲的就是这个事情。最近自己总算是略懂一点了，这里把笔记写下来，作为备忘录。
平衡态分布的展开
当马赫数远小于1时，我们可以对式(1)的平衡态做如下离散：
f
e
q
(
ξ
)
=
ρ
(
2
π
R
T
)
D
/
2
exp
(
−
∣
ξ
∣
2
2
R
T
)
×
[
1
+
ξ
⋅
u
R
T
+
(
ξ
⋅
u
)
2
2
(
R
T
)
2
−
∣
u
∣
2
2
R
T
]
(1-3)
\\begin{aligned} f^{eq} (\\boldsymbol{\\xi}) &= \\frac{\\rho}{ (2 \\pi R T)^{D/2} } \\exp \\left( – \\frac{|\\boldsymbol{\\xi}|^2}{2 R T} \\right)\\\\ & \\quad \\times \\left[ 1 + \\frac{\\boldsymbol{\\xi} \\cdot \\boldsymbol{u}}{R T} + \\frac{(\\boldsymbol{\\xi} \\cdot \\boldsymbol{u})^2}{2 (R T)^2} – \\frac{|\\boldsymbol{u}|^2}{2 R T} \\right] \\end{aligned} \\tag{1-3}
feq(ξ)=(2πRT)D/2ρexp(−2RT∣ξ∣2)×[1+RTξ⋅u+2(RT)2(ξ⋅u)2−2RT∣u∣2](1-3)
式(1-3)中，我们对
u
\\boldsymbol{u}
u 进行展开，并且只截断到2阶。而格子Boltzmann方法（LBM）要做的，便是把下面的这个积分
∫
−
∞
+
∞
ψ
(
ξ
)
f
e
q
(
ξ
)
d
ξ
\\int_{-\\infty}^{+\\infty} \\psi(\\boldsymbol{\\xi}) f^{eq} (\\boldsymbol{\\xi}) \\mathrm{d}\\boldsymbol{\\xi}
∫−∞+∞ψ(ξ)feq(ξ)dξ
用离散的方式近似计算，其中
ψ
(
ξ
)
\\psi(\\boldsymbol{\\xi})
ψ(ξ) 是关于
ξ
\\boldsymbol{\\xi}
ξ 的多项式。下面我们将先以一维的情况为例，不去考虑多维的积分，简单讲讲数值离散的具体操作。
一维情形的离散 —— 以 D1Q3 为例
由于这一节只讨论一维问题，分子速度
ξ
\\xi
ξ 和宏观速度
u
u
u 即是标量值（不需要向量形式的加粗），并且维度
D
=
1
D=1
D=1。所以，式(3)写作
f
e
q
(
ξ
)
=
ρ
(
2
π
R
T
)
1
/
2
exp
(
−
ξ
2
2
R
T
)
×
[
1
+
ξ
u
R
T
+
(
ξ
u
)
2
2
(
R
T
)
2
−
u
2
2
R
T
]
.
(1-4)
\\begin{aligned} f^{eq} (\\xi) &= \\frac{\\rho}{ (2 \\pi R T)^{1/2} } \\exp \\left( – \\frac{\\xi^2}{2 R T} \\right)\\\\ & \\quad \\times \\left[ 1 + \\frac{\\xi u}{R T} + \\frac{(\\xi u)^2}{2 (R T)^2} – \\frac{u^2}{2 R T} \\right]. \\end{aligned} \\tag{1-4}
feq(ξ)=(2πRT)1/2ρexp(−2RTξ2)×[1+RTξu+2(RT)2(ξu)2−2RTu2].(1-4)
积分形式的抽象
式(1-4)的中括号内的四项可以用积分的加法拆成四个积分，并且每一项都可以看成如下形式：
Constant
⋅
∫
−
∞
+
∞
exp
(
−
ξ
2
2
R
T
)
ψ
(
ξ
)
d
ξ
,
\\text{Constant} \\cdot \\int_{-\\infty}^{+\\infty} \\exp \\left( – \\frac{\\xi^2}{2 R T} \\right) \\psi(\\xi) \\mathrm{d}\\xi,
Constant⋅∫−∞+∞exp(−2RTξ2)ψ(ξ)dξ,
其中
ψ
(
ξ
)
\\psi(\\xi)
ψ(ξ) 为关于
ξ
\\xi
ξ 的多项式。那么，我们令
ψ
m
(
ξ
)
=
ξ
m
\\psi_{m}(\\xi) = \\xi^{m}
ψm(ξ)=ξm （
m
m
m 为自然数），将研究目标变为式(1-5)的通式，即：
ρ
(
2
π
R
T
)
1
/
2
⋅
∫
−
∞
+
∞
exp
(
−
ξ
2
2
R
T
)
ψ
m
(
ξ
)
d
ξ
.
(1-5)
\\frac{\\rho}{ (2 \\pi R T)^{1/2} } \\cdot \\int_{-\\infty}^{+\\infty} \\exp \\left( – \\frac{\\xi^2}{2 R T} \\right) \\psi_{m}(\\xi) \\mathrm{d}\\xi. \\tag{1-5}
(2πRT)1/2ρ⋅∫−∞+∞exp(−2RTξ2)ψm(ξ)dξ.(1-5)
为简化这个积分的计算，这里需要引入变量代换
ζ
=
ξ
/
2
R
T
,
\\zeta = \\xi / \\sqrt{2 R T},
ζ=ξ/2RT
,
则
d
ξ
=
2
R
T
d
ζ
.
\\mathrm{d}\\xi = \\sqrt{2 R T} \\mathrm{d}\\zeta.
dξ=2RT
dζ.
代换后，有：
ρ
(
2
π
R
T
)
1
/
2
⋅
∫
−
∞
+
∞
exp
(
−
ζ
2
)
ψ
m
(
2
R
T
ζ
)
2
R
T
d
ζ
=
ρ
(
2
π
R
T
)
1
/
2
⋅
(
2
R
T
)
(
m
+
1
)
/
2
⋅
I
m
=
ρ
π
⋅
(
2
R
T
)
m
/
2
⋅
I
m
\\begin{aligned} & \\frac{\\rho}{ (2 \\pi R T)^{1/2} } \\cdot \\int_{-\\infty}^{+\\infty} \\exp \\left( – \\zeta^2 \\right) \\psi_{m}(\\sqrt{2 R T} \\zeta) \\sqrt{2 R T} \\mathrm{d}\\boldsymbol{\\zeta} \\\\=& \\frac{\\rho}{ (2 \\pi R T)^{1/2} } \\cdot (2 R T)^{(m+1)/2} \\cdot I_m \\\\=& \\frac{\\rho}{\\sqrt{\\pi}} \\cdot (2 R T)^{m/2} \\cdot I_m \\end{aligned}
==(2πRT)1/2ρ⋅∫−∞+∞exp(−ζ2)ψm(2RT
ζ)2RT
dζ(2πRT)1/2ρ⋅(2RT)(m+1)/2⋅Imπ
ρ⋅(2RT)m/2⋅Im
其中
I
m
=
∫
−
∞
+
∞
exp
(
−
ζ
2
)
ζ
m
d
ζ
\\displaystyle I_m = \\int_{-\\infty}^{+\\infty} \\exp \\left( – \\zeta^2 \\right) \\zeta^m \\mathrm{d}\\boldsymbol{\\zeta}
Im=∫−∞+∞exp(−ζ2)ζmdζ，且
m
m
m 为自然数。
所以，下面的积分会被展开为
I
=
∫
−
∞
+
∞
ψ
m
(
ξ
)
f
e
q
(
ξ
)
d
ξ
=
ρ
π
[
(
1
−
u
2
2
R
T
)
(
2
R
T
)
m
/
2
I
m
+
u
R
T
(
2
R
T
)
(
m
+
1
)
/
2
I
m
+
1
+
u
2
2
(
R
T
)
2
(
2
R
T
)
(
m
+
2
)
/
2
I
m
+
2
]
.
(1-6)
\\begin{aligned} I=& \\int_{-\\infty}^{+\\infty} \\psi_m(\\xi) f^{eq} (\\xi) \\mathrm{d}\\xi \\\\=& \\frac{\\rho}{\\sqrt{\\pi}} \\left[ \\left(1 – \\frac{u^2}{2 R T}\\right) (2 R T)^{m/2} I_m + \\frac{u}{RT} (2 R T)^{(m+1)/2} I_{m+1} \\right. \\\\&\\qquad \\left. + \\frac{u^2}{2(RT)^2} (2 R T)^{(m+2)/2} I_{m+2} \\right]. \\end{aligned} \\tag{1-6}
I==∫−∞+∞ψm(ξ)feq(ξ)dξπ
ρ[(1−2RTu2)(2RT)m/2Im+RTu(2RT)(m+1)/2Im+1+2(RT)2u2(2RT)(m+2)/2Im+2].(1-6)
数值积分的构造
假设我们想用一个三点格式的Gauss-Hermite数值积分进行近似，则：
∫
−
∞
+
∞
exp
(
−
ζ
2
)
ϕ
(
ζ
)
d
ζ
=
∑
i
∈
[
−
1
,
0
,
1
]
w
i
ϕ
(
ζ
i
)
.
\\int_{-\\infty}^{+\\infty} \\exp(- \\zeta^2) \\phi(\\zeta) \\mathrm{d} \\zeta = \\sum_{i \\in [-1, 0, 1]} w_{i} \\phi(\\zeta_{i}).
∫−∞+∞exp(−ζ2)ϕ(ζ)dζ=i∈[−1,0,1]∑wiϕ(ζi).
其中：
ζ
0
=
0
,
ζ
1
=
3
/
2
,
ζ
−
1
=
−
3
/
2
,
w
0
=
2
π
/
3
,
w
1
=
π
/
6
,
w
−
1
=
π
/
6.
\\begin{aligned} \\zeta_0 = 0 ,\\quad& \\zeta_{1} = \\sqrt{3/2} ,\\quad& \\zeta_{-1} = -\\sqrt{3/2} ,\\\\ w_0 = 2\\sqrt{\\pi} / 3 ,\\quad& w_{1} = \\sqrt{\\pi} / 6 ,\\quad& w_{-1} = \\sqrt{\\pi} / 6. \\end{aligned}
ζ0=0,w0=2π
/3,ζ1=3/2
,w1=π
/6,ζ−1=−3/2
,w−1=π
/6.
如果用变换关系代换回
ξ
\\xi
ξ，则
ξ
0
=
0
\\xi_0 = 0
ξ0=0,
ξ
−
1
=
−
3
R
T
\\xi_{-1} = -\\sqrt{3RT}
ξ−1=−3RT
和
ξ
1
=
3
R
T
\\xi_{1} = \\sqrt{3RT}
ξ1=3RT
.
所以，积分
I
m
I_m
Im 可以被近似成：
I
m
=
∫
−
∞
+
∞
exp
(
−
ζ
2
)
ζ
m
d
ζ
≈
∑
i
∈
[
−
1
,
0
,
1
]
w
i
ζ
i
m
=
(
2
R
T
)
−
m
/
2
∑
i
∈
[
−
1
,
0
,
1
]
w
i
ξ
i
m
\\begin{aligned} I_m &= \\int_{-\\infty}^{+\\infty} \\exp \\left( – \\zeta^2 \\right) \\zeta^m \\mathrm{d}\\boldsymbol{\\zeta} \\\\&\\approx \\sum_{i \\in [-1, 0, 1]} w_{i} \\zeta_{i}^m \\\\&= (2RT)^{-m/2} \\sum_{i \\in [-1, 0, 1]} w_{i} \\xi_{i}^m \\end{aligned}
Im=∫−∞+∞exp(−ζ2)ζmdζ≈i∈[−1,0,1]∑wiζim=(2RT)−m/2i∈[−1,0,1]∑wiξim
也就是说，
(
2
R
T
)
m
/
2
I
m
≈
∑
i
∈
[
−
1
,
0
,
1
]
w
i
ξ
i
m
(2 R T)^{m/2} I_m \\approx \\sum_{i \\in [-1, 0, 1]} w_{i} \\xi_{i}^m
(2RT)m/2Im≈i∈[−1,0,1]∑wiξim
综上所述，式(1-6)就会被改写为：
I
=
∫
−
∞
+
∞
ψ
m
(
ξ
)
f
e
q
(
ξ
)
d
ξ
≈
∑
i
∈
[
−
1
,
0
,
1
]
ξ
i
m
⋅
{
ρ
w
i
π
⋅
[
1
+
ξ
i
u
R
T
+
(
ξ
i
u
)
2
2
(
R
T
)
2
−
u
2
2
R
T
]
}
.
\\begin{aligned} I=& \\int_{-\\infty}^{+\\infty} \\psi_m(\\xi) f^{eq} (\\xi) \\mathrm{d}\\xi \\\\ \\approx& \\sum_{i \\in [-1, 0, 1]} \\xi_{i}^{m} \\cdot \\left\\{ \\rho \\frac{w_{i}}{\\sqrt{\\pi}} \\cdot \\left[ 1 + \\frac{\\xi_{i} u}{R T} + \\frac{(\\xi_{i} u)^2}{2 (R T)^2} – \\frac{u^2}{2 R T} \\right] \\right\\}. \\end{aligned}
I=≈∫−∞+∞ψm(ξ)feq(ξ)dξi∈[−1,0,1]∑ξim⋅{ρπ
wi⋅[1+RTξiu+2(RT)2(ξiu)2−2RTu2]}.
到这一步，我们就可以定义
ξ
i
\\xi_{i}
ξi 方向的离散平衡态
f
i
e
q
f^{eq}_{i}
fieq 为：
f
i
e
q
=
ρ
W
i
⋅
[
1
+
ξ
i
u
c
s
2
+
(
ξ
i
u
)
2
2
c
s
4
−
u
2
2
c
s
2
]
.
(1-7)
f^{eq}_{i} = \\rho W_{i} \\cdot \\left[ 1 + \\frac{\\xi_{i} u}{c_s^2} + \\frac{(\\xi_{i} u)^2}{2 c_s^4} – \\frac{u^2}{2 c_s^2} \\right]. \\tag{1-7}
fieq=ρWi⋅[1+cs2ξiu+2cs4(ξiu)2−2cs2u2].(1-7)
其中的权重
W
i
=
w
i
/
π
W_{i} = w_i / \\sqrt{\\pi}
Wi=wi/π
。并且，如果令
c
=
3
R
T
c = \\sqrt{3RT}
c=3RT
， 则
c
s
=
R
T
=
c
/
3
c_s = \\sqrt{RT} = c / \\sqrt{3}
cs=RT
=c/3
。
至此，我们就整理出了 D1Q3 模型：
| 0 | 0 | 2/3 |
| -c | -c | 1/6 |
| c | c | 1/6 |
注意：在多数LBM实践中，我们通常将
c
c
c 无量纲化，即令
c
=
1
c=1
c=1。
更高维的扩展
由于各个坐标轴之间是正交的，这种做法同样可以拓展到高维空间。
打个比方，对于三维问题，
ξ
=
[
ξ
x
,
ξ
y
,
ξ
z
]
T
\\boldsymbol{\\xi} = [\\xi_x, \\xi_y, \\xi_z]^{\\mathrm{T}}
ξ=[ξx,ξy,ξz]T， 则
d
ξ
=
d
ξ
x
d
ξ
y
d
ξ
z
\\mathrm{d}\\boldsymbol{\\xi} = \\mathrm{d}\\xi_x \\mathrm{d}\\xi_y \\mathrm{d}\\xi_z
dξ=dξxdξydξz。积分时就只需要将不同方向的分量分离出来就可以了。
对于二维问题，记
ψ
m
,
n
(
ξ
)
=
ξ
x
m
ξ
y
n
\\psi_{m,n}(\\boldsymbol{\\xi})=\\xi_x^{m} \\xi_y^{n}
ψm,n(ξ)=ξxmξyn，则
I
=
∫
ψ
m
,
n
(
ξ
)
f
(
s
t
)
d
ξ
=
ρ
π
(
2
R
T
)
m
+
n
{
(
1
−
u
2
2
R
T
)
I
m
I
n
+
1
2
R
T
2
(
u
x
I
m
+
1
I
n
+
u
y
I
m
I
n
+
1
)
+
1
R
T
(
u
x
2
I
m
+
2
I
n
+
2
u
x
u
y
I
m
+
1
I
n
+
1
+
u
y
2
I
m
I
n
+
2
)
}
(1-8)
\\begin{aligned} I =& \\int\\psi_{m , n} (\\boldsymbol{\\xi}) f^{( \\mathrm{s t} )} \\mathrm{d}\\boldsymbol{\\xi} \\\\=& \\frac{\\rho}{\\pi} ( \\sqrt{2 R T} )^{m+n} \\Biggl\\{\\left( 1-\\frac{u^{2}} {2 R T} \\right) I_{m} I_{n} \\\\&+ \\frac{1} {\\sqrt{2 R T}} 2 ( u_{x} I_{m+1} I_{n}+u_{y} I_{m} I_{n+1} ) \\\\&+ \\frac{1} {R T} (u_{x}^{2} I_{m+2} I_{n}+2 u_{x} u_{y} I_{m+1} I_{n+1}+u_{y}^{2} I_{m} I_{n+2}) \\Biggr\\} \\end{aligned} \\tag{1-8}
I==∫ψm,n(ξ)f(st)dξπρ(2RT
)m+n{(1−2RTu2)ImIn+2RT
12(uxIm+1In+uyImIn+1)+RT1(ux2Im+2In+2uxuyIm+1In+1+uy2ImIn+2)}(1-8)
式(1-8)同样可以整理成：
I
=
∑
i
,
j
∈
[
−
1
,
0
,
1
]
ψ
(
ξ
i
,
j
)
⋅
{
ρ
w
i
w
j
π
[
1
+
ξ
i
,
j
⋅
u
c
s
2
+
(
ξ
i
,
j
⋅
u
)
2
2
c
s
4
−
∥
u
∥
2
2
c
s
2
]
}
(1-9)
I = \\sum_{i,j\\in[-1,0,1]} \\psi(\\boldsymbol{\\xi}_{i,j}) \\cdot \\left\\{ \\rho \\frac{w_i w_j}{\\pi} \\left[ 1 + \\frac{\\boldsymbol{\\xi}_{i,j} \\cdot \\boldsymbol{u}}{c_s^2} + \\frac{(\\boldsymbol{\\xi}_{i,j} \\cdot \\boldsymbol{u})^2}{2 c_s^4} – \\frac{\\|\\boldsymbol{u}\\|^2}{2 c_s^2} \\right] \\right\\} \\tag{1-9}
I=i,j∈[−1,0,1]∑ψ(ξi,j)⋅{ρπwiwj[1+cs2ξi,j⋅u+2cs4(ξi,j⋅u)2−2cs2∥u∥2]}(1-9)
其中
ξ
i
,
j
=
2
R
T
(
ζ
i
,
ζ
j
)
\\boldsymbol{\\xi}_{i,j}=\\sqrt{2RT} (\\zeta_i,\\zeta_j)
ξi,j=2RT
(ζi,ζj)。令
W
i
,
j
=
w
i
w
j
/
π
W_{i,j} = w_i w_j / \\pi
Wi,j=wiwj/π ，便得出了 D2Q9 模型的权重。
类似的，对于三维问题，则可以整理出
I
=
∑
i
,
j
,
k
∈
[
−
1
,
0
,
1
]
ψ
(
ξ
i
,
j
,
k
)
⋅
{
ρ
w
i
w
j
w
k
π
3
/
2
[
1
+
ξ
i
,
j
,
k
⋅
u
c
s
2
+
(
ξ
i
,
j
,
k
⋅
u
)
2
2
c
s
4
−
∥
u
∥
2
2
c
s
2
]
}
(1-10)
I = \\sum_{i,j,k\\in[-1,0,1]} \\psi(\\boldsymbol{\\xi}_{i,j,k}) \\cdot \\left\\{ \\rho \\frac{w_i w_j w_k}{\\pi^{3/2}} \\left[ 1 + \\frac{\\boldsymbol{\\xi}_{i,j,k} \\cdot \\boldsymbol{u}}{c_s^2} + \\frac{(\\boldsymbol{\\xi}_{i,j,k} \\cdot \\boldsymbol{u})^2}{2 c_s^4} – \\frac{\\|\\boldsymbol{u}\\|^2}{2 c_s^2} \\right] \\right\\} \\tag{1-10}
I=i,j,k∈[−1,0,1]∑ψ(ξi,j,k)⋅{ρπ3/2wiwjwk[1+cs2ξi,j,k⋅u+2cs4(ξi,j,k⋅u)2−2cs2∥u∥2]}(1-10)
其中
ξ
i
,
j
,
k
=
2
R
T
(
ζ
i
,
ζ
j
,
ζ
k
)
\\boldsymbol{\\xi}_{i,j,k}=\\sqrt{2RT} (\\zeta_i,\\zeta_j,\\zeta_k)
ξi,j,k=2RT
(ζi,ζj,ζk)。令
W
i
,
j
,
k
=
w
i
w
j
w
k
/
π
3
/
2
W_{i,j,k} = w_i w_j w_k / \\pi^{3/2}
Wi,j,k=wiwjwk/π3/2 ，便得出了 D3Q27 模型的权重。
待定系数法与DnQb模型
上面的内容可以说从数学上给出了明确的证明。然而实际使用中，其实还是用待定系数法来敲定 DnQb 模型更多一些。下面的内容里，为了保持和何雅玲老师的书标号一致，这里将第
α
\\alpha
α 个格子速度向量记作
e
α
\\boldsymbol{e}_{\\alpha}
eα，其权重也记作
ω
α
\\omega_{\\alpha}
ωα。
一般来说，d 维空间下的格子速度的 n 阶张量为：
E
i
1
i
2
.
.
.
i
n
(
n
)
=
∑
α
(
e
α
)
i
1
(
e
α
)
i
2
.
.
.
(
e
α
)
i
n
(2-1)
E_{i_1 i_2 … i_n}^{(n)} = \\sum_{\\alpha} (\\boldsymbol{e}_{\\alpha})_{i_1} (\\boldsymbol{e}_{\\alpha})_{i_2} … (\\boldsymbol{e}_{\\alpha})_{i_n} \\tag{2-1}
Ei1i2…in(n)=α∑(eα)i1(eα)i2…(eα)in(2-1)
假如你选取的离散速度集（共有
M
M
M 个速度）是各向同性的，则应该满足
E
i
1
i
2
.
.
.
i
2
n
+
1
(
2
n
+
1
)
=
0
,
E^{(2n+1)}_{i_1 i_2 … i_{2n+1}} = 0,
Ei1i2…i2n+1(2n+1)=0,
和
E
i
1
i
2
.
.
.
i
2
n
(
2
n
)
=
M
d
(
d
+
2
)
.
.
.
(
d
+
2
n
−
2
)
Δ
i
1
i
2
.
.
.
i
2
n
(
2
n
)
E^{(2n)}_{i_1 i_2 … i_{2n}} = \\frac{M}{d (d+2) … (d+2n-2)} \\Delta^{(2n)}_{i_1 i_2 … i_{2n}}
Ei1i2…i2n(2n)=d(d+2)…(d+2n−2)MΔi1i2…i2n(2n)
下面我们首先从张量算子
Δ
(
2
n
)
\\Delta^{(2n)}
Δ(2n) 开始介绍。详情可见文献 [2] 的第3章，这篇文章在scihub上可以直接下载。
格子张量
张量算子
Δ
(
2
n
)
\\Delta^{(2n)}
Δ(2n) 实际上计算的是将
2
n
2n
2n 个指标全部两两配对的所有方式之和。它的具体形式为：
Δ
i
j
(
2
)
=
δ
i
j
Δ
i
j
k
l
(
4
)
=
δ
i
j
δ
k
l
+
δ
i
k
δ
j
l
+
δ
i
l
δ
j
k
Δ
i
j
k
l
p
q
(
6
)
=
δ
i
j
Δ
i
j
k
l
(
4
)
+
δ
i
k
Δ
j
l
p
q
(
4
)
+
δ
i
l
Δ
j
k
p
q
(
4
)
+
δ
i
p
Δ
j
k
l
q
(
4
)
+
δ
i
q
Δ
j
k
l
p
(
4
)
\\begin{aligned} \\Delta^{(2)}_{ij}&=\\delta_{ij} \\\\ \\Delta^{(4)}_{ijkl}&=\\delta_{ij}\\delta_{kl} + \\delta_{ik}\\delta_{jl} + \\delta_{il}\\delta_{jk} \\\\ \\Delta^{(6)}_{ijklpq}&=\\delta_{ij}\\Delta^{(4)}_{ijkl} + \\delta_{ik}\\Delta^{(4)}_{jlpq} + \\delta_{il}\\Delta^{(4)}_{jkpq} + \\delta_{ip}\\Delta^{(4)}_{jklq} + \\delta_{iq}\\Delta^{(4)}_{jklp} \\end{aligned}
Δij(2)Δijkl(4)Δijklpq(6)=δij=δijδkl+δikδjl+δilδjk=δijΔijkl(4)+δikΔjlpq(4)+δilΔjkpq(4)+δipΔjklq(4)+δiqΔjklp(4)
并存在如下递归关系：
Δ
i
1
i
2
…
i
2
n
(
2
n
)
=
∑
j
=
2
2
n
δ
i
1
,
i
j
⋅
Δ
i
2
…
i
j
−
1
i
j
+
1
…
i
2
n
(
2
n
−
2
)
\\Delta^{(2n)}_{i_1 i_2 \\dots i_{2n}} = \\sum_{j=2}^{2n} \\delta_{i_1, i_{j}} \\cdot \\Delta^{(2n-2)}_{i_2 \\dots i_{j-1} i_{j+1} \\dots i_{2n}}
Δi1i2…i2n(2n)=j=2∑2nδi1,ij⋅Δi2…ij−1ij+1…i2n(2n−2)
其中
δ
i
j
\\delta_{ij}
δij 为克罗内克积。
2
n
2n
2n 个指标分成
n
n
n 对再求和，共有
(
2
n
−
1
)
!
!
(2n−1)!!
(2n−1)!! 项。 计算思路为 ：第一个元素有
(
2
n
−
1
)
(2n−1)
(2n−1) 个选择配对；剩下
(
2
n
−
2
)
(2n−2)
(2n−2) 个元素之中的第一个有
(
2
n
−
3
)
(2n−3)
(2n−3) 个选择；后续依此类推。
(
2
n
)
!
2
n
n
!
=
总积
偶数部分的积
\\frac{(2n)!}{2^n n!} = \\frac{\\text{总积}}{\\text{偶数部分的积}}
2nn!(2n)!=偶数部分的积总积，结果是等价的。 计算思路为 ：
2
n
2n
2n 个元素的全排列，除以每对内部的顺序(
2
n
2^n
2n)，再除以
n
n
n 对之间的顺序(
n
!
n!
n!)。
为了简洁表示
Δ
(
2
n
)
\\Delta^{(2n)}
Δ(2n) 的分量，文献 [2] 引入了 upper simplicial components，其索引是非递增的多重指标（如
(
i
1
,
i
2
,
…
,
i
k
)
(i_1,i_2,\\dots,i_k)
(i1,i2,…,ik)，满足
i
1
≥
i
2
≥
⋯
≥
i
k
i_1 \\geq i_2 \\geq \\dots \\geq i_k
i1≥i2≥⋯≥ik）。 假设我们讨论一下二维空间的
Δ
i
j
k
l
(
4
)
\\Delta^{(4)}_{ijkl}
Δijkl(4) ，其中的
i
,
j
,
k
,
l
i,j,k,l
i,j,k,l 取 1 或 2。 那么下标的所有排列便是：
- 1111（全1，非递增）→ 值3
- 2111（1个2、3个1，非递增）→ 值0
- 2211（2个2、2个1，非递增）→ 值1
- 2221（3个2、1个1，非递增）→ 值0
- 2222（全2，非递增）→ 值3
即
Δ
(
4
)
=
[
3
,
0
,
1
,
0
,
3
]
\\Delta^{(4)}=[3,0,1,0,3]
Δ(4)=[3,0,1,0,3].
在后续的描述中，为了避免写一大堆下标，我们采用如下简记
E
(
n
)
=
E
i
1
i
2
…
i
n
(
n
)
,
Δ
(
2
n
)
=
Δ
i
1
i
2
…
i
2
n
(
2
n
)
,
δ
(
2
n
)
=
δ
i
1
i
2
…
i
2
n
(
2
n
)
,
E^{(n)} = E^{(n)}_{i_1 i_2 \\dots i_{n}}, \\quad \\Delta^{(2n)} = \\Delta^{(2n)}_{i_1 i_2 \\dots i_{2n}}, \\quad \\delta^{(2n)} = \\delta^{(2n)}_{i_1 i_2 \\dots i_{2n}},
E(n)=Ei1i2…in(n),Δ(2n)=Δi1i2…i2n(2n),δ(2n)=δi1i2…i2n(2n),
这里的
Δ
(
2
n
)
\\Delta^{(2n)}
Δ(2n) 就是
{
i
1
,
i
2
,
…
i
2
n
}
\\{i_1, i_2, \\dots i_{2n}\\}
{i1,i2,…i2n} 一系列下标组成的
Δ
i
1
i
2
…
i
2
n
(
2
n
)
\\Delta^{(2n)}_{i_1 i_2 \\dots i_{2n}}
Δi1i2…i2n(2n)。而
δ
(
2
n
)
\\delta^{(2n)}
δ(2n) 则是2n个下标的克罗内克积，当所有下标相等时才等于1，否则为0。
速度矩与权重系数
从
E
(
2
n
)
E^{(2n)}
E(2n) 的定义中，若假设所有
e
α
\\boldsymbol{e}_{\\alpha}
eα 一样长，可得：
1
M
∑
α
(
e
α
⋅
v
)
2
n
=
(
2
n
−
1
)
!
!
d
(
d
+
2
)
.
.
.
(
d
+
2
n
−
2
)
⏟
Q
2
n
∥
v
∥
2
n
\\frac{1}{M} \\sum_{\\alpha} (\\boldsymbol{e}_{\\alpha} \\cdot \\boldsymbol{v})^{2n} = \\underbrace{\\frac{(2n-1)!!}{d (d+2) … (d+2n-2)}}_{Q_{2n}} \\|\\boldsymbol{v}\\|^{2n}
M1α∑(eα⋅v)2n=Q2n
d(d+2)…(d+2n−2)(2n−1)!!∥v∥2n
对于二维情况 (
d
=
2
d=2
d=2)，则
Q
2
=
1
/
2
Q_2 = 1/2
Q2=1/2,
Q
4
=
3
/
8
Q_4 = 3/8
Q4=3/8,
Q
6
=
5
/
16
Q_6 = 5/16
Q6=5/16,
Q
8
=
35
/
128
Q_8 = 35/128
Q8=35/128。对于三维情况 (
d
=
3
d=3
d=3)，则
Q
2
n
=
1
/
(
2
n
+
1
)
Q_{2n} = 1 / (2n+1)
Q2n=1/(2n+1).
同理，拓展到奇数阶，则有：
1
M
∑
α
(
e
α
⋅
v
)
2
n
(
e
α
⋅
v
)
=
Q
2
n
∥
v
∥
2
n
v
\\frac{1}{M} \\sum_{\\alpha} (\\boldsymbol{e}_{\\alpha} \\cdot \\boldsymbol{v})^{2n} (\\boldsymbol{e}_{\\alpha} \\cdot \\boldsymbol{v}) = Q_{2n} \\|\\boldsymbol{v}\\|^{2n} \\boldsymbol{v}
M1α∑(eα⋅v)2n(eα⋅v)=Q2n∥v∥2nv
因此，对于一个所有
e
α
\\boldsymbol{e}_{\\alpha}
eα 一样长的离散速度集来说，它的离散速度矩的权重是确定的一个数 ——
Q
2
n
Q_{2n}
Q2n。
假如
e
α
{\\boldsymbol{e}_{\\alpha}}
eα 中的所有速度并不等长，我们则可以认为：
E
i
1
i
2
.
.
.
i
n
(
n
)
=
∑
α
ω
(
∥
e
α
∥
)
⋅
(
e
α
)
i
1
(
e
α
)
i
2
.
.
.
(
e
α
)
i
n
E_{i_1 i_2 … i_n}^{(n)} = \\sum_{\\alpha} \\omega(\\|\\boldsymbol{e}_{\\alpha}\\|) \\cdot (\\boldsymbol{e}_{\\alpha})_{i_1} (\\boldsymbol{e}_{\\alpha})_{i_2} … (\\boldsymbol{e}_{\\alpha})_{i_n}
Ei1i2…in(n)=α∑ω(∥eα∥)⋅(eα)i1(eα)i2…(eα)in
其中
ω
(
∥
e
α
∥
)
\\omega(\\|\\boldsymbol{e}_{\\alpha}\\|)
ω(∥eα∥) 指代由特征速度的模
∥
e
α
∥
\\|\\boldsymbol{e}_{\\alpha}\\|
∥eα∥ 决定的权重系数。
DnQb 模型的确定
Qian等的 DnQb 模型采用如下的平衡态分布：
f
α
e
q
(
ρ
,
u
;
e
α
,
ω
α
)
=
ω
α
ρ
⋅
[
1
+
e
α
⋅
u
c
s
2
+
(
e
α
⋅
u
)
2
2
c
s
4
−
∥
u
∥
2
2
c
s
2
]
.
f^{eq}_{\\alpha} (\\rho, \\boldsymbol{u}; \\boldsymbol{e}_{\\alpha}, \\omega_{\\alpha}) = \\omega_{\\alpha} \\rho \\cdot \\left[ 1 + \\frac{\\boldsymbol{e}_{\\alpha} \\cdot \\boldsymbol{u}}{c_s^2} + \\frac{(\\boldsymbol{e}_{\\alpha} \\cdot \\boldsymbol{u})^2}{2 c_s^4} – \\frac{\\|\\boldsymbol{u}\\|^2}{2 c_s^2} \\right].
fαeq(ρ,u;eα,ωα)=ωαρ⋅[1+cs2eα⋅u+2cs4(eα⋅u)2−2cs2∥u∥2].
其中
ρ
\\rho
ρ 为流场密度；
u
\\boldsymbol{u}
u 为流场速度，
c
s
2
=
R
T
c_s^2 = RT
cs2=RT。若是为了模拟Navier-Stokes方程，平衡态分布函数需要满足下面的约束：
∑
α
f
α
e
q
=
ρ
,
∑
α
f
α
e
q
e
α
=
ρ
u
,
∑
α
f
α
e
q
e
α
i
e
α
j
=
ρ
u
i
u
j
+
p
δ
i
j
.
(2-2)
\\begin{aligned} \\sum_{\\alpha} f^{eq}_{\\alpha} &= \\rho ,\\\\ \\sum_{\\alpha} f^{eq}_{\\alpha} \\boldsymbol{e}_{\\alpha} &= \\rho \\boldsymbol{u} ,\\\\ \\sum_{\\alpha} f^{eq}_{\\alpha} e_{\\alpha i} e_{\\alpha j} &= \\rho u_i u_j + p \\delta_{ij}. \\end{aligned} \\tag{2-2}
α∑fαeqα∑fαeqeαα∑fαeqeαieαj=ρ,=ρu,=ρuiuj+pδij.(2-2)
在待定系数法的框架下，就是把已经设计好的离散速度集
{
e
α
}
\\{\\boldsymbol{e}_{\\alpha}\\}
{eα} 代入上面的约束中，然后求解对应的权重
ω
α
\\omega_{\\alpha}
ωα。
在这种视角下，我们要怎么确定就是要用 D2Q9、D3Q15、D3Q9 这些形式呢？
首先，正四边形最多能保证三阶格子张量是各向同性的，而到了四阶就变成了
E
i
j
k
l
(
4
)
∣
M
=
4
=
2
δ
i
j
k
l
E_{ijkl}^{(4)}|_{M=4} = 2 \\delta_{ijkl}
Eijkl(4)∣M=4=2δijkl。为了解决这个问题，我们可以在在原有最近邻链接基础上，增加更远邻的链接，扩大速度方向集合
{
e
α
}
\\{\\boldsymbol{e}_{\\alpha}\\}
{eα}。 而正方晶格加对角线的方式，就构成了典型的 D2Q9 模型：此时速度方向等效于正八边形顶点，当
M
=
8
>
4
M=8>4
M=8>4 时，
E
(
4
)
E^{(4)}
E(4) 满足各向同性条件（
M
M
M 不整除4）。
同样的道理，在三维空间中若是使用立方晶格，此时速度方向为立方体顶点（
M
=
8
M=8
M=8，最近邻链接），对称性较低。
E
(
4
)
E^{(4)}
E(4) 也是各向异性的，无法导出标准流体方程。此时也需要把一些对角线方向加进来，让
E
(
4
)
E^{(4)}
E(4) 满足各向同性条件。
Example: 以 D2Q9 为例
首先，我们把 D2Q9 分成三组
| 1 |
( 0 , 0 ) (0,0) (0,0) |
ω 0 \\omega_0 ω0 |
| 2 | cyc:
( ± c , 0 ) (\\pm c,0) (±c,0) |
ω 1 = ω 2 = ω 3 = ω 4 \\omega_1=\\omega_2=\\omega_3=\\omega_4 ω1=ω2=ω3=ω4 |
| 3 |
( ± c , ± c ) (\\pm c,\\pm c) (±c,±c) |
ω 5 = ω 6 = ω 7 = ω 8 \\omega_5=\\omega_6=\\omega_7=\\omega_8 ω5=ω6=ω7=ω8 |
对于组别 2，根据格子张量的表格，我们有
∑
α
∈
[
1
−
4
]
e
α
i
=
0
,
∑
α
∈
[
1
−
4
]
e
α
i
e
α
j
=
2
δ
i
j
⋅
c
2
∑
α
∈
[
1
−
4
]
e
α
i
e
α
j
e
α
k
=
0
,
∑
α
∈
[
1
−
4
]
e
α
i
e
α
j
e
α
k
e
α
l
=
2
δ
i
j
k
l
⋅
c
4
\\begin{aligned} \\sum_{\\alpha\\in[1-4]} e_{\\alpha i} &= 0 ,& \\sum_{\\alpha\\in[1-4]} e_{\\alpha i} e_{\\alpha j} &= 2\\delta_{ij} \\cdot c^2\\\\ \\sum_{\\alpha\\in[1-4]} e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} &= 0 ,& \\sum_{\\alpha\\in[1-4]} e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} e_{\\alpha l} &= 2\\delta_{ijkl} \\cdot c^4\\\\ \\end{aligned}
α∈[1−4]∑eαiα∈[1−4]∑eαieαjeαk=0,=0,α∈[1−4]∑eαieαjα∈[1−4]∑eαieαjeαkeαl=2δij⋅c2=2δijkl⋅c4
对于组别 3，同理可得
∑
α
∈
[
5
−
8
]
e
α
i
=
0
,
∑
α
∈
[
5
−
8
]
e
α
i
e
α
j
=
2
δ
i
j
⋅
(
2
c
)
2
∑
α
∈
[
5
−
8
]
e
α
i
e
α
j
e
α
k
=
0
,
∑
α
∈
[
5
−
8
]
e
α
i
e
α
j
e
α
k
e
α
l
=
(
4
Δ
i
j
k
l
(
4
)
−
8
δ
i
j
k
l
(
4
)
)
⋅
(
2
c
)
4
\\begin{aligned} \\sum_{\\alpha\\in[5-8]} e_{\\alpha i} &= 0 ,& \\sum_{\\alpha\\in[5-8]} e_{\\alpha i} e_{\\alpha j} &= 2\\delta_{ij} \\cdot (\\sqrt{2} c)^2\\\\ \\sum_{\\alpha\\in[5-8]} e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} &= 0 ,& \\sum_{\\alpha\\in[5-8]} e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} e_{\\alpha l} &= (4 \\Delta^{(4)}_{ijkl} – 8 \\delta^{(4)}_{ijkl}) \\cdot (\\sqrt{2} c)^4\\\\ \\end{aligned}
α∈[5−8]∑eαiα∈[5−8]∑eαieαjeαk=0,=0,α∈[5−8]∑eαieαjα∈[5−8]∑eαieαjeαkeαl=2δij⋅(2
c)2=(4Δijkl(4)−8δijkl(4))⋅(2
c)4
代入式(2-2)，依次得：
ρ
(
ω
0
+
4
ω
1
+
4
ω
5
)
+
ρ
∥
u
∥
2
[
ω
1
(
c
2
c
s
4
−
2
c
s
2
)
+
ω
5
(
2
c
2
c
s
4
−
2
c
s
2
)
−
ω
0
2
c
s
2
]
=
ρ
(2-3)
\\rho (\\omega_0 + 4 \\omega_1 + 4 \\omega_5) + \\rho \\|\\boldsymbol{u}\\|^2 \\left[ \\omega_1 \\left(\\frac{c^2}{c_s^4} – \\frac{2}{c_s^2}\\right) \\right.\\\\ \\left. + \\omega_5 \\left(\\frac{2 c^2}{c_s^4} – \\frac{2}{c_s^2}\\right) – \\frac{\\omega_0}{2 c_s^2} \\right] = \\rho \\tag{2-3}
ρ(ω0+4ω1+4ω5)+ρ∥u∥2[ω1(cs4c2−cs22)+ω5(cs42c2−cs22)−2cs2ω0]=ρ(2-3)
ρ
u
(
ω
1
2
c
2
c
s
2
+
ω
5
4
c
2
c
s
2
)
=
ρ
u
(2-4)
\\rho \\boldsymbol{u} \\left( \\omega_1 \\frac{2 c^2}{c_s^2} + \\omega_5 \\frac{4 c^2}{c_s^2} \\right) = \\rho \\boldsymbol{u} \\tag{2-4}
ρu(ω1cs22c2+ω5cs24c2)=ρu(2-4)
ρ
ω
1
∑
α
=
1
4
[
e
α
i
e
α
j
(
1
−
u
2
2
c
s
2
)
+
e
α
i
e
α
j
e
α
k
e
α
l
u
k
u
l
2
c
s
4
]
+
ρ
ω
5
∑
α
=
5
8
[
e
α
i
e
α
j
(
1
−
u
2
2
c
s
2
)
+
e
α
i
e
α
j
e
α
k
e
α
l
u
k
u
l
2
c
s
4
]
=
ρ
u
i
u
j
+
p
δ
i
j
(2-5)
\\begin{aligned} \\rho \\omega_1 \\sum_{\\alpha=1}^{4} \\left[ e_{\\alpha i} e_{\\alpha j} \\left(1 – \\frac{u^2}{2 c_s^2}\\right) + \\frac{e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} e_{\\alpha l} u_k u_l}{2 c_s^4} \\right] &+ \\\\ \\rho \\omega_5 \\sum_{\\alpha=5}^{8} \\left[ e_{\\alpha i} e_{\\alpha j} \\left(1 – \\frac{u^2}{2 c_s^2}\\right) + \\frac{e_{\\alpha i} e_{\\alpha j} e_{\\alpha k} e_{\\alpha l} u_k u_l}{2 c_s^4} \\right] &= \\rho u_i u_j + p \\delta_{ij} \\end{aligned} \\tag{2-5}
ρω1α=1∑4[eαieαj(1−2cs2u2)+2cs4eαieαjeαkeαlukul]ρω5α=5∑8[eαieαj(1−2cs2u2)+2cs4eαieαjeαkeαlukul]+=ρuiuj+pδij(2-5)
由式(2-3)，得：
ω
0
+
4
ω
1
+
4
ω
5
=
1
ω
1
(
c
2
c
s
4
−
2
c
s
2
)
+
ω
2
(
2
c
2
c
s
4
−
2
c
s
2
)
−
ω
0
2
c
s
2
=
0
\\omega_0 + 4 \\omega_1 + 4 \\omega_5 = 1 \\\\ \\omega_1 \\left(\\frac{c^2}{c_s^4} – \\frac{2}{c_s^2}\\right) + \\omega_2 \\left(\\frac{2 c^2}{c_s^4} – \\frac{2}{c_s^2}\\right) – \\frac{\\omega_0}{2 c_s^2} = 0
ω0+4ω1+4ω5=1ω1(cs4c2−cs22)+ω2(cs42c2−cs22)−2cs2ω0=0
由式(2-4)，得：
(
ω
1
+
2
ω
5
)
2
c
2
c
s
2
=
1
(\\omega_1 + 2 \\omega_5) \\frac{2 c^2}{c_s^2} =1
(ω1+2ω5)cs22c2=1
式(2-5)需要用到如下化简：
Δ
i
j
k
l
(
4
)
u
k
u
l
=
δ
i
j
(
δ
k
l
u
k
u
l
)
+
(
δ
i
k
δ
j
l
u
k
u
l
)
+
(
δ
i
l
δ
j
k
u
k
u
l
)
=
δ
i
j
∥
u
∥
2
+
2
u
i
u
j
\\begin{aligned} \\Delta^{(4)}_{ijkl} u_{k} u_{l} &= \\delta_{ij} (\\delta_{kl} u_{k} u_{l}) + (\\delta_{ik}\\delta_{jl}u_{k} u_{l}) + (\\delta_{il}\\delta_{jk} u_{k} u_{l}) \\\\ &=\\delta_{ij} \\|\\boldsymbol{u}\\|^2 + 2 u_i u_j \\end{aligned}
Δijkl(4)ukul=δij(δklukul)+(δikδjlukul)+(δilδjkukul)=δij∥u∥2+2uiuj 其中：
- 第一项仅当
k
=
l
k=l
k=l 时不为0。因此δ
k
l
u
k
u
l
=
∑
k
=
l
u
k
u
l
=
∑
k
u
k
u
k
=
∥
u
∥
2
\\delta_{kl} u_k u_l = \\sum_{k=l} u_k u_l = \\sum_{k} u_k u_k = \\|\\boldsymbol{u}\\|^2
δklukul=∑k=lukul=∑kukuk=∥u∥2。 - 第二项
δ
i
k
\\delta_{ik}
δik 要求k
=
i
k=i
k=i 时才非零，δ
j
l
\\delta_{jl}
δjl 要求l
=
j
l=j
l=j 时才非零。因此对k
,
l
k,l
k,l 求和时，只有k
=
i
k=i
k=i 且l
=
j
l=j
l=j 的项保留，即u
i
u
j
u_i u_j
uiuj。 - 第三项
δ
i
l
\\delta_{il}
δil 要求l
=
i
l=i
l=i，δ
j
k
\\delta_{jk}
δjk 要求k
=
j
k=j
k=j。因此只有l
=
i
l=i
l=i 且k
=
j
k=j
k=j 的项保留，即u
i
u
j
u_i u_j
uiuj。
由于式(2-5)的左侧可被简化为
(
ω
1
−
4
ω
5
)
c
4
c
s
4
ρ
δ
i
j
k
l
u
k
u
l
+
ρ
c
2
δ
i
j
[
(
1
−
∥
u
∥
2
2
c
s
2
)
(
2
ω
1
+
4
ω
5
)
+
2
ω
5
∥
u
∥
2
c
2
c
s
4
]
+
ρ
u
i
u
j
w
5
4
c
2
c
s
4
(2-6)
(\\omega_1 – 4 \\omega_5) \\frac{c^4}{c_s^4} \\rho \\delta_{ijkl} u_k u_l + \\rho c^2 \\delta_{ij} \\left[ (1 – \\frac{\\|\\boldsymbol{u}\\|^2}{2 c_s^2}) (2\\omega_1 + 4\\omega_5) \\right. \\\\ \\left. + 2 \\omega_5 \\frac{\\|\\boldsymbol{u}\\|^2 c^2}{c_s^4} \\right] + \\rho u_i u_j w_5 \\frac{4 c^2}{c_s^4} \\tag{2-6}
(ω1−4ω5)cs4c4ρδijklukul+ρc2δij[(1−2cs2∥u∥2)(2ω1+4ω5)+2ω5cs4∥u∥2c2]+ρuiujw5cs44c2(2-6)
所以
w
5
4
c
2
c
s
4
=
1
ω
1
−
4
ω
5
=
0
(
1
−
∥
u
∥
2
2
c
s
2
)
(
2
ω
1
+
4
ω
5
)
+
2
ω
5
∥
u
∥
2
c
2
c
s
4
=
p
ρ
c
2
\\begin{aligned} w_5 \\frac{4 c^2}{c_s^4} &= 1 \\\\ \\omega_1 – 4 \\omega_5 &= 0 \\\\ \\left(1 – \\frac{\\|\\boldsymbol{u}\\|^2}{2 c_s^2}\\right) (2\\omega_1 + 4\\omega_5) + 2 \\omega_5 \\frac{\\|\\boldsymbol{u}\\|^2 c^2}{c_s^4} &= \\frac{p}{\\rho c^2} \\end{aligned}
w5cs44c2ω1−4ω5(1−2cs2∥u∥2)(2ω1+4ω5)+2ω5cs4∥u∥2c2=1=0=ρc2p
联立这些方程就可以解得：
ω
0
=
4
/
9
\\omega_0 = 4/9
ω0=4/9,
ω
1
=
1
/
9
\\omega_1 = 1/9
ω1=1/9,
ω
5
=
1
/
36
\\omega_5 = 1/36
ω5=1/36,
c
s
=
c
/
3
c_s = c / \\sqrt{3}
cs=c/3
。
附录
（1）对于二维正M边形网格，
e
α
=
e
(
cos
(
2
π
α
M
)
,
sin
(
2
π
α
M
)
)
\\boldsymbol{e}_{\\alpha} = e \\left( \\cos(\\frac{2\\pi\\alpha}{M}), \\sin(\\frac{2\\pi\\alpha}{M})\\right)
eα=e(cos(M2πα),sin(M2πα))。此时阶数
n
n
n 与
M
M
M 的对应关系为：
|
E ( 2 ) E^{(2)} E(2) |
M > 2 M \\gt 2 M>2 |
|
E ( 3 ) E^{(3)} E(3) |
M ≥ 2 , M ≠ 3 M \\ge 2, M \\neq 3 M≥2,M=3 |
|
E ( 4 ) E^{(4)} E(4) |
M > 2 , M ≠ 4 M \\gt 2, M \\neq 4 M>2,M=4 |
|
E ( 5 ) E^{(5)} E(5) |
M ≥ 2 , M ≠ [ 3 , 5 ] M \\ge 2, M \\neq [3,5] M≥2,M=[3,5] |
|
E ( 6 ) E^{(6)} E(6) |
M > 4 , M ≠ 6 M \\gt 4, M \\neq 6 M>4,M=6 |
|
E ( 7 ) E^{(7)} E(7) |
M ≥ 2 , M ≠ [ 3 , 5 , 7 ] M \\ge 2, M \\neq [3,5,7] M≥2,M=[3,5,7] |
（2）文献[2]给出了三维情况下
E
(
n
)
E^{(n)}
E(n) 的各向同性情况，下表用 Y 和 N 分别表示是和否。其中 cyc 表示循环排列,
ϕ
=
(
1
+
5
)
/
2
\\phi=(1+\\sqrt{5})/2
ϕ=(1+5
)/2。
| 四面体 |
( 1 , 1 , 1 ) (1,1,1) (1,1,1), cyc: ( 1 , − 1 , − 1 ) (1,-1,-1) (1,−1,−1) |
4 | Y | N | N | N | N |
| 立方体 | cyc:
( ± 1 , ± 1 , ± 1 ) (\\pm 1,\\pm 1,\\pm 1) (±1,±1,±1) |
8 | Y | Y | N | Y | N |
| 八面体 | cyc:
( ± 1 , 0 , 0 ) (\\pm 1,0,0) (±1,0,0) |
6 | Y | Y | N | Y | N |
| 十二面体 |
( ± 1 , ± 1 , ± 1 ) (\\pm 1,\\pm 1,\\pm 1) (±1,±1,±1), cyc: ( 0 , ± ϕ − 1 , ± ϕ ) (0, \\pm \\phi-1, \\pm \\phi) (0,±ϕ−1,±ϕ) |
20 | Y | Y | Y | Y | N |
| 二十面体 | cyc:
( 0 , ± ϕ , ± 1 ) (0, \\pm \\phi, \\pm 1) (0,±ϕ,±1) |
12 | Y | Y | Y | Y | N |
（3）两种基本正方形格子模型的偶数阶张量
| cyc:
( ± 1 , 0 ) (\\pm 1,0) (±1,0) |
4 |
2 δ ( 2 ) 2 \\delta^{(2)} 2δ(2) |
2 δ ( 4 ) 2 \\delta^{(4)} 2δ(4) |
2 δ ( 6 ) 2 \\delta^{(6)} 2δ(6) |
|
( ± 1 , ± 1 ) (\\pm 1,\\pm 1) (±1,±1) |
8 |
4 δ ( 2 ) 4 \\delta^{(2)} 4δ(2) |
4 Δ ( 4 ) − 8 δ ( 4 ) 4 \\Delta^{(4)} – 8 \\delta^{(4)} 4Δ(4)−8δ(4) |
4 3 Δ ( 6 ) − 16 δ ( 6 ) \\frac{4}{3} \\Delta^{(6)} – 16 \\delta^{(6)} 34Δ(6)−16δ(6) |
（4）几种最对称的三维Bravais晶格的偶数阶张量
| Primitive cubic | cyc:
( ± 1 , 0 , 0 ) (\\pm 1,0,0) (±1,0,0) |
6 |
2 δ ( 2 ) 2 \\delta^{(2)} 2δ(2) |
2 δ ( 4 ) 2 \\delta^{(4)} 2δ(4) |
2 δ ( 6 ) 2 \\delta^{(6)} 2δ(6) |
| Body-centered cubic |
( ± 1 , ± 1 , ± 1 ) (\\pm 1,\\pm 1,\\pm 1) (±1,±1,±1) |
8 |
8 δ ( 2 ) 8 \\delta^{(2)} 8δ(2) |
8 ( Δ ( 4 ) − 2 δ ( 4 ) ) 8 (\\Delta^{(4)} – 2\\delta^{(4)}) 8(Δ(4)−2δ(4)) |
8 ( Δ ( 6 ) − 2 Δ ( 4 , 2 ) + 16 δ ( 6 ) ) 8 (\\Delta^{(6)} – 2 \\Delta^{(4,2)} + 16\\delta^{(6)}) 8(Δ(6)−2Δ(4,2)+16δ(6)) |
| Face-centered cubic | cyc:
( ± 1 , ± 1 , 0 ) (\\pm 1,\\pm 1,0) (±1,±1,0) |
12 |
8 δ ( 2 ) 8 \\delta^{(2)} 8δ(2) |
4 ( Δ ( 4 ) − δ ( 4 ) ) 4 (\\Delta^{(4)} – \\delta^{(4)}) 4(Δ(4)−δ(4)) |
4 ( Δ ( 4 , 2 ) + 13 δ ( 6 ) ) 4 (\\Delta^{(4,2)} + 13\\delta^{(6)}) 4(Δ(4,2)+13δ(6)) |
在文献[2]的原话中：
Δ
(
n
,
m
)
\\Delta^{(n,m)}
Δ(n,m) is the sum of all possible products of pairs of Kronecker delta symbols with
n
n
n and
m
m
m indices。
具体到
Δ
(
4
,
2
)
=
Δ
i
j
k
l
m
n
(
4
,
2
)
\\Delta^{(4,2)} = \\Delta^{(4,2)}_{ijklmn}
Δ(4,2)=Δijklmn(4,2) ，它应该是
C
6
2
=
15
\\mathrm{C}_{6}^{2} = 15
C62=15 个项的和。打个比方：第1个项是
δ
i
j
Δ
k
l
m
n
(
4
)
\\delta_{ij}\\Delta^{(4)}_{klmn}
δijΔklmn(4)；第2个项是
δ
i
k
Δ
j
l
m
n
(
4
)
\\delta_{ik}\\Delta^{(4)}_{jlmn}
δikΔjlmn(4)；后续依次类推，直到把15个组合都轮换出来，然后求和。
参考
[1] HE X, LUO L S. Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation[J]. Physical Review E, 1997, 56(6): 6811-6817. DOI:10.1103/PhysRevE.56.6811. [2] WOLFRAM S. Cellular automaton fluids 1: Basic theory[J]. Journal of Statistical Physics, 1986, 45(3-4): 471-526. DOI:10.1007/BF01021083. [3] 何雅玲, 王勇, 李庆, 童自翔. 格子Boltzmann方法的理论及应用[M]. 1 版. 北京: 高等教育出版社, 2023.
