武汉大学《水力学》梯形渠道恒定非均匀渐变流水面线计算程序 WTSF.FOR

武汉大学《水力学》梯形渠道恒定非均匀渐变流水面线计算程序 WTSF.FORmodule NrType ! 确定常数 implicit none integer,parameter::I2B=selected_int_kind(4) integer,parameter::I4B=selected_int_kind(9) integer,parameter::SP=kind(1.0) ! 单精度 integer,parameter::DP=kind(1.0d0) ! 双精度 real(DP),parameter::ERR=1.0d-6 ! 两实数差值的绝对值小于此常数可认为相等...

module NrType ! 确定常数 implicit none integer,parameter::I2B=selected_int_kind(4) integer,parameter::I4B=selected_int_kind(9) integer,parameter::SP=kind(1.0) ! 单精度 integer,parameter::DP=kind(1.0d0) ! 双精度 real(DP),parameter::ERR=1.0d-6 ! 两实数差值的绝对值小于此常数可认为相等 end module NrType module wtsf use NrType implicit none save ! 模块内全局变量均为静态变量 integer(I2B),parameter::NOUT=2 ! 输出文件号 integer(I2B)::dr ! 水面线计算方向 real(DP),parameter::G=9.8d0 ! 重力加速度 real(DP)::q,i,n,m,b,csn,srm,alfa ! 以上分别为流量,底坡比降,糙率,边坡系数,底宽,底坡余弦值,边坡余割值,动能修正系数 real(DP)::v2,j1,j2,es1,es2,ds ! 流速,水力比降,断面比能,断面间距 real(DP),allocatable::h(:),v(:),s(:) ! 水深,流速,距起始(控制)断面的距离 ! 外部函数 real(DP),external::ErFenFa ! 二分法 real(DP),external::zBrent ! 布伦特法 contains subroutine ReadData() ! 读取数据 ! 使用的全局变量有:q,i,n,m,b,NOUT write(*,*)'键盘输入:流量Q,底坡比降i,糙率n,边坡系数m,底宽b(单位均为国际单位制)' read(*,*)q,i,n,m,b ! 输入流量,底坡比降,糙率,边坡系数,底宽 if(qhk).or.((abs(hd-hk)h0).and.(hd>hk)).or.((h0>hk).and.(hd0梯形时 h1=m*h1/b ! 西北农林科技大学《水力学》书上公式选择判定参数 if(h1<3.0d0)then ! 3.0为分界点 h1=1.3d0*(sqrt(1.0d0+1.54d0*h1)-1.0d0) else h1=1.12d0*h1**0.604d0-0.39d0 end if ! h1<757时,临界水深误差小于0.3%;h1<2560时,临界水深误差小于0.7% h1=h1*b/m ! 作为(牛顿法或迭代法)迭代初值 write(*,'(a,f10.6,a)')'梯形断面临界水深迭代初值:',h1,'(m)' ! 优先选用牛顿法(速度快但不能保证收敛) h2=h1 ! h2为牛顿法迭代初值 do k=1,kMax bb=b+2.0d0*m*h2 ! 水面宽度 aa=(b+m*h2)*h2 ! 过水面积 x=h2+(aa**3/temp-bb)/(2.0d0*m-3.0d0*bb*bb/aa) ! 牛顿法 if(abs(x-h2)>ERR)then ! 精度不够,继续迭代 h2=x else ! 精度满足 write(*,*)'求解临界水深时牛顿法的迭代次数:',k return ! 求出m>0梯形时的临界水深(牛顿法) end if end do ! 牛顿法不收敛而采用迭代法(保证收敛但速度较慢),迭代初值仍为h1 write(*,*)'牛顿法不收敛而采用迭代法!' do k=1,kMax x=(temp*(b+2.0d0*m*h1))**Y/(b+m*h1) if(abs(x-h1)>ERR)then ! 精度不够,继续迭代 h1=x else ! 精度满足 write(*,*)'求解临界水深时迭代法的迭代次数:',k return ! 求出m>0梯形时的临界水深(迭代法) end if end do write(*,'(a,i2,a)')'错误:求解临界水深时迭代次数超过',kMax,'次!' stop end subroutine SubHk subroutine SubH0(x) ! 迭代法求解正常水深 ! 使用的全局变量有:q,n,i,b,srm,m,ERR real(DP),intent(out)::x real(DP)::temp,h1 ! 临时数据 integer(I2B),parameter::kMax=40 ! 最大迭代次数 integer(I2B)::k temp=(q*n/sqrt(i))**0.6d0 ! 迭代时还有用 h1=temp/b**0.6d0 ! 迭代初值 write(*,'(a,f10.6,a)')'正常水深迭代初值:',h1,'(m)' do k=1,kMax x=temp*(b+2.0d0*srm*h1)**0.4d0/(b+m*h1) if(abs(x-h1)>ERR)then ! 精度不够,继续迭代 h1=x else ! 精度满足 write(*,*)'求解正常水深时的迭代次数:',k return ! 求出梯形时的正常水深 end if end do write(*,'(a,i2,a)')'错误:求解正常水深时迭代次数超过',kMax,'次!' stop end subroutine SubH0 function fe(x) ! 水面线计算方程 ! 使用的全局变量有:b,m,q,n,v2,srm,csn,alfa,G,es1,es2,dr,i,j1,j2,ds real(DP),intent(in)::x real(DP)::fe,aa aa=(b+m*x)*x ! 过水面积 v2=q/aa ! 流速 j2=(n*v2)**2/(aa/(b+2.0d0*x*srm))**(4.0d0/3.0d0) ! 水力比降es2=csn*x+alfa*v2*v2/(2.0d0*G) ! 断面比能 fe=es1-es2+dr*(i-(j1+j2)/2.0d0)*ds return end function fe end module wtsf function zBrent(Func,x1,x2,tol) ! 采用布伦特法(有所改动) use NrType real(DP),intent(in)::x1,x2,tol ! tol 表示收敛的精度 real(DP)::zBrent interface function Func(x) use nrtype implicit none real(DP),intent(in)::x real(DP)::Func end function Func end interface integer(I2B),parameter::itMax=100 ! 最大迭代数 real(DP),parameter::eps=3.0d-8 integer(I2B)::iter ! 当前迭代数 real(DP)::a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm a=x1 b=x2 do ! 这个循环是为了保证区间有根 fa=Func(a) fb=Func(b) if(fa*fb<0.0d0)exit ! 为有根区间,退出循环 write(*,*)'(x1,x2)不是有根区间,请重新输入x1,x2' read(*,*)a,b ! 若不是有根区间,用户自己输入确定,直至异号为止end do c=b fc=fb do iter=1,itMax if(fb*fc>0.0d0)then c=a ! 对a,b,c更名并调整区间d fc=fa d=b-a e=d endif if(abs(fc)=tol1.and.abs(fa)>abs(fb))then s=fb/fa ! 欲进行二次反内插 if(abs(a-c)0.0d0)q=-q ! 检查是否在解区间内 p=abs(p) if(2.0d0*ptol1) ! 计算新的试验解fb=func(b) end do write(*,'(a,i2,a)')'错误:布伦特法迭代次数超过',itMax,'次!' stop ! 永远不能执行到此 end function zBrent function ErFenFa(f,x1,x2,eps) ! 采用二分法 use NrType implicit none real(DP),intent(in)::x1,x2,eps real(DP)::ErFenFa,a,b,fa,fb,fm integer(I2B),parameter::ItMax=40 integer(I2B)::iter interface use nrtype implicit none real(DP),intent(in)::x real(DP)::f end function f end interface a=x1 b=x2 do ! 这个循环是为了保证区间有根 fa=f(a) fb=f(b) if(fa*fb<0.0d0)exit ! 有根区间,退出循环 write(*,*)'(x1,x2)不是有根区间,请重新输入x1,x2' read(*,*)a,b ! 若不是有根区间,用户自己输入确定,直至异号为止end do do iter=1,ItMax ErFenFa=(a+b)*0.5d0 if(abs(b-a)

                    本文档为【武汉大学《水力学》梯形渠道恒定非均匀渐变流水面线计算程序 WTSF.FOR】，请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑，
                    图片更改请在作品中右键图片并更换，文字修改请直接点击文字进行修改，也可以新增和删除文档中的内容。 
 该文档来自用户分享，如有侵权行为请发邮件ishare@vip.sina.com联系网站客服，我们会及时删除。

                    [版权声明] 本站所有资料为用户分享产生，若发现您的权利被侵害，请联系客服邮件isharekefu@iask.cn，我们尽快处理。

                    本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权，请谨慎使用。

                    网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传，仅限个人学习分享使用，禁止用于任何广告和商用目的。
                

下载需要：免费已有0 人下载

立即下载

武汉大学《水力学》梯形渠道恒定非均匀渐变流水面线计算程序 WTSF.FOR

你可能还喜欢