计算水力学基础

--本页仅作为文档封面，使用时请直接删除即可--

--内页可以根据需求调整合适字体及大小本页仅作为文档封面，使用时请直接删除即可--

--内页可以根据需求调整合适字体及大小--计算水力学基础(总79页)PAGE计算水力学基础李占松编着郑州大学水利与环境学院

内容简介本讲义是编者根据多年的教学实践，并参考《微机计算水力学》（杨景芳编着，大连理工大学出版社出版，1991年5月第1版）等类似教材，取其精华，编写而成的。目的是使读者掌握通过计算机解水力学问题的方法，为解决更复杂的实际工程问题打下牢固的计算基础。书中内容包括：数值计算基础，偏微分方程式的差分解法，有限单元法；用这些方法解有压管流、明渠流、闸孔出流、堰流、消能、地下水的渗流及平面势流等计算问题。讲义中的用FORTRAN77算法语言编写的计算程序，几乎包括了全部水力学的主要计算问题。另外，结合讲授对象的实际情况，也提供了用VB算法语言编写的计算程序。VB程序编程人员的话为了更好地促进水利水电工程建筑专业的同学学好《微机计算水力学》这门学科，编程员借暑假休息的时间，利用我们专业目前所学的VB中的算法语言部分对水力学常见的计算题型编制成常用程序。希望大家能借此资料更好地学习《微机计算水力学》这门课程。本程序着重程序的可读性，不苛求程序的过分技巧。对水力学中常用的计算题型，用我们现在所学的VB语言编制而成。由于编程员能力有限，程序中缺点和错误在所难免，望老师和同学及时给予批评指正。VB程序编程人员：黄渝桂曹命凯

前言计算水力学的形成与发展计算水力学作为一门新学科，形成于20世纪60年代中期。水力学问题中有比较复杂的紊流、分离、气穴、水击等流动现象，并存在各种界面形式，如自由水面、分层流、交界面等。由各种流动现象而建立的数学模型（由微分方程表示的定解问题），例如连续方程、动量方程等组成的控制微分方程组，多具有非线性和非恒定性，只有少数特定条件下的问题，可根据求解问题的特性对方程和边界条件作相应简化，而得到其解析解。因此长期以来，水力学的发展只得主要藉助于物理模型试验。随着电子计算机和现代计算技术的发展，数值计算已逐渐成为一个重要的研究手段，发展至今，已广泛应用与水利、航运、海洋、流体机械与流体工程等各种技术科学领域。计算水力学的特点是适应性强、应用面广。首先流动问题的控制方程一般是非线性的，自变量多，计算域的几何形状任意，边界条件复杂，对这些无法求得解析解的问题，用数值解则能很好的满足工程需要；其次可利用计算机进行各种数值试验，例如，可选择不同的流动参数进行试验，可进行物理方程中各项的有效性和敏感性试验，以便进行各种近似处理等。它不受物理模型试验模型律的限制，比较省时省钱，有较多的灵活性。但数值计算一是依赖于基本方程的可靠性，且最终结果不能提供任何形式的解析表达式，只是有限个离散点上的数值解，并有一定的计算误差；二是它不像物理模型试验一开始就能给出流动现象并定性地描述，却往往需要由原体观测或物理实验提供某些流动参数，并对建立的数学模型验证；三是程序的编制及资料的收集、整理与正确利用，在很大程度上依赖于经验与技巧。所以计算水力学有自己的原理方法和特点，数值计算与理论分析观测和试验相互联系、促进又不能相互代替，已成为目前解决复杂水流问题的主要手段之一，尤其是在研究流动过程物理机制时，更需要三者有机结合而互相取长补短。近三、四十年来，计算水力学有很大的发展，替代了经典水力学中的一些近似计算法和图解法。例如水面曲线计算；管网和渠系的过水或输沙（排污）能力的计算；有水轮机负荷改变时水力震荡系统的稳定性计算研究；流体机械过流部件的流道计算以及优化设计，还有洪水波、河口潮流计算，以及各种流动条件下，不同排放形式的污染物混合计算等。上世纪70年代中期已从针对个别工程问题建立的单一数学模型，开始建立对整个流域洪泛区已建或规划中的水利水电工程进行系统模拟的系统模型。理论课题的研究中，对扩散问题、传热问题、边界层问题、漩涡运动、紊流等问题的研究也有了很大的发展，并已开始计算非恒定的三维紊流问题。由于离散的基本原理不同，计算水力学可分为两个分支：一是有限差分法，在此基础上发展的有有限分析法；二是有限单元法，在此基础上提出了边界元法和混合元法，另外还有迎风有限元法等。

目录数值计算基础非线性方程式的解法迭代法牛顿-拉普森切线法二等分法线性方程组的解法高斯-塞得尔迭代法高斯列主元消去法插值拟合偏微分方程式的差分数值解法分类、数值解法和差分格式分类数值解法差分格式椭圆型偏微分方程的数值解法抛物型偏微分方程的数值解法有压管流管网的水力计算哈迪—克劳斯法管网的水力计算有限单元法简单管道的水击计算—特征线法短管的水力计算明渠流明渠非恒定流的水力计算—特征线法棱柱形明渠恒定非均匀渐变流水面曲线计算非棱柱形明渠恒定非均匀渐变流水面曲线计算天然河道水面曲线计算堰流和消能宽顶堰流水力计算消力池水力计算

数值计算基础第1节非线性方程式的解法1． 简单迭代法一、基本原理图简单迭代法原理图如图所示，由，选取初值，带入该式得，一直进行下去，则有收敛判别式：（为高阶小量，收敛判别常数）；收敛条件：。二、计算步骤（1）将变为；（2）选取初值；（3）迭代计算：；（4）比较，： 若，迭代结束，； 若，，返回（3）继续迭代计算。例1-1已知梯形断面底宽，边坡系数，，，m3/s，试编写用迭代法求此渠道中正常水深的程序（计算允许误差）。解：若已知正常水深求底宽，则底宽的计算表达式为，其迭代过程，与正常水深的迭代过程类似。程序：简单迭代法求解正常水深表变量说明变量名意义Q，b渠道流量，底宽m，n渠道边坡系数，粗糙系数i渠道底坡h水深PrivateSubCommand1_Click()DimQ!,b!,h!,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^^*(b+2*h0*(1+m^2)^^/(b+m*h0)IfAbs(h-h0)<epsThenExitDoh0=hLoopPrint"该渠道的正常水深为：";Format(h,"");"米"EndSub依次输入下列数据：Q=15，b=8，n=，m=，i=，eps=，输出结果为：“该渠道的正常水深为：米”。2.牛顿-拉普森（Newton-Raphson）切线法计算溢流坝下游收缩断面水深一、基本原理图切线法原理图已知方程，求解该方程的根。，，……，)())()(xfxfxx???收敛判别式：二、计算步骤（1）由函数求其导数；（2）选取初值；（3）迭代计算：；（4）比较：若，迭代结束，；若，，返回（3）继续迭代计算。例1-2．已知某溢流坝上游断面对下游河底的比能，矩形断面河道上，为流速系数，试用牛顿-拉普森方法编写计算的程序。为收缩断面的水深。解：由整理得：（*）令，，（**）（***）将（**）代入（*）式，以及对（***）式求导可得，由牛顿-拉普森切线法，可得收缩断面水深的迭代公式为收敛判别式：。对于本题，也可以利用简单迭代法求解收缩断面的水深，其迭代公式为程序：牛顿-拉普森切线法求解收缩断面水深表变量说明变量名意义Eo坝上游断面的总比能q，x单宽流量，坝面流速系数A，常数B，常数C，常数linjieH引入的函数Hc下游收缩断面水深PublicFunctionlinjieH(x!,q!,hc!)'建立函数Dimh1!,A!,B!,C!,fhc!,fh1c!g=Eo=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/fh1cIfAbs(hc2-hc1)<eThenExitDohc=hc2LoopEndFunctionPrivateSubCommand1_Click()CalllinjieH(x!,q!,hc!)'调用函数PrintFormat(hc,"")EndSub3．二等分法求方程的根。选取区间[a,b]使，则其根在区间[a,b]内。1．2．继续二等分。每二等分一次区间缩短为原来的一半。第k次二等分后区间长度为：收敛判别式：计算步骤：（1）选取；使；（2）令，求；（3）若，则；若，则；（4）若，迭代结束，（或或）；若，返回（2）继续进行二等分。例1-3．梯形渠道底宽，求。解：，选取。其中。程序：变量说明表：程序中：意义：Q，e流量，允许误差Ac，BcHc临界水深hc1，hc2区间（a,b）俩端点水深hc3区间（a,b）的中点值b,m渠底宽度，边坡系数Dimb!,m!,Q!,hc!,g!,e!,Ac!,Bc!,hc1!,x!,hc2!,f1!,f2!,f3!,hc3!,p!PublicFunctionLinjieH(b!,m!,Q!,hc!)g=Ac!=(b+m*hc)*hcBc!=b+2*m*hcp=Ac^3/Bc-Q^2/gEndFunctionPrivateSubCommand1_Click()g=b!=Val(InputBox("请输入矩形断面的底宽b"))m!=Val(InputBox("请输入边坡系数m"))Q!=Val(InputBox("请输入流量Q"))e!=Val(InputBox("请输入计算精度e"))hc1!=x!=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!)Iff1*f3<0Thenf2=f3:hc2=hc3Elsef1=f3:hc1=hc3EndIfLoopUntilAbs(f3)<eCallLinjieH(b!,m!,Q!,hc!)Printhc3EndSub

第2节线性方程组的解法1．高斯-赛德尔迭代法三元一次方程组：初值：雅可比迭代法：赛德尔迭代法：第步：雅可比迭代法：赛德尔迭代法：收敛判别式：对于计算步骤：1． 给定初值：2． 迭代计算：（1）雅可比迭代法：（2）赛德尔迭代法：3．收敛判别式：赛德尔迭代法程序（4）变量说明表：程序中：意义：a(i,j)方程组系数矩阵元素b(i)方程组右端项元素x0(i)，x1(i)方程组待求的未知量N矩阵的维数DimiAsInteger,jAsInteger,nAsInteger,sumAsInteger,s1AsDouble,s2AsDoubleDima()AsInteger,b()AsInteger,x0()AsDouble,x1()AsDouble,EPSAsDouble,t#,z#PublicFunctionGaussSeidel()n=InputBox("输入矩阵的维数：")ReDimPreservea(n,n),b(n),x0(n),x1(n)Fori=1TonForj=1Tona(i,j)=InputBox("输入a(i,j)的值：","a("&i&","&j&")")Printa(i,j);NextjPrintNextiFori=1Tonb(i)=InputBox("输入b(i)的值：","b("&i&")")x0(i)=0Printb(i);NextiPrint'此处判断该系数矩阵是否为严格对角占优阵,若不想判断也可将下10行删除'Fori=1Ton'Forj=1Ton't=Abs(a(i,i))'z=z+Abs(a(i,j))'Nextj'Ift<z-Abs(a(i,i))Then'Print"您输入的方程组的Gauss-Seidel迭代式不收敛！";'EndIf'Ift<z-Abs(a(i,i))ThenExitFunction'Ift>=z-Abs(a(i,i))ThenGoToF'NextiF:DoFori=1TonForj=i+1Tons2=s2+a(i,j)*x0(j)NextjForj=1Toi-1s1=s1+a(i,j)*x1(j)Nextjx1(i)=(b(i)-s1-s2)/a(i,i)EPS=Abs(x1(1)-x0(1))s1=0s2=0NextiFori=1Tonx0(i)=x1(i)Nextisum=sum+1IfEPS<=ThenExitDoLoopPrint"迭代了";sum;"次"Print"x(1)=";Format(x0(1),"");"x(2)=";Format(x0(2),"");"x(3)=";Format(x0(3),"");'声明:此程序是系数矩阵满足严格对角占优阵的情况下,Gauss-Seidel迭代法收敛;对于满足不可约对角占优阵的Gauss-Seidel迭代法收敛情况下,'有兴趣的读者可查看有关资料,自己动手试试3x1+x2+6x3=3x1+x2+2x3=3-x2+2x2-3x3=1x1=-10x2=3x3=5EndFunctionPrivateSubForm_Click()CallGaussSeidelEndSub雅可比迭代法程序DimiAsInteger,jAsInteger,nAsInteger,sumAsInteger,s1AsDouble,s2AsDouble,t#,z#Dima()AsInteger,b()AsInteger,x0()AsDouble,x1()AsDouble,EPSAsDoublePublicFunctionJacobi()n=InputBox("输入矩阵的维数：")ReDimPreservea(n,n),b(n),x0(n),x1(n)Fori=1TonForj=1Tona(i,j)=InputBox("输入a(i,j)的值：","a("&i&","&j&")")Printa(i,j);NextjPrintNextiFori=1Tonb(i)=InputBox("输入b(i)的值：","b("&i&")")x0(i)=0Printb(i);NextiPrintFori=1TonForj=1Tont=Abs(a(i,i))z=z+Abs(a(i,j))NextjIft<z-Abs(a(i,i))ThenPrint"您输入的方程组的Jacobi迭代式不收敛！";EndIfIft<z-Abs(a(i,i))ThenExitFunctionIft>=z-Abs(a(i,i))ThenGoToFNextiF:DoFori=1TonForj=i+1Tons2=s2+a(i,j)*x0(j)NextjForj=1Toi-1s1=s1+a(i,j)*x0(j)Nextjx1(i)=(b(i)-s1-s2)/a(i,i)EPS=Abs(x1(1)-x0(1))s1=0s2=0NextiFori=1Tonx0(i)=x1(i)Nextisum=sum+1IfEPS<=ThenExitDoLoopPrint"迭代了";sum;"次"Print"x(1)=";Format(x0(1),"");"x(2)=";Format(x0(2),"");"x(3)=";Format(x0(3),"");'声明:此程序是系数矩阵满足严格对角占优阵的情况下,Jacobi迭代法收敛;对于满足不可约对角占优阵的Jacobi迭代法收敛情况下,'有兴趣的读者可查看有关资料,自己动手试试'10x1-x2-2x3=72-x1+10x2-2x3=83-x1-x2-5x3=42解为：x1=x2=x3=EndFunctionPrivateSubForm_Click()CallJacobiEndSub2．高斯列主元消去法计算步骤：1．选主元（第k步）（k=1,2,…,n）假设主元在第L行换元：把换为主元（j=k,k+1,…,n,n+1）2．消元（第k步）（k=1,2,…,n）（如第2次消元，即k=2，第3行，即i=3，第4列，即j=4,数据的变化应为：）进行到第n步得：3．回代例：高斯列主元消去法求解线性代数方程组的程序⑸变量说明表：程序中：意义：a(i,j)线性方程组的系数矩阵元素n线性方程组的维数b(i)方程组右端项元素maxV对角元素的最大值maxI，maxJ最大元素在I行J列OptionExplicitPrivateSubCommand1_Click()Dima()AsDouble,b()AsDouble,iAsInteger,nAsInteger,j!n=3ReDima(n,n)AsDouble,b(n)AsDoublen=InputBox("输入矩阵的维数：")ReDimPreservea(n,n),b(n)Fori=1TonForj=1Tona(i,j)=InputBox("输入a(i,j)的值：","a("&i&","&j&")")Printa(i,j);NextjPrintNextiFori=1Tonb(i)=InputBox("输入b(i)的值：","b("&i&")")Printb(i);NextiPrintIfSolve(a(),b())ThenFori=1TonPrint"x(";i;")=";b(i),NextiElsePrint"Thisisasinglematrix!"EndIfPrintEndSubPrivateFunctionSolve(a()AsDouble,b()AsDouble)AsBooleanDimiiAsInteger,iAsInteger,jAsIntegerDimi1AsInteger,d1AsDoubleDimnAsInteger,index()AsInteger,x()AsDoubleDimmaxVAsDouble,maxIAsInteger,maxJAsIntegern=UBound(b)ReDimindex(n)AsInteger,x(n)AsDoubleFori=1Tonindex(i)=iNextiForii=1Ton-1maxV=Abs(a(ii,ii))maxI=iimaxJ=iiFori=iiTonForj=iiTonIfAbs(a(i,j))>maxVThenmaxV=Abs(a(i,j))‘将主对角线上的最大元素找出maxI=i'rowofmaxelementmaxJ=j'columnofmaxelementEndIfNextjNextiIfii<>maxIThen'exchangerowsd1=b(maxI)b(maxI)=b(ii)b(ii)=d1Forj=iiTond1=a(maxI,j)a(maxI,j)=a(ii,j)a(ii,j)=d1NextjEndIfIfii<>maxJThen'exchangecolumnsi1=index(maxJ)index(maxJ)=index(ii)index(ii)=i1Fori=1Tond1=a(i,maxJ)a(i,maxJ)=a(i,ii)a(i,ii)=d1NextiEndIfIfa(ii,ii)=0#ThenExitFunctionEndIfFori=ii+1Ton'eleminationd1=a(i,ii)/a(ii,ii)Forj=ii+1Tona(i,j)=a(i,j)-a(ii,j)*d1Nextjb(i)=b(i)-b(ii)*d1NextiNextiix(n)=b(n)/a(n,n)'backsubstitutionFori=n-1To1Step-1d1=0#Forj=i+1Tond1=d1+a(i,j)*x(j)Nextjx(i)=(b(i)-d1)/a(i,i)NextiFori=1Tonb(index(i))=x(i)NextiSolve=TrueEndFunction

第三节插值例已知消力坎的淹没系数与相对淹没度对应关系如下表：…….…….…….……..hs为下游水位超过坎顶的高度；H0为坎上全水头。试编一分别用线性插值和抛物型插值的程序，计算相对淹没度时的淹没系数值。1．线性插值已知，求，要求满足则令2．拉格朗日插值二次插值（抛物型插值），求。设满足条件：3．n次插值（拉格朗日插值），求线性插值程序⑹变量说明表：程序中：意义：x(i)已知节点坐标y(i)已知节点函数值m待求函数值的节点坐标Dimm!,f1!,f2!,L1x!,L2x!,Lx!PrivateSubCommand1_Click()Dimx(4)AsSingle,y(4)AsSinglex(0)=:x(1)=:x(2)=:x(3)=:x(4)=y(0)=:y(1)=:y(2)=:y(3)=:y(4)=m=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=i+1LoopUntilf1*f2<0PrintFormat(Lx,"")EndSub

第4节 拟合堰流：最小二乘法：偏差若要使取最小值，则①②例：直角三角形薄壁堰，求时，最小二乘法线性拟合直角三角形薄璧堰流量的程序⑺变量说明表：程序中：意义：H直角三角形薄壁堰的堰上水头t1(i)，t2(i)堰上水头，流量，一维数组sum1，sum2公式中AA1，AA2，V中间变量Q欲求流量PrivateSubCommand1_Click()DimH#H=Val(InputBox("请输入直角三角形薄壁堰的堰上水头"))Callnehe(H#)EndSubPublicFunctionnehe(H#)Dimt1(8)AsDouble,t2(8)AsDouble,x(8)AsDouble,y(8)AsDouble,sum1#,sum2#,p1#,p2#,AA1#,AA2#,V#,r1#,r2#t1(0)=:t1(1)=:t1(2)=:t1(3)=:t1(4)=:t1(5)=:t1(6)=:t1(7)=:t1(8)=t2(0)=:t2(1)=:t2(2)=:t2(3)=:t2(4)=:t2(5)=:t2(6)=:t2(7)=:t2(8)=Fori=0To8y(i)=Log(100*t2(i))NextiFori=0To8x(i)=Log(10*t1(i))NextiFori=0To8AA1=x(i)AA2=y(i)sum1=sum1+AA1sum2=sum2+AA2NextiV=sum2-sum1Fori=0To8r1=x(i)*y(i)*V*9r2=x(i)^2*V*9p1=p1+r1p2=p2+r2Nextim=p1/p2b=Abs(sum2-m*sum1)/9Q=Exp(b+m*Log(H*10))PrintQEndFunction

第二章偏微分方程式的差分数值解法第1节分类、数值解法和差分格式1．分类二阶线性偏微分方程式：典型例子：椭圆型偏微分方程式：拉普拉斯方程式抛物型偏微分方程式：扩散方程式双曲型偏微分方程式：波动方程式2．数值解法主要步骤：①划分网格；②构造差分方程式；③求解差分方程式。3．差分格式函数及在处的泰勒级数展开式分别为：（1）向前差分（一阶精度）（2）向后差分（一阶精度）（1）-（2）得：中心差分（二阶精度）（1）+（2）得：中心差分（二阶精度）

第2节 椭圆型偏微分方程式的数值解法基本方程式：（在域A内）边界条件：（1）第一类边界条件：（2） 第二类边界条件：（3）第三类边界条件：有限差分解法：（1）网格划分：（2）构造差分方程式：，五点格式：（3）解差分方程式1．高斯迭代法2．高斯-赛德尔迭代法3．超松弛迭代法ω超松弛因子（ω=1时为高斯-赛德尔迭代法，ω〈1为亚松驰迭代，ω〉1为超松弛迭代）（计算试验表明：ω取~之间数值时迭代次数最小）判别收敛性：例：如图所示为夹在两个无限平行壁间的圆柱绕流，无限远处的来流速度分布均匀，且，二平板间通过的流量为，圆柱半径，试用有限差分法编写求流函数在流场中分布的程序。解：（一）网格划分：，，。（二）构造差分方程式：（在域A内）第一类边界条件：（1）；（2）；；；（3）；第二类边界条件：（4）。五点差分格式：超松弛迭代法：（ω=~）边界条件：（1）（ab边界）；（2）；（3）bc边界的处理：(a)(b)；（4）；（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，sdij1PrivateSubCommand1_Click()Dimpsi(8,5),w,sdij0,sdij1Forj=1To5Fori=1To8psi(i,j)=0NextiNextjForj=2To4psi(1,j)=(j-1)*NextjFori=2To8psi(i,5)=2NextiDow=k=0sdij0=0k=k+1sdij1=0psi(6,2)=psi(5,2)*/psi(7,3)=psi(7,4)*/psi(8,4)=(2*psi(7,4)+2)/4Forj=2To4Fori=2Toj+3dij=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)+dijNextiNextjIfAbs(sdij1-sdij0)<ThenExitDosdij0=sdij1LoopForj=1To5PrintFori=1To8PrintSpace(1),Format(psi(i,j),"");NextiNextjEndSub

第3节抛物型偏微分方程式的数值解法例：两无限平板间油的流动，上平板水平运动，内u=0线性增加到，其后保持不变。求两平板间的速度分布随时间的变化。解：一．基本方程式边界条件：初始条件：二．有限差分解法（一）网格划分取空间步长其中；取时间步长其中。（二）构造差分方程式显式格式：若令，则稳定性条件：隐式格式：若令，则隐式格式无条件稳定。初始条件：边界条件：（三）解差分方程式显格式差分法计算两平板间速度分布随时间变化的程序⑼变量说明表：程序中：公式中：意义：n，Hn，a/m时段数，距离步长MM距离a被等份的份数AA两平板间的距离Vν液体的运动粘滞系数U0平板的最终速度RR显格式的稳定系数u(i,j)点流速Tj中间变量DT△t时间步长PrivateSubCommand1_Click()Dimm,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=:count=0DT=R*H^2/vFori=1Tom+1u(i,1)=0NextiForj=1Ton+1u(1,j)=0t(j)=(j-1)*DTtj=t(j)Iftj<Thenu(m+1,j)=tj*U0/Elseu(m+1,j)=U0EndIfNextjForj=1TonFori=2Tomu(i,j+1)=R*u(i+1,j)+(1-2*R)*u(i,j)+R*u(i-1,j)NextiNextjForj=n+1To1Step-50PrintjFori=1Tom+1count=count+1PrintSpace(1),Format(u(i,j),"");Ifcount=11Thencount=0:PrintEndIfNextiNextjEndSub

第3章有压管流第1节 管网的水力计算哈迪克劳斯法例1：如图所示环状管网，假设。试求：（1） 各管段的流量分配；（2） 各节点的测压管水头。管段12345678910管长600600300600300900300600300600管径沿程阻力系数解：（一）主要公式连续性方程：。n为该点的管段数目，流进节点，流出节点。能量方程式：。环的正方向与该管段流向相同时，否则，，所以：（二）（1）假设各管段的流量分配，使其满足连续性方程；（2）A环路水头残差；设的修正流量为，使；即：所以：，定性分析其正确性：如果初分流量时只是数值误差，不是方向“误差”的话，分母总为正值。则若分子（即残差）为正，需减小。亦即正的水头损失太大了，需减小流量；而负的水头损失太小了，需增加流量。这是正确的。反之亦然。（3）计算步骤：1．根据各管段的计算；2．任意选择各管段中的，满足连续性方程（给定各管段的流动方向）；3．选定各闭合环路的正流动方向（如顺时针方向）来确定各管段的，并构成矩阵，k-环路编号（行），I-管段编号（列）。对于环路中不包含的管段，矩阵元素。当该管段流向与该环方向相同时，；当该管段流向与该环方向相反时，。矩阵：（第k闭合环路--）ik123456789101-1+100000-10+1200+1+10000-1-130000+1-1-1+1+104.计算（对每个闭合环路），扩展为：。不在该环时，；在该环时。5．计算各管段中的流量，，扩展为：。对于非共用管，只有一个不为零；对于平面共用管，只有两个不为零，并且当均规定顺时针方向为环路正方向时，一个为“+1”，一个为“-1”。6．判别计算误差（如的允许误差），是否满足：（其中为总环数），若不满足，返回第4步继续进行各环路修正流量的计算。7．计算各节点的测压管水头（需构造节点及管段号之间的关系数组）。几点解释：（1）只管数值，不管方向。流量方向事先已假定，最终只要不修正为负值，方向就与原假设方向相同。如若修正为负值，则与原假设方向相反。（2）连续性方程中流量方向用符号函数表示。流入为“+1”，流出为“-1”。（3）能量方程式中只表示修正流量的数值。为正值时，即为均按顺时针方向增加流量；为负值时，即为均按逆时针方向增加流量。（4）的意义为：；。哈迪-克劳斯法管网的水力计算的程序⑾变量说明表：程序中：公式中：意义：d(i)，l(i)d,l各管段的管径，管长y(i)λ各管段的沿程阻力系数R(i)R公式中的系数Q(i，j)Q第j次迭代第i个管段的流量M(i，j)表示管段中的流动方向是否与环路方向一致的符号矩阵H(i)H各管段的水头B,A中间变量c(i)，x(i)中间变量PrivateSubForm_Click()Dimd(1To10),l(1To10)AsInteger,y(1To10),R(1To10),Q(1To10,100)AsDouble,w#,M(1To3,1To10)DimQ1(1To3)AsDouble,B,A,c(1To10),H(1To8),x(1To10)‘输入已知数据d(1)=:d(2)=:d(3)=:d(4)=:d(5)=:d(6)=:d(7)=:d(8)=:d(9)=:d(10)=l(1)=600:l(2)=600:l(3)=300:l(4)=600:l(5)=300:l(6)=900:l(7)=300:l(8)=600:l(9)=300:l(10)=600y(1)=:y(2)=:y(3)=:y(4)=:y(5)=:y(6)=:y(7)=:：y(8)=:y(9)=:y(10)=Fori=1To10‘求系数R(i)=8*y(i)*l(i)/*^2*d(i)^5)NextiQ(1,0)=:Q(2,0)=:Q(3,0)=:Q(4,0)=:Q(5,0)=:Q(6,0)=:Q(7,0)=Q(8,0)=:Q(9,0)=:Q(10,0)=M(1,1)=-1:M(1,2)=1:M(1,3)=0:M(1,4)=0:M(1,5)=0:M(1,6)=0:M(1,7)=0:M(1,8)=-1:M(1,9)=0:M(1,10)=1M(2,1)=0:M(2,2)=0:M(2,3)=1:M(2,4)=1:M(2,5)=0:M(2,6)=0:M(2,7)=0:M(2,8)=0:M(2,9)=-1:M(2,10)=-1M(3,1)=0:M(3,2)=0:M(3,3)=0:M(3,4)=0:M(3,5)=1:M(3,6)=-1:M(3,7)=-1:M(3,8)=1:M(3,9)=1:M(3,10)=0A=0:B=0Do‘利用迭代法求各管段的流量T:Fork=1To3Fori=1To10A=A+M(k,i)*R(i)*Q(i,0)^2B=2*Abs(M(k,i))*R(i)*Q(i,0)+BQ1(k)=-A/BNextiNextkFori=1To10Forj=1To3c(i)=0c(i)=c(i)+M(j,i)*Q1(j)Q(i,1)=Q(i,0)+c(i)IfQ(i,1)<>Q(i,0)ThenQ(i,0)=Q(i,1)EndIfNextjNextiFori=1To10Fork=1To3IfAbs(Q1(k)/Q(i,0))<ThenExitDoNextkNextiLoopFori=1To10Print"Q(";i;")=";Format(Q(i,0),"");PrintNextiFori=1To10‘求各节点的水头x(i)=R(i)*Q(i,0)*Abs(Q(i,0))NextiH(1)=40H(7)=H(1)-x(1):H(2)=H(1)-x(2):H(3)=H(2)-x(3):H(4)=H(3)-x(4)H(5)=H(4)-x(5):H(6)=H(5)+x(6):H(8)=H(7)-x(8)Fori=1To8Print"H(";i;")=";Format(H(i),"")NextiEndSub

第2节管网的水力计算有限单元法哈迪—克劳斯法是逐步渐进的经典解法，在大管径小流量情况下不能很快收敛。一、单元方程式水流从令则（流量流出节点为“+”，流入节点为“—”）构成整体方程组（加入流量为“+”，消耗流量为“—”）以节点2为例：单元（1） 单元（3）单元（4）所以i单元，k,j两个节点，ki放到（k,k）,(j,j)位置上，-ki放到（k,j）,(j,k)位置上。单元编号12345始端13223末端21344求解总体方程组（边界条件与解法）为了求解方程组，至少必须已知一个节点的水头值。由于流量只取决于水头差而不是水头本身，如果需要，可以任意规定某个节点的水头值，对流量计算的结果也毫无影响。边界条件引入的两种方法：如本例四个节点的连续性方程，只有三个节点是独立的。即其系数矩阵的行列式值为零，有无穷多个解。只有引入一个边界条件（一节点的水头值），才可求解。这也符合其物理意义。若：第一种处理方法：不减少方程的个数，直接引入。第二种处理方法：减少方程的个数，移项处理。高斯—赛德尔迭代法：当流动在紊流的阻力平方区时，单元特征系数和总体矩阵[k]的各元素都不是常数。它们都是水头的函数。因此总体矩阵方程式是非线性的代数方程组，即解非线性代数方程组可采用高斯-塞德尔迭代法，牛顿-拉普森切线法等。牛顿-拉普森切线法比高斯-塞德尔迭代法收敛速度快。但高斯—赛德尔迭代法编制程序比较简单。若用高斯—赛德尔迭代法计算步骤如下：（1）假设待求各节点的水头值，即给定迭代初值；（2）用初值计算矩阵；（3）用计算待求的各节点的水头值，顶标“1”表示迭代次数；（4）将求得的水头值与初值相比较，如果满足下式：计算结束。否则将符给，重复（2）~（4）步骤，直到满足精度要求为止。（5）求各单元的流量。例1。单元编号12345678910始端1123456742末端7234567888有限单元法管网水力计算程序（12）变量说明表：程序中：公式中：意义：d(i)，l(i)d,l各管段的管径，管长y(i)n各管段的沿程粗糙系数Q(i)Q第j次迭代第i个管段的流量H0(i)，c(i)H0j,cj节点水头的初值，“节点消耗”H1(i)Hj节点的水头n1(i)，n2(i)单元端点编号A(ii,jj)总体矩阵元素，二维数组K(i)单元特征系数，一维数组s1，s2中间变量PrivateSubForm_Click()Dimd(1To10),l(1To10)AsInteger,y(1To10),Q(1To10)AsDoubleDimA(10,10),c(1To8),H0(1To8),H1(1To8),K(10),n1(10),n2(10)‘输入已知数据d(1)=:d(2)=:d(3)=:d(4)=:d(5)=:d(6)=:d(7)=:d(8)=:d(9)=:d(10)=l(1)=600:l(2)=600:l(3)=300:l(4)=600:l(5)=300:l(6)=900:l(7)=300:l(8)=600:l(9)=300:l(10)=600y(1)=:y(2)=:y(3)=:y(4)=:y(5)=:y(6)=:y(7)=:y(8)=:y(9)=:y(10)=c(1)=:c(2)=:c(3)=:c(4)=:c(5)=c(6)=:c(7)=:c(8)=H0(1)=40:H0(2)=33:H0(3)=31:H0(4)=30:H0(5)=28:H0(6)=31:H0(7)=35:H0(8)=n1(1)=1:n1(2)=1:n1(3)=2:n1(4)=3:n1(5)=4n1(6)=6:n1(7)=7:n1(8)=7:n1(9)=8:n1(10)=2n2(1)=7:n2(2)=2:n2(3)=3:n2(4)=4:n2(5)=5n2(6)=5:n2(7)=6:n2(8)=8:n2(9)=4:n2(10)=8Forii=1To8Forjj=1To8A(ii,jj)=0NextjjNextiiFori=1To10‘求总体矩阵元素ii=n1(i)jj=n2(i)K(i)=0K(i)=*d(i)^/(y(i)*Sqr(l(i))*Sqr(Abs(H0(ii)-H0(jj))))A(ii,jj)=A(ii,jj)-K(i)A(jj,ii)=A(jj,ii)-K(i)A(ii,ii)=A(ii,ii)+K(i)A(jj,jj)=A(jj,jj)+K(i)NextiDo‘迭代求水深Fori=1To8Forj=i+1To8H1(1)=H0(1)s2=s2+A(i,j)*H0(j)NextjForj=1Toi-1s1=s1+A(i,j)*H0(j)NextjH1(i)=(c(i)-s1-s2)/A(i,i)EPS=Abs(H1(i)-H0(i))s1=0s2=0NextiFori=1To8H0(i)=H1(i)NextiIfEPS<=ThenExitDoLoopFori=1To8PrintFormat(H1(i),"")NextiFori=1To10‘用已知公式求流量ii=n1(i)jj=n2(i)Q(i)=K(i)*(H1(ii)-H1(jj))PrintFormat(Q(i),"")NextiEndSub

第3节简单管道的水击计算特征线法（1）基本方程式（坡度很小）非恒定流连续性方程：f—沿程阻力系数（λ）运动方程式：所以顺特征线方程：顺特征方程：逆特征线方程：逆特征方程：（2）差分方程式(N个网格)顺特征差分方程：逆特征差分方程：（3）边界条件上游边界条件：上游为水库：下游边界条件：若按孔口出流计算（阀门断面）：为阀门相对开度。特征线法计算简单管路中水击压强的程序⒀例：有一如图所示的水电站的引水钢管，在钢管末端装有自动调节阀门，使阀门按直线规律关闭，钢管长，水库最高水位与管道阀门断面的高差，阀门全开时管中最大流速，水击波速，设阀门完全关闭的时间，试编写程序求：阀门断面处的最大水击压强。解：（一）.主要公式（1）.基本方程式对于坡度较小的有压管道，非恒定流的连续方程式和运动方程式为：式中是沿程阻力系数，而其中波速计算公式：变量说明表：程序中：公式中：意义：DD管径LL钢管的长度H0h0水头计算初值V0ν0流速计算初值CC=c/g水击速度Ts闸门完全关闭时间F沿程阻力系数NN=l/dx管段数Nt时段数H(I,j),v(I,j)管中水头及流速，二维数组T(i),tn(i)t,时间及阀门开度PrivateSubCommand1_Click()Dimh(10,100),v(10,100),t(100),tn(100)Dimd,l,h0,v0,c,ts,f,n,nt,cm#,v1#,v2#,v3#f=Val(InputBox("沿程阻力系数f"))v0=Val(InputBox("流速计算初值v0"))h0=Val(InputBox("水头计算初值h0"))l=Val(InputBox("钢管的长度l"))ts=Val(InputBox("闸门完全关闭时间ts"))c=Val(InputBox("水击速度c"))dx=500n=l/dxnt=40d=2g=dt=l/(n*c)c1=c/gc2=f*dt/(2#*d)c3=c1*c2t(1)=0#tn(1)=1#Fori=1To6‘边界条件赋值h(i,1)=h0-(i-1)*c3*v0^2v(i,1)=v0NextiForj=2Tont+1‘特征法进行运算Fori=2To5h1=h(i-1,j-1)h2=h(i+1,j-1)v1=v(i-1,j-1)v2=v(i+1,j-1)w1=(h1-h2)/c1w2=v1+v2w3=v1*Abs(v1)+v2*Abs(v2)v(i,j)=*(w1+w2-w3*c2)w4=h1+h2w5=c1*(v1-v2)w6=v1*Abs(v1)-v2*Abs(v2)h(i,j)=*(w4+w5-c3*w6)Nexticm=h(2,j-1)-v(2,j-1)*(c1-c3*Abs(v(n,j-1)))cp=h(n,j-1)+v(n,j-1)*(c1-c3*Abs(v(n,j-1)))h(1,j)=h0v(1,j)=(h(1,j)-cm)/c1t(j)=(j-1)*dtIft(j)-ts<>0Thentn(j)=1#-(j-1)*dt/tsc6=(v0*v0*tn(j)*tn(j))/h0c7=c6*c1/2#v(n+1,j)=-c7+Sqr(c7*c7+cp*c6)h(n+1,j)=cp-c1*v(n+1,j)Elsetn(j)=0#v(n+1,j)=0#h(n+1,j)=cp-c1*v(n+1,j)EndIfNextjForj=1TontPrintFormat(t(j),""),Format(tn(j),"")Fori=2Ton+1PrintFormat(h(i,j),"")NextiNextjEndSub

第4节短管的水力计算短管的水力计算程序（10）例：如图所示为新的铸铁管道，已知上下游水位差，管长，总局部阻力系数（不包括出口损失），管壁的当量粗糙度，水温为，其运动粘滞系数，试编写程序求：当流量时的管径；当管径时的流量。解：①管径d的计算：，其中。牛顿切线法：计算过程：⑴假设λ；⑵求A,B,用牛顿切线法求d；⑶求ν,Re,计算λ；⑷比较λ相对误差是否满足精度要求。若满足精度要求则计算结束，否则返回⑵继续循环计算。②已知，求Q：计算过程：⑴假设λ；⑵计算μs,v;⑶计算Re，选取公式，求λ;⑷比较的相对误差是否满足精度要求：若满足，求Q，结束计算；若不满足，则返回⑵循环计算。③已知Q,d,求H：代入，即可求得H。关于沿程阻力系数计算的讨论：上述的第二式：称为布拉休斯公式。此式适用于紊流的光滑区。又有适用于紊流粗糙区的经验公式为：称为希弗林松公式。上述的第三式：可看成是紊流光滑区的布拉休斯公式与粗糙区的希弗林松公式的综合。亦即可看成是适用于紊流三个区的通用公式。如若采用上一套计算公式进行计算，无论计算d或Q，均先假设，把的计算纳入到d的牛顿切线法迭代计算或计算Q的简单迭代法循环中即可。若采用柯列布鲁克—怀特公式：求，则可用二等分法计算。构造关于的函数如下：设，则若，则判别下式是否成立：若成立，计算结束；否则，返回（*）步继续进行二等分。柯列布鲁克—怀特公式是适用于紊流过渡粗糙区的公式，亦可以看成是适用于紊流三个区的通用公式。变量说明表：程序中；公式中：意义：hH水头QQ流量LL管长DD管径Sz∑ζ未出口的局部阻力系数KsKs管径的当量粗糙度Nuν水的运动粘滞系数Fλ沿程阻力系数fEps1,eps2管径和沿程阻力系数的允许的误差j通道管号i循环变量reRe雷诺数PrivateSubCommand1_Click()Dimv#h=Val(InputBox("水位差h"))q=Val(InputBox("流量q"))l=Val(InputBox("管长l"))d=Val(InputBox("管径d"))sz=Val(InputBox("末计出口的局部阻力系数sz"))ks=Val(InputBox("管径的当量粗糙度ks"))nu=Val(InputBox("水的运动粘滞系数nu"))f=Val(InputBox("沿程阻力系数f"))eps1=Val(InputBox("管径和沿程阻力系数的允许的误差eps1"))eps2=Val(InputBox("管径和沿程阻力系数的允许的误差eps2"))j=Val(InputBox("通道管号j"))pi=g=Ifj=1ThenFori=1To100v=4#*q*q/(pi*d*d)re=v*d/nu‘雷诺数的计算Ifre<2300Then‘沿程阻力系数的判断f=64/reElseIfre>2300Andre<100000Thenf=/re^Elsef=*(68/re+ks/d)^EndIfa=8#*q*q*f*l/(g*pi*pi)b=8#*q*q*(1+sz)/(g*pi*pi)fd=a/d^5#+b/d^4#-hfd1=-5#*a/d^6#-4#*b/d^5#d1=d-fd/fd1IfAbs(d1-d)<eps1ThenPrintd=d1PrintFormat(d,"")GoToww‘进行循环迭代Elsed=d1EndIfNextiElsei=1qq:us=1/Sqr(f*l/d+sz+1#)'是求沿程阻力系数的循环v=us*Sqr(2*g*h)re=v*d/nuIfre<2300Thenf1=64#/reElseIfre>2300Andre<100000Thenf1=/re^Elsef1=*(68#/re+ks/d)^EndIfIfAbs(f1-f)<eps2Thenq=*d*d*vf=f1PrintFormat(i,""),Format(f,""),Format(q,"")Elsef=f1i=i+1GoToqqEndIfEndIfww:EndSub

明渠流第1节明渠非恒定流的计算特征线法基本方程组（以为变量，梯形断面明渠）其中特征线方程和特征方程顺特征线方程：顺特征方程：即：逆特征线方程：逆特征方程：即：特征线差分方程和特征差分方程（库朗格式）顺特征线差分方程：顺特征差分方程：即：逆特征线差分方程：逆特征差分方程：即：点L的值据点A和M的值线性插值计算；点R的值据点M和B的值线性插值计算。同理同理上两式相减得：由顺特征差分方程求：或由逆特征差分方程求：边界条件上游：已知，则下游：给定关系曲线。稳定性条件要求时间步长满足下式：例：梯形断面明渠渠中为均匀流，末端断面水深。上游洪水过程如下表：t(s)060120180240300360420Q1517192123252524t(s)480540600660720780840900Q2322212019181715特征线法计算明渠非恒定流的程序⒄如图所示，某梯形断面渠道，下游与一水塘相接，已知：渠道长度，底宽，边坡系数，粗糙系数，底坡，通过流量时渠中为均匀流，其水深，这时渠道末端断面水位与水塘中水位平齐，其后末端断面水深随流量按下式变化：今在渠道的上游流域降一场大雨，在上游端断面测得洪水过程如下表所给：06012018024030036042015171921232525244805406006607207808409002322212019181716试用特征线法编写程序，计算各时刻各断面的水深、流量、流速及波速。变量说明表：程序中公式意义MM边坡系数NN粗糙系数BB宽度S0I底坡H0H渠中水深Ds距离步长Dt时间步长A1A方程式中的系数a1C1C方程式中的系数c1H(l,k),Q(l,k)H,q水深，流量数组V(l,k),C(I,k)v,c=(g*A/B)^(1/2)断面平均流速，波速数组PublicSubCommand1_Click()Dimh(12,25),q(12,25),v(12,25),c(12,25),qu(16)m=Val(InputBox("边坡系数m"))n=Val(InputBox("粗糙系数n"))b=Val(InputBox("b"))s0=Val(InputBox("底坡s0"))h0=Val(InputBox("渠中水深h0"))ds=Val(InputBox("距离步长ds"))dt=Val(InputBox("时间步长dt"))a1=Val(InputBox("方程式中的系数a1"))b1=Val(InputBox("方程式中的系数b1"))c1=Val(InputBox("方程式中的系数c1"))l=11k=25n1=16Forj=1Tok‘边界条件Fori=1Tolh(i,j)=0#q(i,j)=0#v(i,j)=0#c(i,j)=0#NextiNextjFori=1Tolh(i,1)=h0q(i,1)=15#Nextiqu(1)=15:qu(2)=17:qu(3)=19:qu(4)=21:qu(5)=23:qu(6)=25:qu(7)=25:qu(8)=24qu(9)=23:qu(10)=22:qu(11)=21:qu(12)=20:qu(13)=19:qu(14)=18:qu(15)=17:qu(16)=16Forj=1Ton0‘对流量赋值q(i,j)=qu(j)NextjForj=n1+1Tokq(1,j)=15#Nextjg=Forj=2Tok–1‘对内点的p计算Fori=1Tolhm=h(i,j-1)qm=q(i,j-1)bs=b+2*m*hmx=b+2*hm*Sqr(1+m*m)a=(b+m*hm)*hmv(i,j-1)=qm/ac(i,j-1)=Sqr(g*a/bs)vm=v(i,j-1)cm=c(i,j-1)wb=bs*(vm-cm)‘逆水流方向的绝对速度wf=bs*(vm+cm)‘顺水流方向的绝对速度n0=-g*a*s0+g*qm^2*n^2*x^(4#/3)/a^(7#/3)Ifi=1Thenhb=h(i+1,j-1)qb=q(i+1,j-1)qr=qm+dt/ds*(qm-qb)*wb/bshr=hm+dt/ds*(hm-hb)*wb/bsh(i,j)=hr+(q(i,j)-qr+n0*dt)/wfElseIfi>1Andi<1Thenha=h(i-1,j-1)hb=h(i+1,j-1)qa=q(i-1,j-1)qb=q(i+1,j-1)hl=hm+dt/ds*(ha-hm)*wf/bs‘已知时层内l插点的水深ql=qm+dt/ds*(qa-qm)*wf/bs‘已知时层内l插点的流量hr=hm+dt/ds*(hm-hb)*wb/bs‘已知时层内r插点的水深qr=qm+dt/ds*(qm-qb)*wb/bs‘已知时层内r插点的流量h(i,j)=1/(wb-wf)*(qr-ql+wb*hl-wf*hr)q(i,j)=qr+wf*(h(i,j)-hr)-n0*dtElseha=h(i-1,j-1)qa=q(i-1,j-1)hl=hm+dt/ds*(ha-hm)*wf/bsql=qm+dt/ds*(qa-qm)*wf/bsaa=wb*c1bb=wb*b1-1#cc=ql+wb*a1-wb*hl-n0*dtdd=bb^2-4*aa*ccq(i,j)=(-bb-Sqr(dd))/(2*aa)h(i,j)=a1+b1*q(i,j)+c1*q(i,j)^2‘末端断面的水深随流量的变化的函数EndIfN