盐指型双扩散湍流二维直接数值模拟研究

杨延涛

(北京大学 工学院，湍流与复杂系统国家重点实验室，北京 100871)

0 引 言

当流体密度由于温度等标量场的空间不均匀分布出现差异时，流体在重力场的作用下会发生浮力驱动的对流湍流。对流湍流是众多自然界和工程流动中热量和物质混合的关键物理机制。在地球海洋环境中，由于海水密度主要受温度和盐度的影响，温度和盐度的非均匀分布必然会导致浮力驱动的宏观流动发生。这种由温度和盐度共同驱动的对流现象称为双扩散对流湍流。特别地，由于温度的分子扩散速率比盐度快大概两个量级，双扩散湍流表现出单标量场对流湍流所不具有的众多独特动力学过程和丰富的流动机理。

温盐驱动双扩散湍流广泛分布于海洋上层海水中[1-2]。在热带和亚热带海域，由于日照加热和表层海水蒸发，上层海水中特定深度范围内温度和盐度随深度增加而下降，此时盐度梯度驱动流动而温度梯度致稳流动；双扩散湍流处于盐指流态[3]。在高纬度极地海域，由于表层海冰低温和融化，上层海水中存在温度和盐度随深度增加的区域，此时温度梯度驱动流动而盐度梯度致稳流动；双扩散湍流处于扩散流态[4]。无论哪种流态，双扩散湍流导致的温盐台阶形态都具有垂直尺度较小而水平相干性尺度极大的特征，并且对热量和盐度的垂直输运有重要的影响[3-5]。

对于盐指型双扩散湍流，对观测数据的分析表明其可能存在于超过30%的海洋水体中[1]，且能够大大增强当地海域的垂直混合效率[3]。因此理解盐指型双扩散湍流的流动结构和输运规律，并提出参数化描述，具有非常重要的科学意义和应用价值。自盐指型双扩散湍流机制被提出以来[6]，相关研究工作非常丰富，包括实验、计算和理论工作[7-9]。虽然海洋双扩散湍流是温度-盐度体系，但是在实验中，由于实验条件的限制，也有采用如盐度-糖度体系[10-11]、温度-铜离子体系等[12-13]。数值模拟方面，由于盐度场极小的分子扩散率对计算网格分辨率提出非常高的要求，一般为节省计算资源会采用具有较小施密特数的浓度场进行模拟，而不直接采用与海水相同的物性参数[14-15]。

基于自主开发的高效直接数值模拟程序[16]，本课题组针对盐指型双扩散湍流已开展了一系列数值模拟工作，包括真实海水物性参数模拟[17-18]，以及不同物性参数组合体系的模拟[19]。但是，之前的研究只针对固定密度比和少数不同瑞利数展开，关于不同物性参数组合的系统性研究还没有进行[19]。因此，本论文主要针对四种不同物性参数组合，通过变化瑞利数和密度比开展系统性的二维直接数值模拟工作，并在此基础上考察流动结构和传输特性的变化规律。

1 控制方程与数值方法

1.1 控制方程

本文考虑两平行平板间的流体流动。平板间距为H，且重力方向垂直于平板。上下板分别给定恒定的温度和盐度值。本研究采用线性状态方程ρ=ρ0(1 −βT θ+βSs)，其中：ρ为密度，ρ0为参考状态密度，θ为偏离参考状态温度，s为偏离参考状态盐度，βT为热膨胀系数，βS为盐度导致密度变化的线性系数。本文中脚标T代表与温度相关物理量，而脚标S代表与盐度相关物理量。则在Oberbeck-Boussinesq假设下不可压缩流动动量方程和温度、盐度输运方程为：

其中，u为速度矢量，p为压强，υ为运动学黏性系数，g为重力加速度，κT为温度场分子扩散率，κS为盐度场分子扩散率。

温度和盐度用两平板间的恒定温度差ΔT和盐度差ΔS无量纲化，速度用自由落体速度U= (g βSΔSH)1/2无量纲化，长度用板间距H无量纲化。则流动控制参数写为：

本文还同时使用刘易斯数，其为普朗特数和施密特数的组合：

注意Le代表两个组分分子扩散率的比值，而密度比代表两个组分差引起的密度差异的比值。对于盐指型流态，密度比越大，表明致稳流动的温度梯度相对于驱动流动的盐度梯度的强度越大。

1.2 数值方法

本文的流动模型中，上下两个平板处速度场满足无滑移和无穿透边界条件，温度场和盐度场恒定。水平方向所有物理量满足周期边界条件。计算中选取区域水平宽度远大于流动结构的特征水平尺度以保证周期边界条件的有效性。由于计算中涉及大施密特数和大瑞利数情况，本文采用二维模型，即只考虑一个水平方向。以往研究表明，对于盐指型双扩散湍流二维和三维流动非常类似[20]。

直接数值模拟采用自主开发程序[16]。程序采用分步-投影法处理速度-压强耦合。空间离散采用二阶中心有限差分格式。时间方向采用二阶Runge-Kutta格式，其中扩散使用类似Crank-Nicolson格式的半隐式处理方法，非线性对流项采用类似Adams-Bashforth格式的显式处理方法，浮力项处理方法与非线性项一致。针对大施密特数情况，程序特别使用了双分辨率方法。其中盐度场在加密网格上计算，其他物理量在基础网格上计算。基础网格和加密网格分别保证能够解析速度场黏性Kolmogorov尺度和盐度场Batchelor尺度，两者的比值为Sc1/2。速度场从基础网格向加密网格插值采用保证局部无散性的Hermite插值方法。利用该方法，能够有效提高处理高施密特数标量场的模拟效率。注意由于温度场扩散较快，在本文参数范围内其特征尺度一般大于速度场，因此温度场采用与速度场一致的网格分辨率。实际计算中考虑到并行效率，基础网格不能太小，因此加密网格的加密系数要小于Sc1/2，一般单个方向在2～4倍左右。

1.3 模拟参数范围

本研究选取四种不同物性参数组合，对应的普朗特数和施密特数分别为（Pr,Sc） = （7, 21）、（7, 70）、（7, 700）和（700, 2100），其中（Pr,Sc） = （7, 700）对应海洋真实物性参数，（700, 2100）对应实验中使用的盐度-糖度体系[10-11]，（7, 21）为直接数值模拟中常用的数值[14]。为系统研究Le的影响，此处额外计算组合（7,70）。

对于每个Pr和Sc的组合，计算三个不同密度比，Λ= 0.8、1.2、1.6。对于每个密度比，Ra= 1×106～1×1010，计算5个不同瑞利数。总计60个算例。图1给出了Ra-Λ相平面上的模拟参数范围。

图1 Ra-Λ相平面上模拟的参数范围Fig. 1 Simulated parameters on the Ra-Λ plane

2 主要结果与分析

本节首先讨论温盐传输特性随控制参数的变化行为以及相应的幂律模型，然后讨论流动结构及相应特征尺度的行为。

2.1 温盐传输特性

热量和盐度的传输率可由努塞尔特数度量。本文定义努塞尔特数为对流传输率与传导传输率的比值；而湍流运动强度可由均方根速度定义的雷诺数度量。这三个无量纲参数的数学定义如下：

其中，w为垂向速度，尖括号代表全场空间平均和时间平均；urms为速度矢量大小的全场时空均方根值。注意公式（10、11）所定义Nu只考虑对流运动引起的传输率，而没有考虑传导引起的传输。图2～图4展示了上述温盐努塞尔特数和雷诺数随流动控制参数的变化行为。可以看出，对四种不同物性参数组合，三个主要输运参数都随着Ra的增大而增加，且分别满足类似的幂律规律。同时，密度比Λ增加导致努塞尔特数和雷诺数减小；这是因为密度比增加对应致稳温度梯度的相对强度增加，从而导致流动和传输减弱。但是Λ的影响幅度随着Le的增加而降低。尤其是盐度传输率，当Le= 100时NuS几乎不随Λ变化，见图3（c）。特别地，对于不同物性参数组合，努塞尔特数的取值与Le的大小相关，而与Pr和Sc的取值无关。如图2（a、d）和图3（a、d），其Le都为3，但是Pr和Sc各自相差100倍，温度和盐度努塞尔特数大小非常相似。但是随着Pr和Sc的增大，虽然Le保持不变，流动雷诺数显著降低。

图2 温度努塞尔特数随流动参数变化趋势Fig. 2 Heat Nusselt number versus control parameters

图3 盐度努塞尔特数随流动参数变化趋势Fig. 3 Salt Nusselt number versus control parameters

对图2～图4数据的进一步分析表明，两个努赛尔特数基本遵循Ra0.30的规律增加，而雷诺数近似按照Ra0.42幂律变化。由此可假设存在如下幂律标度：

图4 流动雷诺数随流动参数变化趋势Fig. 4 Reynolds number versus control parameters

其中Λ的幂指数取值与物性参数Pr和Sc相关，其数值可根据数据由最小二乘拟合给出。具体结果见表1。可以看到，对于海洋相关的物性参数，即Pr= 7和Sc= 700，盐度雷诺数对应的幂律指数ξ非常接近于零，即与Λ相关性很小。但是对于其他物性参数组合，如盐-糖体系等具有较小Le的情况，传输率受到Λ的显著影响。图5中展示了用相应密度比幂律修正的响应参数随瑞利数的变化规律，可以看到对每组物性参数，数据均重合在一起，表明公式（13～15）的幂律模型确实可以描述响应的数据。

图5 用密度比幂律修正后的相应参数随瑞利数变化趋势Fig. 5 Nusselt numbers and Reynolds numbers rescaled by the power laws of density ratio versus Rayleigh number

表1 不同Pr和Sc对应的Λ幂律指数Table 1 Values of Λ exponents for different combinations of Pr and Sc

2.2 盐指结构

本节讨论典型流动结构和相关特征尺度。对于盐指型双扩散湍流，主要的典型流动结构为盐指结构，即沿垂直方向发展的细长对流胞结构，其水平宽度可远小于竖直高度，如图6所示。

图6 典型盐指结构对应的垂向速度场云图（Pr = 7，Sc = 700，Ra = 1×1010，Λ = 1.2）Fig. 6 Typical finger structures illustrated by the contours of vertical velocity（ Pr = 7, Sc = 700,Ra = 1×1010, and Λ = 1.2）

为提取如图6所示盐指结构的水平宽度，我们在流场中间高度为0.5H的区域内计算垂向速度w的水平自相关系数C(δ)，并以C(δ)第一个零点位置对应δ的两倍作为盐指宽度df。图7中给出盐指宽度df随流动参数的变化趋势。如图所示，对于大部分算例，盐指宽度均满足幂律df～Ra−1/4，该幂律指数与线性稳定性分析给出的结果一致[6]。但是，对于小Le情况，当Ra足够大时，盐指宽度明显偏离该幂律。如图8所示，此时流场中主导结构不再是盐指结构，而是较宽的对流涡结构，因而不再满足盐指结构的幂律规律。注意偏离幂律规律的算例发生在Λ= 0.8，即整体密度梯度处于不稳定分层状态，与经典对流流动类似。对于最大的Le= 100，所有算例结果都满足该幂律。当流场主导结构为盐指时，盐指宽度随着Λ的增加而减小。

图7 盐指宽度随流动参数的变化趋势Fig. 7 Finger width versus control parameters

图8 典型对流涡结构对应的垂向速度场云图（Pr = 7，Sc = 21，Ra = 1×1010，Λ = 0.8）Fig. 8 Typical finger structures illustrated by the contours of vertical velocity（ Pr = 7, Sc = 21, Ra = 1×1010, and Λ = 0.8）

3 结 论

综上所述，本文针对四种不同物性参数组合的双扩散系统进行不同瑞利数和密度比的二维直接数值模拟研究。物性参数组合包括海洋真实物性参数，实验和数值模拟常用参数。模拟瑞利数范围涵盖4个数量级，同时密度比与海洋典型参数相当。

结果表明，热流、盐流和流动雷诺数随瑞利数的整体变化趋势不随物性参数的改变发生显著变化，但其具体数值受密度比的影响。密度比的影响幅度与施密特数和普朗特数的比值，即刘易斯数有显著的关系。对于小刘易斯数，密度比的增大导致努塞尔特数和雷诺数的减小。随着刘易斯数的增加，努塞尔特数和雷诺数的减小幅度下降。对于海水对应的刘易斯数近似为100的情况，温度努塞尔特数和雷诺数仍然随密度比的增大出现明显下降，但是盐度努塞尔特数几乎不随密度比的增大发生变化。基于数据，本研究给出了温盐努塞尔特数和雷诺数随瑞利数和密度比变化的幂律模型描述。并确定了几种实验和计算中常用的双扩散系统的幂律指数。

对于本文中的大多数算例，流动主导结构为沿垂向发展的细长盐指结构，其竖直高度远大于水平宽度。对于小刘易斯数和大瑞利数情况，主导结构变为类似对流涡结构。当主导结构为盐指时，盐指宽度随瑞利数的变化与线性稳定性分析给出的结果类似。对于固定瑞利数，密度比的增加导致盐指宽度减小。

本文给出的结果，如物性参数组合的影响、传输特性的幂律模型、主导结构及尺度的变化趋势，对理解海洋盐指型界面传输特性有重要的价值。同时，不同物性参数结果的对比，对解读实验和模拟数据，以及将实验和模拟结果应用于海洋参数时，有相应的指导意义。