海洋盐指演化过程的数值模拟❋

徐 岩， 刘香全,3,4， 宋仁刚， 岑显荣， 鲁远征， 郭双喜， 周生启

(1.山东科技大学电子通信与物理学院，山东 青岛 266590； 2.中国科学院南海海洋研究所热带海洋环境国家重点实验室，广东 广州 510301；3.中国科学院半导体研究所集成光电子学国家重点实验室，北京 100083； 4.中国科学院大学材料科学与光电技术学院，北京 100049)

双扩散现象是自然界及工业生产中普遍存在的现象，如大气中的热和污染物扩散，海洋中的扩散对流及盐指现象，工业生产中的晶体生长、金属固化等[1]。当流体中的两种或两种以上组分具有不同的扩散速率时，就会产生双扩散现象。在海洋环境中，海水的热扩散率比盐扩散率大两个数量级，当垂向温度和盐度分布对海水的密度梯度具有相反的作用时，在热盐驱动下易形成双扩散对流。双扩散混合是海洋中物质和能量进行垂向输运的重要机制之一。

在亚热带海洋，由于强烈的太阳辐射的作用，使得表层海水的温度较高且蒸发量较大，形成了表层海水高温高盐、次表层海水相对低温低盐的状态，在这种状态下，温度场是稳定的，而盐度场却不稳定[2]。如图1(a)所示，当表层海水是高温高盐的均匀层，次表层海水是低温低盐的均匀层时，若存在一个微小的扰动，使得次表层的低盐海水微团跨过交界面向上运动，由于热扩散速率远大于盐扩散速率，该海水微团就会迅速的在周围吸收热量达到与周围海水微团相同的温度，但此时仍保持原来的盐度不变，这样浮力作用大于其重力的作用，促进其向上运动。同理，对于表层的海水微团跨越交界面向下的扰动，则能促进其向下的运动。可见，在此过程中存在着正反馈机制，使该微团在交界面处于一个相当不稳定的状态，进而产生盐度较大的海水向下、盐度较小海水向上，呈指状分布的现象，即盐指现象[3-11]。

图1 盐指形成机制示意图(a)及海洋盐指模型(b)

在有关海洋盐指的研究中，盐指导致的热盐输运问题受到海洋学家们的广泛关注。盐指中的热盐输运能改变水体的微观结构，引发水体有效位能的下降。特别地，在调节亚热带流涡内区的主温跃层和环流结构方面，从盐指释放出来的能量可以起到决定性的作用。但是，盐指包含了复杂的小尺度动力学过程，这些过程在大尺度环流的数值模式中是难以分辨的，因此很难直接评估盐指的具体贡献。过去，由于海洋观测手段的限制，对于盐指热盐通量的计算主要依赖从室内实验得到的经验公式，例如Turner[12]提出的4/3标度律：FS=Af(Rρ)(ΔS)4/3。Turner的公式表明，盐通量FS由相邻两混合层之间的盐度差ΔS以及密度比Rρ所决定。如今，随着计算机仿真技术的进步，数值模拟逐渐被应用于盐指现象的研究。众多学者对盐指问题作出了重要的贡献，如Yang等[13-16]通过直接数值模拟方法研究了盐指的热盐输运规律和流动结构；李文婷和李永放[3]数值模拟了两层和多层结构中的盐指现象，发现瑞利数(RT、RS)、初始密度稳定率(Rρ)均影响交界面处产生的盐指现象；罗莹莹等[9]利用基于杂交网格的高精度数值方法研究了多孔介质中的盐指现象，分析了孔隙率对盐指传热传质效应的影响; 在现场观测方面，鲁远征等[17]结合微结构数据对南海的盐指及湍流混合过程进行了分析。

过去有关盐指的模拟工作大多是在固定的热瑞利数或盐瑞利数的条件下进行的，鲜有关注盐指在不同的瑞利数组合条件下的演化过程。本文在前人研究的基础上，通过数值方法模拟了盐指的生成和演化过程。通过改变不同的初始条件，模拟分析了不同粘滞系数、热瑞利数和盐瑞利数等因素对盐指形成和演化过程现象的影响。

1 模型建立

1.1 理论基础

用数值计算的手段对盐指现象进行过程模拟具有一定优势，既能够精确地描述流量场及标量场，也能够严格地控制边界条件。1973年，Turner[12]指出盐指双扩散的稳定性可以用密度比进行衡量，其定义如下：

(1)

式中:α=-(1/ρ)∂ρ/∂T和β=-(1/ρ)∂ρ/∂S分别为热膨胀系数和盐收缩系数；ΔT、ΔS分别为盐指结构中上下两层溶液的温度差和盐度差。

到目前为止，对双扩散的理论研究，Boussinesq近似[13,18]得到了广泛的认可。Boussinesq近似只考虑密度对质量力项的影响，密度的线性状态方程为：

(2)

式中:T0、S0、ρ0分别代表参考温度、参考盐度及参考密度；ρ为流体密度；T和S分别为流场的温度和盐度。

在Boussinesq近似假设下，盐指过程的控制方程形式如下：

(3)

为了更好地研究双扩散情况，能够确定可以控制流体流动形态及强度的无量纲参数是非常有必要的。 因此在模拟盐指时，引入了热瑞利数和盐瑞利数[3]，其定义为：

(4)

式中L为流场的特征长度。

1.2 方腔模型的建立

建立盐指模拟的二维方腔模型，如图1(b)所示。设定其长和宽均为L=0.10 m，在直线Y=0.05 m的位置分成上下均匀的两层，上层为高温高盐海水，温度和盐度分别用Ta、Sa表示；下层为低温低盐海水，温度和盐度分别用Tb、Sb表示(Ta>Tb,Sa>Sb)。

ΔT=Ta-Tb, ΔS=Sa-Sb。

(5)

控制方程采用有限体积法求解，从而计算出流场的温度、盐度、速度随时间的演变关系。计算在方形网格上进行，网格数为100×100。为了满足CFL条件，计算的时间步长为0.01 s，每计算1 000步(10 s)输出一组数据。在盐指模拟中，首先设置二维方腔中流体的温度为参考温度(300 K)，浓度为参考浓度(17.0 g/kg)，然后按照表达式(5)及表达式(6)设置上下两层的温度差及浓度差。

(6)

为分析不同参数对盐指现象的影响，分别设计了一系列的参数。为了在较大范围内改变瑞利数，首先通过改变粘滞系数来进行，分别模拟了粘滞系数为1.0×10-7m2·s-1、1.0×10-6m2·s-1、1.0×10-5m2·s-1的情况。然后，在固定粘滞系数为1.0×10-6m2·s-1的情况下，通过改变温度差ΔT和盐度差ΔS以达到改变热瑞利数和盐瑞利数的目的。具体参数设置如表1、2所示。

表1 热瑞利数的参数设置(ΔS=0.1 g/kg)

表2 盐瑞利数的参数设置(ΔT=1 K)

2 结果处理分析

在通过改变粘滞系数η、温度差ΔT和盐度差ΔS来改变瑞利数时，发现在不同的参数下，盐指模拟结果存在着很大的差异，这种差异主要表现在盐指的生成速率及盐指的数目上。为了能够较为准确的得到观察结果，规定形成盐指的标志是大部分指状结构的高度能够达到标准长度(直线Y=0.04 m与直线Y=0.06 m之间的长度)，盐指的数目N是指向上和向下扩散形成的指状结构的总个数，其最小的计数单位是0.5。上下边界处长度达到标准长度，但是宽度为其他指状结构的一半的计为0.5；当大部分指状结构的长度达到标准长度，但某些指状结构未达到标准长度的计为0.5[3]。

2.1 盐指现象模拟情况

取η=1.0×10-6m2·s-1,ΔT=1 K,ΔS=0.1 g/kg，即RT=1.40×107,RS=5.25×106情况下，对盐指的演化过程进行展示。当时间约为150 s时达到标准长度，可以看到当上层高温高盐流体下沉到下层低温低盐流体时，盐度基本保持不变，而温度却迅速衰减(见图2，3)，这一现象与理论(热扩散率比盐扩散率大两个量级)相互印证。在接下来的时间里，指状结构的长度不断增加，并最终延伸到上下边界。

观察中截线(X=0.05 m)上的温度和盐度的变化曲线，在整个过程中温度由上层高温、下层低温的初始状态逐渐趋于均匀，前段时期(约2 260 s之前)温度变化曲线波动较大，在此之后的温度曲线的波动较小，上下层温差逐渐减小，达到一个相对均匀的状态，但是后段时期历时较长。对于盐度曲线而言，模拟开始时由上层高盐下层低盐的初始状态迅速发生振荡，并在约3 380 s时达到一个基本线性分布的状态，在之后的时间里，盐度曲线逐渐向线性分布上修正，当模拟时长为5 h时，该直线上温度和盐度曲线如图4所示。

对于模拟过程中直线上温度和盐度曲线变化的原因给出如下解释：对于盐指现象的初始条件而言，温度分布是处于稳定状态的，而盐度分布却是处于非稳定状态的，在有扰动的条件下，上层水团向下运动，但盐度的扩散速率远小于热扩散率，从而使得密度相对周围的密度大，在重力作用下向下掉，并携带着热量，但当盐度基本达到一个稳定状态时，热量的传递仅仅是依靠分子热传导，因此在此之后上下层温差的减小就耗时较长。

2.2 不同粘滞系数对盐指现象的影响

粘滞系数是与雷诺数、瑞利数相关的参量，且成反比关系，通过改变粘滞系数可以改变这两个参量。模拟在固定ΔT=1 K,ΔS=0.1 g/kg的情况下，分别设置粘滞系数为1.0×10-7m2·s-1、1.0×10-6m2·s-1、1.0×10-5m2·s-1。当盐指达到标准长度时，所对应的结果分别如图3、5、6所示。

通过对比发现，随粘滞系数的增大，对应的雷诺数、瑞利数相应地减小，盐指的生长速度明显变慢，盐指的数目变小且宽度变大。由于流体本身具有粘性，随着粘滞系数的增加，流体更为粘稠，流动更为缓慢，致使出现了这种现象。

2.3 热瑞利数对盐指现象的影响

如表1的参数设置所示，固定η=1.0×10-6m2·s-1，ΔS=0.1 g/kg，即RS=5.25×106，在这种条件下改变盐指模型上下两层的温度差ΔT，以改变热瑞利数RT。ΔT取0.8、1.0、1.2、1.4和1.6 K时，分别观察盐指现象生成的时间及盐指数目的多少。温度图的指状分布可以看成是由盐度的扩散所引起的扰动，是物质携带热量的结果，因此，通常意义上的盐指现象是盐度在溶液中的指状分布。在接下来的讨论中，采用盐度分布图来观察盐指现象，ΔT=1.0 K的情况如图3所示，其余模拟结果如图7所示。盐指现象生成的时间及数目如表3所示。

图2 盐指模拟100 s后的形态

图3 盐指模拟150 s后的形态

图4 盐指模拟5 h后直线x=0.05 m上的温度(a)、盐度(b)曲线图

由表3可以看出，随着热瑞利数的增加(此时初始密度稳定率Rρ0增加)，产生明显盐指现象的时间滞后，盐指的数目增多。时间的滞后意味着盐指现象生成的速率在降低，给出如下解释：由于上下两层温差的增大，热溶液向上流动的趋势对盐水向下流动的阻碍增强，致使盐指现象的生成速率减慢。对于盐指数目的增多，是由于温度差的增加使热瑞利数变大，单个指状结构的宽度变小，盐指的数目增多。

(a)温度分布图；(b)盐度分布图

图6 粘滞系数η=1.0×10-5m2·s-1时的温度(a)、盐度(b)图

图7 上下两层温差不同(热瑞利数不同)情况下，盐指达到标准长度时的盐度分布图

表3 不同热瑞利数情况下出现明显盐 指现象的时间及盐指的数目(ΔS=0.1 g/kg)

2.4 盐瑞利数对盐指现象的影响

如表2的参数设置所示，固定η=1.0×10-6m2·s-1,ΔT=1 K，即RT=1.40×107，在这种条件下改变盐指模型上下两层的盐度差ΔS，以改变盐瑞利数RS。ΔS取0.06、0.08、0.10、0.12、0.14 g/kg时，分别观察盐指现象生成的时间及盐指数目的多少。模拟结果如图8所示。盐指现象生成的时间及数目如表4所示。

表4 不同盐瑞利数情况下出现明显盐 指现象的时间及盐指的数目(Δ T=1 K)Table 4 The time for the appearance of patchy salt fingers and the number of the salt finger with various salinity Rayleigh number (Δ T=1 K)

图8 上下两层盐度差不同(盐瑞利数不同)情况下，盐指达到标准长度时的盐度分布图

从表4中的数据可以看出，在保持热瑞利数不变的情况下，盐瑞利数增加时(此时初始密度稳定率Rρ0减小)，产生明显盐指现象所需的时间显著减少，生成盐指的数目变化却不是很明显，但是缩短的时间跨度逐渐降低，这说明盐指的生成速度虽在增加，但其增加程度是不同的，且是逐渐变慢的。盐指生成速度的增加可以理解为：在上下两层温度差不变的情况下，热效应对盐指向下扩散的阻碍相同，随着盐度差的增加，盐度的不稳定性逐渐增大，向下扩散的趋势逐渐增强，因此盐指的生成速度不断增加。与热瑞利数对盐指数目影响不同的是，盐瑞利数的改变对盐指数目的影响较小，这是海水的热扩散率和盐扩散率不同造成的(热扩散率比盐扩散率大两个量级)，因此热瑞利数能够较大地影响盐指的数目。

3 结语

本文以海洋双扩散对流为背景，在二维方腔模型中对盐指现象进行了模拟研究，展示并解释了盐指的模拟结果。由于正反馈的存在，盐指现象发生的时间很短，达到相对稳定状态之后，主要依靠热扩散进行能量的输运，耗时较长。

通过改变粘滞系数来改变流体的雷诺数、瑞利数，对比观察，得出随粘滞系数的增加(即雷诺数、瑞利数减小)，盐指的演变速度变慢、盐指生成数目减少且宽度增加。建立统一的标准，对不同热瑞利数、盐瑞利数情况下生成明显盐指现象的时间以及盐指的数目进行统计，并对结果进行了对比分析，从而得出了这两个参数对盐指现象的影响效果。研究发现当热瑞利数RT增加时，盐指的生成速度会变慢，盐指的数目会增多；而当盐瑞利数RS增加时，盐指的生成速度会变快。