计算水力学基础

1、计算水力学基础李占松 编著郑州大学水利与环境学院内 容 简 介本讲义是编者根据多年的教学实践，并参考微机计算水力学（杨景芳编著，大连理工大学出版社出版，1991年5月第1版）等类似教材，取其精华，编写而成的。目的是使读者掌握通过计算机解水力学问题的方法，为解决更复杂的实际工程问题打下牢固的计算基础。书中内容包括：数值计算基础，偏微分方程式的差分解法，有限单元法；用这些方法解有压管流、明渠流、闸孔出流、堰流、消能、地下水的渗流及平面势流等计算问题。讲义中的用FORTRAN77算法语言编写的计算程序，几乎包括了全部水力学的主要计算问题。另外，结合讲授对象的实际情况，也提供了用VB算法语言编写的计算

2、程序。VB程序编程人员的话为了更好地促进水利水电工程建筑专业的同学学好微机计算水力学这门学科，编程员借暑假休息的时间，利用我们专业目前所学的VB中的算法语言部分对水力学常见的计算题型编制成常用程序。希望大家能借此资料更好地学习微机计算水力学这门课程。本程序着重程序的可读性，不苛求程序的过分技巧。对水力学中常用的计算题型，用我们现在所学的VB语言编制而成。由于编程员能力有限，程序中缺点和错误在所难免，望老师和同学及时给予批评指正。VB程序编程人员：黄渝桂 曹命凯前 言-计算水力学的形成与发展计算水力学作为一门新学科，形成于20世纪60年代中期。水力学问题中有比较复杂的紊流、分离、气穴、水击等流动

3、现象，并存在各种界面形式，如自由水面、分层流、交界面等。由各种流动现象而建立的数学模型（由微分方程表示的定解问题），例如连续方程、动量方程等组成的控制微分方程组，多具有非线性和非恒定性，只有少数特定条件下的问题，可根据求解问题的特性对方程和边界条件作相应简化，而得到其解析解。因此长期以来，水力学的发展只得主要藉助于物理模型试验。随着电子计算机和现代计算技术的发展，数值计算已逐渐成为一个重要的研究手段，发展至今，已广泛应用与水利、航运、海洋、流体机械与流体工程等各种技术科学领域。计算水力学的特点是适应性强、应用面广。首先流动问题的控制方程一般是非线性的，自变量多，计算域的几何形状任意，边界条件复

4、杂，对这些无法求得解析解的问题，用数值解则能很好的满足工程需要；其次可利用计算机进行各种数值试验，例如，可选择不同的流动参数进行试验，可进行物理方程中各项的有效性和敏感性试验，以便进行各种近似处理等。它不受物理模型试验模型律的限制，比较省时省钱，有较多的灵活性。但数值计算一是依赖于基本方程的可靠性，且最终结果不能提供任何形式的解析表达式，只是有限个离散点上的数值解，并有一定的计算误差；二是它不像物理模型试验一开始就能给出流动现象并定性地描述，却往往需要由原体观测或物理实验提供某些流动参数，并对建立的数学模型验证；三是程序的编制及资料的收集、整理与正确利用，在很大程度上依赖于经验与技巧。所以计算

5、水力学有自己的原理方法和特点，数值计算与理论分析观测和试验相互联系、促进又不能相互代替，已成为目前解决复杂水流问题的主要手段之一，尤其是在研究流动过程物理机制时，更需要三者有机结合而互相取长补短。近三、四十年来，计算水力学有很大的发展，替代了经典水力学中的一些近似计算法和图解法。例如水面曲线计算；管网和渠系的过水或输沙（排污）能力的计算；有水轮机负荷改变时水力震荡系统的稳定性计算研究；流体机械过流部件的流道计算以及优化设计，还有洪水波、河口潮流计算，以及各种流动条件下，不同排放形式的污染物混合计算等。上世纪70年代中期已从针对个别工程问题建立的单一数学模型，开始建立对整个流域洪泛区已建或规划中

6、的水利水电工程进行系统模拟的系统模型。理论课题的研究中，对扩散问题、传热问题、边界层问题、漩涡运动、紊流等问题的研究也有了很大的发展，并已开始计算非恒定的三维紊流问题。由于离散的基本原理不同，计算水力学可分为两个分支：一是有限差分法，在此基础上发展的有有限分析法；二是有限单元法，在此基础上提出了边界元法和混合元法，另外还有迎风有限元法等。目 录第一章 数值计算基础第一节 非线性方程式的解法1 迭代法2 牛顿-拉普森切线法3 二等分法第二节 线性方程组的解法1 高斯-塞得尔迭代法2 高斯列主元消去法第三节 插值第四节 拟合第二章 偏微分方程式的差分数值解法第一节 分类、数值解法和差分格式1 分类

7、2 数值解法3 差分格式第二节 椭圆型偏微分方程的数值解法第三节 抛物型偏微分方程的数值解法第三章 有压管流第一节 管网的水力计算-哈迪克劳斯法第二节 管网的水力计算-有限单元法第三节 简单管道的水击计算特征线法第四节 短管的水力计算第四章 明渠流第一节 明渠非恒定流的水力计算特征线法第二节 棱柱形明渠恒定非均匀渐变流水面曲线计算第三节 非棱柱形明渠恒定非均匀渐变流水面曲线计算第四节 天然河道水面曲线计算第五章 堰流和消能第一节 宽顶堰流水力计算第二节 消力池水力计算第一章 数值计算基础第1节 非线性方程式的解法1简单迭代法一、基本原理图1.1 简单迭代法原理图如图1.2所示，由，选取初值，带

8、入该式得，一直进行下去，则有收敛判别式：（为高阶小量，收敛判别常数）；收敛条件：。二、计算步骤（1）将变为；（2）选取初值；（3）迭代计算：；（4）比较，：若，迭代结束，；若，返回（3）继续迭代计算。例1-1 已知梯形断面底宽，边坡系数，m3/s，试编写用迭代法求此渠道中正常水深的程序（计算允许误差）。解：若已知正常水深求底宽，则底宽的计算表达式为，其迭代过程，与正常水深的迭代过程类似。程序1.1：简单迭代法求解正常水深表1.1 变量说明变量名意义Q，b渠道流量，底宽m，n渠道边坡系数，粗糙系数i渠道底坡h水深Private Sub Command1_Click()Dim Q!, b!, h!

9、, h0!, m!, i!, n!, eps!Q = Val(InputBox(请输入渠道流量Q=)b = Val(InputBox(请输入渠道底宽b=)n = Val(InputBox(请输入渠道粗糙系数n=)m = Val(InputBox(请输入渠道边坡系数m=)i = Val(InputBox(请输入渠道底坡i=)eps = Val(InputBox(请输入精度e=)Doh = (n * Q / i 0.5) 0.6 * (b + 2 * h0 * (1 + m 2) 0.5) 0.4 / (b + m * h0)If Abs(h - h0) eps Then Exit Doh0 =

10、hLoopPrint 该渠道的正常水深为：; Format(h, 0.000); 米End Sub依次输入下列数据：Q =15，b = 8，n = 0.025，m =1.5，i = 0.0009，eps=0.005，输出结果为：“该渠道的正常水深为：1.265米”。2.牛顿-拉普森（Newton-Raphson）切线法计算溢流坝下游收缩断面水深一、基本原理图1.2 切线法原理图已知方程，求解该方程的根。 ，)()(1nnnnxfxfxx-=+ 收敛判别式：二、计算步骤（1）由函数求其导数；（2）选取初值；（3）迭代计算：；（4）比较：若，迭代结束，；若，返回（3）继续迭代计算。例1-2已知某溢

11、流坝上游断面对下游河底的比能，矩形断面河道上，为流速系数，试用牛顿-拉普森方法编写计算的程序。为收缩断面的水深。解：由 整理得： （*）令， （*） （*）将（*）代入（*）式，以及对（*）式求导可得，由牛顿-拉普森切线法，可得收缩断面水深的迭代公式为收敛判别式：。对于本题，也可以利用简单迭代法求解收缩断面的水深，其迭代公式为程序1.2：牛顿-拉普森切线法求解收缩断面水深表1.2 变量说明变量名意义 Eo坝上游断面的总比能q，x单宽流量，坝面流速系数A，常数B，常数C，常数linjieH引入的函数Hc下游收缩断面水深Public Function linjieH(x!, q!, hc!) 建立

12、函数Dim h1!, A!, B!, C!, fhc!, fh1c!g = 9.8Eo = Val(InputBox(请输入断面比能Eo)X = Val(InputBox(请输入流速系数x)q = Val(InputBox(请输入断面流量q)e = Val(InputBox(请输入计算精度e)A = 2 * g * x 2B = 2 * g * x 2 * EoC = q * qhc = 1!Do 迭代过程fhc = A * hc 3 - B * hc 2 + Cfh1c = 3 * A * hc 2 - 2 * B * hchc2 = hc1 - fhc / fh1cIf Abs(hc2 -

13、 hc1) e Then Exit Dohc = hc2LoopEnd FunctionPrivate Sub Command1_Click()Call linjieH(x!, q!, hc!) 调用函数Print Format(hc, 0.000)End Sub3二等分法求方程的根。选取区间a,b使，则其根在区间a,b内。12继续二等分。每二等分一次区间缩短为原来的一半。第k次二等分后区间长度为：收敛判别式：计算步骤：（1）选取；使；（2）令，求；（3）若，则；若，则；（4）若，迭代结束，（或或）； 若，返回（2）继续进行二等分。例1-3梯形渠道底宽，求。解： ， 选取。其中。程序：变量说明

14、表：程序中：意义：Q，e流量，允许误差Ac，BcHc临界水深hc1 ，hc2区间（a,b）俩端点水深 hc3区间（a,b）的中点值b,m渠底宽度，边坡系数Dim b!, m!, Q!, hc!, g!, e!, Ac!, Bc!, hc1!, x!, hc2!, f1!, f2!, f3!, hc3!, p!Public Function LinjieH(b!, m!, Q!, hc!)g = 9.8Ac! = (b + m * hc) * hcBc! = b + 2 * m * hcp = Ac 3 / Bc - Q 2 / gEnd FunctionPrivate Sub Command1

15、_Click()g = 9.8b! = Val(InputBox(请输入矩形断面的底宽b)m! = Val(InputBox(请输入边坡系数m)Q! = Val(InputBox(请输入流量Q)e! = Val(InputBox(请输入计算精度e)hc1! = 0.01x! = Q / bhc2! = (x 2 / g) (1 / 3)Dohc3 = (hc1 + hc2) / 2f1 = LinjieH(b!, m!, Q!, hc1!)f2 = LinjieH(b!, m!, Q!, hc2!)f3 = LinjieH(b!, m!, Q!, hc3!)If f1 * f3 0 Thenf

16、2 = f3: hc2 = hc3Elsef1 = f3: hc1 = hc3End IfLoop Until Abs(f3) eCall LinjieH(b!, m!, Q!, hc!)Print hc3End Sub第2节线性方程组的解法1高斯-赛德尔迭代法三元一次方程组：初值：雅可比迭代法：赛德尔迭代法：第步：雅可比迭代法：赛德尔迭代法：收敛判别式：对于计算步骤：1给定初值：2迭代计算：（1）雅可比迭代法：（2）赛德尔迭代法：3收敛判别式：赛德尔迭代法程序（4）变量说明表：程序中：意义：a(i, j)方程组系数矩阵元素b(i)方程组右端项元素x0(i) ， x1(i)方程组待求的未知量N

17、矩阵的维数Dim i As Integer, j As Integer, n As Integer, sum As Integer, s1 As Double, s2 As DoubleDim a() As Integer, b() As Integer, x0() As Double, x1() As Double, EPS As Double, t#, z#Public Function GaussSeidel()n = InputBox(输入矩阵的维数：)ReDim Preserve a(n, n), b(n), x0(n), x1(n)For i = 1 To nFor j = 1 T

18、o na(i, j) = InputBox(输入a(i,j)的值：, a( & i & , & j & )Print a(i, j);Next jPrintNext iFor i = 1 To nb(i) = InputBox(输入b(i)的值：, b( & i & )x0(i) = 0Print b(i);Next iPrint此处判断该系数矩阵是否为严格对角占优阵,若不想判断也可将下10行删除For i = 1 To nFor j = 1 To nt = Abs(a(i, i)z = z + Abs(a(i, j)Next jIf t z - Abs(a(i, i) ThenPrint 您

19、输入的方程组的Gauss-Seidel迭代式不收敛！;End IfIf t = z - Abs(a(i, i) Then GoTo FNext iF: DoFor i = 1 To nFor j = i + 1 To ns2 = s2 + a(i, j) * x0(j)Next jFor j = 1 To i - 1s1 = s1 + a(i, j) * x1(j)Next jx1(i) = (b(i) - s1 - s2) / a(i, i)EPS = Abs(x1(1) - x0(1)s1 = 0s2 = 0Next iFor i = 1 To nx0(i) = x1(i)Next isu

20、m = sum + 1If EPS = 0. Then Exit DoLoopPrint 迭代了; sum; 次Print x(1)=; Format(x0(1), 0.000); x(2)=; Format(x0(2), 0.000); x(3)=; Format(x0(3), 0.000);声明:此程序是系数矩阵满足严格对角占优阵的情况下,Gauss-Seidel迭代法收敛;对于满足不可约对角占优阵的Gauss-Seidel迭代法收敛情况下,有兴趣的读者可查看有关资料,自己动手试试 3x1+x2+6x3=3 x1+x2+2x3=3 -x2+2x2-3x3=1 x1=-10 x2=3 x3=

21、5End FunctionPrivate Sub Form_Click()Call GaussSeidelEnd Sub雅可比迭代法 程序Dim i As Integer, j As Integer, n As Integer, sum As Integer, s1 As Double, s2 As Double, t#, z#Dim a() As Integer, b() As Integer, x0() As Double, x1() As Double, EPS As DoublePublic Function Jacobi()n = InputBox(输入矩阵的维数：)ReDim Pr

22、eserve a(n, n), b(n), x0(n), x1(n)For i = 1 To nFor j = 1 To na(i, j) = InputBox(输入a(i,j)的值：, a( & i & , & j & )Print a(i, j);Next jPrintNext iFor i = 1 To nb(i) = InputBox(输入b(i)的值：, b( & i & )x0(i) = 0Print b(i);Next iPrintFor i = 1 To nFor j = 1 To nt = Abs(a(i, i)z = z + Abs(a(i, j)Next jIf t z

23、- Abs(a(i, i) ThenPrint 您输入的方程组的Jacobi迭代式不收敛！;End IfIf t = z - Abs(a(i, i) Then GoTo FNext iF: DoFor i = 1 To nFor j = i + 1 To ns2 = s2 + a(i, j) * x0(j)Next jFor j = 1 To i - 1s1 = s1 + a(i, j) * x0(j)Next jx1(i) = (b(i) - s1 - s2) / a(i, i)EPS = Abs(x1(1) - x0(1)s1 = 0s2 = 0Next iFor i = 1 To nx0

24、(i) = x1(i)Next isum = sum + 1If EPS maxV ThenmaxV = Abs(a(i, j) 将主对角线上的最大元素找出maxI = i row of max elementmaxJ = j column of max elementEnd IfNext jNext iIf ii maxI Then exchange rowsd1 = b(maxI)b(maxI) = b(ii)b(ii) = d1For j = ii To nd1 = a(maxI, j)a(maxI, j) = a(ii, j)a(ii, j) = d1Next jEnd IfIf ii

25、 maxJ Then exchange columnsi1 = index(maxJ)index(maxJ) = index(ii)index(ii) = i1For i = 1 To nd1 = a(i, maxJ)a(i, maxJ) = a(i, ii)a(i, ii) = d1Next iEnd IfIf a(ii, ii) = 0# ThenExit FunctionEnd IfFor i = ii + 1 To n eleminationd1 = a(i, ii) / a(ii, ii)For j = ii + 1 To na(i, j) = a(i, j) - a(ii, j)

26、* d1Next jb(i) = b(i) - b(ii) * d1Next iNext iix(n) = b(n) / a(n, n) back substitutionFor i = n - 1 To 1 Step -1d1 = 0#For j = i + 1 To nd1 = d1 + a(i, j) * x(j)Next jx(i) = (b(i) - d1) / a(i, i)Next iFor i = 1 To nb(index(i) = x(i)Next iSolve = TrueEnd Function第三节 插值例 已知消力坎的淹没系数与相对淹没度对应关系如下表：.0.700

27、.720.740.760.78.0.940.930.9130.900.885.hs为下游水位超过坎顶的高度；H0为坎上全水头。试编一分别用线性插值和抛物型插值的程序，计算相对淹没度时的淹没系数值。1线性插值已知，求，要求满足则 令2拉格朗日插值二次插值（抛物型插值），求。设满足条件： 3n次插值（拉格朗日插值），求 线性插值程序变量说明表：程序中：意义：x(i)已知节点坐标y(i)已知节点函数值m待求函数值的节点坐标Dim m!, f1!, f2!, L1x!, L2x!, Lx!Private Sub Command1_Click()Dim x(4) As Single, y(4) As S

28、inglex(0) = 0.7: x(1) = 0.72: x(2) = 0.74: x(3) = 0.76: x(4) = 0.78y(0) = 0.94: y(1) = 0.93: y(2) = 0.915: y(3) = 0.9: y(4) = 0.885m = Val(InputBox(请输入一个数m)i = 0Dof1 = m - x(i)f2 = m - x(i + 1)L1x = (m - x(i + 1) / (x(i) - x(i + 1)L2x = (m - x(i) / (x(i + 1) - x(i)Lx = L1x * y(i) + L2x * y(i + 1)i =

29、 i + 1Loop Until f1 * f2 0Print Format(Lx, 0.0000)End Sub第4节拟合堰流： 最小二乘法：偏差若要使取最小值，则 例：直角三角形薄壁堰，求时， 0.71 0.86 0.90 0.97 1.10 1.25 1.37 1.44 1.53 1.64 1.77 1.871.88 3.13 3.38 3.88 5.14 8.05 9.22 10.50 12.30 14.33 16.85 18.80最小二乘法线性拟合直角三角形薄璧堰流量的程序变量说明表：程序中：意义：H直角三角形薄壁堰的堰上水头t1(i)，t2(i)堰上水头，流量，一维数组sum1，s

30、um2公式中AA1，AA2，V中间变量Q欲求流量Private Sub Command1_Click()Dim H#H = Val(InputBox(请输入直角三角形薄壁堰的堰上水头)Call nehe(H#)End SubPublic Function nehe(H#)Dim t1(8) As Double, t2(8) As Double, x(8) As Double, y(8) As Double, sum1#, sum2#, p1#, p2#, AA1#, AA2#, V#, r1#, r2# t1(0) = 0.71: t1(1) = 0.86: t1(2) = 0.9: t1(3

31、) = 0.97: t1(4) = 1.1: t1(5) = 1.25: t1(6) = 1.37: t1(7) = 1.44: t1(8) = 1.53 t2(0) = 1.88: t2(1) = 3.13: t2(2) = 3.38: t2(3) = 3.88: t2(4) = 5.14: t2(5) = 8.05: t2(6) = 9.22: t2(7) = 10.5: t2(8) = 12.3 For i = 0 To 8 y(i) = Log(100 * t2(i) Next i For i = 0 To 8 x(i) = Log(10 * t1(i) Next i For i =

32、0 To 8 AA1 = x(i) AA2 = y(i) sum1 = sum1 + AA1 sum2 = sum2 + AA2 Next i V = sum2 - sum1 For i = 0 To 8 r1 = x(i) * y(i) * V * 9 r2 = x(i) 2 * V * 9 p1 = p1 + r1 p2 = p2 + r2 Next i m = p1 / p2 b = Abs(sum2 - m * sum1) / 9 Q = Exp(b + m * Log(H * 10) Print QEnd Function第二章 偏微分方程式的差分数值解法第1节 分类、数值解法和差分

33、格式1分类二阶线性偏微分方程式：典型例子：椭圆型偏微分方程式：拉普拉斯方程式 抛物型偏微分方程式：扩散方程式 双曲型偏微分方程式：波动方程式 2数值解法主要步骤：划分网格；构造差分方程式；求解差分方程式。3差分格式函数及在处的泰勒级数展开式分别为： （1）向前差分（一阶精度）（2）向后差分（一阶精度）（1）-（2）得：中心差分（二阶精度）（1）+（2）得：中心差分（二阶精度）第2节 椭圆型偏微分方程式的数值解法基本方程式：（在域A内）边界条件：（1）第一类边界条件：（2）第二类边界条件：（3）第三类边界条件：有限差分解法：（1）网格划分：（2）构造差分方程式： ， 五点格式：（3）解差分方程式

34、1高斯迭代法2高斯-赛德尔迭代法3超松弛迭代法-超松弛因子（=1时为高斯-赛德尔迭代法，1为亚松驰迭代，1为超松弛迭代）（计算试验表明：取1.21.5之间数值时迭代次数最小）判别收敛性： 例：如图所示为夹在两个无限平行壁间的圆柱绕流，无限远处的来流速度分布均匀，且，二平板间通过的流量为，圆柱半径，试用有限差分法编写求流函数在流场中分布的程序。解： （一）网格划分：，。（二）构造差分方程式： （在域A内）第一类边界条件：（1）；（2）；（3）；第二类边界条件：（4）。五点差分格式：超松弛迭代法：（=1.21.5）边界条件：（1）（ab边界）；（2）；（3）bc边界的处理： (a) (b) ；（4

35、）；（ae边界）；（5）cd边界的处理: 的计算：(a)直接转移：(b)直线内插：如上图所示，。则，其中a表示P1与网格线和圆柱壁面交点之间的距离，是（7，3）和（7，4）两网格点之间的距离即y方向的网格步长；b表示P2与j=2网格线和圆柱壁面交点之间的距离，是（6，2）和（5，2）两网格点之间的距离即x方向的网格步长。（三）解差分方程式：逐点迭代计算。每循环一次，进行收敛条件的判别。收敛判别条件：有限差分法计算两个无限平行壁间圆柱绕流流函数在流场中分布的程序变量说明表：程序中：公式中：意义：psi(i, j)节点的流函数值，二维数组w松弛因子dij，sdij0，sdij1Private Su

36、b Command1_Click()Dim psi(8, 5), w, sdij0, sdij1For j = 1 To 5 For i = 1 To 8 psi(i, j) = 0 Next iNext jFor j = 2 To 4 psi(1, j) = (j - 1) * 0.5Next jFor i = 2 To 8 psi(i, 5) = 2Next iDow = 1.35k = 0sdij0 = 0 k = k + 1sdij1 = 0psi(6, 2) = psi(5, 2) * 0.125 / 0.625psi(7, 3) = psi(7, 4) * 0.125 / 0.62

37、5psi(8, 4) = (2 * psi(7, 4) + 2) / 4For j = 2 To 4 For i = 2 To j + 3 dij = w * (psi(i + 1, j) + psi(i - 1, j) + psi(i, j + 1) + psi(i, j - 1) / 4 - psi(i, j) sdij1 = sdij1 + Abs(dij) psi(i, j) = psi(i, j) + dij Next iNext jIf Abs(sdij1 - sdij0) 0.00001 Then Exit Dosdij0 = sdij1LoopFor j = 1 To 5Pri

38、ntFor i = 1 To 8 Print Space(1), Format(psi(i, j), 0.000); Next iNext jEnd Sub第3节 抛物型偏微分方程式的数值解法例：两无限平板间油的流动，上平板水平运动，0.1s内u=0线性增加到，其后保持不变。求两平板间的速度分布随时间的变化。解：一基本方程式 边界条件： 初始条件：二有限差分解法（一）网格划分 取空间步长 其中；取时间步长 其中。（二）构造差分方程式（1） 显式格式：若令，则 稳定性条件：（2） 隐式格式：若令，则 隐式格式无条件稳定。初始条件：边界条件：（三）解差分方程式显格式差分法计算两平板间速度分布随时间

39、变化的程序变量说明表：程序中：公式中：意义：n，Hn，a/m时段数，距离步长MM距离a被等份的份数AA两平板间的距离V液体的运动粘滞系数U0平板的最终速度RR显格式的稳定系数u(i, j)点流速Tj中间变量DTt时间步长Private Sub Command1_Click()Dim m, n, A, v, H, U0, R, DT, u(50, 501), t(501), tj, count%m = 10: n = 500: A = 1: v = 1: H = A / m: U0 = 10: R = 0.1: count = 0DT = R * H 2 / vFor i = 1 To m +

40、1u(i, 1) = 0Next iFor j = 1 To n + 1u(1, j) = 0t(j) = (j - 1) * DTtj = t(j)If tj 0.1 Thenu(m + 1, j) = tj * U0 / 0.1Elseu(m + 1, j) = U0End IfNext jFor j = 1 To nFor i = 2 To mu(i, j + 1) = R * u(i + 1, j) + (1 - 2 * R) * u(i, j) + R * u(i - 1, j)Next iNext jFor j = n + 1 To 1 Step -50Print jFor i =

41、 1 To m + 1count = count + 1Print Space(1), Format(u(i, j), 0.00);If count = 11 Thencount = 0: PrintEnd IfNext iNext jEnd Sub第3章 有压管流第1节管网的水力计算-哈迪克劳斯法例1：如图所示环状管网，假设。试求：（1）各管段的流量分配；（2）各节点的测压管水头。管段12345678910管长600600300600300900300600300600管径0.40.40.30.30.20.20.30.30.30.3沿程阻力系数0.020.020.030.030.040.04

42、0.030.030.030.03解：（一）主要公式连续性方程：。n为该点的管段数目，流进节点，流出节点。能量方程式：。环的正方向与该管段流向相同时，否则 ， ， 所以：（二）（1）假设各管段的流量分配，使其满足连续性方程；（2）A环路水头残差；设的修正流量为，使；即：所以： ，定性分析其正确性：如果初分流量时只是数值误差，不是方向“误差”的话，分母总为正值。则若分子（即残差）为正，需减小。亦即正的水头损失太大了，需减小流量；而负的水头损失太小了，需增加流量。这是正确的。反之亦然。（3）计算步骤：1根据各管段的计算；2任意选择各管段中的，满足连续性方程（给定各管段的流动方向）；3选定各闭合环路的

43、正流动方向（如顺时针方向）来确定各管段的，并构成矩阵，k-环路编号（行），I-管段编号（列）。对于环路中不包含的管段，矩阵元素。当该管段流向与该环方向相同时，；当该管段流向与该环方向相反时，。矩阵：（第k闭合环路- ） ik12345678910 1-1+100000-10+1200+1+10000-1-13 0000+1-1-1+1+104.计算（对每个闭合环路），扩展为：。不在该环时，；在该环时。5计算各管段中的流量，扩展为：。对于非共用管，只有一个不为零；对于平面共用管，只有两个不为零，并且当均规定顺时针方向为环路正方向时，一个为“+1”，一个为“-1”。6判别计算误差（如的允许误差），

44、是否满足：（其中为总环数），若不满足，返回第4步继续进行各环路修正流量的计算。7计算各节点的测压管水头（需构造节点及管段号之间的关系数组）。几点解释：（1）只管数值，不管方向。流量方向事先已假定，最终只要不修正为负值，方向就与原假设方向相同。如若修正为负值，则与原假设方向相反。（2）连续性方程中流量方向用符号函数表示。流入为“+1”，流出为“-1”。（3）能量方程式中只表示修正流量的数值。为正值时，即为均按顺时针方向增加流量；为负值时，即为均按逆时针方向增加流量。（4）的意义为：；。哈迪-克劳斯法管网的水力计算的程序变量说明表：程序中：公式中：意义：d(i)，l(i)d,l各管段的管径，管长y

45、(i)各管段的沿程阻力系数R(i)R公式中的系数Q(i，j)Q第j次迭代第i个管段的流量M(i，j)表示管段中的流动方向是否与环路方向一致的符号矩阵H(i)H各管段的水头B, A中间变量c(i)，x(i)中间变量Private Sub Form_Click()Dim d(1 To 10), l(1 To 10) As Integer, y(1 To 10), R(1 To 10), Q(1 To 10, 100) As Double, w#, M(1 To 3, 1 To 10)Dim Q1(1 To 3) As Double, B, A, c(1 To 10), H(1 To 8), x(1 To 10) 输入已知数据d(1) = 0.4: d(2) = 0.4: d(3) = 0.3: d(4) = 0.3: d(5) = 0.2: d(6) = 0.2: d(7) = 0.3: d(8) = 0.3: d(9) = 0.3: d(10) = 0.3 l(1) = 600: l