水波与淤泥质底床相互作用数值模拟——基于滞后回路特征的淤泥体流变模型

摘要

应用一个滞后回路特征的半经验的淤泥体流变模型描述表面水波作用下底床软泥层的振荡运动，建立了表面水波与淤泥质底床相互作用的垂向二维耦合模型。在数值处理上，为提高数值格式对水－泥界面附近流动的分辨力，较好地再现界面附近速度梯度大的流动特点，采用垂向对数网格处理技巧。模型验证是通过计算水波波高衰减率和底床软泥体输移速度，并与有关的实验测量数据进行比较。

正文

《1 前言》

1 前言

以粘性细颗粒泥沙为主形成的淤泥质河口海岸遍及世界各地。例如, 南美洲亚马孙河口海岸, 美国密西西比河口海岸, 流经印度和孟加拉国的恒河和不拉马普特拉河的河口海岸, 中国主要河流 (黄河、长江、珠江、海河和辽河) 的入海河口及其附近海岸。在风浪、水流等海洋水动力条件下海岸带底床淤泥体易发生输移。但在不同的水动力条件下, 淤泥体输移的主要方式往往不同:在强流弱波条件下以悬沙输移为主;在强波弱流条件下则以底部浮泥层或底泥层的输移为主。由于粘性泥沙对重金属和化学物质的吸附作用, 一方面, 底泥输移对污染物质扩散以及水质和生物生境有重要影响;另一方面, 底泥 (沙) 大量输移常成为风暴期间港口、航道骤淤的重要原因。因而, 水波与淤泥质底床相互作用问题受到国内外学者的普遍关注。

水波在淤泥质海床上传播时, 波浪运动的水体通过向下传递压力波能引起底床浮泥层产生振荡, 导致泥体向波浪传播方向显著地输移和水波波高沿程显著地衰减 ^[1,2,3,4]。

Gade开展了最早的淤泥质底床上浅水波的理论研究, 建立了泥-水两层流体模型, 研究了底泥层对上覆水体波浪运动的影响, 认为波浪大幅衰减的主要原因是泥层运动造成了波能损耗 ^[3]。从此, 以建立表达淤泥动力学特征的本构关系模型为重点的波浪与淤泥质底床相互作用的理论工作大量涌现。Dalrymple等用粘性流体模型表达泥层 ^[5];Mallard等较早提出了淤泥流变学的弹性体模型 ^[6];Bingham流体模型 (粘-塑性体模型) 由Krone提出 ^[7,8];Macpherson等提出了粘-弹性体模型 (Voigt模型) ^[9,10,11];练继建等提出了淤泥流变学的非线性粘弹性体模型 ^[12];Trien 等提出了粘-弹-塑性体模型 ^[13]。最近, Jiang提出了一个半经验的强非线性的流变响应模型 ^[14], 其应力-应变 (及应变率) 关系具有显著的滞后回线特征。尽管淤泥的流变行为相当复杂, 粘-弹-塑性模型也还不是最终的选择, 但是, 这种模型能够用于波浪-淤泥质底床相互作用的研究。

笔者应用半经验的强非线性的流变响应模型, 建立了模拟波浪运动的水体与其底床软泥层之间相互作用的垂向二维数值模型。为了提高水-泥界面附近流动的分辨力, 更好地再现界面附近速度梯度大的流动特点, 计算中采用垂向对数网格的数值处理技巧, 并利用有关的实验室数据进行了验证。

《2 波浪作用下高岭土软泥的经验流变模型》

2 波浪作用下高岭土软泥的经验流变模型

Jiang根据对商业高岭土泥浆体受循环剪切荷载实验数据, 提出了一个经验流变模型 ^[14], 其本构方程为:

$τ = \frac{G_{0} γ}{1 - α | γ |} + \frac{μ_{0} \dot{γ}}{1 + β | \dot{γ} |} (1)$ $τ = \frac{G_{0} γ}{1 - α | γ |} + \frac{μ_{0} \dot{γ}}{1 + β | \dot{γ} |} (1)$

式中, τ, γ和 $\dot{γ}$ $\dot{γ}$ 分别代表剪应力、剪应变及其应变率。G₀, μ₀, α 和β是模型参数。其中, G₀代表初始剪切模量 (γ=0) ;μ₀代表初始粘度 $(\dot{γ} = 0)$ $(\dot{γ} = 0)$ ;α和β是确定骨架曲线形状的两个系数。这些模型参数的经验表达式分别为:

$\begin{array}{l} G_{0} = (5.560 42 \times 10^{5} \cdot W^{- 2.766 31}) \cdot {[- 7.255 26 - \\ 558.832 42 \cdot (1 + \tanh (- 0.578 46 \cdot Τ))] + \\ [81.693 10 + 707.435 01 \cdot (1.0 - \\ \tanh (- 0.471 94 \cdot Τ))] \cdot α^{\frac{1.714}{{1.0 + 9.872 \cdot \exp (- 0.866 \cdot Τ)}}}} (2) \\ μ_{0} = {5.560 42 \times 10^{5} \cdot W^{- 2.766 31}} \cdot \\ {29.033 6 \cdot (1.31^{Τ})} \cdot \ln {[24.714 \cdot \\ \exp (- 0.308 \cdot Τ)] \cdot β + 1.0} (3) \\ α = \frac{1}{γ_{\max}} {1.0 - (0.790 51 - 5.045 71 \times 10^{- 3} \cdot Τ) \cdot \end{array}$ $\begin{array}{l} G_{0} = (5.560 42 \times 10^{5} \cdot W^{- 2.766 31}) \cdot {[- 7.255 26 - \\ 558.832 42 \cdot (1 + \tanh (- 0.578 46 \cdot Τ))] + \\ [81.693 10 + 707.435 01 \cdot (1.0 - \\ \tanh (- 0.471 94 \cdot Τ))] \cdot α^{\frac{1.714}{{1.0 + 9.872 \cdot \exp (- 0.866 \cdot Τ)}}}} (2) \\ μ_{0} = {5.560 42 \times 10^{5} \cdot W^{- 2.766 31}} \cdot \\ {29.033 6 \cdot (1.31^{Τ})} \cdot \ln {[24.714 \cdot \\ \exp (- 0.308 \cdot Τ)] \cdot β + 1.0} (3) \\ α = \frac{1}{γ_{\max}} {1.0 - (0.790 51 - 5.045 71 \times 10^{- 3} \cdot Τ) \cdot \end{array}$

$\begin{array}{l} \exp (- [0.749 46 \cdot \tanh (0.242 983 \cdot Τ)] \cdot γ_{\max})} (4) \\ β = (0.341 35 + 0.088 222 \cdot Τ + \\ 0.011 651 \cdot Τ^{2}) \cdot γ_{\max}^{(1.045 29 - 1.419 46 \cdot Τ^{0.174 215})} (5) \end{array}$ $\begin{array}{l} \exp (- [0.749 46 \cdot \tanh (0.242 983 \cdot Τ)] \cdot γ_{\max})} (4) \\ β = (0.341 35 + 0.088 222 \cdot Τ + \\ 0.011 651 \cdot Τ^{2}) \cdot γ_{\max}^{(1.045 29 - 1.419 46 \cdot Τ^{0.174 215})} (5) \end{array}$

式中, γ_max为最大应变, T为循环剪切荷载 (波浪) 的周期, W为泥体水含量比。

本构方程式 (1) 以及表达各模型参数G₀, μ₀, α 和β的经验公式 (2) - (5) 构成了一个完整的淤泥体流变关系模型。

图1 (a) 和 (b) 分别为根据该经验流变关系计算所得高岭土泥浆体在循环剪切外力作用下泥体中的剪应变, 及其应变率随剪应力的周期性变化而变化的关系。取泥体水含量比W=200%, 循环剪切荷载周期T=6.0 s。图中虚线OC和O′C′表示开始时剪应变随剪应力及剪应变率随剪应力的起始曲线。当剪应力开始减小时, 泥体中的剪应变及其应变率并不沿起始曲线返回, 而是滞后于应力变化, 这就是泥体内的剪应变 (及其应变率) 对循环剪应力变化响应的滞后现象。当剪应力减小为0时, 泥体内的剪应变及其应变率并不为0, 而是分别有一个剩余量γ₀和 $\dot{γ}$ $\dot{γ}$ ₀, 在这里称其为剩余剪应变和剩余剪应变率, 暂留在泥体中。这种受循环剪切外力作用 (波浪荷载) 淤泥体流变关系的滞后响应特征同时也显示了淤泥流变学的强非线性特征。

《图1》

图1剪应力与剪应变 (a) 及其应变率 (b) 的关系

Fig.1 Shear stress vs. shear strain (a) and rate of shear strain (b)

《3 水-泥两层模型运动控制方程》

3 水-泥两层模型运动控制方程

图2为水-泥两层模型定义草图。建立xoz坐标系, x轴水平, z轴向上为正, 以静止时水-泥界面为坐标原点。η和η_m分别代表自由水面和水-泥界面的垂直位移;d_w和d_m分别代表水层和泥层的厚度。波浪向x轴正向传播。假设水层为牛顿流体, 泥层为具有方程 (1) 本构关系的均质流体。于是, 在忽略非线性项的情况下, 这两层流体的运动可统一地用以下形式的Navier-Stokes方程和连续方程表达:

$\begin{array}{l} \frac{\partial u}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{1}{ρ} \frac{\partial τ_{z x}}{\partial z} (6) \\ \frac{\partial w}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial z} + \frac{1}{ρ} \frac{\partial τ_{x z}}{\partial x} (7) \\ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0 (8) \end{array}$ $\begin{array}{l} \frac{\partial u}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial x} + \frac{1}{ρ} \frac{\partial τ_{z x}}{\partial z} (6) \\ \frac{\partial w}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial z} + \frac{1}{ρ} \frac{\partial τ_{x z}}{\partial x} (7) \\ \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0 (8) \end{array}$

式中, p为动压力, u和w分别为速度的x和z分量, ρ为密度。

对于泥层, 通过将剪应变γ和剪应变率 $\dot{γ}$ $\dot{γ}$ 公式 (9a-9b) 代入流变模型方程 (1) 中, 方程 (6) 、 (7) 中剪应力τ_xz和τ_zx可由速度u和w表达,

$γ = \frac{\partial ξ}{\partial z} + \frac{\partial ζ}{\partial x} ； \dot{γ} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} (9 a - 9 b)$ $γ = \frac{\partial ξ}{\partial z} + \frac{\partial ζ}{\partial x} ； \dot{γ} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} (9 a - 9 b)$

式中, ξ和ζ分别代表泥层质点运动位移的x和z分量, 即

$ξ = \int_{0}^{t} u d t; ζ = \int_{0}^{t} w d t (10 a - 10 b)$ $ξ = \int_{0}^{t} u d t; ζ = \int_{0}^{t} w d t (10 a - 10 b)$

方程 (6) — (8) 构成变量u, w和p的封闭的运动方程组。它们有以下形式的解:

$\begin{array}{l} u (x, z, t) = \tilde{u} (z) e^{i (n x - σ t)}; w (x, z, t) = \tilde{w} (z) e^{i (n x - σ t)} ； \\ p (x, z, t) = \tilde{p} (z) e^{i (n x - σ t)}; η (x, t) = a e^{i (n x - σ t)} (11 a - 11 d) \end{array}$ $\begin{array}{l} u (x, z, t) = \tilde{u} (z) e^{i (n x - σ t)}; w (x, z, t) = \tilde{w} (z) e^{i (n x - σ t)} ； \\ p (x, z, t) = \tilde{p} (z) e^{i (n x - σ t)}; η (x, t) = a e^{i (n x - σ t)} (11 a - 11 d) \end{array}$

式中, n是复波数, 其表达式为, n=n_r+iδ;通过其实部n_r可得到波长L=2π/n_r;以其虚部δ代表波浪衰减率。a是波动振幅。

将方程 (11a-11d) 带入方程 (6) — (8) 中, 由方程 (8) 可得

$\tilde{u} = \frac{i {\tilde{w}}^{'}}{n} (12)$ $\tilde{u} = \frac{i {\tilde{w}}^{'}}{n} (12)$

将式 (12) 带入方程 (6) , 得

$\tilde{p} = \frac{ρ_{l} v_{l}}{n^{2}} {\tilde{w}}^{‴} + \frac{ρ_{l} v^{'}_{l}}{n^{2}} {\tilde{w}}^{″} + (\frac{i ρ_{l} σ}{n^{2}} - ρ_{l} v_{l}) {\tilde{w}}^{'} (13)$ $\tilde{p} = \frac{ρ_{l} v_{l}}{n^{2}} {\tilde{w}}^{‴} + \frac{ρ_{l} v^{'}_{l}}{n^{2}} {\tilde{w}}^{″} + (\frac{i ρ_{l} σ}{n^{2}} - ρ_{l} v_{l}) {\tilde{w}}^{'} (13)$

式中, 下标l=1或2。l=1代表水层, v₁表示水的粘度;l=2代表泥层, v₂表示泥的粘度。

用式 (13) 代换方程 (7) 中的 $\tilde{p}$ $\tilde{p}$ 得到关于 $\tilde{w}$ $\tilde{w}$ 的四阶偏微分方程:

$\begin{array}{l} {\tilde{w}}^{'}^{‴} + \frac{v^{'}_{l}}{v_{l}} {\tilde{w}}^{‴} + (- 2 n^{2} + \frac{n^{2} (v^{″}_{l} + i σ)}{v_{l}}) {\tilde{w}}^{″} - \\ \frac{2 n^{2} v^{'}_{l}}{v_{l}} {\tilde{w}}^{'} + (n^{4} - \frac{i n^{2} σ}{v_{l}}) \tilde{w} = 0 (14) \end{array}$ $\begin{array}{l} {\tilde{w}}^{'}^{‴} + \frac{v^{'}_{l}}{v_{l}} {\tilde{w}}^{‴} + (- 2 n^{2} + \frac{n^{2} (v^{″}_{l} + i σ)}{v_{l}}) {\tilde{w}}^{″} - \\ \frac{2 n^{2} v^{'}_{l}}{v_{l}} {\tilde{w}}^{'} + (n^{4} - \frac{i n^{2} σ}{v_{l}}) \tilde{w} = 0 (14) \end{array}$

通过以下边界条件, 方程可解 (包括确定未知变量n和a) 。

《图2》

图2水-泥两层模型定义草图

Fig.2 Sketch of the two-layer water-mud model

《4 边界条件及初始条件》

4 边界条件及初始条件

本文垂向二维问题所应考虑的边界条件如下。

水体自由表面:

$\begin{array}{l} - p + 2 ρ_{1} v_{1} \frac{\partial w}{\partial z} + ρ_{1} g η = 0 ； \\ ρ_{1} v_{1} (\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}) = 0 ； \frac{\partial η}{\partial t} = w 。 (15 a - 15 c) \end{array}$ $\begin{array}{l} - p + 2 ρ_{1} v_{1} \frac{\partial w}{\partial z} + ρ_{1} g η = 0 ； \\ ρ_{1} v_{1} (\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}) = 0 ； \frac{\partial η}{\partial t} = w 。 (15 a - 15 c) \end{array}$

水-泥界面:

$\begin{array}{l} u_{1} = u_{2}; w_{1} = w_{2} ； - p_{1} + 2 ρ_{1} v_{1} \frac{\partial w_{1}}{\partial z} + \\ ρ_{1} g η_{2} = - p_{2} + 2 ρ_{2} v_{2} \frac{\partial w_{2}}{\partial z} + ρ_{2} g η_{2} ； \\ ρ_{1} v_{1} (\frac{\partial u_{1}}{\partial z} + \frac{\partial w_{1}}{\partial x}) = ρ_{2} v_{2} (\frac{\partial u_{2}}{\partial z} + \frac{\partial w_{2}}{\partial x}) ； \\ \frac{\partial η_{1}}{\partial t} = w_{2} 。 (16 a - 16 e) \end{array}$ $\begin{array}{l} u_{1} = u_{2}; w_{1} = w_{2} ； - p_{1} + 2 ρ_{1} v_{1} \frac{\partial w_{1}}{\partial z} + \\ ρ_{1} g η_{2} = - p_{2} + 2 ρ_{2} v_{2} \frac{\partial w_{2}}{\partial z} + ρ_{2} g η_{2} ； \\ ρ_{1} v_{1} (\frac{\partial u_{1}}{\partial z} + \frac{\partial w_{1}}{\partial x}) = ρ_{2} v_{2} (\frac{\partial u_{2}}{\partial z} + \frac{\partial w_{2}}{\partial x}) ； \\ \frac{\partial η_{1}}{\partial t} = w_{2} 。 (16 a - 16 e) \end{array}$

泥层底部:

$u_{2} = 0; w_{2} = 0 (17 a - 17 b)$ $u_{2} = 0; w_{2} = 0 (17 a - 17 b)$

初始条件是假设计算从静止状态开始。

《5 泥体输移速度》

5 泥体输移速度

表面水波的振荡荷载作用迫使底床软泥层发生振荡运动响应, 其结果是一个波浪周期平均下来, 泥质点的平均速度是一个正值, 也即导致泥体向波浪传播方向产生净输移 (越接近泥面, 这种净输移量越大) , 称这种现象为泥体输移, 称泥质点的这种平均速度为泥体输移速度。

泥体输移速度的计算表达式为

$U_{L} = U_{E} + \bar{\frac{\partial u_{m}}{\partial x} \int_{0}^{t} u_{m} d t} + \bar{\frac{\partial u_{m}}{\partial z} \int_{0}^{t} w_{m} d t} (18)$ $U_{L} = U_{E} + \bar{\frac{\partial u_{m}}{\partial x} \int_{0}^{t} u_{m} d t} + \bar{\frac{\partial u_{m}}{\partial z} \int_{0}^{t} w_{m} d t} (18)$

式中, 下标L和E含义分别为拉格朗日和欧拉坐标系。方程式右侧第2项与第3项之和称为Stokes's漂移, 并用U_S代表。

Stokes's漂移的物理含义可解释为, 在Euler空间坐标系里, 质点处在椭圆形轨道路径顶部时, 其速度比它处在底部时的速度稍大一些, 虽然不是很大的差别, 但经过一个波浪周期后质点的轨道路径不再封闭——质点沿波浪传播方向发生了净位移。

平均Euler速度U_E的存在是由于软泥的粘度。Sakakiyama和Bijker推导了一个计算U_E的简化方程:

$μ_{m} \frac{\partial^{2} U_{E}}{\partial z^{2}} = \frac{\partial (ρ_{m} \bar{u_{m}^{2}})}{\partial x} + \frac{\partial (ρ_{m} \bar{u_{m} w_{m}})}{\partial z} (19)$ $μ_{m} \frac{\partial^{2} U_{E}}{\partial z^{2}} = \frac{\partial (ρ_{m} \bar{u_{m}^{2}})}{\partial x} + \frac{\partial (ρ_{m} \bar{u_{m} w_{m}})}{\partial z} (19)$

式中系数μ_m类比于牛顿流体粘度, 可取一个波浪周期内的平均值 ^[15,16]。求解方程还需边界条件:

$U_{E} |_{z = - d_{m}} = 0 ； \frac{\partial U_{E}}{\partial z} |_{z = 0} = 0 (20 a - 20 b)$ $U_{E} |_{z = - d_{m}} = 0 ； \frac{\partial U_{E}}{\partial z} |_{z = 0} = 0 (20 a - 20 b)$

因而, 在获得水层和泥层的速度场之后, 泥体输移速度可通过分别计算Stokes's漂移U_S和平均Euler速度U_E, 然后求和获得。

《6 数值模型验证及泥体输移规律》

6 数值模型验证及泥体输移规律

数值计算采用有限差分方法。垂向网格生成是从水-泥界面起向其两侧区域作对数网格划分。数值试验显示, 对数网格的处理办法较好地再了现界面附近区域水平速度梯度较大的边界层流动特征。图3以速度矢量场再现了使用对数网格系统计算的部分流场, 特别是泥-水界面附近的速度变化。图4对比显示了不同泥密度条件下水平速度振幅沿水深的分布。可看出对数网格能够反映水-泥界面附近水平速度大的剪切变化特征。这里坐标z从泥层底部起算, 泥层厚度为0.10 m, 即水-泥界面平均位置在z=0.10 m处。

《图3》

图3对数网格下水层和泥层界面附近速度场

Fig.3 Velocity field in the water and mud layers using logarithm grid

《图4》

图4对数网格下水平速度振幅沿深度分布

Fig.4 Profiles of horizontal velocity amplitude

波浪运动的水体引起底床软泥层振荡并导致泥体输移和表面水波波高衰减。

首先探究表面水波波高衰减率的一个变化规律, 同时验证笔者的模型。泥床厚度与其Stokes's边界层厚度的比值是无量纲参数, 常被用来研究表面水波波高衰减率的变化 ^[5,17]。设水波引起的泥床振荡为层流运动, 于是, 其波动边界层的厚度可直接套用Stokes's边界层厚度表达式:

$δ_{m} \equiv (2 v_{m} / σ)^{1 / 2} (21)$ $δ_{m} \equiv (2 v_{m} / σ)^{1 / 2} (21)$

式中, v_m定义为淤泥粘度的深度平均值, 即:

$v_{m} \equiv \frac{1}{d_{m}} \int_{0}^{d_{m}} v_{e} (z) d z (22)$ $v_{m} \equiv \frac{1}{d_{m}} \int_{0}^{d_{m}} v_{e} (z) d z (22)$

定义无量纲参数 $\tilde{d}$ $\tilde{d}$ 为泥床厚度d_m与其Stokes's边界层厚度δ_m之比, 该无量纲参数在Dalrymple和Liu及Ng中用来研究波浪衰减率的变化 ^[5,17]:

$\tilde{d} = d_{m} / δ_{m} (23)$ $\tilde{d} = d_{m} / δ_{m} (23)$

图5显示了用本文数值模型计算得到的波高衰减率随 $\tilde{d}$ $\tilde{d}$ 变化的关系。容易看出, 对于一定的底泥密度范围 $(1 200 \sim 1 300 k g / m^{3}) ‚ \tilde{d}$ $(1 200 \sim 1 300 k g / m^{3}) ‚ \tilde{d}$ 值在1.2～1.5范围内, 波高衰减最为显著, 即泥层厚度超过其波动边界层厚度20%～50%的情形最有利于表面水波波高衰减, 这与Dalrymple和Liu及Ng的结论很接近 ^[5,17]。

《图5》

图5波高衰减率随无量纲参数 $\tilde{d}$ $\tilde{d}$ 变化

Fig.5 Rate of decay of surface wave height as function of $\tilde{d}$ $\tilde{d}$

其次计算泥体输移速率。淤泥质底床上有水波传播时, 底床泥层在水波作用下发生振荡响应, 导致向波浪传播方向泥体净输移;越接近水-泥界面, 这种净输移量越大。

为进一步检验所建模型用于泥体输移计算, 按Sakakiyama和Bijker经典实验的工况条件进行了计算 ^[15]。图6显示了数值计算结果与实验数据的比较。图中纵坐标为到泥层底面 (刚性底面) 的距离。Lagrange速度廓线, 即泥体总输移速度的计算结果与实验数据比较可以看出, 泥层上半部计算值与实验值符合得很好, 但泥层下半部符合程度稍差一些, 出现小的偏离。可能由于数值计算垂向使用对数网格划分, 对于提高泥-水界面附近流场分辨率很有帮助, 但这也使得泥层接近刚性底面的部分 (由于网格稀疏) 不能得到更好的描述。

与图6比较, 图7的计算条件除水波周期不同外, 其他均相同。比较发现, 软泥床泥体输移速度随波浪周期增大而呈显著增大趋势。图8除泥密度外, 其他计算条件与图6相同, 用来比较软泥床泥体输移速度随泥床密度变化的趋势。泥体输移速度随泥密度减小略有增大。比较也发现, 软泥床泥体输移速度随表面水波波高的增大而快速增大。张庆河从实验中也得到类似结论 ^[18]。

《图6》

图6模型计算结果与实验测量值比较

Fig.6 Measured and calculated profiles of mud mass transport velocity

《图7》

图7泥体输移速度计算结果

Fig.7 Calculated profiles of mud mass transport velocity

《图8》

图8泥体输移速度计算结果

Fig.8 Calculated profiles of mud mass transport velocity

《7 结论》

7 结论

所建立的表面水波与淤泥质底泥层相互作用耦合数值模型有以下特色:a.应用了半经验的淤泥流变关系模型, 该流变模型表达淤泥体在水波振荡剪切外力作用下其应变 (及应变率) 与应力关系的显著滞后回线特征和非线性特征。b.在数值处理技巧上, 为提高数值格式对水-泥界面附近流动的分辨力, 较好地再现界面附近速度梯度大的流动特点, 垂向网格方法是从水-泥界面向两侧区域使用对数网格。主要结论如下:

1) 在水波振荡外力作用下高岭土软淤泥体流变性质表现为淤泥剪切应变及其应变率对剪切应力的响应呈显著的滞后回线特征, 即当剪切应力减小为0时, 泥体中的剪切应变及其应变率并不为0, 而是分别有一个剩余量, 剩余剪切应变及剩余剪切应变率暂留在泥体中。淤泥流变学的这种滞后响应特征也显示了淤泥流变学的强非线性特征。

2) 对数网格方法对提高数值格式对水-泥界面附近流动的分辨力, 较好地再现界面附近速度梯度大的流动特点是有帮助的。

3) 应用所建立的耦合数值模型计算了水波在淤泥质底床上传播时发生沿程衰减变化及水波引起底床泥体输移响应, 这些数值结果和相关的实验室测量值符合较好, 模型得到了初步验证。

然而, 水波与淤泥质底床相互作用过程是十分复杂的水动力现象, 其中包括许多挑战性的工作。例如, 粘性泥沙的絮凝、沉降和沉积等过程, 泥床的液化以及考虑水流与波浪共存等。从这些方面看, 今后还需更多的基础性工作进一步丰富数值模型的内涵。

展示更多