計(jì)算水力學(xué)基礎(chǔ)

--本頁僅作為文檔封面，使用時(shí)請(qǐng)直接刪除即可--

--內(nèi)頁可以根據(jù)需求調(diào)整合適字體及大小本頁僅作為文檔封面，使用時(shí)請(qǐng)直接刪除即可--

--內(nèi)頁可以根據(jù)需求調(diào)整合適字體及大小--計(jì)算水力學(xué)基礎(chǔ)(總79頁)PAGE計(jì)算水力學(xué)基礎(chǔ)李占松編著鄭州大學(xué)水利與環(huán)境學(xué)院

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

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

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

數(shù)值計(jì)算基礎(chǔ)第1節(jié)非線性方程式的解法1． 簡(jiǎn)單迭代法一、基本原理圖簡(jiǎn)單迭代法原理圖如圖所示，由，選取初值，帶入該式得，一直進(jìn)行下去，則有收斂判別式：（為高階小量，收斂判別常數(shù)）；收斂條件：。二、計(jì)算步驟（1）將變?yōu)?；?）選取初值；（3）迭代計(jì)算：；（4）比較，： 若，迭代結(jié)束，； 若，，返回（3）繼續(xù)迭代計(jì)算。例1-1已知梯形斷面底寬，邊坡系數(shù)，，，m3/s，試編寫用迭代法求此渠道中正常水深的程序（計(jì)算允許誤差）。解：若已知正常水深求底寬，則底寬的計(jì)算表達(dá)式為，其迭代過程，與正常水深的迭代過程類似。程序：簡(jiǎn)單迭代法求解正常水深表變量說明變量名意義Q，b渠道流量，底寬m，n渠道邊坡系數(shù)，粗糙系數(shù)i渠道底坡h水深PrivateSubCommand1_Click()DimQ!,b!,h!,h0!,m!,i!,n!,eps!Q=Val(InputBox("請(qǐng)輸入渠道流量Q="))b=Val(InputBox("請(qǐng)輸入渠道底寬b="))n=Val(InputBox("請(qǐng)輸入渠道粗糙系數(shù)n="))m=Val(InputBox("請(qǐng)輸入渠道邊坡系數(shù)m="))i=Val(InputBox("請(qǐng)輸入渠道底坡i="))eps=Val(InputBox("請(qǐng)輸入精度e="))Doh=(n*Q/i^^*(b+2*h0*(1+m^2)^^/(b+m*h0)IfAbs(h-h0)<epsThenExitDoh0=hLoopPrint"該渠道的正常水深為：";Format(h,"");"米"EndSub依次輸入下列數(shù)據(jù)：Q=15，b=8，n=，m=，i=，eps=，輸出結(jié)果為：“該渠道的正常水深為：米”。2.牛頓-拉普森（Newton-Raphson）切線法計(jì)算溢流壩下游收縮斷面水深一、基本原理圖切線法原理圖已知方程，求解該方程的根。，，……，)())()(xfxfxx???收斂判別式：二、計(jì)算步驟（1）由函數(shù)求其導(dǎo)數(shù)；（2）選取初值；（3）迭代計(jì)算：；（4）比較：若，迭代結(jié)束，；若，，返回（3）繼續(xù)迭代計(jì)算。例1-2．已知某溢流壩上游斷面對(duì)下游河底的比能，矩形斷面河道上，為流速系數(shù)，試用牛頓-拉普森方法編寫計(jì)算的程序。為收縮斷面的水深。解：由整理得：（*）令，，（**）（***）將（**）代入（*）式，以及對(duì)（***）式求導(dǎo)可得，由牛頓-拉普森切線法，可得收縮斷面水深的迭代公式為收斂判別式：。對(duì)于本題，也可以利用簡(jiǎn)單迭代法求解收縮斷面的水深，其迭代公式為程序：牛頓-拉普森切線法求解收縮斷面水深表變量說明變量名意義Eo壩上游斷面的總比能q，x單寬流量，壩面流速系數(shù)A，常數(shù)B，常數(shù)C，常數(shù)linjieH引入的函數(shù)Hc下游收縮斷面水深PublicFunctionlinjieH(x!,q!,hc!)'建立函數(shù)Dimh1!,A!,B!,C!,fhc!,fh1c!g=Eo=Val(InputBox("請(qǐng)輸入斷面比能Eo"))X=Val(InputBox("請(qǐng)輸入流速系數(shù)x"))q=Val(InputBox("請(qǐng)輸入斷面流量q"))e=Val(InputBox("請(qǐng)輸入計(jì)算精度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!)'調(diào)用函數(shù)PrintFormat(hc,"")EndSub3．二等分法求方程的根。選取區(qū)間[a,b]使，則其根在區(qū)間[a,b]內(nèi)。1．2．繼續(xù)二等分。每二等分一次區(qū)間縮短為原來的一半。第k次二等分后區(qū)間長(zhǎng)度為：收斂判別式：計(jì)算步驟：（1）選取；使；（2）令，求；（3）若，則；若，則；（4）若，迭代結(jié)束，（或或）；若，返回（2）繼續(xù)進(jìn)行二等分。例1-3．梯形渠道底寬，求。解：，選取。其中。程序：變量說明表：程序中：意義：Q，e流量，允許誤差A(yù)c，BcHc臨界水深hc1，hc2區(qū)間（a,b）倆端點(diǎn)水深hc3區(qū)間（a,b）的中點(diǎn)值b,m渠底寬度，邊坡系數(shù)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("請(qǐng)輸入矩形斷面的底寬b"))m!=Val(InputBox("請(qǐng)輸入邊坡系數(shù)m"))Q!=Val(InputBox("請(qǐng)輸入流量Q"))e!=Val(InputBox("請(qǐng)輸入計(jì)算精度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．高斯-賽德爾迭代法三元一次方程組：初值：雅可比迭代法：賽德爾迭代法：第步：雅可比迭代法：賽德爾迭代法：收斂判別式：對(duì)于計(jì)算步驟：1． 給定初值：2． 迭代計(jì)算：（1）雅可比迭代法：（2）賽德爾迭代法：3．收斂判別式：賽德爾迭代法程序（4）變量說明表：程序中：意義：a(i,j)方程組系數(shù)矩陣元素b(i)方程組右端項(xiàng)元素x0(i)，x1(i)方程組待求的未知量N矩陣的維數(shù)DimiAsInteger,jAsInteger,nAsInteger,sumAsInteger,s1AsDouble,s2AsDoubleDima()AsInteger,b()AsInteger,x0()AsDouble,x1()AsDouble,EPSAsDouble,t#,z#PublicFunctionGaussSeidel()n=InputBox("輸入矩陣的維數(shù)：")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'此處判斷該系數(shù)矩陣是否為嚴(yán)格對(duì)角占優(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),"");'聲明:此程序是系數(shù)矩陣滿足嚴(yán)格對(duì)角占優(yōu)陣的情況下,Gauss-Seidel迭代法收斂;對(duì)于滿足不可約對(duì)角占優(yōu)陣的Gauss-Seidel迭代法收斂情況下,'有興趣的讀者可查看有關(guān)資料,自己動(dòng)手試試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("輸入矩陣的維數(shù)：")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),"");'聲明:此程序是系數(shù)矩陣滿足嚴(yán)格對(duì)角占優(yōu)陣的情況下,Jacobi迭代法收斂;對(duì)于滿足不可約對(duì)角占優(yōu)陣的Jacobi迭代法收斂情況下,'有興趣的讀者可查看有關(guān)資料,自己動(dòng)手試試'10x1-x2-2x3=72-x1+10x2-2x3=83-x1-x2-5x3=42解為：x1=x2=x3=EndFunctionPrivateSubForm_Click()CallJacobiEndSub2．高斯列主元消去法計(jì)算步驟：1．選主元（第k步）（k=1,2,…,n）假設(shè)主元在第L行換元：把換為主元（j=k,k+1,…,n,n+1）2．消元（第k步）（k=1,2,…,n）（如第2次消元，即k=2，第3行，即i=3，第4列，即j=4,數(shù)據(jù)的變化應(yīng)為：）進(jìn)行到第n步得：3．回代例：高斯列主元消去法求解線性代數(shù)方程組的程序⑸變量說明表：程序中：意義：a(i,j)線性方程組的系數(shù)矩陣元素n線性方程組的維數(shù)b(i)方程組右端項(xiàng)元素maxV對(duì)角元素的最大值maxI，maxJ最大元素在I行J列OptionExplicitPrivateSubCommand1_Click()Dima()AsDouble,b()AsDouble,iAsInteger,nAsInteger,j!n=3ReDima(n,n)AsDouble,b(n)AsDoublen=InputBox("輸入矩陣的維數(shù)：")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))‘將主對(duì)角線上的最大元素找出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é)插值例已知消力坎的淹沒系數(shù)與相對(duì)淹沒度對(duì)應(yīng)關(guān)系如下表：…….…….…….……..hs為下游水位超過坎頂?shù)母叨?；H0為坎上全水頭。試編一分別用線性插值和拋物型插值的程序，計(jì)算相對(duì)淹沒度時(shí)的淹沒系數(shù)值。1．線性插值已知，求，要求滿足則令2．拉格朗日插值二次插值（拋物型插值），求。設(shè)滿足條件：3．n次插值（拉格朗日插值），求線性插值程序⑹變量說明表：程序中：意義：x(i)已知節(jié)點(diǎn)坐標(biāo)y(i)已知節(jié)點(diǎn)函數(shù)值m待求函數(shù)值的節(jié)點(diǎn)坐標(biāo)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("請(qǐng)輸入一個(gè)數(shù)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é) 擬合堰流：最小二乘法：偏差若要使取最小值，則①②例：直角三角形薄壁堰，求時(shí)，最小二乘法線性擬合直角三角形薄璧堰流量的程序⑺變量說明表：程序中：意義：H直角三角形薄壁堰的堰上水頭t1(i)，t2(i)堰上水頭，流量，一維數(shù)組sum1，sum2公式中AA1，AA2，V中間變量Q欲求流量PrivateSubCommand1_Click()DimH#H=Val(InputBox("請(qǐng)輸入直角三角形薄壁堰的堰上水頭"))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

第二章偏微分方程式的差分?jǐn)?shù)值解法第1節(jié)分類、數(shù)值解法和差分格式1．分類二階線性偏微分方程式：典型例子：橢圓型偏微分方程式：拉普拉斯方程式拋物型偏微分方程式：擴(kuò)散方程式雙曲型偏微分方程式：波動(dòng)方程式2．?dāng)?shù)值解法主要步驟：①劃分網(wǎng)格；②構(gòu)造差分方程式；③求解差分方程式。3．差分格式函數(shù)及在處的泰勒級(jí)數(shù)展開式分別為：（1）向前差分（一階精度）（2）向后差分（一階精度）（1）-（2）得：中心差分（二階精度）（1）+（2）得：中心差分（二階精度）

第2節(jié) 橢圓型偏微分方程式的數(shù)值解法基本方程式：（在域A內(nèi)）邊界條件：（1）第一類邊界條件：（2） 第二類邊界條件：（3）第三類邊界條件：有限差分解法：（1）網(wǎng)格劃分：（2）構(gòu)造差分方程式：，五點(diǎn)格式：（3）解差分方程式1．高斯迭代法2．高斯-賽德爾迭代法3．超松弛迭代法ω超松弛因子（ω=1時(shí)為高斯-賽德爾迭代法，ω〈1為亞松馳迭代，ω〉1為超松弛迭代）（計(jì)算試驗(yàn)表明：ω取~之間數(shù)值時(shí)迭代次數(shù)最?。┡袆e收斂性：例：如圖所示為夾在兩個(gè)無限平行壁間的圓柱繞流，無限遠(yuǎn)處的來流速度分布均勻，且，二平板間通過的流量為，圓柱半徑，試用有限差分法編寫求流函數(shù)在流場(chǎng)中分布的程序。解：（一）網(wǎng)格劃分：，，。（二）構(gòu)造差分方程式：（在域A內(nèi)）第一類邊界條件：（1）；（2）；；；（3）；第二類邊界條件：（4）。五點(diǎn)差分格式：超松弛迭代法：（ω=~）邊界條件：（1）（ab邊界）；（2）；（3）bc邊界的處理：(a)(b)；（4）；（ae邊界）；（5）cd邊界的處理:的計(jì)算：(a)直接轉(zhuǎn)移：(b)直線內(nèi)插：如上圖所示，。則，其中a表示P1與網(wǎng)格線和圓柱壁面交點(diǎn)之間的距離，是（7，3）和（7，4）兩網(wǎng)格點(diǎn)之間的距離即y方向的網(wǎng)格步長(zhǎng)；b表示P2與j=2網(wǎng)格線和圓柱壁面交點(diǎn)之間的距離，是（6，2）和（5，2）兩網(wǎng)格點(diǎn)之間的距離即x方向的網(wǎng)格步長(zhǎng)。（三）解差分方程式：逐點(diǎn)迭代計(jì)算。每循環(huán)一次，進(jìn)行收斂條件的判別。收斂判別條件：有限差分法計(jì)算兩個(gè)無限平行壁間圓柱繞流流函數(shù)在流場(chǎng)中分布的程序⑻變量說明表：程序中：公式中：意義：psi(i,j)ψ節(jié)點(diǎn)的流函數(shù)值，二維數(shù)組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é)拋物型偏微分方程式的數(shù)值解法例：兩無限平板間油的流動(dòng)，上平板水平運(yùn)動(dòng)，內(nèi)u=0線性增加到，其后保持不變。求兩平板間的速度分布隨時(shí)間的變化。解：一．基本方程式邊界條件：初始條件：二．有限差分解法（一）網(wǎng)格劃分取空間步長(zhǎng)其中；取時(shí)間步長(zhǎng)其中。（二）構(gòu)造差分方程式顯式格式：若令，則穩(wěn)定性條件：隱式格式：若令，則隱式格式無條件穩(wěn)定。初始條件：邊界條件：（三）解差分方程式顯格式差分法計(jì)算兩平板間速度分布隨時(shí)間變化的程序⑼變量說明表：程序中：公式中：意義：n，Hn，a/m時(shí)段數(shù)，距離步長(zhǎng)MM距離a被等份的份數(shù)AA兩平板間的距離Vν液體的運(yùn)動(dòng)粘滯系數(shù)U0平板的最終速度RR顯格式的穩(wěn)定系數(shù)u(i,j)點(diǎn)流速Tj中間變量DT△t時(shí)間步長(zhǎng)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é) 管網(wǎng)的水力計(jì)算哈迪克勞斯法例1：如圖所示環(huán)狀管網(wǎng)，假設(shè)。試求：（1） 各管段的流量分配；（2） 各節(jié)點(diǎn)的測(cè)壓管水頭。管段12345678910管長(zhǎng)600600300600300900300600300600管徑沿程阻力系數(shù)解：（一）主要公式連續(xù)性方程：。n為該點(diǎn)的管段數(shù)目，流進(jìn)節(jié)點(diǎn)，流出節(jié)點(diǎn)。能量方程式：。環(huán)的正方向與該管段流向相同時(shí)，否則，，所以：（二）（1）假設(shè)各管段的流量分配，使其滿足連續(xù)性方程；（2）A環(huán)路水頭殘差；設(shè)的修正流量為，使；即：所以：，定性分析其正確性：如果初分流量時(shí)只是數(shù)值誤差，不是方向“誤差”的話，分母總為正值。則若分子（即殘差）為正，需減小。亦即正的水頭損失太大了，需減小流量；而負(fù)的水頭損失太小了，需增加流量。這是正確的。反之亦然。（3）計(jì)算步驟：1．根據(jù)各管段的計(jì)算；2．任意選擇各管段中的，滿足連續(xù)性方程（給定各管段的流動(dòng)方向）；3．選定各閉合環(huán)路的正流動(dòng)方向（如順時(shí)針方向）來確定各管段的，并構(gòu)成矩陣，k-環(huán)路編號(hào)（行），I-管段編號(hào)（列）。對(duì)于環(huán)路中不包含的管段，矩陣元素。當(dāng)該管段流向與該環(huán)方向相同時(shí)，；當(dāng)該管段流向與該環(huán)方向相反時(shí)，。矩陣：（第k閉合環(huán)路--）ik123456789101-1+100000-10+1200+1+10000-1-130000+1-1-1+1+104.計(jì)算（對(duì)每個(gè)閉合環(huán)路），擴(kuò)展為：。不在該環(huán)時(shí)，；在該環(huán)時(shí)。5．計(jì)算各管段中的流量，，擴(kuò)展為：。對(duì)于非共用管，只有一個(gè)不為零；對(duì)于平面共用管，只有兩個(gè)不為零，并且當(dāng)均規(guī)定順時(shí)針方向?yàn)榄h(huán)路正方向時(shí)，一個(gè)為“+1”，一個(gè)為“-1”。6．判別計(jì)算誤差（如的允許誤差），是否滿足：（其中為總環(huán)數(shù)），若不滿足，返回第4步繼續(xù)進(jìn)行各環(huán)路修正流量的計(jì)算。7．計(jì)算各節(jié)點(diǎn)的測(cè)壓管水頭（需構(gòu)造節(jié)點(diǎn)及管段號(hào)之間的關(guān)系數(shù)組）。幾點(diǎn)解釋：（1）只管數(shù)值，不管方向。流量方向事先已假定，最終只要不修正為負(fù)值，方向就與原假設(shè)方向相同。如若修正為負(fù)值，則與原假設(shè)方向相反。（2）連續(xù)性方程中流量方向用符號(hào)函數(shù)表示。流入為“+1”，流出為“-1”。（3）能量方程式中只表示修正流量的數(shù)值。為正值時(shí)，即為均按順時(shí)針方向增加流量；為負(fù)值時(shí)，即為均按逆時(shí)針方向增加流量。（4）的意義為：；。哈迪-克勞斯法管網(wǎng)的水力計(jì)算的程序⑾變量說明表：程序中：公式中：意義：d(i)，l(i)d,l各管段的管徑，管長(zhǎng)y(i)λ各管段的沿程阻力系數(shù)R(i)R公式中的系數(shù)Q(i，j)Q第j次迭代第i個(gè)管段的流量M(i，j)表示管段中的流動(dòng)方向是否與環(huán)路方向一致的符號(hào)矩陣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)‘輸入已知數(shù)據(jù)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‘求系數(shù)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é)點(diǎn)的水頭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é)管網(wǎng)的水力計(jì)算有限單元法哈迪—克勞斯法是逐步漸進(jìn)的經(jīng)典解法，在大管徑小流量情況下不能很快收斂。一、單元方程式水流從令則（流量流出節(jié)點(diǎn)為“+”，流入節(jié)點(diǎn)為“—”）構(gòu)成整體方程組（加入流量為“+”，消耗流量為“—”）以節(jié)點(diǎn)2為例：?jiǎn)卧?） 單元（3）單元（4）所以i單元，k,j兩個(gè)節(jié)點(diǎn)，ki放到（k,k）,(j,j)位置上，-ki放到（k,j）,(j,k)位置上。單元編號(hào)12345始端13223末端21344求解總體方程組（邊界條件與解法）為了求解方程組，至少必須已知一個(gè)節(jié)點(diǎn)的水頭值。由于流量只取決于水頭差而不是水頭本身，如果需要，可以任意規(guī)定某個(gè)節(jié)點(diǎn)的水頭值，對(duì)流量計(jì)算的結(jié)果也毫無影響。邊界條件引入的兩種方法：如本例四個(gè)節(jié)點(diǎn)的連續(xù)性方程，只有三個(gè)節(jié)點(diǎn)是獨(dú)立的。即其系數(shù)矩陣的行列式值為零，有無窮多個(gè)解。只有引入一個(gè)邊界條件（一節(jié)點(diǎn)的水頭值），才可求解。這也符合其物理意義。若：第一種處理方法：不減少方程的個(gè)數(shù)，直接引入。第二種處理方法：減少方程的個(gè)數(shù)，移項(xiàng)處理。高斯—賽德爾迭代法：當(dāng)流動(dòng)在紊流的阻力平方區(qū)時(shí)，單元特征系數(shù)和總體矩陣[k]的各元素都不是常數(shù)。它們都是水頭的函數(shù)。因此總體矩陣方程式是非線性的代數(shù)方程組，即解非線性代數(shù)方程組可采用高斯-塞德爾迭代法，牛頓-拉普森切線法等。牛頓-拉普森切線法比高斯-塞德爾迭代法收斂速度快。但高斯—賽德爾迭代法編制程序比較簡(jiǎn)單。若用高斯—賽德爾迭代法計(jì)算步驟如下：（1）假設(shè)待求各節(jié)點(diǎn)的水頭值，即給定迭代初值；（2）用初值計(jì)算矩陣；（3）用計(jì)算待求的各節(jié)點(diǎn)的水頭值，頂標(biāo)“1”表示迭代次數(shù)；（4）將求得的水頭值與初值相比較，如果滿足下式：計(jì)算結(jié)束。否則將符給，重復(fù)（2）~（4）步驟，直到滿足精度要求為止。（5）求各單元的流量。例1。單元編號(hào)12345678910始端1123456742末端7234567888有限單元法管網(wǎng)水力計(jì)算程序（12）變量說明表：程序中：公式中：意義：d(i)，l(i)d,l各管段的管徑，管長(zhǎng)y(i)n各管段的沿程粗糙系數(shù)Q(i)Q第j次迭代第i個(gè)管段的流量H0(i)，c(i)H0j,cj節(jié)點(diǎn)水頭的初值，“節(jié)點(diǎn)消耗”H1(i)Hj節(jié)點(diǎn)的水頭n1(i)，n2(i)單元端點(diǎn)編號(hào)A(ii,jj)總體矩陣元素，二維數(shù)組K(i)單元特征系數(shù)，一維數(shù)組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)‘輸入已知數(shù)據(jù)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é)簡(jiǎn)單管道的水擊計(jì)算特征線法（1）基本方程式（坡度很?。┓呛愣鬟B續(xù)性方程：f—沿程阻力系數(shù)（λ）運(yùn)動(dòng)方程式：所以順特征線方程：順特征方程：逆特征線方程：逆特征方程：（2）差分方程式(N個(gè)網(wǎng)格)順特征差分方程：逆特征差分方程：（3）邊界條件上游邊界條件：上游為水庫：下游邊界條件：若按孔口出流計(jì)算（閥門斷面）：為閥門相對(duì)開度。特征線法計(jì)算簡(jiǎn)單管路中水擊壓強(qiáng)的程序⒀例：有一如圖所示的水電站的引水鋼管，在鋼管末端裝有自動(dòng)調(diào)節(jié)閥門，使閥門按直線規(guī)律關(guān)閉，鋼管長(zhǎng)，水庫最高水位與管道閥門斷面的高差，閥門全開時(shí)管中最大流速，水擊波速，設(shè)閥門完全關(guān)閉的時(shí)間，試編寫程序求：閥門斷面處的最大水擊壓強(qiáng)。解：（一）.主要公式（1）.基本方程式對(duì)于坡度較小的有壓管道，非恒定流的連續(xù)方程式和運(yùn)動(dòng)方程式為：式中是沿程阻力系數(shù)，而其中波速計(jì)算公式：變量說明表：程序中：公式中：意義：DD管徑LL鋼管的長(zhǎng)度H0h0水頭計(jì)算初值V0ν0流速計(jì)算初值CC=c/g水擊速度Ts閘門完全關(guān)閉時(shí)間F沿程阻力系數(shù)NN=l/dx管段數(shù)Nt時(shí)段數(shù)H(I,j),v(I,j)管中水頭及流速，二維數(shù)組T(i),tn(i)t,時(shí)間及閥門開度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("沿程阻力系數(shù)f"))v0=Val(InputBox("流速計(jì)算初值v0"))h0=Val(InputBox("水頭計(jì)算初值h0"))l=Val(InputBox("鋼管的長(zhǎng)度l"))ts=Val(InputBox("閘門完全關(guān)閉時(shí)間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‘特征法進(jìn)行運(yùn)算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é)短管的水力計(jì)算短管的水力計(jì)算程序（10）例：如圖所示為新的鑄鐵管道，已知上下游水位差，管長(zhǎng)，總局部阻力系數(shù)（不包括出口損失），管壁的當(dāng)量粗糙度，水溫為，其運(yùn)動(dòng)粘滯系數(shù)，試編寫程序求：當(dāng)流量時(shí)的管徑；當(dāng)管徑時(shí)的流量。解：①管徑d的計(jì)算：，其中。牛頓切線法：計(jì)算過程：⑴假設(shè)λ；⑵求A,B,用牛頓切線法求d；⑶求ν,Re,計(jì)算λ；⑷比較λ相對(duì)誤差是否滿足精度要求。若滿足精度要求則計(jì)算結(jié)束，否則返回⑵繼續(xù)循環(huán)計(jì)算。②已知，求Q：計(jì)算過程：⑴假設(shè)λ；⑵計(jì)算μs,v;⑶計(jì)算Re，選取公式，求λ;⑷比較的相對(duì)誤差是否滿足精度要求：若滿足，求Q，結(jié)束計(jì)算；若不滿足，則返回⑵循環(huán)計(jì)算。③已知Q,d,求H：代入，即可求得H。關(guān)于沿程阻力系數(shù)計(jì)算的討論：上述的第二式：稱為布拉休斯公式。此式適用于紊流的光滑區(qū)。又有適用于紊流粗糙區(qū)的經(jīng)驗(yàn)公式為：稱為希弗林松公式。上述的第三式：可看成是紊流光滑區(qū)的布拉休斯公式與粗糙區(qū)的希弗林松公式的綜合。亦即可看成是適用于紊流三個(gè)區(qū)的通用公式。如若采用上一套計(jì)算公式進(jìn)行計(jì)算，無論計(jì)算d或Q，均先假設(shè)，把的計(jì)算納入到d的牛頓切線法迭代計(jì)算或計(jì)算Q的簡(jiǎn)單迭代法循環(huán)中即可。若采用柯列布魯克—懷特公式：求，則可用二等分法計(jì)算。構(gòu)造關(guān)于的函數(shù)如下：設(shè)，則若，則判別下式是否成立：若成立，計(jì)算結(jié)束；否則，返回（*）步繼續(xù)進(jìn)行二等分?？铝胁剪斂恕獞烟毓绞沁m用于紊流過渡粗糙區(qū)的公式，亦可以看成是適用于紊流三個(gè)區(qū)的通用公式。變量說明表：程序中；公式中：意義：hH水頭QQ流量LL管長(zhǎng)DD管徑Sz∑ζ未出口的局部阻力系數(shù)KsKs管徑的當(dāng)量粗糙度Nuν水的運(yùn)動(dòng)粘滯系數(shù)Fλ沿程阻力系數(shù)fEps1,eps2管徑和沿程阻力系數(shù)的允許的誤差j通道管號(hào)i循環(huán)變量reRe雷諾數(shù)PrivateSubCommand1_Click()Dimv#h=Val(InputBox("水位差h"))q=Val(InputBox("流量q"))l=Val(InputBox("管長(zhǎng)l"))d=Val(InputBox("管徑d"))sz=Val(InputBox("末計(jì)出口的局部阻力系數(shù)sz"))ks=Val(InputBox("管徑的當(dāng)量粗糙度ks"))nu=Val(InputBox("水的運(yùn)動(dòng)粘滯系數(shù)nu"))f=Val(InputBox("沿程阻力系數(shù)f"))eps1=Val(InputBox("管徑和沿程阻力系數(shù)的允許的誤差eps1"))eps2=Val(InputBox("管徑和沿程阻力系數(shù)的允許的誤差eps2"))j=Val(InputBox("通道管號(hào)j"))pi=g=Ifj=1ThenFori=1To100v=4#*q*q/(pi*d*d)re=v*d/nu‘雷諾數(shù)的計(jì)算Ifre<2300Then‘沿程阻力系數(shù)的判斷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‘進(jìn)行循環(huán)迭代Elsed=d1EndIfNextiElsei=1qq:us=1/Sqr(f*l/d+sz+1#)'是求沿程阻力系數(shù)的循環(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é)明渠非恒定流的計(jì)算特征線法基本方程組（以為變量，梯形斷面明渠）其中特征線方程和特征方程順特征線方程：順特征方程：即：逆特征線方程：逆特征方程：即：特征線差分方程和特征差分方程（庫朗格式）順特征線差分方程：順特征差分方程：即：逆特征線差分方程：逆特征差分方程：即：點(diǎn)L的值據(jù)點(diǎn)A和M的值線性插值計(jì)算；點(diǎn)R的值據(jù)點(diǎn)M和B的值線性插值計(jì)算。同理同理上兩式相減得：由順特征差分方程求：或由逆特征差分方程求：邊界條件上游：已知，則下游：給定關(guān)系曲線。穩(wěn)定性條件要求時(shí)間步長(zhǎng)滿足下式：例：梯形斷面明渠渠中為均勻流，末端斷面水深。上游洪水過程如下表：t(s)060120180240300360420Q1517192123252524t(s)480540600660720780840900Q2322212019181715特征線法計(jì)算明渠非恒定流的程序⒄如圖所示，某梯形斷面渠道，下游與一水塘相接，已知：渠道長(zhǎng)度，底寬，邊坡系數(shù)，粗糙系數(shù)，底坡，通過流量時(shí)渠中為均勻流，其水深，這時(shí)渠道末端斷面水位與水塘中水位平齊，其后末端斷面水深隨流量按下式變化：今在渠道的上游流域降一場(chǎng)大雨，在上游端斷面測(cè)得洪水過程如下表所給：06012018024030036042015171921232525244805406006607207808409002322212019181716試用特征線法編寫程序，計(jì)算各時(shí)刻各斷面的水深、流量、流速及波速。變量說明表：程序中公式意義MM邊坡系數(shù)NN粗糙系數(shù)BB寬度S0I底坡H0H渠中水深Ds距離步長(zhǎng)Dt時(shí)間步長(zhǎng)A1A方程式中的系數(shù)a1C1C方程式中的系數(shù)c1H(l,k),Q(l,k)H,q水深，流量數(shù)組V(l,k),C(I,k)v,c=(g*A/B)^(1/2)斷面平均流速，波速數(shù)組PublicSubCommand1_Click()Dimh(12,25),q(12,25),v(12,25),c(12,25),qu(16)m=Val(InputBox("邊坡系數(shù)m"))n=Val(InputBox("粗糙系數(shù)n"))b=Val(InputBox("b"))s0=Val(InputBox("底坡s0"))h0=Val(InputBox("渠中水深h0"))ds=Val(InputBox("距離步長(zhǎng)ds"))dt=Val(InputBox("時(shí)間步長(zhǎng)dt"))a1=Val(InputBox("方程式中的系數(shù)a1"))b1=Val(InputBox("方程式中的系數(shù)b1"))c1=Val(InputBox("方程式中的系數(shù)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‘對(duì)流量賦值q(i,j)=qu(j)NextjForj=n1+1Tokq(1,j)=15#Nextjg=Forj=2Tok–1‘對(duì)內(nèi)點(diǎn)的p計(jì)算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)‘逆水流方向的絕對(duì)速度wf=bs*(vm+cm)‘順?biāo)鞣较虻慕^對(duì)速度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‘已知時(shí)層內(nèi)l插點(diǎn)的水深ql=qm+dt/ds*(qa-qm)*wf/bs‘已知時(shí)層內(nèi)l插點(diǎn)的流量hr=hm+dt/ds*(hm-hb)*wb/bs‘已知時(shí)層內(nèi)r插點(diǎn)的水深qr=qm+dt/ds*(qm-qb)*wb/bs‘已知時(shí)層內(nèi)r插點(diǎn)的流量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‘末端斷面的水深隨流量的變化的函數(shù)EndIfN