LES和DNS在平坡异重流数值模拟中的比较研究

Comparative Study of LES and DNS on Numerical Simulation of Gravity Current over a Flat Bed

摘要: 本文利用FLUENT软件模拟了开闸式异重流的演变特性。使用二维大涡模拟(LES)和直接数值模拟(DNS)，对异重流的形态变化、瞬时掺混系数、势能转变，动能剖面变化等进行了相关的对比分析。结果表明，LES和DNS对于异重流头部位置的预测都具有较好的准确性，DNS的结果略微精确，但需要较多的计算网格数。此外，在异重流与环境流体的交界面发展过程中，相比于LES模拟的结果，DNS的结果具有较少的K-H不稳定性产生。DNS模拟的异重流基本上集中于头部，而且DNS模拟的结果与实验异重流的形态较为接近。LES和DNS模拟计算得出的异重流与环境流体的掺混程度相近，但DNS的交界面较为光滑。LES和DNS对于异重流沿途势能变化的结果相近，但DNS的结果具有较小的背景势能。对于特征断面处动能逐时变化，相较于LES，DNS模拟较少出现多次波峰的震荡现象。本文的结果可为后续模拟异重流运动时，提供选择数值模型的参考依据。

Abstract: In this paper, the FLUENT software was used to simulate the evolution characteristics of the lock-exchange gravity currents over a flat bed. Two-dimensional large eddy simulation (LES) and direct numerical simulation (DNS) were utilized to compare the morphological change, instantane-ous mixing coefficient, potential energy transformation, and kinetic energy at the characteristic section. The results show that both LES and DNS have good accuracy in predicting the head position of the gravity currents, but the results of DNS are slightly more accurate with more computational grids required. For simulating the development process of gravity currents, the DNS results show less K-H instability occurring at the interface of current and ambient fluid than the LES simulations do. For DNS simulations, heavy fluid basically concentrates on the current head. In addition, the DNS results are closer to the gravity current morphology from experimental observations. The mixing degree of density flow with ambient fluid estimated by LES and DNS simulations is similar, but the result of DNS shows smoother current-ambient fluid interface. LES and DNS results obtain similar potential energy changes along the course of density flow, but the results of DNS have smaller background potential energy. For the temporal variations of kinetic energy at characteristic cross-sections, compared with LES, DNS performs fewer oscillating peaks. The results of this paper can provide a reference for the selection of numerical models on gravity current simulations in the future.

1. 引言

异重流作为由密度差导致的重流体入侵轻流体的流动现象，在自然界和工程界中广泛存在，是常见的流体流动形式 [1]。水库泥沙冲淤、大坝开闸泄洪、沙尘暴、海底浊流等现象都与异重流有密切联系 [2] [3]。

近年来，随着计算机技术的快速发展和各种数值方法的改进，数值模拟逐渐成为了研究异重流的重要手段。根据所采用湍流模式的不同，数值模拟可分为雷诺平均模拟(Reynolds-Averaged Navier-Stokes Simulation, RANS)、大涡模拟(Large-Eddy Simulation, LES)和直接数值模拟(Direct Numerical Simulation, DNS) [4] [5] [6]。RANS可以计算高雷诺数的复杂流动，但给出的是平均运动的结果，不能反映流场紊动的细节信息 [7]；大涡模拟通过将原有控制方程进行过滤，只计算大于过滤尺度的湍流，而将小于过滤尺度的湍流信息用亚网格应力加以刻画，以得到模拟异重流运动的大涡模拟模型 [4]；而DNS通过综合考虑流场的所有细节，可以获得湍流场的精确信息，但是会耗费极高的计算机资源 [8]；

通过模拟在高格拉斯霍夫数(Grashof Number)时开闸式异重流运动形态的时变特征，Ooi等人发现二维LES模型能够模拟包括头部结构、K-H不稳定性等在内的大部分动力特性 [9]。通过对比异重流沿坡运动的二维和三维LES模拟结果，Nourazar和Safavi证明了二维模型可以捕捉异重流的某些动力过程 [10]。Härtel等人 [11] [12] 通过对开闸式无颗粒异重流在平坡上运动的进行三维DNS模拟，着重分析了异重流头部的耳垂和裂缝结构以及边界条件对异重流头部特性的影响。利用二维DNS模型，Birman等人 [13] 模拟了非布辛奈斯克(Non-Boussinesq)无颗粒异重流的运动过程，并研究了边界条件、雷诺数以及初始时刻流体之间的密度差异对异重流头部速度和能量变化的影响。为了研究边界条件和初始闸门相对高度对异重流运动特性的影响，Liu和Jiang [14] 通过开源软件Open FOAM构建了相应的DNS模型，并发现边界条件对异重流运动的影响会受到初始闸门相对高度的影响。

总体而言，在异重流数值研究中，前人的研究主要集中于LES和DNS各自的优缺点，较少有关注它们之间的具体差异。本文将使用FLUENT软件，比较LES和DNS在二维平坡水槽异重流数值模拟中的具体差别，并研究这两种模型下异重流运动时的瞬时掺混系数、能量变化等特性参数。本文结果可为后续异重流模拟中，湍流模型的选取提供参考依据。

2. 数值模型

本研究参照林颖典等人开闸式异重流实验 [15]，分别使用LES和DNS模拟异重流的运动过程，其物理模型如图1所示。平坡水槽闸门右侧为密度较大的盐水，闸室宽 $x_0$ = 0.1 m，密度为 ${\rho }_d$，水深为 $h_0$ = 0.15 m；左侧为环境水体，水体密度为 ${\rho }_a$，水深为H = 0.15 m，水槽全长为L = 2 m。

Figure 1. Schematic diagram of physical model

图1. 物理模型示意图

2.1. 数值模型

本文的数值研究对象为如图1所示的开闸式异重流，其计算区域分为初始异重流区域(图1中红色区域)及环境流体区域(图1中蓝色区域)。

2.1.1. 网格及求解方法

采用笛卡尔坐标系，并结合矩型结构网格对计算区域进行剖分。在LES模型中，亚格子模型(Subgrid-Scale Model)采用Smagorinsky-Lily假设，DNS模型中选用层刘默写(Laminar flow)，多相流模型都使用Mixture模型，压力采用PRESTO!算法，压力与速度的耦合采用SIMPLE(压力速度耦合方程组的半隐式算法)，动量方程采用有界中心差分格式(Bounded Central Differencing)，体积分数采用一阶迎风格式。当各计算变量残差小于0.001时可认为计算收敛。

2.1.2. 边界条件

顶部边界(图1绿色区域)选用剪力为0的边界假设 [15]，左、右和下部壁面(图1黑色区域)采用无滑移边界条件。异重流与环境流体分界面(图1黄色区域)采用interior (内部边界)条件。

2.2. 模型验证

首先，将数值模型的结果和林颖典等人实验研究结果 [15] 进行对比，以验证模型的准确性。这里采用其文章中的典型工况N2进行验证。对于LES模型，对比了网格数量为7.5万(2 mm × 2 mm)、30万(1 mm × 1 mm)、120万(0.5 mm × 0.5 mm)的数值结果，发现当网格数量大于30万结果已经收敛。因此，出于准确性和节省计算资源考虑，LES模型选用网格大小为1 mm × 1 mm，共计数量为30万。而对于DNS模型，由于其对网格的要求比较高，因此选用1 mm × 1 mm (30万)、0.8 mm × 0.8 mm (约45万)、0.7 mm × 0.7 mm (约60万)、0.58 mm × 0.58 mm (约90万)总计4种网格大小，其对比结果如图2所示。

Figure 2. Comparison of temporal head position of gravity current between numerical and experimental results

图2. LES和DNS中异重流头部位置随时间变化与实验结果 [5] 对比图

图中，无量纲头部位置X及无量纲时间T为：

$X=\frac{x_f-x_0}{H}$ (1)

$T=\frac{t\sqrt{g^{\prime }H}}{H}$ (2)

其中， $x_f$ 为异重流的头部位置，t为实际时间， $g^{\prime }$ 为约化的重力加速度。

从图中可以看出，对于DNS模型，当网格数量超过45万后，其结果之间已经非常接近，可认为收敛。且DNS和LES和实验结果都具有很好的一致性，但DNS比LES模型具有更高的精确度，尤其是在异重流发展的后期。因此，在下文中出于对准确性与节省计算资源的综合考量，LES模型选用1 mm × 1 mm网格总计数量30万，DNS选用0.8 m × 0.8 mm网格总计数量约45万。

2.3. 特征参数及组别

本研究工况分为LES工况和DNS工况，除所采用模型不同以外，其余参数保持一致，具体如表1所示。

Table 1. Parameters of numerical cases

表1. 数值工况参数

3. 结果分析

3.1. 流态分析

图3为LES模拟(LES-N3)和DNS模拟DNS-N3)的异重流发展过程，并做了相应的验证实验进行比较，分别都选取了6 s、9 s、14 s和27 s五个典型时刻。对比图3(a)和图3(b)可知：LES模型模拟的结果和DNS模拟的结果在异重流整体运动的趋势上接近，但LES模型的结果会产生更多的K-H波，尤其是发展的后期更明显；且从图中可以很容易的发现，相较于LES，DNS模拟的异重流所占有的体积较少，且绝大部分体积被头部占有，尾部只占很少的部分。对于异重流头部形态而言，DNS模拟的结果与实验异重流的形态较为接近。

Figure 3. Development processes of gravity current in the case N3: (a) LES; (b) DNS; (c) Experimental results

图3. N3工况异重流的发展过程：(a) LES；(b) DNS；(c) 实验结果

3.2. 卷吸和掺混

异重流发展时会由于K-H不稳定性的存在导致自身与环境流体不断的卷吸和掺混，并使得自身占有的体积不断增大。本文参考Ottolenghi等人对瞬时卷吸系数(Instantaneous Entrainment Coefficient)的定义 [5]，使用下式衡量异重流运动的任意i时刻，异重流与环境流体的卷吸掺混程度，表示如下：

$E_i=\frac{2\left(M_i-M_{i-1}\right)}{\left(x_{f,i}+x_{f,i-1}\right)\left(U_{b,i}+U_{b,i-1}\right)\left(t_i-t_{i-1}\right)}$ (3)

其中， $M_i$ 为i时刻下等密度线ρ^* ≥ 0.02时所包围的异重流的侧向面积 [5]， $x_{f,i}$ 为i时刻的异重流头部位置， $U_{b,i}$ 为0至i时刻的异重流的平均速度， $t_i$ 为计算开始至i时刻的时间。

图4为异重流的掺混系数随头部位置变化关系图。从图中可以看出，不论是LES模拟的结果还是DNS模拟的结果，瞬时卷吸系数都在不断地变小，这是因为异重流在发展过程中导致掺混的驱动浮力不

Figure 4. The variation of entrainment parameter with front position of gravity current: (a) N1; (b) N2; (c) N3

图4. 异重流的掺混系数随头部位置变化关系图：(a) N1；(b) N2；(c) N3

断减弱。在LES模拟的结果中，异重流发展的后期，瞬时卷吸系数有一个较为明显的波动，说明在这个阶段异重流与环境流体掺混不稳定。其原因为：一方面，开闸式异重流往前传播时，由于头部没有得到及时重流体的补充，导致密度不断被稀释；另一方面，由于K-H涡使一部分重流体被抬升，导致与环境流体的接触面积增大，从而使得卷吸作用增强。DNS的结果中，瞬时卷吸系数和LES模拟的结果在数值大小上具有较好的一致性，但DNS的瞬时卷吸系数的结果随时间变化较为平滑，这是因为DNS模拟的异重流产生的K-H涡较LES少，相比之下，密度不断被稀释对掺混影响的结果占优，导致瞬时卷吸系数较少有波动。

3.3. 势能转换

流体势能(Potential Energy)是异重流的初始动力源。异重流运动时，总势能被分成来了可用势能和背景势能。可用势能用于维持水平前进的动能，不可逆势能转化为异重流与环境流体交界面上不可逆的掺混耗散能，因此，背景势能是衡量异重流掺混强烈程度的有效参数 [5]。根据Winters等人对流体能量的定义，异重流体系的总势能、背景势能和可用势能可以分别表示为 [16]：

$E_p\left(t\right)=g{\displaystyle \underset{V}{\int }\langle \rho \left(x,z,t\right)\rangle z\text{d}V}$ (4)

$E_b\left(t\right)=g{\displaystyle \underset{V}{\int }\langle \tilde{\rho }\left(x,z,t\right)\rangle z\text{d}V}$ (5)

$E_a\left(t\right)=E_p\left(t\right)-E_b\left(t\right)$ (6)

其中 $\langle \rho \left(x,z,t\right)\rangle$ 为流场平均密度；V为整个流场流体总体积(包含异重流和环境流体)； $\langle \tilde{\rho }\left(x,z,t\right)\rangle$ 是流场在绝热条件(与外界无热交换)下重新排布为势能最小状态(流场完全稳定水平分层)时的密度场。基于此，可用势能 $E_a\left(t\right)$ 可以理解为 $\langle \rho \left(x,z,t\right)\rangle$ 绝热转变为 $\langle \tilde{\rho }\left(x,z,t\right)\rangle$ 释放的能量。

图5为N1工况下，异重流总势能 $E_p$，背景势能 $E_b$ 和可用势能 $E_a$ 随异重流沿程头部位置变化的关系图，并利用系统初始时刻总势能 $E_p\left(0\right)$ 对能量参数无量纲化。从图中可以看出，不论是LES的结果还是DNS的结果都具有如下的共通性：在异重流速度恒定的初始坍塌阶段 [17]，总势能 $E_p$ 的变化趋势与可用势能 $E_a$ 变化趋势高度吻合，而背景势能几乎为0，这说明此阶段总势能完全转变为动能，异重流与环境流体几乎不掺混。背景势能 $E_b$ 随后逐渐增大，也就意味着发生在异重流与环境流体交界面不可逆的掺混逐渐增强，表现为K-H不稳定性主导的涡旋。总势能曲线与可用势能曲线逐渐分离，直到背景势能与可用势能曲线慢慢接近。但相较于LES的结果，DNS的结果中背景势能 $E_b$ 所占的比重更小，这是因为DNS模拟的结果中，异重流的不可逆掺混程度较低，导致描述异重流与环境流体不可逆掺混的间接度量 $E_b$ 也更小一些。

Figure 5. The variation of the dimensionless total potential energy, background potential energy, and available potential energy with head position of gravity current in N1: (a) LES; (b) DNS

图5. N1工况下异重流无量纲总势能、背景势能和可用势能随头部位置变化关系图：(a) LES；(b) DNS

3.4. 动能峰值

动能(Kinetic Energy)是异重流水平运动强烈程度的直接度量。此处选取了三个特征断面，分别是距离闸室0.2 m处(P1)，0.8 m处(P2)和1.4 m处(P3)，给出各个特征断面处的动能随时间变化的图像(图6)，并利用下式对动能进行计算 [18]：

$K=\frac{1}{2}{\displaystyle \underset{\Omega }{\int }{\bar{u}}^2+{\bar{w}}^2\text{d}\Omega }$ (7)

其中， $\Omega$ 表示积分断面， ${\bar{u}}^2$ 和 ${\bar{w}}^2$ 分别表示某时刻单位质量流体在主流方向和垂直方向上瞬时速度平均值的平方。其中，平均值是指某一计算时刻前后各0.05 s (即0.1 s内)的平均。

图6比较了N3工况时，LES和DNS在各个特征断面处动能的变化情况。从图中可以看出，LES和DNS模拟的结果基本上保持一致。在各个特征断面处，DNS的峰值都比LES模拟的结果略大。但他们有如下区别 LES模拟的结果(图6(a))和DNS模拟的结果(图6(b))在P1和P2断面处容易出现多个波峰，而P3则不会。这是因为当异重流的头部运动至P1和P2断面时，使得该断面的流体动能增大，但由于和环境流体的掺混导致浮力损失，速度下降，即出现第一次波峰。随后，异重流身体后方的补充流体到达该区域时，再次激活该断面处流体的动能，于是出现第二个波峰。但随着异重流流经更远的P3断面时，由于和环境流体掺混较为完全，后方的流体无法及时补充头部流体，于是便不会出现多个波峰。但相对于LES模拟的特征断面处动能的结果(图6(a))，DNS模拟的结果(图6(b))在各个特征断面处都具有较为光滑的变化，这是因为DNS的结果中绝大部分的能量都集中在异重流头部区域，后方补充的流体很少因此较难形成多次波峰。

Figure 6. Kinetic energy profile at the characteristic section of density currents of the case N3: (a) LES; (b) DNS

图6. N3工况下异重流特征断面处动能峰值图：(a) LES；(b) DNS

4. 结语

本文通过选取LES和DNS湍流模型对开闸式异重流进行比较研究，得出如下结论：

1) LES和DNS对于异重流头部位置的预测都具有较好的准确性，DNS的结果略微精确，但所需网格数较多；

2) 在异重流与环境流体的交界面发展过程中，相比于LES模拟的结果，DNS的结果具有较少的K-H不稳定性产生，而DNS模拟的结果与实验异重流的形态较为接近；

3) LES和DNS模拟计算得出的异重流与环境流体的掺混程度相近，但DNS的交界面较为光滑；

4) LES和DNS对于异重流沿途势能变化的计算结果相近，但DNS的结果具有较小的背景势能；

5) 对于特征断面处动能逐时变化，相较于LES，DNS模拟较少出现多次波峰的震荡现象。

本文的结果提供了LES和DNS模型模拟平坡异重流的特性及各项参数的比较，也对实验异重流的头部位置和运动形态进行验证，确定不同模型的准确性，可为后续进行异重流模拟时，湍流模型的选取提供参考依据。

基金项目

国家重点研发计划(2017YFC0405205)；浙江省自然科学基金资助项目(LY20A020009)；中央高校基本科研业务费专项资金资助(2020QNA4038)。

NOTES

^*通讯作者。