本文同步发表于知乎 和CSDN
我在阅读其他文献时,看到有文献提到 Kupershtokh 的一篇关于 LBM 外力项的文章,他提出的格式在其他文章中又被称为 Exact difference method(EDM) 。
本想搜索这篇文章的PDF进行阅读,但最后只找到一个网站 介绍了部分信息。这份粗略笔记的内容,都是基于其他文献中对该方法的介绍而整理的。
介绍
Boltzmann-BGK方程的近似
密度分布函数 f f f 的 Boltzmann-BGK方程写作:
∂ f ∂ t + ξ ⋅ ∇ f + a ⋅ ∇ ξ f = − 1 τ f [ f − f ( e q ) ] (1) \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}
∂ t ∂ f + ξ ⋅ ∇ f + a ⋅ ∇ ξ f = − τ f 1 [ f − f ( e q ) ] ( 1 )
其中:
f ( e q ) = ρ ( 2 π R g T ) D / 2 exp ( − ( ξ − u ) 2 2 R g T ) (2) 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}
f ( e q ) = ( 2 π R g T ) D / 2 ρ exp ( − 2 R g T ( ξ − u ) 2 ) ( 2 )
符号含义分别为
ξ \bold{\xi} ξ : 粒子速度;
ρ \rho ρ : 流体密度;
u \bold{u} u : 流场宏观速度;
R g R_g R g : 气体常数;
T T T : 温度;
a = d u / d t \bold{a}=\mathrm{d}\bold{u}/\mathrm{d}t a = d u / d t : 流场加速度。
当Knudsen数较小 时,Boltzmann方程的体力项可以基于f ( e q ) f^{(eq)} f ( e q ) 进行近似,即
∇ ξ f ≈ ∇ ξ f ( e q ) \nabla_{\bold{\xi}} f \approx \nabla_{\bold{\xi}} f^{(eq)}
∇ ξ f ≈ ∇ ξ f ( e q )
根据式(2), f ( e q ) f^{(eq)} f ( e q ) 是分子随机速度 C = ξ − u \bold{C} = \bold{\xi} - \bold{u} C = ξ − u 的指数函数,则:
∂ f ( e q ) ∂ u = − ∂ f ( e q ) ∂ ξ = − ∇ ξ f ( e q ) \frac{\partial f^{(eq)}}{\partial \bold{u}} = -\frac{\partial f^{(eq)}}{\partial \bold{\xi}} = -\nabla_{\bold{\xi}} f^{(eq)}
∂ u ∂ f ( e q ) = − ∂ ξ ∂ f ( e q ) = − ∇ ξ f ( e q )
若假设密度 ρ \rho ρ 为常数 ,则 d ρ d t = 0 \frac{\mathrm{d} \rho}{\mathrm{d} t}=0 d t d ρ = 0 ,那么 f ( e q ) f^{(eq)} f ( e q ) 的全导数为:
d f ( e q ) d t = ∂ f ( e q ) ∂ u ⋅ d u d t + ∂ f ( e q ) ∂ ρ ⋅ d ρ d t ≈ ∂ f ( e q ) ∂ u ⋅ d u d t = − a ⋅ ∇ ξ f ( e q ) \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}
d t d f ( e q ) = ∂ u ∂ f ( e q ) ⋅ d t d u + ∂ ρ ∂ f ( e q ) ⋅ d t d ρ ≈ ∂ u ∂ f ( e q ) ⋅ d t d u = − a ⋅ ∇ ξ f ( e q )
也就是说,Kupershtokh 等近似的 Boltzmann-BGK方程写作:
∂ f ∂ t + ξ ⋅ ∇ f = − 1 τ f [ f − f ( e q ) ] + d f ( e q ) d t (3) \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}
∂ t ∂ f + ξ ⋅ ∇ f = − τ f 1 [ f − f ( e q ) ] + d t d f ( e q ) ( 3 )
近似方程的离散化
将式(3)在[ t , t + δ t ] [t,t+\delta_t] [ t , t + δ t ] 内沿着e α \bold{e}_\alpha e α 方向积分,可得
f α ( x + e α δ t , t + δ t ) − f α ( x , t ) = − 1 τ [ f α ( x , t ) − f α ( e q ) ( ρ , u ∗ ) ] + [ f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) ] (4) \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}
f α ( x + e α δ t , t + δ t ) − f α ( x , t ) = − τ 1 [ f α ( x , t ) − f α ( e q ) ( ρ , u ∗ ) ] + [ f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) ] ( 4 )
其中:
δ t \delta_t δ t 为时间步长;
τ = τ f / δ t \tau = \tau_f / \delta_t τ = τ f / δ t 为无量纲松弛时间;
x \bold{x} x 为网格点位置,ρ , u ∗ \rho, \bold{u^*} ρ , u ∗ 分别为该点在 t t t 时刻的密度和速度(与宏观速度 u \bold{u} u 有区别);
f α , f α ( e q ) f_\alpha, f_\alpha^{(eq)} f α , f α ( e q ) 分别为 e α \bold{e}_\alpha e α 方向的分布函数 和平衡态 ;
f α ( e q ) = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + ( e α ⋅ u ) 2 2 c s 4 − ∣ u ∣ 2 2 c s 2 ] = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + 1 2 c s 4 ( u u ) : ( e α e α − c s 2 [ I ] ) ] (4-1) \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}
f α ( e q ) = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 ( e α ⋅ u ) 2 − 2 c s 2 ∣ u ∣ 2 ] = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 1 ( u u ) : ( e α e α − c s 2 [ I ] ) ] ( 4 - 1 )
这里 w α w_\alpha w α 为 e α \bold{e}_\alpha e α 的权重系数, [ I ] [\bold{I}] [ I ] 为单位矩阵;
两个同样为m*n大小的矩阵 [ A ] [\bold{A}] [ A ] 和 [ B ] [\bold{B}] [ B ] ,则 [ A ] : [ B ] = ∑ i = 1 m ∑ j = 1 n A i j B i j \displaystyle [\bold{A}] : [\bold{B}] = \sum^{m}_{i=1}\sum^{n}_{j=1} A_{ij} B_{ij} [ A ] : [ B ] = i = 1 ∑ m j = 1 ∑ n A i j B i j ;
F α F_\alpha F α 为外力项,其表达式基于外力 F = ρ d u ∗ d t = ρ a \bold{F} = \rho \dfrac{\mathrm{d} \bold{u^*}}{\mathrm{d} t} = \rho \bold{a} F = ρ d t d u ∗ = ρ a 进行构建。
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) F_\alpha = f_\alpha^{(eq)} (\bold{u^*}(t+\delta_t)) - f_\alpha^{(eq)} (\bold{u^*}(t))
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) )
并且
u ∗ ( t ) = ( ∑ α f α e α ) / ( ∑ α f α ) u ∗ ( t + δ t ) = u ∗ ( t ) + ∫ t t + δ t d u ∗ d t d t = u ∗ ( t ) + ∫ t t + δ t F ρ d t ≈ u ∗ ( t ) + F ρ δ t (5) \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}
u ∗ ( t ) u ∗ ( t + δ t ) = ( α ∑ f α e α ) / ( α ∑ f α ) = u ∗ ( t ) + ∫ t t + δ t d t d u ∗ d t = u ∗ ( t ) + ∫ t t + δ t ρ F d t ≈ u ∗ ( t ) + ρ F δ t ( 5 )
这意味着加速度 a = F / ρ \bold{a}=\bold{F}/\rho a = F / ρ 在 [ t , t + δ t ] [t,t+\delta_t] [ t , t + δ t ] 时间段内被视为常数。
Kupershtokh等在这篇文章说明,在使用该源项时,流体宏观量的计算应表示为:
ρ = ∑ α f α , ρ u = ∑ i e α f α + δ t F 2 (6) \rho = \sum_{\alpha} f_\alpha ,\quad \rho \bold{u} = \sum_{i} \bold{e}_\alpha f_\alpha + \frac{\delta_t \bold{F}}{2}
\tag{6}
ρ = α ∑ f α , ρ u = i ∑ e α f α + 2 δ t F ( 6 )
Li 等通过 Chapman-Enskog 展开证明:该格式向宏观方程引入的误差项为:
E r r o r E D M = − δ t 2 ∇ ⋅ ( F F 4 ρ ) (7) \mathbf{Error}_{\mathrm{EDM}} = -\delta_t^2 \nabla\cdot\left( \frac{\bold{FF}}{4 \rho} \right)
\tag{7}
E r r o r E D M = − δ t 2 ∇ ⋅ ( 4 ρ F F ) ( 7 )
EDM 源项的展开式
为便于分析,这里采用如下书写形式的平衡态:
f α ( e q ) = w α ρ ⋅ [ 1 + e α ⋅ u c s 2 + 1 2 c s 4 ( u u ) : ( e α e α − c s 2 [ I ] ) ] 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]
f α ( e q ) = w α ρ ⋅ [ 1 + c s 2 e α ⋅ u + 2 c s 4 1 ( u u ) : ( e α e α − c s 2 [ I ] ) ]
将其带入 EDM 源项,记 Δ u = F ρ δ t \Delta\bold{u} = \frac{\bold{F}}{\rho} \delta_t Δ u = ρ F δ t ,整理得:
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) = ρ w α { e α ⋅ Δ u c s 2 + 1 2 c s 4 [ ( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ ] : ( e α e α − c s 2 [ I ] ) } \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}
F α = f α ( e q ) ( u ∗ ( t + δ t ) ) − f α ( e q ) ( u ∗ ( t ) ) = ρ w α { c s 2 e α ⋅ Δ u + 2 c s 4 1 [ ( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ ] : ( e α e α − c s 2 [ I ] ) }
这里涉及矩阵 ( u ∗ + Δ u ) ( u ∗ + Δ u ) (\bold{u}^* + \Delta\bold{u})(\bold{u}^* + \Delta\bold{u}) ( u ∗ + Δ u ) ( u ∗ + Δ u ) 和 u ∗ u ∗ \bold{u}^* \bold{u}^* u ∗ u ∗ 的相减,由于:
( u ∗ + Δ u ) ( u ∗ + Δ u ) = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) , u ∗ u ∗ = u i ∗ u j ∗ (\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 ∗ + Δ u ) ( u ∗ + Δ u ) = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) , u ∗ u ∗ = u i ∗ u j ∗
其中 u i ∗ u_i^* u i ∗ 和 Δ u i = F i δ t / ρ \Delta u_i = F_i \delta_{t} / \rho Δ u i = F i δ t / ρ 分别为 u ∗ \bold{u}^* u ∗ 和 Δ u \Delta\bold{u} Δ u 在 i 轴上的分量,所以:
( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) − u i ∗ u j ∗ = δ t ρ ( u i ∗ F j + u j ∗ F i ) + δ t 2 ρ 2 F i F j = δ t ρ ( u i ∗ + F i δ t 2 ρ ) F j + δ t ρ ( u j ∗ + F j δ t 2 ρ ) F 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}
( u ∗ + Δ u ) ( u ∗ + Δ u ) − u ∗ u ∗ = ( u i ∗ + Δ u i ) ( u j ∗ + Δ u j ) − u i ∗ u j ∗ = ρ δ t ( u i ∗ F j + u j ∗ F i ) + ρ 2 δ t 2 F i F j = ρ δ t ( u i ∗ + 2 ρ F i δ t ) F j + ρ δ t ( u j ∗ + 2 ρ F j δ t ) F i
这里令 v E D M = u ∗ + F 2 ρ δ t \bold{v}_{\rm{EDM}} = \bold{u}^{*} + \frac{\bold{F}}{2 \rho} \delta_{t} v E D M = u ∗ + 2 ρ F δ t ,则 EDM 源项的展开式化简为:
F α = δ t w α [ e α ⋅ F c s 2 + 1 2 c s 4 ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ] (8) 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}
F α = δ t w α [ c s 2 e α ⋅ F + 2 c s 4 1 ( v E D M F + F v E D M ) : ( e α e α − c s 2 [ I ] ) ] ( 8 )
上式即为 EDM 源项的展开结果。