[Ocean Modelling for Begineers] Ch4. Long Waves in a Channel

Ch4. Long Waves in a Channel

简介

本章主要介绍明渠中分层流体模拟。练习包括浅水表面波，风暴潮、内波和分层流体模拟。

4.1 有限差分法详细介绍

4.1.1 泰勒公式

4.1.2 前差，后差，中心差分

4.1.3 二阶微分格式

4.1.4 截断误差

考虑函数：

[f(x) = A sin(2pi x/lambda) ]

方程的导数为：

[frac{df}{dx} = 2 pi A/ lambda cos(2 pi x/lambda) ]

如果我们使用中心差分近似一阶导数，可以得到

[frac{f(x+Delta x) - f(x-Delta x)}{2 Delta x} = frac{A sin[2 pi (x+Delta x) / lambda] - A sin[2 pi (x-Delta x) / lambda]}{2 Delta x} = 2 pi A / lambda cos(2 pi x / lambda) cdot [1-epsilon] ]

相对误差与绝对误差的比值

[epsilon(Delta x) = 1-frac{sin(2 pi Delta x/ lambda)}{2 pi Delta x / lambda} ]

使用有限差分法模拟波动过程时，只有波长范围内有10个以上节点时才能准确描述波动过程。

4.2 表面重力长波

4.2.1 单个过程分析

NS方程描述了在不同长度尺度和时间尺度上同时发生的流体运动过程。在某些假定下，我们可以取出单独过程进行研究。比如说，当周期远远小于惯性周期时，我们可以忽略科氏力的作用而对控制方程进行简化。

4.2.2 浅水过程

4.2.3 浅水模型

我们做以下假定：

波周期相比惯性周期很小，可以忽略科氏力作用
假设摩擦力为一阶近似
忽略非线性项，即波速远远超过颗粒运动速度

4.2.4 控制方程

4.2.5 波动解析解

4.2.6 动画演示

4.2.7 数值计算网格

使用交错网格，流速节点位于相邻水位节点中间。

4.2.8 有限差分格式

4.2.9 稳定准则

CFL准则：

[lambda = frac{Delta t}{Delta x}sqrt{g h_{max}} le 1 ]

其中(h_max)为模型中最大水深。

4.2.10 一阶Shapiro过滤器

为了去掉数值振荡现象，使用Shapiro过滤算子（Shapiro, 1970）：

[eta\_k^{n+1} = (1-epsilon)eta\_k^* + 0.5 epsilon (eta\_{k-1}^* + eta\_{k+1}^*) ]

其中(epsilon)为光滑系数。

4.2.11 陆地与海岸边界

没有流动可以穿过陆地与海岸线，除非特殊处理的计算格式（4.4节介绍）。由于交叉网格设置，海岸线需要特殊的条件，若(h_{k+1}le 0)时(u_k)为0。

4.2.12 侧向陆边界条件

模型计算网格从(k=1)到(k=nx)，边界处单元(k=0)与(k=(n+1)x)也需要赋值才能计算。一种方法是令边界闭合，横向对流流速为0：

[u_0^n = 0 quad ext{and} quad u_{(n+1)x}^n = 0quad ]

0梯度边界条件用来消除物理量在边界处扩散流量，边界条件为：

[C_0^n = C_1^n quad ext{and} quad C_{(n+1)x}^n = C_{nx}^nquad ]

周期边界条件表示为：

[C_0^n = C_{(n+1)x}^n quad ext{and} quad C_{(n+1)x}^n = C_1^nquad ]

4.2.13 模块化Fortran脚本

4.2.14 Fortran代码结构

4.3 练习5：明渠中长波运动

4.3.1 目标

4.3.2 说明

水渠长1km，网格步长10m，设置101个节点单元，包括最后一个边界节点，计算时间步长0.1s，满足CFL准则。

溃坝算例
平衡状态水深为10m，在中心110m宽长度上初始高程提升1m
波动算例
在中点处放置一个造波板，以1m振幅，20s周期造波

4.3.3 简单代码及动画示例

4.3.4 结果

4.4 练习6：干湿算法

4.4.1 目标

可以模拟洪水传播到陆地过程。例如由潮水或风暴潮淹没海岸情形。

4.4.2 重新定义干湿

湿节点定义为水深超过某一阀值(h_{min})，通常为几厘米左右。这个阀值作用主要是避免水体退去后湿单元出现负水深，从而使模型崩溃。干单元通常定义为(hle h_{min})

4.4.3 淹没干单元

4.4.4 淹没边坡

(h_0)为静水位置距底部距离，在有水地方为正值，而在陆地为负值
初始高程在湿单元处为0，在干单元与陆地高程相同(eta = -min(0, h_0))
真正水深为(h = h_0 + eta)，这个水深用来判别干湿节点

4.4.5 终极测试（？Ultimate Crash Tests）

两组测试算例，长1km，网格步长10m，设置101个节点单元，前后增加一个边界节点，计算时间步长0.1s，满足CFL准则。

小岛最高位于水面0.5m，左侧200m有一风暴潮，初始高程1m，两侧为陆边界
斜坡上有一高程为4m水坑

代码几处要点

过滤算子只在湿节点作用，跳过干节点
流速更新只考虑两种情形，两侧为湿节点；由湿节点流向干节点
干湿阀值水深为0.1m，水深越小越精确（随着阀值减小，水体下流速度变快？）。但是此值不能为0，否则水体前进速度会变得特别快，并不符合实际。

存在问题

干湿水深阀值如何确定？

4.5 多层浅水模型

4.5.1 基础

在浅水假定条件下，可以构造每一层都是密度为常数多层浅水模型。第i层的动量方程为

[frac{partial u_i}{partial t} = -frac{1}{ ho_i} frac{partial P_i}{partial x} ]

其中，(i)为层数，每层水体动压为

[egin{array}{ll} P_1 = ho_1 g eta_1 cr P_2 = P_1 + ( ho_2 - ho_1) g eta_1 cr cdots end{array}]

更一般形式为

[P_i = P_{i-1} + ( ho_i - ho_{i-1})g eta_i quad ext{for} quad i = 1,2,3,... ]

连续方程为

[frac{partial h_i}{partial t} = - frac{partial u_i h_i}{partial x} ]

其中第i层水体厚度为

[h_i = h_{i, 0} + eta_i - eta_{i+1} ]

这里(h_{i, 0})为第i层水体未扰动时的厚度。

4.6 练习7：长波在分层流体中运动

4.6.1 目标

模拟重力波在不同密度分层水体内运动情况

4.6.2 算例描述

顶层水体密度为(1025 kg/m^3)，第二层至底层密度分别从(1026 kg/m^3)到(1026.5 kg/m^3)逐渐增加。模型由左侧正弦振荡自由表面驱动，振幅为1s，周期分别为10s与2h。两侧为封闭的固边界。计算时间取10倍波动周期，时间步0.25s。

4.6.3 代码

4.6.4 结果

4.6.5 内波相速度

在两层流体模型中，内部长波的相速度为：

[c_{iw} = sqrt{g'h^*} ]

这里(g')为相对重力加速度，(h^* = h_1 h_2 /(h_1 + h_2))为相对水深，对(h_2gg h_1)时，(h*approx h_1)。因此内部波运动速度远小于表面重力波，与表面重力波一样，在静压模型中，无法模拟内波破碎现象。

4.6.6 封闭水体自由振荡

在波节点处，流体只进行水平运动而没有垂向运动。与之相反，在波腹处流体只有垂向运动。