孙孟喆，胡 鹏, 2，李有为，刘奇峰, 曹志先
(1. 浙江大学 海洋学院，浙江 舟山 316021; 2. 浙江大学 舟山海洋研究中心，浙江 舟山 316021; 3.长江航道规划设计研究院，湖北 武汉 430040; 4. 武汉大学 水利水电学院，湖北 武汉 430072)
福姜沙水道位于长江口感潮河段，上起江阴鹅鼻嘴，下迄护槽港，全长约30 km。水动力受上游径流和下游潮流的共同作用，泥沙输移和河床演变复杂，是航道整治的难点河段[1-4]。福姜沙水道在平面上呈现“两级分汊、三汊并存”的形态：主流经福姜沙分为左汊和右汊(福南水道)，左汊水流经双涧沙分为福北水道和福中水道，形成长江下游典型的多级分汊多级汇流的复杂河势格局。福北水道和福中水道为主航道，分别为12.5 m深水航道的上行通航分道和上下行通道。历史上双涧沙的冲淤变化以及滩面横流、窜沟发展对福中、福北水道的航行条件产生了不利影响[5-6]，为了限制双涧沙的不稳定对各汊道航行条件的不利影响，2009年至2011年期间实施了双涧沙守护工程，分为双涧沙头部潜堤、北顺堤、南顺堤3部分。工程虽然达到“封堵窜沟、固滩护沙”的整治目的，但由于工程规模有限，滩面流未完全遏制，两侧深槽发展仍不能有效控制[7-9]，甚至对航道产生淤积等不良影响[4]。为解决这一问题，将长江口深水航道进一步向上延伸，航道部门对福姜沙水道实施了深水航道二期工程，在福姜沙左缘建设4道丁坝，在双涧沙北侧建设4道丁坝、南侧建设8道丁坝，优化双涧沙头部潜堤结构。2017年3月二期工程完工后，由于底沙下泄和弯道环流输沙影响，仍需对航槽边缘及航槽内淤浅部位进行维护以达到12.5 m深水航道通航尺度[10-11]。研究福姜沙水道的冲淤特征，并揭示其在航道工程作用下的滩槽演变特性显得尤为重要和紧迫。
众多学者基于实测资料分析、物理模型试验、数值模拟的方法，从河流动力学和河流地貌学的角度对福姜沙水道的滩槽演变特性进行了研究[1, 3-4, 6-14]。研究表明，水流动力轴线摆动是河床调整变化的主要动力因子[12]，边心滩变形、航道回淤与特征流量级的大小和持续时间有关[14]，河道整体冲淤与上下游河道演变密切相关[15]。已有研究成果通常基于定性的描述，且研究背景多集中于二期工程实施前，对航道工程实施后的河道演变特性研究较为匮乏。基于实测地形、水文资料分析表明，航道工程实施前后河道的演变特征发生了改变[16]，正在向与工程相适应的河势调整，福姜沙左缘边滩和福姜沙沙体体积增大，福中水道深槽扩容，福北水道深槽淤浅[11]。过去十余年，水沙数学模型被大量应用于福姜沙水道的整治方案设计、航槽回淤、滩漕演变规律[12]等问题的研究[4, 12, 17-21]。基于TVD—有限体积法的水沙模型发展已有二十余年的历史[22-26]，但其在野外尺度的应用非常少见。这主要是因为：第一，TVD—有限体积法为显式格式，计算时间步长受限于CFL稳定性条件；
第二，为了编程的方便，以往数学模型习惯上采用满足CFL条件的全局最小时间步长。当计算网格尺度或水流状态不均匀时，整体时间步长被最小时间步长限制，大大降低模型的计算效率[27-28]。对于清水模型，局部时间步长方法(LTS)可以使每个单元的变量采用满足稳定性条件尽可能大的时间步长进行更新，从而降低计算耗时[29-33]。进一步采用混合局部时间步长和整体最大时间步长方法(Hybrid LTS/GMaTS)将局部时间步长技术推广至动床模型中，已应用于长江中游水道的冲淤数值计算[34]，但尚未应用于径潮流共同作用的水道。
采用基于Hybrid LTS/GMaTS方法的平面二维水—沙—床耦合数学模型，以福姜沙水道径潮流和河床冲淤过程为例，验证其计算效率和精度。通过8组数值算例，分析二期工程实施后，福姜沙水道在不同流量特征下的水动力特征和河床冲淤特点，尝试从洲滩特征长度变化的角度定量揭示洲滩演变对径流特征的响应过程，在此基础上建立航道工程影响下航道冲淤量与洲滩演变、径流特征的数学关系式，以加深对福姜沙水道滩槽演变规律的认识，为航道工程提供参考。
1.1 控制方程
水—沙—床耦合数学模型的控制方程[35]为：
(1)
(2)
(3)
式中：t为时间；
x和y为笛卡尔坐标系下的空间坐标；
h为水深；
u和v分别为x和y方向深度积分平均的流速分量；
g为重力加速度；
ck为第k粒径组泥沙的深度平均体积含沙量；
zb为床面高程；
δ为床面活动层厚度；
fa,k为活动层第k粒径组泥沙体积百分比；
Ek、Dk为床面与水流之间第k粒径组泥沙交换的上扬和沉降通量，DT=∑Dk，ET=∑Ek；
p0为泥沙孔隙率；
Sbx和Sby分别为x、y方向的底坡坡度；
Sfx和Sfy分别为x、y方向的阻力坡度；
fs,k为底床存储层和活动层交界面处第k粒径组泥沙体积百分比；
η=zb-δ为活动层下界面的高程。
1.2 经验公式
阻力坡度采用曼宁公式计算：
(4)
(5)
其中，n为曼宁糙率系数。
底床与水流之间的泥沙上扬和沉降通量按式(6)～(7)计算：
Dk=αkckωk(1-αkck)m
(6)
Ek=αkce,kωk(1-αkce,k)m
(7)
式中：αk表示近底泥沙含量与深度平均含沙量之间的比值，考虑到该河段引起底床冲淤变化主要是临底层的高浓度泥沙[13]，αk取3.0；
ωk为第k粒径组单颗粒泥沙沉速，采用张瑞瑾统一沉速公式计算；
ce,k为第k粒径组泥沙的水流挟沙力。采用张瑞瑾水流挟沙力计算公式计算[36]：
(8)
(9)
式中：υ为水流运动黏滞系数，取υ=1.2×10-6m2/s;dk为第k粒径组泥沙的特征粒径；
s=(ρs-ρw)/ρw，ρs、ρw分别为清水和泥沙密度；
K0、m0为经验参数，本文取K0=0.12，m0=0.92。
1.3 基于Hybrid LTS/GMaTS的数值算法
采用非结构三角形网格如图1所示。图1(a)为计算区域某单元以及三个相邻单元；
图1(b)为某个角点以及其所属的三角形单元；
图1(c)为某个界面及其左右两侧的单元。单元变量使用Hybrid LTS/GMaTS进行更新，对控制方程采用有限体积法进行离散，挟沙水流的质量与动量守恒方程、随流泥沙质量守恒方程[式(10)]采用局部时间步长(LTS)求解，底床所有组分泥沙总质量守恒方程[式(11)]、床沙分粒径组的质量守恒方程[式(12)]采用整体最大时间步长(GMaTS)求解[34]：
图1 非结构三角形网格Fig. 1 Sketches of the unstructured triangular meshes
(10)
(11)
(12)
2.1 模型设置
计算区域如图2所示。模型采用非结构化网格，对航道工程处进行局部加密，网格尺度最大为260 m，最小为56 m，平均为135 m。选取2017年2月实测地形作为初始地形，采用实测床沙级配数据。悬沙粒径分为6个组别：0.042 mm、0.07 mm、0.085 mm、0.11 mm、0.16 mm和0.35 mm，各个粒径的组分根据2017年8月福姜沙水道所有断面悬移质实测数据加权平均得到。河道糙率采用经验公式n=0.01+0.01/H计算[38]。上游边界取在江阴河道内，采用大尺度模型为上游边界提供流量过程，来沙采用大通站实测含沙量，下游边界位于徐六泾水文站断面，采用实测水位过程。
图2 计算区域和初始地形Fig. 2 Location of the study area and initial bed topography
2.2 Hybrid LTS/GMaTS方法对计算结果的影响
使用均方根误差(εRMSE)量化Hybrid LTS/GMaTS方法与传统最小时间步长方法之间计算结果的相对误差，使用全局泥沙质量的相对误差ε量化泥沙计算的守恒性：
(13)
(14)
式中：f表示水位、流速、河床高程变化等物理量；
N为网格节点数；
fLG+GMa和fGMi分别表示Hybrid LTS/GMaTS方法和传统最小时间步长方法的计算结果；
V(t2)和V(t1)分别表示计算区域内t2和t1两个时刻水体中泥沙的总体积；
Vin和Vout分别表示两个时刻之间从边界流入或流出的泥沙体积；
VE和VD分别表示两个时刻之间河床与水体泥沙交换过程中上扬和沉降的泥沙体积。
模拟2017年8月15日至18日的水动力和河床冲淤过程，考虑muser=1～7的七个工况(muser是使用者给出的参数，用于限制LTS层级，muser越大，整体最大时间步长能达到的值越大)，为方便计算误差，将最小时间步长设置为1 s。Hybrid LTS/GMaTS方法计算效率和相对误差如表1所示。当muser=0时得到传统整体时间步长方法，即每个网格的时间步长取全局最小值。Hybrid LTS/GMaTS方法可以大幅度提高计算效率，计算耗时随着muser的增大而减小。Hybrid LTS/GMaTS方法(muser=7时)相比于传统算法提升近十倍。传统方法计算结果与实测数据的相对误差为0.09 m(水位)、0.12 m/s(流速)、1.64 m(床面高程变化，2017年2月至2018年2月)，而Hybrid LTS/GMaTS方法相对于传统方法的最大误差仅12.93×10-4m(水位)、9.67×10-4m/s(流速)、7.01×10-4m(床面高程变化)，远小于传统方法计算结果与实测值之间的相对误差。该方法在大幅度减少计算成本的基础上带来的计算误差可忽略不计。Hybrid LTS/GMaTS方法(muser=7)的全局泥沙质量相对误差量级仅为10-9，说明泥沙计算具有良好的守恒性。
表1 Hybrid LTS/GMaTS方法计算效率和相对误差Tab. 1 Calculation cost and relative error using Hybrid LTS/GMaTS method
2.3 水动力及含沙量验证
根据2017年洪季水文测验结果，验证福姜沙水道至通州沙水道沿岸四个潮位站(2L、4R、6R、任港)的水位(图3)，四个测点(FZ-C、FZ1-B、JLG-A、TSG-B)的流速、流向和含沙量(图4、图5)。测站和测点位置如图2所示。对比结果显示，模型计算的水位、流速、流向、含沙量与实测过程线趋势基本一致；
水位计算结果与实测值的均方根误差为0.07 m，计算水位与实测值的偏差均在0.1 m内；
流速计算结果的均方根误差为0.07 m/s，考虑到实测平均流速为0.75 m/s，平均流速偏差在10%以内；
流向计算结果的均方根误差为9°，平均流向偏差在10°以内；
含沙量计算结果的均方根误差为0.02 kg/m3，考虑到实测平均含沙量为0.10 kg·m3，平均含沙量偏差在30%以内。本模型的水动力和含沙量验证满足《水运工程模拟试验技术规范》(JTS/T 231—2021)的规定。
图3 计算水位与实测值对比Fig. 3 Comparison of calculated and measured water levels in different stations
图4 计算垂向平均流速流向与实测值对比Fig. 4 Comparison of calculated and measured vertical average current speeds and directions
2.4 河床冲淤验证
图6分别为福姜沙水道2017年2月至2018年2月实测和计算冲淤分布。对比结果表明，模型成功复演了本河段2017年2月至2018年2月的滩漕冲淤过程的主要特征：双涧沙整体淤涨，两侧深槽冲刷，靖江边滩冲刷后退，计算冲淤部位和冲淤趋势与实测结果相似。该河段实测冲淤量分别为4 953万m3和4 391万m3，计算冲淤量分别为4 296万m3和4 152万m3，冲淤量的偏差分别为13%和5%，冲淤量误差均控制在20%内，满足《海岸及河口潮流泥沙模拟技术规程》(JTS/T 231-2—2010)，说明本模型的泥沙模块满足精度要求。
图6 2017年2月至2018年2月地形冲淤变化Fig. 6 Topographic changes from February 2017 to February 2018
3.1 工况设置
长江径流流量存在洪枯季，且每时都发生变化，河道内滩槽演变过程是对不断变化的径流流量的综合响应。径流过程对于滩槽演变的综合作用大致相当于某一特征流量级[39]。如表2所示，选取大通站2003年后多年流量累计频率为5%(59 800 m3/s)、25%(40 500 m3/s)、50%(25 000 m3/s)、75%(15 000 m3/s)作为洪、中、枯水期的特征流量，选取大通站2003年后多年洪、中、枯水期平均含沙量(分别为0.24 kg/m3、0.14 kg/m3、0.09 kg/m3)作为来流特征含沙量。下游边界条件考虑径流和潮汐的影响，将水位过程分解为日平均水位和潮波水位[40]：特征流量与基准流量之差作为径流增量对日平均水位进行预报[41]，徐六泾站实测潮位序列做调和分析对潮波水位进行预报。
表2 工况设置Tab. 2 Simulation cases
选取福姜沙水道2017年2月实测地形作为模型计算初始地形。一方面，此时二期工程已接近竣工，考虑航道工程对水道滩槽演变的影响较为合理。另一方面，水流切割下来的滩体淤积在六助港附近，福北水道航道条件有恶化趋势[42]，历史上靖江边滩多次发生切割过程，切割下来的沙体呈现形成、冲刷下移、并入福北水道及双涧沙、再形成的周期性过程[15]，该地形下河道演变过程对航道条件的影响具有代表性。每组工况中设置两个算例，在第一个算例中考虑航道工程措施，包括整治建筑物和基建性疏浚工程，如图7(a)所示。疏浚区域见图7(a)，设计挖槽深度为理论最低潮面以下12.5 m，超深取0.5 m；
第二个不考虑航道工程(即自然状态)作为对比算例。8组工况的模拟时长均为90 d。
3.2 不同流量级条件下水动力及冲淤特征
流速大小反映了水流动力作用的强弱，河流断面或沿程流速分布不均，导致水流挟沙力差异，使河床产生冲淤变化，从而影响河流的滩槽演变特性。水流经江阴水道流入福姜沙水道时，干流受到左岸炮台圩节点和右岸鹅鼻嘴节点共同控制[42]。该河段内落潮流是塑造主槽河床的主要动力[4]，图8是不同流量特征下各断面大潮落潮流流速分布，断面位置见图7(a)。枯水期(流量为15 000 m3/s)大潮落潮流主流位于左汊深槽处，左汊深槽处水流流速大于左岸浅滩，靖江边滩位于缓流区而淤积。随着流量的升高，水流逐渐脱离炮台圩节点摆动至右岸，对鹅鼻嘴节点的顶冲作用增大，节点挑流作用增强，水动力轴线逐渐向左岸边滩摆动。当流量大于40 500 m3/s时，左岸浅滩处水流流速大于左汊深槽，大潮落潮流主流偏向靖江边滩处，直接作用在靖江边滩上，使靖江边滩冲淤加剧，如图7(c)、(d)所示。因此，该河段进口干流具有“大水趋直，小流坐弯”[4]的水流特性：枯水时主流在边滩束水作用下走深槽，中洪水时主流线偏靠边滩而取直，特别是在大流量(59 800 m3/s)下，边滩呈现头部冲刷，中部淤积，尾部冲刷的冲淤特点，说明边滩在主流直接作用下出现冲刷后退的现象，有进一步被切割成心滩的风险。以上水流动力轴线在不同流量级条件下的摆动造成靖江边滩洪冲枯淤的演变特性。而各流量级下航道工程对于福姜沙入口处水动力影响较小，在福姜沙左缘边滩丁坝群作用下，深槽及左岸浅滩处流速略有增加，最大增加约0.1 m/s，右岸浅滩处流速有所减小，减小幅度约为0.05 m/s。
图7 各工况下地形冲淤计算结果Fig. 7 Calculated results of topographic erosion and deposition
图8 各断面大潮落潮流流速分布Fig. 8 Cross-section ebb velocity distribution of spring tide
主流经双涧沙分为福北、福中水道两汊，主流进入福中水道，分流比约为70%[17]，如图7(b)～(d)所示，在自然状态下，流量为15 000 m3/s时，福中水道主流位于靠近福姜沙左缘边滩的深槽内，福北水道主流于入口处位于偏左岸的深槽内，而后偏向双涧沙一侧。随着流量增大，入口处主流摆至靠近双涧沙沙头，当流量超过40 500 m3/s时，主流进一步摆动至双涧沙上部，而福中水道水流动力轴线摆向福姜沙北缘浅滩、福北水道水流动力轴线则摆动至双涧沙左缘。这种主流摆动特征与两汊河道呈现洪季“凸岸冲刷、凹岸淤积”、枯季“凹岸冲刷、凸岸淤积”[42]的演变特征一致，进而导致双涧沙沙头“洪冲枯淤”、双涧沙两侧浅滩在不同流量级下冲淤交替的演变特征，带来双涧沙滩体左右摆动的不利影响。在福姜沙左缘丁坝和双涧沙丁坝群的共同作用下，即使流量超过40 500 m3/s，福中水道主流也控制在主槽内，仅在丁坝挑流作用下向右岸浅滩一侧轻微摆动。各流量级下福北水道入口处水动力条件均得到显著改善，航槽刷深显著。流量为15 000 m3/s时，福中水道中段总体呈现淤积趋势，特别是福中、福南水道交汇处受支汊入流顶托效应影响，淤浅尤为显著。然而随着上游来流流量增大，深槽处水动力增强，福中水道总体上转淤为冲，呈现出“洪冲枯淤”的特点。福北水道入口处流速相比于自然状态下反而有所降低，受到双涧沙北缘丁坝群挑流作用，福北水道主流稳定在深槽内，随着流量增大，水流动力轴线偏向左岸，是靠近双涧沙一侧深槽转冲为淤[图7(a)～(d)]的主要原因。上游来流流量增大、入口处水动力作用显著增强的情况下，淤积反而增加，说明水动力作用不是影响福北水道冲淤变化的主要因素。作用于双涧沙滩体上的水动力条件减弱，说明滩面流在航道工程作用下得到控制，且位于丁坝坝田区的滩体在各流量级下均呈现淤积趋势，特别是洪季双涧沙滩头冲刷趋势得到抑制，双涧沙滩体形态趋于稳定，对于两汊的发展产生有利影响。
3.3 不同流量级条件下滩槽演变特征
3.3.1 滩槽演变特征
将12.5 m等深线(基于理论最低潮面)作为划分滩槽的依据，绘制各流量级下滩槽分布图，如图9所示。流量为15 000 m3/s时，在自然状态下跃进港附近边滩冲刷萎缩，而蟛蜞港附近边滩淤积展宽，说明在水流作用下，靖江边滩整体上呈现冲刷后退的趋势，并随着流量的增大，冲刷趋势增强。当流量达到59 800 m3/s时，靖江边滩尾部深入主流的低矮滩体被水流切割脱离边滩，进入航槽形成心滩，符合靖江边滩在高于60 000 m3/s流量下完成切割过程的历史演变特征[13]。中枯水时期，双涧沙沙舌形态几乎不发生变化，但当流量大于40 500 m3/s，双涧沙沙舌受水流冲击后退。福姜沙左缘边滩和双涧沙右缘浅滩在各流量级下滩面形态几乎没有发生大的变化，仅在福姜沙尾部有向福南水道左岸浅滩淤积展宽的趋势，因此福中水道航槽条件无明显影响。初始地形下福北水道夏仕港附近12.5 m深槽有中断，随着流量增大，中断长度有所减小。但福北水道深槽在洪季(流量大于40 500 m3/s)发生较大程度的变形，左岸浅滩向内侧挤压航槽，而双涧沙左缘浅滩则受冲刷后退，福北水道航槽整体有向南摆动的趋势，对于航槽稳定性及左岸港口水深产生了不利影响。
图9 各工况下滩槽格局变化Fig. 9 Change in morphology of bars and channels
在航道工程作用下，靖江边滩尾部和福北水道两侧浅滩的滩槽演变趋势产生了根本性变化。在一次性疏浚条件下，福北水道12.5 m航槽被打通。夏仕港至焦港航槽在各流量级下均保持稳定，且在双涧沙左缘丁坝的控制下，12.5 m深槽有向左岸展宽的趋势。焦港至如皋港段水流脱离双涧沙航道工程控制，双涧沙一侧浅滩向内挤压深槽，随着流量增大，深槽内回淤程度加剧。双涧沙沙舌向福北水道航槽内淤涨，当流量达到40 500 m3/s时，双涧沙沙舌与淤积在六助港附近的浅滩相连，依靠疏浚打通的12.5 m深水航槽被阻断。当流量达到59 800 m3/s时靖江边滩尾部深入主流的低矮滩体没有被水流切割，而是继续向航槽内发育，在福北水道入口处形成浅滩，进一步影响福北水道的航行条件。说明航道工程建成后，靖江边滩的冲刷下移和双涧沙沙舌的横向发育是影响福北水道入口处航槽稳定的主要因素。
3.3.2 洲滩演变及特征流量对航道条件的影响
将靖江边滩垂直流向最宽处至炮台圩节点的直线距离作为靖江边滩特征尺度，用于表征靖江边滩冲刷后退的程度，绘制靖江边滩特征尺度关于特征流量的散点图，如图10(a)所示。将双涧沙沙舌垂直流向最宽处作为双涧沙沙舌的特征尺度，用于表征双涧沙沙舌横向发育的强弱，绘制双涧沙沙舌特征长度关于特征流量的散点图，如图10(b)所示。靖江边滩特征长度和流量级呈现正相关关系，且流量大于40 500 m3/s时增长速率显著提高，靖江边滩会出现较大规模的下移，航道工程下靖江边滩下移速率大于自然状态。双涧沙沙舌在自然状态下与流量级相关较弱，航道工程下特征尺度随着流量增大而快速增长。通过以上分析，一方面说明年内流量过程中洪季流量对靖江边滩和双涧沙沙舌变形起到主导作用，另一方面，航道工程实施后，加快了靖江边滩冲刷后退、双涧沙沙舌横向发育的速率，福北水道入口处水动力条件减弱，不利于冲刷下移的泥沙通过福北水道向下游输移，造成局部区域形成碍航浅滩。
图10 靖江边滩及双涧沙沙舌特征尺度随流量级变化关系Fig. 10 Relationship between discharge and characteristic scale of point /channel bar deformation
靖江边滩—双涧沙沙舌—福北水道是主要的输沙通道，基于2014年7月至2016年11月的实测数据，王建军等[13]建立了福北水道航槽冲淤量V福北水道航槽与靖江边滩冲淤量V靖江边滩、双涧沙沙舌冲淤量V双涧沙沙舌及径流量V径流量的关系，见式(15)。该拟合式的实际意义是：靖江边滩和双涧沙沙舌冲刷下来的泥沙进入福北水道航槽并淤积，且随着径流量的增加航槽淤积量增多。将计算结果代入式(10)进行验证，计算结果总体上符合拟合值的变化趋势(相关性系数R=0.932 1)；
而航道工程影响下的计算值与拟合结果偏差较大，且与拟合值呈现负相关(相关性系数R=-0.804 3)，反映了航道工程完工后，靖江边滩、双涧沙沙舌、福北水道的滩槽演变关系和对径流特征响应过程发生根本性变化。式(15)不适用于描述航道工程影响下福北水道航槽冲淤量与靖江边滩和双涧沙沙舌、径流特征的关系。
根据算例A1、B1、C1、D1(考虑航道工程)的计算结果，可以得到航道工程影响下的福北水道航槽冲淤量与靖江边滩和双涧沙沙舌冲淤量、径流量的拟合关系式(16)(R=0.849 8)。在靖江边滩、双涧沙沙舌冲淤及特征流量的综合作用下，靖江边滩和双涧沙沙舌每冲刷约100万 m3泥沙，约11.9万 m3会淤积在福北水道航槽内；
径流量每减少1 000亿 m3，福北水道航槽淤积量会增加约17万 m3。如图11所示，航道工程作用下福北水道的淤积量随着流量的增长而减小，特别是洪水期转淤为冲，而维护性疏浚量反而增大，说明对航道整体而言淤积程度放缓，洪季航槽大规模回淤的现象得到显著改善，但受到靖江边滩冲刷下移、双涧沙沙舌横向淤宽和局部区域水动力减弱的不利影响，随着流量增长，福北水道入口处区域的航道条件发生恶化，应作为深水航道维护的重点关注区域。
图11 福北水道整体冲淤量和维护性疏浚量Fig. 11 Total erosion and deposition volume and maintenance dredging volume of Fubei channel
V福北水道航槽=-0.348(V靖江边滩+V双涧沙沙舌)+0.033V径流量+93.82
(15)
V福北水道航槽=-0.119(V靖江边滩+V双涧沙沙舌)-0.017V径流量+29.383
(16)
采用基于Hybrid LTS/GMaTS方法的水沙—床耦合数学模型对工程河道冲淤过程进行数值模拟研究。模型使用非结构三角形网格，对航道工程区域进行局部加密，通过Hybrid LTS/GMaTS方法提高了计算效率。Hybrid LTS/GMaTS方法和传统整体最小时间步长方法的水位、流速、含沙量、床面高程变化的计算结果相对误差较小，但该方法大幅度提高了计算效率，在本文计算工况下提高了近10倍。模型复演了2017年7月17日至18日的径潮流水动力过程和2017年2月至2018年2月的河床冲淤过程，说明该模型可以应用于长江感潮河段的水沙数值模拟。基于8组数值算例，得到以下结论：
1)水动力作用是福姜沙水道洲滩变形的主要动力，体现在洲滩整体冲淤和局部形态改变两方面。一方面，靖江边滩受主流摆动影响冲淤交替，枯水期处于缓流区而淤积、洪水期处于主流区而冲刷；
双涧沙在航道工程作用下水动力条件减弱，整体上呈现淤涨趋势。另一方面，靖江边滩受水动力作用向下游移动，双涧沙沙舌受水力输沙作用横向淤宽，洪季流量对靖江边滩和双涧沙沙舌变形起到主导作用。
2)航道工程作用下，福中水道入口处水动力条件明显改善，各特征流量级下均呈现刷深趋势；
航道中段呈现“洪冲枯淤”的冲淤特征；
航道工程作用下福中水道的主流被归顺在深槽内，仅在丁坝挑流作用下轻微向福姜沙左缘边滩摆动，因此随着流量增长，航槽向南摆动。
3)边心滩演变和水动力作用共同影响福北水道航道冲淤变化。航道工程建成后，靖江边滩、双涧沙沙舌、福北水道的滩槽演变关系和对径流特征响应过程发生根本性变化。建立了航道工程影响下的福北水道航槽冲淤量与靖江边滩和双涧沙沙舌、径流特征的拟合关系式。在靖江边滩、双涧沙沙舌冲淤及特征流量的综合作用下，靖江边滩和双涧沙沙舌每冲刷约100万 m3泥沙，约11.9万 m3泥沙会淤积在福北水道航槽内；
径流量每减少1 000亿 m3，福北水道航槽淤积量会增加约17万m3。航道工程实施后，对福北水道整体而言，其淤积程度放缓，特别是洪季航槽大规模回淤得到显著改善。然而受到靖江边滩冲刷下移、双涧沙沙舌横向淤宽和局部区域水动力减弱的不利影响，随着流量增长，福北水道入口处航道条件将发生恶化，是重点疏浚维护区域。