一种经gpu加速的浸没边界-格子玻尔兹曼流固耦合模拟方法
【技术领域】
[0001] 本发明属于计算流体动力学和计算固体力学的交叉领域，特别设及高性能计算在 上述领域的应用。
【背景技术】
[0002] 流体和可变形的结构之间的相互作用在自然界普遍存在，特别在生物组织和器官 中更是随处可见，典型的例子包括昆虫翅膀、鱼罐、人类屯、脏瓣膜和声带。虽然上述结构属 于解创学或生理学领域的研究对象，但是他们在=维空间的运动学特征和随时变化的空间 形状对生物体完成指定的功能至关重要。因为结构所经历的剧烈变形和周围流场的内在复 杂性，所W迄今为止正确模拟=维条件下的流固禪合问题仍然是学术界的一大挑战之一。
[0003] 传统的用来模拟流固禪合问题的方法通常基于适体网格，典型的例子包括任意拉 格朗日一欧拉法和有限元法。上述方法的劣势在于底层流体网格需要根据结构位置来进行 调整，而且为了防止网格的剧烈变形需要进行特殊处理。总的来说网格生成是代价高昂的 过程，特别是当结构形状复杂时为了更新网格拓扑结构需要大量耗时的计算，使得将传统 方法推向实际应用存在困难。
[0004] 因此，非常有必要针对效率问题开发一种新的流固禪合模拟方法，能够在不W降 低计算精度为代价的前提下，缩短计算时间，为实际工程提供可靠的指导建议。
[0005] 浸没边界法作为一种专口模拟流固禪合问题的数值算法自提出W来已经得到了 广泛和成功的应用。在此法中，流体和结构体分别采用欧拉法和拉格朗日法描述，且分别用 来离散流体和结构体的正交笛卡尔网格和非结构化网格没有必要重合在一起，而该正是浸 没边界法相比于传统方法的一大优势，解决了每一时间步都需要根据结构体位置和形态重 新构造网格的难点。流体和结构体的相互作用是通过狄拉克函数实现的，利用此函数能够 将定义在结构体网格节点上的作用力扩散至周边流体网格点，并且也能够将流体速度插值 到结构体网格点上，从而更新结构体的位置。除了上述步骤外，无需额外的附加处理即可非 常简洁地完成流固禪合的模拟。
[0006] 在浸没边界法中，最基本的一点在于如何求解粘性不可压缩NS方程。该方程的求 解方法非常丰富，比如说快速傅里叶变换、投影方法或者是有限体积法。此外新近发展起来 的非常受欢迎的格子玻尔兹曼方法也是很好的选择。与传统方法求解的NS方法不一样，采 用介观视角的格子玻尔兹曼方法求解的是分布函数所满足的玻尔兹曼方程。在过去的十年 间，格子玻尔兹曼方法已经被确认为一种非常强大的计算流体动力学工具，并且在其分支 领域都已得到了成功应用。和传统方法比起来，格子玻尔兹曼方法的易于并行W及额外物 理效应的方便引入等优点得到越来越多学者的关注。
[0007] 自从2003年W来，半导体工业的微处理器设计开始朝两个主要方向发展。首先多 核路径依靠强化串行代码的执行效率和多核技术提高性能。多核开始于两核模型，并且随 着每一次更新换代核的数量加倍，比如Intel的Core i7处理器有四核。与此相反，众核路 径则将关注的重点放在了提高并行代码执行的吞吐量上面，和中央处理器比较，众核处理 器将包含数量更多但是小得多的核屯、，比如NVIDIA公司的Ge化rceGTX280图形处理器包 含240个核。作为比较，2009年图形处理器和中央处理器的峰值运算性能的比值是10:1， 2014年NVIDIA公司的Geforce780Ti所提供的峰值吞吐量为5500GFL0PS，而Intel公司 的Ivy化idge架构的数据为700G化OPS。如此大的运算性能差距已经让越来越多的科研人 员将它们程序中的计算密集部分转移到图形处理单元来执行。
[000引因为比如像化en化和Direct3D技术需要直接针对巧片编程，而直到2006年程 序员必须通过应用程序编程接口来访问其核屯、，所W图形处理器的使用非常困难。该些应 用程序编程接口限制了能够为巧片编写的应用程序的类型，所W只有小部分程序员掌握了 利用图形处理器运行特定程序的必要技术，进而限制了该项技术的广泛传播。一切随着 NVIDIA公司在2007年推出CUDA技术而改变。因为规避了传统的应用程序编程接口，通过 CUDA，现在程序员们能够用他们所熟悉的类C/C++语言来直接针对巧片编程，从而使得并 行编程技术即便对于初学者来说也非常容易上手，大大地推进了该项技术在科学领域的应 用。
[0009] 因为格子玻尔兹曼方法是定义在正交网格上的显格式，并且只需要交换相邻节点 的信息，所W它非常适合于在图形处理器上并行化。发展至今，发表的关于格子玻尔兹曼方 法并行化的文章已非常丰富，并且相比于中央处理器，运行速度一般能提高两个数量级。虽 然格子玻尔兹曼方法的并行化已经日臻完善，但是将其与浸没边界法结合起来模拟流固禪 合问题的研究仍然没有起步，更没有全面透彻的分析浸没边界法潜在数据并行性的文章发 表。
【发明内容】
[0010] 针对现有技术存在的不足，本发明提出了一种基于GPU的高效流固禪合模拟方 法。
[0011] 为了解决上述技术问题，本发明采用如下的技术方案：
[0012] 本发明的具体操作步骤如下：
[0013] 1)在专业前处理软件中建立待测结构体及其内部流体的浸没边界体型，并划分边 界网格W供程序读取；
[0014] 2)程序读取边界体型数据，设置计算区域的边界条件和程序模型参数；
[0015] 3)根据边界当前时刻构型和参考构型利用有限元法计算边界施加在流体上的作 用力；
[0016] 4)将边界的作用力经由狄拉克函数扩散至周边流体；
[0017] 5)利用格子玻尔兹曼方法求解携带外力项的纳维-斯托克斯方程；
[0018] 6)同样经由狄拉克函数将流体速度插值到边界上得到边界的移动速度，进而更新 边界位置；
[0019] 7)重复步骤3-6,直至达到计算终点。
[0020] 所述专业前处理软件为Gambit或ICEM或ProE。
[0021] 本发明中，详细过程如下：
[0022] 一、采用浸没边界法处理流体和结构体之间的相互作用。
[0023] 具体来说包含W下几方面的特征：
[0024] 1)二维条件下，采用一维线网格离散结构体。S维条件下，采用非结构化S角形网 格离散结构体。
[0025] 2)二维条件下，同时考虑结构体抵抗拉伸和弯曲的能力，定义在节点上的抗拉和 抗弯力密度分别按下式计算：
【主权项】
1. 一种经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，其特征在于：步骤 如下： (1) 在专业前处理软件中建立待测结构体及其内部流体的浸没边界体型，并划分边界 网格以供程序读取； (2) 程序读取边界体型数据，设置计算区域的边界条件和程序模型参数； (3) 根据边界当前时刻构型和参考构型利用有限元法计算边界施加在流体上的作用 力； (4) 将边界的作用力经由狄拉克函数扩散至周边流体； (5) 利用格子玻尔兹曼方法求解携带外力项的纳维-斯托克斯方程； (6) 同样经由狄拉克函数将流体速度插值到边界上得到边界的移动速度，进而更新边 界位置； (7) 重复步骤3-6,直至达到计算终点。
2. 如权利要求1所述的经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，其 特征在于：所述专业前处理软件为Gambit或ICEM或ProE。
3. 如权利要求1所述的经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，其 特征在于：所述步骤3中： (1) 当浸没边界体型为二维时，采用一维线网格离散结构体，同时考虑结构体抵抗拉伸 和弯曲的能力，定义在节点上的抗拉和抗弯力密度分别按下式计算：
上式中的（Fs)1和（Fb)别表示定义在编号为1节点上的抗拉和抗弯作用力线密度； KJP Kb分别为抗拉系数和弯曲刚度；Λ a 1为编号1的节点的原始弧长；X 和Xm分别 为编号m+1、m-Ι和m的节点坐标；δ ml为克罗内克符号，当m = 1时其取值为1，其余情况 下取值为〇 ;nf为用来离散结构体的总节点数； (2) 当浸没边界体型为三维时，采用非结构化三角形网格离散结构体， 结构体在编号为m节点上的作用力密度计算公式如下：
其中：W为结构的总能量，3为偏微分符号，Xm为节点m的坐标；结构总能量的计算式 为： W = / wsdA 式中的dA为面元面积，而ws为结构的拉能密度，其计算式为：
，其中E为弹性模量，IjP 12分别为第一和第二应变不变量，可以 证明应变不变量能够按下式计算：
其中λ Jp λ 2是法平面内沿一对垂直方向的拉伸比，而a、b、c的计算式为：
上式中几何参数的定义如下：在结构体参考构型和变形后构型上选取两块相应的三角 形面元，分别在上述三角形中任意选取两条边，在参考构型上其边长分别为^和Γ C1，而在 变形后构型上其边长分别为1和Γ，灼和炉分别为上述两条边在参考构型和变形后构型 上的夹角。
4. 如权利要求1所述的经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，其 特征在于：所述步骤4和步骤6中的作用力扩散方程为：
上式中的f(x，t)为t时刻点X处流体受到的作用力密度，FmS编号为m节点上的作用 力密度，分别针对二维和三维条件，△ Sk分别是定义在编号为k的结构体节点上的弧长和面 积，而S (x-Xm)即为狄拉克函数，一般采取如下具体形式：
为了能够得到结构体形状的变化过程，我们需要根据结构体网格点的速度更新结构体 的构型，在已知流场速度分布条件下，结构体速度同样将利用狄拉克函数通过从流场插值 得到：
上式中的Uni分别为结构体网格点的坐标和速度，u(x，t)是t时刻点X处流体的 速度，Ah是网格步长。
5. 如权利要求1所述的经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，其 特征在于：所述步骤5中，格子玻尔兹曼方程为：
式中，ga (x，t)表示t时刻X处沿离散速度方向ξ α的分布函数，而g α (χ+ξ α St，t+ St)则为t+St时刻χ+ξ aSt处沿离散速度方向ξ α的分布函数，g】q(x，0是t时刻x处 沿离散速度方向ξ a的平衡态分布函数，τ是无量纲的松弛时间，与分子间的粘度系数有 关，St为时间步长。
【专利摘要】本发明提供一种经GPU加速的浸没边界-格子玻尔兹曼流固耦合模拟方法，步骤如下：在专业前处理软件中建立待测结构体及其内部流体的浸没边界体型，并划分边界网格以供程序读取；程序读取边界体型数据，设置所需测量点的边界条件和程序模型参数；根据边界当前时刻构型和参考构型利用有限元法计算边界施加在流体上的作用力；将边界的作用力经由狄拉克函数扩散至周边流体；利用格子玻尔兹曼方法求解携带外力项的纳维-斯托克斯方程；同样经由狄拉克函数将流体速度插值到边界上得到边界的移动速度，进而更新边界位置；重复步骤3-6，直至达到计算终点。本发明能够处理任意形状的结构体，计算过程中不需要重构底层流体网格，计算效率高。
【IPC分类】G06F17-50
【公开号】CN104866695
【申请号】CN201510355061
【发明人】吴家阳, 程永光, 张春泽, 刁伟
【申请人】武汉大学
【公开日】2015年8月26日
【申请日】2015年6月24日