瓯江河口潮流水动力场数值模拟

Numerical Simulation of Tide Hydrodynamic Field in Oujiang Estuary

摘要: 针对瓯江口的复杂地形和强潮条件，建立了正交曲线坐标系下的3D、非恒定瓯江口潮流水动力学模型，研究了瓯江口在径流与潮汐共同作用下的水动力学特点。模型区域梅岱至洞头长约80 Km，计算网格8693个。模型应用2010年大、中和小潮条件，对8个站点潮位及6个站点的流场进行了验证，模拟了不同潮汐和径流条件下的潮位及水动力场变化。模型值与实测值吻合较好，模型有较好适应性，能满足波流作用下的瓯江口水动力研究，并可用于模拟建筑物对水动力场的影响。

Abstract: Considering the complex topography and strong tide conditions of Oujiang Estuary, the 3D, time- dependent hydrodynamic model under orthogonal curvilinear coordinate system was established, and the characteristics of the interaction of runoff and tide in Oujiang Estuary was studied. The length of the model area from Mei’ao to Dongtou is about 80 km, with 8693 computational grids. The model was applied to verify the flow field of 8 stations and 6 stations under the conditions of spring, middle and neap tides in 2010, and simulated the changes of tidal level and hydrodynamic field under different tidal and runoff conditions. The model values are in good agreement with the measured values, and the model has good adaptability. It can be used to simulate the influence of buildings on hydrodynamic field.

1. 研究背景

瓯江口位于浙江东南部、北纬27˚~28˚、东经119˚~121˚之间，是浙江第二大河流瓯江的入海口，也是重要的沿海经济带(如图1)。近年来，随着地区经济发展和航运业发展，研究瓯江口水动力条件和波流特征对开展航运基础设施建设具有十分重要的意义，并取得了一些研究成果 [1] [2] [3]。影响瓯江河口水动力场的因素很多，主要包括：1) 河口区域的复杂地形。口门往上形成了灵昆岛、江心屿、七都岛等江心洲，河道分汊呈藕节状。2) 潮汐条件和径流变化。受瓯江干流、楠溪江径流及潮流耦合叠加影响。3) 盐度变化和风应力因素。受台风及咸水入侵作用，对泥沙絮凝淤积有一定影响。4) 人工建筑物对水动力场的影响，沿江分布有众多码头、船厂、丁顺坝、桥梁等水工建筑物 [4]。本文针对瓯江口复杂地形和强潮条件，建立了瓯江口3D、非恒定水动力模型，模型基于正交曲线网格，研究瓯江口在径流与潮汐共同作用下的水动力学特点，并研究了不同潮汐和径流条件下水动力场分布。

Figure 1. Sketch of Oujiang estuary (a) and site distribution (b)

图1. 瓯江口地形图(a)和站点分布图(b)

2. 数学模型的建立

2.1. 模型方程

垂向σ坐标变化有：

$z=\left(z^{*}+h\right)/\left(\zeta +h\right)$ (1)

这里：*代表原始物理坐标，h和 $\zeta$ 是河底高程和水位。

在静水压、Boussinesq、垂向正压条件下，基本方程由不可压、变密度流的动量方程、连续方程和温度、盐度输运方程组成 [5] [6]：

${\partial }_t\left(m_x{m}_{y}H\right)+{\partial }_x\left(m_{y}Hu\right)+{\partial }_y\left(m_{x}Hv\right)+{\partial }_z\left(m_x{m}_{y}w\right)=Q_H$ (2)

$\begin{array}{l}{\partial }_t\left(m_x{m}_{y}Hv\right)+{\partial }_x\left(m_{y}Huv\right)+{\partial }_y\left(m_{x}Hvv\right)+{\partial }_z\left(m_x{m}_{y}wv\right)+m_x{m}_y{f}_{e}Hu\\ =-m_{x}H{\partial }_y\left(p+g\eta \right)-m_x\left({\partial }_{y}h-z{\partial }_{y}H\right){\partial }_{z}p+{\partial }_z\left(m_x{m}_y{H}^{-1}{A}_v{\partial }_{z}u\right)+Q_v\end{array}$ (3)

$\begin{array}{l}{\partial }_t\left(m_x{m}_{y}Hu\right)+{\partial }_x\left(m_{y}Huu\right)+{\partial }_y\left(m_{x}Hvu\right)+{\partial }_z\left(m_x{m}_{y}wu\right)-m_x{m}_y{f}_{e}Hv\\ =-m_{y}H{\partial }_x\left(p+g\eta \right)-m_y\left({\partial }_{x}h-z{\partial }_{x}H\right){\partial }_{z}p+{\partial }_z\left(m_x{m}_y{H}^{-1}{A}_v{\partial }_{z}u\right)+Q_u\end{array}$ (4)

${\partial }_{z}p=-gHb=-gH\left(\rho -{\rho }_o\right){\rho }_o^{-1}$ (5)

${\partial }_t\left(m\xi \right)+{\partial }_x\left(m_{y}H{\displaystyle \underset{0}{\overset{1}{\int }}u\text{d}z}\right)+{\partial }_y\left(m_{x}H{\displaystyle \underset{0}{\overset{1}{\int }}v\text{d}z}\right)=0$ (6)

$\rho =\rho \left(S,T\right)$ (7)

${\partial }_t\left(mHS\right)+{\partial }_x\left(m_{y}HuS\right)+{\partial }_y\left(m_{x}HvS\right)+{\partial }_z\left(mwS\right)={\partial }_z\left(m{H}^{-1}{A}_b{\partial }_{z}S\right)+Q_S$ (8)

${\partial }_t\left(mHT\right)+{\partial }_x\left(m_{y}HuT\right)+{\partial }_y\left(m_{x}HvT\right)+{\partial }_z\left(mwT\right)={\partial }_z\left(m{H}^{-1}{A}_b{\partial }_{z}T\right)+Q_T$ (9)

式中，u，v，w是三个方向流速， $H=\zeta +h$ 是水深， $m_x$， $m_y$ 是坐标变化张量， $m=m_x{m}_y$，f是Coriolis系数， $m_x{m}_y{f}_e\equiv m_x{m}_{y}f-u{\partial }_y{m}_x+v{\partial }_x{m}_y$，压力 $p={\rho }_{0}gH\left(1-z\right)$， ${\rho }_0$ 是参考密度；密度 $\rho$ 是温度T和盐度S的函数，在不可压、滞弹性假定下(Mellor，1991, Clark和Hall，1991)，也可以是压力的弱函数；浮力b定义为垂向密度与参考密度 ${\rho }_0$ 的相对偏离。A_v，A_b是垂向紊动粘度和扩散率，联立求解(2)~(8)式可得变量u，v，w，P， $\zeta$， $\rho$，S，T。

垂向紊动动量扩散率(A_v)和质量传输系数(A_b)，可采用M-Y (1982)方程求解。方程定义湍流强度 $q^2$ 、湍流尺度l和Richardson数 $R_q$。

$A_v={\varphi }_{v}ql=0.4{\left(1+36R_q\right)}^{-1}{\left(1+6R_q\right)}^{-1}\left(1+8R_q\right)ql$ (10)

$A_b={\varphi }_{b}ql=0.5{\left(1+36R_q\right)}^{-1}ql$ (11)

$R_q=\frac{gH{\partial }_{z}b}{q^2}\frac{l}{H^2}$ (12)

${\varphi }_b$， ${\varphi }_v$ 称为稳定函数，用于在稳定和非稳定密度分层下，减小或增强垂向混合和输运。

正交曲线坐标系下，垂向无量纲化，采用σ坐标：

$\begin{array}{l}{\partial }_t\left(mH{q}^2\right)+{\partial }_x\left(m_{y}Hu{q}^2\right)+{\partial }_y\left(m_{x}Hv{q}^2\right)+{\partial }_z\left(mw{q}^2\right)\\ ={\partial }_z(m{H}^{-1}{A}_q{\partial }_z{q}^2)-2mH{\left(B_{1}l\right)}^{-1}{q}^3+2m\left(H^{-1}{A}_v\left({\left({\partial }_{z}u\right)}^2+{\left({\partial }_{z}v\right)}^2\right)+{\eta }_p{c}_p{D}_p{\left(u^2+v^2\right)}^{3/2}+g{K}_v{\partial }_{z}b\right)+Q_q\end{array}$ (13)

$\begin{array}{l}{\partial }_t\left(mH{q}^{2}l\right)+{\partial }_x\left(m_{y}Hu{q}^{2}l\right)+{\partial }_y\left(m_{x}Hv{q}^{2}l\right)+{\partial }_z\left(mw{q}^{2}l\right)\\ ={\partial }_z\left(m{H}^{-1}{A}_q{\partial }_z{q}^{2}l\right)-mH{B}_1^{-1}{q}^3\left(1+E_2{\left(\frac{l}{\kappa Hz}\right)}^2+E_3{\left(\frac{I}{\kappa H\left(1-Z\right)}\right)}^2\right)\\ +2m{H}^{-1}{E}_{1}l{A}_v\left({\left({\partial }_{z}u\right)}^2+{\left({\partial }_{z}v\right)}^2+{\eta }_p{c}_p{D}_p{\left(u^2+v^2\right)}^{3/2}+g{K}_v{\partial }_{z}b\right)+Q_l\end{array}$ (14)

式中，经验常数 $\left(E_1,E_2,E_3\right)=\left(1.8,1.33,0.25\right)$。

2.2. 边界条件

平面边界条件：进出边界采用开边界，河岸及固定建筑物采用固壁边界。

风剪应力有下式：

$\left({\tau }_{xz},{\tau }_{yz}\right)=\left({\tau }_{sx},{\tau }_{sy}\right)=c_s\sqrt{U_W^2+V_W^2}\left(U_W,V_W\right)$ (15)

剪切系数

$c_s=0.001\frac{{\rho }_a}{{\rho }_w}\left(0.8+0.065\sqrt{U_W^2+V_W^2}\right)$ (16)

U_W和V_W是水面风速，ρ_a，ρ_w是大气和水的密度。

河床底部的剪切力：

$\left({\tau }_{xz},{\tau }_{yz}\right)=\left({\tau }_{bx},{\tau }_{by}\right)=c_b\sqrt{u_1^2+v_1^2}\left(u_1,v_1\right)$ (17)

下标1表示底层水体流速，假定近床流场剖面满足对数分布，则

$c_b={\left(\frac{\kappa }{\ln \left({\Delta }_1/2z_o\right)}\right)}^2$ (18)

式中κ是von Karman常数，Δ₁是底层厚度， $z_o=z_o^{*}/H$ 是无量纲糙度高度。

2.3. 计算网格

模型区域如图1，上起瓯江梅岱，下至洞头岛，全长约80 km。模型计算网格采用正交贴体网格，i×j为233 × 97，共采用了8693个格点，如图2。

Figure 2. The computational grid (I × J = 233 × 97)

图2. 计算网格(I × J = 233 × 97)

3. 模型验证

模型验证采用2010年大、中、小潮条件下，以温州、状元、龙湾、乌牛、海思、黄华、洞头和坎门等8个站点潮位与1#~6#共6个站点的流场观测值与模拟值对比得出。站点位置见图1(b)。

潮位校准采用了温州、状元、龙湾、乌牛、海思、黄华、洞头和坎门等8个站点潮位进行比较。如下图3~5，圆点代表观测值，实线代表模拟值。大潮条件下，潮位观测值与模拟值之差平均值最大为0.1 m，最小为0.01 m，均方根误差RSME最大为0.15 m，最小为0.03 m；中潮条件下，潮位观测值与模拟值之差平均值最大为0.08 m，最小为0.01 m，均方根误差RSME最大为0.11 m，最小为0.04 m；小潮条件下，潮位观测值与模拟值之差平均值最大为0.1 m，最小为0.01 m，均方根误差RSME最大为0.15 m，最小为0.03 m。

流速观测值与模拟值之差平均值最大为0.12 m/s，最小为0.03 m/s，计算结果与观察值吻合良好。

Figure 3. Comparisons of observed and simulated values at different stations under high-tide

图3. 大潮下的各站点观测值和模拟值比较图

Figure 4. Comparisons of observed and simulated values at different stations under mid-tide

图4. 中潮下的各站点观测值和模拟值比较图

Figure 5. Comparisons of observed and simulated values at different stations under low-tide

图5. 小潮下的各站点观测值和模拟值比较图

4. 模型应用

4.1. 潮差

对模型结果进行分析，可以得出瓯江口河段在不同径流和潮汐状态下的潮差，如下表1，潮位和潮差是研究河口水动力的基础，对于认识其河口水动力状况具有重要意义。

Table 1. Statistics on tidal range of each station (unit: m)

表1. 各站点潮差统计表(单位：m)

4.2. 潮流场特征

大潮涨急、落急工况下流速场分布见图6、图7。根据计算结果，工程海域流场特征表现为，瓯江河段潮流受河道深槽地形控制，基本上作往复运动。涨潮流向为W~NW向，落潮流向为SE~E向。瓯江河口处在为强潮海区，涨落潮流速较大，且总体呈现落潮大于涨潮特点。七都岛段水域，大潮情况下涨潮平均流速约介于0.7~1.0 m/s，涨潮最大流速约介于1.1~1.6 m/s，落潮平均流速约介于0.9~1.3 m/s，落潮最大流速约介于1.2~1.9 m/s。七都岛北汊水域在大潮期间其平均流速约介于0.8~1.2 m/s之间，最大流速约介于1.1~1.7 m/s之间。七都南汊水域在大潮期间其平均流速约介于0.7~1.0 m/s之间，最大流速约介于1.0~1.5 m/s之间。从平面分布来看，七都岛北汊水域流速强度较大于南汊水域流速，北汊水流流速更强劲些。从垂向分布来看，表层流速最大，0.6 H层流速其次，底层流速最小。表层流速与底层流速的倍数关系多介于1.5~1.8倍。

Figure 6. Distribution of velocity field in spring tide

图6. 大潮涨急流速场分布

Figure 7. Distribution of Velocity field in spring ebb

图7. 大潮落急流速场分布

5. 结论

瓯江口属于强潮河口，并具有复杂地形条件，为研究该河口的水动力特点，生成了正交曲线坐标网格，并建立了3D、非恒定瓯江口水动力学模型，研究了瓯江口在径流与潮汐共同作用下的水动力学特点。模拟了梅岱至洞头河段，采用2010年的大、中和小潮条件对8个站点潮位进行了验证，模拟值与实测值吻合较好，表明模型有较好适应性，能满足波流作用下的瓯江口水动力研究。应用于该河口段的水动力场的流场和潮差研究，结果表明：从上游至下游潮差逐渐减小，上游江心寺大潮潮差6.12 m，小潮潮差4.85 m；坎门大潮潮差4.93 m，小潮潮差3.74 m。瓯江河口涨落潮流速较大，且总体呈现落潮大于涨潮特点。涨潮平均流速约0.7~1.0 m/s，落潮平均流速约介于0.9~1.3 m/s。七都岛北汊水域流速强度较大于南汊水域流速。

基金项目

国家重点研发计划(2016YFC0402006)。