河流模擬課程設(shè)計(jì)

1、水庫(kù)一維泥沙淤積計(jì)算課程設(shè)計(jì)目錄一、目的與要求1二、基本原理11、基本方程12、方程離散13、公式補(bǔ)充2三、計(jì)算步驟2四、計(jì)算框圖3五、計(jì)算結(jié)果41、歷年輸沙量特征值42、各年淤積總量53、各年水位庫(kù)容關(guān)系64、水面線的變化75、深泓變化86、壩前斷面變化9六、結(jié)果分析121、剖面形態(tài)分析122、庫(kù)容損失合理性分析12七、計(jì)算程序13一、目的與要求通過(guò)課程設(shè)計(jì)，初步掌握一維數(shù)學(xué)模型建立數(shù)學(xué)模型的基本過(guò)程和計(jì)算方法，具備一定的解決實(shí)際問(wèn)題的能力。以水流、泥沙方程為基礎(chǔ)，構(gòu)建恒定流條件下的河道一維水沙數(shù)學(xué)模型，并編制出完整的計(jì)算程序，并以某個(gè)水庫(kù)為實(shí)例，進(jìn)行水庫(kù)泥沙淤積計(jì)算。水流條件：恒定非均勻流

2、。泥沙條件：包括懸移質(zhì)，推移質(zhì)的均勻沙模型，推移質(zhì)計(jì)算模式為飽和輸沙，懸移質(zhì)計(jì)算模式為不飽和輸沙，水流泥沙方程采用非耦合解?；驹?、基本方程水流連續(xù)方程:Q2A水流運(yùn)動(dòng)方程hgAgAioifx或QQ2gAZ22Qn0txAg4AR3泥沙連續(xù)方程SAQS(SS*)tx河床變形方程y。0tGbx(SS*)推移質(zhì)平衡輸沙方程Mayer-_Peter公式,G=G水流挾沙力公式采用張瑞瑾公式，推移質(zhì)輸沙率公式采用MAYER-PETER公式中的能坡J按均勻流曼寧公式近似計(jì)算（每個(gè)斷面不同）2、方程離散方程在恒定流情況下有&0，離散為：Q=constx方程變形為Q2n2A2R43Q2-72x2g

3、A22Qn24AR3上式離散為2gQj2A2Zj1zj2xnQj4-R3j1乎4(1AjRj3方程S（4）去掉時(shí)間項(xiàng)得到Sj1Sj1SjS*jexpS*jSji-1expxBGbx對(duì)2號(hào)斷面以下，上式可以離散為:由方程（4-5）可得QS0By。BGbBGbQSQSxxBjy°j(10)Bj1y°j1t對(duì)于進(jìn)口斷面，推移質(zhì)不考慮,懸移質(zhì)采用單點(diǎn)離散（SS）q該方程的解析解為：方程（5）可離散為:y。(S1Sm)t3、公式補(bǔ)充3kgR0.124,m取1.05，干密度o取1.3恢復(fù)飽和系數(shù)0.25均勻沙粒徑為d=0.041mm（懸移質(zhì)），d=2mm（推移質(zhì)）二、計(jì)算步驟1、輸入河

4、床地形糙率等數(shù)據(jù)求得斷面面積與水位的關(guān)系（AZ）,進(jìn)而求得斷面平均流速uQ，水力學(xué)半徑RAA2、讀入一個(gè)時(shí)段的水沙數(shù)據(jù)(特別注意，不要一次性將數(shù)據(jù)全部讀入)讀入第一時(shí)段(Q,S)值3、計(jì)算水面線，同時(shí)得到各斷面的水力要素求得各個(gè)斷面的河寬、斷面面積、水深、平均流速等值計(jì)算前要注意在壩前輸入水位，各斷面均應(yīng)對(duì)流量賦值4、計(jì)算懸移質(zhì)水流挾沙力3mSk-gRK取0.124,m取1.05。5、計(jì)算推移質(zhì)輸沙率(采用mayer-peter公式)3/2n3/2()hJ0.047(s)dngb0.1251/2()gs6、計(jì)算各斷面含沙量公式Sj1S*j1SjSjexpLS*jSj11expqx7、計(jì)算各斷面

5、沖淤厚度進(jìn)口斷面y0一(Si&)t08、修改水各斷面水下河床高程HH°Y9、重新進(jìn)入(2)進(jìn)行下一循環(huán)10、計(jì)算10年河床變形，計(jì)算時(shí)段為一天，單位為秒(s)11、淤積總量年輸出一次，其余每?jī)赡贻敵鲆淮斡?jì)算結(jié)果四、計(jì)算框圖讀入地形資料和糙率i=i+1計(jì)算水流子程序計(jì)算S和Gb計(jì)算S(j)和dy(j)Y-J<npnxt重新計(jì)算S(j)修改河床地形五、計(jì)算結(jié)果1。（其中推移質(zhì)輸沙量是按照1、歷年輸沙量特征值計(jì)算各年懸疑質(zhì)輸沙量、推移質(zhì)輸沙量、年均流量見表懸疑質(zhì)輸沙量的1.5%計(jì)算的）。表1龍開口壩址歷年流量輸沙量特征值代表年年平均懸移質(zhì)輸沙量（萬(wàn)t）年平均推移質(zhì)輸沙量（萬(wàn)t

6、）年平均流量（萬(wàn)t）1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900取大值645696.841900最小值198229.731290多年平均值428764.3116402、各年淤積總量表2歷年淤積總量的變化年懸移質(zhì)輸沙量（萬(wàn)噸）年均流量（萬(wàn)噸）總淤積量（億立方米）000011821.494845560.14

7、25291.20510994180.387316993.5419658070.879420027.3424697751.085521668.3327656171.204623853.9830875691.368730155.0537460621.701834879.5943121571.987936754.1846403382.1171039289.6650197622.243LA2515aTK-210漫積最億立方米024b圖1歷年淤積總量的變化bio岸拍3、各年水位庫(kù)容關(guān)系表3各年水位庫(kù)容關(guān)系水位（g0年2年4年6年8年10年2732.5742.2071.5721.3130.7650.569

8、2712.3852.0181.3841.1320.6160.4382692.2021.8351.2020.980.4930.3332672.0241.6611.0440.8590.3980.262651.8521.510.9280.7730.3460.232631.6881.3720.830.6980.3050.2042611.5321.2450.7450.6280.2690.1792591.3841.1260.6680.5630.2370.1562571.2471.0140.5940.5010.2070.1332551.1170.9070.5230.4420.1790.1122530.995

9、0.8080.4570.3870.1540.0922510.880.7170.3990.3340.1310.0742490.7710.6320.3440.2870.110.0592470.670.5530.2930.2440.0890.0452450.5760.4780.2440.2040.0710.0342430.4920.4080.1990.1680.0530.0252410.4170.3470.1590.1370.0380.017Jj75p"0511.522.5容（10"nnO圖2各年水位庫(kù)容關(guān)系曲線4、水面線的變化表4水面線變化年距壩里程(m)水位(m)年距壩里程(

10、m)水位(m)年距壩里程(m)水位(m)034920267.6434920267.88834920274.43033810267.16433810267.7683381027401431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267424050267.07824050271.5021846267421846267.02821846270.53019896267419896267819896269.580182702674182

11、70267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距壩里程(m)水位(m)年距壩里程(m)水位(m)年距壩里程(m)水位(m)234920267011034920273.99233810267.12633810271.751033810

12、273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.391019896271.07218270267618270267.111018270270.482157902676157902671015790269.721317026761317026710

13、13170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267圖3水面線變化曲線5、深泓變化表5逐年深泓變化年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)034920264.1434920264.77834920271.75033810261.7433810264.12833810270.46031350253.3431350262.29831350268

14、.95029280249.4429280261.92829280268.6027000246.6427000263.23827000269.21024050241.3424050261.99824050267.71021846239.2421846261.5821846267.87019896237.9419896261.498198962681418270252.83818270266.84015790230.1415790250.82815790264.62013170228.2413170240.94813170264.74010170224.441017023

15、9.22810170262.5907820222.147820230.9887820253.7604950219.844950228.2984950237.2802100215.9242100219.9782100233.9900213.1740221.0480224.22年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)年距壩里程(m)深泓(m)234920262.79634920269.961034920271.76233810262.41633810269.081033810270.93231350259.92631350267.251031350269.37229280259.536

16、29280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242941019896269.7218270242.84618270264.711018270268.61215790233.11615790258.711015790267.03213170232.41613170249.521013170266.88210170226.496101702

17、42.18101017026527820224.1167820234.07107820264.7924950220.6564950228.12104950245.6322100216.8662100222.71102100236.5420213.7860219.7100229.24U年澇苴T-的謙泓年澇眼十6年避泓8年配泓.10年滁第10年水司凌2001CQD0200QQ000040000距以里程(m)圖4深泓變化曲線6、壩前斷面變化該水庫(kù)為三角形淤積，最靠近壩前斷面淤積并不會(huì)特別嚴(yán)重，為使斷面變化情況表現(xiàn)的更加直觀，這里采用IP-17斷面作為典型斷面并分析其淤積情況。表6壩前斷面變化年起點(diǎn)距

18、(m)高程(m)年起點(diǎn)距(m)高程(m)年起點(diǎn)距(m)高程(m)00300.220300.240300.2010.5272.29210.5272.29410.5272.29015.8266.93215.8267415.8267.01018.2260.24218.2261.18418.2264.29018.9259.62218.9260.56418.9263.67019.4259.88219.4260.82419.4263.93024.8259.91224.8260.85424.8263.96026.4259.82226.4260.76426.4263.87044.2246.15244.2247

19、.09444.2250.2061.2234.47261.2235.41461.2238.52078.6218.17278.6219.11478.6222.22091.5215.92291.5216.86491.5219.97097.3216.1297.3217.04497.3220.150106.8216.982106.8217.924106.8221.030121.6217.472121.6218.414121.6221.520132.4218.052132.4218.994132.4222.10145.6218.452145.6219.394145.6222.50173.8218.9621

20、73.8219.94173.8223.010191.7219.432191.7220.374191.7223.480203.6230.782203.6231.724203.6234.830216.5242.882216.5243.824216.5246.930228.7255.042228.7255.984228.7259.090241.1267.442241.1267.444241.1267.440253.7279.812253.7279.814253.7279.810265.3291.22265.3291.24265.3291.20279.1304.522279.1304.524279.1

21、304.520293.1300.32293.1300.34293.1300.3年起點(diǎn)距(m)高程(m)年起點(diǎn)距(m)高程(m)年起點(diǎn)距(m)高程(m)60300.280300.2100300.2610.5272.29810.5272.291010.5272.29615.8267.01815.8267.011015.8267.01618.2267818.22671018.2267618.9266.41818.92671018.9267619.4266.67819.42671019.4267624.8266.7824.82671024.8267626.4266.61826.42671026.426

22、7644.2252.94844.2264.221044.2267661.2241.26861.2252.541061.2255.09678.6224.96878.6236.241078.6238.79691.5222.71891.5233.991091.5236.54697.3222.89897.3234.171097.3236.726106.8223.778106.8235.0510106.8237.66121.6224.268121.6235.5410121.6238.096132.4224.848132.4236.1210132.4238.676145.6225.248145.6236.

23、5210145.6239.076173.8225.758173.8237.0310173.8239.586191.7226.228191.7237.510191.7240.056203.6237.578203.6248.8510203.6251.46216.5249.678216.5260.9510216.5263.56228.7261.838228.726710228.72676241.1267.448241.1267.4410241.1267.446253.7279.818253.7279.8110253.7279.816265.3291.28265.3291.210265.3291.26

24、279.1304.528279.1304.5210279.1304.526293.1300.38293.1300.310293.1300.3220320高3UO程m24(*200050100150200250300孟點(diǎn)也W)圖5壩前斷面變化曲線六、結(jié)果分析1、剖面形態(tài)分析A、縱斷面淤積情況分析：由深泓變化曲線圖可以看出，隨著時(shí)間的增加，河道深泓點(diǎn)高程成增高趨勢(shì)。該水庫(kù)運(yùn)行10年后，在距壩約10000米處，深泓點(diǎn)高程增加40多米，計(jì)算結(jié)果表明，該水庫(kù)泥沙以淤積為主，致使水庫(kù)庫(kù)容逐年損失，且情況較為嚴(yán)重。8壩前斷面變化曲線分析：由壩前斷面變化曲線圖可以看出，隨著時(shí)間的增加，壩前斷面各點(diǎn)的總體高程也

25、成增高趨勢(shì)。該水庫(kù)運(yùn)行10后，該斷面河床抬高20多米，表明在該處河道逐年以淤積為主，且淤積情況較為嚴(yán)重。2、庫(kù)容損失合理性分析由水位庫(kù)容曲線可以看出，水庫(kù)在運(yùn)行10年后，經(jīng)計(jì)算可知，在273m水位條件下，庫(kù)容的損失率為77.89%,在241m的水位條件下，庫(kù)容的損失率高達(dá)95.92%,計(jì)算結(jié)果表明，經(jīng)過(guò)十年的運(yùn)行，水庫(kù)庫(kù)容基本成全部淤滿狀態(tài)，所以在10年后將不能再繼續(xù)使用。結(jié)果表明，該水庫(kù)的調(diào)節(jié)能力較差。上述計(jì)算結(jié)果與現(xiàn)實(shí)不太相符，其原因是在計(jì)算該水庫(kù)庫(kù)容的泥沙損失時(shí)，計(jì)算沖往下游的泥沙量不合理。在給定條件下，該水庫(kù)并沒(méi)有考慮水位的變化及壩前沖沙情況。該水庫(kù)的計(jì)算條件為正常蓄水位保持不變，在高

26、水位運(yùn)行條件下，必然會(huì)導(dǎo)致水庫(kù)的嚴(yán)重淤積，而實(shí)際情況是，在洪水來(lái)臨時(shí)，壩前水位要降至防洪限制水位，以增大防洪庫(kù)容，該條件下，水流的速度及水流挾沙力較大，會(huì)造成庫(kù)區(qū)泥沙處于沖刷狀態(tài)，以增加水庫(kù)庫(kù)容。七、計(jì)算程序programmainparameter(nn=16,mm=60,ndisp=0,npxt=16,ND=3653,NY=10,T=15)paramete(rs=26.5,r=10.0,d=0.000041,failev=0.5,faidy=0.5,dt=86400)dimensionx(mm,2,nn),rough(nn),npoint(nn),B(nn),A(nn),Gb(nn),dy(

27、nn),Nday(NY)dimensiondx(nn),Xw(nn),dxa(nn),alow(nn),alow0(nn),Q(nn),Sx(nn),S(nn),Zlevel(nn)data(Nday(m),m=1,10)/366,731,1096,1461,1827,2192,2557,2922,3288,3653/!*打開輸出文件*open(10,file='河床地形.txt',status='old')open(11,file='流量沙量.txt',status='old')open(12,file='逐年累計(jì)淤積量

28、.txt',status='unknown')open(13,file='逐年深泓變化.txt',status='unknown')open(14,file='逐年水面線變化.txt',status='unknown')open(15,file='逐年壩前斷面淤積變化.txt',status='unknown')open(16,file='逐年水位庫(kù)容關(guān)系變化.txt',status='unknown')!*讀入地形數(shù)據(jù)*callQbed(mm

29、,npxt,x,dxa,alow0,npoint)*計(jì)算沉降速度*W=FW(T,r,rs,d,0)write(*,*)W!*計(jì)算斷面間距*doi=1,npxt-1dx(i)=dxa(i)-dxa(i+1)enddo!*計(jì)算初始水面線、深泓、庫(kù)容水位關(guān)系*dWQST=0.0dWQS=0.0dWGb=0.0AVQ=0.0doI=1,NPXTrough(I)=0.035enddoQ(1)=295.doI=2,NPXTQ(I)=Q(I-1)enddoZlevel(npxt)=267.0CALLlevel(x,rough,nn,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev

30、,NDISP)!初始水面線!*輸出到文件*CALLFILE(12,13,14,15,16,0,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dWQST,AVQ,dWQS,dWGb,ndisp)!*計(jì)算各個(gè)時(shí)段水庫(kù)的淤積情況*DO40,I=1,NDWRITE(*,*)'讀入天數(shù):',I!*讀入流量、沙量數(shù)據(jù)*READ(11,*)Q(1),S(1)IF(Q(1)=0)THENGOTO41ENDIFDOJ=2,NPXTQ(J)=Q(1)ENDDO!*CALLlevel(x,rough,nn,zlevel,dx,Q,npoint,b,a,xw,nn,mm,fail

31、ev,0)!*計(jì)算各斷面的水流挾沙力和推移質(zhì)輸沙率*CALLSGB(nn,a,b,Q,Sx,Gb,rough,w,ndisp)!*計(jì)算各斷面的含沙量和沖淤厚度*CALLDEPO(nn,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alow0,npoint,faidy,ndisp)!*檢驗(yàn)質(zhì)量守恒*CALLCHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp)41CONTINUE!*計(jì)算懸移質(zhì)、推移質(zhì)年均流量*WQS1=Q(1)*S(1)Gb(1)=WQS1*0.015WQSN=Q(nn)*S(nn)dWQST=dWQST+WQ

32、S1*dt/10000000!懸移質(zhì)累計(jì)輸沙量dWQS=dWQS+(WQS1-WQSN)*dt/1300.0/100000000!懸移質(zhì)累計(jì)沖淤量dWGb=dWGb+(Gb(1)-Gb(nn)*dt/1300.0/100000000!推移質(zhì)累計(jì)淤積量AVQ=AVQ+Q(1)*dt/10000!累計(jì)流量!*輸出相關(guān)文件*DOm=1,NYIF(I=Nday(m)thenCALLFILE(12,13,14,15,16,m,x,npoint,Zlevel,alow,dxa,dx,mm,nn,dWQST,AVQ,dWQS,dWGb,ndisp)ENDIFENDDO40CONTINUEclose(10)c

33、lose(11)close(12)close(13)close(14)close(15)close(16)end!*讀入地形數(shù)據(jù)*subroutineQbed(mm,nn,x,dxa,alow0,npoint)dimensionx(mm,2,nn),dxa(nn),alow0(nn),npoint(nn)doi=1,nnread(10,*)read(10,*)dxa(i),npoint(i)read(10,*)doj=1,npoint(i)read(10,*)x(j,1,i),x(j,2,i)enddo!*計(jì)算深泓*alow0(i)=x(1,2,i)dom=2,npoint(i)if(alow

34、0(i)>x(m,2,i)thenalow0(i)=x(m,2,i)endifenddoenddoend!*泥沙沉降速度的計(jì)算，采用張瑞瑾公式*FUNCTIONFW(T,gama,gamas,d,ndisp)IF(NDISPEQ.-1)WRITE(*,*)'INTOFWcallVISCOS(T,CMU)A1=CMU*13.95/dA2=1.09*9.8*d*(gamas-gama)/gamaFW=(A1*2.0+A2)*0.5-A1RETURNENDSUBROUTINEVISCOS(T,CMU)X=1.775E-06A=1+0.0337*T+0.000221*T*TCMU=X/A

35、RETURNEND!*水面線及各斷面水力要素計(jì)算函數(shù)*subroutinelevel(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP)dimensionx(mm,2,nn),rough(npxt),dx(npxt),zlevel(npxt)dimensionq(npxt),npoint(npxt),b(npxt),a(npxt),xw(npxt)IF(NDISP=-1)WRITE(*,*)'INTOLEVEL'nc=1000dz=0.5dz1=0.1317dz2=2.079callarea(npoint(npx

36、t),x(1,1,npxt),x(1,2,npxt),b(npxt),a(npxt),xw(npxt),zlevel(npxt),ndisp)doip=npxt-1,1,-1NR=0zmin=zlevel(ip+1)+dz30callarea(npoint(ip),x(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(npxt),zmin,ndisp)if(A(ip)<=0)thenzmin=zmin+0.173goto30endiffmin=flevel(zlevel(ip+1),zmin,dx(ip),q(ip+1),q(ip),rough(ip),b(ip+1),b

37、(ip),a(ip+1),a(ip),failev,ndisp)if(fmin>0)thenfr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip)if(fr<1)thenzmax=zmin+dz!*計(jì)算斷面要素*20callarea(npoint(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),failev,ndisp)if(fm

38、ax>0)thenzmax=zmax+dzgoto20endifelsezmin=zmin+dz2NR=NR+1IF(NR>NC)THENzmin=zmin-dz1write(*,*)'theloopisdeath,pause'WRITE(*,*)'nr,ip,zmin,fmin',nr,ip,zmin,fminread(*,*)endifgoto30endifelsezmin=zmin-dz1goto30endifcallbisec(zmin,zmax,fmin,zlevel(ip),npoint(ip),X(1,1,ip),x(1,2,ip),b

39、(ip),a(ip),xw(ip),q(ip),ndisp,zlevel(ip+1),b(ip+1),a(ip+1),q(ip+1),dx(ip),rough(ip),failev)enddoreturnendFUNCTIONflevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev,ndisp)IF(NDISPEQ.-1)WRITE(*,*)'INTOFLEVEL'HLOWER=ALOWER/BLOWERH=A/BAA=ZLELO-ZLEVELB1=failev*(QLOWER/BLOWER)*2.0/HLOW

40、ER*(10/3.0)B2=(1-failev)*(Q/B)*2.0/H*(10/3.0)BB=DX*ROUGH*2.0*(B1+B2)C1=(QLOWER/ALOWER)*2.0C2=(Q/A)*2.0CC=(C1-C2)/(2*9.8)FLEVEL=AA+BB+CCRETURNEND!*計(jì)算斷面要素*SUBROUTINEAREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)DIMENSIONX(NPOINT),Z(NPOINT)IF(NDISP<=-1)WRITE(*,*)'INTOAREA'IF(NPOINT<=2)THENWRITE(*,*

41、)'INAREATHEINPUTDATAAREFALSEN=',NpointSTOPENDIFIF(Z(1)<ZLEVEL.OR.Z(NPOINT)<ZLEVEL)THENWRITE(*,*)'INAREATHEWATERLEVELISTOOHIGHORTHEHEIGHTOFTHESECTIONATEDGEISTOOLOWSTOPENDIFA=0.0B=0.0XW=0.0DOI=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1)IF(ZMIN<ZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1)IF(ZMAX<ZLEV

42、EL)THENDB=X(I+1)-X(I)DH=ZLEVEL-0.5*(Z(I)+Z(I+1)DX=(Z(I)-Z(I+1)*2.0+DB*DB)*0.5ELSEDB=(ZLEVEL-ZMIN)/(ZMAX-ZMIN)*(X(I+1)-X(I)DH=0.5*(ZLEVEL-ZMIN)DX=(2*DH)*2.0+DB*DB)*0.5ENDIFIF(DB<0.0)THENWRITE(*,*)'THEDISTANCEFORNODE',I,'AND',I+1,'AREFALSE'WRITE(*,*)'INAREA',I,I+1,X(

43、I),X(I+1)STOPENDIFDA=DB*DHB=B+DBA=A+DAXW=XW+DXENDIFENDDORETURNEND!*二分法計(jì)算*subroutinebisec(zmin,zmax,fmin,zlevel,npoint,x,z,b,a,xw,q,ndisp,zlelo,blower,alower,qlower,dx,rough,failev)dimensionx(npoint),z(npoint)IF(NDISPEQ.-1)WRITE(*,*)'INTOBISEC'ERR=0.00110ddz=zmax-zminzlevel=0.5*(zmin+zmax)cal

44、larea(npoint,x,z,b,a,xw,zlevel,ndisp)f=flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev,ndisp)if(ddz>err)thenff=f*fminif(ff>0)thenzmin=zlevelelsezmax=zlevelendifddz=zmax-zmingoto10endifreturnend!懸移質(zhì)挾沙力和推移質(zhì)輸沙率計(jì)算，懸移質(zhì)采用張瑞謹(jǐn)挾沙力公式!推移質(zhì)采用mayer-peter公式SUBROUTINESGB(nn,a,b,Q,Sx,Gb,rough,w,

45、ndisp)DIMENSIONA(nn),B(nn),Q(nn),Sx(nn),Gb(nn),rough(nn)IF(ndisp=-1)WRITE(*,*)'INTOSGB'DOI=1,nnH=A(I)/B(I)U=Q(I)/A(I)VA1=U*3.0VA2=9.8*H*WVA3=(VA1VA2)*1.05Sx(I)=VA3*0.124!懸移質(zhì)挾沙力!推移質(zhì)輸沙率，斷面輸沙率Gb(I)=FGb(rough(I),10.0,H,Q(I),A(I),26.5,0.002,1.0,2.65,0)*B(i)ENDDOENDFUNCTIONFGb(rn,r,H,Q,A,rs,d,rou,

46、rous,ndisp)IF(ndisp=-1)WRITE(*,*)'INTOFGb'A1=(rn*2.0)*(Q*2.0)A2=(A*2.0)*(H*(4/3.0)RJ=A/A2rnn=(d*(1/6.0)/26.0AA=(rnn/rn)*1.5)*rs*H*RJBB=0.047*(rs-r)*dCC=(rous-rou)/rous)*9.8*(rou*0.5)*0.125IF(AA-BB)<0)THENFGB=0ELSEFGB=(AA-BB)*1.5)/CCENDIFRETURNEND!*斷面沖淤厚度計(jì)算*SUBROUTINEDEPO(nn,mm,w,x,S,Sx,Gb

47、,dt,B,Q,dx,dy,alow,alow0,npoint,faidy,ndisp)DIMENSIONx(mm,2,nn),npoint(nn),B(nn),Gb(nn),dy(nn),dx(nn),alow(nn),alow0(nn),Q(nn),Sx(nn),S(nn)IF(ndisp=-1)WRITE(*,*)'INTODEPO'rou0=1300.0dy(1)=0.25*w*(S(1)-Sx(1)*dt/rou0!進(jìn)口斷面河床變形方程單點(diǎn)離散!*上游第二斷面開始離散求沖淤量*DOJ=2,NNqq=Q(J-1)/(B(J-1)+B(J)*0.5)FAC=-(0.25*

48、w*dx(j-1)/qqSA=(S(J-1)-Sx(J-1)*EXP(FAC)SB=(Sx(J-1)-Sx(J)*(1-EXP(FAC)/FACS(J)=Sx(J)+SA-SBdGb=(Gb(j)-Gb(j-1)/dx(j-1)dQS=(Q(J)*S(J)-Q(J-1)*S(J-1)/dx(j-1)dy1=-(dGb+dQS)*dt/rou0dy2=-B(J-1)*dy(J-1)*(1-faidy)dy3=B(J)*faidydy(J)=(dy1+dy2)/dy3!*判斷庫(kù)尾沖刷問(wèn)題*alow(J)=x(1,2,J)DOL=2,npoint(J)if(alow(J)>x(L,2,J)al

49、ow(J)=x(L,2,J)ENDDOIF(alow(J)+dy(J)<alow0(J)THENdy(J)=0temp=rou0*B(j-1)*dy(J-1)*(1-faidy)/dtS(J)=-(temp+dGb)*dx(J-1)/Q(J-1)+S(J-1)ENDIFENDDOEND!*檢驗(yàn)?zāi)嗌呈欠袷睾?SUBROUTINECHECK(dy,B,S,Gb,dt,npxt,Zlevel,x,npoint,Q,dx,ndisp)DIMENSIONx(60,2,npxt),npoint(npxt),B(npxt),Gb(npxt),dy(npxt),dx(npxt),Q(npxt),S(npxt),Zlevel(npxt)IF(ndisp=-1)WRITE(*,*)'INTOCHECK'DVV=0.0DOM=2,npxtDVV=DVV+(B(M-1)*dy(M-1)+B(M)*dy(M)*dx(M-1)/2.0TOL仁DVV*1300TOL2=(Q(1)*(S(1)-S(M