剖面二维非恒定悬移质泥沙扩散方程的数值方法论文.docVIP

　　剖面二维非恒定悬移质泥沙扩散方程的数值方法论文 摘要：通过讨论剖面二维非恒定泥沙扩散方程的数值方法，建立了一种用于求解含沙量分布沿程变化的差分格式(Z-C格式)并通过一个具体的数值例子说明了计算的方法步骤。 关键词：扩散方程 差分格式 精度 稳定性 1 引言 数学模拟方法正在成为研究河流泥沙问题的重要手段。目前，一维数学模型发展较成熟，已广泛应用于模拟长河段的长期变形，但它只能给出河段平均冲淤深度的沿程变化，如需了解短河段的河床变形细节，则要采用二维以至三维数学模型。不论是一维数学模型还是平面二维维数学模型，都不能反映含沙量沿垂线的分布状况.freel (3) 其中εm为动量传递系数，β为修正值。由勃兰特尔掺长理论可得 εm=κu*y/(1-y/h) (4) 式中κ为卡门常数，u*为摩阻流速。 对于u，我们取卡曼-勃兰特尔对数流速分布公式 (umax-u)/u*=(1/κ)ln(h/y) (5) 令时,D1Sn+1/2j,M-1+D2Sn+1/2j,M+D3Sn+1/2j,M+1=C3LxSnj,M-Snj,M 令Hl=-Snj,l+C3LxSnj,l (1 l M) H1=-Snj,1+C3LxSnj,1-D1Sn+1/2j,0 HM=-Snj,M+C3LxSnj,M-D3Sn+1/2j,M+1 其中Sn+1/2j,0和Sn+1/2j,M+1由边界条件给出，则用矩阵形式表示为 （26） 第二步，再对(25)式用追赶法求第n+1层的值 令Fj=E1δ2ySn+1/2j,l+E2LySn+1/2j,l+Sn+1/2j,l (1 j N) F1=E1δ2ySn+1/21,l+E2LySn+1/21,l+Sn+1/21,l-C1Sn+10,l FN=E1δ2ySn+1/2N,l+E2LySn+1/2N,l+Sn+1/2N,l-C3Sn+1N+1,l 其中Sn+10,l和Sn+1N+1，l由边界条件给出，则同理可得矩阵方程 (27) 这样，按此步骤一层层地计算。 4.3 数值模拟合理性分析 受所掌握的实测资料的限制，目前尚无法对本文提出的算法与含沙量沿垂线分布的实测值进行对比。我们用库里·阿雷克沉沙池５的实测资料作了垂线平均值沿程变化的比较。该沉沙池的主要数据为：池深h=1.53m;平均流速u=0.12m/s;泥沙沉速ω=0.0176cm/s;悬浮指标Z=0.01。计算时取卡门常数κ=0.4，a=0.05h。表1给出了计算值和实测值，结果表明，计算值和实测值比较符合。 表1 断面平均含沙量验证(单位：kg/m3) Verifications of cross-sectional average sediment concentrations 距离(m) 0 200 400 600 800 1000 1200 计算值 3.012 2.573 2.071 1.639 1.209 0.849 0.539 实测值 3.012 2.650 1.810 1.520 1.141 0.822 0.561 为了进一步分析含沙量垂线分布计算结果的合理性，我们对另一组较粗的泥沙(ω=0.616cm/s,Z=0.4)进行了对比计算。图1/,图2分别为两组沙的计算结果。从图中可以看出，计算结果符合含沙量沿垂线分布的一般规律，粗沙分布不均匀，细沙分布较均匀；近底浓度相对较大，水面浓度相对较小，不存在Rouse公式中水面含沙量为0的缺陷；含沙量沿程衰减的特性较为明显。图3为较粗一组泥沙的相对含沙量沿垂线分布的沿程变化情况。图3表明，尽管进口断面按Rouse公式给出了含沙量沿垂线的分布，但由于该断面实际处于不平衡输沙状态，这种分布并不是稳定的。在距进口200m处，泥沙的分布调整到一种不平衡输沙状态，随着泥沙的沿程淤积，水流输沙向平衡方向发展，垂线平均含沙量趋向于水流挟沙力，而含沙量沿垂线分布向平衡时的分布状态(Rouse公式)发展。由于这种发展是趋向于稳定状态，因此愈接近下游，分布愈靠近Rouse公式。计算结果表明，本文提出的计算方法是合理可行的。 图1 含沙量垂线分布的沿程变化(Z=0.01) 图2 含沙量垂线分布的沿程变化(Z=0.4) Changes of vertical distributions of sediment concentrations(Z=0.01) Changes of vertical distributions of sediment concentrations(Z=0.4) 5 结语 本文建立了求解二维非恒定泥沙扩散方程的一种差分格式(Z-C格式)。这种格式具有如下特点： 1.精度较高(具有二阶精度)。 2.稳定性好(无条件稳定)。 3.计算较方便(每一时段利用两次追赶法即可)。 图3