河口海岸水動力模擬技術講義PPT課件

1、第一章第一章 緒論緒論v海岸：是海陸相互作用的重要地帶，也是海、陸、氣交互作用的重要空間，這種表現(xiàn)在： 岸線演變（自然和人為） 颶風（臺風）帶來的災難性破壞； 海洋潮汐環(huán)境的變化。v河口：海岸常伴隨有江河湖泊的出海口，通常稱為河口。第1頁/共173頁v海岸河口問題： 潮流問題 波浪問題 徑流、異重流（密度流）、污染物（COD)擴散。v研究海岸河口問題的方法 物理模型（水力學比尺模型） 數(shù)學模型（數(shù)值模擬）第2頁/共173頁沿岸過程動力因素物質過程流（潮流）波（風浪）鹽水入侵泥沙輸移污染物擴散波流相互作用海水入侵控制反饋流載波波生流第3頁/共173頁v數(shù)值模擬：一門綜合性的模擬技術，它采用數(shù)學模

2、型來模擬某中物理現(xiàn)象，并通過計算機用數(shù)值計算法進行近似求解，籍以復演自然演變過程的總稱。v水力學、泥沙數(shù)值模擬：以水力學和泥沙動力學為理論基礎，并結合具體工程的一門新型實用科學。第4頁/共173頁v水動力泥沙數(shù)值模擬：以微分方程為理論，并通過微分方程的離散，變成代數(shù)方程，最后采用計算機進行近似求解。v數(shù)值模擬的特點： （1）一般以線性理論為基礎，但實際自然現(xiàn)象和描述這些現(xiàn)象的微分方程均為非線性的； （2）需要豐富的經(jīng)驗，現(xiàn)場資料和一定的技巧； （3）數(shù)值模擬不僅僅是一種近似計算，可以作為一種實驗或研究及預測方法。第5頁/共173頁v數(shù)值模擬的優(yōu)點： （1）實驗費用少； （2）速度快、周期短；

3、（3）可以模擬多種因素相互作用的復雜物理過程。如可以模擬水（潮）流、風、柯氏力等多種因素共同作用下的多種泥沙及地形演變的復雜過程。 （4）可以完全控制流體的物理性質（如密度、容重、粘度、含沙量等） （5）模型建成后，長期保存、隨時調用修改。 （6）無法模擬微分方程不能描述的物理現(xiàn)象。第6頁/共173頁v數(shù)值模擬工作的基本步驟（1）建立數(shù)學模型和編制源程序 建立或選擇的微分方程； 根據(jù)模擬域邊界條件選擇合適的網(wǎng)格； 按一定的格式離散方程，得到代數(shù)方程和采用合適的數(shù)值方法求解代數(shù)方程； 編制源程序求解代數(shù)方程。 數(shù)值模擬分析（收斂性、穩(wěn)定性、相容性、誤差程度等）第7頁/共173頁（2）調試源程序（

4、3）模型驗證 調整模型中有關參數(shù)（糙率、紊動動量摻混系數(shù)等），使模型有良好的穩(wěn)定性和收斂性，并與現(xiàn)場資料有良好的吻合；（4）正式方案試驗 第8頁/共173頁v河口、海岸水動力模擬的發(fā)展方向 1、資料同化將是河口數(shù)值模型發(fā)展和結合的一個新技術切入點,也是帶動河口動力數(shù)字模擬技術革新的一種重要方法.2、數(shù)字河口動力模型第9頁/共173頁1、河口模型四維資料同化v四維同化在河口中的一類主要應用是為河口數(shù)值模式提供優(yōu)化的初始狀態(tài),從而校正不合理的邊界條件所帶來的偏差,提高模擬精度.v變分同化技術還可應用于確定河口模式中的未知參數(shù). 例如確定了潮汐河口的摩擦系數(shù),通常這種系數(shù)是通過經(jīng)驗獲得的.v另外,河

5、口模式的外部強迫場(如風場、海-氣界面的熱通量等)都可以通過這種技術反演出來.第10頁/共173頁v資料同化的關鍵在同化模型在過濾和插值觀測資料的機制方面有高效的預測能力。模型應對前期信息在時間上能向前插值。v動力場的模型同化熟路也被證明對陸-氣系統(tǒng)內部環(huán)流的科學研究極有價值v同化是唯一能產(chǎn)生復雜非線性過長的場的手段，并且在物理學與動力學的機制上協(xié)調一致。第11頁/共173頁2、數(shù)字河口動力模型v河口水動力特征及外界強迫作用因子構建數(shù)字河口動力模型(模型計算要素象溫、鹽、流、泥沙的特征均以數(shù)字矩陣形式記載),若再與數(shù)字河網(wǎng)模型嵌套聯(lián)結,最終可以獲得區(qū)域河網(wǎng)河口的動力-沉積-地貌的機制解釋,揭示

6、河口物理規(guī)律和解決工程實際問題.v在數(shù)字河口模型的構架下,可以將人類已經(jīng)擁有的河口科學理論、知識與具有較強物理概念的水動力學模型集成于一體,為河口數(shù)值模型計算提供良好的平臺。第12頁/共173頁數(shù)字河口動力模型具有許多優(yōu)勢:首先,數(shù)字河口模型是基于數(shù)字區(qū)域地形構建而成的,地形要素可自動生成,無需手工操作,大大提高了工作效率;其次,數(shù)字模型不僅能輸出傳統(tǒng)模型的結果,而且能夠十分方便地給出河口水文要素和水文狀態(tài)變量的空間分布場,這些對近岸河口動力科學研究與河口、港口、航道工程都有著廣闊的應用前景.數(shù)字河口模型研究的最終目的就是利用已有的河口基礎科學理論和知識,在數(shù)字區(qū)域地形的基礎之上將觀測點的水文

7、信息拓展、同化至區(qū)域平面上乃至區(qū)域三維立體上的信息,并形成數(shù)字成品。第13頁/共173頁v參考文獻：vKoutitar 著“Mathematical Model in Coastal Engineering”1)模型簡單易懂2）附有Basic程序，而且有驗證的算例3）介紹各種數(shù)值處理技術v曹祖德、王運洪”水動力泥沙數(shù)值模擬第14頁/共173頁第二章第二章 水動力數(shù)值模擬的理論基礎水動力數(shù)值模擬的理論基礎)()()()(022222222zVAzyVxVAygfUdtDVzUAzyUxUAxgfVdtDUzWyVxUzHzH ),()tyxzyVxUtW2.1 基本方程自由面運動學邊界條件：底部

8、運動學邊界條件：),()(yxhzyhVxhUW第15頁/共173頁zWyVxUtdtDU,V,W為x,y,z 方向上的流速分量。（，）為距平均海平面的自由表面水位。（，）為平均海平面距底部邊界的水深。為水平擴散系數(shù)。為垂直渦動系數(shù)。第16頁/共173頁0)()()()(222222yQxQtQDCVUgygDfQQVyQUxtQQDCVUgxgDfQQVyQUxtQyxyxyyyxyxxx初始條件),()0 ,(0yxuyxu),()0 ,(0yxvyxv),()0 ,(0yxyx邊界條件岸邊界：法向流速為零。水邊界：給定潮位過程。第17頁/共173頁0)()()(0oxbxbbBxAZhg

9、xZhgAxuQtQxQtASaint Venant 方程A 過水斷面面積，Q為流量，Q=Bhu，B為河寬，h為水深，u為斷面平均流速。 分別為水面剪切力與水底剪切力Zb 為水底的豎向坐標位置X軸為沿水流縱向方向。第18頁/共173頁三、二、一維方程的定解條件三、二、一維方程的定解條件v初始條件u,v,w,|t=0=u0,v0,w0,0邊界條件開邊界：計算域水體與外部水體相接處。(u,v,w)=(u(t),v(t),w(t)=(t)固邊界：計算域與陸地或建筑物接壤處無滑動：u,v,w=0有滑動： 垂直邊界的速度為0。0nV第19頁/共173頁2.22.2數(shù)值計算數(shù)值計算v在計算水動力、泥沙數(shù)值

10、模擬時，大都將基本方程組離散成代數(shù)方程組，最后求解代數(shù)方程組，此處介紹微分方程組的離散技術有限差分法和線性代數(shù)方程組的數(shù)值解法。第20頁/共173頁有限差分法有限差分法v有限差分法是工程中常用的一種離散技術，將計算域分成有限個網(wǎng)格，通過差分法求網(wǎng)格結點的微分方程的近似值，也稱網(wǎng)格法。v將網(wǎng)格結點上的函數(shù)f(x,y,z,t)表示成 ，i,j,k分別表示x,y,z方向的坐標位置，n表示時間。nkjif,第21頁/共173頁1、工程中常用的幾種差分和微分的關系（一維）（1）一階向前差分10 ,)(2)()()(1122xxxxxfxxxfxxfxxf10 ,)(2)()()(2222xxxxxfxx

11、xxfxfxxf（2）一階向后差分第22頁/共173頁（3）一階中心差分10 , 2/10 , 2/)()(48)2()2()(4443333433332xxxxxxxxfxxfxxxxfxxfxxf（4）二階中心差分10 , 2/10 , 2/)()(96)2/()(2)2()2()(6665554644542222xxxxxxxxfxxfxxxfxxfxxfxxf第23頁/共173頁2、幾種常見的差分格式以一維熱傳導方程為例：022xftfxti-1 i i+1n+1 n第24頁/共173頁（1）古典顯式格式0)(222111xtOxffftffninininini022111xffftf

12、fninininini2111/),()21 (xtrffrfrfnininini可由已知值直接求解未知值，且只涉及二層，稱為雙層顯式差分格式第25頁/共173頁（2）古典隱式格式0)(222111111xtOxffftffninininini022111111xffftffninininini21111/,)21 (xtrfrffrrfnininini不能直接求解未知值，且只涉及二層，稱為雙層隱式差分格式第26頁/共173頁（3）六點格式（Crank-Nicolson），雙層六點隱式格式在x點和n+n/2時層，對t和x均采用中心差分0)(2221222112111111xtOxfffxfff

13、tffnininininininini21111111/,22)1 (2/)1 (2/xtrfrfrfrfrfrfrnininininini022212112111111xfffxffftffnininininininini第27頁/共173頁（4）Richardson格式，三層顯式格式在x點和n時層，對t和x均采用中心差分0)(222221111xtOxffftffninininini21111/),2(2xtrfffrffninininini02221111xffftffninininini第28頁/共173頁（5）加權六點格式，隱式格式在x點和n+n時層，01，對t和x均采用中心差分2/

14、10)(2)1 (222112111111xtOxfffxffftffnininininininini21111111/)1 ()1 ()1 (21 )21 (xtrrfrffrrffrrfnininininini02)1 (22112111111xfffxffftffnininininininini第29頁/共173頁線性方程組的數(shù)值解線性方程組的數(shù)值解v有限差分法是工程中常用的一種離散技術，將計算域分成有限個網(wǎng)格，通過差分法求網(wǎng)格結點的微分方程的近似值，也稱網(wǎng)格法。v將網(wǎng)格結點上的函數(shù)f(x,y,z,t)表示成 ，i,j,k分別表示x,y,z方向的坐標位置，n表示時間。nkjif,第30頁

15、/共173頁1、解線性方程組的兩種方法:v直接法:通過有限步算術運算直接求出方程組的精確解,最常用的是消元結合代入的方法.v實際上除非是采用無窮位精度計算,一般都得不到精確解.v直接方法適用于解低階稠密矩陣方程組.第31頁/共173頁v迭代法 類似于方程求根的迭代法,用一個迭代過程逐步逼近方程組的解.v迭代有可能不收斂,或雖然收斂,但收斂速度慢.v迭代法適用于求解高階稀疏矩陣方程組.v稀疏矩陣:矩陣非零元素較少,且在固定的位置上.v稀疏矩陣一般是人為構造的,例如36頁三轉角插值時方程組(8.12),(8.15)的系數(shù)矩陣.第32頁/共173頁GaussGauss消去法消去法( (第一次消元第一

16、次消元) )v考慮方程組A(1)x=b(1)11111111221111112112222211111122( )( )( )( )( )( )( )( )( )( )( )( )nnnnnnnnnnaxaxaxbaxaxaxbaxaxaxb v第一次消元用第一個方程將后面方程的x1消去.111111( )( )iiama v計算乘數(shù)v條件:a11(1)0v用-mi1乘以第一個方程加到第i個(i=1,n)方程上,則消去了第i個方程中的x1.第33頁/共173頁GaussGauss消去法消去法( (第一次消元第一次消元) )v經(jīng)過上述過程,得到方程組A(2)x=b(2),111111121112

17、2222222222200( )( )( )( )( )( )( )( )( )( )nnnnnnnaaabxxaabxbaa v其中2111 1( )( )( ),ijijijaam a2111 1( )( )( ),iiibbm b111111( )( )iiama 2 3( , , )i jn 第34頁/共173頁GaussGauss消去法消去法( (第第k k次消元次消元) )v假設已完成k-1次消元,得到方程組A(k)x=b(k).111111112111122222221221111111000000000( )( )( )( )( ),( )( )( )( ),()()(),(

18、)( )( )( )kknkknkkkkkkkknkkkkknkknknnaaaaaaaaaaaaaaaa v第k次消元的目的是將akk(k) (稱為主元)下面的元素變?yōu)?.第35頁/共173頁GaussGauss消去法消去法( (第第k次消元次消元) )v對A(k)右下角的矩陣( )( )( )( )kkkkknkknknnaaaa( )( )kikikkkkama v計算乘數(shù)v條件:akk(k)0v用-mik乘以第k個方程加到第i個(i=k+1,n)方程上,則消去了第i個方程中的xk,得到方程組A(k+1)x=b(k+1).第36頁/共173頁GaussGauss消去法消去法( (第第k次

19、消元次消元) )v第一步消元的計算公式2111 1( )( )( ),ijijijaam a2111 1( )( )( ),iiibbm b111111( )( )iiama v類似可以得到第k步消元的計算公式1()( )( ),kkkijijikkjaam a 1()( )( ),kkkiiikkbbm b ( )( )kikikkkkama 1( , )i jkn第37頁/共173頁GaussGauss消去法消去法v消去法完成后最終得到與原方程組等價的三角形方程組A(n)x=b(n).1111111211122222222000( )( )( )( )( )( )( )( )( )nnnn

20、nnnnaaabxxaabxba v一共需進行 ? 步n-1第38頁/共173頁GaussGauss消去法消去法( (算法算法) )1,1kn( )( )/kkikikkkmaa 1,ikn1()( )( )kkkijijikkjaam a 1,i jkn1()( )( )kkkiiikkbbm b 1,ikn( )( )/nnnnnnxba 1,1kn( )( )( )1()/nkkkkkkjjkkj kxbaxa 第39頁/共173頁追趕法求解三對角方程組追趕法求解三對角方程組v上面的方程組可以利用追趕法求解(P185).v對于下面形式的方程組11112222211111nnnnnnnnn

21、fbcxabcxfabcxfabxf v將系數(shù)矩陣進行三角分解第40頁/共173頁v比較兩邊對應元素可以得到11223111111nnnnrrr 11222111nnnnnbcabcabcab11,b 111,c 2(, ),iiar in12(, ),iiiibrin 21(,).iiicin 第41頁/共173頁v因此有11,b 111,c 2(, ),iiar in12(, ),iiiibrin 21(,).iiicin iiic 1iiiicbr 121(,)iiiicinba v又11111ccb v因此所有i的可遞推求出,進一步可求出i,ri.第42頁/共173頁v在得到系數(shù)矩陣的

22、分解后,原方程組轉化為vLUx=f.v先求解Ly=f11122223111nnnnnnnyfryfryfryf v顯然有y1=f1/1,vyi=(fi-riyi-1)/i=(fi-iyi-1)/(bi-aii-1)(i=2,n)第43頁/共173頁v再求解Ux=y,111221111111nnnnnyxxyxyxy v顯然有xn=yn, xi=yi-bixi+1(i=n-1,1)第44頁/共173頁迭代法迭代法v在處理一元方程f(x)=0時,我們將其轉化為x=j(x)的形式,然后用不動點迭代的方法進行求解.v對于線性方程組Ax=b,我們也可以將其轉化為類似的形式: x=Bx+f,v任取初始向量

23、x(0),令x(k+1)=Bx(k)+f(k=0,1,),則得到一個向量的序列x(k).v若該序列收斂于向量x*,對x(k+1)=Bx(k)+f 兩邊取極限得到x*=Bx*+f,即x*是方程組的解.第45頁/共173頁JacobiJacobi迭代法與迭代法與Gauss-SeidelGauss-Seidel迭代法迭代法v對于方程組1231231238322041133631236xxxxxxxxx v我們將其改寫為1232133121322081433111633612()()()xxxxxxxxx 第46頁/共173頁JacobiJacobi迭代法迭代法v寫成矩陣的形式為x=B0 x+f,其中

24、03208841011116301212B 20833113612f 第47頁/共173頁JacobiJacobi迭代法迭代法v利用x(k+1)=Bx(k)+f 進行迭代,得到結果如下kx1(k)x2(k)x3(k)000012.53.03.03.022.87500000 2.36363636 1.000000002.083.00020012 2.00063786 0.999830513.30e-393.00028157 1.99991182 0.999740487.26e-410 3.00003181 1.99987402 0.999881262.50e-41( )()|kkxx 第48頁/

25、共173頁JacobiJacobi迭代法迭代法v從上表可以看出,迭代序列逐步逼近方程組的精確解(3,2,1)T.v注:在迭代中,我們不可能得到x(k)和精確解之間的誤差,一般我們用|x(k)-x(k-1)|(通常用無窮范數(shù))的值來判斷是否終止迭代.v在上面的例子中,我們將第i個方程變形為左邊是xi,右邊是其它分量和常數(shù)的線性組合,然后進行迭代,這一方法稱位Jacobi迭代.第49頁/共173頁JacobiJacobi迭代法迭代法v一般的,對于方程組Ax=b,設A非奇異且aii0(i=1,2,n),將A改寫為A=D L U,其中112233nnaaDaa 2131321230000nnnaLaa

26、aaa 1213123230000nnnaaaaaUa 第50頁/共173頁JacobiJacobi迭代法迭代法v將方程組改寫為vDx=(L+U)x+bvx=D1(L+U)x+D1bv令B0=D1(L+U)(稱位Jacobi迭代矩陣),f=D1b,上式簡記為x=B0 x+f.v我們得到Jacobi迭代公式 x(k+1)=B0 x(k)+f.111()( )()nkkiiijjiijj ixba xa v寫成分量的形式為第51頁/共173頁Gauss-SeidelGauss-Seidel迭代法迭代法v在前面的例子中,我們計算x1(k+1),用的是第k步的x2,x3;112312131312132

27、2081433111633612()( )( )()( )( )()( )( )()()()kkkkkkkkkxxxxxxxxx v計算x2(k+1),用的是第k步的x1,x3,我們有理由認為已經(jīng)計算出的第k+1步的x1比第k步的“好”.因此,我們應該用第k+1步的x1和第k步的x3來計算x2.11()kx v類似地,我們也應該用新信息計算x3.11()kx 12()kx 第52頁/共173頁Gauss-SeidelGauss-Seidel迭代法迭代法v我們可以將上面一般的Jacobi迭代公式改寫為111()( )()nkkiiijjiijj ixba xa 111111()()( )()in

28、kkkiiijjijjiijj ixba xa xa v這一迭代方法稱為Gauss-Seidel迭代.第53頁/共173頁Gauss-SeidelGauss-Seidel迭代法迭代法( (算例算例) )v其Gauss-Seidel迭代公式為v對于方程組1231231238322041133631236xxxxxxxxx 1123112131113121322081433111633612()( )( )()()( )()()()()()()kkkkkkkkkxxxxxxxxx 第54頁/共173頁Gauss-SeidelGauss-Seidel迭代法迭代法( (算例算例) )v同樣取x(0)=

29、(0,0,0)T,迭代結果如下kx1(k)x2(k)x3(k)000012.50000000 2.09090909 1.227272732.522.97272727 2.02892562 1.004132230.47733.00981405 1.99680691 0.995891253.25e-242.99982978 1.99968838 1.000163029.98e-352.99984239 2.00007213 1.000060773.84e-41( )()|kkxx 第55頁/共173頁超松弛迭代超松弛迭代(SOR)(SOR)方法方法v沿著從xi(k)到xi (k+1) (G)的方向

30、再向前走,就得到超松弛迭代(SOR)方法.( )kix1()()ki Gx 11( )()()()kkii Gxx v假設已知第k步的迭代向量x(k)以及第k+1步迭代向量x(k+1)的前i1個分量已知,Gauss-Seidel迭代法取111111()()( )()inkkkiiijjijjiijj ixba xa xa 第56頁/共173頁超松弛迭代方法超松弛迭代方法v我們定義新的xi(k+1)為xi(k)與 的加權平均.1()kix 1111()( )()( )()( )()()kkkkkkiiiiiixxxxxx111( )()( )()inkkkiiijjijjiijj ixba xa

31、 xa v在=1時,上述方法就是Gauss-Seidel方法,1時稱為超松弛法(有時不管的范圍,統(tǒng)稱為超松弛方法).第57頁/共173頁超松弛迭代方法超松弛迭代方法( (算例算例) )v對于方程組123441111141111141111141xxxx v松弛方法迭代格式為(1)( )( )( )( )( )111234(1)( )(1)( )( )( )221234(1)( )(1)(1)( )( )331234(1)( )(1)(1)(1)( )441234(14)/4(14)/4(14)/4(14)/4kkkkkkkkkkkkkkkkkkkkkkkkxxxxxxxxxxxxxxxxxxx

32、xxxxx 第58頁/共173頁超松弛迭代方法超松弛迭代方法( (算例算例) )v取x(0)=0,=1.3,終止準則為|x(k)x(k1)|10-5.kx1(k)x2(k)x3(k)x4(k)000001-0.32500000-0.43062500-0.57057813-0.756016020.7562-0.79858622-0.88649937-0.94718783-0.953687310.47410-1.00000717-0.99999179-1.00000289-1.000001703.45e-511-0.99999667-1.00000287-0.99999954-0.99999919

33、1.11e-512-1.00000152-0.99999922-1.00000012-1.000000524.85e-61( )()|kkxx 第59頁/共173頁超松弛迭代方法超松弛迭代方法( (算例算例) )v我們來觀察松弛因子對收斂速度的影響.0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0步數(shù)301 156 104765947383126211.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0步數(shù)1712121518243555114*v步數(shù)表示|x(k)x(k1)|10-5時的迭代步數(shù),=2.0時,500步以內不收斂.第60頁/

34、共173頁超松弛迭代方法超松弛迭代方法( (矩陣表示矩陣表示) )v超松弛迭代格式可以寫為1111()( )()( )()inkkkkiiiijjijjiijj ixxba xa xa 111111()( )()( )()()inkkkkiiiiiiiijjijjjj ia xa xba xa x v用矩陣可以表示為111()( )()( )()()kkkkDxDxbLxUx11()( )()()kkDL xDU xb 第61頁/共173頁第三章第三章 二維水動力數(shù)值模擬二維水動力數(shù)值模擬一、二維水動力數(shù)值模擬系統(tǒng)的分類1、按差分網(wǎng)格分：三角形、正方形、矩形、四邊形、多邊形、曲線坐標網(wǎng)格以及各

35、種形狀網(wǎng)格的組合2、按計算方法分：顯式法、隱式法、顯隱混合法3、按模擬格式分：三角元法、ADI法、破開算子法、單元體積法、MADI法、準分析法、貼體坐標法。第62頁/共173頁二、平面二維水動力數(shù)學模型的一般形式000)()(00hhfuygyvvxvutvhhfvxgyuvxuutuyvhxuhtsybysxbx第63頁/共173頁wwwassywwwassxbybxvvufuvufhvvunhuvun22223/12223/1222第64頁/共173頁v定解條件0( , ,0)0( , ,0)0( , ,0)( , )0u x yv x yx yx y岸邊界：法向流速為 ；水邊界：給定潮位

36、過程；第65頁/共173頁三、三、ADIADI法法v交叉方向隱式格式是對x,y方向交叉使用隱式格式（同時也交叉使用了顯式格式），使得求解過程簡化。222221/221/21211/221/2112212222u()()/ 2()/ 2nnnnnnnnuuatxyuuatuuatuxuy是的差分算子；是的差分算子；第66頁/共173頁1、網(wǎng)格正方形或矩形，變量u,v,分別交錯布置于網(wǎng)格中心和兩側。2、ADI基本思想（1）分步（2）交錯顯隱第67頁/共173頁ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位 、水深、水深uv第68頁/共173頁1/21/2,1/2,3/

37、2,1/2,1/2,11/2,11/21/2,1/2,1,1/2,1/2,1/2,1/2,0t/222( )bxsxnnnnnnijijijijijijnnijijnniji jbx ijnsxijijiuuuuvgfvtxyxhhuuuuuuuvxygfvxhh1/322222 1/32 1/321/21/2,1/2,3/2,1/2,1/2,11/2,11/21/2,1/2,21/2,1,1/2,0g=t/222jbxbxnnnnnnijijijijijijnnijijnnnijiji jnijnVuhnVun uv uuv uhh hh hChuuuuuuuvxyg uvgfvx21/21

38、/2,1/2,1/2,21/2,1/2,1/2,1/2,( )0()()nnij ijsx ijnnijijijijuC hh3、差分格式X向運動方程在（i+1/2,j)點離散第69頁/共173頁3、差分格式X向運動方程在（i+1/2,j)點離散1/21/21/2,1/2,1,221/2,1/2,3/2,1/2,201/2,1/2,11/2,11/2,;2()1();42();2()22nnnii jiijiijiinnijijnniijijnijinnijijniijabucdg taxuvttbuugxChg tcxuutdufy 其中）1/2,1,1/2,1/2,1/2,1/21,1/2

39、1,1/20001.1,1()41()21()2ijiji jnijnnnnniji ji jijijnnnijiji jvvvvvvhhh第70頁/共173頁001/21/21/2,01/2,1/2,1/2,01/2,1/2,1/2,0 ,1/2,1/2,1/20 ,1/2,1/2,1/2()()0()()/ 2()()0nnnnnni ji jijijijijijijnnnni ji ji ji ji ji jhuhvtxyhuhutxhvhvy連續(xù)方程在（i,j)點離散第71頁/共173頁1/21/21/21/2,1/2,01/2,01/2,1/20,1/2,1/20,1/2();21;

40、();2()()2nnniijii jiijiniijiniijnnnnnii ji ji ji ji jaubcudtahxbtchxtdvhvhy 其中第72頁/共173頁聯(lián)立上面的式子得到下面的線性方程組injiinjiinjiidcuba2/1, 12/1, 2/12/1, injiinjiinjiiducbua2/1, 2/12/1,2/1, 2/1第73頁/共173頁其中，u與存在如下關系iiiiiiiiiiiiiiiiiiiiiiiiiniiniiniinibEaFadHbEacGbGaHadFbGacEFuEHGu2/12/12/12/12/12/12/12/12/12/12/

41、12/12/1;其中第74頁/共173頁Y向運動方程在（i,j1/2)點離散)(412)(2)(4;)()()(2)(411, 2/11, 2/1, 2/1, 2/12/1,2/12/1,1,2/1, 12/1, 12/12/1,2/1,2/1,0222/1,22/12/1,2/1,2/3,2/12/1,jijijijijinjinjinjinjinjinjinjiinjinjinjinjinjiiinjiiuuuuuu ftytgvvxtuvMhCvugtvvytNMvN其中第75頁/共173頁v在后半個時間步長內，按上述同樣原理，在y向掃描。因此只須將上述各計算公式做如下的變換：xy,xy

42、;ij,ij; uv,uv;即可求出vi,j+1/2n+1, i,jn+1,然后顯式求ui+1/2,jn+1。第76頁/共173頁四、分步全隱式格式將x向動量方程與連續(xù)方程聯(lián)立，求出un+1,n+1/2,將y向動量方程與連續(xù)方程聯(lián)立，求出vn+1,n+1,第77頁/共173頁ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位水深水深uv第78頁/共173頁X向動量方程在（i+1/2,j)點離散0)()()()()()()()()(, 2/1, 2/1, 2/1, 2/12, 2/12, 2/12, 2/12/12, 2/12, 2/11, 2/12/1,2/1, 12

43、/1, 2/12/1, 2/1, 2/1, 11, 2/1, 2/11, 2/1njinjinjisxnjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihv fhCvuugxgyuuvxuuutuu 第79頁/共173頁連續(xù)方程在（i,j)點離散0)()()()(2/2/1,2/1,12/1,2/1,2/1,2/1, 2/1, 2/11, 2/1, 2/1, 2/11, 2/1,2/1,yhvhvxhuhutnjinjinjinjinjinjinjinjinjinjinjinjinjinji 第80頁/共173頁)(21)(1)(41)(21

44、)(21)(21)(21);(21, 1, 2/16/1, 2/1, 2/1, 2/12/1, 12/1, 12/1,2/1, 2/1, 2/11, 2/12/1, 2/11, 2/1, 2/12/1, 2/1, 2/11, 2/12/1, 2/1, 2/1, 2/1, 2/3, 2/1, 1njinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihnCvvvvvuuuuuuuuuuuuuuu 第81頁/共173頁injiinjiinjiidcuba2/1, 11, 2/12/1, injiin

45、jiinjiiducbua1, 2/12/1,1, 2/1 以上離散后的公式整理后如下:同理 在y方向，y向動量方程和連續(xù)方程聯(lián)立得如下：injiinjiinjiidcvba11,12/1,1, injiinjiinjiidvcbva12/1,1,12/1, 12/1,nnu 11,nnv 第82頁/共173頁五、移步雙向交替顯、隱式交錯法MADI法：將水深、水位、流速等變化量均布置在同一網(wǎng)格節(jié)點上，由此將基本方程離散成新的差分代數(shù)方程組，并建立一種新的解法，這種解法既吸收了原有傳統(tǒng)的ADI法的優(yōu)點，又有較高的穩(wěn)定性、收斂性和精度。第83頁/共173頁方程離散時，時間導數(shù)項采用前差表示，空間采

46、用中心差分，將t分成兩等分。v當nt(n+1/2)t,x向動量方程和連續(xù)方程建立差分方程iniiniiniiininiiniidcubadcbua2/112/12/112/112/12/11第84頁/共173頁2/112/11211111,00niniiiiiiiiiiiiiudduucbacba2/12/111211111,00niniiiiiiiiiiiiiudduucbacba第85頁/共173頁)()(2)(21)(2)2(22)()()(2)(4121,01,1,01, 10, 10,1,1,02,2,2, 1, 1njinjinjinjinjiinjiiinjiinjinjinji

47、njiiinjijinjinjinjinjiiihvhvytdhxtcbhxtavfyuutudxtgchCvutguuxtbxtga第86頁/共173頁)()()()()()(1 ()(1(1)(1 ()(1(22/1,02,2,22/1,1,1,) 1,() 1,(2/1, 1, 12/1,2/1,vsignKusignEhCvugvvvKvvKytNygfuvvEvvEuxtvMNMvvunjijinjinjinjinjinjivnjinjivjinjibynjibynjinjinjiunjinjiunjinjijijijinjiv當（n+1/2）t(n+1)t,y向動量方程和連續(xù)方程建

48、立差分方程,求解v,顯式求u，過程同上。第87頁/共173頁六、三角元法由于采用矩形網(wǎng)格而形成鋸齒形岸線以及鋸齒形堤、壩必然給數(shù)值模擬結果帶來不良影響，因此有限差分法的矩形網(wǎng)格，難以準確模擬不規(guī)則曲線形岸線，即使采用空間變步長網(wǎng)格，也不能完全模擬復雜的曲線形邊界。為了解決這些問題，采用了三角形網(wǎng)格。這種不規(guī)則三角形網(wǎng)格有以下優(yōu)點：v可以隨意加密計算網(wǎng)點形成不規(guī)則三角形，從而較準確地模擬出復雜的邊界，如：岸線、建筑物輪廓及航道，v可根據(jù)模擬區(qū)域的重要性，期F好網(wǎng)格的疏密程度及漸變程度。重要地區(qū)，網(wǎng)格密些，不重要地區(qū)，網(wǎng)格疏些；第88頁/共173頁000)()(hfuygyvvxvutvhfvxg

49、yuvxuutuyhvxhutbybx 基本方程第89頁/共173頁離散后方程0)()()()(0)()()()(0)()(11111111111nibynininininininininibxninininininininininininihfuygyvvxvutvvhfvxgyuvxuutuuyhvxhut 第90頁/共173頁iiiyaxaaF210yaxaayxF210),(jjjyaxaaF210kkkyaxaaF210第91頁/共173頁ekkkkjjjjiiiisFycxbaFycxbaFycxbayxF2/ )()()(),(ekkkkejjjjeiiiisycxbaNsycx

50、baNsycxbaN2/ )(2/ )(2/ )(kkjjiiFNFNFNyxF),(第92頁/共173頁kkjjiikkjjiiFcFcFcyFFbFbFbxFekkjjiiSFbFbFbMyF2/ )()(ekkjjiiSFcFcFcMxF2/ )()(eeSS第93頁/共173頁出于該格式采用顯式求解。同時受到復雜地形的影響，在計算過程中會產(chǎn)生一些擾。當這些擾動擴大并傳播易使計算失敗為消除這些擾動的影響，采用濾波公式mjmjjjiiiLLFFF1122*)()1( 第94頁/共173頁七、邊界處理邊界處理合適與否影響數(shù)值模擬的成敗。實際工程中，計算域的邊界常由不規(guī)則的曲線組成，如何利用有

51、限差分的矩形網(wǎng)格來模擬曲線邊界，以及如何選取邊界值，這是邊界處理的重要內容。第95頁/共173頁1、邊界類型（1）Dirichlet邊界（2）Neumann邊界（3）淺灘活動邊界第96頁/共173頁2、 Dirichlet邊界在邊界處有已知的函數(shù)值。在實際情況下，水邊界通常屬于這類邊界，該處的實測水位、流速，可作為已知函數(shù)值。1ff第97頁/共173頁v如果邊界網(wǎng)格結點正好在實測點上，則該結點的邊界值可直接采用已知值。v如果邊界網(wǎng)格結點不在實測點上，則分以下不同情況分別處理。v邊界結點與實測點很靠近，邊界值直接等于實測值；v如果實測點與邊界點較遠，但與虛擬外結點較近，則令虛擬外結點直接等于實測

52、值，然后引入邊界條件求出邊界值。第98頁/共173頁v邊界網(wǎng)格結點或虛擬結點與實測點均不靠近，不能直接引用實測值，這時，可根據(jù)前移時間步長所得的結果，用線性差值可得邊界值。xiicxjijixiicjijijixxfffxxffff )(, 1,)(, 1, 1第99頁/共173頁3、Neumann邊界在實際情況下，固定邊界屬于此類邊界0()Sfffn第100頁/共173頁v邊界正好通過網(wǎng)格或邊界與網(wǎng)格結點很靠近v邊界的外法線方向與坐標軸平行，直接得到邊界值；以A點位例：對于c點，可得第101頁/共173頁v邊界的外法線方向與坐標軸不平行，考慮外法線與坐標軸的夾角，帶入邊界條件后離散得到邊界值

53、。第102頁/共173頁v邊界與網(wǎng)格結點距離較大第103頁/共173頁4、淺灘活動邊界（1）開挖法：將灘地開挖至可能出現(xiàn)的最低水位之下，為使水量平衡，將岸邊界向水域內移動，并增大開挖部位的糙率以求得動量上的平衡。這種方法一般只適用于潮灘問題，而且灘地面積只占整個海區(qū)較小的情況，而對于水位變化小，坡度較緩的地形，這種處理易失真。第104頁/共173頁（2）凍結法根據(jù)水深判斷網(wǎng)格結點是否露出水面，對談度單元取糙率系數(shù)為一接近于無窮大的正數(shù)（糙率和水位均布置在網(wǎng)格中心），使單元四周的流速為一接近于0的無窮小量，這樣的處理結果相當于使該單元的潮位在計算過程中北凍結不變這種方法適用于寬淺，坡度較坦的露灘

54、問題，而對于潮灘相間的海岸、河口海域，則因水量和動量的過分凍結而失真。第105頁/共173頁（3）切削法稱水位判別法或薄層水法，該法對露灘單元并不凍結，而引入一個富裕水深（相當于灘地上存在很薄的水層）以保證計算過程的完整，相當于將原始地形切削降低，而一旦判斷實際水深大于富裕水深時，恢復原始地形，和開挖法具有相同的局限性。第106頁/共173頁（4）窄縫法假想在岸灘的每個網(wǎng)格上存在一條很窄的縫隙，它的深度和岸灘前的水深一致，根據(jù)水量平衡，將窄縫內的水量平鋪到岸灘上，把計算邊界設在岸灘的窄縫內，成為具有一定水深的固定邊界。系數(shù)不易確定，只適用于岸邊界的露灘問題。第107頁/共173頁bbzzssz

55、zBzzeBBBzBb0)(0)()( B(z) 為單位淺灘長度內的窄縫系數(shù)為單位淺灘長度內的窄縫系數(shù)第108頁/共173頁第109頁/共173頁2)()(2zBzBA第110頁/共173頁v 窄縫法的優(yōu)點是在計算中、不必隨時變化邊界計算點只礙事先定出一個足夠大的固定計算邊界，就可以動態(tài)地模擬出計算城內的水位、流速變化。由于計算域事先已固定，求解時的系數(shù)矩陣大小不變，有利于計算的穩(wěn)定性。又由于窄縫寬度很小。續(xù)內流速也可人為加以衰減，故對精度無大影響。缺點是窄縫參數(shù)不易確定。第111頁/共173頁（5）干濕法根據(jù)每步計算結果判斷每個單元干、濕？濕單元參加方程的計算。第112頁/共173頁二維淺水

56、水流的一種三角形網(wǎng)格二維淺水水流的一種三角形網(wǎng)格FVM計算格式計算格式對開闊寬淺型水域的水流運動, 可以守恒型二維非恒定淺水方程描述:ByGxFtQTbysybxsxTTTyZghfhuxZghfhvByvhghhhvvyuhhuvhvGxvhhuvxuhghhhuuhuFhvhuQ), 0()21,(),21,(),( 第113頁/共173頁式中: x , y 空間坐標; t 時間坐標; u, v 在x , y 方向沿水深積分平均流速分量; 潮位; h 水深; f 柯氏力系數(shù); 水的密度; Z 地形高程; x , y 方向沿水深平均的紊動粘性系數(shù); sx， sy 沿x , y 方向的風應力;

57、 bx， by沿x , y 方向的河床底應力.閉邊界( 岸邊界) : 采用流動法向通量為零開邊界( 水流邊界) 一般在河道較順直段選取已知流速或水位過程。第114頁/共173頁對計算區(qū)域中的任意一個三角形單元, 各水力變量P ( x , y ) 的計算值均布置在節(jié)點上, 并假定單元內變量分布呈一階近似:P( x , y ) = a + bx + cy第115頁/共173頁v基本方程的離散及求解對基本方程在控制體積( 圖1) 上積分:0)()(GdxFdydxdyBtQ 由于引入一階近似假定, 單元內各水力變量, h, u, v 均可表示為節(jié)點處變量值的線性函數(shù), 從而可直接求解積分方程. 離散

58、結果表示為如下形式:第116頁/共173頁183116bhStjjji()(buwuSCvugthjjjuijjiiiij 38311083222)()(bvwvSCvugthjjjvijjiiiij 第117頁/共173頁第118頁/共173頁)()(4)()()(4)(11112011111301jjjijjjijjjjjijjjijjjjjxxxxyyyyShhxxxxyyyyShhw第119頁/共173頁第120頁/共173頁第四章第四章 河口三維流體動力學模型河口三維流體動力學模型v在寬闊且較深的海岸河口地區(qū)，研究水流運動，海岸演變及泥沙運動時，通常二位數(shù)值

59、模擬就不能滿足要求了，此外，像疏浚拋泥、油膜運動、水質污染擴散等一些專門課題，不是二位數(shù)值模擬能解決的，必須采用三位數(shù)值模擬技術。第121頁/共173頁vPOM（Princeton Ocean Model)模型由Blumberg和Mellor1978提出，經(jīng)多年的改進，已成為比較廣泛使用的海洋模式。vPOM在淺水海域水深小于3米時，退潮時，模擬灘地退出水面遇到困難，計算不穩(wěn)定。vPOM模型采用模式分離技術，三維控制方程組及其定解條件構成模型的內模式，而外模由全積分內模式方程得到。第122頁/共173頁第四章第四章 河口三維流體動力學模型河口三維流體動力學模型v連續(xù)方程：0zwyvxuv動量方程

60、：xMFzuKzxpfvzuwyuvxuutu)(10yMFzvKzypfuzvwyvvxvutv)(10gzpMK為動量垂向渦粘系數(shù) 第123頁/共173頁v對于三維斜壓模型，還需要同時考慮溫度、鹽度的擴散過程，其控制方程為：THFzTKzzTwyTvxTutT)(SHFzSKzzSwySvxSutS)(),(pTS T為勢溫（對河口及近岸地區(qū)可為現(xiàn)場實際溫度），S為鹽度，HK為反映溫度、鹽度垂向紊動混合的垂向擴散系數(shù)。第124頁/共173頁xvyuAyxuAxFMMx2xvyuAxyvAyFMMy2ySTAyxSTAxFHHST, xFyFTFSF是對模型網(wǎng)格無法分辨的所謂次網(wǎng)格運動過程用