河口、海岸水動(dòng)力模擬技術(shù)_2014_12

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

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

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

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

5、值觀測(cè)資料的機(jī)制方面有高效的預(yù)測(cè)能力。模型應(yīng)對(duì)前期信息在時(shí)間上能向前插值。v動(dòng)力場(chǎng)的模型同化熟路也被證明對(duì)陸-氣系統(tǒng)內(nèi)部環(huán)流的科學(xué)研究極有價(jià)值v同化是唯一能產(chǎn)生復(fù)雜非線性過長(zhǎng)的場(chǎng)的手段，并且在物理學(xué)與動(dòng)力學(xué)的機(jī)制上協(xié)調(diào)一致。2、數(shù)字河口動(dòng)力模型v河口水動(dòng)力特征及外界強(qiáng)迫作用因子構(gòu)建數(shù)字河口動(dòng)力模型(模型計(jì)算要素象溫、鹽、流、泥沙的特征均以數(shù)字矩陣形式記載),若再與數(shù)字河網(wǎng)模型嵌套聯(lián)結(jié),最終可以獲得區(qū)域河網(wǎng)河口的動(dòng)力-沉積-地貌的機(jī)制解釋,揭示河口物理規(guī)律和解決工程實(shí)際問題.v在數(shù)字河口模型的構(gòu)架下,可以將人類已經(jīng)擁有的河口科學(xué)理論、知識(shí)與具有較強(qiáng)物理概念的水動(dòng)力學(xué)模型集成于一體,為河口數(shù)值模型

6、計(jì)算提供良好的平臺(tái)。數(shù)字河口動(dòng)力模型具有許多優(yōu)勢(shì):首先,數(shù)字河口模型是基于數(shù)字區(qū)域地形構(gòu)建而成的,地形要素可自動(dòng)生成,無需手工操作,大大提高了工作效率;其次,數(shù)字模型不僅能輸出傳統(tǒng)模型的結(jié)果,而且能夠十分方便地給出河口水文要素和水文狀態(tài)變量的空間分布場(chǎng),這些對(duì)近岸河口動(dòng)力科學(xué)研究與河口、港口、航道工程都有著廣闊的應(yīng)用前景.數(shù)字河口模型研究的最終目的就是利用已有的河口基礎(chǔ)科學(xué)理論和知識(shí),在數(shù)字區(qū)域地形的基礎(chǔ)之上將觀測(cè)點(diǎn)的水文信息拓展、同化至區(qū)域平面上乃至區(qū)域三維立體上的信息,并形成數(shù)字成品。v參考文獻(xiàn)：vKoutitar 著“Mathematical Model in Coastal Engin

7、eering”1)模型簡(jiǎn)單易懂2）附有Basic程序，而且有驗(yàn)證的算例3）介紹各種數(shù)值處理技術(shù)v曹祖德、王運(yùn)洪”水動(dòng)力泥沙數(shù)值模擬第二章第二章 水動(dòng)力數(shù)值模擬的理論基礎(chǔ)水動(dòng)力數(shù)值模擬的理論基礎(chǔ))()()()(022222222zVAzyVxVAygfUdtDVzUAzyUxUAxgfVdtDUzWyVxUzHzH ),()tyxzyVxUtW2.1 基本方程自由面運(yùn)動(dòng)學(xué)邊界條件：底部運(yùn)動(dòng)學(xué)邊界條件：),()(yxhzyhVxhUWzWyVxUtdtDU,V,W為x,y,z 方向上的流速分量。（，）為距平均海平面的自由表面水位。（，）為平均海平面距底部邊界的水深。為水平擴(kuò)散系數(shù)。為垂直渦動(dòng)系數(shù)。

8、0)()()()(222222yQxQtQDCVUgygDfQQVyQUxtQQDCVUgxgDfQQVyQUxtQyxyxyyyxyxxx初始條件),()0 ,(0yxuyxu),()0 ,(0yxvyxv),()0 ,(0yxyx邊界條件岸邊界：法向流速為零。水邊界：給定潮位過程。0)()()(0oxbxbbBxAZhgxZhgAxuQtQxQtASaint Venant 方程A 過水?dāng)嗝婷娣e，Q為流量，Q=Bhu，B為河寬，h為水深，u為斷面平均流速。 分別為水面剪切力與水底剪切力Zb 為水底的豎向坐標(biāo)位置X軸為沿水流縱向方向。 三、二、一維方程的定解條件三、二、一維方程的定解條件v初始

9、條件u,v,w,|t=0=u0,v0,w0,0邊界條件開邊界：計(jì)算域水體與外部水體相接處。(u,v,w)=(u(t),v(t),w(t)=(t)固邊界：計(jì)算域與陸地或建筑物接壤處無滑動(dòng)：u,v,w=0有滑動(dòng)： 垂直邊界的速度為0。0nV2.22.2數(shù)值計(jì)算數(shù)值計(jì)算v在計(jì)算水動(dòng)力、泥沙數(shù)值模擬時(shí)，大都將基本方程組離散成代數(shù)方程組，最后求解代數(shù)方程組，此處介紹微分方程組的離散技術(shù)有限差分法和線性代數(shù)方程組的數(shù)值解法。.1有限差分法有限差分法v有限差分法是工程中常用的一種離散技術(shù)，將計(jì)算域分成有限個(gè)網(wǎng)格，通過差分法求網(wǎng)格結(jié)點(diǎn)的微分方程的近似值，也稱網(wǎng)格法。v將網(wǎng)格結(jié)點(diǎn)上的函數(shù)f(x,

10、y,z,t)表示成 ，i,j,k分別表示x,y,z方向的坐標(biāo)位置，n表示時(shí)間。nkjif,1、工程中常用的幾種差分和微分的關(guān)系（一維）（1）一階向前差分10 ,)(2)()()(1122xxxxxfxxxfxxfxxf10 ,)(2)()()(2222xxxxxfxxxxfxfxxf（2）一階向后差分（3）一階中心差分10 , 2/10 , 2/)()(48)2()2()(4443333433332xxxxxxxxfxxfxxxxfxxfxxf（4）二階中心差分10 , 2/10 , 2/)()(96)2/()(2)2()2()(6665554644542222xxxxxxxxfxxfxxxf

11、xxfxxfxxf2、幾種常見的差分格式以一維熱傳導(dǎo)方程為例：022xftfxti-1 i i+1n+1 n（1）古典顯式格式0)(222111xtOxffftffninininini022111xffftffninininini2111/),()21 (xtrffrfrfnininini可由已知值直接求解未知值，且只涉及二層，稱為雙層顯式差分格式（2）古典隱式格式0)(222111111xtOxffftffninininini022111111xffftffninininini21111/,)21 (xtrfrffrrfnininini不能直接求解未知值，且只涉及二層，稱為雙層隱式差分格式（

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

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

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

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

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

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

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

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

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

21、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時(shí),我們將其轉(zhuǎn)化為x=j(x)的形式,然后用不動(dòng)點(diǎn)迭代的方法進(jìn)行求解.v對(duì)于線性方程組Ax=b,我們也可以將其轉(zhuǎn)化為類似的形式: x=Bx+f,v任取初始向量x(0),令x(k+1)=Bx(k)+f(k=0,1,),則得到一個(gè)向量的序列x(k).v若該序列收斂于向量x*,對(duì)x(k+1)=Bx(k)+f 兩邊取

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

23、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ù))的值來判

24、斷是否終止迭代.v在上面的例子中,我們將第i個(gè)方程變形為左邊是xi,右邊是其它分量和常數(shù)的線性組合,然后進(jìn)行迭代,這一方法稱位Jacobi迭代.JacobiJacobi迭代法迭代法v一般的,對(duì)于方程組Ax=b,設(shè)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,上式簡(jiǎn)記為x=B0 x+

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

26、2.11()kx v類似地,我們也應(yīng)該用新信息計(jì)算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對(duì)于方程組1231231238322041133631236xxxxxxxxx 11231

27、12131113121322081433111633612()( )( )()()( )()()()()()()kkkkkkkkkxxxxxxxxx Gauss-SeidelGauss-Seidel迭代法迭代法( (算例算例) )v同樣取x(0)=(0,0,0)T,迭代結(jié)果如下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

28、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假設(shè)已知第k步的迭代向量x(k)以及第k+1步迭代向量x(k+1)的前i1個(gè)分量已知,Gauss-Seidel迭代法取111111()()( )()inkkkiiijjijjiijj ixba xa xa 超松弛迭代方法超松弛迭代方法v我們定義新的xi(

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

30、1234(1)( )(1)(1)(1)( )441234(14)/4(14)/4(14)/4(14)/4kkkkkkkkkkkkkkkkkkkkkkkkxxxxxxxxxxxxxxxxxxxxxxxx 超松弛迭代方法超松弛迭代方法( (算例算例) )v取x(0)=0,=1.3,終止準(zhǔn)則為|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-

31、0.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我們來觀察松弛因子對(duì)收斂速度的影響.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.

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

33、值模擬一、二維水動(dòng)力數(shù)值模擬系統(tǒng)的分類1、按差分網(wǎng)格分：三角形、正方形、矩形、四邊形、多邊形、曲線坐標(biāo)網(wǎng)格以及各種形狀網(wǎng)格的組合2、按計(jì)算方法分：顯式法、隱式法、顯隱混合法3、按模擬格式分：三角元法、ADI法、破開算子法、單元體積法、MADI法、準(zhǔn)分析法、貼體坐標(biāo)法。二、平面二維水動(dòng)力數(shù)學(xué)模型的一般形式000)()(00hhfuygyvvxvutvhhfvxgyuvxuutuyvhxuhtsybysxbxwwwassywwwassxbybxvvufuvufhvvunhuvun22223/12223/1222v定解條件0( , ,0)0( , ,0)0( , ,0)( , )0u x yv x

34、yx yx y岸邊界：法向流速為 ；水邊界：給定潮位過程；三、三、ADIADI法法v交叉方向隱式格式是對(duì)x,y方向交叉使用隱式格式（同時(shí)也交叉使用了顯式格式），使得求解過程簡(jiǎn)化。222221/221/21211/221/2112212222u()()/ 2()/ 2nnnnnnnnuuatxyuuatuuatuxuy是的差分算子；是的差分算子；1、網(wǎng)格正方形或矩形，變量u,v,分別交錯(cuò)布置于網(wǎng)格中心和兩側(cè)。2、ADI基本思想（1）分步（2）交錯(cuò)顯隱ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位 、水深、水深uv1/21/2,1/2,3/2,1/2,1/2,11/2

35、,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/2,1/2,1/2,21/

36、2,1/2,1/2,1/2,( )0()()nnij ijsx ijnnijijijijuC hh3、差分格式X向運(yùn)動(dòng)方程在（i+1/2,j)點(diǎn)離散3、差分格式X向運(yùn)動(dòng)方程在（i+1/2,j)點(diǎn)離散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/21,1/20001.1,1()41()21()2

37、ijiji jnijnnnnniji ji jijijnnnijiji jvvvvvvhhh001/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)點(diǎn)離散1/21/21/21/2,1/2,01/2,01/2,1/20,1/2,1/20,1/2();21;();2()()2nnniijii jiijiniijiniijnnnnnii ji j

38、i ji ji jaubcudtahxbtchxtdvhvhy 其中聯(lián)立上面的式子得到下面的線性方程組injiinjiinjiidcuba2/1, 12/1, 2/12/1, injiinjiinjiiducbua2/1, 2/12/1,2/1, 2/1其中，u與存在如下關(guān)系iiiiiiiiiiiiiiiiiiiiiiiiiniiniiniinibEaFadHbEacGbGaHadFbGacEFuEHGu2/12/12/12/12/12/12/12/12/12/12/12/12/1;其中Y向運(yùn)動(dòng)方程在（i,j1/2)點(diǎn)離散)(412)(2)(4;)()()(2)(411, 2/11, 2/1,

39、 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其中v在后半個(gè)時(shí)間步長(zhǎng)內(nèi)，按上述同樣原理，在y向掃描。因此只須將上述各計(jì)算公式做如下的變換：xy,xy;ij,ij; uv,uv;即可求出vi,j+1/2n+1, i,jn+1,然后顯式求ui+1/2,jn+1。四、分步全隱式格式將x向動(dòng)量方程與連續(xù)方程聯(lián)立，求出un+

40、1,n+1/2,將y向動(dòng)量方程與連續(xù)方程聯(lián)立，求出vn+1,n+1,ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位水深水深uvX向動(dòng)量方程在（i+1/2,j)點(diǎn)離散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/1, 2/12/1, 2/1, 2/1, 11, 2/1, 2/11, 2/1njinjinjisxnjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihv fhCvuugxg

41、yuuvxuuutuu 連續(xù)方程在（i,j)點(diǎn)離散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 )(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

42、, 2/3, 2/1, 1njinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihnCvvvvvuuuuuuuuuuuuuuu injiinjiinjiidcuba2/1, 11, 2/12/1, injiinjiinjiiducbua1, 2/12/1,1, 2/1 以上離散后的公式整理后如下:同理 在y方向，y向動(dòng)量方程和連續(xù)方程聯(lián)立得如下：injiinjiinjiidcvba11,12/1,1, injiinjiinjiidvcbva12/1,1,12/1, 12/1,nnu 11,nn

43、v 五、移步雙向交替顯、隱式交錯(cuò)法MADI法：將水深、水位、流速等變化量均布置在同一網(wǎng)格節(jié)點(diǎn)上，由此將基本方程離散成新的差分代數(shù)方程組，并建立一種新的解法，這種解法既吸收了原有傳統(tǒng)的ADI法的優(yōu)點(diǎn)，又有較高的穩(wěn)定性、收斂性和精度。方程離散時(shí)，時(shí)間導(dǎo)數(shù)項(xiàng)采用前差表示，空間采用中心差分，將t分成兩等分。v當(dāng)nt(n+1/2)t,x向動(dòng)量方程和連續(xù)方程建立差分方程iniiniiniiininiiniidcubadcbua2/112/12/112/112/12/112/112/11211111,00niniiiiiiiiiiiiiudduucbacba2/12/111211111,00niniiiii

44、iiiiiiiiudduucbacba)()(2)(21)(2)2(22)()()(2)(4121,01,1,01, 10, 10,1,1,02,2,2, 1, 1njinjinjinjinjiinjiiinjiinjinjinjinjiiinjijinjinjinjinjiiihvhvytdhxtcbhxtavfyuutudxtgchCvutguuxtbxtga)()()()()()(1 ()(1(1)(1 ()(1(22/1,02,2,22/1,1,1,) 1,() 1,(2/1, 1, 12/1,2/1,vsignKusignEhCvugvvvKvvKytNygfuvvEvvEuxtvM

45、NMvvunjijinjinjinjinjinjivnjinjivjinjibynjibynjinjinjiunjinjiunjinjijijijinjiv當(dāng)（n+1/2）t(n+1)t,y向動(dòng)量方程和連續(xù)方程建立差分方程,求解v,顯式求u，過程同上。六、三角元法由于采用矩形網(wǎng)格而形成鋸齒形岸線以及鋸齒形堤、壩必然給數(shù)值模擬結(jié)果帶來不良影響，因此有限差分法的矩形網(wǎng)格，難以準(zhǔn)確模擬不規(guī)則曲線形岸線，即使采用空間變步長(zhǎng)網(wǎng)格，也不能完全模擬復(fù)雜的曲線形邊界。為了解決這些問題，采用了三角形網(wǎng)格。這種不規(guī)則三角形網(wǎng)格有以下優(yōu)點(diǎn)：v可以隨意加密計(jì)算網(wǎng)點(diǎn)形成不規(guī)則三角形，從而較準(zhǔn)確地模擬出復(fù)雜的邊界，如：岸

46、線、建筑物輪廓及航道，v可根據(jù)模擬區(qū)域的重要性，期F好網(wǎng)格的疏密程度及漸變程度。重要地區(qū)，網(wǎng)格密些，不重要地區(qū)，網(wǎng)格疏些；000)()(hfuygyvvxvutvhfvxgyuvxuutuyhvxhutbybx 基本方程離散后方程0)()()()(0)()()()(0)()(11111111111nibynininininininininibxninininininininininininihfuygyvvxvutvvhfvxgyuvxuutuuyhvxhut iiiyaxaaF210yaxaayxF210),(jjjyaxaaF210kkkyaxaaF210ekkkkjjjjiiiisFyc

47、xbaFycxbaFycxbayxF2/ )()()(),(ekkkkejjjjeiiiisycxbaNsycxbaNsycxbaN2/ )(2/ )(2/ )(kkjjiiFNFNFNyxF),(kkjjiikkjjiiFcFcFcyFFbFbFbxFekkjjiiSFbFbFbMyF2/ )()(ekkjjiiSFcFcFcMxF2/ )()(eeSS出于該格式采用顯式求解。同時(shí)受到復(fù)雜地形的影響，在計(jì)算過程中會(huì)產(chǎn)生一些擾。當(dāng)這些擾動(dòng)擴(kuò)大并傳播易使計(jì)算失敗為消除這些擾動(dòng)的影響，采用濾波公式mjmjjjiiiLLFFF1122*)()1( 七、邊界處理邊界處理合適與否影響數(shù)值模擬的成敗。實(shí)際

48、工程中，計(jì)算域的邊界常由不規(guī)則的曲線組成，如何利用有限差分的矩形網(wǎng)格來模擬曲線邊界，以及如何選取邊界值，這是邊界處理的重要內(nèi)容。1、邊界類型（1）Dirichlet邊界（2）Neumann邊界（3）淺灘活動(dòng)邊界2、 Dirichlet邊界在邊界處有已知的函數(shù)值。在實(shí)際情況下，水邊界通常屬于這類邊界，該處的實(shí)測(cè)水位、流速，可作為已知函數(shù)值。1ffv如果邊界網(wǎng)格結(jié)點(diǎn)正好在實(shí)測(cè)點(diǎn)上，則該結(jié)點(diǎn)的邊界值可直接采用已知值。v如果邊界網(wǎng)格結(jié)點(diǎn)不在實(shí)測(cè)點(diǎn)上，則分以下不同情況分別處理。v邊界結(jié)點(diǎn)與實(shí)測(cè)點(diǎn)很靠近，邊界值直接等于實(shí)測(cè)值；v如果實(shí)測(cè)點(diǎn)與邊界點(diǎn)較遠(yuǎn)，但與虛擬外結(jié)點(diǎn)較近，則令虛擬外結(jié)點(diǎn)直接等于實(shí)測(cè)值，然后

49、引入邊界條件求出邊界值。v邊界網(wǎng)格結(jié)點(diǎn)或虛擬結(jié)點(diǎn)與實(shí)測(cè)點(diǎn)均不靠近，不能直接引用實(shí)測(cè)值，這時(shí)，可根據(jù)前移時(shí)間步長(zhǎng)所得的結(jié)果，用線性差值可得邊界值。xiicxjijixiicjijijixxfffxxffff )(, 1,)(, 1, 13、Neumann邊界在實(shí)際情況下，固定邊界屬于此類邊界0()Sfffnv邊界正好通過網(wǎng)格或邊界與網(wǎng)格結(jié)點(diǎn)很靠近v邊界的外法線方向與坐標(biāo)軸平行，直接得到邊界值；以A點(diǎn)位例： 對(duì)于c點(diǎn)，可得 v邊界的外法線方向與坐標(biāo)軸不平行，考慮外法線與坐標(biāo)軸的夾角，帶入邊界條件后離散得到邊界值。 v邊界與網(wǎng)格結(jié)點(diǎn)距離較大4、淺灘活動(dòng)邊界（1）開挖法：將灘地開挖至可能出現(xiàn)的最低水位

50、之下，為使水量平衡，將岸邊界向水域內(nèi)移動(dòng)，并增大開挖部位的糙率以求得動(dòng)量上的平衡。這種方法一般只適用于潮灘問題，而且灘地面積只占整個(gè)海區(qū)較小的情況，而對(duì)于水位變化小，坡度較緩的地形，這種處理易失真。（2）凍結(jié)法根據(jù)水深判斷網(wǎng)格結(jié)點(diǎn)是否露出水面，對(duì)談度單元取糙率系數(shù)為一接近于無窮大的正數(shù)（糙率和水位均布置在網(wǎng)格中心），使單元四周的流速為一接近于0的無窮小量，這樣的處理結(jié)果相當(dāng)于使該單元的潮位在計(jì)算過程中北凍結(jié)不變這種方法適用于寬淺，坡度較坦的露灘問題，而對(duì)于潮灘相間的海岸、河口海域，則因水量和動(dòng)量的過分凍結(jié)而失真。（3）切削法稱水位判別法或薄層水法，該法對(duì)露灘單元并不凍結(jié)，而引入一個(gè)富裕水深（相

51、當(dāng)于灘地上存在很薄的水層）以保證計(jì)算過程的完整，相當(dāng)于將原始地形切削降低，而一旦判斷實(shí)際水深大于富裕水深時(shí)，恢復(fù)原始地形，和開挖法具有相同的局限性。（4）窄縫法假想在岸灘的每個(gè)網(wǎng)格上存在一條很窄的縫隙，它的深度和岸灘前的水深一致，根據(jù)水量平衡，將窄縫內(nèi)的水量平鋪到岸灘上，把計(jì)算邊界設(shè)在岸灘的窄縫內(nèi)，成為具有一定水深的固定邊界。系數(shù)不易確定，只適用于岸邊界的露灘問題。bbzzsszzBzzeBBBzBb0)(0)()( B(z) 為單位淺灘長(zhǎng)度內(nèi)的窄縫系數(shù)為單位淺灘長(zhǎng)度內(nèi)的窄縫系數(shù)2)()(2zBzBAv 窄縫法的優(yōu)點(diǎn)是在計(jì)算中、不必隨時(shí)變化邊界計(jì)算點(diǎn)只礙事先定出一個(gè)足夠大的固定計(jì)算邊界，就可以

52、動(dòng)態(tài)地模擬出計(jì)算城內(nèi)的水位、流速變化。由于計(jì)算域事先已固定，求解時(shí)的系數(shù)矩陣大小不變，有利于計(jì)算的穩(wěn)定性。又由于窄縫寬度很小。續(xù)內(nèi)流速也可人為加以衰減，故對(duì)精度無大影響。缺點(diǎn)是窄縫參數(shù)不易確定。（5）干濕法根據(jù)每步計(jì)算結(jié)果判斷每個(gè)單元干、濕？濕單元參加方程的計(jì)算。二維淺水水流的一種三角形網(wǎng)格二維淺水水流的一種三角形網(wǎng)格FVM計(jì)算格式計(jì)算格式對(duì)開闊寬淺型水域的水流運(yùn)動(dòng), 可以守恒型二維非恒定淺水方程描述:ByGxFtQTbysybxsxTTTyZghfhuxZghfhvByvhghhhvvyuhhuvhvGxvhhuvxuhghhhuuhuFhvhuQ), 0()21,(),21,(),( 式中

53、: x , y 空間坐標(biāo); t 時(shí)間坐標(biāo); u, v 在x , y 方向沿水深積分平均流速分量; 潮位; h 水深; f 柯氏力系數(shù); 水的密度; Z 地形高程; x , y 方向沿水深平均的紊動(dòng)粘性系數(shù); sx， sy 沿x , y 方向的風(fēng)應(yīng)力; bx， by沿x , y 方向的河床底應(yīng)力.閉邊界( 岸邊界) : 采用流動(dòng)法向通量為零開邊界( 水流邊界) 一般在河道較順直段選取已知流速或水位過程。對(duì)計(jì)算區(qū)域中的任意一個(gè)三角形單元, 各水力變量P ( x , y ) 的計(jì)算值均布置在節(jié)點(diǎn)上, 并假定單元內(nèi)變量分布呈一階近似:P( x , y ) = a + bx + cyv基本方程的離散及求

54、解對(duì)基本方程在控制體積( 圖1) 上積分:0)()(GdxFdydxdyBtQ 由于引入一階近似假定, 單元內(nèi)各水力變量, h, u, v 均可表示為節(jié)點(diǎn)處變量值的線性函數(shù), 從而可直接求解積分方程. 離散結(jié)果表示為如下形式:183116bhStjjji()(buwuSCvugthjjjuijjiiiij 38311083222)()(bvwvSCvugthjjjvijjiiiij )()(4)()()(4)(11112011111301jjjijjjijjjjjijjjijjjjjxxxxyyyyShhxxxxyyyyShhw第四章第四章 河口三維流體動(dòng)力學(xué)模型河口

55、三維流體動(dòng)力學(xué)模型v在寬闊且較深的海岸河口地區(qū)，研究水流運(yùn)動(dòng)，海岸演變及泥沙運(yùn)動(dòng)時(shí)，通常二位數(shù)值模擬就不能滿足要求了，此外，像疏浚拋泥、油膜運(yùn)動(dòng)、水質(zhì)污染擴(kuò)散等一些專門課題，不是二位數(shù)值模擬能解決的，必須采用三位數(shù)值模擬技術(shù)。vPOM（Princeton Ocean Model)模型由Blumberg和Mellor1978提出，經(jīng)多年的改進(jìn)，已成為比較廣泛使用的海洋模式。vPOM在淺水海域水深小于3米時(shí)，退潮時(shí)，模擬灘地退出水面遇到困難，計(jì)算不穩(wěn)定。vPOM模型采用模式分離技術(shù)，三維控制方程組及其定解條件構(gòu)成模型的內(nèi)模式，而外模由全積分內(nèi)模式方程得到。第四章第四章 河口三維流體動(dòng)力學(xué)模型河口三

56、維流體動(dòng)力學(xué)模型v連續(xù)方程：0zwyvxuv動(dòng)量方程：xMFzuKzxpfvzuwyuvxuutu)(10yMFzvKzypfuzvwyvvxvutv)(10gzpMK為動(dòng)量垂向渦粘系數(shù) v對(duì)于三維斜壓模型，還需要同時(shí)考慮溫度、鹽度的擴(kuò)散過程，其控制方程為：THFzTKzzTwyTvxTutT)(SHFzSKzzSwySvxSutS)(),(pTS T為勢(shì)溫（對(duì)河口及近岸地區(qū)可為現(xiàn)場(chǎng)實(shí)際溫度），S為鹽度，HK為反映溫度、鹽度垂向紊動(dòng)混合的垂向擴(kuò)散系數(shù)。xvyuAyxuAxFMMx2xvyuAxyvAyFMMy2ySTAyxSTAxFHHST, xFyFTFSF是對(duì)模型網(wǎng)格無法分辨的所謂次網(wǎng)格運(yùn)

57、動(dòng)過程用水平紊動(dòng)擴(kuò)散過程參數(shù)化后的產(chǎn)生項(xiàng)，分別為：對(duì)基本方程簡(jiǎn)化時(shí)采用的假定與近似對(duì)基本方程簡(jiǎn)化時(shí)采用的假定與近似 （1）靜壓假定：在河口、近岸淺水地區(qū)，垂向速度的時(shí)間變率（即：垂向加速度）與重力加速度相比甚微，可略去不計(jì)，因此垂向動(dòng)量方程可簡(jiǎn)化為：，即壓強(qiáng)沿水深的變化符合靜水壓強(qiáng)分布。（2）Boussinesq近似：海水密度為時(shí)均值（參考密度）和脈動(dòng)值之和，將其代入動(dòng)量方程后，除在重力加速度的前面保留外，其余各項(xiàng)的均略去。（3）Boussinesq假定：由于在時(shí)均運(yùn)動(dòng)方程中包含了較難處理的雷諾應(yīng)力張量，Boussinesq在1877年提出了關(guān)于可以將水流紊動(dòng)應(yīng)力類比于層流粘性應(yīng)力的假定，即用

58、層流粘性應(yīng)力的形式對(duì)紊動(dòng)應(yīng)力進(jìn)行參數(shù)化。v坐標(biāo)系中的動(dòng)力學(xué)方程： Hz010tyDVxDUxMFUDKdxDDxgDxgDfVDUyUVDxDUtUD0022yMFVDKdyDDygDygDfUDVyDVxUVDtVD0022ttDyyDVxxDUWv坐標(biāo)系中的溫度、鹽度擴(kuò)散方程分別為 THFTDKTyTVDxTUDtTD)(SHFSDKSySVDxSUDtSD)(v紊流閉合模型：得到動(dòng)量、標(biāo)量垂向擴(kuò)散系數(shù)KM、KH qHMqFlBDqKgVUDKqDKqyDVqxDUqtDq1302222222222lHMqFWBDqKgEVUDKlElqDKlqylDVqxlDUqtlDq13032212

59、22220yDVxDUtddxDxDgDGxgDDVfFyDVUxDUtDUxbxsxx01002ddyDyDgDGygDDUfFyDVxDVUtDVybysyy01002v外模式方程外模式方程v邊界條件v1、自由表面邊界條件sysxMVUDK,0SHSTDKH,023212sUBq02lq0 yxz,0 v2 2、水體底部邊界條件、水體底部邊界條件yxHz,1bybxMVUDK,023212bUBq0,ST02lq0 v3 3、側(cè)向閉邊界條件、側(cè)向閉邊界條件0nU0,STn v4、側(cè)向開邊界條件v對(duì)水位開邊界，通常用實(shí)測(cè)的水位資料或者用更大范圍數(shù)學(xué)模型計(jì)算的水位值作為強(qiáng)迫水位控制條件。v對(duì)于

60、溫、鹽等守恒性標(biāo)量物質(zhì)的開邊界條件，通常分為入流邊界和出流邊界兩種情況。對(duì)于入流邊界，一般采用開邊界實(shí)測(cè)的溫、鹽數(shù)據(jù)；對(duì)于出流邊界，采用對(duì)流型開邊界形式：v數(shù)值求解方法（1）模式分裂技術(shù) 將垂向積分的運(yùn)動(dòng)方程（外模態(tài)）從反映流速垂向結(jié)構(gòu)的運(yùn)動(dòng)方程（內(nèi)模態(tài)）中分離出來，用較少的計(jì)算量通過求解外模態(tài)得到自由表面，然后通過求解內(nèi)模態(tài)得到水流的垂向結(jié)構(gòu)，這稱為模式分裂技術(shù)。 （2）時(shí)間分步法 將三維內(nèi)模態(tài)的計(jì)算分為計(jì)算垂向擴(kuò)散過程時(shí)步和計(jì)算對(duì)流和水平擴(kuò)散時(shí)步兩部分。分步的目的是為了提高計(jì)算效率和提高垂向分辨率兩大要求。前一部分為了適應(yīng)垂向網(wǎng)格的高分辨率需要單獨(dú)采用隱格式，后者則單獨(dú)采用顯格式。 ECO