沙运动规律数学模型研究及清淤方案优选

1、1 基本理论和方法1.1 模型的基本方程 水流模型方程水流连续方程：(1)水流运动方程：(2)(3)其中：x、y为空间坐标，t为时间,h为水深，z为水位z=z0+h,z0为床面高程，u,v为流速在x、y方向上的分量，f为柯氏力系数f=2sin,为地球自转速度，为纬度，Cs为谢才系数，曼宁公式Cs=h1/6/n,n为糙率，g为重力加速度，为涡动粘滞系数，拉普拉斯算子。 泥沙模型方程    悬沙连续方程：(4)    河床冲淤演变过程：(5)其中：C为悬移质含沙量，Kx,Ky分别为x、y方向的泥沙扩散系数，P为水流作用下泥沙的起悬量，D为

2、水流作用下泥沙的沉降量，Cbx、Cby分别为x、y方向的河口浮泥层推移输沙率，z0为床面高程，Cm为床面极限含沙量。1.2 泥沙模型方程中物理量的确定 泥沙沉降量D的确定    泥沙沉降量D与悬移含量C和泥沙沉速*有关，可表示为：D=1C*(6)其中：1为校正值。 泥沙起悬量P的确定    泥沙起悬量P与饱和挟沙状态下的含沙量C*和泥沙沉速*有关，可表示为：P=2C*(7)    C*通过LaneKalinske公式确定为：C*=5.55I00.5(u*T/)exp(-(*/u*T)2)1.61(8)式中

3、u*T为潮流与波浪合成摩阻流速，u2*T=(u2+v2)/C2s+fwu20/4(9)u0=H/Tsh(2h/L)，波浪底部最大轨迹流速H波高，T波周期，L波长，fw波浪阻力系数。    ab=H/2sh(2h/L)-1，近底波浪质点振幅，r*=120d50床面粗度，d50床面泥沙中径；    ,泥沙的平均沉速，1，2分别为泥沙最大(d1)与最小(d2)颗粒的沉速；I0为床面活动泥沙含量。    这样，将床面层附近泥沙起悬量和沉降量区分开来处理的方法，具有明确的物理意义，可以清楚地反映出悬移泥沙与床面层内

4、运动泥沙及床面活动泥沙发生相互交换过程中矛盾的两个方面。在波浪掀沙和潮流输沙过程中，平衡只是相对的、暂时的，而不平衡则是绝对的、长期的。泥沙起悬量作为一个起悬函数(Pickup Function)的形式出现，有助于反映泥沙在床面层附近的交换机理，同时采用Lane-Kalinske公式，其优越之处在该公式在建立时，充分考虑了底部水流湍流脉动特性，反映了一定的力学机理，从目前流体力学研究来看，所不完善的是尚未区分大尺度拟序结构和小尺度随机脉动的作用，这方面我们正进行深入研究，但作为工程应用，该公式已经过了许多实验室和现场天然实测资料的验证，证明是可靠的1,2。总之，这样选取，有助于逐步弄清泥沙在床

5、面层附近的交换机理，也有助于从力学机理上反映掀沙和输沙过程的不平衡特性。1.3 模型方程的数值求解方法 水流方程的数值解    (1)杂交有限元方法    杂交有限元方法是三角形单元和等参四边形单元相交错结合的有限元方法，这种方法对于具有复杂边界和长宽尺寸悬殊的水域是非常合适的3,4。    设单元中的近似函数取：(10)其中i(x,y)为基函数，Z、U、V取相同的基函数；如果单元为三角形，i(x,y)为线性三角元基函数；如果单元是等参四边形，则i(x,y)为四结点等参基函数。  &#

6、160; 将单元插值近似函数(10)代入方程(1)、(2)、(3)的Galerkin积分式中，可得杂交有限元方程：(11)(12)(13)    式中A、B1、C1、D1、G、B2、C2、D2、B4、C4分别为系数矩阵，其元素分别为：Aij=(ij),，B3的元素与B1相同；C3的元素与C2相同。()代表在三角形单元或四边形单元内积分。对于不同的区域三角形单元或等参四边形单元；i,j,k分别为1，2，3或1，2，3，4.    (2)质量集中方法    在过行长系列时间过程量模拟时，为节省计算时间和计算机

7、内存，提高计算效率，按以下三种单元组合情况(图1)，分别采用集中质量法对单元系数矩阵进行化简，具体方法可参见文献3。(a)三角形单元(b)等参四边形单元(c)三角形-等参四边形单元图1    (3)时间方向的差分离散和半隐格式    为计算稳定，采用半隐差分格式，先利用Un、Vn和Zn通过离散式(12)和(13)，计算出Un+1和Vn+1值：AUn+1=(A-tB1-tG)Un-(C1·t-tfA)Vn-tD1Zn(14)AVn+1=(-tfA-tB2)Un-(C2·t-A+tG)Vn-tD2Zn(15)然后通过离

8、散式(11)，利用Un+1、Vn+1和Zn，求得Zn+1值：AZn+1=AZn-(tB3+tC3)Zn-(tB4Un+1+tC4Vn+1)(16) 泥沙方程的数值解 泥沙扩散方程和河床演变方程采用集中质量迎风有限元方法，具体如下：   (1)泥沙浓度场的数值模拟    为保证含沙度场的计算精度，同时考虑涨落潮流的实际背景，采用集中质量迎风有限元素方法。下面主要以三角形单元为例加以说明，等参四边形方法与之类似。    取插值函数为：(17)其中i(x,y)为基函数；如果单元为三角形，i(x,y)为线性三角元基函数；如

9、果单元是等参四边形，则i(x,y)为四结点等参基函数。    采用Garlerkin方法，方程(4)可写为：(18)    然后在计算节点Pi的相关单元中寻找其迎风单元，并将dcdt在迎风单元上离散，设(x1,y1),(xi,yi),(xi,xi)为节点Pi迎风单元ei的三个顶点坐标，i为迎风单元ei的面积。    其中：，这样，由式(18)可得：(19)    Pik(k=1,2,ni)为Pi的相关节点，其中：aiik=aii·iki+aii·iki，这

10、里，E代表Pi节点的相关单元点数。    (2)海床演变的数值模拟    利用向前差分格式离散床面随时间的变化项，然后代入方程(6)，整理得：(20)2 独流减河口潮流泥沙数学模型的研究2.1模型域的确定及剖分 河口海域模型域：以工农兵闸及以远10m等深线为东西相界，距离15km；南北以工农兵闸为起点，向北3km，向南4km，共7km，105km2的区域，河口区域网格剖分考虑清淤及导堤方案的需求，采用嵌套模型，即：在河口5km以内，工农兵闸附近，采取加密模型网格；该小模型网格为200m×100m,200m×50m等

11、网格，而在河口距工农兵闸5km以远，采用稀疏的大网格，为1000m×1000m的大模型网格，在大网格与加密小网格之间，采去变网格，使之大小相互匹配，从而使小模型形成统一的整体(参见图2).图2 独流减河河口区网格剖分2.2 河口水流数学模型的调试验证 河口模型采用1997年6月13日14日，17日18日，21日22日全潮观测资料，以6月13日14日中潮资料作为调试资料，以大潮17日18日，小潮21日22日作为验证资料，河口地形取1996年10月地形，各测点平面分布如图2，共8个全潮观测点。以图3 水位过程验证曲线(3号点)图4 流速模过程验证曲线(3号点)  &#

12、160; 从以上水位、流速过程线计算及实测图可以看出，水位计算值与实测值非常吻合；流速模吻合也较好；同时涨落潮流速椭圆图5也基本上符合实际情况，因此计算结果正确，潮流模型是可靠的，可进一步用于规划工程方案的预报。2.3 泥沙模型的调试与验证 采用非平衡泥沙数学模型，结合97年6月全潮观测的泥沙过程资料，对模型进行调试与验证。模型验证主要通过测点含沙量过程和床面冲淤变化情况来确定。    (1)含沙量验证：由于篇幅所限，仅给出一组点的验证结果(图5).    (2)床面冲淤变化情况：在含沙量验证正确的基础上，对床面演变情况进行了验证计算

13、，采用6月1314日中潮作为典型潮型，并考虑波浪作用，取多年平均波高，以两个月为一计时段，即每计算两个潮段，一天的冲淤量，使其乘以60，与原来地形叠加，再进行计算，如此累积，从96年10月计算到98年3月，其中考虑97年3月挖泥情况，2km内泥沙回淤计算结果见表1，床面地形变化实测与计算结果见图6.2.4 河口泥沙回淤计算结果分析 从表1可以看出，96年10月挖泥后，头两个月回淤量17万m3；之后逐月减少；而再次挖泥后，在97年5月至6月之间，回淤量又有所增加，以后又呈逐减趋势，直至最后达到平衡，但泥沙总的回淤量随时间的增长亦仍呈增长趋势，只是增长幅度有所减缓。表1 回淤量及总回淤量(万m3)

14、时间回淤量/两个月累积回淤量1996.1096.1117.317.31996.1297.116.033.31997.297.314.547.8挖泥清淤-407.81997.597.616.6924.491997.797.812.7937.281997.997.108.3245.501997.1197.125.5851.081998.198.24.4255.42图5 含沙量过程验证曲线(3号点)    利用1998年3月与1996年10月地形对照，2km范围内，泥沙回淤总量为63.9万m3，计算值(55.42万m3)与实测值(63.9万m3)较为接近，泥沙冲淤分布计

15、算与实测分布基本一致，含沙量和床面演变计算表明泥沙数学模型计算结果是正确的，所建模型是可靠的，可用于其它问题的研究和工程方案的预报计算。图6(a) 独流减河口98年3月实测地形(单位：m)图6(b) 独流减河口98年3月计算地形(单位：m)3 河口清淤方案泄流能力回淤量的计算及清淤方案的优选3.1 清淤方案泄流能力及回淤量计算 根据历年清淤情况，经海委规划设计科研处筛选，提出共12个清淤方案(见表2)，然后针对不同的清淤方案，分别进行河口泄流能力和泥沙回淤能力计算，从中筛选推荐清淤方案。利用河口大模型嵌套加密小模型进行回淤计算，以2km内各清淤方案的回淤总量作比较，优选最佳清淤方案。 

16、;   参照海河口大、中、小潮出现的机率为19.5%，48.5%，32%，通过各种潮型的搭配，计算泥沙回淤，最终得出各清淤方案下河口的泥沙年回淤量。为明显起见，将各清淤方案下河口的清淤量、年回淤量、清淤的工程效益(年回淤量与清淤量之比值)、河口日均最大泄流量列于表2.表2 清淤方案总表方案部位(m)清宽(m)清深(m)挖方量(万m3)年回淤量(万m3)日均泄流量(m3/s)年回淤量与清淤量之比10+000-0+260200-2.620+260-0+520200-150-2.624038.3365795.8%0+520-0+840150-2.6220+000-0+260200

17、-2.00+260-0+520200-150-2.035.430.9338088.2%0+520-0+840150-2.020+000-0+2602000+260-0+520200-150-3.044.548.23784108%0+520-0+84015040+000-1+1002000+100+0+200200-100-2.531.140.373280130%0+200-1+00010050+000-0+1002000+100-0+200200-100-2.536.243.83450120.9%0+200-1+20010060+000-0+1002000+100-0+200200-100-2

18、.540.3745.93560113.6%0+200-1+40010070+00-0+1002000+100-0+200200-80-2.525.3638.83225152.9%0+200-1+0008080+00-0+1002000+100-0+200200-60-2.519.635.03190178.6%0+200-1+0006090+00-0+1001500+100-0+200150-100-2.416.030.03100187.5%0+200-0+700100100+000-0+260200-2.00+260-0+520200-150-2.041.435.0361584.5%0+520

19、-1+000150-2.0110+000-0+1002000+100-0+200200-100-2.542.3640.5376395.6%0+200-1+500100120+00-0+1002000+100-0+200200-80-2.534.8938.03698109.5%0+200-1+500803.2 清淤方案优选    (1)比较相同清淤长度、相同清淤宽度、不同清淤深度的清淤方案(1、2、3)，可以看出，清淤深度越深(方案3)，河口日均泄量也越大，回淤量也越大。而清淤深度越浅的方案(方案2)，回淤量也越小。在这种清淤方式下，清淤深度越小，清淤量较小，回淤率也较小，相应地在满足泄流要求时