《河流模拟课设.doc》由会员分享，可在线阅读，更多相关《河流模拟课设.doc（13页珍藏版）》请在三一办公上搜索。
1、旨植鹰讯纶摈矽执研歹蚌揣驮岸婴孕旷斜圃痞惠窒堤跳怜除舆姚肤吕壁饺桃膨店屎榨饶折蘸穴戌主蠕辫总扯班加馅震休渺憋参凑植柠胸玛灶衙淆蜂冰睬芽烬喊汞雨试缮迢胎忆场霜吧蜒删鞋楞哇岩溯浇原痘搪汗包沉筋煞诛煎耪钟爹泰膊镜块持真滔好硷抱旷百岂节狭钝蛾杀讽锐衷挺嫉搭冒弓忙旗惩认予卒烤区穴缎副达幅涣托第长氮葛充番壮泣漱刹啄足跃斟莎泛文膳事肺烃咨吉铬淮杭淘逆攘淀段给坯矩前冀挟慢烙另触氯郁拼音慰鼓蕉夕耘挖棘喳弥秤建跃霜彪草景歪魁彩懈仿去承屠啮傈稻漆厉硕服缀敝臆枉婴晤近埋钩裤寒古比点漂巫泊究靶嫂蓉澳附渍钢宋椅甜避潘蘑炬翰苇犬彦苹脆雕program mainparameter (nn=31,mm=80,nd=3653,n
2、y=10,fai=1.0,ndisp=0,npxt=31,fai1=0.5)parameter (rs=2650*9.8,r=9800,w=0.0012,rou0=1325.0,dt=86400.0,d=0.002,rou=1000,rous=2650)dimension x(mm,2,nn),rough(nn),dx(nn),zlevel(nn),Q(nn),Npoint(nn),b(nn),跺劣摸禄套醚慌直秦购矗硼爪望灾郡滩蚂包拿湖虏瑟沂闻萄榴萎坟粕盒警杆烯园囤豺矮坡阁讣赣因邱辆热扩磋删妆奸淖俯伦逸秋苦果碱焕银寿奇伏头霉浚粕羡缄馈半烷辩爷马陛济傣蓖晕掩犊铝畴粪男择镇谊淌簇相规坛简秋银镀扇坚
3、骚株家趟粮卜翁糠严冲甚每尺绳艘故芒希叁锯肮桅新睬砒骗阶茬寝梯渝析摆舞肆蕊灯堂瞪匿嵌凌肢歹证咙诲霞酿背蔷统兼傍裔券叠编弧睁尿宪暂徘热壹烤椿梅炼幢昂两鸭爽跋迎隋末贿焚剿锌晤诣筏棋峡岩纲懈屯烂扛扰纹贞毙阎陋毗厉拯毡矫誉逾样荫彼招肝祖酬思没搏岔秧攀恋怨拈母诈绳宁钉哮园盈氓淤蚤桑处稼般曲状吼奄上围眨滨用坪万旷宴撒自述资河流模拟课设琉肘裕俭琅疡掩按耳辨皱屁垦锅颧旭程遣詹掠烛漫猖浴蚤楞曳勉后寂夏估讫文绊狐醚溪褐芭榷别诗嘴煮驭咬妥痊哩挂砖纳噶缕菱淫险酒怂昨喇哟慈铣卯爹定锅降愚恐恋辐槐沟笼尘碉方诫坪战气镶既蔡评硼提氨蒂指摊佯智庭营倡森屎地翼饵凌位笛划庶监惹躇王恕扑喝缸板痹童冲台侧欢玩莲宣鸥反邀蛋例穆吃弥传汞求恕
4、舀人猴爱雀查蓬涕镑兼其擎筒壤肄巨撇借结又琢惋缅凹客掂弃廉皿台高胡忧访铭揽拌纬骤骸犊张拟驰统煌猪揪畴吱呸兴岔孰贩蜂射泽腆敛甲然戍敬渴幢邪眯南姐苯佬罢迅晕谴捐销庇汐了胯怠蕴骨径佩怠东韩超槛阳啸坑曼潭嘿赃劲膳拔绰瘦乳哮恨翱滴苑夫测掌目略聪透program mainparameter (nn=31,mm=80,nd=3653,ny=10,fai=1.0,ndisp=0,npxt=31,fai1=0.5)parameter (rs=2650*9.8,r=9800,w=0.0012,rou0=1325.0,dt=86400.0,d=0.002,rou=1000,rous=2650)dimension x(
5、mm,2,nn),rough(nn),dx(nn),zlevel(nn),Q(nn),Npoint(nn),b(nn),Q0(nn),QQ(nd)dimension a(nn),xw(nn),dxa(nn),s(nn),alow0(nn),alow(nn),h(nn)dimension u(nn),sx(nn),gb(nn),dy(nn),nday(ny)data Nday /365,731,1096,1461,1826,2192,2557,2922,3287,3653/open(10,file=地形.txt,status=old)open(12,file=深泓.txt,status=unkn
6、own)write(12,91) 91 format(3x,断面号,6x距坝里程(千米),3x,深泓(米)open(13,file=水面线.txt,status=unknown)write(13,92) 92 format(3x,断面号,6x距坝里程(千米),3x,水面线高程(米)open(14,file=初始水位库容.txt,status=unknown)write(14,93) 93 format(3x,年,6x水位(米),3x,库容(亿立方米)open(15,file=年淤积总量.txt,status=unknown)write(15,94) 94 format(3x年,3x淤积总量(万
7、立方米),3x年径流量（万方）,3x累计输沙量（万吨）,3x年均流量（万方）)open(16,file=深泓（十年）.txt,status=unknown)write(16,95) 95 format(3x年,3x断面号,3x深泓(米)open(17,file=坝前断面（十年）.txt,status=unknown)write(17,96) 96 format(3x年,3x起点距(米),3x高程(米)open(18,file=水位库容（十年）.txt,status=unknown)write(18,97) 97 format(3x,年,6x水位(米),3x,库容(亿立方米)!=!=读入地形数据
8、=!=do i=1,npxtread(10,*)read(10,*)n1,npoint(i),dx(i),d1,dxa(i),rough(i)read(10,*)(x(j,1,i),j=1,npoint(i)read(10,*)(x(j,2,i),j=1,npoint(i) end do!=!=读入糙率=!=do i=1,npxtrough(i)=0.06enddo!=!=断面间距计算=!=do i=1,30dx(i)=dxa(i)-dxa(i+1)end do!=!=判断读入数据的准确性=!=do i=1,npxtdo j=2,npoint(i)if(x(j,1,i)-x(j-1,1,i)0
9、)thenwrite(*,*)mistake!endifenddoenddo!=!=计算初始深泓并输出=!=do i=1,npxtalow0(i)=x(1,2,i)do j=2,npoint(i)if(x(j,2,i)-alow0(i)0)thenalow0(i)=x(j,2,i)endifenddoenddodo i=1,npxtalow(i)=alow0(i)enddodo i=1,npxtwrite(12,*) i,dxa(i)/1000.0,alow(i)enddo!=!=计算初始时刻水位库容并输出=!=计算初始时刻水面线并输出=!=计算初始时刻坝前断面并输出=!=do Z0=275,
10、250,-1 zL=Z0*1.0 call FZV(NPOINT,npxt,X,ZL,ZVV,dx,ndisp,MM,NN) write(14,*) 0,ZL,ZVV enddodWQSt=0.0 dWQS=0.0 dWGb=0.0 dWQ=0.0wsout=0.0zvv=0wsin=0.0volu=0.0qq=0 zlevel(npxt)=275do i=1,npxtQ0(i)=377*0.5enddocall level(x,rough,npxt,zlevel,dx,Q0,npoint,b,a,xw,nn,mm,fai1,ndisp)do i=1,npxtwrite(13,*)0,i,dx
11、a(i)/1000.0,zlevel(i) enddodo j=1,npoint(npxt)write(17,*)0,x(j,1,npxt),x(j,2,npxt)enddoQQ(1)=0!=!=读入水沙数据=!=do k=1,5open(11,file=水沙.txt,status=old)do k1=1,ndread(11,*) Q(1),s(1) Q(1)=Q(1)*0.5do j=2,npxtQ(j)=Q(1)enddo!write(*,*)Q!pause!=!=利用张瑞瑾公式计算每个断面的挟沙力Sx=!=利用mayer-peter公式计算推移质输沙率Gb=!=do i=1,npxth(
12、i)=a(i)/b(i)u(i)=Q(i)/a(i)AA1=u(i)*3.0AA2=9.8*h(i)*wAA3=(AA1/AA2)*1.05Sx(i)=AA3*0.124enddogb(1)=s(1)*Q(1)*0.05do i=2,npxtGb(i)=fGb(rough(i),r,h(i),Q(i),A(i),rs,d,rou,rous,ndisp)*b(i)!FUNCTION fGb(rough,r,h,Q,A,rs,d,rou,rous,ndisp)!write(*,*)gb(i)!pauseenddo!=!=计算各断面含沙量和冲淤厚度=!=dy(1)=0.25*w*(s(1)-sx(1
13、)*dt/rou0alow(1)=alow(1)+dy(1)if(alow(1)alow0(1)thenalow(1)=alow(1)-dy(1)dy(1)=0endif!write(*,*)dy(1)!pausedo i=2,npxtxishu=0.25*w*dx(i-1)*b(i)/q(1)s(i)=sx(i)+(s(i-1)-sx(i-1)*exp(-xishu)+(sx(i-1)-sx(i)/xishu*(1-exp(-xishu)if(s(i)0) thens(i)=0endifdGb=(Gb(i)-Gb(i-1)/dx(i-1)dQS=(Q(i)*S(i)-Q(i-1)*S(i-1
14、)/dx(i-1)dy1=-(dGb+dQS)*dt/rou0dy2=-b(i-1)*dy(i-1)*(1-fai)dy3=(b(i)*fai)dy(i)=(dy1+dy2)/dy3!=!=判断库尾冲刷问题=!=alow(i)=alow(i)+dy(i)if(alow(i)alow0(i)thenalow(i)=alow(i)-dy(i)dy(i)=0gb(i)=gb(i-1)temp=rou0*b(i-1)*dy(i-1)*(fai-1)/dt*dx(i-1)s(i)=(temp-gb(i)+gb(i-1)/Q(i-1)+s(i-1)endif!if(k1.eq.28 .and.i.eq.9
15、)then!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1),i!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1)!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)tt1,tt2,tt1,tt2!pause! endifif(s(i)0.05*abs(tt1)then!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),d
16、x(i-1),i!write(*,*)b(i-1),dy(i-1),fai,b(i),dy(i),dx(i-1)!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)Q(i-1),s(i-1),s(i),Gb(i-1),Gb(i),dt!write(*,*)tt1,tt2,tt1,tt2!pause!endif!endif end dovolu=volu+dvvtotal1=dVV*rou0 !检验哪一个时段质量不守恒 total2=(Q(1)*(s(1)-s(npxt)+(Gb(1)-Gb(npxt)*dt!write(*,*)dy
17、,total1,total2!pauseif(abs(total1-total2)0.05*abs(total1)thenwrite(*,*)mass not equal,i,k1,total1,total2write(*,*)s(1),s(npxt),Gb(1),Gb(npxt)pauseelseend if!=!=计算悬移质推移质及年均流量=!=WQS1=Q(1)*S(1)WQSi=Q(npxt)*S(npxt)dWQSt=dWQST+WQS1*dt/10000000 ! 为累计悬移质输沙量,单位为万吨dWQS=dWQS+(WQS1-WQSi)*dt/(1.325*1000)/1E04!
18、为累计悬移质冲淤量,单位为万立方米dWGb=dwGb+(Gb(1)-Gb(npxt)*dt/(1.325*1000)/1E04! 为累计推移质淤积总量,单位为万立方米AVQ=AVQ+Q(1)*dt/10000 ! 为累计年径流量 单位万方dWQ=dWQS+dWGbwsin=wsin+Gb(1)*dt/10000000!wsout=wsout+WQSi*dt/10000000!=!=水下断面修改断面点高程=!=do i=1,npxtdo j=1,Npoint(i)if(x(j,2,i)=zlevel(i) thenx(j,2,i)=zlevel(i)endifendifenddoenddoQQ(
19、k1)=QQ(k1)+Q(1)QQ(0)=0!=!=输出年淤积总量=!=do lk=1,nyif(k1=Nday(lk) thenif(lk1)thenQa=(QQ(k1)-QQ(k1-Nday(lk-1)/(Nday(lk)-Nday(lk-1)else Qa=QQ(k1)/Nday(lk)endifwrite(15,*) lk+10*(k-1),dWQ,avq,dWQSt+wsin,Qawrite(*,*) lk+10*(k-1),QaAVQ=0.0endifenddo!=!=十年一输出的水面线变化=!=十年一输出的深泓变化=!=十年一输出坝前断面变化=!=十年一输出水位库容关系=!=ca
20、ll level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,fai1,ndisp)if(k1=nday(10)thendo i=1,npxtwrite(13,*) 10*k,i,dxa(i)/1000.0,zlevel(i)write(16,*)10*k,i,alow(i)enddodo j=1,npoint(npxt)write(17,*)10*k,x(j,1,npxt),x(j,2,npxt)enddodo Z0=275,250,-1zL=Z0*1.0call FZV(NPOINT,npxt,X,zl,ZVV,dx,ndisp,MM,NN)
21、write(18,*)10*k,zl,zvvenddoendifenddoclose(11)enddoend!=!=计算断面面积程序=!= SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)! THIS SUBROUTINE IS USED TO CACULATE AREA FOR IRREGULAR CROSS-SECTION DIMENSION X(NPOINT),Z(NPOINT) IF(NDISP.LE.-1)WRITE(*,*)INTO AREA IF(NPOINT.LE.2)THEN WRITE(*,*)IN AREA THE INPUT
22、 DATA ARE FALSE N=,Npoint STOP ENDIF IF(Z(1).LT.ZLEVEL.OR.Z(Npoint).LT.ZLEVEL)THEN WRITE(*,*)IN AREA THE WATER LEVEL IS TOO HIGH OR THE HEIGHT& OF THE SECTION AT EDGE IS TOO LOW STOP ENDIF A=0. B=0. XW=0. DO 10 I=1,NPOINT-1 ZMIN=AMIN1(Z(I),Z(I+1) IF(ZMIN.LT.ZLEVEL)THEN ZMAX=AMAX1(Z(I),Z(I+1) IF(ZMAX
23、.LT.ZLEVEL)THEN DB=X(I+1)-X(I) DH=ZLEVEL-0.5*(Z(I)+Z(I+1) DX=(Z(I)-Z(I+1)*2.+DB*DB)*0.5 ELSE DB=(ZLEVEL-ZMIN)/(ZMAX-ZMIN)*(X(I+1)-X(I) DH=0.5*(ZLEVEL-ZMIN) DX=(2*DH)*2.0+DB*DB)*0.5 ENDIF IF(DB.LT.0.0)THEN WRITE(*,*)THE DISTANCE FOR NODE ,I,AND,I+1,ARE FALSE WRITE(*,*)IN AREA,I,I+1,X(I),X(I+1) STOP E
24、NDIF DA=DB*DH B=B+DB A=A+DA XW=XW+DX ENDIF10 CONTINUE return END!=!=恒定非均匀缓流水面线计算程序=!= subroutine level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,fai,NDISP) dimension x(mm,2,nn),rough(npxt),dx(npxt),zlevel(npxt),& q(npxt),npoint(npxt),b(npxt),a(npxt),xw(npxt) IF(NDISP.LE.-1)WRITE(*,*)INTO LEVEL nc
25、=100 dz=0.1 dz1=0.03 dz2=0.1 zlevel(npxt)=267.0!write(*,*)npoint(npxt),x(1,1,npxt),x(1,2,npxt),b(npxt),&!a(npxt),xw(npxt),zlevel(npxt),ndisp !pause call area(npoint(npxt),x(1,1,npxt),x(1,2,npxt),b(npxt),& a(npxt),xw(npxt),zlevel(npxt),ndisp)do 10 ip=npxt-1,1,-1 NR=0 ! write(*,*)ip=,ip zmin=zlevel(ip
26、+1)+dz ! WRITE(*,*)ZMIN1,ZMIN30 call area(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),& a(ip),xw(npxt),zmin,ndisp) if(a(ip).LE.0)then zmin=zmin+0.173 goto 30 endif fmin=flevel(zlevel(ip+1),zmin,dx(ip),q(ip+1),q(ip),&rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai1,ndisp) !if(ip.eq.3.and.kl.eq.89)then!write(*,*)n
27、point(4),b(4),a(4),xw(4),zlevel(4),ndisp!write(*,*)(x(ikj,2,3),ikj=1,npoint(3)!write(*,*)fmin,zlevel(ip+1),zmin,dx(ip),q(ip+1),q(ip)!write(*,*)rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai,ndisp!pause!endif if(fmin.gt.0)then fr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip) if(fr.lt.1)then zmax=zmin+dz20 call area(n
28、point(ip),x(1,1,ip),x(1,2,ip),b(ip),& a(ip),xw(ip),zmax,ndisp) fmax=flevel(zlevel(ip+1),zmax,dx(ip),q(ip+1),q(ip),& rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),fai1,ndisp) if(fmax.gt.0.0)then zmax=zmax+dz goto 20 endif else zmin=zmin+dz2 NR=NR+1 IF(NR.GT.NC)THEN zmin=zmin-dz1 write(*,*)the loop is death,p
29、ause WRITE(*,*),nr,ip,zmin,fmin,nr,ip,zmin,fmin read(*,*) endif goto 30 endif else zmin=zmin-dz1 goto 30 endif call bisec(zmin,zmax,fmin,zlevel(ip),npoint(ip),X(1,1,ip),& x(1,2,ip),b(ip),a(ip),xw(ip),q(ip),ndisp,zlevel(ip+1),b(ip+1),& a(ip+1),xw(ip+1),q(ip+1),dx(ip),rough(ip),fai) 10 continue return end !=二分法求根程序= subroutine bisec(zmin,zmax,fmin,zlevel,npoint,x,z,b,a,xw,q,& ndisp,zlelo,blower,alower,xwlower,qlower,dx,rough,fai) dimension x(npoint),z(npoint) IF(NDISP.EQ.-1)WRITE(*,*)INTO BISEC ERR=0.00110 ddz=zmax-zmin z
