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

1基本理论和方法1.1模型的基本方程1.1.1水流模型方程水流连续方程：(1)水流运动方程：(2)(3)其中：x、y为空间坐标，t为时间,h为水深，z为水位z=z0+h,z0为床面高程，u,v为流速在x、y方向上的分量，f为柯氏力系数f=2ωsinφ,ω为地球自转速度，φ为纬度，Cs为谢才系数，曼宁公式Cs=h1/6/n,n为糙率，g为重力加速度，λ为涡动粘滞系数，，拉普拉斯算子。1.1.2泥沙模型方程悬沙连续方程：(4)河床冲淤演变过程：(5)其中：C为悬移质含沙量，Kx,Ky分别为x、y方向的泥沙扩散系数，P为水流作用下泥沙的起悬量，D为水流作用下泥沙的沉降量，Cbx、Cby分别为x、y方向的河口浮泥层推移输沙率，z0为床面高程，Cm为床面极限含沙量。1.2泥沙模型方程中物理量的确定1.2.1泥沙沉降量D的确定泥沙沉降量D与悬移含量C和泥沙沉速ω*有关，可表示为：D=Φ1Cω*(6)其中：Φ1为校正值。1.2.2泥沙起悬量P的确定泥沙起悬量P与饱和挟沙状态下的含沙量C*和泥沙沉速ω*有关，可表示为：P=Φ2C*ω*(7)C*通过Lane—Kalinske公式确定为：C*=5.55I0[0.5(u*T/ω)exp(-(ω*/u*T)2)]1.61(8)式中u*T为潮流与波浪合成摩阻流速，u2*T=(u2+v2)/C2s+fwu20/4(9)u0=πH/[Tsh(2πh/L)]，波浪底部最大轨迹流速H—波高，T—波周期，L—波长，fw—波浪阻力系数。ab=H/2[sh(2πh/L)]-1，近底波浪质点振幅，r*=120d50床面粗度，d50床面泥沙中径；,泥沙的平均沉速，ω1，ω2分别为泥沙最大(d1)与最小(d2)颗粒的沉速；I0为床面活动泥沙含量。这样，将床面层附近泥沙起悬量和沉降量区分开来处理的方法，具有明确的物理意义，可以清楚地反映出悬移泥沙与床面层内运动泥沙及床面活动泥沙发生相互交换过程中矛盾的两个方面。在波浪掀沙和潮流输沙过程中，平衡只是相对的、暂时的，而不平衡则是绝对的、长期的。泥沙起悬量作为一个起悬函数(PickupFunction)的形式出现，有助于反映泥沙在床面层附近的交换机理，同时采用Lane-Kalinske公式，其优越之处在该公式在建立时，充分考虑了底部水流湍流脉动特性，反映了一定的力学机理，从目前流体力学研究来看，所不完善的是尚未区分大尺度拟序结构和小尺度随机脉动的作用，这方面我们正进行深入研究，但作为工程应用，该公式已经过了许多实验室和现场天然实测资料的验证，证明是可靠的[1,2]。总之，这样选取，有助于逐步弄清泥沙在床面层附近的交换机理，也有助于从力学机理上反映掀沙和输沙过程的不平衡特性。1.3模型方程的数值求解方法1.3.1水流方程的数值解(1)杂交有限元方法杂交有限元方法是三角形单元和等参四边形单元相交错结合的有限元方法，这种方法对于具有复杂边界和长宽尺寸悬殊的水域是非常合适的[3,4]。设单元中的近似函数取：(10)其中Φi(x,y)为基函数，Z、U、V取相同的基函数；如果单元为三角形，Φi(x,y)为线性三角元基函数；如果单元是等参四边形，则Φi(x,y)为四结点等参基函数。将单元插值近似函数(10)代入方程(1)、(2)、(3)的Galerkin积分式中，可得杂交有限元方程：(11)(12)(13)式中A、B1、C1、D1、G、B2、C2、D2、B4、C4分别为系数矩阵，其元素分别为：Aij=(ΦiΦj),，，，，，，，，，B3的元素与B1相同；C3的元素与C2相同。()代表在三角形单元或四边形单元内积分。对于不同的区域——三角形单元或等参四边形单元；i,j,k分别为1，2，3或1，2，3，4.(2)质量集中方法在过行长系列时间过程量模拟时，为节省计算时间和计算机内存，提高计算效率，按以下三种单元组合情况(图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)然后通过离散式(11)，利用Un+1、Vn+1和Zn，求得Zn+1值：AZn+1=AZn-(ΔtB3+ΔtC3)Zn-(ΔtB4Un+1+ΔtC4Vn+1)(16)1.3.2泥沙方程的数值解泥沙扩散方程和河床演变方程采用集中质量迎风有限元方法，具体如下：(1)泥沙浓度场的数值模拟为保证含沙度场的计算精度，同时考虑涨落潮流的实际背景，采用集中质量迎风有限元素方法。下面主要以三角形单元为例加以说明，等参四边形方法与之类似。取插值函数为：(17)其中Φi(x,y)为基函数；如果单元为三角形，Φi(x,y)为线性三角元基函数；如果单元是等参四边形，则Φ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γ，这里，E代表Pi节点的相关单元点数。(2)海床演变的数值模拟利用向前差分格式离散床面随时间的变化项，然后代入方程(6)，整理得：(20)2独流减河口潮流泥沙数学模型的研究2.1模型域的确定及剖分河口海域模型域：以工农兵闸及以远10m等深线为东西相界，距离15km；南北以工农兵闸为起点，向北3km，向南4km，共7km，105km2的区域，河口区域网格剖分考虑清淤及导堤方案的需求，采用嵌套模型，即：在河口5km以内，工农兵闸附近，采取加密模型网格；该小模型网格为200m×100m,200m×50m等网格，而在河口距工农兵闸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个全潮观测点。以2.4.6.7.8测点作为边界点，以内点1.3.5作为验证点，由于篇幅所限，仅给出一组点的大潮计算验证结果(图3、图4)图3水位过程验证曲线(3号点)图4流速模过程验证曲线(3号点)从以上水位、流速过程线计算及实测图可以看出，水位计算值与实测值非常吻合；流速模吻合也较好；同时涨落潮流速椭圆图[5]也基本上符合实际情况，因此计算结果正确，潮流模型是可靠的，可进一步用于规划工程方案的预报。2.3泥沙模型的调试与验证采用非平衡泥沙数学模型，结合97年6月全潮观测的泥沙过程资料，对模型进行调试与验证。模型验证主要通过测点含沙量过程和床面冲淤变化情况来确定。(1)含沙量验证：由于篇幅所限，仅给出一组点的验证结果(图5).(2)床面冲淤变化情况：在含沙量验证正确的基础上，对床面演变情况进行了验证计算，采用6月13～14日中潮作为典型潮型，并考虑波浪作用，取多年平均波高，以两个月为一计时段，即每计算两个潮段，一天的冲淤量，使其乘以60，与原来地形叠加，再进行计算，如此累积，从96年10月计算到98年3月，其中考虑97年3月挖泥情况，2km内泥沙回淤计算结果见表1，床面地形变化实测与计算结果见图6.2.4河口泥沙回淤计算结果分析从表1可以看出，96年10月挖泥后，头两个月回淤量17万m3；之后逐月减少；而再次挖泥后，在97年5月至6月之间，回淤量又有所增加，以后又呈逐减趋势，直至最后达到平衡，但泥沙总的回淤量随时间的增长亦仍呈增长趋势，只是增长幅度有所减缓。表1回淤量及总回淤量(万m3)时间回淤量/两个月累积回淤量1996.10～96.1117.317.31996.12～97.116.033.31997.2～97.314.547.8挖泥清淤-407.81997.5～97.616.6924.491997.7～97.812.7937.281997.9～97.108.3245.501997.11～97.125.5851.081998.1～98.24.4255.42图5含沙量过程验证曲线(3号点)利用1998年3月与1996年10月地形对照，2km范围内，泥沙回淤总量为63.9万m3，计算值(55.42万m3)与实测值(63.9万m3)较为接近，泥沙冲淤分布计算与实测分布基本一致，含沙量和床面演变计算表明泥沙数学模型计算结果是正确的，所建模型是可靠的，可用于其它问题的研究和工程方案的预报计算。图6(a)独流减河口98年3月实测地形(单位：m)图6(b)独流减河口98年3月计算地形(单位：m)3河口清淤方案泄流能力回淤量的计算及清淤方案的优选3.1清淤方案泄流能力及回淤量计算根据历年清淤情况，经海委规划设计科研处筛选，提出共12个清淤方案(见表2)，然后针对不同的清淤方案，分别进行河口泄流能力和泥沙回淤能力计算，从中筛选推荐清淤方案。利用河口大模型嵌套加密小模型进行回淤计算，以2km内各清淤方案的回淤总量作比较，优选最佳清淤方案。参照海河口大、中、小潮出现的机率为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-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.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-1+000150-2.0110+000-0+1002000+100-0+200200-100-2.542.3640.