泥沙参数选择方法及在悬沙浓度研究中的应用

董 佳，张宁川

（大连理工大学海岸和近海工程国家重点实验室，大连116024）

目前，在港口和海岸工程工可研阶段，论证工程对泥沙运动及海床演变影响多采用数值模拟方法［1-3］。该方法中需要确定诸多影响泥沙运移的基本参数，如泥沙沉降速度、泥沙扩散系数、控制侵蚀速度的比例值、软底床的侵蚀系数、硬底床的侵蚀指数、床面糙率、河床淤积物的干密度、床层间的泥沙转换率等。上述参数选择恰当与否将直接影响数值计算的可靠性。通常情况，研究者通过实测资料验证或经验获取得，然而当实测资料缺乏或经验不足时，往往难以得到符合实际的计算结果。

对于泥沙沉降速度可通过物理模型试验或较为成熟的经验公式计算得到较为可靠的数值；对于以潮流为主要环境动力的淤泥质海床而言，经过数值计算表明，泥沙扩散系数，控制侵蚀速度的比例值，软底床的侵蚀系数，硬底床的侵蚀指数等对悬沙浓度的影响较小。因而以潮流为主要动力的淤泥质海域悬沙浓度计算中，床面糙率、河床淤积物的干密度、床层间的泥沙转换率3个参数的最佳取值需要进行详细探讨。

首先建立二维水动力模型，在通过实测资料验证的基础上，进行悬沙浓度计算，定量地讨论了每个参数的选择对悬沙浓度的影响程度。在床面糙率、河床淤积物的干密度、床层间的泥沙转换率可变范围内，分别调整各个参数的取值，以计算得到的悬沙浓度变化与实测资料吻合为目标，搜寻得到上述3个参数的最优取值。在此基础上，对洋山港全水域中各特征点的悬沙浓度进行了模拟并和实测值吻合较好，验证了所选择参数的可靠性。

1 水流泥沙数学模型

1.1 控制方程

水流泥沙数学模型控制方程包括潮流控制方程和悬沙浓度控制方程。

对于潮流控制方程，采用包含二阶紊动项的平面二维非恒定流体运动方程

式中：S为悬沙浓度；Sm为泥沙的沉降和再悬浮浓度；Cb为底层泥沙浓度；τb为瞬时床底剪切应力；τcd为淤积临界剪切应力；τce为侵蚀临界剪切应力；ε为泥沙水平扩散系数；α为软底床的侵蚀系数；n为硬底床的侵蚀指数；E为控制侵蚀速度的比例因子。

1.2 定解条件

（1）初始条件：在计算的初始时刻给出潮位、悬沙浓度等初始值，其选取范围较为宽松，因为初始条件的误差可以在正确的边界条件控制下很快消失。

（2）潮流模型开边界条件：在海域开边界，海面水位的边界条件由边界处主要分潮的调和常数计算可得

式中：Em为该点相对于平均海平面的水位；ai、wi、φi分别是第i格分潮的振幅、频率和迟角。

（3）潮流模型固体侧边界条件：在研究区域的固体侧边界处，一般假定垂直于固体海岸的法向速度为零，即。

（4）悬沙浓度计算模型开边界条件：在入流时段S为给定值，在出流时段由下式计算得

（5）悬沙浓度计算模型固体边界条件：一般假定固边界处，法向泥沙通量为零，即。

1.3 数值计算方法与网格生成方法

数值计算采用中心有限体积法。考虑到非结构网格具有节点随意编号、网格随意加密以提高精度、可以对复杂的固岸边界和水下地形进行比较精确的模拟、少量的节点数可得到满意的精度等优点，故在数值计算中网格采用非结构网格。

1.4 水动力模型验证

以洋山港工程水域为水动力模型验证对象。洋山深水港位于上海市南汇嘴东南海域的崎岖列岛海域，该工程近十几年积累了丰富且可靠性较高水动力、泥沙现场同步观测资料。模型范围为东经121.8°～122.3°，北纬30.3°～30.8°，该范围计算网格共剖分为35 394个。

水动力模型验证包括计算水域潮位、流速及流向验证等内容。图1给出了潮位计算值与实测值的比较示例，图1结果表明，不同位置潮位实测值和计算值吻合良好。图2给出了流速、流向计算值与实测值的比较示例，图2结果表明，不同位置流速及流向的实测值和计算值的吻合程度很高。模拟计算结果较好地再现了洋山港区非正规半日潮、水流流动往复性较强的特点。

2 泥沙参数取值及对悬沙浓度的影响

2.1 悬沙浓度计算的主要影响参数筛选

如前所述，影响泥沙运动的参数包括泥沙沉降速度、泥沙扩散系数、控制侵蚀速度的比例值、软底床的侵蚀系数、硬底床的侵蚀指数、床面糙率、河床淤积物的干密度、床层间的泥沙转换率。

对泥沙沉降速度，理论上与泥沙粒度有关。对于泥沙粒径d<0.01 mm的泥沙，将发生絮凝现象，絮凝沉速在0.000 3～0.000 6 m/s，可直接取絮凝沉速为0.000 4 m/s；而泥沙粒径在0.03

表1 泥沙沉降速度Tab.1 Velocity of sediment

对于二维模型中泥沙水平扩散系数ε，一般是基于涡粘系数类比公式计算而得，即ε=λεt（εt为涡粘系数），涡粘系数比例因子在0.5～1.5。

在泥沙的沉降速度已确定的条件下，分别计算了λ为0.5、1.0、1.5时的悬沙浓度值，结果几乎没有什么变化，表明ε对悬沙浓度计算影响很小，数值计算时可取涡粘系数比例因子λ为1.0。

对于控制侵蚀速度的比例值E，硬泥层一般取0.000 1（kg·m-2）/s，软泥层可0.000 005～0.000 02（kg·m-2）/s取值，分别取软底床的为0.000 005（kg·m-2）/s、0.000 01（kg·m-2）/s、0.000 015（kg·m-2）/s、0.000 02（kg·m-2）/s 4 组值进行悬沙浓度的计算，结果表明悬沙浓度变化很小，即表明E对悬沙浓度计算的影响不大，可直接取值0.000 02（kg·m-2）/s。

对于硬底床的侵蚀指数n，一般取10；对于软底床的侵蚀系数α通常在4～20变化，对于淤积偏重的地区可选择下限值4，侵蚀偏重的地区取20，而对于冲淤相当的可在中间挑选一个值即可。

综上所述，泥沙沉降速度ω、泥沙扩散系数ε、控制侵蚀速度的比例值E、软底床的侵蚀系数α、硬底床的侵蚀指数n，在悬沙浓度的数值计算中分别按照上述方法取值，对结果的影响不大。由此筛选出对悬沙浓度计算影响最大的3个参数分别为：床面糙率、河床淤积物的干密度以及床层间的泥沙转换率。

下面分别讨论床面糙率、河床淤积物的干密度、床层间的泥沙转换率对悬沙浓度计算结果的影响。

2.2 床面糙率k

床面糙率是描述海床对水流阻力的参数，可采用谢才系数、曼宁系数、尼古拉兹糙率等表示。当采用尼古拉兹糙率时，床面糙率依赖于床面形状和底床泥沙颗粒的粒径。进行含沙量计算时，可将计算范围的局地床面形状均值化，即床面糙率仅依赖底床泥沙颗粒的粒径，其值可按n倍泥沙颗粒的粒径计算，即k=nd。

在扩散系数ε等其他5个参数固定不变、且按照2.1节所述取值的条件下，分别选取床面糙率k=1.0 d、2.5 d、4.0 d、5.5 d，分析床面糙率k的取值变化对悬沙浓度的影响。

根据经验，当d<0.5 mm时，有些学者直接取k=0.001 m。对于洋山港工程水域，底质平均中值粒为0.037 mm。因此，在讨论床面糙率对悬沙浓度的影响时，也将k=0.001 m作为一种处理方法，与k=nd的处理方法进行比较。图3给出了采用不同的床面糙率时悬沙浓度的计算结果。

由图3可见，对于床面糙率，若直接取0.001 m，计算得到的悬沙浓度随时间的变化过程与实测结果差别较大。若床面糙率按照粒径的倍数来计算，当糙率取不同的粒径倍数时，计算的悬沙浓度和实测浓度的误差也将随之变化（表2）。

将表2的不同粒径倍数情况下所计算得到的平均悬沙浓度和实测的平均悬沙浓度进行对比，将粒径倍数值与误差关系绘出曲线（图4）发现，当选取k≈2 d时，计算的平均悬沙浓度和实测浓度的误差最小，即床面糙率的最佳取值应为k≈2 d。

表2 床面糙率对平均悬沙浓度值的影响汇总Tab.2 Influence of bed roughness on value of suspended sediment concentration

2.3 河床淤积物的干密度

河床淤积物的干密度是一个重要的泥沙参数，对于冲淤体积和重量的换算、泥沙起动、河床冲刷及水流挟沙能力等均有相当重要的影响［6］，干密度是不随时间变化的。河床的最上层为最弱层，主要是浮泥或新淤积的泥沙，称之为软泥；其下面层的密度及强度都不断递增，称之为硬泥。

河床淤积物的干容量与孔隙率成反比，孔隙率与泥沙粒径成反比，河床淤积物的干容重随着时间的变化很小，根据大量的实验得出：软泥的干密度［7］会在100～400 kg/m3浮动，具体的数值是取决于各自的环境，和新淤积泥沙的湿容重也有密不可分的关系，而硬泥则是由底质粒径按照公式γd=1 750×d0.183［8］计算，其中粒径的单位为mm。

港区的海床分为两层：上层为软底床，即泥沙最新输运落淤的床层，分别计算了干密度为150 kg/m3、250 kg/m3、400 kg/m3这3种不同的密度量级（图5）。区域中新落淤的泥沙干密度取250 kg/m3时计算的泥沙浓度和实测值拟合较好，其他的误差在10%～60%不等；下层为硬底床，按照底层粒径和公式γd=1 750×d0.183计算干密度。

2.4 床层间的泥沙转换率

在模拟计算时，通常将海床分为可动层和固结层两层，层间的泥沙在动力作用下的泥沙交换量的多少采用泥沙转换率描述。该参数对于研究泥沙运动机理、考虑床沙和运动泥沙的关系、计算悬沙浓度及海床冲淤强度等均为最重要影响参数之一。其取值范围一般在0.1～0.001（kg·m-2）/s，该范围上、下限量值上相差100倍，各个取值对悬沙浓度计算结果的影响是显而易见的。

在此，与分析床面糙率k对悬沙浓度影响一样，在探讨泥沙转换率对悬沙浓度影响时，也将扩散系数ε等其他 5 个参数固定不变（按照 2.1 节所述取值），分别选取 m 为 0.000 1（kg·m-2）/s、0.001（kg·m-2）/s、0.01（kg·m-2）/s、0.1（kg·m-2）/s，计算了 4 组不同的转换率对应的悬沙浓度随时间的变化过程，与实测结果的比较汇总于图6。

由图 6 可知，当 m 取值为 0.01（kg·m-2）/s和 0.1（kg·m-2）/s时，即固结率较高时，由于大部分悬沙落於到床面后固结，使得悬沙浓度大幅度降低。换言之，悬沙浓度模拟计算时，如果过高估计了泥沙转换率，将使得计算悬沙浓度小于实际悬沙浓度。反之，计算悬沙浓度将大于实际悬沙浓度。计算表明，当泥沙转换率为0.001（kg·m-2）/s时，模拟的泥沙浓度和实测的泥沙浓度值拟合最好。即最佳床层间的泥沙转换率为0.001（kg·m-2）/s。

3 悬沙浓度计算在洋山港实际工程中的验证

对于影响悬沙浓度计算的各个泥沙参数，采用前节泥沙参数的取值方法，分别取泥沙沉降速度：絮凝沉速ω=0.000 4 m/s，非絮凝沉速按表1取值；泥沙扩散系数：涡粘系数比例因子λ=1.0；控制侵蚀速度的比例值：硬泥层 E=0.000 1（kg·m-2）/s，软泥层 E=0.000 02（kg·m-2）/s；硬底床的侵蚀指数：n=10；软底床的侵蚀系数：对于港区南部（淤积）α=4，对于港区北部（冲刷）α=20；床面糙率：k=2 d；泥沙干密度：对软底床γd=250 kg/m3；对硬底床γd=1 750×d0.183；床层间的泥沙转换率：m=0.001（kg·m-2）/s。

选定上述泥沙参数后，对洋山港全水域各特征代表点的悬沙浓度分别进行数值计算，并与实测结果进行对比。选取的代表点位置分布在洋山港外南北附近、港区、汊道等典型水域，代表点的具体位置见图7。图8给出了各代表点垂线平均悬沙浓度计算和实测结果比较示例。

图 8-a 为小洋山北部 N2 点（北纬°39′26.54″，东经 2°04′04.85″）大潮期间垂线平均含沙量时间变化过程实测结果和计算结果的比较示例；图 8-b～图 8-f依次为大洋山南部 S1 点（30°33′19.14″，122°04′00.86″）、三期码头中部 M11 点（30°36′39.82″，122°04′30.31″）、西口门中部 W2 点（30°37′32.10″，122°00′13.05″）、大山塘—大贴饼汊道 SC3 点（30°35′30.58″，122°02′25.58″）、以及小洋山一颗珠山汊道 K2 点（30°38′41.28″，122°02′23.38″）的大潮期间垂线平均含沙量时间变化过程计算结果和实测结果的比较示例。表3给出了上述各代表点大潮和小潮期间的时均垂线平均含沙量计算与实测值的比较。

表3 洋山港工程水域典型泥沙测点时均垂线平均悬沙浓度值的计算值与实测值比较Tab.3 Comparison between computed vertical suspended sediment concentration in the average time and measured results at representative point in the Yangshan Port

由图8和表3的比较可知，在洋山港周围水域，按照前述推荐的各参数取值，计算得到的含沙量浓度随时间变化过程及在一个潮周期内的时均值与实测结果的吻合程度是令人满意的。

4 结语

基于经过实测资料验证的二维水动力模型，在较准确模拟潮流动力的基础上，对于影响悬沙浓度计算有关的8个参数通过数值计算的方法：首先筛选出影响最大的3个参数为床面糙率、河床淤积物的干密度、床层间的泥沙转换率；然后在其各个参数取值的可变范围内，分别调整各个参数的取值，计算对比得到上述3个参数的最优取值。最后，就洋山港实际工程进行了悬沙浓度模拟计算，验证了所选择参数的可靠性。

研究表明，对于以潮流为主要环境动力的淤泥质海床而言，泥沙沉降速度可通过物理模型试验或较为成熟的经验公式计算得到较为可靠的数值；泥沙扩散系数，控制侵蚀速度的比例值，软底床的侵蚀系数，硬底床的侵蚀指数等对悬沙浓度的影响较小。对于床面糙率、软底床淤积物的干密度、床层间的泥沙转换率，最佳取值分别为 2 d、250 kg/m3、0.001（kg·m-2）/s。

［1］ZUO S H，ZHANG N C，LI B，et al.Numerical Simulations of Tidal Current&Sediment and Sea Bed Erosion and Deposition for Yangshan Deep-water Harbor of Shanghai［J］.International Journal of Sediment Research，2009，24（3）：287-298.

［2］左书华，张宁川，张征.岛群海域环境下泥沙运动及地形冲淤变化数值模拟研究［J］.泥沙研究，2011（2）：1-8.ZUO S H，ZHANG N C，ZHANG Z.Numerical modeling of sediment transport and seabed erosion and siltation in archipelago sea area［J］.Journal of Sediment Research，2011（2）：1-8.

［3］李文丹，李孟国，庞启秀.台山核电站取排水工程潮流泥沙数值模拟研究［J］.水道港口，2011，32（2）：94-100.LI W D，LI M G，PANG Q X.The tide and sediment numerical simulation study of Taishan nuclear power plant project［J］.Journal of Waterway and Harbor，2011，32（2）：94-100.

［4］XIN W J.Computational Techniques of 2D Tidal Flow in Estuaries and Bays［J］.China Ocean Engineering，1995，9（4）：395-404.

［5］曹祖德，杨树森，杨华.粉沙质海岸的界定及其泥沙运动特点［J］.水运工程，2003，352（5）：15-20.CAO Z D，YANG S S，YANG H.Definition of silt-sandy beach and its characteristics of sediment movement［J］.Port&Waterway Engineering，2003，352（5）：15-20.

［6］潘庆燊.三峡工程泥沙问题研究［M］.北京：中国水利水电出版社，1999.

［7］韩其为，王玉成，向熙珑.淤积物的初期干容重［J］.泥沙研究，1981（1）：1-6.HAN Q W，WANG Y C，XIANG X L.Initial specific weight of deposits［J］.Journal of Sediment Research，1981（1）：1-6.

［8］罗肇森.潮汐通道口拦门沙航道的淤积计算［J］.海洋工程，1992（3）：10-18.LUO Z S.Prediction of sedimentation for the navigation channel of a tidal inlet with mouth bay［J］.Ocean Engineering，1992（3）：10-18.