1、河口、海岸水动力模河口、海岸水动力模拟技术拟技术第一章第一章 绪论绪论v海岸：是海陆相互作用的重要地带，也是海、陆、气交互作用的重要空间，这种表现在： 岸线演变（自然和人为） 飓风（台风）带来的灾难性破坏； 海洋潮汐环境的变化。v河口：海岸常伴随有江河湖泊的出海口，通常称为河口。v海岸河口问题： 潮流问题 波浪问题 径流、异重流（密度流）、污染物（COD)扩散。v研究海岸河口问题的方法 物理模型（水力学比尺模型） 数学模型（数值模拟）沿岸过程动力因素物质过程流（潮流）波（风浪）盐水入侵泥沙输移污染物扩散波流相互作用海水入侵控制反馈流载波波生流v数值模拟：一门综合性的模拟技术，它采用数学模型来模
2、拟某中物理现象，并通过计算机用数值计算法进行近似求解，籍以复演自然演变过程的总称。v水力学、泥沙数值模拟：以水力学和泥沙动力学为理论基础，并结合具体工程的一门新型实用科学。v水动力泥沙数值模拟：以微分方程为理论，并通过微分方程的离散，变成代数方程，最后采用计算机进行近似求解。v数值模拟的特点： （1）一般以线性理论为基础，但实际自然现象和描述这些现象的微分方程均为非线性的； （2）需要丰富的经验，现场资料和一定的技巧； （3）数值模拟不仅仅是一种近似计算，可以作为一种实验或研究及预测方法。v数值模拟的优点： （1）实验费用少； （2）速度快、周期短； （3）可以模拟多种因素相互作用的复杂物理过
3、程。如可以模拟水（潮）流、风、柯氏力等多种因素共同作用下的多种泥沙及地形演变的复杂过程。 （4）可以完全控制流体的物理性质（如密度、容重、粘度、含沙量等） （5）模型建成后，长期保存、随时调用修改。 （6）无法模拟微分方程不能描述的物理现象。v数值模拟工作的基本步骤（1）建立数学模型和编制源程序 建立或选择的微分方程； 根据模拟域边界条件选择合适的网格； 按一定的格式离散方程，得到代数方程和采用合适的数值方法求解代数方程； 编制源程序求解代数方程。 数值模拟分析（收敛性、稳定性、相容性、误差程度等）（2）调试源程序（3）模型验证 调整模型中有关参数（糙率、紊动动量掺混系数等），使模型有良好的稳
4、定性和收敛性，并与现场资料有良好的吻合；（4）正式方案试验 v河口、海岸水动力模拟的发展方向 1、资料同化将是河口数值模型发展和结合的一个新技术切入点,也是带动河口动力数字模拟技术革新的一种重要方法.2、数字河口动力模型1、河口模型四维资料同化v四维同化在河口中的一类主要应用是为河口数值模式提供优化的初始状态,从而校正不合理的边界条件所带来的偏差,提高模拟精度.v变分同化技术还可应用于确定河口模式中的未知参数. 例如确定了潮汐河口的摩擦系数,通常这种系数是通过经验获得的.v另外,河口模式的外部强迫场(如风场、海-气界面的热通量等)都可以通过这种技术反演出来.v资料同化的关键在同化模型在过滤和插
5、值观测资料的机制方面有高效的预测能力。模型应对前期信息在时间上能向前插值。v动力场的模型同化熟路也被证明对陆-气系统内部环流的科学研究极有价值v同化是唯一能产生复杂非线性过长的场的手段，并且在物理学与动力学的机制上协调一致。2、数字河口动力模型v河口水动力特征及外界强迫作用因子构建数字河口动力模型(模型计算要素象温、盐、流、泥沙的特征均以数字矩阵形式记载),若再与数字河网模型嵌套联结,最终可以获得区域河网河口的动力-沉积-地貌的机制解释,揭示河口物理规律和解决工程实际问题.v在数字河口模型的构架下,可以将人类已经拥有的河口科学理论、知识与具有较强物理概念的水动力学模型集成于一体,为河口数值模型
6、计算提供良好的平台。数字河口动力模型具有许多优势:首先,数字河口模型是基于数字区域地形构建而成的,地形要素可自动生成,无需手工操作,大大提高了工作效率;其次,数字模型不仅能输出传统模型的结果,而且能够十分方便地给出河口水文要素和水文状态变量的空间分布场,这些对近岸河口动力科学研究与河口、港口、航道工程都有着广阔的应用前景.数字河口模型研究的最终目的就是利用已有的河口基础科学理论和知识,在数字区域地形的基础之上将观测点的水文信息拓展、同化至区域平面上乃至区域三维立体上的信息,并形成数字成品。v参考文献：vKoutitar 著“Mathematical Model in Coastal Engin
7、eering”1)模型简单易懂2）附有Basic程序，而且有验证的算例3）介绍各种数值处理技术v曹祖德、王运洪”水动力泥沙数值模拟第二章第二章 水动力数值模拟的理论基础水动力数值模拟的理论基础)()()()(022222222zVAzyVxVAygfUdtDVzUAzyUxUAxgfVdtDUzWyVxUzHzH ),()tyxzyVxUtW2.1 基本方程自由面运动学边界条件：底部运动学边界条件：),()(yxhzyhVxhUWzWyVxUtdtDU,V,W为x,y,z 方向上的流速分量。（，）为距平均海平面的自由表面水位。（，）为平均海平面距底部边界的水深。为水平扩散系数。为垂直涡动系数。
8、0)()()()(222222yQxQtQDCVUgygDfQQVyQUxtQQDCVUgxgDfQQVyQUxtQyxyxyyyxyxxx初始条件),()0 ,(0yxuyxu),()0 ,(0yxvyxv),()0 ,(0yxyx边界条件岸边界：法向流速为零。水边界：给定潮位过程。0)()()(0oxbxbbBxAZhgxZhgAxuQtQxQtASaint Venant 方程A 过水断面面积，Q为流量，Q=Bhu，B为河宽，h为水深，u为断面平均流速。 分别为水面剪切力与水底剪切力Zb 为水底的竖向坐标位置X轴为沿水流纵向方向。 三、二、一维方程的定解条件三、二、一维方程的定解条件v初始
9、条件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数值计算数值计算v在计算水动力、泥沙数值模拟时，大都将基本方程组离散成代数方程组，最后求解代数方程组，此处介绍微分方程组的离散技术有限差分法和线性代数方程组的数值解法。2.2.12.2.1有限差分法有限差分法v有限差分法是工程中常用的一种离散技术，将计算域分成有限个网格，通过差分法求网格结点的微分方程的近似值，也称网格法。v将网格结点上的函数f(x,
10、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()(4443333433332xxxxxxxxfxxfxxxxfxxfxxf（4）二阶中心差分10 , 2/10 , 2/)()(96)2/()(2)2()2()(6665554644542222xxxxxxxxfxxfxxxf
11、xxfxxfxxf2、几种常见的差分格式以一维热传导方程为例：022xftfxti-1 i i+1n+1 n（1）古典显式格式0)(222111xtOxffftffninininini022111xffftffninininini2111/),()21 (xtrffrfrfnininini可由已知值直接求解未知值，且只涉及二层，称为双层显式差分格式（2）古典隐式格式0)(222111111xtOxffftffninininini022111111xffftffninininini21111/,)21 (xtrfrffrrfnininini不能直接求解未知值，且只涉及二层，称为双层隐式差分格式（
12、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(2xtrfffrffninininini0
13、2221111xffftffninininini（5）加权六点格式，隐式格式在x点和n+n时层，01，对t和x均采用中心差分2/10)(2)1 (222112111111xtOxfffxffftffnininininininini21111111/)1 ()1 ()1 (21 )21 (xtrrfrffrrffrrfnininininini02)1 (22112111111xfffxffftffnininininininini2.2.22.2.2线性方程组的数值解线性方程组的数值解v有限差分法是工程中常用的一种离散技术，将计算域分成有限个网格，通过差分法求网格结点的微分方程的近似值，也称网格法
14、。v将网格结点上的函数f(x,y,z,t)表示成 ，i,j,k分别表示x,y,z方向的坐标位置，n表示时间。nkjif,1、解线性方程组的两种方法:v直接法:通过有限步算术运算直接求出方程组的精确解,最常用的是消元结合代入的方法.v实际上除非是采用无穷位精度计算,一般都得不到精确解.v直接方法适用于解低阶稠密矩阵方程组.v迭代法 类似于方程求根的迭代法,用一个迭代过程逐步逼近方程组的解.v迭代有可能不收敛,或虽然收敛,但收敛速度慢.v迭代法适用于求解高阶稀疏矩阵方程组.v稀疏矩阵:矩阵非零元素较少,且在固定的位置上.v稀疏矩阵一般是人为构造的,例如36页三转角插值时方程组(8.12),(8.1
15、5)的系数矩阵.GaussGauss消去法消去法( (第一次消元第一次消元) )v考虑方程组A(1)x=b(1)11111111221111112112222211111122( )( )( )( )( )( )( )( )( )( )( )( )nnnnnnnnnnaxaxaxbaxaxaxbaxaxaxb v第一次消元用第一个方程将后面方程的x1消去.111111( )( )iiama v计算乘数v条件:a11(1)0v用-mi1乘以第一个方程加到第i个(i=1,n)方程上,则消去了第i个方程中的x1.GaussGauss消去法消去法( (第一次消元第一次消元) )v经过上述过程,得到方程
16、组A(2)x=b(2),1111111211122222222222200( )( )( )( )( )( )( )( )( )( )nnnnnnnaaabxxaabxbaa v其中2111 1( )( )( ),ijijijaam a 2111 1( )( )( ),iiibbm b 111111( )( )iiama 2 3( , , )i jn GaussGauss消去法消去法( (第第k k次消元次消元) )v假设已完成k-1次消元,得到方程组A(k)x=b(k).111111112111122222221221111111000000000( )( )( )( )( ),( )( )
17、( )( ),()()(),( )( )( )( )kknkknkkkkkkkknkkkkknkknknnaaaaaaaaaaaaaaaa v第k次消元的目的是将akk(k) (称为主元)下面的元素变为0.GaussGauss消去法消去法( (第第k次消元次消元) )v对A(k)右下角的矩阵( )( )( )( )kkkkknkknknnaaaa( )( )kikikkkkama v计算乘数v条件:akk(k)0v用-mik乘以第k个方程加到第i个(i=k+1,n)方程上,则消去了第i个方程中的xk,得到方程组A(k+1)x=b(k+1).GaussGauss消去法消去法( (第第k次消元次消
18、元) )v第一步消元的计算公式2111 1( )( )( ),ijijijaam a 2111 1( )( )( ),iiibbm b 111111( )( )iiama v类似可以得到第k步消元的计算公式1()( )( ),kkkijijikkjaam a 1()( )( ),kkkiiikkbbm b ( )( )kikikkkkama 1( , )i jkn GaussGauss消去法消去法v消去法完成后最终得到与原方程组等价的三角形方程组A(n)x=b(n).1111111211122222222000( )( )( )( )( )( )( )( )( )nnnnnnnnaaabxxa
19、abxba v一共需进行 ? 步n-1GaussGauss消去法消去法( (算法算法) )1,1kn ( )( )/kkikikkkmaa 1,ikn1()( )( )kkkijijikkjaam a 1,i jkn 1()( )( )kkkiiikkbbm b 1,ikn( )( )/nnnnnnxba 1,1kn ( )( )( )1()/nkkkkkkjjkkj kxbaxa 追赶法求解三对角方程组追赶法求解三对角方程组v上面的方程组可以利用追赶法求解(P185).v对于下面形式的方程组11112222211111nnnnnnnnnfbcxabcxfabcxfabxf v将系数矩阵进行三
20、角分解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在得到系数矩阵的分解后,原方程组转化为vLUx=f.v先求解Ly=f11122223111nnnnnnnyfryfryfryf v显
21、然有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 两边取极限
22、得到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(k)
23、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)|(通常用无穷范数)的值来判断是
24、否终止迭代.v在上面的例子中,我们将第i个方程变形为左边是xi,右边是其它分量和常数的线性组合,然后进行迭代,这一方法称位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+f.
25、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,我们有理由认为已经计算出的第k+1步的x1比第k步的“好”.因此,我们应该用第k+1步的x1和第k步的x3来计算x2.
26、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 1123112
27、131113121322081433111633612()( )( )()()( )()()()()()()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.0001
28、63029.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(k+
29、1)为xi(k)与 的加权平均.1()kix 1111()( )()( )()( )()()kkkkkkiiiiiixxxxxx 111( )()( )()inkkkiiijjijjiijj ixba xa xa v在=1时,上述方法就是Gauss-Seidel方法,1时称为超松弛法(有时不管的范围,统称为超松弛方法).超松弛迭代方法超松弛迭代方法( (算例算例) )v对于方程组123441111141111141111141xxxx v松弛方法迭代格式为(1)( )( )( )( )( )111234(1)( )(1)( )( )( )221234(1)( )(1)(1)( )( )3312
30、34(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.
31、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步数301 156 104765947383126211.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
32、2.0步数1712121518243555114*v步数表示|x(k)x(k1)|10-5时的迭代步数,=2.0时,500步以内不收敛.超松弛迭代方法超松弛迭代方法( (矩阵表示矩阵表示) )v超松弛迭代格式可以写为1111()( )()( )()inkkkkiiiijjijjiijj ixxba xa xa 111111()( )()( )()()inkkkkiiiiiiiijjijjjj ia xa xba xa x v用矩阵可以表示为111()( )()( )()()kkkkDxDxbLxUx 11()( )()()kkDL xDU xb 第三章第三章 二维水动力数值模拟二维水动力数值模
33、拟一、二维水动力数值模拟系统的分类1、按差分网格分：三角形、正方形、矩形、四边形、多边形、曲线坐标网格以及各种形状网格的组合2、按计算方法分：显式法、隐式法、显隐混合法3、按模拟格式分：三角元法、ADI法、破开算子法、单元体积法、MADI法、准分析法、贴体坐标法。二、平面二维水动力数学模型的一般形式000)()(00hhfuygyvvxvutvhhfvxgyuvxuutuyvhxuhtsybysxbxwwwassywwwassxbybxvvufuvufhvvunhuvun22223/12223/1222v定解条件0( , ,0)0( , ,0)0( , ,0)( , )0u x yv x yx
34、 yx y岸边界：法向流速为 ；水边界：给定潮位过程；三、三、ADIADI法法v交叉方向隐式格式是对x,y方向交叉使用隐式格式（同时也交叉使用了显式格式），使得求解过程简化。222221/221/21211/221/2112212222u()()/ 2()/ 2nnnnnnnnuuatxyuuatuuatuxuy是的差分算子；是的差分算子；1、网格正方形或矩形，变量u,v,分别交错布置于网格中心和两侧。2、ADI基本思想（1）分步（2）交错显隐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,1
35、1/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/2,
36、1/2,1/2,1/2,( )0()()nnij ijsx ijnnijijijijuC hh3、差分格式X向运动方程在（i+1/2,j)点离散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/21,1/20001.1,1()41()21()2ij
37、iji 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连续方程在（i,j)点离散1/21/21/21/2,1/2,01/2,01/2,1/20,1/2,1/20,1/2();21;();2()()2nnniijii jiijiniijiniijnnnnnii ji ji
38、ji ji jaubcudtahxbtchxtdvhvhy 其中联立上面的式子得到下面的线性方程组injiinjiinjiidcuba2/1, 12/1, 2/12/1, injiinjiinjiiducbua2/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
39、/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在后半个时间步长内，按上述同样原理，在y向扫描。因此只须将上述各计算公式做如下的变换：xy,xy;ij,ij; uv,uv;即可求出vi,j+1/2n+1, i,jn+1,然后显式求ui+1/2,jn+1。四、分步全隐式格式将x向动量方程与连续方程联立，求出un+1,
40、n+1/2,将y向动量方程与连续方程联立，求出vn+1,n+1,ji+1/2i+1i-1i-1/2j+1j+1/2ij-1/2j-1水位水位水深水深uvX向动量方程在（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/1, 2/12/1, 2/1, 2/1, 11, 2/1, 2/11, 2/1njinjinjisxnjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjihv fhCvuugxgyu
41、uvxuuutuu 连续方程在（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 )(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向动量方程和连续方程联立得如下：injiinjiinjiidcvba11,12/1,1, injiinjiinjiidvcbva12/1,1,12/1, 12/1,nnu 11,nnv
43、五、移步双向交替显、隐式交错法MADI法：将水深、水位、流速等变化量均布置在同一网格节点上，由此将基本方程离散成新的差分代数方程组，并建立一种新的解法，这种解法既吸收了原有传统的ADI法的优点，又有较高的稳定性、收敛性和精度。方程离散时，时间导数项采用前差表示，空间采用中心差分，将t分成两等分。v当nt(n+1/2)t,x向动量方程和连续方程建立差分方程iniiniiniiininiiniidcubadcbua2/112/12/112/112/12/112/112/11211111,00niniiiiiiiiiiiiiudduucbacba2/12/111211111,00niniiiiiii
44、iiiiiiudduucbacba)()(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,vsignKusignEhCvugvvvKvvKytNygfuvvEvvEuxtvMNM
45、vvunjijinjinjinjinjinjivnjinjivjinjibynjibynjinjinjiunjinjiunjinjijijijinjiv当（n+1/2）t(n+1)t,y向动量方程和连续方程建立差分方程,求解v,显式求u，过程同上。六、三角元法由于采用矩形网格而形成锯齿形岸线以及锯齿形堤、坝必然给数值模拟结果带来不良影响，因此有限差分法的矩形网格，难以准确模拟不规则曲线形岸线，即使采用空间变步长网格，也不能完全模拟复杂的曲线形边界。为了解决这些问题，采用了三角形网格。这种不规则三角形网格有以下优点：v可以随意加密计算网点形成不规则三角形，从而较准确地模拟出复杂的边界，如：岸线、
46、建筑物轮廓及航道，v可根据模拟区域的重要性，期F好网格的疏密程度及渐变程度。重要地区，网格密些，不重要地区，网格疏些；000)()(hfuygyvvxvutvhfvxgyuvxuutuyhvxhutbybx 基本方程离散后方程0)()()()(0)()()()(0)()(11111111111nibynininininininininibxninininininininininininihfuygyvvxvutvvhfvxgyuvxuutuuyhvxhut iiiyaxaaF210yaxaayxF210),(jjjyaxaaF210kkkyaxaaF210ekkkkjjjjiiiisFycxb
47、aFycxbaFycxbayxF2/ )()()(),(ekkkkejjjjeiiiisycxbaNsycxbaNsycxbaN2/ )(2/ )(2/ )(kkjjiiFNFNFNyxF),(kkjjiikkjjiiFcFcFcyFFbFbFbxFekkjjiiSFbFbFbMyF2/ )()(ekkjjiiSFcFcFcMxF2/ )()(eeSS出于该格式采用显式求解。同时受到复杂地形的影响，在计算过程中会产生一些扰。当这些扰动扩大并传播易使计算失败为消除这些扰动的影响，采用滤波公式mjmjjjiiiLLFFF1122*)()1( 七、边界处理边界处理合适与否影响数值模拟的成败。实际工程
48、中，计算域的边界常由不规则的曲线组成，如何利用有限差分的矩形网格来模拟曲线边界，以及如何选取边界值，这是边界处理的重要内容。1、边界类型（1）Dirichlet边界（2）Neumann边界（3）浅滩活动边界2、 Dirichlet边界在边界处有已知的函数值。在实际情况下，水边界通常属于这类边界，该处的实测水位、流速，可作为已知函数值。1ffv如果边界网格结点正好在实测点上，则该结点的边界值可直接采用已知值。v如果边界网格结点不在实测点上，则分以下不同情况分别处理。v边界结点与实测点很靠近，边界值直接等于实测值；v如果实测点与边界点较远，但与虚拟外结点较近，则令虚拟外结点直接等于实测值，然后引入
49、边界条件求出边界值。v边界网格结点或虚拟结点与实测点均不靠近，不能直接引用实测值，这时，可根据前移时间步长所得的结果，用线性差值可得边界值。xiicxjijixiicjijijixxfffxxffff )(, 1,)(, 1, 13、Neumann边界在实际情况下，固定边界属于此类边界0()Sfffnv边界正好通过网格或边界与网格结点很靠近v边界的外法线方向与坐标轴平行，直接得到边界值；以A点位例： 对于c点，可得 v边界的外法线方向与坐标轴不平行，考虑外法线与坐标轴的夹角，带入边界条件后离散得到边界值。 v边界与网格结点距离较大4、浅滩活动边界（1）开挖法：将滩地开挖至可能出现的最低水位之下
50、，为使水量平衡，将岸边界向水域内移动，并增大开挖部位的糙率以求得动量上的平衡。这种方法一般只适用于潮滩问题，而且滩地面积只占整个海区较小的情况，而对于水位变化小，坡度较缓的地形，这种处理易失真。（2）冻结法根据水深判断网格结点是否露出水面，对谈度单元取糙率系数为一接近于无穷大的正数（糙率和水位均布置在网格中心），使单元四周的流速为一接近于0的无穷小量，这样的处理结果相当于使该单元的潮位在计算过程中北冻结不变这种方法适用于宽浅，坡度较坦的露滩问题，而对于潮滩相间的海岸、河口海域，则因水量和动量的过分冻结而失真。（3）切削法称水位判别法或薄层水法，该法对露滩单元并不冻结，而引入一个富裕水深（相当于
