河口海岸泥沙输移的分组数学模型研究.pdf

河口海岸泥沙输移的分组数学模型研究 杨 晨 1,路 宽2 (1.中国石油大学 石油工程学院,北京 102249；2.国家海洋技术中心,天津 300112) 摘要：根据河口海岸水沙输移的特点,建立了一个新的二维分组数学模型,用来预测该区域的 水沙输移过程。该模型耦合了水动力模块、泥沙输移模块和床面演变模块。其中水动力模块基于 浅水方程组,综合考虑了柯氏力、床面切应力以及表面风应力的影响,引入干湿判断法处理动边 界。 泥沙输移模块首先将泥沙按照粒径分组,针对不同泥砂性质,对各组泥沙分别进行建模求 解。 床面演变模块基于质量守恒方程,实时更新床面高程以及床沙级配变化,并传递给水动力模 块,更新底部边界。 该模型被应用在了英国塞汶(Severn)河口,其预测的泥沙浓度和实测数据以 及不分组的模型的预测结果进行了比较,结果显示,文中建立的分组模型预测的结果要明显好 于不分组模型。 关键词：河口海岸；泥沙输移；分组模型；地形演变；水动力 中图分类号：TV14文献标志码：A文章编号：1003-2029(2015)02-0101-05 收稿日期： 2014-09-29 基金项目： 国家自然科学基金资助项目 （51304227） ； 清华大学水沙与水利水电工程国家重点实验室开放基金资助项目 （sklhse-2013-B-03） 作者简介： 杨晨 （1981-） 男， 博士， 讲师， 主要研究方向为水沙多相流的数值模拟与实验研究。E-mail： yangchen 河口海岸地区是陆地与海洋相互作用的交流 地带， 它既是海产品、 石油和天然气等资源的产地， 又为海洋化工产品的生产提供原料；另一方面， 由 于河口独特的地理位置， 又通常会形成大面积的天 然湿地， 不仅可以调蓄洪水、 改善气候， 而且可以维 持生物多样性， 维持区域生态平衡， 对地球的环境 和生态健康都有着巨大的影响。河口海岸地区一般 水深较浅， 水流和泥沙的相互作用对一些水环境和 水生态问题都有至关重要的影响， 比如重金属的输 移、 细菌和微生物的输移和衰减等。所以， 掌握河口 海岸地区水沙输移规律， 准确预测河口海岸地区的 泥沙浓度， 一直是广大学者重点关注的科学问题。 河口海岸地区的泥沙粒径分布并不均匀， 粗细 之间会相差几十倍以上， 而且粘性沙和非粘性沙都 同时存在1。所以， 对这个区域泥沙输移进行数值模 拟时， 有必要将泥砂区别对待， 对不同性质的泥沙 分别建模求解。泥沙分组模型的想法最先由 Ein- stein2提出， 随后 Misri 等3建立了分组泥沙输移的 概念模型， Samage 等在概念模型的基础上进一步完 善， 在大量观测数据的基础上， 提出了分组泥沙输沙 率计算的经验公式4。James 等利用商业软件 EFDC 对 180 度转弯处不均匀泥沙输移实验进行了模拟， 计算结果与实测结果相近5。吴挺峰等建立了基于 斜压水动力模型的垂向二维泥沙模型， 该模型可以 独立求解粘性沙和非粘性沙的沉降速度及水土界 面的交换通量， 模型对富春江水库的水沙输移情况 做了预测， 模拟结果与实测结果接近6。 Yang 等建立 了考虑地形演变的分组泥沙输移模型， 在黄河口应 用并取得了良好的效果7。 本文在前人研究的基础上， 提出了可以模拟河 口海岸泥砂输移的分组数学模型。该模型的水动力 部分考虑了柯氏力、 床面切应力以及表面风应力的 影响， 引入干湿判断法处理动边界问题。泥沙输移 部分将泥沙按照粒径进行了分组模拟， 并且综合考 虑了粘性沙和非粘性沙不同的运动规律。床面演变 部分实现了对床面高度和床沙级配变化的模拟， 并 实时传递给水动力模块来更新计算边界。模型被应 用在了英国的塞汶 （Severn） 河口， 本模型预测的泥 第 34 卷 第 2 期 2015 年 4 月 海 洋 技 术 学 报 JOURNAL OF OCEAN TECHNOLOGY Vol. 34,No.2 Apr,2015 第 34 卷海 洋 技 术 学 报 沙浓度和实测数据与不分组的模型的预测结果进 行了比较， 比较结果显示本文建立的分组模型的预 测结果要明显优于不分组模型。 1数学模型 河口海岸地区由于水体垂向尺度相对于横向 尺度较小， 沿着水深可以有静水压强假定， 纵向速 度和加速度可以忽略， 所以把 Navier-Stokes 方程沿 着水深积分， 得到如下浅水方程组: 鄣 鄣t + 鄣p 鄣x + 鄣q 鄣y =qm（1） 鄣p 鄣t +鄣pU 鄣x +鄣pV 鄣y =fq-gH鄣 鄣x + a CwWxW 2 x+W 2 y - gpp2+q2 H 2C2 + 2鄣 2p 鄣x2 +鄣 2p 鄣y2 + 鄣 2q 鄣x鄣y （2） 鄣q 鄣t +鄣qU 鄣x +鄣qV 鄣y =-fp-gH鄣 鄣y + a CwWyW 2 x+W 2 y - gpp2+q2 H 2C2 + 鄣2q 鄣x2 +2鄣 2q 鄣y2 + 鄣2p 鄣x鄣y （3） 式中： U， V 分别为 x 和 y 方向深度平均速度； p 和 q 分别为 x 和 y 方向的单宽流量； qm为单位面积 上的源流量； 为速度纵剖面上的动量修正因子； f 为柯氏力参数； H 为总水深； 为基准面之上的水面 抬升高度； h 为基准面之下的水深； a为空气密度； 为流体密度； C 为谢才系数，可用 Manning 公式计 算； Cw为空气与水交界面的阻力系数； 为深度平 均的涡粘性； Wx和 Wy分别为水面上 10 m 处的 x 和 y 方面的流速8。 对于动边界问题， 可采用干湿判断法处理9。将 计算网格分为干网格和湿网格，对水位降落过程， 进行干网格检索。用 来描述粗糙高度或沙波高 度， 检索过程如下： （1）当单元网格任意边中心点水深小于 时， 其水深和流速都为 0， 没有流动经过这个单元格。 （2）当单元网格周围四点水深都小于 时， 那 么该网格为干网格，将这个网格从计算区域删除， 本网格的水位保留之前半个时段的计算数值。 （3）一个湿网格周边四个边中心水深的平均值 小于 ， 该网格也设为干网格。 （4）如果一个湿网格四个边中心水深的平均值 不小于 ， 但小于 2， 该网格假设为可能干网格： 如 果网格周边四个中心水深至少一个介于 和 2 之 间， 且流动是从该网格发散出去， 那么该网格可以 设定为干网格； 相反该网格不能从计算区域删除。 对于水位上升过程， 在下列情况下干网格需要 重新进入计算区域： （1）当一个干网格周围水位高于干网格的水位 时， 该干网格变为湿网格， 恢复到计算区域中。 （2）当干网格与周围湿网格之间的交界面水深 大于 时， 该干网格变为湿网格。 （3）干网格中心的水深 H 时。 河口和海岸地区的泥沙粒径分布十分不均匀， 所以首先将泥沙按照粒径分组，对于每组泥沙， 其 在水中的输移可用深度平均的带源项的对流扩散 方程进行描述10： 鄣i 鄣t + 鄣Ui 鄣x + 鄣Vi 鄣y = 1 H 鄣 鄣x DxxH 鄣i 鄣x +DxyH 鄣i 鄣y + 1 H 鄣 鄣y DyxH 鄣i 鄣x +DyyH 鄣i 鄣y +Ssi+Sai（4） 式中： i是第 i 组泥沙深度平均的浓度； Ssi是第 i 组泥沙包含沉降和冲刷的源项； Dxx和 Dyy分别是 x 和 y 方向水深积分的横向扩散系数， 附加源项 Sai可 以通过 Sai=i 鄣U 鄣x + 鄣V 鄣y ()计算得到。 对于非粘性沙， 方程 （4） 中的源项 Ssi可由以下 公式计算11： Ssi=si(aei-ai) /H（5） 式中： wi是第 i 组泥沙的沉降速度； aei是第 i 组泥沙在参考点处的平衡浓度； ai是第 i 组泥沙在 参考点附近的实际泥沙浓度8。 对于粘性沙， 源项 Ssi可表示为： Ssi= （Ei-Di） /H（6） 式中： Di是沉降泥沙通量； Ei是冲刷泥沙通量8。 床面高程变化的控制方程其实就是泥沙的质 量守恒方程，对于做悬移质运动的每组粘性泥沙， 相应的床面高程变化为： s Zbi t - 1 1-p0 (Di-Ei)=0（7） 式中： p0是床沙的孔隙率。 对于做悬移质运动的每组非粘性泥沙， 相应的 床面高程变化为： s Zbi t - 1 1-p0 i(i-ei)=0（8） 102 在一个时间步长内， 总的床面高程变化为： Zb= N i=1 Zbi（9） 在泥沙的输移过程中， 床沙级配的变化很大， 假 设在一个时间步内，床底高程由 HT变化到 HT+Zb， 则第 i 组泥沙在床沙中所占的体积百分比变为： P f i= HTP b i+Zbi HT+Zb （10） 式中： Pi是第 i 组泥沙在床沙中的体积分数， 上 标 b 和 f歩分别代表时间之前和之后的值。 2模型应用 为了验证模型的适用性， 模型被应用在了英国 的塞汶 （Severn） 河口。如图 1 所示， Severn 河口位于 英国的威尔士和英格兰之间的西南海， 该区域有着 世界第二高的潮差， 泥沙浓度有时超过 1 000 mg/L， 是英国最大的河口之一。 因为在研究区域内，一维和二维流动同时存 在， 所以模型采用了两套网格体系， 在河道内采用 一维网格 （1D） ， 在河口区域采用二维网格 （2D） 。如 图 2 所示， 二维部分被划分为了 242168 个正方形 网格， 每个网格均为 600600 m， 一维区域 （从 Sev- ern Bridge 到 Haw Bridge） 被划分为了 351 个网格， 平均每个网格长度为 240 m。二维区域的向海边界 设为水位开边界， 潮位数据可由 Proudman Oceano- graphic Laboratory (POL) 的潮流模型获得。一维区 域的上游边界 （Haw Bridge） 设置为流速边界条件。 其中，二维模型和一维模型是交互耦合的关系， 二 维模型为一维模型提供水位数据， 一维模型为二维 模型提供流量数据。 应用本文建立的模型对区域内的水沙输移进 行了数值模拟，模拟时间为 300 h，包括了一次涨 潮、 一次小潮以及一次中潮， 时间步长设为 25 s。模 型中的动量修正因子取为 1.016， 2D 区域的床底糙 率高度设定为 35 mm， 1D 区域的 Manning 粗糙系数 取为 0.022。 模型对整个区域的的水动力过程进行了模拟， 并在 S Wales 观测站处与观测数据进行了对比（如 图 2 所示） 。通过比较可以看出，本模型计算的水 深、 流速以及流向都与观测数据吻合得非常好。 水动力模块被验证之后， 应用本模型对区域内 的泥沙输移过程进行了数值模拟， Severn 河口的泥 沙粒径的几个特征值 D16， D50， D84， D90分别为 0.026 mm， 0.058 mm， 0.126 mm 以及 0.15 mm13， 粘性沙的 平均直径介于 0.01 mm 和 0.063 mm 之间。利用这 些特征值将泥沙分成 5 组，每组泥沙的性质如表 1 所示。由分组结果可以看出， 第 5 组泥沙的平均粒 径比第 1 组泥沙粒径的 10 倍还要大，这也证明了 采用分组模型的必要性。 图 1 研究区域及各观测点位置示意图12 图 2 本模型计算的水深， 流速以及流向与实测数据比较 （S Wale 观测点, 2001 年 7 月 24 日） 时间/h 时间/h 时间/h 杨 晨,等：河口海岸泥沙输移的分组数学模型研究第 2 期 103 第 34 卷海 洋 技 术 学 报 模拟结果和现场实测数据在 Southerndown 观 测点进行了对比，为了验证本文分组模型的优势， 同时也和 Gao 等人利用不分组模型的计算结果12做 了对比。对比结果如图 3 所示， 本模型的计算结果 和实测结果更加吻合， 而不分组模型计算的泥沙浓 度要明显小于实测浓度， 这主要是因为不分组模型 用 D50作为代表粒径来模拟泥沙运动，所以粗沙和 细沙的运动实际上被隐藏了， 而水体里细沙的浓度 是占主要部分的， 由于这部分被隐藏， 所以会导致 计算结果的偏低。如果要依据模拟结果制定控沙政 策， 这样明显偏低的结果将是十分危险的， 须引起 足够的重视。 另外， 本模型的计算结果会比不分组模型结果 滞后， 这可能是因为不分组模型无法考虑床沙级配 的变化， 使得床沙被 “过快” 冲起所导致的。 3结 论 本文建立了可以模拟河口海岸水沙输移过程 的二维数学模型。该模型的水动力模块综合考虑了 柯氏力、 床面切应力以及表面风应力的影响， 引入 干湿判断法处理动边界。泥沙运动模块应将泥沙按 照粒径进行分组， 考察每一组泥沙的性质然后分别 建模求解。床面演变模块可以模拟床面高程的变化 和床沙的级配变化。 模型被应用在了英国的 Severn 河口， 通过和实 测数据以及不分组模型结果的比较， 体现了本模型 在水沙数输移值模拟方面的优越性。模型可以被用 于研究河口海岸水沙输移的科学， 也可以为管理者 提供制定政策的依据。 表 1 各组泥沙的平均粒径以及质量百分比 组别12345 平均粒径/mm0.0180.0420.0920.1380.190 百分比/%163434610 图 3本模型计算的泥沙浓度与实测结果以及不分组模型 结果比较 (Southerndown 观测点) 时间/h 参考文献： 1 Yang C, Jiang C, Lin B. Modelling Graded Sediment Transport and Bed Evolutionin a Tidal Harbor J. Journal of Coastal Research, 2013, 29(3): 736-744. 2 Einstein H A. The Bed-load Function for Sediment Transport in Open Channel Flows R. Washington, DC:Department of Agriculture, Technical Bulletin 1026， 1950. 3 Misri R L, Ranga Raju K G, Garde R J. Bed Loadtransport of Coarse Nonuniform Sediments J. Journal of Hydraulic Engineering, 1984,110(3): 312-328. 4 Samaga B R, Ranga Raju K G, Garde R J. Suspended Load Transport Rate of Sediment Mixture J. Journal of Hydraulic Engineering, 1986, 11(11-13): 1019-1038. 5 James S C, Shrestha P L, Jesse D,et al. Modeling Noncohesive Sediment Transport Using Multiple Sedimentsize Classes J. Journal of Coastal Research, 2006,22(5): 1125-1132. 6 吴挺峰, 罗潋葱, 崔广柏, 等. 河流型水库垂向二维水沙数学模型J. 水科学进展, 2009, 20(2): 216-221. 7 Yang, C, Jiang C. A New Model for Predicting Bed Evolution in Estuarine Area and Its Application in Yellow River Delta J. Journal of Hydrodynamics, 2011, 23(4): 457-465. 8 杨晨. 地形演变模型与湿地健康评价方法及其在黄河口的应用D. 北京: 清华大学, 2011. 9 王光谦, 胡春宏.泥沙研究进展M. 北京: 中国水利水电出版社, 2006. 10 Wu Y, Falconer R A. Refined Two-dimensional Ultimatequickest Scheme for Conservative Solute Transport Modelling C/ Proceed- ings of Third International Conference on Hydro-Scienceand Engineering,Cottbus, Germany, the University of Mississippi, 1998, 1 (143): 1-13. 104 11 Lin B, Falconer R A. Modelling Sediment Fluxes Inestuarine Waters Using a Curvilinear Co-ordinate Grid System J. Estuarine, Coastal, and Shelf Science, 1995, 41: 413-428. 12 Gao G, Falconer R A, Lin B. Numerical Modelling Sediment-bacteria Interaction Processes in the Severn EstuaryJ. Journal of Water Resource and Protection, 2011, 3(1): 22-31. 13 Stapleton C M, Wyer M D, Kay D, et al. Fate and Transport of Particles in EstuariesR. Bristol, UK: Environment Agency, 2007. Study on the Fractional Modeling of Sediment Transport in Estuarine and Coastal Waters YANG Chen1, LU Kuan2 1. College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China; 2. National Ocean Technology Center, Tianjin 300112, China Abstract：A new two -dimensional fractional numerical model is established in this paper to simulate the hydrodynamics and graded sediment transport process in estuarine and coastal waters. The model is coupled by the hydrodynamic module, sediment transport module and bed evolution module. The Corlioli effect, bed shear and wind shear are introduced in the hydrodynamic module and the wetting-drying method is adopted to deal with the moving bo