《数值模拟及数值试验-华东师范大学数学系.pdf》由会员分享，可在线阅读，更多相关《数值模拟及数值试验-华东师范大学数学系.pdf（19页珍藏版）》请在三一文库上搜索。
1、1 ADI 方法数值离散海洋原始方程组并在理想湖泊风生环流中的应用 数学系薛鹏飞B00111624 指导导师：朱建荣 摘要： 随着海洋科学的发展，对海洋现象的研究越来越倚重于定量和预测。计算机速度，容 量和计算方法的飞速发展，使得海洋数值模式的应用越来越广泛，海洋数值模式在定量和预 测海洋动力过程的研究和应用中已起着不可替代的作用。本文首先应用ADI 方法对海洋原 始方程组数值离散，再在理想湖泊的假设前提下，结合物理海洋学，有限差分法，利用 MATLAB数学软件，进行三维数值模拟，以期对海洋数值模式做出最基础的实现，并从动 力机制上分析，验证其模拟结果的近似准确性。 Abstract: Wit
2、h the scientific development of ocean, the study on marine phenomenon relies on more and more the ration and predicts. The development at full speed of capacity and computing technology and the speed of the computer, making the application of the marine numerical model more and more extensive, and m
3、arine numerical model with predict marine power research of course and already play an irreplaceable role of using in ration. This text use ADI to dispersed equation group at first, on the premise of assumption of the ideal lake, combining physical oceangraphy, finite difference , MATLAB ,carry on t
4、hree dimension numerical simulation, expect to make the most basic realization to the marine numerical model, analyse , prove its simulation accuracy of result from motive force mechanism 关键词：海洋控制方程组，f 平面近似，三维数值模拟，理想湖泊，半动量格式，ADI 方法， MATLAB Key words:Control the equation group in the ocean , f-level
5、approximate , three dimension value simulation, ideal lake, half a momentum form , ADI method , MATLAB 一、海洋运动控制方程组 海洋运动可以通过求解一组数学方程来描写。这些方程包括（1）运动方程， （2）连续 方程以及（ 3）热量和盐量守恒方程。可利用质量守恒定律和牛顿第二定律导出这些方程。 因本文为了进行简单的三维数值模拟，设计了一个理想湖泊，形状为一个长方体，四周 都是固边界， 底形无起伏， 湖泊里的水盐度为0，纯净， 无泥沙等杂质， 为均质不可压流体， 即密度为常数，所以只需要考虑动力学
6、因素。 控制方程组由动量方程，连续方程组成： )()()( 1 z u K zy u A yx u A xx P fv z u w y u v x u u t u myx )()()( 1 z v K zy v A yx v A xy P fu z v w y v v x v u t v myx 2 g z P 0 z w y v x u 由于密度取为常数,再由连续方程,利用边界条件，可以得到如下一组方程： )()()( z u K zy u A yx u A xx gfv z u w y u v x u u t u myx )()()( z v K zy v A yx v A xy gfu
7、z v w y v v x v u t v myx 0 t vdz y udz x hh 0| z h z h zx vdz y udz x w 其中wvu,分别为x,y,z 方向的速度， myx KAA,分别为水位波动，水平湍流系数和 垂向湍流系数。 上述海洋方程组为时间，空间的变量， 作为数学物理的适定问题，还必须给出初始条件和边 界条件。 初始条件： 因海洋的动力（流场和水位）过程调整较快，初值取为0 u(x,y,z,0)=0, v(x,y,z,0)=0, w(x,y,z,0)=0, (x,y,0)=0 。 边界条件： 海洋模式的边界包括垂向的海表面和海底，水平的固边界（岸界）和开边界（
8、水界），根 据本模式的简单假设，边界条件取为： 海表面边界条件： 运动学边界条件： w(x,y,0,t)=0 动力学边界条件： zat z v K z u K yzmxzm |,|（湖面） x ， y 分别为风应力矢量在x,y 方向的分量。 海底边界条件： 运动学边界条件： w（x,y,-h,t） =0， 动力学边界条件： 3 hzat z v K z u K byhzmbxhzm |,|（湖底） bx,by分别为底摩擦力矢量在 x,y 方向的分量。 岸边界条件： 0 t v， 0 n v。 t 表示岸边界的切向，n 表示法向。 2方程离散及求解思路 接下来对方程组进行离散，变量空间配置采用A
9、rakawa C 格式， 解释如下： 取一个小立方 体，在其中心求水位，在其前后左右四个面的中心求X，Y 方向速度U，V，上下两个面中 心求 Z 方向速度W。平面图如下所示： 其中 r 表示水位， u 表示 x 方向速度， v 表示 y 方向速度。 整个模式数值计算方法的基本框架式欧拉前差，空间中央差，水位方程 （连续方程）的求 解采用隐式，整理成三对角方程组，矩阵形式后由MATLAB求解。模式中的动量方程分三 步做时间积分。第一步： fv z uu w y uu v x uu u t uu kjikji kji kjikji kji kjikji kji n kji n kji, 1, ,
10、, 1, , , 1 , , 1 , 1 2 , 1,1, 2 ,1, 1 22 y uuu A x uuu A kjikjikji y kjikjikji x fu z vv w y vv v x vv u t vv kjikji kji kjikji kji kjikji kji n kji n kji,1, , , 1, , ,1 , , 1 , 1 2 , 1, 1, 2 ,1, 1 22 y vvv A x vvv A kjikjikji y kjikjikji x 在第一步中，仅考虑非线性项，科氏力项，水平涡动粘滞项，所有这些项均作显示处理。 第二步： x g t uuji n j
11、i nn , 1 , 1 11 12 4 y g t vvji n ji nn , 1 1, 11 12 在第二步中， 仅考虑外重力波产生的水位梯度力，这个快过程作隐式处理，把它代入连续 方程，用 ADI 方法求解水位分布，消除了CFL 判据的严格限制，提高了计算效率。 第三步： 2 1,1, 11 2 23 z uuu K t uukjikjikji m nn 2 1,1, 11 2 23 z vvv K t vvkjikjikji m nn 注意：在上下边界处， z u Km 用 x代替， z v Km 用 y 代替。 在第三步中，仅考虑垂向涡动粘滞项，并作隐式处理，这样在垂向可大大提高分
12、辨率。由 第三步可直接得出三维速度场的分布。循环计算，可得到每一时间层的流场和水位。 下面，在方程的某些项作一个说明： 2.1 科氏力项： 由前面介绍知道，科氏参数f为与纬度有关的一个变数，考虑到所取湖泊尺度不是很大， 将科氏参数取为一个定值，即取定纬度为固定纬度，这就是通常的f平面近似。 2.2 下面介绍一下海表面风应力 的计算（上面讨论假设是已知的，但在实际计算中必须给 出实际表达式） 。 ada VWC，z z V km o 处 a 为空气密度，1.2 10 -3 kgm -3 ； d C为大气对海面的拖曳系数； a V 为风矢量， jViUV aaa ；W为风速， 22 aa VUW。
13、关于 Cd给出，有不同的表达式。Large and Pond (1981) 给出为 Cd = 1.2 10 -3 W11 m/s Cd =(0.49+0.065W) 10 -3 W 11 m/s 在台风、热带气旋、飓风情况下，一般采用表达式（Miller, 1964; Smiter, 1980） ： Cd =(0.73+0.069W) 10 -3 在目前情况下，Cd应取随风速的函数，不取为常数。风大浪高，可提升风对海 水的作用 2.3 下面介绍海底摩擦力 b的计算 adb VWC， 5 底应力拖曳系数由近海底 ab z处的流速呈对数分布计算， 0025.0 ,)ln(/max 2 0 2 z
14、z kC ab d ， 其中k为卡门常数， ab z一般取为垂向分层底层厚度的一半， 0 z为海底粗糙度，一般取为 0.0010.002m。 2.4 垂向湍流扩散系数变化较大，尤其是在河口浅海，它应是时间，空间的函数，本模型 中，为了简单起见，还是取为了常数。 下面我们对ADI （Alternating Direction Implicit）格式的稳定性作一简单的分析： 交叉方向隐式格式是对x,y 方向交叉使用隐式格式（同时也交叉使用了显示格式），使得求 解过程简化。现以下面方程为例，在概念，方法上作一介绍。 )( 2 2 2 2 2 y u x u t u 写成差分形式： n n u n u
15、DuD t uu 2 2 1 1 2 1 2/ ， （2.1） 1 2 2 1 1 2 1 1 2/ n n n n uDuD t uu （2.2） 其中 2, 1D D分别是 2 2 2 x u ， 2 2 2 x u 的差分算子。 ADI 思想是在一个时间步长 t内，分成 两步。第一步在x 方向隐式， y 方向显式求解（2.1）式，第二步在第一步的基础上，x 方向 显式， y 方向隐式求解（2.2）式。这就是交叉方向隐式的意思，把一个二维线性代数方程组 问题转化为一维线性代数方程组问题，求解有多种方法，如追赶法。 再考虑其稳定性： )(ykxujinn eBU （2.3） 将（ 2.3）式
16、带入（ 2.1）式， )2cos2()2cos2( 2 2 2 2 2 2 1 2 1 y y a Bx x a B t BB n n n n (2.4) 记， 2 sin 22 2 2 1 xu x ta a， 2 sin 22 2 2 2 x y ta ， 改写式（ 2.4） 6 21 2 1 2 1 n n n n BBBB )1()1 ( 21 2 1 n n BB 1 2 2 1 1 1 n n B B ， 将式（ 2.3）代入式（ 2.2） ，重复以上步骤可得 2 1 2 1 1 1 1 n n B B ， 这样，加辐因子 AF )1)(1( )1)(1( 12 12 1 n n
17、B B 由此我们可以得出对任意时间步长t，|AF|1,即 ADI 方法是无条件绝对稳定的。所以 ADI 方法无论是在计算精度（证明从略），计算量，计算稳定性都是相当令人满意的，在海洋模 式中已得到了广泛的应用。 在模式的第一步中，我们使用了显式求解，而且所用到的方程组是复杂的非线性偏微分 方程，对相应线性化方程的稳定性分析只能给出非线性差分方程计算稳定性的必要条件，即 使满足这个条件，也可能会产生虚假的能量转移，扩大， 甚至产生计算湍流，而造成计算不 稳定。 因非线性作用产生的不稳定称为非线性不稳定。非线性不稳定常常是突变的，在满足 线性稳定性判据的条件下，非线性差分方程可以稳定的计算一段时间
18、，随后计算结果突然迅 速增长而产生不稳定，这种不稳定不能用缩短时间步长来克服。 克服非线性不稳定的方法： （1）进行空间和时间平滑，滤去短波分量。 （2）构造守恒的差分格式，如能量，动量，涡度拟能守恒。 所以，为了提高第一步中（仅考虑非线性项，科氏力项，水平涡动粘滞项，所有这些项均作 显示处理）的稳定性，采用了著名的半动量格式。该数值差分格式为： y y y x x x vuv（2.48） 它即为处理非线性项的一，二次守恒格式。 一次守恒对正负数值变化很大是，总的效应可以 体现不出来，而二次守恒对每点都不能变化很大，因此二次守恒更能保证计算的稳定性。 二：程序设计方法及相关命令简介： 模式的计
19、算机实现选择了优秀的数学软件MATLAB ，首先采用均匀长方形网格分层，即 在 x,y,z 三个方向分别是等间距的，这样，就可以把所有网格点上的值与一个三维数组相对 应，我们知道，MATLAB的基本运算单位就是矩阵，所以基本的想法就是把所有网格点上 的赋值给一个长方体矩阵（这是MATLAB特有的），为了便于书写，取行增加的方向为x 方向，列增加的方向为y 方向，页增加的方向为z的负方向。 值得强调的一点是，由于采用的是Arakawa C 格式， 如果 x 方向分层数为a，y 方向分层 数为 b，垂向分层数为c，则 uu（a,b-1,c-1） ，v=v(a-1,b,c-1) ，w=w(a-1,b
20、-1,c) ，r=r(a-1,b-1), 再整个循环计算结束后，再通过求平均，在水位点，也即每个小立方体的中心求得相应的 u,v,w，分别将其赋值给矩阵uu=uu(a-1,b-1,c-1),vv=vv(a-1,b-1,c-1),ww=ww(a-1,b-1,c-1)，以便 7 于最后数据的可视化。 另外，我们知道，ADI 方法需要求解两个一维的线性方程组,通常是三对角线性方程组， 这里将其整理为矩阵表达形式，方程求解可以直接利用MATLAB内部多种函数命令直接求 解， （例如共轭梯度法，松弛迭代法等等） 最后当整个求解过程完毕以后，为了给出直观的理解，还需要进行数据可视化处理，也 利用 MATL
21、AB的画矢量图的命令quiver（用来画流场分布） ，以及画等值线的命令contour （用来画水位分布） ，下面给出所需命令的使用帮助。 在本模型中，理想湖泊的长（length） ，宽（ width） ，高（ hight） ，x 方向分层数a，y 方向 分层数 b，垂向分层数c，x 方向风速windx ，y 方向风速windy ，湖泊纬度，以及风吹时间 均为可自由调节量，以便于模拟不同风速风向下，不同大小的理想湖泊的流场和水位分布。 实际运行当中， 取了长宽均为1000 公里，深为 20 米的理想湖泊， 在 5 米每秒的定常北风 （这 里为模型里的x 方向， 即矩阵行增加的方向） ，循环计算
22、10800 步，即为恒定风吹30 个小时 以后的，流场，水位分布。 三计算可视化结果及其动力机制分析 分析：由于定常北风和科氏力的共同作用，使得表层流场在向南走的同时，有略微西偏 的趋势。 8 第二层流场主要是由于南岸水位堆积产生向北的正压压强梯度力，以及表层水流向南产 生的垂向湍流剪切力， 在这两个力的共同作用下产生了上述流场，并且由图我们也可以看出， 第二层流速相对较小。 9 由于风应力对第三层流场影响较小，而水位梯度力（向北）不变，所以产生了向北的回 流。 底层流场情况大致与第三层流场相似，只是由于底摩擦作用使得流速略小于第三层流速。 10 垂向切面图可以清楚的显示流场的连续性 在北风作
23、用下水体向南运动，沿南岸堆积，导致水位上升，科氏力的作用使西侧水位高 于东侧水位 11 结论：这个理想湖泊中，由于风应力（北风）的驱动，使得表层水流向南流动，同时在 科氏力的作用下，使得表层流场略西偏，由于南岸水位堆积产生向北的正压压强梯度力，以 及表层水流向南产生的垂向湍流剪切力，在这两个力的共同作用下产生第二层流场，到第三 层，随着深度的增加，垂向湍流剪切力减小，由于正压压强梯度力不随深度变化，使得合力 偏向北， 产生回流， 底层基本受力基本与第三层一致，只是由于底摩擦力使得流速略有减小， 同时由于补偿性质（连续性和质量守恒），北岸近岸产生上升流，南岸产生下降流。从水位 图整个水位分布呈现
24、南高北低，西高东低的情况也是符合动力机制的。 四：海洋数值模式简介： 随着海洋科学的发展，对海洋现象的研究越来越倚重于定量和预测，所以海洋数值模式 的应用越来越广泛，已在水动力， 污染物和泥沙的输运和生态等方面取得了一系列的重大研 究成果。 水动力是物质输运的基础，流场对温度，盐度，污染物，泥沙，营养盐，浮游植物和浮 游动物等的输运起着决定性的作用，它们之间相互作用和影响，如温度和盐度的变化改变密 度场， 从而影响流场， 也影响着浮游植物和浮游动物的生长；泥沙媳妇和释放污染物和营养 盐，又阻碍太阳透射，从而影响光合作用，又影响叶绿素的生成，再影响浮游植物，浮游植 物再影响浮游动物。因此，海洋现
25、象之间相互作用，相互关联。 目前水动力模式已较为成熟， 在此基础上联合污染物，泥沙，营养盐，浮游植物，浮游动物等方程，采用三维高分辨率网 格，设计质量， 能量和涡度拟能守恒的差分格式，建立一个海洋整体模式，是海洋模式发展 的趋势。 基于海洋原始动力方程组的海洋数值模式在定量，预测海洋现象上有着独特的优势。在 科学研究上， 为揭示海洋变化的内在原因，可不受物理模型相似性理论的限制，单独和组合 考虑各种动力因子，通过数值模拟揭示海洋变化的动力机制，如海洋环流，河口盐水入侵， 河口入海物质输运， 全球变化对河口海洋的影响等。在社会发展和经济建设上，如三峡工程， 南水北调对长江河口及黄海，东海水文，
26、环境和生态等的影响，长江口深水航道工程对流场， 分流比， 泥沙输运的影响等，数值模式可在工程开工前尚无实测资料的情况下，定量预测它 们对河口海洋的影响， 减少负面影响和工程的失误。河口盐水入侵研究涉及水库建设的选址， 容量和水资源等问题。在军事上，舰艇的登陆需要预测波浪，潮滩的漫滩情况，潮流等，潜 艇航行需要预测密度跃层的分布和变化，以避免跟踪和窃听，涉及国家和军事的安全。因此 发展海洋数值模式，对社会发展，经济建设和军事活动等具有重要的意义。 目前国际上先进的和使用较为广泛的河口海洋数值模式有美国普林斯顿大学的POM （Princeton Ocean Model ）和 ECOM （Estra
27、ry，Cvoast and Ovcean Mvodel ） ，佛罗里达大学 的等密度面模式，荷兰的Delft 模式，德国的汉堡模式和丹麦模式。我国缺少自己研制的在 国际上有影响力的模式，使用较多的是POM 和 ECOM 模式，但一般的使用着缺乏了解模 式的数值计算方法，如遇到实际问题，很难找出原因和改进模式。这就是我国在海洋数值模 式上与国际先进水平的主要差距。 参考文献： 1 朱建荣2003 海洋数值计算方法和数值模式海洋出版社 2 叶安乐李凤歧1991 物理海洋学青岛海洋大学出版社 3 李荣华1995 微分方程数值解法高等教育出版社 4 张宜华1998 精通 MATLAB5 清华大学出版社
28、 附录：模型源代码 close,clear all 湖泊长宽高赋值 length=100000; 12 width=100000; hight=20; x,y,z 方向分层数 a=25; b=25; c=5; x,y,z 方向空间步长 dtax=length/(a-1); dtay=width/(b-1); dtaz=hight/(c-1); 科氏力计算 cor=30; f=2*7.292*10(-5)*sin(cor*pi*2/360); 水平，垂向涡动系数 km=0.001; kn=0.0005; 风速赋值 windx=5; windy=0; 风应力大小 windxy=sqrt(windx
29、2+windy2); if windxy11 cd=1.2*10(-3); else cd=(0.49+0.065*windxy)*10(-3); end 底摩擦力大小 z0=0.001; kk=(0.41/log(dtaz/2/z0)2; cb=max(0.0025,kk); tausx=1.2*10(-3)*cd*windxy*windx; tausy=1.2*10(-3)*cd*windxy*windy; 时间步长 dtat=10; 循环次数 n=10800; 所需常量赋值 g=9.8; u=zeros(a+2,b+1,c+1);v=zeros(a+1,b+2,c+1);w=zeros(
30、a+1,b+1,c+2);r=zeros(a-1,b-1);uu=zero s(a-1,b-1,c-1);vv=zeros(a-1,b-1,c-1);taubx=r;tauby=r;uu=zeros(a-1,b-1,c-1);vv=zeros(a-1,b-1,c-1); u3=u; vvu=u3; wwu=vvu; 13 v3=v; uuv=v3; wwv=uuv; vvw=w; uuw=w; q1=(dtat/4/dtax); q2=(dtat/4/dtay); q3=(dtat/4/dtaz); e1=(dtaz*dtat2*g*(c-1)/4/dtax2); e2=(dtaz*dtat2
31、*g*(c-1)/4/dtay2); ee1=(dtat*dtaz/dtax/2); ee2=(dtat*dtaz/dtay/2); e3=km*dtat/(dtaz2); 开始循环计算 for tt=1:n 显示循环次数 if mod(tt,60)=0 times=tt end 第一步，显式部分计算 dx=uu(:,:,(c-1); dy=vv(:,:,(c-1); taubx=cb*dx.*sqrt(dx.2+dy.2); tauby=cb*dy.*sqrt(dx.2+dy.2); for k=2:c for j=2:b for i=3:a vvu(i,j,k)=(v(i,j,k)+v(i
32、,j+1,k)+v(i-1,j,k)+v(i-1,j+1,k)/4; wwu(i,j,k)=(w(i-1,j,k)+w(i-1,j,k+1)+w(i,j,k)+w(i,j,k+1)/4; end end end for k=2:c for j=2:b for i=3:a sx11=u(i+1,j,k)+u(i,j,k); sx12=u(i+1,j,k)-u(i,j,k); sx21=u(i,j,k)+u(i-1,j,k); sx22=u(i,j,k)-u(i-1,j,k); sy11=vvu(i,j+1,k)+vvu(i,j,k); 14 sy12=u(i,j+1,k)-u(i,j,k); s
33、y21=vvu(i,j,k)+vvu(i,j-1,k); sy22=u(i,j,k)-u(i,j-1,k); sz11=wwu(i,j,k+1)+wwu(i,j,k); sz12=u(i,j,k+1)-u(i,j,k); sz21=wwu(i,j,k)+wwu(i,j,k-1); sz22=u(i,j,k)-u(i,j,k-1); u3(i,j,k)=u(i,j,k)-q1*(sx11*sx12+sx21*sx22)-q2*(sy11*sy12+sy21*sy22)-q3*(sz11*sz12+sz21*s z22)+f*vvu(i,j,k)+(kn/dtax2)*(u(i+1,j,k)+u(
34、i-1,j,k)-2*u(i,j,k)+(kn/dtay2)*(u(i,j+1,k)+u(i,j-1,k)- 2*u(i,j,k); end end end for k=2:c for i=2:a for j=3:b uuv(i,j,k)=(u(i,j,k)+u(i+1,j,k)+u(i+1,j-1,k)+u(i,j-1,k)/4; wwv(i,j,k)=(w(i,j,k)+w(i,j,k+1)+w(i,j-1,k)+w(i,j-1,k+1)/4; end end end for k=2:c for i=2:a for j=3:b sx11=uuv(i+1,j,k)+uuv(i,j,k); s
35、x12=v(i+1,j,k)-v(i,j,k); sx21=uuv(i,j,k)+uuv(i-1,j,k); sx22=v(i,j,k)-v(i-1,j,k); sy11=v(i,j+1,k)+v(i,j,k); sy12=v(i,j+1,k)-v(i,j,k); sy21=v(i,j,k)+v(i,j-1,k); sy22=v(i,j,k)-v(i,j-1,k); sz11=wwv(i,j,k+1)+wwv(i,j,k); sz12=v(i,j,k+1)-v(i,j,k); sz21=wwv(i,j,k)+wwv(i,j,k-1); sz22=v(i,j,k)-v(i,j,k-1); v3(
36、i,j,k)=v(i,j,k)-q1*(sx11*sx12+sx21*sx22)-q2*(sy11*sy12+sy21*sy22)-q3*(sz11*sz12+sz21*s z22)-f*uuv(i,j,k)+(kn/dtax2)*(v(i+1,j,k)+v(i-1,j,k)-2*v(i,j,k)+(kn/dtay2)*(v(i,j+1,k)+v(i,j-1,k)- 15 2*v(i,j,k); end end end for i=1:a for j=1:b-1 for k=1:c-1 u1(i,j,k)=u3(i+1,j+1,k+1); end end end for j=1:b for i
37、=1:a-1 for k=1:c-1 v1(i,j,k)=v3(i+1,j+1,k+1); end end end ，第二步，ADI 方法计算 mm,nn,pp=size(r); gg=ones(mm,3); gg(:,1)=-e1; gg(:,2)=2*e1+1; gg(:,3)=-e1; aa=spdiags(gg,-1 0 1,mm,nn);aa(1)=aa(1)-e1;aa(mm*nn)=aa(mm*nn)-e1; for j=1:b-1 for i=1:a-1 sv=(dtaz/dtay)*(sum(v1(i,j+1,:)-sum(v1(i,j,:); bb(i)=r(i,j)-ee
38、1*(sum(u1(i+1,j,:)-sum(u1(i,j,:)-(dtat/2)*(sv+r(i,j)*(u1(i+1,j,1)-u1(i,j,1)/dtax+(v 1(i,j+1,1)-v1(i,j,1)/dtay); end r(:,j)=aa(bb); %r(:,j)=bicg(aa,bb); end for k=1:c-1 for j=1:b-1 for i=2:a-1 u2(i,j,k)=u1(i,j,k)-(g*dtat/2/dtax)*(r(i,j)-r(i-1,j); end end end 16 mm,nn,pp=size(r); gg=ones(mm,3); gg(:,1
39、)=-e2; gg(:,2)=2*e2+1; gg(:,3)=-e2; aa=spdiags(gg,-1 0 1,mm,nn);aa(1)=aa(1)-e2;aa(mm*nn)=aa(mm*nn)-e2; for i=1:a-1 for j=1:b-1 su=(dtaz/dtax)*(sum(u2(i+1,j,:)-sum(u2(i,j,:); bb(j)=r(i,j)-ee2*(sum(v1(i,j+1,:)-sum(v1(i,j,:)-(dtat/2)*(su+r(i,j)*(u2(i+1,j,1)-u2(i,j,1)/dtax+(v 1(i,j+1,1)-v1(i,j,1)/dtay);
40、 end r(i,:)=(aa(bb); %r(i,:)=bicg(aa,bb); end for k=1:c-1 for j=1:b-1 for i=2:a-1 u1(i,j,k)=u1(i,j,k)-(g*dtat/dtax)*(r(i,j)-r(i-1,j); end end end for k=1:c-1 for i=1:a-1 for j=2:b-1 v1(i,j,k)=v1(i,j,k)-(g*dtat/dtay)*(r(i,j)-r(i,j-1); end end end 第三步，垂向隐式计算 mm,nn,pp=size(r); gg=ones(mm,3); gg(:,1)=-e
41、3; gg(:,2)=2*e3+1; gg(:,3)=-e3; aa=spdiags(gg,-1 0 1,mm,nn);aa(1)=aa(1)-e3;aa(mm*nn)=aa(mm*nn)-e3; for i=2:a-1 for j=1:b-1 for k=1:c-1 bb(k)=u1(i,j,k); end bb(1)=bb(1)+(dtat*tausx/dtaz); 17 bb(c-1)=bb(c-1)-(dtat*taubx(i,j)/dtaz); qq=(aabb); %qq=bicg(aa,bb); for k=1:c-1 u(i+1,j+1,k+1)=qq(k); end end
42、end for j=2:b-1 for i=1:a-1 for k=1:c-1 bb(k)=v1(i,j,k); end bb(1)=bb(1)+(dtat*tausy/dtaz); bb(c-1)=bb(c-1)-(dtat*tauby(i,j)/dtaz); qq=aabb; %qq=bicg(aa,bb); for k=1:c-1 v(i+1,j+1,k+1)=qq(k); end end end 求最后三维速度 for i=2:a for j=2:b for k=3:c uuw(i,j,k)=(u(i+1,j,k-1)+u(i,j,k-1)+u(i+1,j,k)+u(i,j,k)/4;
43、 vvw(i,j,k)=(v(i,j+1,k)+v(i,j,k)+v(i,j+1,k-1)+v(i,j,k-1)/4; end end end for i=2:a for j=2:b for k=3:c w(i,j,k)=-dtaz*(sum(uuw(i+1,j,(k:c)-sum(uuw(i,j,(k:c)/dtax-dtaz*(sum(vvw(i,j+1,(k:c)-sum( vvw(i,j,(k:c)/dtay; end end end for i=1:a 18 for j=1:b-1 for k=1:c-1 u3(i+1,j+1,k+1)= u1(i,j,k); end end end
44、 for j=1:b for i=1:a-1 for k=1:c-1 v3(i+1,j+1,k+1)=v1(i,j,k); end end end uu=zeros(a-1,b-1,c-1);vv=uu;ww=vv; for k=1:c-1 for j=1:b-1 for i=1:a-1 uu(i,j,k)=(u(i+1,j+1,k+1)+u(i+2,j+1,k+1)/2; end end end for k=1:c-1 for i=1:a-1 for j=1:b-1 vv(i,j,k)=(v(i+1,j+2,k+1)+v(i+1,j+1,k+1)/2; end end end for i=1
45、:a-1 for j=1:b-1 for k=1:c-1 ww(i,j,k)=(w(i+1,j+1,k+2)+w(i+1,j+1,k+1)/2; end end end end 数据可视化处理 X=ones(size(r); Y=X; for i=1:size(Y,1) 19 Y(i,:)=Y(i,:)*(-i); end for j=1:size(X,1) X(:,j)=X(:,j)*j; end save wind1 for i=1:c-1 scale=(sum(sum(abs(uu(:,:,i)+sum(sum(abs(vv(:,:,i)/(sum(sum(abs(uu(:,:,1)+s
46、um(sum(abs (vv(:,:,1); figure(i),quiver(X,Y ,vv(:,:,i),-uu(:,:,i),scale),pause; end %figure(1),quiver(X,Y ,vv(:,:,1),-uu(:,:,1),pause; %figure(2),quiver(X,Y ,vv(:,:,2),-uu(:,:,2),pause; %figure(3),quiver(X,Y ,vv(:,:,3),-uu(:,:,3),pause; %figure(4),quiver(X,Y ,vv(:,:,4),-uu(:,:,4) figure(c),contour(X,Y ,r,10)