高雷诺数湍流风场大涡模拟的并行直接求解方法1)

包 芸 习令楚

(中山大学航空航天学院，广州 510275)

引言

土地沙漠化是当今人类面临的一个重要生态环境问题，风沙流是风沙环境力学的核心研究内容.风沙流发生在高雷诺数大气边界层中，同时具有湍流和多相流两个流体力学最具挑战的复杂问题[1].在这类复杂流体力学问题中，风场是沙尘运动和输运的动力源和基础.兰州大学在过去的二十年里，建立了世界上顶尖的野外风沙特性观察站，获得了大量实际野外风沙特性的第一手资料，并发性在大气边界层中存在超大型尺度湍流结构[2].同时他们还展开了风沙特性的风洞实验和数值模拟工作.在风沙特性包括风沙运输过程、电磁作用、沙尘暴等多方位的研究都取得大量的成果[3-9].风沙流等研究成果很丰富[10-12].

采用高雷诺数的壁湍流模型研究风场及其对沙尘特性的影响是当今研究风沙流特性的主要手段之一.计算机数值模拟是研究壁湍流及槽道流流动特性的重要方法[13-16].采用大涡模拟方法[17]，计算大气边界层的流动特性，已有很多研究工作和成果，并应用到城市大气污染扩散和其他相关大气风场流动等研究[18-25].

对于自然环境大气边界层，数值模拟研究用于风沙流的高雷诺数湍流风场，需要足够大的计算规模，才能有效的反映出实际野外环境中沙尘颗粒在气流作用下的运动规律.目前对高雷诺数湍流流动的数值模拟研究仍受到计算规模的限制.因此，充分利用我国的超级计算机资源，发展高效率并行计算技术，对高雷诺数湍流风场以及进一步的风沙流数值模拟研究有重要的意义.

在对高瑞利数湍流热对流进行DNS 数值模拟的研究中，创建了高效的并行直接求解方法[26]，完成了大量的高和极高瑞利数的湍流热对流DNS模拟，得到了丰富的计算数据，并在此基础上取得了系列湍流热对流物理特性的研究成果[27-30].在高雷诺数的壁湍流风场的数值模拟研究中，大涡模拟(LES)是有效的计算工具.本文的工作是将用于湍流热对流DNS 模拟的高效并行直接求解方法，扩展到三维壁湍流风场的大涡模拟计算中，建立新的用于三维壁湍流风场大涡模拟的高效并行直接求解方法，为大涡模拟计算研究高雷诺数的壁湍流风场提供有力的工具.

1 大涡模拟控制方程和亚网格模式

大气边界层中的空气流动，在Bussinesqe 假定下可认为是分层的不可压缩流体.风场大涡模拟控制方程是不可压NS 方程.无量纲化的不可压NS 方程为

其中，u为速度，P为压力，Re为雷诺数，τij是亚网格应力.

LES 方程和NS 方程的不同之处是LES 方程包含了亚网格通量.这些亚网格通量需要建立封闭模式，是大涡模拟的关键问题之一.关于复杂湍流中的亚网格通量模式仍在继续研究和发展中，目前最常用的模式是涡黏和涡扩散模型.Smagorinsky 涡黏模式是最早提出的亚网格模式，并被工业界广泛应用.涡黏模式的公式为

其中Smagorinsky 的涡黏模式

Cs是模式系数，∆是网格尺度，,而是可解尺度速度变形率张量.各项同性湍流理论可确定模式系数Cs=0.10.但在实际应用中发现，Smagorinsky 涡黏模式在底边界壁面附近耗散过大，需要调整模式系数.本文将采用Smagorinsky 涡黏模式作为亚网格模式进行大涡模拟计算.

大涡模拟的控制方程特性与NS方程DNS模拟在求解过程中没有本质的区别，可以用相同的数值求解方法进行数值计算.

高雷诺数湍流风场的数值模拟研究受到计算规模的限制，一直是湍流数值模拟研究工作中的难题之一.因此，在超级计算机上研究和应用高效的并行计算技术，对开展大规模高雷诺数湍流风场的模拟以及为后续沙尘颗粒等运动特性的联立求解，有着现实的意义.

2 不可压流动LES 模拟的并行直接求解方法

在过去的八年里，我们开展了极高瑞利数湍流热对流的DNS 模拟研究.由于极高瑞利数湍流热对流问题的特殊性，要求DNS 模拟的计算规模非常之大，使改进和发展新的并行模拟求解方法的工作成为必然.2016 年终于在规模并行计算技术上有了突破性进展，创建了不可压流动湍流热对流的并行直接求解方法(parallel direct method of DNS,PDMDNS)，并在“天河二号”超级计算机上实现了系列高瑞利数湍流热对流的DNS模拟计算，并得到了系列的高瑞利数湍流热对流特性的研究成果.

本文拟将湍流热对流DNS 模拟的并行直接求解方法，扩展到湍流风场的LES 模拟中，希望能显著有效的提高壁湍流风场LES 模拟的计算效率，由此可进一步扩大壁湍流风场的计算规模，提高计算壁湍流风场的雷诺数.

2.1 不可压流动LES 的并行直接求解方法

由于该计算方法只限用于水平采用等距网格的矩形计算区域，适用于风场的LES 模拟计算，如图1 所示的风场示意图.

图1 风沙流中的风场计算域示意图Fig.1 The wind field calculation domain schematic

采用LES 模拟不可压NS方程，本文采用投影法分步骤计算.动量方程中预估速度的计算采用容易实现并行的显示格式.联立压力驱动速度项和连续方程，得到压力泊松方程，需要全流场联立求解，是数值求解过程中计算工作量最大的部分，同时全场联立求解也给大规模并行造成困难.大规模并行计算的压力泊松方程高效求解技术，是解决高雷诺数湍流风场大规模并行计算的关键.

在二维湍流热对流DNS 模拟的并行直接求解方法PDM-DNS[26]中，压力泊松方程是二维的，因此仅需要在水平方向进行FFT 变换解耦泊松方程，将方程变为三对角方程，再利用PPD 并行算法，综合建立高效并行计算方法.

在“天河二号”超级计算机上，大规模的并行计算可以进行MPI 和OpenMP 混编，MPI 并行计算需要交换分区边界的数据.

类似二维湍流热对流DNS 的直接求解方法，在三维湍流风场的LES 模拟中将矩形的风场计算区域沿水平面对z进行分割.相邻的MPI 分割区域在并行计算中需要数据通讯，区域内部可用OpenMP 并行，无需数据通讯.在此基础上进行压力泊松方程的水平面FFT 变化时将不用并行计算，而是在纵向三对角方程求解时使用PDD 技术进行并行计算.因而大大的提高了泊松方程并行求解的计算效率.

由此可见，二维流动的并行直接求解方法扩展到三维计算，泊松方程的求解过程仅是在FFT 解耦的过程上有所不同.二维流动计算中用一维FFT，三维流动计算中用平面二维FFT 变换，其它求解过程基本一致.

对三维湍流风场的LES 模拟，同样并行计算的困难在于求解压力泊松方程.三维泊松方程的直接求解方法的具体做法为，在三维空间的xoy水平平面上使用二维平面FFT 变换，将空间3 个方向上都要求联立求解的压力泊松方程解耦，使泊松方程变换成为只在z方向上的三对角方程.求解三对角方程后通过FFT 反变换就可得到全场泊松方程解.

三维压力泊松方程为

其对应的上下压力无梯度边界条件为

沿流向及展向方向的压力边界条件为周期边界条件.

在x和y方向使用等距网格，z方向可使用变距网格.在点(i,j,k)对应的3 个方向网格长度分别为∆x,∆y和∆zk，压力泊松方程的二阶精度中心差分离散格式可写成如下形式

根据以上离散泊松方程的特点，使用xy平面上的二维离散余弦傅里叶变换，可以对泊松方程的离散方程进行解耦.为自动满足压力的边界条件，选用二维离散余弦傅里叶变换，表达式为

M对应x方向的网格数Nx，N对应y方向的网格数Ny.以上变换可以通过使用FFTW 软件包实现.将上式变换应用到压力泊松方程中，可得到沿z方向联立的三对角方程

对于给定的p和q，变换后的压力泊松方程在需要xy和z方向上解耦，变为系列的三对角方程组.完成所有三对角方程组的求解后，通过对应的变换即可完成压力泊松方程的求解.

利用以上高效的压力泊松方程并行直接求解方法，联合其它易并行的动量方程等计算，本文将已创建的二维湍流热对流并行直接求解方法，扩展到三维不可压湍流风场的LES 模拟中，建立LES 模拟的并行直接求解方法，为大规模高效并行计算高雷诺数湍流流动，提供全新的数值计算技术和计算方法.

2.2 并行计算规模和计算效率

实际计算研究中，壁湍流风场的大涡模拟已有许多成熟的研究成果.但在实际自然灾害问题的风场研究中，空间尺度都很大，从而导致风场具有很高的雷诺数，需要计算规模巨大.如何实现大规模高雷诺数风场数值计算仍是一个瓶颈问题.本论文的研究特点主要在采用新的并行计算技术，建立了一个可进行大规模高效并行的计算方法，可以显著的提高风场大涡模拟的计算效率，为实现高Re风场以及由风带来的物质输运的计算研究提供有价值的计算工具.

采用半槽道壁湍流大涡模拟计算，对新方法PDM-LES 的并行效率进行验证.

计算域为三维矩形.作为验证计算，计算网格在(x,y,z)方向均取256×256×256.以x作为流动方向，y方向为与x方向垂直的另一水平方向，z方向为垂直地面方向.边界条件为x方向上下游及y方向为周期边界，z方向下底面为无滑移边界，上边界为对称边界条件.在计算过程中，采用定流量计算方法，以保证流动持续进行，即在每个时间迭代计算中对由于粘性摩擦损失的流量进行修正.

计算在“天河二号”超级计算机上进行.每个计算机节点可以有32 线程无需交换数据的并行计算，采用OpenMP 方式进行并行计算.节点间在并行计算是需要交换数据，采用MPI方式进行并行计算.

表1 杆不同计算节点的计算效率Table 1 The computational efficiency of different compute nodes

表1 给出了采用不同节点数的计算结果，分别采用1，2 和4 个节点进行了计算效率的对比计算.以1 个节点的计算作为基础计算工作量，可以看到，采用1 个节点每天可迭代71.6 万个时间步长，当用4 个节点共128 核的并行计算，加速比为3.59，计算效率达到90%.

由此可见，本文建立的针对风场的壁湍流大涡模拟的并行直接求解方法(PDM-LES)，具有较高的并行计算效率，可用于规模风场的数值计算研究.

3 三维风场的规模并行计算

采用建立的PDM-LES 方法，对三维壁湍流风场进行大涡模拟计算.取来流Re=10000 进行试算，检验壁湍流计算结果的合理性.在以高度为1 的无量纲计算域中，计算域选取流动流向、展向和垂向长度分别为8×2×1，网格分别为Nx⊆Ny⊆Nz=1024 ⊆192 ⊆160.水平采用等距网格，垂向采用变距网格，计算网格规模约为3.800×107.

图2 给出了来流Re=10 000 时三维瞬时速度分布.三维瞬时速度分布图中可以看到，湍流流动的速度分布存在明显的脉动特性.计算结果得出，在现有计算条件下Reτ≈4000.

图2 三维瞬时速度云图Fig.2 3D instantaneous speed cloud map

图3 为不同高度xy平面的瞬时速度分布.图3(a)为近壁面平行于xy平面上的瞬时速度分布(z+=12)，可以看到壁湍流近壁面区域分布着众多的沿流动方向分布的低速条纹和高速条纹.图3(b)为离壁面距离略高的xy平面上的瞬时速度分布(z+=200)，仍然可以看到明显的条带状速度分布.随着高度的增加，条带状速度分布的尺度增大.这些条纹反映出壁面附近的流动存在沿流向方向的条带状拟序结构.瞬时速度条带状拟序结构是湍流流动的最主要流动特征，并且起到控制湍流流动特性的重要作用.

图3 不同高度xy 平面的瞬时速度分布Fig.3 Instant velocity distribution of xy planes at different altitudes

计算湍流平均速度特性.选用多个不同时刻的瞬时速度场取平均值，再对水平平面的速度进行平均计算，得到壁湍流流动在时间和空间上的平均速度场湍流特性纵向分布.

图4 中给出了平均速度场U+和雷诺应力的纵向分布特性.图4(a)为速度U+的纵向分布，纵向距离z+取对数.可以看到，第一个数值点的距离小于1 个z+，表明近壁面网格距离基本满足计算的要求，本算例的计算网格选取是合理的.U+的纵向分布在z+≈20 处分为两段，近壁面区域U+的分布曲线是弯曲的，离开壁面区域U+随z+的分布在单对数坐标中呈现直线关系，表明速度分布存在对数律.图4(b) 是雷诺应力的纵向分布情况.雷诺应力的分布特性符合湍流边界层的特性规律，但雷诺应力的最大值偏小.

图4 壁湍流平均速度场特性Fig.4 Average velocity field characteristics of wall turbulence

图4 壁湍流平均速度场特性(续)Fig.4 Average velocity field characteristics of wall turbulence(continued)

总体上说，本文计算的壁湍流风场的湍流特性基本合理.

下一步工作是，在“天河二号”超级计算机上开展高雷诺数的壁湍流大涡模拟计算，增加计算网格规模，同时也增加有效的并行计算线程数，使得能在可接受的时间内完成高雷诺数的壁湍流大涡模拟计算，为风沙流的模拟计算等研究服务.

4 结论

利用大涡模拟数值研究三维高雷诺数风场特性，超大规模的数值模拟计算是十分必要的.依靠超级计算机技术的发展，建立高效并行的计算方法，对更深入研究计算实际野外风沙流的物理和流动特性具有很重要的意义.

(1)本文将高效的并行直接求解方法扩展到三维壁湍流风场的大涡模拟计算中，建立了三维壁湍流风场大涡模拟的高效并行直接求解方法PDM-LES.

(2)对并行直接求解方法的并行计算效率进行的验证计算结果表明，在“天河二号”超级计算机上进行的1 至4 个节点，最大用到128 核的并行计算效率达到90%.

(3)计算Re=10 000 的结果显示，湍流风场瞬时速度分布在近壁区存在条带状的拟序结构，并远离壁面条带结构尺度增加.平均场的速度分布符合湍流速度分布特性，即近壁面在z+≈20 之内U+为线性分布，在外区U+满足对数律分布.计算得到的风场湍流特性基本合理.