河口、海岸水動力模擬技術

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

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

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

4、定性和收斂性，并與現(xiàn)場資料有良好的吻合；（4）正式方案試驗 v河口、海岸水動力模擬的發(fā)展方向1、河口模型四維資料同化v四維同化在河口中的一類主要應用是為河口數(shù)值模式提供優(yōu)化的初始狀態(tài),從而校正不合理的邊界條件所帶來的偏差,提高模擬精度.v變分同化技術還可應用于確定河口模式中的未知參數(shù). 例如確定了潮汐河口的摩擦系數(shù),通常這種系數(shù)是通過經(jīng)驗獲得的.v另外,河口模式的外部強迫場(如風場、海-氣界面的熱通量等)都可以通過這種技術反演出來.資料同化將是河口數(shù)值模型發(fā)展和結合的一個新技術切入點,也是帶動河口動力數(shù)字模擬技術革新的一種重要方法.2、數(shù)字河口動力模型v河口水動力特征及外界強迫作用因子構建數(shù)字

5、河口動力模型(模型計算要素象溫、鹽、流、泥沙的特征均以數(shù)字矩陣形式記載),若再與數(shù)字河網(wǎng)模型嵌套聯(lián)結,最終可以獲得區(qū)域河網(wǎng)河口的動力-沉積-地貌的機制解釋,揭示河口物理規(guī)律和解決工程實際問題.v在數(shù)字河口模型的構架下,可以將人類已經(jīng)擁有的河口科學理論、知識與具有較強物理概念的水動力學模型集成于一體,為河口數(shù)值模型計算提供良好的平臺。數(shù)字河口動力模型具有許多優(yōu)勢:首先,數(shù)字河口模型是基于數(shù)字區(qū)域地形構建而成的,地形要素可自動生成,無需手工操作,大大提高了工作效率;其次,數(shù)字模型不僅能輸出傳統(tǒng)模型的結果,而且能夠十分方便地給出河口水文要素和水文狀態(tài)變量的空間分布場,這些對近岸河口動力科學研究與河口

6、、港口、航道工程都有著廣闊的應用前景.總而言之,數(shù)字河口模型研究的最終目的就是利用已有的河口基礎科學理論和知識,在數(shù)字區(qū)域地形的基礎之上將觀測點的水文信息拓展、同化至區(qū)域平面上乃至區(qū)域三維立體上的信息,并形成數(shù)字成品。v參考文獻：vKoutitar 著“Mathematical Model in Coastal Engineering”1)模型簡單易懂2）附有Basic程序，而且有驗證的算例3）介紹各種數(shù)值處理技術v曹祖德、王運洪”水動力泥沙數(shù)值模擬第二章第二章 水動力數(shù)值模擬的理論基礎水動力數(shù)值模擬的理論基礎)()()()(022222222zVAzyVxVAygfUdtDVzUAzyUxU

7、AxgfVdtDUzWyVxUzHzH ),()tyxzyVxUtW2.1 基本方程自由面運動學邊界條件：底部運動學邊界條件：),()(yxhzyhVxhUWzWyVxUtdtDU,V,W為x,y,z 方向上的流速分量。（，）為距平均海平面的自由表面水位。（，）為平均海平面距底部邊界的水深。為水平擴散系數(shù)。為垂直渦動系數(shù)。0)()()()(222222yQxQtQDCVUgygDfQQVyQUxtQQDCVUgxgDfQQVyQUxtQyxyxyyyxyxxx初始條件),()0 ,(0yxuyxu),()0 ,(0yxvyxv),()0 ,(0yxyx邊界條件岸邊界：法向流速為零。水邊界：給定

8、潮位過程。0)()()(0oxbxbbBxAZhgxZhgAxuQtQxQtASaint Venant 方程三、二、一維方程的定解條件三、二、一維方程的定解條件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。0nV2.22.2數(shù)值計算數(shù)值計算v在計算水動力、泥沙數(shù)值模擬時，大都將基本方程組離散成代數(shù)方程組，最后求解代數(shù)方程組，此處介紹微分方程組的離散技術有限差分法和線性代數(shù)方程組的數(shù)值解法。.

9、1有限差分法有限差分法v有限差分法是工程中常用的一種離散技術，將計算域分成有限個網(wǎng)格，通過差分法求網(wǎng)格結點的微分方程的近似值，也稱網(wǎng)格法。v將網(wǎng)格結點上的函數(shù)f(x,y,z,t)表示成 ，i,j,k分別表示x,y,z方向的坐標位置，n表示時間。nkjif,1、工程中常用的幾種差分和微分的關系（一維）（1）一階向前差分10 ,)(2)()()(1122xxxxxfxxxfxxfxxf10 ,)(2)()()(2222xxxxxfxxxxfxfxxf（2）一階向后差分（3）一階中心差分10 , 2/10 , 2/)()(48)2()2()(4443333433332xxxxxxxxfxxfxxxx

10、fxxfxxf（4）二階中心差分10 , 2/10 , 2/)()(96)2/()(2)2()2()(6665554644542222xxxxxxxxfxxfxxxfxxfxxfxxf2、幾種常見的差分格式以一維熱傳導方程為例：022xftf（1）古典顯式格式0)(222111xtOxffftffninininini022111xffftffninininini2111/),()21 (xtrffrfrfnininini（2）古典隱式格式0)(222111111xtOxffftffninininini022111111xffftffninininini21111/,)21 (xtrfrffrr

11、fnininini（3）六點格式（Crank-Nicolson），雙層六點隱式格式在x點和n+n/2時層，對t和x均采用中心差分0)(2221222112111111xtOxfffxffftffnininininininini21111111/,22)1 (2/)1 (2/xtrfrfrfrfrfrfrnininininini022212112111111xfffxffftffnininininininini（4）Richardson格式，三層顯式格式在x點和n時層，對t和x均采用中心差分0)(222221111xtOxffftffninininini21111/),2(2xtrfffrffn

12、inininini02221111xffftffninininini（5）加權六點格式，隱式格式在x點和n+n時層，01，對t和x均采用中心差分2/10)(2)1 (222112111111xtOxfffxffftffnininininininini21111111/)1 ()1 ()1 (21 )21 (xtrrfrffrrffrrfnininininini02)1 (22112111111xfffxffftffnininininininini.2線性方程組的數(shù)值解線性方程組的數(shù)值解v有限差分法是工程中常用的一種離散技術，將計算域分成有限個網(wǎng)格，通過差分法求網(wǎng)格結點的微分方程

13、的近似值，也稱網(wǎng)格法。v將網(wǎng)格結點上的函數(shù)f(x,y,z,t)表示成 ，i,j,k分別表示x,y,z方向的坐標位置，n表示時間。nkjif,1、解線性方程組的兩種方法:v直接法:通過有限步算術運算直接求出方程組的精確解,最常用的是消元結合代入的方法.v實際上除非是采用無窮位精度計算,一般都得不到精確解.v直接方法適用于解低階稠密矩陣方程組.v迭代法 類似于方程求根的迭代法,用一個迭代過程逐步逼近方程組的解.v迭代有可能不收斂,或雖然收斂,但收斂速度慢.v迭代法適用于求解高階稀疏矩陣方程組.v稀疏矩陣:矩陣非零元素較少,且在固定的位置上.v稀疏矩陣一般是人為構造的,例如36頁三轉角插值時方程組(

14、8.12),(8.15)的系數(shù)矩陣.GaussGauss消去法消去法( (第一次消元第一次消元) )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.GaussGauss消去法消去法( (第一次消元第一次消元) )v經(jīng)

15、過上述過程,得到方程組A(2)x=b(2),1111111211122222222222200( )( )( )( )( )( )( )( )( )( )nnnnnnnaaabxxaabxbaa v其中2111 1( )( )( ),ijijijaam a2111 1( )( )( ),iiibbm b111111( )( )iiama 2 3( , , )i jn GaussGauss消去法消去法( (第第k k次消元次消元) )v假設已完成k-1次消元,得到方程組A(k)x=b(k).111111112111122222221221111111000000000( )( )( )( )(

16、),( )( )( )( ),()()(),( )( )( )( )kknkknkkkkkkkknkkkkknkknknnaaaaaaaaaaaaaaaa v第k次消元的目的是將akk(k) (稱為主元)下面的元素變?yōu)?.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).GaussGauss消去法消去法( (

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

18、abxxaabxba v一共需進行 ? 步n-1GaussGauss消去法消去法( (算法算法) )1,1kn( )( )/kkikikkkmaa 1,ikn1()( )( )kkkijijikkjaam a 1,i jkn1()( )( )kkkiiikkbbm b 1,ikn( )( )/nnnnnnxba 1,1kn( )( )( )1()/nkkkkkkjjkkj kxbaxa 追趕法求解三對角方程組追趕法求解三對角方程組v上面的方程組可以利用追趕法求解(P185).v對于下面形式的方程組11112222211111nnnnnnnnnfbcxabcxfabcxfabxf v將系數(shù)矩陣進

19、行三角分解v比較兩邊對應元素可以得到11223111111nnnnrrr 11222111nnnnnbcabcabcab11,b 111,c 2(, ),iiar in12(, ),iiiibrin 21(,).iiicin v因此有11,b 111,c 2(, ),iiar in12(, ),iiiibrin 21(,).iiicin iiic 1iiiicbr 121(,)iiiicinba v又11111ccb v因此所有i的可遞推求出,進一步可求出i,ri.v在得到系數(shù)矩陣的分解后,原方程組轉化為vLUx=f.v先求解Ly=f11122223111nnnnnnnyfryfryfryf

20、v顯然有y1=f1/1,vyi=(fi-riyi-1)/i=(fi-iyi-1)/(bi-aii-1)(i=2,n)v再求解Ux=y,111221111111nnnnnyxxyxyxy v顯然有xn=yn, xi=yi-bixi+1(i=n-1,1)迭代法迭代法v在處理一元方程f(x)=0時,我們將其轉化為x=j(x)的形式,然后用不動點迭代的方法進行求解.v對于線性方程組Ax=b,我們也可以將其轉化為類似的形式: x=Bx+f,v任取初始向量x(0),令x(k+1)=Bx(k)+f(k=0,1,),則得到一個向量的序列x(k).v若該序列收斂于向量x*,對x(k+1)=Bx(k)+f 兩邊取

21、極限得到x*=Bx*+f,即x*是方程組的解.JacobiJacobi迭代法與迭代法與Gauss-SeidelGauss-Seidel迭代法迭代法v對于方程組1231231238322041133631236xxxxxxxxx v我們將其改寫為1232133121322081433111633612()()()xxxxxxxxx JacobiJacobi迭代法迭代法v寫成矩陣的形式為x=B0 x+f,其中03208841011116301212B 20833113612f JacobiJacobi迭代法迭代法v利用x(k+1)=Bx(k)+f 進行迭代,得到結果如下kx1(k)x2(k)x3(

22、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 JacobiJacobi迭代法迭代法v從上表可以看出,迭代序列逐步逼近方程組的精確解(3,2,1)T.v注:在迭代中,我們不可能得到x(k)和精確解之間的誤差,一般我們用|x(k)-x(k-1)|(通常用無窮范數(shù))的值來判

23、斷是否終止迭代.v在上面的例子中,我們將第i個方程變形為左邊是xi,右邊是其它分量和常數(shù)的線性組合,然后進行迭代,這一方法稱位Jacobi迭代.JacobiJacobi迭代法迭代法v一般的,對于方程組Ax=b,設A非奇異且aii0(i=1,2,n),將A改寫為A=D L U,其中112233nnaaDaa 2131321230000nnnaLaaaaa 1213123230000nnnaaaaaUa JacobiJacobi迭代法迭代法v將方程組改寫為vDx=(L+U)x+bvx=D1(L+U)x+D1bv令B0=D1(L+U)(稱位Jacobi迭代矩陣),f=D1b,上式簡記為x=B0 x+

24、f.v我們得到Jacobi迭代公式 x(k+1)=B0 x(k)+f.111()( )()nkkiiijjiijj ixba xa v寫成分量的形式為Gauss-SeidelGauss-Seidel迭代法迭代法v在前面的例子中,我們計算x1(k+1),用的是第k步的x2,x3;1123121313121322081433111633612()( )( )()( )( )()( )( )()()()kkkkkkkkkxxxxxxxxx v計算x2(k+1),用的是第k步的x1,x3,我們有理由認為已經(jīng)計算出的第k+1步的x1比第k步的“好”.因此,我們應該用第k+1步的x1和第k步的x3來計算x

25、2.11()kx v類似地,我們也應該用新信息計算x3.11()kx 12()kx Gauss-SeidelGauss-Seidel迭代法迭代法v我們可以將上面一般的Jacobi迭代公式改寫為111()( )()nkkiiijjiijj ixba xa 111111()()( )()inkkkiiijjijjiijj ixba xa xa v這一迭代方法稱為Gauss-Seidel迭代.Gauss-SeidelGauss-Seidel迭代法迭代法( (算例算例) )v其Gauss-Seidel迭代公式為v對于方程組1231231238322041133631236xxxxxxxxx 11231

26、12131113121322081433111633612()( )( )()()( )()()()()()()kkkkkkkkkxxxxxxxxx Gauss-SeidelGauss-Seidel迭代法迭代法( (算例算例) )v同樣取x(0)=(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.00

27、0163029.98e-352.99984239 2.00007213 1.000060773.84e-41( )()|kkxx 超松弛迭代超松弛迭代(SOR)(SOR)方法方法v沿著從xi(k)到xi (k+1) (G)的方向再向前走,就得到超松弛迭代(SOR)方法.( )kix1()()ki Gx 11( )()()()kkii Gxx v假設已知第k步的迭代向量x(k)以及第k+1步迭代向量x(k+1)的前i1個分量已知,Gauss-Seidel迭代法取111111()()( )()inkkkiiijjijjiijj ixba xa xa 超松弛迭代方法超松弛迭代方法v我們定義新的xi(

28、k+1)為xi(k)與 的加權平均.1()kix 1111()( )()( )()( )()()kkkkkkiiiiiixxxxxx111( )()( )()inkkkiiijjijjiijj ixba xa xa v在=1時,上述方法就是Gauss-Seidel方法,1時稱為超松弛法(有時不管的范圍,統(tǒng)稱為超松弛方法).超松弛迭代方法超松弛迭代方法( (算例算例) )v對于方程組123441111141111141111141xxxx v松弛方法迭代格式為(1)( )( )( )( )( )111234(1)( )(1)( )( )( )221234(1)( )(1)(1)( )( )331

29、234(1)( )(1)(1)(1)( )441234(14)/4(14)/4(14)/4(14)/4kkkkkkkkkkkkkkkkkkkkkkkkxxxxxxxxxxxxxxxxxxxxxxxx 超松弛迭代方法超松弛迭代方法( (算例算例) )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

30、.99999179-1.00000289-1.000001703.45e-511-0.99999667-1.00000287-0.99999954-0.999999191.11e-512-1.00000152-0.99999922-1.00000012-1.000000524.85e-61( )()|kkxx 超松弛迭代方法超松弛迭代方法( (算例算例) )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

31、 2.0步數(shù)1712121518243555114*v步數(shù)表示|x(k)x(k1)|10-5時的迭代步數(shù),=2.0時,500步以內不收斂.超松弛迭代方法超松弛迭代方法( (矩陣表示矩陣表示) )v超松弛迭代格式可以寫為1111()( )()( )()inkkkkiiiijjijjiijj ixxba xa xa 111111()( )()( )()()inkkkkiiiiiiiijjijjjj ia xa xba xa x v用矩陣可以表示為111()( )()( )()()kkkkDxDxbLxUx11()( )()()kkDL xDU xb 第三章第三章 二維水動力數(shù)值模擬二維水動力數(shù)值模

32、擬一、二維水動力數(shù)值模擬系統(tǒng)的分類1、按差分網(wǎng)格分：三角形、正方形、矩形、四邊形、多邊形、曲線坐標網(wǎng)格以及各種形狀網(wǎng)格的組合2、按計算方法分：顯式法、隱式法、顯隱混合法3、按模擬格式分：三角元法、ADI法、破開算子法、單元體積法、MADI法、準分析法、貼體坐標法。二、平面二維水動力數(shù)學模型的一般形式000)()(00hhfuygyvvxvutvhhfvxgyuvxuutuyvhxuhtsybysxbxwwwassywwwassxbybxvvufuvufhvvunhuvun22223/12223/1222v定解條件；水邊界：給定潮位過程；岸邊界：法向流速為0),()0 ,(0)0 ,(0)0 ,

33、(0yxyxyxvyxu三、三、ADIADI法法1、網(wǎng)格正方形或矩形，變量u,v,分別交錯布置于網(wǎng)格中心和兩側。2、ADI基本思想（1）分步（2）交錯顯隱ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位 、水深、水深uv3、差分格式X向運動方程在（i+1/2,j)點離散)(21)(21)(41)2(2;2;)()(2)(41;2, 1. 100002/1, 12/1, 12/1,2/1, 2/1, 2/11, 2/11, 2/1, 2/1, 2/1022, 2/12, 2/1, 2/1, 2/32/1, 12/1, 2/12/1, 1,2/1njinjinjinnj

34、injinjinjijijinjinjinjiiinjinjinjinjinjiiiinjiinjiinjiijijijihhhvvvvvvfyuutudxtgchCvutguuxtbxtgadcuba）其中連續(xù)方程在（i,j)點離散)()(2;)(2; 1;)(22/1,02/1,2/1,02/1, 2/10, 2/102/1, 2/12/1,2/1, 2/1njinjinjinjinjiinjiiinjiiinjiinjiinjiihvhvytdhxtcbhxtaducbua其中聯(lián)立上面的式子得到下面的線性方程組injiinjiinjiidcuba2/1, 12/1, 2/12/1, in

35、jiinjiinjiiducbua2/1, 2/12/1,2/1, 2/1其中，u與存在如下關系iiiiiiiiiiiiiiiiiiiiiiiiiniiniiniinibEaFadHbEacGbGaHadFbGacEFuEHGu2/12/12/12/12/12/12/12/12/12/12/12/12/1;其中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,jijijijijin

36、jinjinjinjinjinjinjiinjinjinjinjinjiiinjiiuuuuuu ftytgvvxtuvMhCvugtvvytNMvN其中v在后半個時間步長內，按上述同樣原理，在y向掃描。因此只須將上述各計算公式做如下的變換：xy,xy;ij,ij; uv,uv;即可求出vi,j+1/2n+1, i,jn+1,然后顯式求ui+1/2,jn+1。四、分步全隱式格式將x向動量方程與連續(xù)方程聯(lián)立，求出un+1,n+1/2,將y向動量方程與連續(xù)方程聯(lián)立，求出vn+1,n+1,ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位水深水深uvX向動量方程在（i+1/

37、2,j)點離散0)()()()()()()()()(, 2/1, 2/1, 2/1, 2/12, 2/12, 2/12, 2/12/12, 2/12, 2/12/1, 2/12/1,2/1, 12/1, 2/12/1, 2/1, 2/1, 11, 2/1, 2/11, 2/1njinjinjisxnjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihv fhCvuugxgyuuvxuuutuu 連續(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,

38、 2/1,2/1,yhvhvxhuhutnjinjinjinjinjinjinjinjinjinjinjinjinjinji )(21)(1)(41)(21)(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, 1njinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinji

39、njihnCvvvvvuuuuuuuuuuuuuuu injiinjiinjiidcuba2/1, 11, 2/12/1, injiinjiinjiiducbua1, 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 五、移步雙向交替顯、隱式交錯法MADI法：將水深、水位、流速等變化量均布置在同一網(wǎng)格節(jié)點上，由此將基本方程離散成新的差分代數(shù)方程組，并建立一種新的解法，這種解法既吸收了原有

40、傳統(tǒng)的ADI法的優(yōu)點，又有較高的穩(wěn)定性、收斂性和精度。方程離散時，時間導數(shù)項采用前差表示，空間采用中心差分，將t分成兩等分。v當nt(n+1/2)t,x向動量方程和連續(xù)方程建立差分方程iniiniiniiininiiniidcubadcbua2/112/12/112/112/12/112/112/11211111,00niniiiiiiiiiiiiiudduucbacba2/12/111211111,00niniiiiiiiiiiiiiudduucbacba)()(2)(21)(2)2(22)()()(2)(4121,01,1,01, 10, 10,1,1,02,2,2, 1, 1njinji

41、njinjinjiinjiiinjiinjinjinjinjiiinjijinjinjinjinjiiihvhvytdhxtcbhxtavfyuutudxtgchCvutguuxtbxtga)()()()()()(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(

42、n+1)t,y向動量方程和連續(xù)方程建立差分方程,求解v,顯式求u，過程同上。六、三角元法由于采用矩形網(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)格疏些；000)()(hfuygyvvxvutvhfvxgyu

43、vxuutuyhvxhutbybx 基本方程離散后方程0)()()()(0)()()()(0)()(11111111111nibynininininininininibxninininininininininininihfuygyvvxvutvvhfvxgyuvxuutuuyhvxhut iiiyaxaaF210yaxaayxF210),(jjjyaxaaF210kkkyaxaaF210ekkkkjjjjiiiisFycxbaFycxbaFycxbayxF2/ )()()(),(ekkkkejjjjeiiiisycxbaNsycxbaNsycxbaN2/ )(2/ )(2/ )(kkjjiiF

44、NFNFNyxF),(kkjjiikkjjiiFcFcFcyFFbFbFbxFekkjjiiSFbFbFbMyF2/ )()(ekkjjiiSFcFcFcMxF2/ )()(eeSS出于該格式采用顯式求解。同時受到復雜地形的影響，在計算過程中會產(chǎn)生一些擾。當這些擾動擴大并傳播易使計算失敗為消除這些擾動的影響，采用濾波公式七、邊界處理邊界處理合適與否影響數(shù)值模擬的成敗。實際工程中，計算域的邊界常由不規(guī)則的曲線組成，如何利用有限差分的矩形網(wǎng)格來模擬曲線邊界，以及如何選取邊界值，這是邊界處理的重要內容。1、邊界類型（1）Dirichlet邊界（2）Neumann邊界（3）淺灘活動邊界2、 Diric

45、hlet邊界在邊界處有已知的函數(shù)值。在實際情況下，水邊界通常屬于這類邊界，該處的實測水位、流速，可作為已知函數(shù)值。1ffv如果邊界網(wǎng)格結點正好在實測點上，則該結點的邊界值可直接采用已知值。v如果邊界網(wǎng)格結點不在實測點上，則分以下不同情況分別處理。v邊界結點與實測點很靠近，邊界值直接等于實測值；v如果實測點與邊界點較遠，但與虛擬外結點較近，則令虛擬外結點直接等于實測值，然后引入邊界條件求出邊界值。v邊界網(wǎng)格結點或虛擬結點與實測點均不靠近，不能直接引用實測值，這時，可根據(jù)前移時間步長所得的結果，用線性差值可得邊界值。xiicxjijixiicjijijixxfffxxffff )(, 1,)(,

46、1, 13、Neumann邊界在實際情況下，固定邊界屬于此類邊界2)(ffnfSv邊界正好通過網(wǎng)格或邊界與網(wǎng)格結點很靠近v邊界的外法線方向與坐標軸平行，直接得到邊界值；v邊界的外法線方向與坐標軸不平行，考慮外法線與坐標軸的夾角，帶入邊界條件后離散得到邊界值。v邊界與網(wǎng)格結點距離較大4、淺灘活動邊界（1）開挖法：將灘地開挖至可能出現(xiàn)的最低水位之下，為使水量平衡，將岸邊界向水域內移動，并增大開挖部位的糙率以求得動量上的平衡。這種方法一般只適用于潮灘問題，而且灘地面積只占整個海區(qū)較小的情況，而對于水位變化小，坡度較緩的地形，這種處理易失真。（2）凍結法根據(jù)水深判斷網(wǎng)格結點是否露出水面，對談度單元取糙

47、率系數(shù)為一接近于無窮大的正數(shù)（糙率和水位均布置在網(wǎng)格中心），使單元四周的流速為一接近于0的無窮小量，這樣的處理結果相當于使該單元的潮位在計算過程中北凍結不變這種方法適用于寬淺，坡度較坦的露灘問題，而對于潮灘相間的海岸、河口海域，則因水量和動量的過分凍結而失真。（3）切削法稱水位判別法或薄層水法，該法對露灘單元并不凍結，而引入一個富裕水深（相當于灘地上存在很薄的水層）以保證計算過程的完整，相當于將原始地形切削降低，而一旦判斷實際水深大于富裕水深時，恢復原始地形，和開挖法具有相同的局限性。（4）窄縫法假想在岸灘的每個網(wǎng)格上存在一條很窄的縫隙，它的深度和岸灘前的水深一致，根據(jù)水量平衡，將窄縫內的水量

48、平鋪到岸灘上，把計算邊界設在岸灘的窄縫內，成為具有一定水深的固定邊界。系數(shù)不易確定，只適用于岸邊界的露灘問題。bbzzsszzBzzeBBBzBb0)(0)()( B(z) 為單位淺灘長度內的窄縫系數(shù)為單位淺灘長度內的窄縫系數(shù)2)()(2zBzBAv 窄縫法的優(yōu)點是在計算中、不必隨時變化邊界計算點只礙事先定出一個足夠大的固定計算邊界，就可以動態(tài)地模擬出計算城內的水位、流速變化。由于計算域事先已固定，求解時的系數(shù)矩陣大小不變，有利于計算的穩(wěn)定性。又由于窄縫寬度很小。續(xù)內流速也可人為加以衰減，故對精度無大影響。缺點是窄縫參數(shù)不易確定。（5）干濕法根據(jù)每步計算結果判斷每個單元干、濕？濕單元參加方程的

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

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

51、8311083222)()(buwuSCvugthjjjuijjiiiij 38311083222)()(bvwvSCvugthjjjvijjiiiij )()(4)()()(4)(11112011111301jjjijjjijjjjjijjjijjjjjxxxxyyyyShhxxxxyyyyShhw第四章第四章 河口三維流體動力學模型河口三維流體動力學模型v在寬闊且較深的海岸河口地區(qū)，研究水流運動，海岸演變及泥沙運動時，通常二位數(shù)值模擬就不能滿足要求了，此外，像疏浚拋泥、油膜運動、水質污染擴散等一些專門課題，不是二位數(shù)值模擬能解決的，必須采用三位數(shù)值模擬技術。vPOM（Princeton O

52、cean Model)模型由Blumberg和Mellor1978提出，經(jīng)多年的改進，已成為比較廣泛使用的海洋模式。vPOM在淺水海域水深小于3米時，退潮時，模擬灘地退出水面遇到困難，計算不穩(wěn)定。vPOM模型采用模式分離技術，三維控制方程組及其定解條件構成模型的內模式，而外模由全積分內模式方程得到。第四章第四章 河口三維流體動力學模型河口三維流體動力學模型v連續(xù)方程：0zwyvxuv動量方程：xMFzuKzxpfvzuwyuvxuutu)(10yMFzvKzypfuzvwyvvxvutv)(10gzpMK為動量垂向渦粘系數(shù) v對于三維斜壓模型，還需要同時考慮溫度、鹽度的擴散過程，其控制方程為：

53、THFzTKzzTwyTvxTutT)(SHFzSKzzSwySvxSutS)(),(pTS T為勢溫（對河口及近岸地區(qū)可為現(xiàn)場實際溫度），S為鹽度，HK為反映溫度、鹽度垂向紊動混合的垂向擴散系數(shù)。xvyuAyxuAxFMMx2xvyuAxyvAyFMMy2ySTAyxSTAxFHHST, xFyFTFSF是對模型網(wǎng)格無法分辨的所謂次網(wǎng)格運動過程用水平紊動擴散過程參數(shù)化后的產(chǎn)生項，分別為：對基本方程簡化時采用的假定與近似對基本方程簡化時采用的假定與近似 （1）靜壓假定：在河口、近岸淺水地區(qū)，垂向速度的時間變率（即：垂向加速度）與重力加速度相比甚微，可略去不計，因此垂向動量方程可簡化為：，即壓強

54、沿水深的變化符合靜水壓強分布。（2）Boussinesq近似：海水密度為時均值（參考密度）和脈動值之和，將其代入動量方程后，除在重力加速度的前面保留外，其余各項的均略去。（3）Boussinesq假定：由于在時均運動方程中包含了較難處理的雷諾應力張量，Boussinesq在1877年提出了關于可以將水流紊動應力類比于層流粘性應力的假定，即用層流粘性應力的形式對紊動應力進行參數(shù)化。v坐標系中的動力學方程： Hz010tyDVxDUxMFUDKdxDDxgDxgDfVDUyUVDxDUtUD0022yMFVDKdyDDygDygDfUDVyDVxUVDtVD0022ttDyyDVxxDUWv坐標系

55、中的溫度、鹽度擴散方程分別為 THFTDKTyTVDxTUDtTD)(SHFSDKSySVDxSUDtSD)(v紊流閉合模型：得到動量、標量垂向擴散系數(shù)KM、KH qHMqFlBDqKgVUDKqDKqyDVqxDUqtDq1302222222222lHMqFWBDqKgEVUDKlElqDKlqylDVqxlDUqtlDq1303221222220yDVxDUtddxDxDgDGxgDDVfFyDVUxDUtDUxbxsxx01002ddyDyDgDGygDDUfFyDVxDVUtDVybysyy01002v外模式方程外模式方程v邊界條件v1、自由表面邊界條件sysxMVUDK,0SHSTD

56、KH,023212sUBq02lq0 yxz,0 v2 2、水體底部邊界條件、水體底部邊界條件yxHz,1bybxMVUDK,023212bUBq0,ST02lq0 v3 3、側向閉邊界條件、側向閉邊界條件0nU0,STn v4、側向開邊界條件v對水位開邊界，通常用實測的水位資料或者用更大范圍數(shù)學模型計算的水位值作為強迫水位控制條件。v對于溫、鹽等守恒性標量物質的開邊界條件，通常分為入流邊界和出流邊界兩種情況。對于入流邊界，一般采用開邊界實測的溫、鹽數(shù)據(jù)；對于出流邊界，采用對流型開邊界形式：v數(shù)值求解方法（1）模式分裂技術 將垂向積分的運動方程（外模態(tài)）從反映流速垂向結構的運動方程（內模態(tài)）中分離出來，用較少的計算量通過求解外模態(tài)得到自由表面，然后通過求解內模態(tài)得到水流的垂向結構，這稱為模式分裂