CN103366045A - 基于格子Boltzmann的流体可视化仿真方法 - Google Patents

Abstract

本发明涉及流体实时模拟仿真技术领域，具体涉及一种基于格子Boltzmann的流体可视化仿真方法，对流体速度进行归一化处理，依据特征值划分区间，对不同区间内的数据施加不同的权，用三种颜色强度及其混合值分别表示不同流速的流体，与传统的加权绘制方法相比，本发明减少了计算量，提高了仿真的实时性，同时保证了仿真的精度和稳定性，实现了清晰准确的流体仿真动画，且绘制结果显示可清晰、准确的表现流体的动力学特性。

Description

基于格子Boltzmann的流体可视化仿真方法

技术领域

本发明涉及流体实时模拟仿真技术领域，具体涉及一种基于格子Boltzmann的流体可视化仿真方法。 

背景技术

在科研和工程领域中，流体的计算和可视化具有十分广泛的应用，如何生动、实时的对流体进行可视化是计算机图形学中的一个重要问题。 

由于真实流体速度变化范围大，且流体速度为矢量值，难以直接转化为可用于计算机绘制图像的像素值。 

E.W.Llewellin在文献《LBflow:An extensible lattice Boltzmann framework for the simulation of geophysical》中提出的LBM（Lattice Boltzmann Method）仿真实现，具有参数可控，可扩展性等优点，但其计算过程和可视化过程分别在两个程序中实现，只能先进行数值计算，然后通过读入计算结果进行可视化显示，不能对流体实时数据进行可视化显示，而且没有给出流体可视化的绘制方法。 

Jose Guilherme Mayworm和Sicilia Ferreira Judice等在文献《Computational System for Visualization and Lattice Boltzmann Fluid Simulation》开发了一个用于流体数据计算的可视化模拟系统，可实现流体动画的实时动态显示，其主要是介绍基于Lattice Boltzmann的流体计算方法，没有介绍其流体可视化方法，而且其流体可视化效果不佳，流体颜色过渡突兀，成像效果模糊，难以辨别低速流场。 

传统的流体绘制思路是首先由流体速度求得流体速率，即速度的模，然后乘以权系数使之达到图形学绘制范围，进而通过不同像素值的颜色信息表示流体速 度的快慢，但是该方法存在一个显著的缺点：即权系数的选择难以把握，例如，对于流体速度变化范围大的同一流场，选择小权系数可以使高速流场绘制清晰，但低速流场处理后的数据难以达到像素值要求，颜色显示就会不明显，选择大权系数可以保证低速流场绘制准确，但是容易造成高速流场内数据溢出，易造成高速流场绘制失真。因此，一种准确、简单的流体可视化方法的提出具有十分重要的应用价值。 

发明内容

本发明的目的是通过可视化方法表征流体的动力学特性，解决流体可视化的实时性和逼真性绘制问题，克服了现有技术的不足。 

本发明的技术方案是： 

一种基于格子Boltzmann的流体可视化仿真方法，包括以下步骤： 

步骤1：根据格子Boltzmann方法构建流体力学模型，选择格子Boltzmann方程初始化分布函数，所述格子Boltzmann方程初始化分布函数为： 

$f_{\mathrm{\α}}(r+e_{\mathrm{\α}}{\mathrm{\δ}}_t,t+{\mathrm{\δ}}_t)-f_{\mathrm{\α}}(r,t)=-\frac{1}{\mathrm{\τ}}[f_{\mathrm{\α}}(r,t)-f_{\mathrm{\α}}^{\mathrm{eq}}(r,t)]$

其中τ定义为 $\mathrm{\τ}=\frac{v}{c_s^2\·{\mathrm{\δ}}_t}+\frac{1}{2}$

式中格子声速c_s取

，δt为时间步长，

(r,t)为平衡态分布函数，r为空间位置矢量，t为时间； 

ν为运动粘度系数，其与流场尺寸和流体速度成正比，与雷诺系数Re成反比，即： $v=\frac{\mathrm{Re}}{L\·U};$

式中：Ｌ为方腔的高度，Ｕ为顶盖流驱动速度； 

步骤2：设置初始值，即：流体初始速度、雷诺系数、流体密度和流场尺寸，对 流场进行网格划分； 

步骤3：设定区间[0.0,1.0]作为值域A，选取0.25、0.5、0.75作为特征值将值域A划分为[0.0,0.25]、[0.25,0.5]、[0.5,0.75]和[0.75,1.0]四个区间； 

步骤4：用步骤3中得到的四个区间分别表示蓝色、绿色、红色和黄色的强度值，即：区间[0.0,0.25]表示蓝色值强度，其定义为区间B；区间[0.25,0.5]表示绿色值强度，其定义为区间G；区间[0.5,0.75]表示红色值强度，其定义为区间R；区间[0.75,1.0]表示黄色值强度，其定义为区间RG； 

步骤5：步骤1中的格子Boltzmann方程初始化分布函数向前演进一个时间步长，图像处理单元接收流场内的节点数据，即读取速度信息，所述速度信息包含x方向的速度分量和y方向的速度分量； 

步骤6：将步骤5中所读取的速度信息转化为无量纲数值信息，并进行归一化处理，所述归一化处理所采用的公式为： 

$u_x=\frac{1}{\mathrm{\ρ}}\underset{\mathrm{\α}}{\mathrm{\Σ}}{f}_{\mathrm{\α}}/e_{\mathrm{\α}}|$

$u_y=\frac{1}{\mathrm{\ρ}}\underset{\mathrm{\α}}{\mathrm{\Σ}}{f}_{\mathrm{\α}}/e_{\mathrm{\α}}|$

$m=\frac{(\sqrt{{u_x}^2+{u_y}^2}-u_{\min })}{u_{\max }-u_{\min }}$

其中u_x和u_y代表速度在坐标轴x、y方向上的分量；ρ代表流体宏观密度；f_α代表局部平衡态分布函数；α代表离散速度方向；e_α代表各个方向上的离散速度；u_max为流场内流体速率的最大值；u_min为流场内流体速率的最小值；m表示归一化处理后的数据值； 

步骤7：将步骤6中归一化处理所得到的数据值m与步骤4中得到的四个区间作 比较，如果m∈B，设定权系数w_B∈[0.0,4.0]，w_G=w_R=0；如果m∈G，设定权系数w_B∈[0,2.0]，w_G∈[0.0,2.0]，w_R=0；如果m∈R，设定权系数w_B=0，w_G∈[0,4/3]，w_R∈[0,4/3]；如果m∈RG，设定权系数w_B=0，w_G∈[0,1.0]，w_R=1.0； 

步骤8：基于OpenGL中的RGB颜色显示模式确定R、G、B的强度值： 

glColor3f(floatr,floatg,floatb) 

其中，r=m×w_R，g=m×w_G，b=m×w_B，利用此颜色配比方案，绘制边长为0.5个像素的正方形，表示当前格点处流体速度值： 

glBegin(GL_POLYGON); 

    glVertex2f(i-0.5,j-0.5); 

    glVertex2f(i-0.5,j+0.5); 

    glVertex2f(i+0.5,j+0.5); 

    glVertex2f(i+0.5,j-0.5); 

glEnd(); 

其中i、j表示当前格点的坐标； 

步骤9：遍历整个流场，绘制出整个流场内节点的颜色信息； 

步骤10：重复步骤5至步骤9，直到达到收敛条件，停止图像绘制此时图像即为流场稳定后的成像效果。 

步骤5中所述的时间步长为流场节点每进行一次演进的时间，具体指程序中屏幕每刷新一帧就进行一次演进的过程，所述的一帧即为一个时间步长。 

步骤10中所述的收敛条件为： 

$\mathrm{Error}=\frac{\sqrt{{\mathrm{\&Sigma;}}_{i,j}\left\{\right[u_x(i,j,t+{\mathrm{\&delta;}}_t)-u_x{(i,j,t)]}^2+[u_y(i,j,t+{\mathrm{\&delta;}}_t)-u_y{(i,j,t)]}^2\}}}{\sqrt{{\mathrm{\&Sigma;}}_{i,j}[u_x{(i,j,t+{\mathrm{\&delta;}}_t)}^2+u_y{(i,j,t+{\mathrm{\&delta;}}_t)}^2]}}<\mathrm{\&epsiv;}$

其中u_x和u_y为速度在x、y方向上的分量；δ_t为时间步长；ε为一个无穷小量， 当error小于这个量时，近似认为流场稳定，程序收敛。 

本发明的有益积极效果是： 

1.本发明提出的基于格子Boltzmann的流体可视化仿真方法，操作简便，适用范围广泛，可以对流速不同、种类不同的多种流体进行有效的可视化仿真。 

2.本发明对流体速度进行归一化处理，依据特征值划分区间，对不同区间内的数据施加不同的权，用三种颜色强度及其混合值分别表示不同流速的流体，与传统的加权绘制方法相比，本发明减少了计算量，提高了仿真的实时性，同时保证了仿真的精度和稳定性，实现了清晰准确的流体仿真动画，且绘制结果显示可清晰、准确的表现流体的动力学特性。 

附图说明

图1为本发明的流程示意图； 

图2为通过本发明绘制的流体可视化仿真效果图； 

图3为通过现有技术绘制的流体可视化仿真效果图。 

具体实施方式

参见图1至图3所示，一种基于格子Boltzmann的流体可视化仿真方法，包括以下步骤： 

步骤1：根据格子Boltzmann方法构建流体力学模型，选择格子Boltzmann方程初始化分布函数，所述格子Boltzmann方程初始化分布函数为： 

$f_{\mathrm{\&alpha;}}(r+e_{\mathrm{\&alpha;}}{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t)-f_{\mathrm{\&alpha;}}(r,t)=-\frac{1}{\mathrm{\&tau;}}[f_{\mathrm{\&alpha;}}(r,t)-f_{\mathrm{\&alpha;}}^{\mathrm{eq}}(r,t)]$

其中τ定义为： $\mathrm{\&tau;}=\frac{v}{c_s^2\&CenterDot;{\mathrm{\&delta;}}_t}+\frac{1}{2}$

式中格子声速c_s取，δt为时间步长，

(r,t)为平衡态分布函数，r为空间位置矢量，t为时间； 

ν为运动粘度系数，其与流场尺寸和流体速度成正比，与雷诺系数Re成反比，即： $v=\frac{\mathrm{Re}}{L\&CenterDot;U};$

式中：Ｌ为方腔的高度，Ｕ为顶盖流驱动速度； 

步骤2：设置初始值，即：流体初始速度、雷诺系数、流体密度和流场尺寸，对流场进行网格划分； 

步骤3：设定区间[0.0,1.0]作为值域A，选取0.25、0.5、0.75作为特征值将值域A划分为[0.0,0.25]、[0.25,0.5]、[0.5,0.75]和[0.75,1.0]四个区间； 

步骤4：用步骤3中得到的四个区间分别表示蓝色、绿色、红色和黄色的强度值，即：区间[0.0,0.25]表示蓝色值强度，其定义为区间B；区间[0.25,0.5]表示绿色值强度，其定义为区间G；区间[0.5,0.75]表示红色值强度，其定义为区间R；区间[0.75,1.0]表示黄色值强度，其定义为区间RG； 

步骤5：步骤1中的格子Boltzmann方程初始化分布函数向前演进一个时间步长，图像处理单元接收流场内的节点数据，即读取速度信息，所述速度信息包含x方向的速度分量和y方向的速度分量； 

步骤6：将步骤5中所读取的速度信息转化为无量纲数值信息，并进行归一化处理，所述归一化处理所采用的公式为： 

$u_x=\frac{1}{\mathrm{\&rho;}}\underset{\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{f}_{\mathrm{\&alpha;}}/e_{\mathrm{\&alpha;}}|$

$u_y=\frac{1}{\mathrm{\&rho;}}\underset{\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{f}_{\mathrm{\&alpha;}}/e_{\mathrm{\&alpha;}}|$

$m=\frac{(\sqrt{{u_x}^2+{u_y}^2}-u_{\min })}{u_{\max }-u_{\min }}$

其中u_x和u_y代表速度在坐标轴x、y方向上的分量；ρ代表流体宏观密度；f_α代 表局部平衡态分布函数；α代表离散速度方向；e_α代表各个方向上的离散速度；u_max为流场内流体速率的最大值；u_min为流场内流体速率的最小值；m表示归一化处理后的数据值； 

步骤7：将步骤6中归一化处理所得到的数据值m与步骤4中得到的四个区间作比较，如果m∈B，设定权系数w_B∈[0.0,4.0]，w_G=w_R=0；如果m∈G，设定权系数w_B∈[0,2.0]，w_G∈[0.0,2.0]，w_R=0；如果m∈R，设定权系数w_B=0，w_G∈[0,4/3]，w_R∈[0,4/3]；如果m∈RG，设定权系数w_B=0，w_G∈[0,1.0]，w_R=1.0； 

步骤8：基于OpenGL中的RGB颜色显示模式确定R、G、B的强度值： 

glColor3f(floatr,floatg,floatb) 

其中，r=m×w_R，g=m×w_G，b=m×w_B，利用此颜色配比方案，绘制边长为0.5个像素的正方形，表示当前格点处流体速度值： 

glBegin(GL_POLYGON); 

    glVertex2f(i-0.5,j-0.5); 

    glVertex2f(i-0.5,j+0.5); 

    glVertex2f(i+0.5,j+0.5); 

    glVertex2f(i+0.5,j-0.5); 

glEnd(); 

其中i、j表示当前格点的坐标； 

步骤9：遍历整个流场，绘制出整个流场内节点的颜色信息； 

步骤10：重复步骤5至步骤9，直到达到收敛条件，停止图像绘制此时图像即为流场稳定后的成像效果。 

步骤5中所述的时间步长为流场节点每进行一次演进的时间，具体指程序中屏幕每刷新一帧就进行一次演进的过程，所述的一帧即为一个时间步长。 

步骤10中所述的收敛条件为： 

$\mathrm{Error}=\frac{\sqrt{{\mathrm{\&Sigma;}}_{i,j}\left\{\right[u_x(i,j,t+{\mathrm{\&delta;}}_t)-u_x{(i,j,t)]}^2+[u_y(i,j,t+{\mathrm{\&delta;}}_t)-u_y{(i,j,t)]}^2\}}}{\sqrt{{\mathrm{\&Sigma;}}_{i,j}[u_x{(i,j,t+{\mathrm{\&delta;}}_t)}^2+u_y{(i,j,t+{\mathrm{\&delta;}}_t)}^2]}}<\mathrm{\&epsiv;}$

其中u_x和u_y为速度在x、y方向上的分量；δ_t为时间步长；ε为一个无穷小量，当error小于这个量时，近似认为流场稳定，程序收敛。 

其中附图3为文献《ComputationalSystemforVisualizationandLattice BoltzmannFluidSimulation》绘制的流体可视化仿真效果图，通过与本发明所绘制的流体可视化仿真效果图相比较，可直观的看出本发明的成像效果更为清洗，颜色过渡平滑，实现了流体可视化并提高了仿真的实时性，同时保证了仿真的精度和稳定性，可广泛应用于工程施工、地址灾害仿真、科学实验等领域的流体仿真和分析。 

Claims (3)

1.一种基于格子Boltzmann的流体可视化仿真方法，其特征在于：包括以下步骤：步骤1：根据格子Boltzmann方法构建流体力学模型，选择格子Boltzmann方程初始化分布函数，所述格子Boltzmann方程初始化分布函数为：

$f_{\mathrm{\&alpha;}}(r+e_{\mathrm{\&alpha;}}{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t)-f_{\mathrm{\&alpha;}}(r,t)=-\frac{1}{\mathrm{\&tau;}}[f_{\mathrm{\&alpha;}}(r,t)-f_{\mathrm{\&alpha;}}^{\mathrm{eq}}(r,t)]$

其中τ定义为 $\mathrm{\&tau;}=\frac{v}{c_s^2\&CenterDot;{\mathrm{\&delta;}}_t}+\frac{1}{2}$

式中格子声速c_s取

δ_t为时间步长，

为平衡态分布函数，r为空间位置矢量，t为时间；

ν为运动粘度系数，其与流场尺寸和流体速度成正比，与雷诺系数Re成反比，即： $v=\frac{\mathrm{Re}}{L\&CenterDot;U};$

式中：Ｌ为方腔的高度，Ｕ为顶盖流驱动速度；

步骤2：设置初始值，即：流体初始速度、雷诺系数、流体密度和流场尺寸，对流场进行网格划分；

步骤3：设定区间[0.0,1.0]作为值域A，选取0.25、0.5、0.75作为特征值将值域A划分为[0.0,0.25]、[0.25,0.5]、[0.5,0.75]和[0.75,1.0]四个区间；

步骤4：用步骤3中得到的四个区间分别表示蓝色、绿色、红色和黄色的强度值，即：区间[0.0,0.25]表示蓝色值强度，其定义为区间B；区间[0.25,0.5]表示绿色值强度，其定义为区间G；区间[0.5,0.75]表示红色值强度，其定义为区间R；区间[0.75,1.0]表示黄色值强度，其定义为区间RG；

步骤5：步骤1中的格子Boltzmann方程初始化分布函数向前演进一个时间步长，图像处理单元接收流场内的节点数据，即读取速度信息，所述速度信息包含x方向的速度分量和y方向的速度分量；

步骤6：将步骤5中所读取的速度信息转化为无量纲数值信息，并进行归一化处理，所述归一化处理所采用的公式为：

$u_x=\frac{1}{\mathrm{\&rho;}}\underset{\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{f}_{\mathrm{\&alpha;}}\left|e_{\mathrm{\&alpha;}}\right|$

$u_y=\frac{1}{\mathrm{\&rho;}}\underset{\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{f}_{\mathrm{\&alpha;}}\left|e_{\mathrm{\&alpha;}}\right|$

$m=\frac{(\sqrt{{u_x}^2+{u_y}^2}-u_{\min })}{u_{\max }-u_{\min }}$

其中u_x和u_y代表速度在坐标轴x、y方向上的分量；ρ代表流体宏观密度；f_α代表局部平衡态分布函数；α代表离散速度方向；e_α代表各个方向上的离散速度；u_max为流场内流体速率的最大值；u_min为流场内流体速率的最小值；m表示归一化处理后的数据值；

步骤7：将步骤6中归一化处理所得到的数据值m与步骤4中得到的四个区间作比较，如果m∈B，设定权系数w_B∈[0.0,4.0]，w_G=w_R=0；如果m∈G，设定权系数w_B∈[0,2.0]，w_G∈[0.0,2.0]，w_R=0；如果m∈R，设定权系数w_B=0，w_G∈[0,4/3]，w_R∈[0,4/3]；如果m∈RG，设定权系数w_B=0，w_G∈[0,1.0]，w_R=1.0；

步骤8：基于OpenGL中的RGB颜色显示模式确定R、G、B的强度值：

glColor3f(floatr,floatg,floatb)

其中，r=m×w_R，g=m×w_G，b=m×w_B，利用此颜色配比方案，绘制边长为0.5个像素的正方形，表示当前格点处流体速度值：

glBegin(GL_POLYGON);

glVertex2f(i-0.5,j-0.5);

glVertex2f(i-0.5,j+0.5);

glVertex2f(i+0.5,j+0.5);

glVertex2f(i+0.5,j-0.5);

glEnd();

其中i、j表示当前格点的坐标；

步骤9：遍历整个流场，绘制出整个流场内节点的颜色信息；

步骤10：重复步骤5至步骤9，直到达到收敛条件，停止图像绘制此时图像即为流场稳定后的成像效果。

2.根据权利要求1所述的基于格子Boltzmann的流体可视化仿真方法，其特征在于：步骤5中所述的时间步长为流场节点每进行一次演进的时间，具体指程序中屏幕每刷新一帧就进行一次演进的过程，所述的一帧即为一个时间步长。

3.根据权利要求1所述的基于格子Boltzmann的流体可视化仿真方法，其特征在于：步骤10中所述的收敛条件为：

$\mathrm{Error}=\frac{\sqrt{{\mathrm{\&Sigma;}}_{i,j}\{{[u_x(i,j,t+{\mathrm{\&delta;}}_t)-u_x(i,j,t)]}^2+{[u_y(i,j,t+{\mathrm{\&delta;}}_t)-u_y(i,j,t)]}^2\}}}{\sqrt{{\mathrm{\&Sigma;}}_{i,j}[u_x{(i,j,t+{\mathrm{\&delta;}}_t)}^2+u_y{(i,j,t+{\mathrm{\&delta;}}_t)}^2]}}<\mathrm{\&epsiv;}$

其中u_x和u_y为速度在x、y方向上的分量；δ_t为时间步长；ε为一个无穷小量，当error小于这个量时，近似认为流场稳定，程序收敛。

Images