对壁面流动进行了直接数值模拟（DNS）和大涡模拟（LES） 使用格子Boltzmann方法（LBM）和多个GPU（图形处理单元）。 在DNS中，采用了8个K20M GPU。 最大网格数为 ，这将导致无量纲网格大小为 对于整个解决方案域。 GPU-LBM解算器需要24小时才能进行模拟 LBM步骤。 测试了分辨率域的纵横比，以获得DNS的准确结果。 因此，平均速度和湍流变量，如雷诺应力和速度波动，都与Kim等人（1987）的结果完全一致，当时顺流方向和展向方向的纵横比分别为8和2。 对于大涡模拟，测试并使用了局部网格细化技术。 使用 网格和Smagorinsky常数 ，取得了良好的效果。 验证了LBM模拟湍流的能力和有效性。
1.简介 在过去的二十年中，格子Boltzmann方法被认为是模拟复杂物理流体流动的一种很有前景的替代方法。 与基于宏观连续体方程离散化的传统数值格式不同，LBM基于微观模型和介观动力学方程。 它有许多优点，包括易于实现边界条件、易于编程和完全并行算法[ 1 ]. 因此，LBM在生物流体、多孔介质等领域得到了广泛的应用[ 2 ]. 但其在湍流直接数值模拟和大涡模拟中的可行性和有效性仍有待进一步验证[ 三 , 4 ]。
充分发展的紊流壁面流动是一个经典问题，看似简单却又复杂。 它上面的DNS和LES，包括边界层流动，已经成为研究湍流机制的重要工具。 长期以来，简单的几何形状吸引了许多关于壁附近复杂湍流相互作用的实验和理论研究。 其中，Kim等人进行了直接数值模拟[ 5 , 6 ]被引用最多的是进行比较。 他们的DNS结果 ( 和 分别是摩擦速度和通道半宽） 已被视为数值模拟的基准。 由于DNS和LES的巨大计算成本，尽管近十年来高性能计算（HPC）的性能越来越高，但已发表的关于边界流的相关结果却很少。
目前，GPU的出现为大规模DNS和LES湍流预测带来了可能。 安装在计算机图形板上的GPU最初用于快速处理大型图形数据集。 最近，随着计算平台的出现和发展，例如CUDA[ 7 ]和OpenCL，使用GPU加速非图形计算已引起越来越多的关注[ 8 , 9 ]. 由于其浮点算术运算的高性能、宽内存带宽和更好的可编程性，GPU被进一步用于通用目的，即GPGPU。 基于网格的计算流体动力学（CFD）是一个重要的应用。 根据我们的经验，在单个GPU上使用标记和单元（MAC）求解器求解不可压缩Navier-Stokes方程来模拟流体流动是 比高度优化的基于CPU的实现快倍。 基于LBE GPU的计算甚至可以获得100多个加速比[ 10 – 13 ]. 青木集团[ 14 ]自2006年以来，一直致力于GPU高性能计算。 团队成员使用数百甚至数千个GPU进行了各种数值模拟，如天气预报、晶体生长过程和城市气流，PFLOPS取得了良好的性能[ 15 – 18 ]. 目前，GPU在能源、环境、航空航天等领域的工程应用成为热点。
在目前的工作中，DNS和LES在 使用晶格玻尔兹曼方法的D3Q19模型和多个GPU进行。 有两个目标。 一是提出了一种基于LBM和多GPU的高性能计算方法，为湍流预测提供了可能的途径。 二是验证LBM模拟湍流的有效性和能力。
本文的结构如下。 在节中 2 描述了带强迫项的DNS和LES的格子Boltzmann方程（LBE）。 章节 三 介绍了CUDA-MPI中仿真和数据传输的模型系统。 DNS和LES对充分发展的湍流壁面流动的结果在章节中进行了讨论 4 和 5 分别为。 总结和结论见第节 6 .
2.1. 用于DNS的LBE 在LBM中，首先，在离散晶格上求解以下离散速度分布的Boltzmann方程[ 19 ]: 哪里 是粒子速度分布函数， 是粒子在 th方向，以及 是速度的数量。 对于LBM，按格对不同方法进行分类的方法是DnQm方案。 这里“ “代表” 尺寸，“while” “代表” 速度。” 对于2D模型，9速LBM是最常见的。 对于3D模型，有几种立方晶格模型，如D3Q13、D3q15、D3Q19和D3Q27( 13、15、19或27）。 是碰撞操作符。 使用Boltzmann-BGK近似[ 20 ], 我们得到了 哪里 是局部均衡分布 是放松时间。 为了得到Navier-Stokes方程，必须仔细选择平衡分布函数形式。 在9速方形晶格和13、15、19和27速立方晶格中，提出了一个合适的平衡分布函数，如下所示：[ 2 ] 在本研究中，采用D3Q19模型 (0, 0, 0), (1, 0, 0), (−1, 0, 0), (0, 1, 0), (0, −1, 0), (0, 0, 1), (0, 0, −1), (1, 1, 0), (−1, 1, 0), (1, −1, 0), (−1, −1, 0), (1, 0, 1), (−1, 0, 1), (1, 0, −1), (−1, 0, −1), (0, 1, 1), (0, −1, 1), （0，1，−1），和 (0, −1, −1). 图 1 显示了我们在这里使用的D3Q19模型。 相应的权重因子为 , 、和 . 和 u个 是宏观量，可以评估为[ 2 ] 声音的速度是 ，状态方程是理想气体的状态方程，如下所示： 方程式( 三 )可以在物理空间中进一步离散化 和时间 .的完全离散形式( 三 )是 方程式( 7 )是著名的LBGK模型，其中 是无量纲弛豫时间。 宏观Navier-Stokes方程中的粘度可以由( 7 )作为 方程式( 7 )通常用标准形式求解，假设 按照以下两个步骤进行。 碰撞步骤： 流式处理步骤： 哪里 和 分别表示分布函数的冲突前和冲突后状态。
2.2. LES的LBE LES的LBE过滤形式定义为[ 21 – 23 ] 哪里 和 分别表示解析尺度的分布函数和平衡分布函数。 此外， 和 是对应于分子粘度的弛豫时间 和涡流粘度 分别为。 因此， 由提供[ 21 – 23 ] 涡流粘度 由标准Smagorinsky模型产生的可以用滤波器长度标度计算 和滤波应变率张量 作为[ 23 ] 在这里， 是Smagorinsky常数。 此外，本文采用有限差分近似计算应变率张量。
2.3. 外力 对于外部强制项，有许多实现方法。 在所提出的壁面流模拟中，改进的速度法[ 24 ]应用的是使用修改后的速度计算平衡分布函数。 具体描述如下： 哪里 是外力。 由于此处模拟的壁面流由恒定的物体力驱动，即均匀压降，我们定义 作为 在这里， 和 分别表示摩擦速度和通道半宽。
然后计算平衡分布函数 具有 而不是 .所以( 4 )可以作为
2.4. 边界条件 关于边界条件 ，在本工作的流向和翼展方向上使用了周期性边界条件。对于墙壁，应用了以下边界条件： 哪里 由宏观值计算 和 根据( 4 ). 不能直接计算，并且被假定等于其相邻内部节点的值。 和 由计算得出( 9 )和( 4 )分别使用相应的 和 .
3.模型系统和参数 图 2 显示了模型系统。 计算域长度为 , 、和 与流向相对应 ，正常 、和沿span方向 方向。 是墙法线方向的半宽，用于定义雷诺数。 标准雷诺数可由下式给出 平均中心线速度。 此外，通常使用“壁单元”或粘性长度来定义此类流动模型中的雷诺数 表面摩擦力计算如下 相应的“墙单位”时间为 式中，摩擦雷诺数为 这意味着 通常，上标（+）表示以“墙单位”计量的数量 和 .在雷诺数瞬态范围以外充分发展的湍流壁面流中， .
为了对湍流进行数值模拟，必须首先定义描述问题的无量纲物理参数。 对于壁面流动，计算域的空间范围对获得准确的结果起着重要作用。 特别是，纵横比为 和 .参考中的研究[ 26 ]，Spasov等人[ 27 ]进行了直接的数值模拟 在 和 。与Kim等人获得的结果相比，平均速度预测过高[ 5 ]. 此外，雷诺应力也不太一致。 在节中 4 ，我们验证了不同组的纵横比以执行DNS。
对于DNS，网格大小不能大于本地Kolmogorov长度 .在充分发展的湍流中， 在靠近墙壁的区域[ 28 ]。 随着离墙距离的增加而增加。 由于LBM中使用了统一的网格，因此整个领域的网格规模应满足 对于LBM， , ，这导致 因此，根据本研究中的摩擦雷诺数 ，半宽 在计算域中应该大于120。
至于初始条件， 和 由对数摩擦定律得出： ，其中von Karman常数 , 、和 设置为0.1。 速度场初始化如下： , 、和 ，其中波动项来自 , , 、和 基于墙附近通用平均速度剖面的常用工程近似 无量纲速度 以“墙单位”计量。选择的特定值， 对于速度剖面拐点， Van Karman常数，以及 ，预计不会对用于计算已开发湍流统计的最终流量剖面产生影响。
当前的并行计算是使用8个GPU通过1D域划分进行的。 图 三 显示了中的1D分区 -本研究中执行的指导。 数据与高层和底层进行交换。 有关数据传输的详细信息，请参阅我们以前的工作[ 12 ]. GPU之间的数据通信通过消息传递接口（MPI）和cudaMemcpy完成。 首先，数据通过CUDA API cudaMemcpy（…，cudaMemcpyDeviceToHost）从GPU复制到CPU。 然后使用MPI等并行工具在相应的CPU之间交换数据。 最后，通过cudaMemcpy（…，cudaMemcpyHostToDevice）将交换的数据从CPU复制到GPU。
4.边界流DNS DNS是在湍流壁面流上执行的 以及各种纵横比。 相应的无量纲网格大小为 对于整个解决方案域。 图 4 显示了（a）恒定流向位置的速度模式和（b）速度梯度张量第二不变量的等值面 当纵横比为 , .
4.1. 纵横比的影响 图 5 显示法向平均速度的剖面 -纵横比方向（a） , ，（b） , ，（c） , 和（d） , 其中以实线和虚线显示的DNS结果表示Kim等人获得的DNS结果[ 5 ]. 很明显 , 与Kim等人的结果完全吻合。虽然沿展向使用了周期性边界条件，但对于沿流向的相同尺寸，较大的展向比使结果更符合图 5（a） 带数字 5（b） 以及比较图 5（c） 带数字 5（d） 然而，顺流方向的尺寸，即， ，对平均速度几乎没有影响（对比图 5（a） 带数字 5（c） 和图 5（b） 带数字 5（d） ).
图 6 介绍了LBM-DNS计算的湍流统计数据，并与Kim等人的DNS数据进行了比较[ 5 ]对于（i）壁面法线方向上的雷诺应力，（ii）流向、展向和壁面法向速度波动中的均方根（rms）分量，at（a） , ，（b） , 分别为。 可以看出，我们的LBM-DNS结果和Kim的数据[ 5 ]对于纵横比的两种情况，它们基本上是一致的。 在 , 尤其是由 _rms（均方根）。 与平均速度和雷诺应力相比，速度波动的DNS结果对模拟所选的流向长度更为敏感。
4.2. GPU的计算性能 对于这项DNS工作，模拟大约需要24小时 LBM步骤 带有8个K20M GPU的LBM网格。 性能达到2333 MLUPS； 也就是说， 网格每秒处理一次。 与中的类似DNS相比[ 27 ]由36个CPU完成，其中性能为5.6MLUPS，GPU的当前性能大约高出416倍。 使用LBM-GPU，可以轻松实现高性能，而无需任何编程困难。
5.壁面流LES 由于DNS方法解决了所有相关的空间和时间尺度，它可以高保真地预测所有可能的流体运动，但成本高昂。 在这一部分中，我们通过将标准的Smagorinsky模型应用于大涡模拟来模拟壁面流。 此外，测试并使用了LBM的局部网格细化技术。
5.1. 局部网格优化 局部网格细化通常应用于预计解决方案会发生较大变化的区域。 对于LBM，局部细化的面片叠加到全局粗网格上，在一定程度上节省了计算时间，并在需要时实现了高分辨率。 网格细化是通过将空间步长除以细化因子来执行的 在LBGK模型的框架中定义的运动粘度取决于具有 .在粗格栅上保持恒定粘度和相同雷诺数( )和精细网格( )，放松时间 英寸( 8 )必须重新定义为[ 29 ] 在这里 和 分别表示细网格和粗网格上的松弛时间。 考虑到当 非常接近2，LBGK方案变得不稳定 ，我们估计( 20 )的上限 大约50岁。 在本研究中，细化因子的值 为2。 因此，精细网格上的时间步长减少了 ，其中 是粗网格上的时间步长。由于水动力变量及其导数在两个网格界面上的连续性，从( 4 ), ( 9 )、和( 20 )可以分别描述如下： 带细化因子的数值实现策略 如图所示 7 首先，对粗网格进行仿真 然后，根据粗网格的结果，在空间和时间上使用二阶插值将解映射到精细网格 ，根据( 22 ). 最后，根据得到的精细网格边界条件，及时提出粗网格和细网格上的因变量。
通过模拟在 使用或不使用局部网格优化。 使用局部网格细化，在中间部分使用粗网格（64×32节点），在1/4顶部和底部计算区域使用细网格（16×127节点）。 Ghia等人的NSE结果之间的比较[ 25 ]我们在这个问题上使用和不使用网格细化的LBE结果由图中法线方向中线上的速度剖面表示 8 显然，采用网格细化的LBE结果与Ghia等人的结果吻合较好[ 25 ]在边界附近，网格比均匀粗网格系统精细。
5.2. 结果和讨论 对于LES，首先，基于DNS，每个方向的网格数减少了4倍，纵横比保持为 , 也就是说，LES的总网格是DNS的1/64。 图 9 显示法线中平均速度的轮廓 -在各种情况下使用LBM-LES进行指导 对于网格系统 ，其网格比例为 .不同条件下的结果 是相同的，这意味着所使用的子网格模型在如此大的网格范围内不起作用。 同时考虑网格规模和计算成本 顶墙和底墙的厚度增加了一倍。 图 10 显示了在不同的 在局部网格优化系统中计算，其中 与墙壁相邻的顶部和底部零件被改进为 ( )中间部分保持在 ( ). 正如所观察到的，当 为0.13，小于 用于NS-LES[ 28 ]。
图 11 显示了法向雷诺应力的剖面 -使用LBM-LES的方向 对于上述相同的局部网格细化系统。 可以观察到，数值偏离了壁附近的LBM-DNS结果，这意味着网格不够精细，无法捕捉边界层中非常精细的涡流。 此外，常数的SGS模型 靠近墙壁时效率较低。
6.结论 在目前的工作中，DNS和LES是在 使用格子Boltzmann方法和多个GPU。 在DNS中，采用了8个K20M GPU。 验证了纵横比的影响，网格数为 以获得最大的纵横比。 结果表明，虽然沿展向和流向使用周期性边界，但两个方向的尺寸对DNS结果都有影响。 对于以下情况 , 平均速度和湍流统计变量，如雷诺应力和速度波动，与Kim等人的结果完全一致[ 5 ]. 同时，获得了2333 MLUPS的高性能。 对于大涡模拟，测试并使用了局部网格细化技术。 使用 网格和 ，在平均值上获得了良好的结果。 然而，湍流统计值与墙壁附近的LBM-DNS数据存在偏差。 这可能意味着应在墙壁附近使用更多网格，并应利用动态Smagorinsky常数来获得更准确的结果。 此外，多GPU具有超强的计算能力，与LBM的良好并行性完美匹配，在这项工作中表现出令人惊讶的高性能。
利益冲突 作者声明，本论文的出版不存在利益冲突。
致谢 作者感谢国家自然科学基金（编号：11302165和11242010）对这项工作的资助。 此外，这项研究还得到了日本科学促进会（JSPS）KAKENHI、科学研究补助金（B）编号23360046和青年科学家补助金（B）编号25870226以及日本科学技术署（JST）进化科学与技术核心研究（CREST）的部分支持 “后大规模计算的高效高性能应用程序框架”研究项目
