計算水力學基礎

--本頁僅作為文檔封面，使用時請直接刪除即可--

--內頁可以根據需求調整合適字體及大小本頁僅作為文檔封面，使用時請直接刪除即可--

--內頁可以根據需求調整合適字體及大小--計算水力學基礎(總79頁)PAGE計算水力學基礎李占松編著鄭州大學水利與環(huán)境學院

內容簡介本講義是編者根據多年的教學實踐，并參考《微機計算水力學》（楊景芳編著，大連理工大學出版社出版，1991年5月第1版）等類似教材，取其精華，編寫而成的。目的是使讀者掌握通過計算機解水力學問題的方法，為解決更復雜的實際工程問題打下牢固的計算基礎。書中內容包括：數值計算基礎，偏微分方程式的差分解法，有限單元法；用這些方法解有壓管流、明渠流、閘孔出流、堰流、消能、地下水的滲流及平面勢流等計算問題。講義中的用FORTRAN77算法語言編寫的計算程序，幾乎包括了全部水力學的主要計算問題。另外，結合講授對象的實際情況，也提供了用VB算法語言編寫的計算程序。VB程序編程人員的話為了更好地促進水利水電工程建筑專業(yè)的同學學好《微機計算水力學》這門學科，編程員借暑假休息的時間，利用我們專業(yè)目前所學的VB中的算法語言部分對水力學常見的計算題型編制成常用程序。希望大家能借此資料更好地學習《微機計算水力學》這門課程。本程序著重程序的可讀性，不苛求程序的過分技巧。對水力學中常用的計算題型，用我們現在所學的VB語言編制而成。由于編程員能力有限，程序中缺點和錯誤在所難免，望老師和同學及時給予批評指正。VB程序編程人員：黃渝桂曹命凱

前言計算水力學的形成與發(fā)展計算水力學作為一門新學科，形成于20世紀60年代中期。水力學問題中有比較復雜的紊流、分離、氣穴、水擊等流動現象，并存在各種界面形式，如自由水面、分層流、交界面等。由各種流動現象而建立的數學模型（由微分方程表示的定解問題），例如連續(xù)方程、動量方程等組成的控制微分方程組，多具有非線性和非恒定性，只有少數特定條件下的問題，可根據求解問題的特性對方程和邊界條件作相應簡化，而得到其解析解。因此長期以來，水力學的發(fā)展只得主要藉助于物理模型試驗。隨著電子計算機和現代計算技術的發(fā)展，數值計算已逐漸成為一個重要的研究手段，發(fā)展至今，已廣泛應用與水利、航運、海洋、流體機械與流體工程等各種技術科學領域。計算水力學的特點是適應性強、應用面廣。首先流動問題的控制方程一般是非線性的，自變量多，計算域的幾何形狀任意，邊界條件復雜，對這些無法求得解析解的問題，用數值解則能很好的滿足工程需要；其次可利用計算機進行各種數值試驗，例如，可選擇不同的流動參數進行試驗，可進行物理方程中各項的有效性和敏感性試驗，以便進行各種近似處理等。它不受物理模型試驗模型律的限制，比較省時省錢，有較多的靈活性。但數值計算一是依賴于基本方程的可靠性，且最終結果不能提供任何形式的解析表達式，只是有限個離散點上的數值解，并有一定的計算誤差；二是它不像物理模型試驗一開始就能給出流動現象并定性地描述，卻往往需要由原體觀測或物理實驗提供某些流動參數，并對建立的數學模型驗證；三是程序的編制及資料的收集、整理與正確利用，在很大程度上依賴于經驗與技巧。所以計算水力學有自己的原理方法和特點，數值計算與理論分析觀測和試驗相互聯系、促進又不能相互代替，已成為目前解決復雜水流問題的主要手段之一，尤其是在研究流動過程物理機制時，更需要三者有機結合而互相取長補短。近三、四十年來，計算水力學有很大的發(fā)展，替代了經典水力學中的一些近似計算法和圖解法。例如水面曲線計算；管網和渠系的過水或輸沙（排污）能力的計算；有水輪機負荷改變時水力震蕩系統的穩(wěn)定性計算研究；流體機械過流部件的流道計算以及優(yōu)化設計，還有洪水波、河口潮流計算，以及各種流動條件下，不同排放形式的污染物混合計算等。上世紀70年代中期已從針對個別工程問題建立的單一數學模型，開始建立對整個流域洪泛區(qū)已建或規(guī)劃中的水利水電工程進行系統模擬的系統模型。理論課題的研究中，對擴散問題、傳熱問題、邊界層問題、漩渦運動、紊流等問題的研究也有了很大的發(fā)展，并已開始計算非恒定的三維紊流問題。由于離散的基本原理不同，計算水力學可分為兩個分支：一是有限差分法，在此基礎上發(fā)展的有有限分析法；二是有限單元法，在此基礎上提出了邊界元法和混合元法，另外還有迎風有限元法等。

目錄數值計算基礎非線性方程式的解法迭代法牛頓-拉普森切線法二等分法線性方程組的解法高斯-塞得爾迭代法高斯列主元消去法插值擬合偏微分方程式的差分數值解法分類、數值解法和差分格式分類數值解法差分格式橢圓型偏微分方程的數值解法拋物型偏微分方程的數值解法有壓管流管網的水力計算哈迪—克勞斯法管網的水力計算有限單元法簡單管道的水擊計算—特征線法短管的水力計算明渠流明渠非恒定流的水力計算—特征線法棱柱形明渠恒定非均勻漸變流水面曲線計算非棱柱形明渠恒定非均勻漸變流水面曲線計算天然河道水面曲線計算堰流和消能寬頂堰流水力計算消力池水力計算

數值計算基礎第1節(jié)非線性方程式的解法1． 簡單迭代法一、基本原理圖簡單迭代法原理圖如圖所示，由，選取初值，帶入該式得，一直進行下去，則有收斂判別式：（為高階小量，收斂判別常數）；收斂條件：。二、計算步驟（1）將變?yōu)?；?）選取初值；（3）迭代計算：；（4）比較，： 若，迭代結束，； 若，，返回（3）繼續(xù)迭代計算。例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）繼續(xù)迭代計算。例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．二等分法求方程的根。選取區(qū)間[a,b]使，則其根在區(qū)間[a,b]內。1．2．繼續(xù)二等分。每二等分一次區(qū)間縮短為原來的一半。第k次二等分后區(qū)間長度為：收斂判別式：計算步驟：（1）選取；使；（2）令，求；（3）若，則；若，則；（4）若，迭代結束，（或或）；若，返回（2）繼續(xù)進行二等分。例1-3．梯形渠道底寬，求。解：，選取。其中。程序：變量說明表：程序中：意義：Q，e流量，允許誤差Ac，BcHc臨界水深hc1，hc2區(qū)間（a,b）倆端點水深hc3區(qū)間（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節(jié)線性方程組的解法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'此處判斷該系數矩陣是否為嚴格對角占優(yōu)陣,若不想判斷也可將下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),"");'聲明:此程序是系數矩陣滿足嚴格對角占優(yōu)陣的情況下,Gauss-Seidel迭代法收斂;對于滿足不可約對角占優(yōu)陣的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),"");'聲明:此程序是系數矩陣滿足嚴格對角占優(yōu)陣的情況下,Jacobi迭代法收斂;對于滿足不可約對角占優(yōu)陣的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

第三節(jié)插值例已知消力坎的淹沒系數與相對淹沒度對應關系如下表：…….…….…….……..hs為下游水位超過坎頂的高度；H0為坎上全水頭。試編一分別用線性插值和拋物型插值的程序，計算相對淹沒度時的淹沒系數值。1．線性插值已知，求，要求滿足則令2．拉格朗日插值二次插值（拋物型插值），求。設滿足條件：3．n次插值（拉格朗日插值），求線性插值程序⑹變量說明表：程序中：意義：x(i)已知節(jié)點坐標y(i)已知節(jié)點函數值m待求函數值的節(jié)點坐標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節(jié) 擬合堰流：最小二乘法：偏差若要使取最小值，則①②例：直角三角形薄壁堰，求時，最小二乘法線性擬合直角三角形薄璧堰流量的程序⑺變量說明表：程序中：意義：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節(jié)分類、數值解法和差分格式1．分類二階線性偏微分方程式：典型例子：橢圓型偏微分方程式：拉普拉斯方程式拋物型偏微分方程式：擴散方程式雙曲型偏微分方程式：波動方程式2．數值解法主要步驟：①劃分網格；②構造差分方程式；③求解差分方程式。3．差分格式函數及在處的泰勒級數展開式分別為：（1）向前差分（一階精度）（2）向后差分（一階精度）（1）-（2）得：中心差分（二階精度）（1）+（2）得：中心差分（二階精度）

第2節(jié) 橢圓型偏微分方程式的數值解法基本方程式：（在域A內）邊界條件：（1）第一類邊界條件：（2） 第二類邊界條件：（3）第三類邊界條件：有限差分解法：（1）網格劃分：（2）構造差分方程式：，五點格式：（3）解差分方程式1．高斯迭代法2．高斯-賽德爾迭代法3．超松弛迭代法ω超松弛因子（ω=1時為高斯-賽德爾迭代法，ω〈1為亞松馳迭代，ω〉1為超松弛迭代）（計算試驗表明：ω取~之間數值時迭代次數最?。┡袆e收斂性：例：如圖所示為夾在兩個無限平行壁間的圓柱繞流，無限遠處的來流速度分布均勻，且，二平板間通過的流量為，圓柱半徑，試用有限差分法編寫求流函數在流場中分布的程序。解：（一）網格劃分：，，。（二）構造差分方程式：（在域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方向的網格步長。（三）解差分方程式：逐點迭代計算。每循環(huán)一次，進行收斂條件的判別。收斂判別條件：有限差分法計算兩個無限平行壁間圓柱繞流流函數在流場中分布的程序⑻變量說明表：程序中：公式中：意義：psi(i,j)ψ節(jié)點的流函數值，二維數組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節(jié)拋物型偏微分方程式的數值解法例：兩無限平板間油的流動，上平板水平運動，內u=0線性增加到，其后保持不變。求兩平板間的速度分布隨時間的變化。解：一．基本方程式邊界條件：初始條件：二．有限差分解法（一）網格劃分取空間步長其中；取時間步長其中。（二）構造差分方程式顯式格式：若令，則穩(wěn)定性條件：隱式格式：若令，則隱式格式無條件穩(wěn)定。初始條件：邊界條件：（三）解差分方程式顯格式差分法計算兩平板間速度分布隨時間變化的程序⑼變量說明表：程序中：公式中：意義：n，Hn，a/m時段數，距離步長MM距離a被等份的份數AA兩平板間的距離Vν液體的運動粘滯系數U0平板的最終速度RR顯格式的穩(wěn)定系數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節(jié) 管網的水力計算哈迪克勞斯法例1：如圖所示環(huán)狀管網，假設。試求：（1） 各管段的流量分配；（2） 各節(jié)點的測壓管水頭。管段12345678910管長600600300600300900300600300600管徑沿程阻力系數解：（一）主要公式連續(xù)性方程：。n為該點的管段數目，流進節(jié)點，流出節(jié)點。能量方程式：。環(huán)的正方向與該管段流向相同時，否則，，所以：（二）（1）假設各管段的流量分配，使其滿足連續(xù)性方程；（2）A環(huán)路水頭殘差；設的修正流量為，使；即：所以：，定性分析其正確性：如果初分流量時只是數值誤差，不是方向“誤差”的話，分母總為正值。則若分子（即殘差）為正，需減小。亦即正的水頭損失太大了，需減小流量；而負的水頭損失太小了，需增加流量。這是正確的。反之亦然。（3）計算步驟：1．根據各管段的計算；2．任意選擇各管段中的，滿足連續(xù)性方程（給定各管段的流動方向）；3．選定各閉合環(huán)路的正流動方向（如順時針方向）來確定各管段的，并構成矩陣，k-環(huán)路編號（行），I-管段編號（列）。對于環(huán)路中不包含的管段，矩陣元素。當該管段流向與該環(huán)方向相同時，；當該管段流向與該環(huán)方向相反時，。矩陣：（第k閉合環(huán)路--）ik123456789101-1+100000-10+1200+1+10000-1-130000+1-1-1+1+104.計算（對每個閉合環(huán)路），擴展為：。不在該環(huán)時，；在該環(huán)時。5．計算各管段中的流量，，擴展為：。對于非共用管，只有一個不為零；對于平面共用管，只有兩個不為零，并且當均規(guī)定順時針方向為環(huán)路正方向時，一個為“+1”，一個為“-1”。6．判別計算誤差（如的允許誤差），是否滿足：（其中為總環(huán)數），若不滿足，返回第4步繼續(xù)進行各環(huán)路修正流量的計算。7．計算各節(jié)點的測壓管水頭（需構造節(jié)點及管段號之間的關系數組）。幾點解釋：（1）只管數值，不管方向。流量方向事先已假定，最終只要不修正為負值，方向就與原假設方向相同。如若修正為負值，則與原假設方向相反。（2）連續(xù)性方程中流量方向用符號函數表示。流入為“+1”，流出為“-1”。（3）能量方程式中只表示修正流量的數值。為正值時，即為均按順時針方向增加流量；為負值時，即為均按逆時針方向增加流量。（4）的意義為：；。哈迪-克勞斯法管網的水力計算的程序⑾變量說明表：程序中：公式中：意義：d(i)，l(i)d,l各管段的管徑，管長y(i)λ各管段的沿程阻力系數R(i)R公式中的系數Q(i，j)Q第j次迭代第i個管段的流量M(i，j)表示管段中的流動方向是否與環(huán)路方向一致的符號矩陣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‘求各節(jié)點的水頭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節(jié)管網的水力計算有限單元法哈迪—克勞斯法是逐步漸進的經典解法，在大管徑小流量情況下不能很快收斂。一、單元方程式水流從令則（流量流出節(jié)點為“+”，流入節(jié)點為“—”）構成整體方程組（加入流量為“+”，消耗流量為“—”）以節(jié)點2為例：單元（1） 單元（3）單元（4）所以i單元，k,j兩個節(jié)點，ki放到（k,k）,(j,j)位置上，-ki放到（k,j）,(j,k)位置上。單元編號12345始端13223末端21344求解總體方程組（邊界條件與解法）為了求解方程組，至少必須已知一個節(jié)點的水頭值。由于流量只取決于水頭差而不是水頭本身，如果需要，可以任意規(guī)定某個節(jié)點的水頭值，對流量計算的結果也毫無影響。邊界條件引入的兩種方法：如本例四個節(jié)點的連續(xù)性方程，只有三個節(jié)點是獨立的。即其系數矩陣的行列式值為零，有無窮多個解。只有引入一個邊界條件（一節(jié)點的水頭值），才可求解。這也符合其物理意義。若：第一種處理方法：不減少方程的個數，直接引入。第二種處理方法：減少方程的個數，移項處理。高斯—賽德爾迭代法：當流動在紊流的阻力平方區(qū)時，單元特征系數和總體矩陣[k]的各元素都不是常數。它們都是水頭的函數。因此總體矩陣方程式是非線性的代數方程組，即解非線性代數方程組可采用高斯-塞德爾迭代法，牛頓-拉普森切線法等。牛頓-拉普森切線法比高斯-塞德爾迭代法收斂速度快。但高斯—賽德爾迭代法編制程序比較簡單。若用高斯—賽德爾迭代法計算步驟如下：（1）假設待求各節(jié)點的水頭值，即給定迭代初值；（2）用初值計算矩陣；（3）用計算待求的各節(jié)點的水頭值，頂標“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節(jié)點水頭的初值，“節(jié)點消耗”H1(i)Hj節(jié)點的水頭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節(jié)簡單管道的水擊計算特征線法（1）基本方程式（坡度很?。┓呛愣鬟B續(xù)性方程：f—沿程阻力系數（λ）運動方程式：所以順特征線方程：順特征方程：逆特征線方程：逆特征方程：（2）差分方程式(N個網格)順特征差分方程：逆特征差分方程：（3）邊界條件上游邊界條件：上游為水庫：下游邊界條件：若按孔口出流計算（閥門斷面）：為閥門相對開度。特征線法計算簡單管路中水擊壓強的程序⒀例：有一如圖所示的水電站的引水鋼管，在鋼管末端裝有自動調節(jié)閥門，使閥門按直線規(guī)律關閉，鋼管長，水庫最高水位與管道閥門斷面的高差，閥門全開時管中最大流速，水擊波速，設閥門完全關閉的時間，試編寫程序求：閥門斷面處的最大水擊壓強。解：（一）.主要公式（1）.基本方程式對于坡度較小的有壓管道，非恒定流的連續(xù)方程式和運動方程式為：式中是沿程阻力系數，而其中波速計算公式：變量說明表：程序中：公式中：意義：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節(jié)短管的水力計算短管的水力計算程序（10）例：如圖所示為新的鑄鐵管道，已知上下游水位差，管長，總局部阻力系數（不包括出口損失），管壁的當量粗糙度，水溫為，其運動粘滯系數，試編寫程序求：當流量時的管徑；當管徑時的流量。解：①管徑d的計算：，其中。牛頓切線法：計算過程：⑴假設λ；⑵求A,B,用牛頓切線法求d；⑶求ν,Re,計算λ；⑷比較λ相對誤差是否滿足精度要求。若滿足精度要求則計算結束，否則返回⑵繼續(xù)循環(huán)計算。②已知，求Q：計算過程：⑴假設λ；⑵計算μs,v;⑶計算Re，選取公式，求λ;⑷比較的相對誤差是否滿足精度要求：若滿足，求Q，結束計算；若不滿足，則返回⑵循環(huán)計算。③已知Q,d,求H：代入，即可求得H。關于沿程阻力系數計算的討論：上述的第二式：稱為布拉休斯公式。此式適用于紊流的光滑區(qū)。又有適用于紊流粗糙區(qū)的經驗公式為：稱為希弗林松公式。上述的第三式：可看成是紊流光滑區(qū)的布拉休斯公式與粗糙區(qū)的希弗林松公式的綜合。亦即可看成是適用于紊流三個區(qū)的通用公式。如若采用上一套計算公式進行計算，無論計算d或Q，均先假設，把的計算納入到d的牛頓切線法迭代計算或計算Q的簡單迭代法循環(huán)中即可。若采用柯列布魯克—懷特公式：求，則可用二等分法計算。構造關于的函數如下：設，則若，則判別下式是否成立：若成立，計算結束；否則，返回（*）步繼續(xù)進行二等分?？铝胁剪斂恕獞烟毓绞沁m用于紊流過渡粗糙區(qū)的公式，亦可以看成是適用于紊流三個區(qū)的通用公式。變量說明表：程序中；公式中：意義：hH水頭QQ流量LL管長DD管徑Sz∑ζ未出口的局部阻力系數KsKs管徑的當量粗糙度Nuν水的運動粘滯系數Fλ沿程阻力系數fEps1,eps2管徑和沿程阻力系數的允許的誤差j通道管號i循環(huán)變量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‘進行循環(huán)迭代Elsed=d1EndIfNextiElsei=1qq:us=1/Sqr(f*l/d+sz+1#)'是求沿程阻力系數的循環(huán)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節(jié)明渠非恒定流的計算特征線法基本方程組（以為變量，梯形斷面明渠）其中特征線方程和特征方程順特征線方程：順特征方程：即：逆特征線方程：逆特征方程：即：特征線差分方程和特征差分方程（庫朗格式）順特征線差分方程：順特征差分方程：即：逆特征線差分方程：逆特征差分方程：即：點L的值據點A和M的值線性插值計算；點R的值據點M和B的值線性插值計算。同理同理上兩式相減得：由順特征差分方程求：或由逆特征差分方程求：邊界條件上游：已知，則下游：給定關系曲線。穩(wěn)定性條件要求時間步長滿足下式：例：梯形斷面明渠渠中為均勻流，末端斷面水深。上游洪水過程如下表：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