《河口、海岸水动力模拟技术海岸：是海陆相互作用的重要地带，也是海、陆、气交互作用的重要空间》由会员分享，可在线阅读，更多相关《河口、海岸水动力模拟技术海岸：是海陆相互作用的重要地带，也是海、陆、气交互作用的重要空间（159页珍藏版）》请在金锄头文库上搜索。
1、河口、海岸水动力模拟技术,第一章 绪论,海岸：是海陆相互作用的重要地带，也是海、陆、气交互作用的重要空间，这种表现在：岸线演变（自然和人为）飓风（台风）带来的灾难性破坏；海洋潮汐环境的变化。 河口：海岸常伴随有江河湖泊的出海口，通常称为河口。,海岸河口问题：潮流问题波浪问题径流、异重流（密度流）、污染物（COD)扩散。 研究海岸河口问题的方法物理模型（水力学比尺模型）数学模型（数值模拟）,沿岸过程,动力因素,物质过程,流（潮流）,波（风浪）,盐水入侵,泥沙输移,污染物扩散,波流相互作用,海水入侵,控制,反馈,流载波,波生流,数值模拟：一门综合性的模拟技术，它采用数学模型来模拟某中物理现象，并通
2、过计算机用数值计算法进行近似求解，籍以复演自然演变过程的总称。水力学、泥沙数值模拟：以水力学和泥沙动力学为理论基础，并结合具体工程的一门新型实用科学。,水动力泥沙数值模拟：以微分方程为理论，并通过微分方程的离散，变成代数方程，最后采用计算机进行近似求解。 数值模拟的特点：（1）一般以线性理论为基础，但实际自然现象和描述这些现象的微分方程均为非线性的；（2）需要丰富的经验，现场资料和一定的技巧；（3）数值模拟不仅仅是一种近似计算，可以作为一种实验或研究及预测方法。,数值模拟的优点：（1）实验费用少；（2）速度快、周期短；（3）可以模拟多种因素相互作用的复杂物理过程。如可以模拟水（潮）流、风、柯氏
3、力等多种因素共同作用下的多种泥沙及地形演变的复杂过程。（4）可以完全控制流体的物理性质（如密度、容重、粘度、含沙量等）（5）模型建成后，长期保存、随时调用修改。（6）无法模拟微分方程不能描述的物理现象。,数值模拟工作的基本步骤 （1）建立数学模型和编制源程序建立或选择的微分方程；根据模拟域边界条件选择合适的网格；按一定的格式离散方程，得到代数方程和采用合适的数值方法求解代数方程；编制源程序求解代数方程。数值模拟分析（收敛性、稳定性、相容性、误差程度等）,（2）调试源程序 （3）模型验证调整模型中有关参数（糙率、紊动动量掺混系数等），使模型有良好的稳定性和收敛性，并与现场资料有良好的吻合； （4
4、）正式方案试验,河口、海岸水动力模拟的发展方向 1、河口模型四维资料同化 2、数字河口动力模型 数字河口动力模型具有许多优势: 首先,数字河口模型是基于数字区域地形构建而成的,地形要素可自动生成,无需手工操作,大大提高了工作效率; 其次,数字模型不仅能输出传统模型的结果,而且能够十分方便 地给出河口水文要素和水文状态变量的空间分布场,这些对近岸河口动力科学研究与河口、港口、航道工程都有着广阔的应用前景. 总而言之,数字河口模型研究的最终目的就是利用已有的河口基础科学理论和知识,在数字区域地形的基础之上将观测点的水文信息拓展、同化至区域平面上乃至区域三维立体上的信息,并形成数字成品,为国家宏观决
5、策和国民经济各行各业服务。,参考文献： Koutitar 著“Mathematical Model in Coastal Engineering” 1)模型简单易懂 2）附有Basic程序，而且有验证的算例 3）介绍各种数值处理技术 曹祖德、王运洪”水动力泥沙数值模拟,第二章 水动力数值模拟的理论基础,2.1 基本方程,自由面运动学边界条件：,底部运动学边界条件：,U,V,W为x,y,z 方向上的流速分量。 （，）为距平均海平面的自由表面水位。 （，）为平均海平面距底部边界的水深。 为水平扩散系数。 为垂直涡动系数。,初始条件,边界条件 岸边界：法向流速为零。,水边界：给定潮位过程。,Sain
6、t Venant 方程,三、二、一维方程的定解条件,初始条件 u,v,w,|t=0=u0,v0,w0,0 边界条件 开边界：计算域水体与外部水体相接处。 (u,v,w)=(u(t),v(t),w(t) =(t) 固边界：计算域与陆地或建筑物接壤处 无滑动：u,v,w=0 有滑动： 垂直边界的速度为0。,2.2数值计算,在计算水动力、泥沙数值模拟时，大都将基本方程组离散成代数方程组，最后求解代数方程组，此处介绍微分方程组的离散技术有限差分法和线性代数方程组的数值解法。,2.2.1有限差分法,有限差分法是工程中常用的一种离散技术，将计算域分成有限个网格，通过差分法求网格结点的微分方程的近似值，也称
7、网格法。 将网格结点上的函数f(x,y,z,t)表示成 ，i,j,k分别表示x,y,z方向的坐标位置，n表示时间。,1、工程中常用的几种差分和微分的关系（一维） （1）一阶向前差分,（2）一阶向后差分,（3）一阶中心差分,（4）二阶中心差分,2、几种常见的差分格式 以一维热传导方程为例：,（1）古典显式格式,（2）古典隐式格式,（3）六点格式（Crank-Nicolson），双层六点隐式格式 在x点和n+n/2时层，对t和x均采用中心差分,（4）Richardson格式，三层显式格式 在x点和n时层，对t和x均采用中心差分,（5）加权六点格式，隐式格式 在x点和n+n时层，01，对t和x均采用
8、中心差分,2.2.2线性方程组的数值解,有限差分法是工程中常用的一种离散技术，将计算域分成有限个网格，通过差分法求网格结点的微分方程的近似值，也称网格法。 将网格结点上的函数f(x,y,z,t)表示成 ，i,j,k分别表示x,y,z方向的坐标位置，n表示时间。,1、解线性方程组的两种方法: 直接法:通过有限步算术运算直接求出方程组的精确解,最常用的是消元结合代入的方法. 实际上除非是采用无穷位精度计算,一般都得不到精确解. 直接方法适用于解低阶稠密矩阵方程组.,迭代法 类似于方程求根的迭代法,用一个迭代过程逐步逼近方程组的解. 迭代有可能不收敛,或虽然收敛,但收敛速度慢. 迭代法适用于求解高阶
9、稀疏矩阵方程组. 稀疏矩阵:矩阵非零元素较少,且在固定的位置上. 稀疏矩阵一般是人为构造的,例如36页三转角插值时方程组(8.12),(8.15)的系数矩阵.,Gauss消去法(第一次消元),考虑方程组A(1)x=b(1),第一次消元用第一个方程将后面方程的x1消去.,计算乘数,条件:a11(1)0,用-mi1乘以第一个方程加到第i个(i=1,n)方程上,则消去了第i个方程中的x1.,Gauss消去法(第一次消元),经过上述过程,得到方程组A(2)x=b(2),其中,Gauss消去法(第k次消元),假设已完成k-1次消元,得到方程组A(k)x=b(k).,第k次消元的目的是将akk(k) (称
10、为主元)下面的元素变为0.,Gauss消去法(第k次消元),对A(k)右下角的矩阵,计算乘数,条件:akk(k)0,用-mik乘以第k个方程加到第i个(i=k+1,n)方程上,则消去了第i个方程中的xk,得到方程组A(k+1)x=b(k+1).,Gauss消去法(第k次消元),第一步消元的计算公式,类似可以得到第k步消元的计算公式,Gauss消去法,消去法完成后最终得到与原方程组等价的三角形方程组A(n)x=b(n).,一共需进行 ? 步,n-1,Gauss消去法(算法),追赶法求解三对角方程组,上面的方程组可以利用追赶法求解(P185).,对于下面形式的方程组,将系数矩阵进行三角分解,比较两
11、边对应元素可以得到,因此有,又,因此所有bi的可递推求出,进一步可求出ai,ri.,在得到系数矩阵的分解后,原方程组转化为 LUx=f.,先求解Ly=f,显然有y1=f1/a1, yi=(fi-riyi-1)/ai=(fi-aiyi-1)/(bi-aibi-1)(i=2,n),再求解Ux=y,显然有xn=yn, xi=yi-bixi+1(i=n-1,1),迭代法,在处理一元方程f(x)=0时,我们将其转化为x=j(x)的形式,然后用不动点迭代的方法进行求解. 对于线性方程组Ax=b,我们也可以将其转化为类似的形式: x=Bx+f, 任取初始向量x(0),令x(k+1)=Bx(k)+f(k=0,
12、1,),则得到一个向量的序列x(k). 若该序列收敛于向量x*,对x(k+1)=Bx(k)+f 两边取极限得到x*=Bx*+f,即x*是方程组的解.,Jacobi迭代法与Gauss-Seidel迭代法,对于方程组,我们将其改写为,Jacobi迭代法,写成矩阵的形式为x=B0x+f,其中,Jacobi迭代法,利用x(k+1)=Bx(k)+f 进行迭代,得到结果如下,Jacobi迭代法,从上表可以看出,迭代序列逐步逼近方程组的精确解(3,2,1)T. 注:在迭代中,我们不可能得到x(k)和精确解之间的误差,一般我们用|x(k)-x(k-1)|(通常用无穷范数)的值来判断是否终止迭代. 在上面的例子
13、中,我们将第i个方程变形为左边是xi,右边是其它分量和常数的线性组合,然后进行迭代,这一方法称位Jacobi迭代.,Jacobi迭代法,一般的,对于方程组Ax=b,设A非奇异且aii0(i=1,2,n),将A改写为A=D L U,其中,Jacobi迭代法,将方程组改写为 Dx=(L+U)x+b x=D1(L+U)x+D1b 令B0=D1(L+U)(称位Jacobi迭代矩阵),f=D1b,上式简记为x=B0x+f. 我们得到Jacobi迭代公式 x(k+1)=B0x(k)+f.,写成分量的形式为,Gauss-Seidel迭代法,在前面的例子中,我们计算x1(k+1),用的是第k步的x2,x3;,
14、计算x2(k+1),用的是第k步的x1,x3,我们有理由认为已经计算出的第k+1步的x1比第k步的“好”.因此,我们应该用第k+1步的x1和第k步的x3来计算x2.,类似地,我们也应该用新信息计算x3.,Gauss-Seidel迭代法,我们可以将上面一般的Jacobi迭代公式改写为,这一迭代方法称为Gauss-Seidel迭代.,Gauss-Seidel迭代法(算例),其Gauss-Seidel迭代公式为,对于方程组,Gauss-Seidel迭代法(算例),同样取x(0)=(0,0,0)T,迭代结果如下,超松弛迭代(SOR)方法,沿着从xi(k)到xi (k+1) (G)的方向再向前走,就得到
15、超松弛迭代(SOR)方法.,假设已知第k步的迭代向量x(k)以及第k+1步迭代向量x(k+1)的前i1个分量已知,Gauss-Seidel迭代法取,超松弛迭代方法,我们定义新的xi(k+1)为xi(k)与 的加权平均.,在w=1时,上述方法就是Gauss-Seidel方法,w1时称为超松弛法(有时不管w的范围,统称为超松弛方法).,超松弛迭代方法(算例),对于方程组,松弛方法迭代格式为,超松弛迭代方法(算例),取x(0)=0,w=1.3,终止准则为|x(k)x(k1)|10-5.,超松弛迭代方法(算例),我们来观察松弛因子w对收敛速度的影响.,步数表示|x(k)x(k1)|10-5时的迭代步数
16、,w=2.0时,500步以内不收敛.,超松弛迭代方法(矩阵表示),超松弛迭代格式可以写为,用矩阵可以表示为,第三章 二维水动力数值模拟,一、二维水动力数值模拟系统的分类 1、按差分网格分： 三角形、正方形、矩形、四边形、多边形、曲线坐标网格以及各种形状网格的组合 2、按计算方法分： 显式法、隐式法、显隐混合法 3、按模拟格式分： 三角元法、ADI法、破开算子法、单元体积法、MADI法、准分析法、贴体坐标法。,二、平面二维水动力数学模型的一般形式,定解条件,三、ADI法,1、网格 正方形或矩形，变量u,v,分别交错布置于网格中心和两侧。 2、ADI基本思想 （1）分步 （2）交错显隐,水位,、水深,u,v,3、差分格式 X向运动方程在（i+1/2,j)点离散,
