CN104268943B - 一种基于欧拉‑拉格朗日耦合方法的流体仿真方法 - Google Patents

Abstract

本发明公开了一种基于欧拉‑拉格朗日耦合方法的流体仿真方法。其以实现高效的流体仿真为目的，结合欧拉法和拉格朗日法求解N‑S方程进行流体仿真的核心思想，在SPH和LBM的流体仿真基础上，围绕多方法耦合的计算方法展开研究。首先利用基于欧拉网格的LBM方法进行大规模流场建模，构建流体仿真的主体部分，然后利用基于拉格朗日粒子思想的SPH方法仿真诸如浪花、水珠等流体细节，并通过设计耦合算法整合LBM流体与SPH流体，实现逼真、快速、兼顾细节的大规模流体仿真。

Description

一种基于欧拉-拉格朗日耦合方法的流体仿真方法

技术领域

本发明涉及一种基于欧拉-拉格朗日耦合方法的流体仿真技术。

背景技术

目前，流体仿真朝着多学科，多角度，多平台的方向飞速发展，但是碍于现有硬件水平的限制和仿真算法的局限性，对流体模型进行精确、迅速、逼真的还原还存在诸多困难。近年来，满足工程和科研要求的真实流体仿真备受国内外学者和科研机构的关注。基于物理的真实感流体仿真研究对于从事流体相关的科研工作者和计算机图形学的研究者来说无疑面临着巨大的挑战，也具有巨大的吸引力。

基于物理的流体仿真方法主要分为两种：基于网格思想的欧拉方法和基于粒子思想的拉格朗日方法。拉格朗日方法中SPH方法多用于流体的细节彰显和模拟小规模流体，由于流体细节的精细刻画、流体仿真的速度和流体仿真规模这三者具有不可调和的矛盾，在进行精细流体细节渲染的同时很难保证实现大规模流体仿真。

LBM方法是基于欧拉思想的典型网格方法，具有天然的并行性，且边界条件容易处理，适合于大尺度、低精度、规则边界的流体仿真，如山洪地质灾害、泥石流、大规模水面等。通过精细网格划分或非规则网格划分可以拓展LBM流体仿真范围，但是由此带来的稳定性损失和效率下降将使LBM失去其固有优势。因此，与其他方法相结合，发挥LBM方法天然优势的混合方法正逐渐引起越来越多学者的关注。

为了解决上述问题，本发明提出一种基于欧拉-拉格朗日耦合方法的流体仿真技术，以实现兼顾速度、规模和细节的流体仿真为目标，围绕基于物理的流体仿真方法展开研究，首先用LBM方法进行大规模的流场建模，并利用PLSM方法提取和追踪流体表面，然后在已有LBM仿真流体的基础上，应用SPH方法刻画流体诸如水花、泡沫的细节，通过设计耦合算法将SPH流体和LBM流体集成到同一场景中，实现了兼顾细节的大规模交互场景的实时逼真绘制。

发明内容

本发明解决的技术问题是：克服了现有的流体仿真方法的效率问题，提供了一种基于欧拉-拉格朗日耦合方法的流体仿真技术，解决了流体仿真规模与细节的矛盾，满足了对流体运动过程进行实时模拟的需要。

本发明采用的技术方案为：一种基于欧拉-拉格朗日耦合方法的流体仿真方法，包括以下四个步骤：

步骤(1)、大规模流体的物理仿真：利用求解N-S方程的LBGK控制方程，采用DnQm多维离散网格模型，进行大规模流场建模，这是流体仿真的主体部分。

步骤(2)、流体自由表面追踪：应用PLSM方法构建流体自由表面。输出步骤(1)中得到的流体网格，用PLSM方法标记流体运动界面在LBM网格内部的位置，追踪流体界面，以此作为流体的近似表面。

步骤(3)、流体细节仿真：通过SPH方法来对流体细节进行物理仿真，模拟流体运动过程中产生的诸如浪花、水珠等细节。

步骤(4)、集成仿真场景、进行算法耦合。把步骤(1)、(2)、(3)中LBM流体和SPH仿真的流体集成到同一个场景中，利用耦合算法进行LBM和SPH算法的参数传递，例如密度、速度、位置等。根据耦合算法控制流体仿真过程中SPH粒子的生成时机。

本发明的原理在于：

(1)通过LBM方法模拟大规模流体运动，利用空间网格对流体的运动状态进行计算，以此发挥LBM方法适合模拟大规模、低精度流体的特点，并能在保证流体仿真规模的同时保证速度优势。

(2)为了快速精确的追踪流体的表面，本发明利用PLSM方法通过距离函数来近似代表流体液面的空间位置。该方法通过使用PLSM追踪的流体表面在网格节点中的位置，然后分析PLSM方法进行液面追踪时造成的耗散误差，通过补偿粒子的方法弥补流体液面处被错误平滑掉的流体信息。接着用快速行进算法得到流体液面，采用OpenGL对流体液面进行绘制，保持了边界的锐度并且可以直接以此作为大规模流体表面的近似代表。这样得到的流体表面，保持了良好的真实性，并具有非常高的效率。

(3)为了对流体进行逼真的细节刻画，本发明提出了一种基于LBM和SPH的耦合流体仿真方法。该方法在LBM流体的基础上，通过SPH方法模拟诸如浪花、泡沫等流体细节，并设计LBM-SPH的耦合算法控制细节生成位置和生成时机，达到良好的流体真实视觉表现。

本发明与现有技术相比的优点在于：

1、本发明提出的高效数据组织方式，流体仿真过程中的流体参数可在各个方法间无障碍的传递，保证了算法的鲁棒性。

2、高效率的网格方法，本发明利用PLSM对流体液面进行追踪，由于LBM和PLSM同属网格方法，LBM的空间网格可以被PLSM方法重用，因此有效地减少了计算时间。

3、结合不同流体仿真方法的优势，既发挥了LBM方法速度快、可仿真大规模流场的优势，又发挥了SPH方法擅长刻画流体细节的优势。

4、快速的液面追踪和绘制，利用PLSM方法可对流体液面进行追踪和补偿，提高了算法效率。

附图说明

图1基于欧拉-拉格朗日耦合方法的流体仿真流程图；

图2二维LBM示意图；

图3算法划分示意图；

图4修正网格示意图；

图5 PLSM与LBM耦合后仿真的水滴下落过程；

图6仿真结果与真实照片对比；

图7渲染后耦合方法的仿真结果。

具体实施方式

图1给出了基于欧拉-拉格朗日耦合方法的流体仿真的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于欧拉-拉格朗日耦合方法的流体仿真技术，具体实施为LBM-SPH耦合的流体仿真，主要步骤介绍如下：

1.大规模流体的物理仿真

本发明采用LBM方法来模拟大规模流体的运动，LBM由上世纪70年代提出并发展的格子气自动机(Lattice Gas Automata，LGA)发展而来，因此可把LBM看做是描述流体运动的一个物理模型。LBM方法改善了传统的LGA方法中诸如噪声、指数复杂性等问题。在LBM中，用粒子分布函数替代LGA中粒子本身进行演化，直接采用格子Boltzmann方程(LBE)作为演化方程，根据其分布函数计算流体的密度和速度。

连续的Boltzmann方程表示为：

ξ是粒子速度，Ω(f)是表示碰撞过程的碰撞因子，▽是哈密顿算子。f是连续方程的平衡态分布函数，平衡态分布函数保证了流体局部的质量与动量守恒。f线性依赖于流体密度ρ和速度u。根据上式，可通过求数值解的方法模拟流体的宏观运动。然而碰撞项Ω(f)并非分布函数的线性方程，只与分子作用力有关，因此求解起来具有难度，这使得Boltzmann方程在实际应用中受到了很大限制。改进了的Bhatnager-Gross-Krook(BGK)模型假设流体运动具有平衡态，分子之间的碰撞过程会促使分布函数f趋近于平衡态。由于碰撞导致的变化量和f与平衡态的差值成正比，即：

τ_c为松弛时间，1/τ_c称之为平均碰撞频率。

从宏观运动的角度，粒子速度ξ和分布函数f也相应的被离散到n维空间，这意味着在t时刻，某个格点x处将有沿各个离散方向的分布函数。格子Boltzmann方程(LBE)是BGK方程的特殊离散格式，从空间、时间和速度的角度将连续矢量离散成沿各个特定的方向的标量值。以二维流场为例，LBM模型示意图如图2所示

此模型中，流体运动被看作为两个步骤——碰撞(Collision)和迁移(Stream)。不含外力项的LBGK方程定义为：

f_q是粒子速度分布函数，其中q代表离散方向(0，1，2，3…，q)。是粒子的平衡态分布函数。此方程可简化为两部分，将碰撞过程和迁移过程分开表示，这样的好处在于方便计算，并使碰撞项完全线性化，用单松弛时间逼近近似表现碰撞过程，使得系统趋于平衡。

碰撞：

迁移：

在上式中，表示碰撞后的分布函数。平衡态分布函数为：

ω_q表示粒子速度方向上的权重，e_q为离散速度方向，ρ是宏观密度。c是格子声速，其值为δx/δt，与声速有关，δx和δt分别代表时间步长和网格步长。在均匀网格计算中通常取δt＝δx。在宏观上，LBM对粒子分部函数进行统计，计算出宏观压力和宏观流速等变量：

b表示离散速度数量。

此外，在LBM模型中时空并不是完全独立的，而是通过粒子的离散速度e联系起来：

δx＝e_qδt

粒子在一个时间步长内由一个节点运动到另一个节点，进行“迁移”过程，同时其他节点运动到此节点的物理量通过“碰撞”过程进行更新。

应用于流体仿真的LBM方法有许多离散速度模型，Lattice Bhatnager-Gross-Krook(LBGK)模型是迄今为止应用最为广泛的模型，本发明使用应用广泛的DnQm模型。其中n代表空间维数，m代表离散速度个数。不同模型的平衡态分布函数中权系数ω和格子声速c也不同，本发明应用D3Q15模型进行大规模流体建模，即在三维空间内，每个格子节点有15个离散的子速度。

2.基于PLSM的流体液面追踪

a)流体分区域仿真算法分配如下：

1)用LBM方法模拟流体主体，发挥LBM方法善于计算大尺度流体的优点；

2)在需要精细刻画和细节渲染的液面附近结合SPH和LBM方法，渲染流体细节，增强显示效果；

3)用水平集方法对液面进行跟踪绘制流体表面，保证流体的真实感绘制。

如图3所示，在进行大规模流体仿真的区域，相对平稳的流体占据绝大部分空间，因此采用无需刻画细节且具有速度优势的LBM方法建模；Coupling Band表示PLSM作用区域，在此区域内应用LBM和PLSM偶和算法进行流体液面跟踪；同时在液面和空气中利用SPH生成粒子表示液面诸如浪花、飞沫、气泡等细节增强真实感。

b)网格优化

流体仿真根据区域的不同，对精度的要求也不同，LBM求解区域对精度要求较低，速度较快，可用粗糙网格求解，PLSM和SPH求解区域对精度要求较高，可应用精细网格求解。基于这一思想，本文对网格进行重新划分，修正网格尺寸，在保证算法效率的同时尽可能提高算法精确度。重新划分后的网格示意图如图4图示。

粗细格点之间数据传递遵循以下方程：

上式分别表示表示精细网格的分布函数和粗糙网格的分布函数。其中，和表示粗细网格的平衡态分布函数，τ_c、τ_f分别表示粗糙网格和精细网格的松弛时间，具有如下关系：

在构建PLS区域的同时设置粒子和网格的耦合区域，通过耦合算法实现粒子仿真流体与网格仿真流体区域之间的衔接。在液面附近设置过渡区域，区域以上的流体(主要是浪花)完全采用SPH方法；区域以下的流体主要采用LBM方法进行流体仿真，由于冲撞、挤压等物理过程产生的运动至流体内部的气泡采用SPH方法渲染；在耦合区域的泡沫粒子采用SPH方法。浪花、气泡和泡沫这三种粒子通过不同的图形学绘制手段区分。根据粒子类型选择参数不同的光滑核函数，以体现不同粒子属性；渲染流体粒子；对LBM网格格点属性进行插值，交换格点和粒子属性，通过耦合方程实现两种方法的数据传递。

3.LBM方法与PLSM的耦合

PLSM方法通过追踪流体运动构建液面网格，而流体运动由LBM方法计算得到，因此需要构建耦合算法，使LBM方法求得的流体参数准确传递给PLSM方法，进而对液面进行追踪。在本文中，流体仿真算法LBM和液面跟踪技术PLSM同属欧拉网格方法，因此在参数传递、模型合并、方法耦合等方面具有优势。

为实现模型结合后算法的正确耦合，需要对LBM的算法形式加以变换，PLSM的网格尺寸与LBM网格尺寸相同，这使得LBM速度场能够直接传递给PLSM的距离函数。在本发明中，应用半拉格朗日对流项传递速度参数。半拉格朗日方法的基本思想是：已知给定的LevelSet方程和速度场u，格点x处的速度每Δt时间进行一次更新，液面位置通过Φ(x)进行跟踪：

Φ(x)＝Φ(x-u(x)Δt)

在液面跟踪过程中，需要对网格进行分类：节点处没有速度值的网格标记为空网格(空气网格)；节点处全部具有速度信息的网格标记为液体网格；部分节点具有速度信息，部分节点为空的网格标记为液面网格。其中，流体网格用LBM算法进行计算，空网格所有节点处属性置零，因为只需对液面网格的算法进行改进即可。根据距离函数标记的网格类型，气体网格标记为-1，液体网格标记为+1，因此同时具有气体网格和液体网格邻居的网格必为液面网格。液面网格的演进过程遵照LBM演进规则，碰撞过程按照如下方程进行修正：

x表示液面网格，i是液面网格相邻的气体网格方向，与i方向相反，ρ_A是密度，u是x处的速度。

算法表示的步骤为：

1)流体部分依据LBM算法模型初始化流场区域，根据LBGK方程计算分布函数；

2)时间步长Δt开始计数，各节点依次进行碰撞和演进两个过程，平衡态分布函数及分布函数同时更新；

3)根据流场内网格节点属性判断节点类型，液体节点标记为1，空节点标记为0，介于液体节点和空节点之间的节点标记为液面节点；

4)得到液面位置信息，修正碰撞后的分布函数；

5)根据修正后的分布函数计算液面节点处的速度信息，并传递到PLS算法中，重新初始化距离函数；

6)根据修正的距离函数跟踪液面新位置信息；

7)根据得到的液面位置信息判断当前格点x的类型，在下一个时间步长内反馈给LBM算法；

8)LBM算法根据当前格点类型重新初始化整个流场网格，重新初始化分布函数；

9)时间步长Δt+1，进行新一轮迭代。

耦合后的模型示意图如图5所示。

4.LBM方法与SPH方法的耦合

网格划分过程中，远离流体表面的区域划分为稀疏网格以节省资源，提高运算效率，在液面耦合区域(即PLS区域)为了提高仿真的精确性，实现细节的逼真刻画，需要划分精细的网格。区域内LBM和SPH仿真的流体之间的速度和密度传导通过以下方程实现：

Xc表示耦合位置，通过以上两个方程，确保了SPH粒子可以准确得到初始化属性信息(速度、密度)，并通过局部求解将运算结果反馈给LBM算法，实现了全局数值传递，从而达到了两种算法的耦合。

5.液面绘制

最后，运用快速行进算法(Fast Marching Mehod，FMM)对PLSM追踪出来的网格进行遍历，并用三角面片填充，SPH渲染的流体细节直接用OpenGL绘制成粒子球。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于欧拉-拉格朗日耦合方法的流体仿真方法，其特征在于包括以下5个步骤：

步骤(1)、大规模流体的物理仿真：利用求解纳维-斯托克斯方程Navier-Stokes，N-S的LatticeBhatnagar-gross-krook即LBGK控制方程，采用n维离散空间的m速度DnQm多维离散网格模型，进行大规模流场建模，这是流体仿真的主体部分；

步骤(2)、流体自由表面追踪：应用Particle Level Set Method即PLSM方法构建流体自由表面；输出步骤(1)中得到的流体网格，用PLSM方法标记流体运动界面在LatticeBoltzmann Method即LBM网格内部的位置，追踪流体界面，以此作为流体的近似表面；

步骤(3)、流体细节仿真：通过Smoothed Particle Hydrodynamics即SPH方法来对流体细节进行物理仿真，首先建立流体仿真物理模型，包括流体求解区域和边界条件，其中边界条件为PLSM追踪到的流体界面；然后设定SPH初始条件，确立粒子生成位置，初始化粒子属性信息；再然后，将控制流体运动的N-S方程离散化，用离散了的控制方程近似表达粒子属性变量及其导数形式；接下来，根据时间步长积分形式Δt开始累加，更新求解区域的粒子位置、速度、密度信息，计算粒子属性在一个时间步长内的变化量；紧接着Δt+1，计算下一时刻的速度变化量和新的位置信息，更新SPH粒子系统；最后，用新的速度和位置信息计算瞬时变量；重复以上各步骤，直到指定结束时刻和计算达到收敛条件，仿真结束；以此模拟流体运动过程中产生的浪花、水珠细节；

步骤(4)、集成仿真场景、进行算法耦合；把步骤(1)、(2)、(3)中LBM流体和SPH仿真的流体集成到同一个场景中，利用耦合算法进行LBM和SPH算法的参数传递，包括密度、速度、位置；根据耦合算法控制流体仿真过程中SPH粒子的生成时机；

步骤(5)、运用快速行进算法Fast Marching Method，FMM对PLSM追踪出来的网格进行遍历，并用三角面片填充，SPH渲染的流体细节直接用OpenGL绘制成粒子球。

2.根据权利要求1所述的基于欧拉-拉格朗日耦合方法的流体仿真方法，其特征在于：步骤(1)中所述的大规模流体的物理仿真方法，该方法使用空间场内的均匀网格表示流场，利用DnQm网格模型表示流体在网格格点处的属性，将纳维－斯托克斯Navier-Stokes，N-S方程转化为只与时间相关的离散形式的LBGK方程形式求解，具体过程为假设流体运动具有平衡态，分子之间的碰撞过程会促使流体分布函数趋近于平衡态，在流体趋于平衡态的过程中，从空间、时间和速度的角度将连续矢量离散成沿各个特定的方向的标量值，将连续的流体运动看作为两个离散步骤——碰撞(Collision)和迁移(Stream)，通过离散时间Δt驱动 流体在网格格点处的碰撞或迁移过程以实现流体模拟过程。

3.根据权利要求1所述的基于欧拉-拉格朗日耦合方法的流体仿真方法，其特征在于：步骤(3)中所述的流体细节仿真，该方法使用一组离散粒子表示流体，利用平滑核函数近似表示粒子的物理参数，将纳维－斯托克斯Navier-Stokes，N-S方程转化为只与时间相关的离散形式的微分方程，再利用积分方法求解；描述流体运动的控制方程可以由一系列拉格朗日型偏微分方程表示，这些方程遵循能量守恒定律；N-S方程可用时间全导数形式表示如下：

动量方程为：

能量方程：

其中，α，β表示坐标方向，ρ表示流体密度，t表示时间，u为流体速度，x为空间位置，σ为总应力张量，对于给定的粒子i，宏观密度可根据SPH方法用支持域内粒子密度累加得到：

其中，光滑核函数W_ij为：

W_ij＝W(|x_i-x_j|,h)＝W(R_ij,h)

R_ij表示粒子i和j之间的相对距离；r_ij/h，r_ij是两个粒子的绝对距离；h是光滑长度，直接影响SPH方法的求解精度。