长江口洪水期潮波变形数值模拟研究

路川藤，罗小峰，韩玉芳

(南京水利科学研究院 水文水资源与水利工程科学国家重点实验室，江苏 南京 210029)

研究潮波变形的方法主要有三种。一为现场水文观测研究，这是研究河口潮波变形的最直接手段，结论可信实用，但成本较高。如Kosuth 等［2］指出特定径流与潮汐情况下，亚马逊河口可同时出现3 个波峰。二为半经验解析模型，是估算潮波传播比较实用的工具，但计算结果误差较大。该方法基于的理论方程为圣维南方程以及潮波能量输运方程，通过不同的离散格式，离散理论方程，在概化地形的基础上，考虑底摩擦［3］、河口边界形状［4-6］、径流对潮差、潮波能量沿程变化的影响。三为数学模型研究，是目前研究潮波传播的主要手段。数学模型采用不同方法［7-8］离散N-S 方程，在计算机上实现程序运算，计算成本低，可操作性强，可视化好，但数学模型计算结果的优劣取决于计算格式、网格的疏密［9］及计算参数［10］等众多因素。

在近海或者河口地区，水深相对较浅，潮波向前传播过程中，高潮位潮波传播速度大于低潮位，即潮波传播过程中，潮波变形会逐渐加剧。然据长江口2005 ～2010年实测水文资料(图1 所示)分析知，洪水期，潮波上溯过程中，潮波变形具有先加剧后趋缓的特点，枯水期这一现象不明显，这有别于一般的河口潮波传播认识。为深入认识洪水期长江口潮波变形特点，建立大通至外海的长江口整体数学模型，分析长江口洪水期潮波变形特征。

图1 长江口潮波变形沿程变化Fig.1 The tidal wave deformation in the Yangtze estuary

1 控制方程

1.1 控制方程及离散

水流运动方程向量形式:

1.2 方程求解

采用三角形网格对计算区域进行离散，并将单一的网格单元作为控制单元，水深布置在网格顶点，其他物理变量配置在每个单元的中心，如图2 所示。

将第i 号控制元记为Ωi，在Ωi上对向量式的基本方程组(1)进行积分，并利用Green 公式将面积分化为线积分，得

方程(2)求解主要分三部分，一为对流项求解;二为紊动项求解;三为底坡项处理。对流项数值通量可采用Roe 格式的近似Riemann 解，紊动项采用单元交界面的平均值计算通过该界面紊动粘性项的数值通量，有限体积法底坡项若不加任何处理，则会造成静水的伪流动现象，如图3 所示。

图2 物理量布置示意Fig.2 The arrangement of physical quantities

图3 有限体积法伪流动现象Fig.3 The phenomenon of false flow in FVM

为避免水流伪流动现象，记某个控制体三个

顶点坐标为xi，yi，zi(i =1，2，3)，三个顶点按逆时针方向排序，由这三个顶点确定的平面方程为:

当Dz≠0 时，式(4)可写为:

1.3 边界条件

1.3.1 水位开边界

开边界又分为急流开边界和缓流开边界，因这里所建模型为缓流模型，故只给出缓流开边界的处理方法。Roe 利用间断左右的物理量UR、UL构造常数矩阵代替非线性雅可比矩阵A(Un)，并提出，当UR，UL→Un，即:

式中:cL和cR表示单元左右静水波传播速度。

则:

式中:Zd为边界上通量积分点处的底高程。

1.3.2 闭边界

采用镜像法［11］处理。在闭边界外侧虚拟一个单元，边界上两侧的法向流速相反，切向流速相同，即DR= DL，un，R= － un，L，uτ，R= － uτ，L，un、uτ表示单元法向和切向流速。

1.3.3 动边界

动边界是数学模型的难点，动边界的处理精确与否直接影响模型计算的质量。文中动边界技术采用干湿法，当单元水深小于hmin(动边界水深设定值)时，认为该单元为干单元，不参与数值计算，当水深大于hmin时，该单元重新参与数值计算。

2 数学模型

2.1 模型范围

长江潮区界位于安徽大通，大通以上水域水位基本不受潮波影响，作为模型的上边界;长江口外－50 m等深线处受径流影响可忽略不计，作为模型外边界，模型总长700 多公里。模型北至江苏吕四港南侧，南至金山嘴，如图4 所示。模型网格总数114 489 个，最小边长128 m，时间步长4 s，紊动粘性系数取常数0.1，动边界水深0.02 m。

图4 数学模型范围Fig.4 The range of mathematical model

2.2 模型验证

徐六泾以上地形选用2003年长江下游地形测量数据，徐六泾至绿华山地形数据选用2005年8月长江口实测地形数据，外海地形选用大范围海图地形数据。模型计算时间为2005年8月18日～2005年8月25日，上游边界采用大通实时流量控制，外海潮位边界由东中国海模型提供，验证点位置如图4 所示。由表1及图5 ～6 知，各站高、低潮位值偏差大都均在10 cm 之内，个别点高低潮位值偏差超过10 cm，高低潮位相位偏差在30 分钟内。CS6 点位于航道边缘，流速偏差较大，因网格较粗，模型地形与实际地形有偏差，其他各点偏差基本都在10%之内，验证良好，满足规程要求。

表2 为各汊道分流比情况验证，长江口北支日益萎缩，至2001年时，其落潮分流比仅占1% ～4%，南支分流比96% ～99%［12］，南北港及南北槽分流比选用2006年8月大潮实测资料验证，本数学模型计算分流比值与实测值基本吻合。验证说明本数学模型具有复演长江口潮波传播的能力。

表1 潮位误差分析Tab.1 Tidal level error analysis

表2 落潮分流比验证Tab.2 Ebb flow ratio verification

图5 模型潮位验证Fig.5 The tidal level verification

图6 模型流速验证Fig.6 The tidal current verification

3 潮波变形影响因素分析

影响河口潮波变形的因素主要有河床底摩擦、河口边界形状、径流、水深及边滩反射等。河床底摩擦及径流削弱潮波能量，而河口边界形状收缩与边滩反射使潮波能量增强。对于某一特定河口，边界形状、浅滩、水深等因素短时间内相对稳定，不会发生大的格局变化，其对潮波传播的影响亦相对稳定。

河床底摩擦在数学模型上体现为糙率，糙率是表征河床底部和岸壁影响水流阻力的综合因素系数，与河床底质构成及水深密切相关。径流对潮波传播的影响反映在径流量的大小对潮波传播阻力影响不同。因此，年际内对潮波变形影响较大的因素为河床底摩擦与径流。下文基于数学模型，研究河床底摩擦与径流对潮波变形的影响。

3.1 计算条件

3.1.1 糙率

本数学模型空间尺度大，自大通至外海总长约700 km，上下游糙率取值差别大。表3 列举了前人在长江口数学模型中糙率的取值，由表知，大通附近最大糙率取值0.032，最小0.018，口外海滨段取值约为0.012。

表3 长江感潮河段糙率选取Tab.3 Roughness selection in the Yangtze estuary

文中数学模型糙率采用附加糙率方法，即:n = n0+n＇/h，基本糙率n0自大通至外海线性插值，附加糙率n＇计算过程中固定不变，基本糙率分四种工况计算讨论，见表4 所示。

表4 数学模型中大通至外海4 种基本糙率n0的取值Tab.4 The selection of basic roughness n0 from Datong to seas in mathematic model

3.1.2 径流

模型上边界位于长江口潮区界安徽大通，大通流量即为径流量。据大通站实测资料分析［21］，大通站多年平均流量为28 518 m3/s，多年最小月平均流量10 400 m3/s，多年最大月平均流量48 600 m3/s，最大月平均洪峰流量84 000 m3/s。为充分分析径流对潮波传播的影响，分别取10 000、30 000、60 000和90 000 m3/s四种流量代表枯水流量、年平均流量、洪水流量和特大洪水流量。

3.2 潮波变形对糙率的响应

谢才公式:

式中:v 为流速，R 为水力半径，J 为水面比降。

由上式可知，糙率的变化直接影响流速的大小。潮汐河口中，潮量由流速和水位共同反应，当潮量不变时，流速的变化必然影响潮位。由图7 知，随着糙率的增大，河口平均水位逐渐升高，越往上游，水位升高幅度越大，南京站最大糙率与最小糙率水位差接近0.6 m，换言之，糙率的增大可增大河口平均水面比降(如表5 所示)。徐六泾以下水域受潮汐影响大，且水面开阔，糙率对水位影响相对较小。

图8 为潮波变形对糙率变化的响应。徐六泾站以下水域，糙率的变化对潮波变形影响相对较小。徐六泾以上水域，随着糙率的增大，涨落潮历时比值逐渐增大，潮波变形加剧。自下游至上游，落潮历时先增大后减小的趋势不变，说明糙率不是洪季潮波变形拐点形成的主要原因。

图7 平均水位水面线变化Fig.7 The variation of mean water level

表5 糙率变化条件下长江感潮河段平均水位水面比降Tab.5 Mean water surface slope of different roughnesses in the Yangtze tidal reach (×10 －5)

图8 潮波变形对糙率变化的响应Fig.8 Tidal wave deformation caused by roughness changing

3.3 潮波变形对径流的响应

径流量在水位上的反映为水位随径流量的增大而升高，越往上游该现象越明显(如图9 所示)，因越往上游径流作用越大，潮汐作用越弱，南京站径流为10 000 m3/s 与90 000 m3/s 时，平均水位差超过7 m，北槽口外，潮汐作用为主，径流作用微弱，故径流变化后，该处平均水位变化微小。由此可见，径流量的增大改变了长江河口平均水位的水面比降(表6)，越往上游，水面坡降增加越大。

图9 平均水位水面线变化Fig.9 The variation of mean water level

图10 为潮波变形对径流变化的响应。当径流量为10 000 m3/s 时，自下游至上游，落潮历时呈逐渐增大的趋势，由于径流量较小，径流对潮波传播阻力小，落潮历时整体相对较短。径流量增大到30 000 m3/s 时，落潮历时整体大于径流量为10 000 m3/s 时，且自下游至上游，落潮历时出现先增大后减小的现象，拐点位于扬中附近。径流量增大到60 000 m3/s 时，自下游至上游，落潮历时出现先增大后减小的现象明显，且转折点向下游偏移至江阴上游，当径流量增大到90 000 m3/s 时，拐点位于江阴附近。据路川藤［20］的研究，径流量为30 000、60 000、90 000 m3/s 时，潮流界的位置位于炮子洲、天生港、白茆河处，潮波变形拐点均位于潮流界上游。潮流界下游，随着径流量的增大，潮波上溯阻力增大，落潮历时延长。潮流界上游，随着径流量的增大，落潮历时出现先增大后减小的现象愈趋显著，拐点随径流量的增大向下游偏移。

表6 径流变化条件下平均水位水面比降Tab.6 Mean water surface slope under different runoffs (×10 －5)

图10 潮波变形对径流量变化的响应Fig.10 The tidal wave deformation caused by changing runoff

3.4 潮波变形分析

前文研究了特定河道中，径流及糙率对潮波变形的影响。经分析知，糙率对潮波变形的影响远小于径流，径流量的大小在河道中主要反映为平均水位水面比降(图9、表6)的变化，即平均水位水面比降的变化对潮波变形影响较大。

洪水期，潮波上溯过程中，潮波变形存在明显的转折点。对于转折点以下河段，径流量的增大致使潮波上溯阻力增大，潮波能量损耗多，潮波变形加剧。转折点上游河段，潮波变形逐渐趋缓的现象难以用潮波理论解释，下文试从力学角度解释该现象。

假设径流与潮波二者相互独立，潮波叠加在径流水面之上，潮波传播在河道上反应为水面的升降。将平均水面线以上高潮位水体看作一个整体，高潮位水体在径流水面向上游运动，其受自身重力和径流向下游运动而产生的水流拖曳力及水面支持力，如图11 所示。重力和水流拖曳力是高潮位向上游运动的阻力，可表示为:

式中:FH为高潮位水体向上游运动的阻力，f 为水流拖曳力，G 为重力，α 为平均水面与水平方向的夹角。

低潮位向上游传播表现为水位的下降，若将平均水面线以下低潮位看作一个整体，则低潮位向上游运动只受径流向下游运动而产生的水流拖曳力(设涨落潮过程中水流拖曳力不变):

式中:FL为低潮位水体向上游运动的阻力。

当α →0 时，FH→FL，说明平均水面比降越小，高潮位重力对潮波向上游传播影响越小，相反，平均水面比降越大，高潮位重力对潮波向上游传播影响越大。

洪水期，潮流界以上河段水面比降明显大于潮流界以下河段(表6)，潮波变形程度受高低潮位潮波传播速度差异和高潮位自身重力两个因素影响。高、低潮位潮波传播速度差异致使越往上游潮波变形越剧烈，高潮位自身重力致使越往上游潮波变形越趋缓;当高潮位自身重力对潮波传播的影响大于高低潮位潮波传播速度差异时，潮波变形趋缓;当高潮位自身重力对潮波传播的影响小于高低潮位潮波传播速度差异时，潮波变形继续加剧。

图11 潮波传播力学分析示意Fig.11 The force analysis of tidal wave propagation

4 结 语

1)据长江口实测水文资料分析，洪水期，潮波上溯过程中潮波变形先加剧后趋缓，枯水期这一现象不明显。

2)基于非结构网格FVM 方法，研究了潮波变形对糙率和径流量变化的响应。研究认为径流量对潮波变形的影响远大于糙率，径流量对潮波变形的影响主要原因在于平均水位水面比降的变化。

3)据潮波上溯力学分析，认为潮流界以上水域，高潮位上溯传播的阻力较大，当阻力对潮波传播的影响大于潮差引起的高低潮位潮波传播速度差异时，便造成了洪水期潮波变形先加剧后趋缓的特点。

［1］黄胜，卢启苗.河口动力学［M］.北京:水利电力出版社，1992.(HUANG Sheng，LU Qimiao.Estuarine dynamics［M］.Beijing:Water Power Press，1992.(in Chinese))

［2］KOSUTH P，CALLEDE J，LARAQUE A.Ocean tide waves progagation along downstream amazon river:measuring the amazon discharge at the estuary［J］.Building Partnerships，2000(310):1-10.

［3］VAN RIJN C.Analytical and numerical analysis of tides and salinities in estuaries;part I:tidal wave propagation in convergent estuaries ocean dynamics［J］.Ocean Dynamics，2011，61(11):1719-1741.

［4］HUBERT H，SABENIJE G.Lagrangian solution of st.Venant’s equations for alluvial estuary［J］.Journal of Hydraulic Engineering，1992，118(8):1153-1163.

［5］HUBERT H，SABENIJE G.Determination of estuary parameters on basis of lagrangian analysis［J］.Journal of Hydraulic Engineering，1993，119(5):628-642.

［6］SABENIJE G，HUBERT H.Analytical expression for tidal damping in alluvial estuaries［J］.Journal of Hydraulic Engineering，1998，124(6):615-618.

［7］PAN Cunhong，LIN Bingyao，MAO Xianzhong.Case study:numerical modeling of the tidal bore on the Qiantang river，China［J］.Journal of Hydraulic Engineering，2007，133(2):130-138.

［8］BAO W，ZHANG X，YU Z，et al.Real-time equivalent conversion correction on river stage forecasting with manning’s formula［J］.Journal of Hydrologic Engineering，2011(1)，16:1-9.

［9］HASAN G，DIRK VAN M，CHEONG H，et al.Improving hydrodynamic modeling of an estuary in a mixed tidal regime by grid refining and aligning［J］.Ocean Dynamics，2011(11):1-15.

［10］路川藤.长江口潮波传播［D］.南京:南京水利科学研究院，2009.(LU Chuanteng.The tidal propagation in Yangtze estuary［D］.Nanjing:Nanjing Hydraulic Research Institute，2009.(in Chinese))

［11］陈丕翔.基于有限体积法的二维水流水质模拟［D］.南京:河海大学，2007.(CHEN Pixiang.Water quanlity simulation based on 2D finite volume method［D］.Nanjing:Hohai University，2007.(in Chinese))

［12］恽才兴.图说长江河口演变［M］.北京:海洋出版社，2010.(YUN Caixing.The Yangtze river estuary evolution［M］.Beijing:Ocean Press，2010.(in Chinese))

［13］肖庆华，岳志远，刘怀汉，等。长江下游二维浅水非恒定流数值模拟［J］.水力发电学报，2013，32(5):115-121.(XIAO Qinghua，YUE Zhiyuan，LIU Huanhan，et al.Numerical simulation of unsteady 2D shallow water flow of the lower Yangtze river［J］.Journal of Hydroelectric Engineering，2013，32(5):115-121.(in Chinese))

［14］夏云峰.一种简便的非交错曲线网格三维水流数值模型［J］.水动力学研究与进展，2002，17(5):638-646.(XIA Yunfeng.A simplified 3D flow model for natural river with curvilinear collocated grids［J］.Journal of Hydrodynamics，2002，17(5):638-646.(in Chinese))

［15］梁婷.长江河口洲滩高强度开发水动力响应研究［D］.南京:南京水利科学研究院，2011.(LIANG Ting.The study of hydrodynamic response on Yangtze estuary segment beach high-intensity development［D］.Nanjing:Nanjing Hydraulic Research Institute，2011.(in Chinese))

［16］李健镛.长江大通—徐六泾河段水沙特征及河床演变研究［D］.南京:河海大学，2007.(LI Jianyong.Study on the evolution and sediment characteristics in the Datong-Xuliujing River［D］.Nanjing:Hohai University，2007.(in Chinese))

［17］曾小辉，李国杰，姜昱.长江感潮河段平面二维潮流数值模拟［J］.水运工程，2012(4):12-16.(ZENG Xiaohui，LI Guojie，JIANG Yu.Numerical solution of 2-D tidal flow of the estuary in the Yangtze river［J］.Port ＆Waterway Engineering，2012(4):12-16.(in Chinese))

［18］窦希萍，李褆来，窦国仁.长江口全沙数学模型研究［J］.水利水运科学研究，1999(2):136-145.(DOU Xiping，LI Tilai，DOU Guoren.Numerical model study on total sediment transport in the Yangtze River estuary［J］.Hydro-Science and Engineering，1999(2):136-145.(in Chinese))

［19］马进荣，陈志昌.长江口风暴潮流场计算［J］.水利水运工程学报，2003(1):35-39.(MA Jinrong，CHEN Zhichang.Simulation of storm surge current in the Yangtz Estuary［J］.Hydro-Science and Engineering，2003(1):35-39.(in Chinese))

［20］路川藤.长江口潮波传播模拟研究及主要影响因素分析［D］.南京:南京水利科学研究院，2013.(LU Chuanteng.The simulation of tidal propagation in Yangtze estuary and the analysis of its main influencing factors［D］.Nanjing:Nanjing Hydraulic Research Institute，2013.(in Chinese))

［21］沈焕庭，李九发.长江河口水沙输运［M］.北京:海洋出版社，2011.(SHEN Huanting，LI Jiufa.Water and sediment transport in the Yangtze estuary［M］.Beijing:Ocean Press，2011.(in Chinese ))