基于二维水动力模型的密度驱动流模拟

01 研究背景

气候变化带来的海平面上升可能给河流的入海口带来巨大的改变。原本盐度较低的河水和盐度较高的海水之间维持的盐平衡可能被打破。盐度不同带来的密度差异将驱动海水进一步入侵到陆地的河流体系中。针对入海口处的仿真计算必须考虑这一现象带来的影响。

在二维模型中，密度在垂向上的分布是假设均匀的，不能体现入海口处的密度分层现象，只能考虑密度在水平方向上分布不均引起的效应。为了理解和评估二维模型模拟密度驱动流的效果，该案例进行了三次仿真测试。其中两个为矩形河道和梯形河道的理论情形算例，余下一个是地中海尼罗（Nile）河入海口的实际情形算例。这些算例的计算都将在通用水动力仿真软件里进行。

图 1 河流入海口处的盐度空间变化示意图

02 测试算例及结果讨论

 

每个算例都选择了四个工况进行计算，分别是：

1. A.    忽略扩散，只考虑水平方向上的密度差异的静止水

2. B.    忽略水平方向上的密度差异，只考虑扩散的静止水

3. C.    同时考虑扩散和水平方向上的密度差异的静止水

4. D.    同时考虑扩散和水平方向上的密度差异的流动水

通过这四个工况的计算，我们可以评估密度差异对结果的影响，并与扩散带来的影响进行对比。在尼罗河入海口算例中，流动的水更符合真实情况，因此静止水都将替换成流动水，只存在三个工况。

1. 矩形河道

矩形河道是横截面为矩形，深5 m，宽200 m，长1000 m的河道。底部高程设置为0。整体网格由10 m的三角形网格构成，在中间区域进行了加密（与初始盐度分布的设定相匹配），加密后的网格尺寸为4 m。节点数为3303，单元数为6344。

图2 矩形河道的网格

该河道拥有两个开放边界。上游边界在河流入流方向，为固定流量边界，静止水和流动水两种情况下的流量分别是0和80 m3/s，该边界上的盐度为零。下游边界在近海方向，为固定水位边界，水位5 m，盐度35 mg/l。曼宁摩擦系数选取0.022。粘度为0.001 m2/s。

在靠近下游边界的一半区域内初始盐度为35 mg/l，在另一半区域里与上游边界的盐度一致。时间步为5 s，仿真时长为100 days。

图3 矩形河道算例的初始状态盐度分布

矩形河道算例的计算结果如图4所示。对比工况A和工况B，可以发现扩散的影响要比密度驱动的小。进一步增加10 days的仿真时长，工况B（只有扩散）的盐度分布几乎没有变化。

图 4 四种工况下的矩形河道算例结果

（工况A：只考虑扩散，静止水；

工况B：只考虑密度驱动，静止水；

工况C：同时考虑扩散和密度驱动，静止水；

工况D：同时考虑扩散和密度驱动，流动水）

2.  梯形河道

梯形河道的深度、宽度、长度与矩形河道一致，其横截面如下图所示。

图5 梯形河道的横截面

网格参数的设置，以及其余数值和物理参数的设置与矩形河道算例保持一致。节点数为3325，单元数为6388。网格和底部高程如下图所示。

图6 梯形河道的网格

此时，如图7所示，四种工况下的盐度分布变化速度比矩形河道算例要快很多。只经历10天的仿真时间后，就能看到很明显的差异。

在工况A（只有密度驱动），基本上整条河道都变成了盐水区，这表明盐度上的差异会驱使下游的含盐水入侵到原本的淡水区。另外，盐度随深度的变化也体现了出来，高盐度的水集中在深水区。

在河道边缘，深度的变化将导致底部摩擦和湍流的改变。同时，在重力的影响下，密度驱动流也更倾向于流向深水区。除此之外，密度梯度引起的驱动压随水深的提高而提高，故而深水区的流动要比浅水区更加剧烈一些。

在工况B中，扩散的影响仍然很有限，只局限在中间区域。

对比工况C与工况A，添加扩散与否的区别也比较小。比较工况D下的矩形河道算例和梯形河道算例，可以看到，在上游有流量时，矩形河道的盐度边界与梯形河道的有较大不同，这反映了地形对盐度输运的影响。

图 7四种工况下的梯形河道算例结果

（工况A：只考虑扩散，静止水；

工况B：只考虑密度驱动，静止水；

工况C：同时考虑扩散和密度驱动，静止水；

工况D：同时考虑扩散和密度驱动，流动水）

3. 尼罗河入海口

该案例选择尼罗河的罗塞塔（Rosetta）分支进行研究。其范围为北纬31.32°到31.45°，东经30.34°到30.53°，总长约35 km。上游边界为流量边界，位于伊德费纳坝（Edfina barrage）处，下游边界为水位边界，位于地中海入海口。

图8 所研究的河流段的区域地图

罗塞塔河年平均径流为83.6 m3/s 下游的平均水位为0.37 m。该研究选取这两个平均值作为上游和下游的边界条件。

全河段的水深变化范围从2.30 m跨越到26.5 m，平均河宽约500 m。整体网格尺寸为40 m，在底部高程剧烈变化的和河宽较小的地方，则采用了20 m的网格。节点数为19448，单元数为36669。

图9 尼罗河入海口算例的局部网格

该研究进行了热启动计算。热启动的全场初始速度为零，初始水位为0.37 m，时间步为5 s，仿真时长为3 days。其结果作为输运模拟的初始状态。曼宁系数为0.022，粘度为0.01 m2/s。输运模拟的初始盐度的由下式计算：

其中Li 为含盐水侵入到河流的长度，Reo 为入海口处的雷诺数，Frdo 为考虑密度的弗劳德数。计算得到的侵入长度为15329 m。入海口处的盐度为38.5 kg/m3，假设盐度沿河道的变化是均匀的，那么初始盐度分布如下图所示。

图10 尼罗河入海口算例的初始盐度分布图

图11为仿真时间一天后的各工况盐度分布图。可以发现，含盐水的侵入长度在一天时间内很快变小了。但对比工况A和B，工况A下仍然保持着一个更大的入侵深度。这说明密度差异带来的流动更能抵抗上游来水的影响。对于工况C，其结果与工况A的比较相似。这些分析与前两个理论情形算例得到的结论是一致的。

图 11 仿真时间一天后的各工况盐度分布图。

（工况A：只考虑密度驱动；

工况B：只考虑扩散；

工况C：同时考虑扩散和密度驱动）

03 案例总结

 

该案例使用二维水动力数值仿真研究了二维模型下密度分布不均带来的影响，对比了扩散和密度驱动两种作用的效果，发现密度驱动带来的效应比扩散作用更加明显。

对密度驱动来说，河道地形也会对结果造成影响。在该案例中，梯形河道的盐度分布变化速度要比矩形河道更快。而对于扩散而言，两种河道的差异很小。此外，深水区的盐度要比浅水区要高。在尼罗河入海口算例中，密度差异带来的流动有效地抵抗了上游来水的影响，使得含盐水的侵入长度没有出现大幅的减小。

04 小结

 

在气候变化的大背景下，海平面上升已成事实，入海口处的盐平衡将被打破，或对这些区域的生态等方面造成较大的影响。该案例展现了通用二维水动力软件在模拟存在盐度分布差异的水动力仿真能力，验证了二维水动力在评估海平面上升对入海口区域的影响的可能性。

文章来源：远算云仿真