对流换热仿真-主题076_格子玻尔兹曼方法-格子玻尔兹曼方法
主题076：格子玻尔兹曼方法
目录
引言
1.1 从微观到宏观的桥梁
格子玻尔兹曼方法（Lattice Boltzmann Method, LBM）是一种基于介观（mesoscopic）描述的流体模拟方法。它不同于传统的宏观方法（如有限体积法求解Navier-Stokes方程），也不同于完全微观的分子动力学模拟，而是建立了一种独特的"介观"视角。
LBM的核心思想：
- 不直接追踪单个分子的运动
- 不直接求解宏观的守恒方程
- 而是追踪粒子分布函数在离散速度空间中的演化
这种方法具有独特的优势：
- 物理直观：基于粒子运动图像，易于理解
- 计算简单：主要是局部碰撞和迁移操作，适合并行计算
- 边界处理灵活：易于处理复杂几何边界
- 多物理场耦合：自然扩展到多相流、传热等问题
1.2 LBM的历史发展
格子玻尔兹曼方法的发展经历了几个重要阶段：
早期格子气自动机（Lattice Gas Automata, LGA）：
- 1986年，Frisch、Hasslacher和Pomeau提出了FHP模型
- 使用布尔变量描述粒子存在与否
- 存在统计噪声问题
格子玻尔兹曼方法的诞生：
- 1988年，McNamara和Zanetti引入分布函数代替布尔变量
- 解决了统计噪声问题
- 奠定了现代LBM的基础
BGK近似的引入：
- 1992年，Chen和Qian等独立提出单松弛时间模型
- 大大简化了碰撞算子
- 使LBM成为实用的工程工具
多松弛时间（MRT）模型：
- 提高了数值稳定性
- 增强了可调参数的自由度
1.3 本教程内容概述
本教程将系统介绍格子玻尔兹曼方法的理论基础和实际应用：
- LBM的基本原理和数学框架
- 各种碰撞模型（BGK、MRT、TRT、MRT等）
- 边界条件的实现方法
- 多相流和热流动的模拟技术
- 5个完整的Python仿真案例，包含GIF动画
格子玻尔兹曼方法基础理论
2.1 玻尔兹曼方程
格子玻尔兹曼方法源于经典的玻尔兹曼输运方程：
∂ f ∂ t + v ⋅ ∇ f = Ω ( f ) \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f = \Omega(f) ∂t∂f+v⋅∇f=Ω(f)
其中：
- f ( x , v , t ) f(\mathbf{x}, \mathbf{v}, t) f(x,v,t)：粒子分布函数，表示在位置 x \mathbf{x} x、速度 v \mathbf{v} v、时间 t t t处找到粒子的概率密度
- v ⋅ ∇ f \mathbf{v} \cdot \nabla f v⋅∇f：对流项，描述粒子的自由运动
- Ω ( f ) \Omega(f) Ω(f)：碰撞项，描述粒子间的相互作用
宏观量的计算：
密度： ρ = ∫ f d v \rho = \int f d\mathbf{v} ρ=∫fdv
动量： ρ u = ∫ v f d v \rho \mathbf{u} = \int \mathbf{v} f d\mathbf{v} ρu=∫vfdv
能量： ρ e = 1 2 ∫ ∣ v − u ∣ 2 f d v \rho e = \frac{1}{2} \int |\mathbf{v} - \mathbf{u}|^2 f d\mathbf{v} ρe=21∫∣v−u∣2fdv
2.2 离散化
LBM的核心创新在于对速度空间的离散化：
速度空间离散化：
将连续的速度空间离散为有限个速度方向：
v ∈ { e 0 , e 1 , . . . , e Q − 1 } \mathbf{v} \in \{\mathbf{e}_0, \mathbf{e}_1, ..., \mathbf{e}_{Q-1}\} v∈{e0,e1,...,eQ−1}
常用的格子模型：
| 模型 | 维度 | 速度数 | 描述 |
|---|---|---|---|
| D1Q3 | 1D | 3 | 一维三速 |
| D2Q5 | 2D | 5 | 二维五速 |
| D2Q9 | 2D | 9 | 二维九速（最常用） |
| D3Q15 | 3D | 15 | 三维十五速 |
| D3Q19 | 3D | 19 | 三维十九速（最常用） |
| D3Q27 | 3D | 27 | 三维二十七速 |
D2Q9模型：
速度方向定义：
- e 0 = ( 0 , 0 ) \mathbf{e}_0 = (0, 0) e0=(0,0)，静止粒子
- e 1 = ( 1 , 0 ) , e 2 = ( 0 , 1 ) , e 3 = ( − 1 , 0 ) , e 4 = ( 0 , − 1 ) \mathbf{e}_1 = (1, 0), \mathbf{e}_2 = (0, 1), \mathbf{e}_3 = (-1, 0), \mathbf{e}_4 = (0, -1) e1=(1,0),e2=(0,1),e3=(−1,0),e4=(0,−1)，轴向
- e 5 = ( 1 , 1 ) , e 6 = ( − 1 , 1 ) , e 7 = ( − 1 , − 1 ) , e 8 = ( 1 , − 1 ) \mathbf{e}_5 = (1, 1), \mathbf{e}_6 = (-1, 1), \mathbf{e}_7 = (-1, -1), \mathbf{e}_8 = (1, -1) e5=(1,1),e6=(−1,1),e7=(−1,−1),e8=(1,−1)，对角
权重系数：
w i = { 4 / 9 i = 0 1 / 9 i = 1 , 2 , 3 , 4 1 / 36 i = 5 , 6 , 7 , 8 w_i = \begin{cases} 4/9 & i = 0 \\ 1/9 & i = 1,2,3,4 \\ 1/36 & i = 5,6,7,8 \end{cases} wi=⎩
⎨
⎧4/91/91/36i=0i=1,2,3,4i=5,6,7,8
2.3 格子玻尔兹曼方程
离散化后的格子玻尔兹曼方程：
f i ( x + e i Δ t , t + Δ t ) − f i ( x , t ) = Ω i ( x , t ) f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = \Omega_i(\mathbf{x}, t) fi(x+eiΔt,t+Δt)−fi(x,t)=Ωi(x,t)
或等价地写成两步：
碰撞步：
f i ∗ ( x , t ) = f i ( x , t ) + Ω i ( x , t ) f_i^*(\mathbf{x}, t) = f_i(\mathbf{x}, t) + \Omega_i(\mathbf{x}, t) fi∗(x,t)=fi(x,t)+Ωi(x,t)
迁移步：
f i ( x + e i Δ t , t + Δ t ) = f i ∗ ( x , t ) f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) = f_i^*(\mathbf{x}, t) fi(x+eiΔt,t+Δt)=fi∗(x,t)
2.4 宏观量恢复
从分布函数恢复宏观量：
密度：
ρ = ∑ i = 0 Q − 1 f i \rho = \sum_{i=0}^{Q-1} f_i ρ=i=0∑Q−1fi
速度：
ρ u = ∑ i = 0 Q − 1 e i f i \rho \mathbf{u} = \sum_{i=0}^{Q-1} \mathbf{e}_i f_i ρu=i=0∑Q−1eifi
压力：
p = c s 2 ρ p = c_s^2 \rho p=cs2ρ
其中 c s c_s cs是格子声速，对于D2Q9模型： c s = 1 / 3 c_s = 1/\sqrt{3} cs=1/3。
2.5 平衡分布函数
平衡分布函数是Maxwell-Boltzmann分布在低马赫数下的二阶展开：
f i e q = w i ρ [ 1 + e i ⋅ u c s 2 + ( e i ⋅ u ) 2 2 c s 4 − u ⋅ u 2 c s 2 ] f_i^{eq} = w_i \rho \left[1 + \frac{\mathbf{e}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{e}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2}\right] fieq=wiρ[1+cs2ei⋅u+2cs4(ei⋅u)2−2cs2u⋅u]
物理意义：
- 第一项：密度贡献
- 第二项：动量贡献（一阶）
- 第三、四项：动能贡献（二阶）
展开到二阶的重要性：
- 一阶展开只能恢复欧拉方程
- 二阶展开才能恢复Navier-Stokes方程
碰撞模型
3.1 BGK模型（单松弛时间）
Bhatnagar-Gross-Krook（BGK）模型是最简单的碰撞模型：
Ω i = − 1 τ ( f i − f i e q ) \Omega_i = -\frac{1}{\tau}(f_i - f_i^{eq}) Ωi=−τ1(fi−fieq)
其中 τ \tau τ是松弛时间，与运动粘度相关：
ν = c s 2 ( τ − 1 2 ) Δ t \nu = c_s^2 (\tau - \frac{1}{2}) \Delta t ν=cs2(τ−21)Δt
优点：
- 形式简单，计算效率高
- 物理意义清晰
缺点：
- 数值稳定性受限（ τ > 0.5 \tau > 0.5 τ>0.5）
- 普朗特数固定为1
3.2 多松弛时间（MRT）模型
MRT模型在矩空间中进行松弛：
Ω = − M − 1 S ( m − m e q ) \mathbf{\Omega} = -\mathbf{M}^{-1} \mathbf{S} (\mathbf{m} - \mathbf{m}^{eq}) Ω=−M−1S(m−meq)
其中：
- M \mathbf{M} M：变换矩阵，将分布函数投影到矩空间
- S \mathbf{S} S：对角松弛矩阵
- m = M f \mathbf{m} = \mathbf{M} \mathbf{f} m=Mf：矩向量
D2Q9的矩向量：
m = ( ρ , e , ϵ , j x , q x , j y , q y , p x x , p x y ) T \mathbf{m} = (\rho, e, \epsilon, j_x, q_x, j_y, q_y, p_{xx}, p_{xy})^T m=(ρ,e,ϵ,jx,qx,jy,qy,pxx,pxy)T
各矩的物理意义：
- ρ \rho ρ：密度
- e e e：能量
- ϵ \epsilon ϵ：能量平方
- j x , j y j_x, j_y jx,jy：动量分量
- q x , q y q_x, q_y qx,qy：热流
- p x x , p x y p_{xx}, p_{xy} pxx,pxy：应力张量分量
松弛参数的选择：
S = diag ( s 0 , s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 ) \mathbf{S} = \text{diag}(s_0, s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8) S=diag(s0,s1,s2,s3,s4,s5,s6,s7,s8)
通常：
- s 0 = s 3 = s 5 = 0 s_0 = s_3 = s_5 = 0 s0=s3=s5=0（守恒量）
- s 7 = s 8 = 1 / τ s_7 = s_8 = 1/\tau s7=s8=1/τ（粘度相关）
- 其他参数调节稳定性
3.3 双松弛时间（TRT）模型
TRT模型是对称和反对称分量的分别松弛：
f i + = f i + f i ˉ 2 , f i − = f i − f i ˉ 2 f_i^+ = \frac{f_i + f_{\bar{i}}}{2}, \quad f_i^- = \frac{f_i - f_{\bar{i}}}{2} fi+=2fi+fiˉ,fi−=2fi−fiˉ
其中 i ˉ \bar{i} iˉ是 i i i的反向方向。
碰撞：
f i ∗ = f i − 1 τ + ( f i + − f i e q + ) − 1 τ − ( f i − − f i e q − ) f_i^* = f_i - \frac{1}{\tau^+}(f_i^+ - f_i^{eq+}) - \frac{1}{\tau^-}(f_i^- - f_i^{eq-}) fi∗=fi−τ+1(fi+−fieq+)−τ−1(fi−−fieq−)
优点：
- 比MRT计算量小
- 可以独立调节粘度和边界精度
3.4 正则化（Regularized）模型
正则化模型通过投影消除非物理的高阶矩：
f i ∗ = f i e q + ( 1 − 1 τ ) ( f i n e q − f i ( 1 ) ) f_i^* = f_i^{eq} + (1 - \frac{1}{\tau})(f_i^{neq} - f_i^{(1)}) fi∗=fieq+(1−τ1)(fineq−fi(1))
其中 f i ( 1 ) f_i^{(1)} fi(1)是通过Chapman-Enskog展开得到的一阶非平衡项。
优点：
- 提高数值稳定性
- 减少网格依赖性
边界条件
4.1 无滑移边界（壁面）
反弹格式（Bounce-Back）：
最简单的边界处理方法：
f i ˉ ( x w , t + Δ t ) = f i ∗ ( x w , t ) f_{\bar{i}}(\mathbf{x}_w, t + \Delta t) = f_i^*(\mathbf{x}_w, t) fiˉ(xw,t+Δt)=fi∗(xw,t)
其中 x w \mathbf{x}_w xw是壁面节点， i ˉ \bar{i} iˉ是 i i i的反向方向。
特性：
- 实现简单
- 二阶精度
- 壁面位置在格点中间
改进的反弹格式：
对于运动壁面：
f i ˉ = f i − 2 w i ρ c s 2 e i ⋅ u w f_{\bar{i}} = f_i - \frac{2 w_i \rho}{c_s^2} \mathbf{e}_i \cdot \mathbf{u}_w fiˉ=fi−cs22wiρei⋅uw
4.2 速度边界（入口）
Zou-He格式：
基于反弹思想的非平衡反弹格式：
- 已知速度 u i n \mathbf{u}_{in} uin
- 计算密度 ρ i n \rho_{in} ρin
- 计算平衡分布 f i e q f_i^{eq} fieq
- 非平衡部分反弹
优点：
- 严格满足边界条件
- 数值稳定性好
4.3 压力边界（出口）
密度/压力给定：
给定密度 ρ o u t \rho_{out} ρout，计算速度：
u o u t = − 1 + f 0 + f 3 + f 4 + 2 ( f 1 + f 5 + f 8 ) ρ o u t u_{out} = -1 + \frac{f_0 + f_3 + f_4 + 2(f_1 + f_5 + f_8)}{\rho_{out}} uout=−1+ρoutf0+f3+f4+2(f1+f5+f8)
4.4 周期性边界
对于进出口周期性流动：
f i ( x i n ) = f i ( x o u t − e i Δ t ) f_i(\mathbf{x}_{in}) = f_i(\mathbf{x}_{out} - \mathbf{e}_i \Delta t) fi(xin)=fi(xout−eiΔt)
4.5 复杂几何处理
浸没边界法（IBM）：
将边界力作为体积力加入LBE：
f i ( x + e i Δ t , t + Δ t ) − f i ( x , t ) = − 1 τ ( f i − f i e q ) + F i f_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\tau}(f_i - f_i^{eq}) + F_i fi(x+eiΔt,t+Δt)−fi(x,t)=−τ1(fi−fieq)+Fi
曲边界处理：
使用插值或浸入边界方法处理非对齐网格的曲面。
多相流模拟
5.1 颜色梯度模型（Color Gradient Model）
Gunstensen等人提出的原始多相LBM：
两种分布函数：
- f i R f_i^R fiR：红色流体
- f i B f_i^B fiB：蓝色流体
碰撞步骤：
- 分别对两种流体进行BGK碰撞
- 添加颜色梯度引起的扰动
- 重新着色（recoloring）
颜色梯度：
G ( x ) = ∑ i e i [ ρ R ( x + e i ) − ρ B ( x + e i ) ] \mathbf{G}(\mathbf{x}) = \sum_i \mathbf{e}_i [\rho_R(\mathbf{x} + \mathbf{e}_i) - \rho_B(\mathbf{x} + \mathbf{e}_i)] G(x)=i∑ei[ρR(x+ei)−ρB(x+ei)]
5.2 伪势模型（Pseudo-Potential Model）
Shan-Chen模型是最流行的多相LBM：
相互作用力：
F ( x ) = − G ψ ( x ) ∑ i w i ψ ( x + e i ) e i \mathbf{F}(\mathbf{x}) = -G \psi(\mathbf{x}) \sum_i w_i \psi(\mathbf{x} + \mathbf{e}_i) \mathbf{e}_i F(x)=−Gψ(x)i∑wiψ(x+ei)ei
其中：
- G G G：相互作用强度
- ψ ( ρ ) \psi(\rho) ψ(ρ)：有效密度（伪势）
常用的伪势形式：
ψ ( ρ ) = ρ 0 [ 1 − exp ( − ρ / ρ 0 ) ] \psi(\rho) = \rho_0 \left[1 - \exp(-\rho/\rho_0)\right] ψ(ρ)=ρ0[1−exp(−ρ/ρ0)]
力项的处理：
通过速度修正实现：
u e q = u + τ F ρ \mathbf{u}^{eq} = \mathbf{u} + \frac{\tau \mathbf{F}}{\rho} ueq=u+ρτF
5.3 自由能模型（Free Energy Model）
基于热力学自由能的相场LBM：
Cahn-Hilliard方程：
∂ ϕ ∂ t + ∇ ⋅ ( ϕ u ) = ∇ ⋅ ( M ∇ μ ) \frac{\partial \phi}{\partial t} + \nabla \cdot (\phi \mathbf{u}) = \nabla \cdot (M \nabla \mu) ∂t∂ϕ+∇⋅(ϕu)=∇⋅(M∇μ)
化学势：
μ = δ F δ ϕ \mu = \frac{\delta F}{\delta \phi} μ=δϕδF
压力张量：
P α β = [ p − κ ϕ ∇ 2 ϕ − κ 2 ∣ ∇ ϕ ∣ 2 ] δ α β + κ ∂ α ϕ ∂ β ϕ P_{\alpha\beta} = \left[p - \kappa \phi \nabla^2 \phi - \frac{\kappa}{2}|\nabla \phi|^2\right] \delta_{\alpha\beta} + \kappa \partial_\alpha \phi \partial_\beta \phi Pαβ=[p−κϕ∇2ϕ−2κ∣∇ϕ∣2]δαβ+κ∂αϕ∂βϕ
5.4 多相流的界面处理
表面张力：
在伪势模型中自然出现表面张力：
σ ∼ ∫ − ∞ ∞ ( p N − p T ) d n \sigma \sim \int_{-\infty}^{\infty} \left(p_N - p_T\right) dn σ∼∫−∞∞(pN−pT)dn
接触角处理：
通过调节壁面处的相互作用强度实现不同的润湿性。
热格子玻尔兹曼方法
6.1 双分布函数法
最常用的热流动模拟方法：
流体分布函数 f i f_i fi：描述质量、动量守恒
温度分布函数 g i g_i gi：描述能量守恒
温度演化方程：
g i ( x + e i Δ t , t + Δ t ) − g i ( x , t ) = − 1 τ g ( g i − g i e q ) g_i(\mathbf{x} + \mathbf{e}_i \Delta t, t + \Delta t) - g_i(\mathbf{x}, t) = -\frac{1}{\tau_g}(g_i - g_i^{eq}) gi(x+eiΔt,t+Δt)−gi(x,t)=−τg1(gi−gieq)
温度平衡分布：
g i e q = w i T [ 1 + e i ⋅ u c s 2 ] g_i^{eq} = w_i T \left[1 + \frac{\mathbf{e}_i \cdot \mathbf{u}}{c_s^2}\right] gieq=wiT[1+cs2ei⋅u]
宏观温度：
T = ∑ i g i T = \sum_i g_i T=i∑gi
热扩散系数：
α = c s 2 ( τ g − 1 2 ) Δ t \alpha = c_s^2 (\tau_g - \frac{1}{2}) \Delta t α=cs2(τg−21)Δt
6.2 Boussinesq近似
自然对流模拟：
浮力项：
F = − ρ g β ( T − T 0 ) \mathbf{F} = -\rho \mathbf{g} \beta (T - T_0) F=−ρgβ(T−T0)
其中：
- g \mathbf{g} g：重力加速度
- β \beta β：热膨胀系数
- T 0 T_0 T0：参考温度
实现方法：
将浮力作为体积力加入流体分布函数的碰撞步。
6.3 温度依赖的物性
对于大温差问题，需要考虑粘度、热导率的温度依赖性：
变松弛时间：
τ ( T ) = τ 0 + τ 1 ( T − T 0 ) \tau(T) = \tau_0 + \tau_1 (T - T_0) τ(T)=τ0+τ1(T−T0)
Sutherland定律：
μ ( T ) = μ 0 ( T T 0 ) 3 / 2 T 0 + S T + S \mu(T) = \mu_0 \left(\frac{T}{T_0}\right)^{3/2} \frac{T_0 + S}{T + S} μ(T)=μ0(T0T)3/2T+ST0+S
Python仿真案例
案例1：泊肃叶流动
本案例模拟二维管道中的层流流动。
"""
案例1：泊肃叶流动
模拟二维管道中的层流流动
"""
import numpy as np
import matplotlib.pyplot as plt
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')
# 网格参数
Nx, Ny = 200, 41
lx, ly = Nx-1, Ny-1
# D2Q9参数
c = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
[1, 1], [-1, 1], [-1, -1], [1, -1]])
w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
cs2 = 1/3
# 流动参数
Re = 100
u_max = 0.1
nu = u_max * ly / Re
tau = 3 * nu + 0.5
print("案例1：泊肃叶流动")
print(f"Re = {Re}, u_max = {u_max}, tau = {tau:.3f}")
# 初始化
f = np.ones((9, Nx, Ny))
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
def equilibrium(rho, ux, uy):
"""计算平衡分布函数"""
feq = np.zeros((9, Nx, Ny))
uu = ux**2 + uy**2
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
feq[i] = w[i] * rho * (1 + 3*eu + 4.5*eu**2 - 1.5*uu)
return feq
# 主循环
n_steps = 10000
for step in range(n_steps):
# 碰撞
feq = equilibrium(rho, ux, uy)
f = f - (f - feq) / tau
# 迁移
f_streamed = np.zeros_like(f)
for i in range(9):
f_streamed[i] = np.roll(np.roll(f[i], c[i, 0], axis=0), c[i, 1], axis=1)
f = f_streamed
# 边界条件：入口（Zou-He）
ux[0, 1:-1] = u_max * (1 - ((np.arange(1, Ny-1) - ly/2) / (ly/2))**2)
uy[0, 1:-1] = 0
rho[0, 1:-1] = 1 / (1 - ux[0, 1:-1]) * (
f[0, 0, 1:-1] + f[2, 0, 1:-1] + f[4, 0, 1:-1] +
2*(f[3, 0, 1:-1] + f[6, 0, 1:-1] + f[7, 0, 1:-1])
)
# 边界条件：出口（零梯度）
f[:, -1, :] = f[:, -2, :]
# 边界条件：壁面（反弹）
for i in range(9):
if i == 0: continue
f[i, :, 0] = f[[0,3,4,1,2,7,8,5,6][i], :, 0]
f[i, :, -1] = f[[0,3,4,1,2,7,8,5,6][i], :, -1]
# 计算宏观量
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
if step % 1000 == 0:
print(f"Step {step}/{n_steps}")
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# 速度场
ax = axes[0]
speed = np.sqrt(ux**2 + uy**2)
im = ax.imshow(speed.T, origin='lower', cmap='RdYlBu_r', aspect='auto')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度幅值')
plt.colorbar(im, ax=ax)
# 速度矢量
ax = axes[1]
X, Y = np.meshgrid(range(0, Nx, 5), range(0, Ny, 2))
ax.quiver(X, Y, ux[::5, ::2].T, uy[::5, ::2].T, scale=5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度矢量')
ax.set_aspect('equal')
# 中心线速度分布
ax = axes[2]
y_center = np.arange(Ny)
ax.plot(ux[Nx//2, :], y_center, 'b-o', markersize=4, label='LBM')
# 理论解
y_theory = np.linspace(0, ly, 100)
u_theory = u_max * (1 - ((y_theory - ly/2) / (ly/2))**2)
ax.plot(u_theory, y_theory, 'r--', linewidth=2, label='理论解')
ax.set_xlabel('u')
ax.set_ylabel('y')
ax.set_title('中心线速度分布')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题076_格子玻尔兹曼方法\\仿真结果\\case1_poiseuille.png', dpi=150, bbox_inches='tight')
plt.close()
print("案例1完成！")
案例2：圆柱绕流
本案例模拟圆柱绕流，观察卡门涡街。
"""
案例2：圆柱绕流
模拟圆柱绕流，观察卡门涡街现象
"""
import numpy as np
import matplotlib.pyplot as plt
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')
# 网格参数
Nx, Ny = 400, 100
lx, ly = Nx-1, Ny-1
# D2Q9参数
c = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
[1, 1], [-1, 1], [-1, -1], [1, -1]])
w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
# 流动参数
Re = 200
u_max = 0.1
D = 20 # 圆柱直径
cy = Ny // 2
cx = Ny
nu = u_max * D / Re
tau = 3 * nu + 0.5
print("案例2：圆柱绕流")
print(f"Re = {Re}, D = {D}, tau = {tau:.3f}")
# 圆柱掩码
def make_cylinder_mask():
mask = np.zeros((Nx, Ny), dtype=bool)
for i in range(Nx):
for j in range(Ny):
if (i - cx)**2 + (j - cy)**2 <= (D/2)**2:
mask[i, j] = True
return mask
cylinder_mask = make_cylinder_mask()
# 初始化
f = np.ones((9, Nx, Ny))
rho = np.sum(f, axis=0)
ux = np.zeros((Nx, Ny))
uy = np.zeros((Nx, Ny))
ux[0, :] = u_max
def equilibrium(rho, ux, uy):
feq = np.zeros((9, Nx, Ny))
uu = ux**2 + uy**2
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
feq[i] = w[i] * rho * (1 + 3*eu + 4.5*eu**2 - 1.5*uu)
return feq
# 主循环
n_steps = 50000
save_interval = 500
vorticity_history = []
print("开始计算...")
for step in range(n_steps):
# 碰撞
feq = equilibrium(rho, ux, uy)
f = f - (f - feq) / tau
# 迁移
f_streamed = np.zeros_like(f)
for i in range(9):
f_streamed[i] = np.roll(np.roll(f[i], c[i, 0], axis=0), c[i, 1], axis=1)
f = f_streamed
# 入口：均匀流
ux[0, :] = u_max
uy[0, :] = 0
rho[0, :] = 1
feq_in = equilibrium(rho[0, :], ux[0, :], uy[0, :])
f[:, 0, :] = feq_in
# 出口：零梯度
f[:, -1, :] = f[:, -2, :]
# 上下壁面：反弹
for i in range(9):
if i == 0: continue
f[i, :, 0] = f[[0,3,4,1,2,7,8,5,6][i], :, 0]
f[i, :, -1] = f[[0,3,4,1,2,7,8,5,6][i], :, -1]
# 圆柱表面：反弹
for i in range(9):
if i == 0: continue
f[i, cylinder_mask] = f[[0,3,4,1,2,7,8,5,6][i], cylinder_mask]
# 计算宏观量
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
# 圆柱内部速度置零
ux[cylinder_mask] = 0
uy[cylinder_mask] = 0
if step % save_interval == 0 and step > 20000:
# 计算涡量
vorticity = np.roll(uy, -1, axis=0) - np.roll(uy, 1, axis=0) \
- np.roll(ux, -1, axis=1) + np.roll(ux, 1, axis=1)
vorticity_history.append(vorticity.copy())
print(f"Step {step}/{n_steps}, 保存涡量场")
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 速度场
ax = axes[0, 0]
speed = np.sqrt(ux**2 + uy**2)
im = ax.imshow(speed.T, origin='lower', cmap='RdYlBu_r', aspect='auto', vmin=0, vmax=u_max*1.5)
ax.add_patch(plt.Circle((cx, cy), D/2, color='black', fill=True))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度幅值')
plt.colorbar(im, ax=ax)
# 涡量场
ax = axes[0, 1]
vorticity = np.roll(uy, -1, axis=0) - np.roll(uy, 1, axis=0) \
- np.roll(ux, -1, axis=1) + np.roll(ux, 1, axis=1)
im = ax.imshow(vorticity.T, origin='lower', cmap='RdBu_r', aspect='auto', vmin=-0.1, vmax=0.1)
ax.add_patch(plt.Circle((cx, cy), D/2, color='black', fill=True))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('涡量')
plt.colorbar(im, ax=ax)
# 压力场
ax = axes[1, 0]
p = rho / 3
im = ax.imshow(p.T, origin='lower', cmap='RdYlBu_r', aspect='auto')
ax.add_patch(plt.Circle((cx, cy), D/2, color='black', fill=True))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('压力')
plt.colorbar(im, ax=ax)
# 速度矢量
ax = axes[1, 1]
X, Y = np.meshgrid(range(0, Nx, 10), range(0, Ny, 5))
ax.quiver(X, Y, ux[::10, ::5].T, uy[::10, ::5].T, scale=10)
ax.add_patch(plt.Circle((cx, cy), D/2, color='black', fill=True))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度矢量')
ax.set_aspect('equal')
ax.set_xlim(0, Nx)
ax.set_ylim(0, Ny)
plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题076_格子玻尔兹曼方法\\仿真结果\\case2_cylinder_flow.png', dpi=150, bbox_inches='tight')
plt.close()
print("案例2完成！")
案例3：瑞利-贝纳德对流
本案例模拟底部加热引起的自然对流。
"""
案例3：瑞利-贝纳德对流
模拟底部加热引起的自然对流
"""
import numpy as np
import matplotlib.pyplot as plt
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')
# 网格参数
Nx, Ny = 100, 50
# D2Q9参数
c = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
[1, 1], [-1, 1], [-1, -1], [1, -1]])
w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
# 物理参数
Ra = 10000 # 瑞利数
Pr = 0.71 # 普朗特数
# 计算松弛时间
g = 0.001 # 重力加速度
beta = 0.01 # 热膨胀系数
nu = np.sqrt(Pr * g * beta * (1.0 - 0.0) * Ny**3 / Ra)
alpha = nu / Pr
tau_nu = 3 * nu + 0.5
tau_alpha = 3 * alpha + 0.5
print("案例3：瑞利-贝纳德对流")
print(f"Ra = {Ra}, Pr = {Pr}")
print(f"nu = {nu:.6f}, alpha = {alpha:.6f}")
print(f"tau_nu = {tau_nu:.3f}, tau_alpha = {tau_alpha:.3f}")
# 初始化
f = np.ones((9, Nx, Ny)) # 流体分布函数
g = np.zeros((9, Nx, Ny)) # 温度分布函数
g[0] = 0.5 # 初始温度
rho = np.sum(f, axis=0)
ux = np.zeros((Nx, Ny))
uy = np.zeros((Nx, Ny))
T = np.sum(g, axis=0)
def equilibrium_f(rho, ux, uy):
feq = np.zeros((9, Nx, Ny))
uu = ux**2 + uy**2
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
feq[i] = w[i] * rho * (1 + 3*eu + 4.5*eu**2 - 1.5*uu)
return feq
def equilibrium_g(T, ux, uy):
geq = np.zeros((9, Nx, Ny))
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
geq[i] = w[i] * T * (1 + 3*eu)
return geq
# 主循环
n_steps = 100000
for step in range(n_steps):
# 计算宏观量
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
T = np.sum(g, axis=0)
# 添加浮力（Boussinesq近似）
Fy = -g * beta * (T - 0.5)
uy += tau_nu * Fy / rho
# 流体碰撞
feq = equilibrium_f(rho, ux, uy)
f = f - (f - feq) / tau_nu
# 温度碰撞
geq = equilibrium_g(T, ux, uy)
g = g - (g - geq) / tau_alpha
# 流体迁移
f_streamed = np.zeros_like(f)
for i in range(9):
f_streamed[i] = np.roll(np.roll(f[i], c[i, 0], axis=0), c[i, 1], axis=1)
f = f_streamed
# 温度迁移
g_streamed = np.zeros_like(g)
for i in range(9):
g_streamed[i] = np.roll(np.roll(g[i], c[i, 0], axis=0), c[i, 1], axis=1)
g = g_streamed
# 温度边界条件
g[:, :, 0] = equilibrium_g(1.0, ux[:, 0], uy[:, 0]) # 底部：T = 1
g[:, :, -1] = equilibrium_g(0.0, ux[:, -1], uy[:, -1]) # 顶部：T = 0
# 左右周期性边界
g[:, 0, :] = g[:, -1, :]
# 壁面反弹
for i in range(9):
if i == 0: continue
f[i, :, 0] = f[[0,3,4,1,2,7,8,5,6][i], :, 0]
f[i, :, -1] = f[[0,3,4,1,2,7,8,5,6][i], :, -1]
if step % 10000 == 0:
print(f"Step {step}/{n_steps}")
# 最终计算
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
T = np.sum(g, axis=0)
# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 温度场
ax = axes[0, 0]
im = ax.imshow(T.T, origin='lower', cmap='RdYlBu_r', aspect='auto')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('温度场')
plt.colorbar(im, ax=ax)
# 速度矢量
ax = axes[0, 1]
X, Y = np.meshgrid(range(0, Nx, 4), range(0, Ny, 2))
ax.quiver(X, Y, ux[::4, ::2].T, uy[::4, ::2].T, scale=2)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度矢量')
ax.set_aspect('equal')
# 流线图
ax = axes[1, 0]
stream = np.zeros((Nx, Ny))
for i in range(1, Nx):
stream[i, :] = stream[i-1, :] + uy[i, :]
ax.contour(stream.T, levels=20, colors='black', linewidths=0.5)
im = ax.imshow(T.T, origin='lower', cmap='RdYlBu_r', aspect='alpha=0.7')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('流线图')
# 温度剖面
ax = axes[1, 1]
ax.plot(T[Nx//2, :], range(Ny), 'b-', linewidth=2, label='中心线')
ax.plot(np.mean(T, axis=0), range(Ny), 'r--', linewidth=2, label='平均')
ax.set_xlabel('温度')
ax.set_ylabel('y')
ax.set_title('温度剖面')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题076_格子玻尔兹曼方法\\仿真结果\\case3_rayleigh_benard.png', dpi=150, bbox_inches='tight')
plt.close()
print("案例3完成！")
案例4：液滴撞击液膜的伪势模型
本案例演示使用Shan-Chen伪势模型模拟多相流。
"""
案例4：液滴撞击液膜的伪势模型
使用Shan-Chen模型模拟多相流动
"""
import numpy as np
import matplotlib.pyplot as plt
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')
# 网格参数
Nx, Ny = 200, 100
# D2Q9参数
c = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
[1, 1], [-1, 1], [-1, -1], [1, -1]])
w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
# 伪势参数
G = -5.0 # 相互作用强度（负值表示吸引）
rho_l = 2.0 # 液相密度
rho_g = 0.3 # 气相密度
tau = 1.0
print("案例4：液滴撞击液膜的伪势模型")
print(f"G = {G}, rho_l = {rho_l}, rho_g = {rho_g}")
def psi(rho):
"""伪势函数"""
return rho_0 * (1 - np.exp(-rho/rho_0))
rho_0 = 1.0
def equilibrium(rho, ux, uy):
feq = np.zeros((9, Nx, Ny))
uu = ux**2 + uy**2
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
feq[i] = w[i] * rho * (1 + 3*eu + 4.5*eu**2 - 1.5*uu)
return feq
# 初始化：液滴+液膜
rho = np.ones((Nx, Ny)) * rho_g
# 底部液膜
rho[:, :20] = rho_l
# 液滴
drop_x, drop_y = Nx//2, 70
drop_r = 15
for i in range(Nx):
for j in range(Ny):
if (i - drop_x)**2 + (j - drop_y)**2 < drop_r**2:
rho[i, j] = rho_l
# 初始化分布函数
ux = np.zeros((Nx, Ny))
uy = np.zeros((Nx, Ny))
uy[:, :20] = 0 # 液膜静止
f = equilibrium(rho, ux, uy)
# 主循环
n_steps = 5000
for step in range(n_steps):
# 计算宏观量
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
# 计算相互作用力（Shan-Chen）
psi_val = psi(rho)
Fx = np.zeros((Nx, Ny))
Fy = np.zeros((Nx, Ny))
for i in range(1, 9):
psi_neighbor = np.roll(np.roll(psi_val, c[i, 0], axis=0), c[i, 1], axis=1)
Fx += G * w[i] * c[i, 0] * psi_neighbor
Fy += G * w[i] * c[i, 1] * psi_neighbor
Fx *= psi_val
Fy *= psi_val
# 添加重力
Fy -= 0.001 * rho
# 速度修正
ux += tau * Fx / rho
uy += tau * Fy / rho
# 碰撞
feq = equilibrium(rho, ux, uy)
f = f - (f - feq) / tau
# 迁移
f_streamed = np.zeros_like(f)
for i in range(9):
f_streamed[i] = np.roll(np.roll(f[i], c[i, 0], axis=0), c[i, 1], axis=1)
f = f_streamed
# 边界条件
for i in range(9):
if i == 0: continue
f[i, :, 0] = f[[0,3,4,1,2,7,8,5,6][i], :, 0]
f[i, :, -1] = f[[0,3,4,1,2,7,8,5,6][i], :, -1]
if step % 500 == 0:
print(f"Step {step}/{n_steps}")
# 最终计算
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# 密度场
ax = axes[0]
im = ax.imshow(rho.T, origin='lower', cmap='RdYlBu_r', aspect='auto', vmin=rho_g, vmax=rho_l)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('密度场')
plt.colorbar(im, ax=ax)
# 速度场
ax = axes[1]
speed = np.sqrt(ux**2 + uy**2)
im = ax.imshow(speed.T, origin='lower', cmap='viridis', aspect='auto')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度幅值')
plt.colorbar(im, ax=ax)
# 速度矢量
ax = axes[2]
X, Y = np.meshgrid(range(0, Nx, 5), range(0, Ny, 3))
ax.quiver(X, Y, ux[::5, ::3].T, uy[::5, ::3].T, scale=5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('速度矢量')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题076_格子玻尔兹曼方法\\仿真结果\\case4_droplet_impact.png', dpi=150, bbox_inches='tight')
plt.close()
print("案例4完成！")
案例5：顶盖驱动方腔流的GIF动画
本案例模拟顶盖驱动方腔流，生成GIF动画。
"""
案例5：顶盖驱动方腔流的GIF动画
模拟经典CFD测试案例：顶盖驱动方腔流
"""
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')
# 网格参数
Nx, Ny = 100, 100
# D2Q9参数
c = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
[1, 1], [-1, 1], [-1, -1], [1, -1]])
w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
# 流动参数
Re = 1000
u_lid = 0.1
nu = u_lid * Nx / Re
tau = 3 * nu + 0.5
print("案例5：顶盖驱动方腔流的GIF动画")
print(f"Re = {Re}, u_lid = {u_lid}, tau = {tau:.3f}")
def equilibrium(rho, ux, uy):
feq = np.zeros((9, Nx, Ny))
uu = ux**2 + uy**2
for i in range(9):
eu = c[i, 0] * ux + c[i, 1] * uy
feq[i] = w[i] * rho * (1 + 3*eu + 4.5*eu**2 - 1.5*uu)
return feq
# 初始化
f = np.ones((9, Nx, Ny))
rho = np.sum(f, axis=0)
ux = np.zeros((Nx, Ny))
uy = np.zeros((Nx, Ny))
# 存储帧
frames = []
save_interval = 2000
n_steps = 60000
print("生成动画帧...")
for step in range(n_steps):
# 计算宏观量
rho = np.sum(f, axis=0)
ux = np.sum(f * c[:, 0].reshape(9, 1, 1), axis=0) / rho
uy = np.sum(f * c[:, 1].reshape(9, 1, 1), axis=0) / rho
# 顶盖速度
ux[:, -1] = u_lid
uy[:, -1] = 0
# 碰撞
feq = equilibrium(rho, ux, uy)
f = f - (f - feq) / tau
# 迁移
f_streamed = np.zeros_like(f)
for i in range(9):
f_streamed[i] = np.roll(np.roll(f[i], c[i, 0], axis=0), c[i, 1], axis=1)
f = f_streamed
# 边界条件：反弹
for i in range(9):
if i == 0: continue
# 底边
f[i, :, 0] = f[[0,3,4,1,2,7,8,5,6][i], :, 0]
# 左边
f[i, 0, :] = f[[0,3,4,1,2,7,8,5,6][i], 0, :]
# 右边
f[i, -1, :] = f[[0,3,4,1,2,7,8,5,6][i], -1, :]
# 顶盖（Zou-He）
rho[:, -1] = 1 / (1 + u_lid) * (
f[0, :, -1] + f[1, :, -1] + f[3, :, -1] +
2*(f[2, :, -1] + f[5, :, -1] + f[6, :, -1])
)
f[4, :, -1] = f[2, :, -1] - 2/3 * rho[:, -1] * u_lid
f[7, :, -1] = f[5, :, -1] + 0.5*(f[1, :, -1] - f[3, :, -1]) - 1/6 * rho[:, -1] * u_lid
f[8, :, -1] = f[6, :, -1] - 0.5*(f[1, :, -1] - f[3, :, -1]) - 1/6 * rho[:, -1] * u_lid
# 保存帧
if step % save_interval == 0 and step > 10000:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 速度场
ax = axes[0]
speed = np.sqrt(ux**2 + uy**2)
im = ax.imshow(speed.T, origin='lower', cmap='RdYlBu_r',
aspect='auto', vmin=0, vmax=u_lid)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'速度场 (step={step})')
plt.colorbar(im, ax=ax)
# 涡量
ax = axes[1]
vorticity = np.roll(uy, -1, axis=0) - np.roll(uy, 1, axis=0) \
- np.roll(ux, -1, axis=1) + np.roll(ux, 1, axis=1)
im = ax.imshow(vorticity.T, origin='lower', cmap='RdBu_r',
aspect='auto', vmin=-0.1, vmax=0.1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('涡量')
plt.colorbar(im, ax=ax)
plt.tight_layout()
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=80, bbox_inches='tight')
buf.seek(0)
frames.append(Image.open(buf))
plt.close()
print(f"帧 {len(frames)}/25, step={step}")
# 保存GIF
print("保存GIF动画...")
frames[0].save('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题076_格子玻尔兹曼方法\\仿真结果\\case5_lid_driven_cavity.gif',
save_all=True, append_images=frames[1:], duration=200, loop=0)
print("案例5完成！")
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
0
0- 0
已为社区贡献264条内容
所有评论(0)