計算水力學基礎

1、計算水力學基礎李占松 編著鄭州大學水利與環(huán)境學院內 容 簡 介本講義是編者根據(jù)多年的教學實踐，并參考微機計算水力學（楊景芳編著，大連理工大學出版社出版，1991年5月第1版）等類似教材，取其精華，編寫而成的。目的是使讀者掌握通過計算機解水力學問題的方法，為解決更復雜的實際工程問題打下牢固的計算基礎。書中內容包括：數(shù)值計算基礎，偏微分方程式的差分解法，有限單元法；用這些方法解有壓管流、明渠流、閘孔出流、堰流、消能、地下水的滲流及平面勢流等計算問題。講義中的用FORTRAN77算法語言編寫的計算程序，幾乎包括了全部水力學的主要計算問題。另外，結合講授對象的實際情況，也提供了用VB算法語言編寫的計算

2、程序。VB程序編程人員的話為了更好地促進水利水電工程建筑專業(yè)的同學學好微機計算水力學這門學科，編程員借暑假休息的時間，利用我們專業(yè)目前所學的VB中的算法語言部分對水力學常見的計算題型編制成常用程序。希望大家能借此資料更好地學習微機計算水力學這門課程。本程序著重程序的可讀性，不苛求程序的過分技巧。對水力學中常用的計算題型，用我們現(xiàn)在所學的VB語言編制而成。由于編程員能力有限，程序中缺點和錯誤在所難免，望老師和同學及時給予批評指正。VB程序編程人員：黃渝桂 曹命凱前 言-計算水力學的形成與發(fā)展計算水力學作為一門新學科，形成于20世紀60年代中期。水力學問題中有比較復雜的紊流、分離、氣穴、水擊等流動

3、現(xiàn)象，并存在各種界面形式，如自由水面、分層流、交界面等。由各種流動現(xiàn)象而建立的數(shù)學模型（由微分方程表示的定解問題），例如連續(xù)方程、動量方程等組成的控制微分方程組，多具有非線性和非恒定性，只有少數(shù)特定條件下的問題，可根據(jù)求解問題的特性對方程和邊界條件作相應簡化，而得到其解析解。因此長期以來，水力學的發(fā)展只得主要藉助于物理模型試驗。隨著電子計算機和現(xiàn)代計算技術的發(fā)展，數(shù)值計算已逐漸成為一個重要的研究手段，發(fā)展至今，已廣泛應用與水利、航運、海洋、流體機械與流體工程等各種技術科學領域。計算水力學的特點是適應性強、應用面廣。首先流動問題的控制方程一般是非線性的，自變量多，計算域的幾何形狀任意，邊界條件復

4、雜，對這些無法求得解析解的問題，用數(shù)值解則能很好的滿足工程需要；其次可利用計算機進行各種數(shù)值試驗，例如，可選擇不同的流動參數(shù)進行試驗，可進行物理方程中各項的有效性和敏感性試驗，以便進行各種近似處理等。它不受物理模型試驗模型律的限制，比較省時省錢，有較多的靈活性。但數(shù)值計算一是依賴于基本方程的可靠性，且最終結果不能提供任何形式的解析表達式，只是有限個離散點上的數(shù)值解，并有一定的計算誤差；二是它不像物理模型試驗一開始就能給出流動現(xiàn)象并定性地描述，卻往往需要由原體觀測或物理實驗提供某些流動參數(shù)，并對建立的數(shù)學模型驗證；三是程序的編制及資料的收集、整理與正確利用，在很大程度上依賴于經驗與技巧。所以計算

5、水力學有自己的原理方法和特點，數(shù)值計算與理論分析觀測和試驗相互聯(lián)系、促進又不能相互代替，已成為目前解決復雜水流問題的主要手段之一，尤其是在研究流動過程物理機制時，更需要三者有機結合而互相取長補短。近三、四十年來，計算水力學有很大的發(fā)展，替代了經典水力學中的一些近似計算法和圖解法。例如水面曲線計算；管網(wǎng)和渠系的過水或輸沙（排污）能力的計算；有水輪機負荷改變時水力震蕩系統(tǒng)的穩(wěn)定性計算研究；流體機械過流部件的流道計算以及優(yōu)化設計，還有洪水波、河口潮流計算，以及各種流動條件下，不同排放形式的污染物混合計算等。上世紀70年代中期已從針對個別工程問題建立的單一數(shù)學模型，開始建立對整個流域洪泛區(qū)已建或規(guī)劃中

6、的水利水電工程進行系統(tǒng)模擬的系統(tǒng)模型。理論課題的研究中，對擴散問題、傳熱問題、邊界層問題、漩渦運動、紊流等問題的研究也有了很大的發(fā)展，并已開始計算非恒定的三維紊流問題。由于離散的基本原理不同，計算水力學可分為兩個分支：一是有限差分法，在此基礎上發(fā)展的有有限分析法；二是有限單元法，在此基礎上提出了邊界元法和混合元法，另外還有迎風有限元法等。目 錄第一章 數(shù)值計算基礎第一節(jié) 非線性方程式的解法1 迭代法2 牛頓-拉普森切線法3 二等分法第二節(jié) 線性方程組的解法1 高斯-塞得爾迭代法2 高斯列主元消去法第三節(jié) 插值第四節(jié) 擬合第二章 偏微分方程式的差分數(shù)值解法第一節(jié) 分類、數(shù)值解法和差分格式1 分類

7、2 數(shù)值解法3 差分格式第二節(jié) 橢圓型偏微分方程的數(shù)值解法第三節(jié) 拋物型偏微分方程的數(shù)值解法第三章 有壓管流第一節(jié) 管網(wǎng)的水力計算-哈迪克勞斯法第二節(jié) 管網(wǎng)的水力計算-有限單元法第三節(jié) 簡單管道的水擊計算特征線法第四節(jié) 短管的水力計算第四章 明渠流第一節(jié) 明渠非恒定流的水力計算特征線法第二節(jié) 棱柱形明渠恒定非均勻漸變流水面曲線計算第三節(jié) 非棱柱形明渠恒定非均勻漸變流水面曲線計算第四節(jié) 天然河道水面曲線計算第五章 堰流和消能第一節(jié) 寬頂堰流水力計算第二節(jié) 消力池水力計算第一章 數(shù)值計算基礎第1節(jié) 非線性方程式的解法1簡單迭代法一、基本原理圖1.1 簡單迭代法原理圖如圖1.2所示，由，選取初值，帶

8、入該式得，一直進行下去，則有收斂判別式：（為高階小量，收斂判別常數(shù)）；收斂條件：。二、計算步驟（1）將變?yōu)?；?）選取初值；（3）迭代計算：；（4）比較，：若，迭代結束，；若，返回（3）繼續(xù)迭代計算。例1-1 已知梯形斷面底寬，邊坡系數(shù)，m3/s，試編寫用迭代法求此渠道中正常水深的程序（計算允許誤差）。解：若已知正常水深求底寬，則底寬的計算表達式為，其迭代過程，與正常水深的迭代過程類似。程序1.1：簡單迭代法求解正常水深表1.1 變量說明變量名意義Q，b渠道流量，底寬m，n渠道邊坡系數(shù)，粗糙系數(shù)i渠道底坡h水深Private Sub Command1_Click()Dim Q!, b!, h!

9、, h0!, m!, i!, n!, eps!Q = Val(InputBox(請輸入渠道流量Q=)b = Val(InputBox(請輸入渠道底寬b=)n = Val(InputBox(請輸入渠道粗糙系數(shù)n=)m = Val(InputBox(請輸入渠道邊坡系數(shù)m=)i = Val(InputBox(請輸入渠道底坡i=)eps = Val(InputBox(請輸入精度e=)Doh = (n * Q / i 0.5) 0.6 * (b + 2 * h0 * (1 + m 2) 0.5) 0.4 / (b + m * h0)If Abs(h - h0) eps Then Exit Doh0 =

10、hLoopPrint 該渠道的正常水深為：; Format(h, 0.000); 米End Sub依次輸入下列數(shù)據(jù)：Q =15，b = 8，n = 0.025，m =1.5，i = 0.0009，eps=0.005，輸出結果為：“該渠道的正常水深為：1.265米”。2.牛頓-拉普森（Newton-Raphson）切線法計算溢流壩下游收縮斷面水深一、基本原理圖1.2 切線法原理圖已知方程，求解該方程的根。 ，)()(1nnnnxfxfxx-=+ 收斂判別式：二、計算步驟（1）由函數(shù)求其導數(shù)；（2）選取初值；（3）迭代計算：；（4）比較：若，迭代結束，；若，返回（3）繼續(xù)迭代計算。例1-2已知某溢

11、流壩上游斷面對下游河底的比能，矩形斷面河道上，為流速系數(shù)，試用牛頓-拉普森方法編寫計算的程序。為收縮斷面的水深。解：由 整理得： （*）令， （*） （*）將（*）代入（*）式，以及對（*）式求導可得，由牛頓-拉普森切線法，可得收縮斷面水深的迭代公式為收斂判別式：。對于本題，也可以利用簡單迭代法求解收縮斷面的水深，其迭代公式為程序1.2：牛頓-拉普森切線法求解收縮斷面水深表1.2 變量說明變量名意義 Eo壩上游斷面的總比能q，x單寬流量，壩面流速系數(shù)A，常數(shù)B，常數(shù)C，常數(shù)linjieH引入的函數(shù)Hc下游收縮斷面水深Public Function linjieH(x!, q!, hc!) 建立

12、函數(shù)Dim h1!, A!, B!, C!, fhc!, fh1c!g = 9.8Eo = Val(InputBox(請輸入斷面比能Eo)X = Val(InputBox(請輸入流速系數(shù)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 / fh1cIf Abs(hc2 -

13、 hc1) e Then Exit Dohc = hc2LoopEnd FunctionPrivate Sub Command1_Click()Call linjieH(x!, q!, hc!) 調用函數(shù)Print Format(hc, 0.000)End Sub3二等分法求方程的根。選取區(qū)間a,b使，則其根在區(qū)間a,b內。12繼續(xù)二等分。每二等分一次區(qū)間縮短為原來的一半。第k次二等分后區(qū)間長度為：收斂判別式：計算步驟：（1）選??；使；（2）令，求；（3）若，則；若，則；（4）若，迭代結束，（或或）； 若，返回（2）繼續(xù)進行二等分。例1-3梯形渠道底寬，求。解： ， 選取。其中。程序：變量說明

14、表：程序中：意義：Q，e流量，允許誤差Ac，BcHc臨界水深hc1 ，hc2區(qū)間（a,b）倆端點水深 hc3區(qū)間（a,b）的中點值b,m渠底寬度，邊坡系數(shù)Dim b!, m!, Q!, hc!, g!, e!, Ac!, Bc!, hc1!, x!, hc2!, f1!, f2!, f3!, hc3!, p!Public Function LinjieH(b!, m!, Q!, hc!)g = 9.8Ac! = (b + m * hc) * hcBc! = b + 2 * m * hcp = Ac 3 / Bc - Q 2 / gEnd FunctionPrivate Sub Command1

15、_Click()g = 9.8b! = Val(InputBox(請輸入矩形斷面的底寬b)m! = Val(InputBox(請輸入邊坡系數(shù)m)Q! = Val(InputBox(請輸入流量Q)e! = Val(InputBox(請輸入計算精度e)hc1! = 0.01x! = 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!)If f1 * f3 0 Thenf

16、2 = f3: hc2 = hc3Elsef1 = f3: hc1 = hc3End IfLoop Until Abs(f3) eCall LinjieH(b!, m!, Q!, hc!)Print hc3End Sub第2節(jié)線性方程組的解法1高斯-賽德爾迭代法三元一次方程組：初值：雅可比迭代法：賽德爾迭代法：第步：雅可比迭代法：賽德爾迭代法：收斂判別式：對于計算步驟：1給定初值：2迭代計算：（1）雅可比迭代法：（2）賽德爾迭代法：3收斂判別式：賽德爾迭代法程序（4）變量說明表：程序中：意義：a(i, j)方程組系數(shù)矩陣元素b(i)方程組右端項元素x0(i) ， x1(i)方程組待求的未知量N

17、矩陣的維數(shù)Dim i As Integer, j As Integer, n As Integer, sum As Integer, s1 As Double, s2 As DoubleDim a() As Integer, b() As Integer, x0() As Double, x1() As Double, EPS As Double, t#, z#Public Function GaussSeidel()n = InputBox(輸入矩陣的維數(shù)：)ReDim Preserve a(n, n), b(n), x0(n), x1(n)For i = 1 To nFor j = 1 T

18、o na(i, j) = InputBox(輸入a(i,j)的值：, a( & i & , & j & )Print a(i, j);Next jPrintNext iFor i = 1 To nb(i) = InputBox(輸入b(i)的值：, b( & i & )x0(i) = 0Print b(i);Next iPrint此處判斷該系數(shù)矩陣是否為嚴格對角占優(yōu)陣,若不想判斷也可將下10行刪除For i = 1 To nFor j = 1 To nt = Abs(a(i, i)z = z + Abs(a(i, j)Next jIf t z - Abs(a(i, i) ThenPrint 您

19、輸入的方程組的Gauss-Seidel迭代式不收斂！;End IfIf t = z - Abs(a(i, i) Then GoTo FNext iF: DoFor i = 1 To nFor j = i + 1 To ns2 = s2 + a(i, j) * x0(j)Next jFor j = 1 To i - 1s1 = s1 + a(i, j) * x1(j)Next jx1(i) = (b(i) - s1 - s2) / a(i, i)EPS = Abs(x1(1) - x0(1)s1 = 0s2 = 0Next iFor i = 1 To nx0(i) = x1(i)Next isu

20、m = sum + 1If EPS = 0. Then Exit DoLoopPrint 迭代了; sum; 次Print x(1)=; Format(x0(1), 0.000); x(2)=; Format(x0(2), 0.000); x(3)=; Format(x0(3), 0.000);聲明:此程序是系數(shù)矩陣滿足嚴格對角占優(yōu)陣的情況下,Gauss-Seidel迭代法收斂;對于滿足不可約對角占優(yōu)陣的Gauss-Seidel迭代法收斂情況下,有興趣的讀者可查看有關資料,自己動手試試 3x1+x2+6x3=3 x1+x2+2x3=3 -x2+2x2-3x3=1 x1=-10 x2=3 x3=

21、5End FunctionPrivate Sub Form_Click()Call GaussSeidelEnd Sub雅可比迭代法 程序Dim i As Integer, j As Integer, n As Integer, sum As Integer, s1 As Double, s2 As Double, t#, z#Dim a() As Integer, b() As Integer, x0() As Double, x1() As Double, EPS As DoublePublic Function Jacobi()n = InputBox(輸入矩陣的維數(shù)：)ReDim Pr

22、eserve a(n, n), b(n), x0(n), x1(n)For i = 1 To nFor j = 1 To na(i, j) = InputBox(輸入a(i,j)的值：, a( & i & , & j & )Print a(i, j);Next jPrintNext iFor i = 1 To nb(i) = InputBox(輸入b(i)的值：, b( & i & )x0(i) = 0Print b(i);Next iPrintFor i = 1 To nFor j = 1 To nt = Abs(a(i, i)z = z + Abs(a(i, j)Next jIf t z

23、- Abs(a(i, i) ThenPrint 您輸入的方程組的Jacobi迭代式不收斂！;End IfIf t = z - Abs(a(i, i) Then GoTo FNext iF: DoFor i = 1 To nFor j = i + 1 To ns2 = s2 + a(i, j) * x0(j)Next jFor j = 1 To i - 1s1 = s1 + a(i, j) * x0(j)Next jx1(i) = (b(i) - s1 - s2) / a(i, i)EPS = Abs(x1(1) - x0(1)s1 = 0s2 = 0Next iFor i = 1 To nx0

24、(i) = x1(i)Next isum = sum + 1If EPS maxV ThenmaxV = Abs(a(i, j) 將主對角線上的最大元素找出maxI = i row of max elementmaxJ = j column of max elementEnd IfNext jNext iIf ii maxI Then exchange rowsd1 = b(maxI)b(maxI) = b(ii)b(ii) = d1For j = ii To nd1 = a(maxI, j)a(maxI, j) = a(ii, j)a(ii, j) = d1Next jEnd IfIf ii

25、 maxJ Then exchange columnsi1 = index(maxJ)index(maxJ) = index(ii)index(ii) = i1For i = 1 To nd1 = a(i, maxJ)a(i, maxJ) = a(i, ii)a(i, ii) = d1Next iEnd IfIf a(ii, ii) = 0# ThenExit FunctionEnd IfFor i = ii + 1 To n eleminationd1 = a(i, ii) / a(ii, ii)For j = ii + 1 To na(i, j) = a(i, j) - a(ii, j)

26、* d1Next jb(i) = b(i) - b(ii) * d1Next iNext iix(n) = b(n) / a(n, n) back substitutionFor i = n - 1 To 1 Step -1d1 = 0#For j = i + 1 To nd1 = d1 + a(i, j) * x(j)Next jx(i) = (b(i) - d1) / a(i, i)Next iFor i = 1 To nb(index(i) = x(i)Next iSolve = TrueEnd Function第三節(jié) 插值例 已知消力坎的淹沒系數(shù)與相對淹沒度對應關系如下表：.0.700

27、.720.740.760.78.0.940.930.9130.900.885.hs為下游水位超過坎頂?shù)母叨?；H0為坎上全水頭。試編一分別用線性插值和拋物型插值的程序，計算相對淹沒度時的淹沒系數(shù)值。1線性插值已知，求，要求滿足則 令2拉格朗日插值二次插值（拋物型插值），求。設滿足條件： 3n次插值（拉格朗日插值），求 線性插值程序變量說明表：程序中：意義：x(i)已知節(jié)點坐標y(i)已知節(jié)點函數(shù)值m待求函數(shù)值的節(jié)點坐標Dim m!, f1!, f2!, L1x!, L2x!, Lx!Private Sub Command1_Click()Dim x(4) As Single, y(4) As S

28、inglex(0) = 0.7: x(1) = 0.72: x(2) = 0.74: x(3) = 0.76: x(4) = 0.78y(0) = 0.94: y(1) = 0.93: y(2) = 0.915: y(3) = 0.9: y(4) = 0.885m = Val(InputBox(請輸入一個數(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 =

29、 i + 1Loop Until f1 * f2 0Print Format(Lx, 0.0000)End Sub第4節(jié)擬合堰流： 最小二乘法：偏差若要使取最小值，則 例：直角三角形薄壁堰，求時， 0.71 0.86 0.90 0.97 1.10 1.25 1.37 1.44 1.53 1.64 1.77 1.871.88 3.13 3.38 3.88 5.14 8.05 9.22 10.50 12.30 14.33 16.85 18.80最小二乘法線性擬合直角三角形薄璧堰流量的程序變量說明表：程序中：意義：H直角三角形薄壁堰的堰上水頭t1(i)，t2(i)堰上水頭，流量，一維數(shù)組sum1，s

30、um2公式中AA1，AA2，V中間變量Q欲求流量Private Sub Command1_Click()Dim H#H = Val(InputBox(請輸入直角三角形薄壁堰的堰上水頭)Call nehe(H#)End SubPublic Function nehe(H#)Dim t1(8) As Double, t2(8) As Double, x(8) As Double, y(8) As Double, sum1#, sum2#, p1#, p2#, AA1#, AA2#, V#, r1#, r2# t1(0) = 0.71: t1(1) = 0.86: t1(2) = 0.9: t1(3

31、) = 0.97: t1(4) = 1.1: t1(5) = 1.25: t1(6) = 1.37: t1(7) = 1.44: t1(8) = 1.53 t2(0) = 1.88: t2(1) = 3.13: t2(2) = 3.38: t2(3) = 3.88: t2(4) = 5.14: t2(5) = 8.05: t2(6) = 9.22: t2(7) = 10.5: t2(8) = 12.3 For i = 0 To 8 y(i) = Log(100 * t2(i) Next i For i = 0 To 8 x(i) = Log(10 * t1(i) Next i For i =

32、0 To 8 AA1 = x(i) AA2 = y(i) sum1 = sum1 + AA1 sum2 = sum2 + AA2 Next i V = sum2 - sum1 For i = 0 To 8 r1 = x(i) * y(i) * V * 9 r2 = x(i) 2 * V * 9 p1 = p1 + r1 p2 = p2 + r2 Next i m = p1 / p2 b = Abs(sum2 - m * sum1) / 9 Q = Exp(b + m * Log(H * 10) Print QEnd Function第二章 偏微分方程式的差分數(shù)值解法第1節(jié) 分類、數(shù)值解法和差分

33、格式1分類二階線性偏微分方程式：典型例子：橢圓型偏微分方程式：拉普拉斯方程式 拋物型偏微分方程式：擴散方程式 雙曲型偏微分方程式：波動方程式 2數(shù)值解法主要步驟：劃分網(wǎng)格；構造差分方程式；求解差分方程式。3差分格式函數(shù)及在處的泰勒級數(shù)展開式分別為： （1）向前差分（一階精度）（2）向后差分（一階精度）（1）-（2）得：中心差分（二階精度）（1）+（2）得：中心差分（二階精度）第2節(jié) 橢圓型偏微分方程式的數(shù)值解法基本方程式：（在域A內）邊界條件：（1）第一類邊界條件：（2）第二類邊界條件：（3）第三類邊界條件：有限差分解法：（1）網(wǎng)格劃分：（2）構造差分方程式： ， 五點格式：（3）解差分方程式

34、1高斯迭代法2高斯-賽德爾迭代法3超松弛迭代法-超松弛因子（=1時為高斯-賽德爾迭代法，1為亞松馳迭代，1為超松弛迭代）（計算試驗表明：取1.21.5之間數(shù)值時迭代次數(shù)最?。┡袆e收斂性： 例：如圖所示為夾在兩個無限平行壁間的圓柱繞流，無限遠處的來流速度分布均勻，且，二平板間通過的流量為，圓柱半徑，試用有限差分法編寫求流函數(shù)在流場中分布的程序。解： （一）網(wǎng)格劃分：，。（二）構造差分方程式： （在域A內）第一類邊界條件：（1）；（2）；（3）；第二類邊界條件：（4）。五點差分格式：超松弛迭代法：（=1.21.5）邊界條件：（1）（ab邊界）；（2）；（3）bc邊界的處理： (a) (b) ；（4

35、）；（ae邊界）；（5）cd邊界的處理: 的計算：(a)直接轉移：(b)直線內插：如上圖所示，。則，其中a表示P1與網(wǎng)格線和圓柱壁面交點之間的距離，是（7，3）和（7，4）兩網(wǎng)格點之間的距離即y方向的網(wǎng)格步長；b表示P2與j=2網(wǎng)格線和圓柱壁面交點之間的距離，是（6，2）和（5，2）兩網(wǎng)格點之間的距離即x方向的網(wǎng)格步長。（三）解差分方程式：逐點迭代計算。每循環(huán)一次，進行收斂條件的判別。收斂判別條件：有限差分法計算兩個無限平行壁間圓柱繞流流函數(shù)在流場中分布的程序變量說明表：程序中：公式中：意義：psi(i, j)節(jié)點的流函數(shù)值，二維數(shù)組w松弛因子dij，sdij0，sdij1Private Su

36、b Command1_Click()Dim psi(8, 5), w, sdij0, sdij1For j = 1 To 5 For i = 1 To 8 psi(i, j) = 0 Next iNext jFor j = 2 To 4 psi(1, j) = (j - 1) * 0.5Next jFor i = 2 To 8 psi(i, 5) = 2Next iDow = 1.35k = 0sdij0 = 0 k = k + 1sdij1 = 0psi(6, 2) = psi(5, 2) * 0.125 / 0.625psi(7, 3) = psi(7, 4) * 0.125 / 0.62

37、5psi(8, 4) = (2 * psi(7, 4) + 2) / 4For j = 2 To 4 For i = 2 To j + 3 dij = 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) + dij Next iNext jIf Abs(sdij1 - sdij0) 0.00001 Then Exit Dosdij0 = sdij1LoopFor j = 1 To 5Pri

38、ntFor i = 1 To 8 Print Space(1), Format(psi(i, j), 0.000); Next iNext jEnd Sub第3節(jié) 拋物型偏微分方程式的數(shù)值解法例：兩無限平板間油的流動，上平板水平運動，0.1s內u=0線性增加到，其后保持不變。求兩平板間的速度分布隨時間的變化。解：一基本方程式 邊界條件： 初始條件：二有限差分解法（一）網(wǎng)格劃分 取空間步長 其中；取時間步長 其中。（二）構造差分方程式（1） 顯式格式：若令，則 穩(wěn)定性條件：（2） 隱式格式：若令，則 隱式格式無條件穩(wěn)定。初始條件：邊界條件：（三）解差分方程式顯格式差分法計算兩平板間速度分布隨時間

39、變化的程序變量說明表：程序中：公式中：意義：n，Hn，a/m時段數(shù)，距離步長MM距離a被等份的份數(shù)AA兩平板間的距離V液體的運動粘滯系數(shù)U0平板的最終速度RR顯格式的穩(wěn)定系數(shù)u(i, j)點流速Tj中間變量DTt時間步長Private Sub Command1_Click()Dim m, 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 = 0.1: count = 0DT = R * H 2 / vFor i = 1 To m +

40、1u(i, 1) = 0Next iFor j = 1 To n + 1u(1, j) = 0t(j) = (j - 1) * DTtj = t(j)If tj 0.1 Thenu(m + 1, j) = tj * U0 / 0.1Elseu(m + 1, j) = U0End IfNext jFor j = 1 To nFor i = 2 To mu(i, j + 1) = R * u(i + 1, j) + (1 - 2 * R) * u(i, j) + R * u(i - 1, j)Next iNext jFor j = n + 1 To 1 Step -50Print jFor i =

41、 1 To m + 1count = count + 1Print Space(1), Format(u(i, j), 0.00);If count = 11 Thencount = 0: PrintEnd IfNext iNext jEnd Sub第3章 有壓管流第1節(jié)管網(wǎng)的水力計算-哈迪克勞斯法例1：如圖所示環(huán)狀管網(wǎng)，假設。試求：（1）各管段的流量分配；（2）各節(jié)點的測壓管水頭。管段12345678910管長600600300600300900300600300600管徑0.40.40.30.30.20.20.30.30.30.3沿程阻力系數(shù)0.020.020.030.030.040.04

42、0.030.030.030.03解：（一）主要公式連續(xù)性方程：。n為該點的管段數(shù)目，流進節(jié)點，流出節(jié)點。能量方程式：。環(huán)的正方向與該管段流向相同時，否則 ， ， 所以：（二）（1）假設各管段的流量分配，使其滿足連續(xù)性方程；（2）A環(huán)路水頭殘差；設的修正流量為，使；即：所以： ，定性分析其正確性：如果初分流量時只是數(shù)值誤差，不是方向“誤差”的話，分母總為正值。則若分子（即殘差）為正，需減小。亦即正的水頭損失太大了，需減小流量；而負的水頭損失太小了，需增加流量。這是正確的。反之亦然。（3）計算步驟：1根據(jù)各管段的計算；2任意選擇各管段中的，滿足連續(xù)性方程（給定各管段的流動方向）；3選定各閉合環(huán)路的

43、正流動方向（如順時針方向）來確定各管段的，并構成矩陣，k-環(huán)路編號（行），I-管段編號（列）。對于環(huán)路中不包含的管段，矩陣元素。當該管段流向與該環(huán)方向相同時，；當該管段流向與該環(huán)方向相反時，。矩陣：（第k閉合環(huán)路- ） ik12345678910 1-1+100000-10+1200+1+10000-1-13 0000+1-1-1+1+104.計算（對每個閉合環(huán)路），擴展為：。不在該環(huán)時，；在該環(huán)時。5計算各管段中的流量，擴展為：。對于非共用管，只有一個不為零；對于平面共用管，只有兩個不為零，并且當均規(guī)定順時針方向為環(huán)路正方向時，一個為“+1”，一個為“-1”。6判別計算誤差（如的允許誤差），

44、是否滿足：（其中為總環(huán)數(shù)），若不滿足，返回第4步繼續(xù)進行各環(huán)路修正流量的計算。7計算各節(jié)點的測壓管水頭（需構造節(jié)點及管段號之間的關系數(shù)組）。幾點解釋：（1）只管數(shù)值，不管方向。流量方向事先已假定，最終只要不修正為負值，方向就與原假設方向相同。如若修正為負值，則與原假設方向相反。（2）連續(xù)性方程中流量方向用符號函數(shù)表示。流入為“+1”，流出為“-1”。（3）能量方程式中只表示修正流量的數(shù)值。為正值時，即為均按順時針方向增加流量；為負值時，即為均按逆時針方向增加流量。（4）的意義為：；。哈迪-克勞斯法管網(wǎng)的水力計算的程序變量說明表：程序中：公式中：意義：d(i)，l(i)d,l各管段的管徑，管長y

45、(i)各管段的沿程阻力系數(shù)R(i)R公式中的系數(shù)Q(i，j)Q第j次迭代第i個管段的流量M(i，j)表示管段中的流動方向是否與環(huán)路方向一致的符號矩陣H(i)H各管段的水頭B, A中間變量c(i)，x(i)中間變量Private Sub Form_Click()Dim d(1 To 10), l(1 To 10) As Integer, y(1 To 10), R(1 To 10), Q(1 To 10, 100) As Double, w#, M(1 To 3, 1 To 10)Dim Q1(1 To 3) As Double, B, A, c(1 To 10), H(1 To 8), x(1 To 10) 輸入已知數(shù)據(jù)d(1) = 0.4: d(2) = 0.4: d(3) = 0.3: d(4) = 0.3: d(5) = 0.2: d(6) = 0.2: d(7) = 0.3: d(8) = 0.3: d(9) = 0.3: d(10) = 0.3 l(1) = 600: l