数值模拟及数值试验-华东师范大学数学系

ADI 方法数值离散海洋原始方程组并在理想湖泊风生环流中的应用数学系 薛鹏飞 B00111624 指导导师：朱建荣摘要：随着海洋科学的发展，对海洋现象的研究越来越倚重于定量和预测。

计算机速度，容量和计算方法的飞速发展，使得海洋数值模式的应用越来越广泛，海洋数值模式在定量和预测海洋动力过程的研究和应用中已起着不可替代的作用。

本文首先应用ADI 方法对海洋原始方程组数值离散，再在理想湖泊的假设前提下，结合物理海洋学，有限差分法，利用MATLAB 数学软件，进行三维数值模拟，以期对海洋数值模式做出最基础的实现，并从动力机制上分析，验证其模拟结果的近似准确性。

Abstract: With 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 marine 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 three 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 方法，MATLABKey words:Control the equation group in the ocean , f-level approximate , three dimension value simulation, ideal lake, half a momentum form , ADI method , MATLAB一、海洋运动控制方程组海洋运动可以通过求解一组数学方程来描写。

这些方程包括（1）运动方程，（2）连续方程以及（3）热量和盐量守恒方程。

可利用质量守恒定律和牛顿第二定律导出这些方程。

因本文为了进行简单的三维数值模拟，设计了一个理想湖泊，形状为一个长方体，四周都是固边界，底形无起伏，湖泊里的水盐度为0，纯净，无泥沙等杂质，为均质不可压流体，即密度为常数，所以只需要考虑动力学因素。

控制方程组由动量方程，连续方程组成：)()()(1z u K z y u A y x u A x x P fv z u w y u v x u u t u m y x ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂-=-∂∂+∂∂+∂+∂∂ρ )()()(1zv K z y v A y x v A x y P fu z v w y v v x v u t v m y x ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂-=+∂∂+∂∂+∂∂+∂∂ρg zPρ-=∂∂ 0=∂∂+∂∂+∂∂zw y v x u 由于密度取为常数,再由连续方程,利用边界条件，可以得到如下一组方程：)()()(z u K z y u A y x u A x x g fv z u w y u v x u u t u m y x ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂-=-∂∂+∂∂+∂+∂∂ς )()()(z v K z y v A y x v A x y g fu z v w y v v x v u t v m y x ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂-=+∂∂+∂∂+∂∂+∂∂ς 0=∂∂+∂∂+∂∂⎰⎰--t vdz y udz x h h ςςς 0|=∂∂+∂∂+⎰⎰--=zhz h z x vdz y udz x w 其中w v u ,,分别为x,y,z 方向的速度，m y x K A A ,,,ς分别为水位波动，水平湍流系数和垂向湍流系数。

上述海洋方程组为时间，空间的变量，作为数学物理的适定问题，还必须给出初始条件和边界条件。

初始条件：因海洋的动力（流场和水位）过程调整较快，初值取为0u(x,y,z,0)=0, v(x,y,z,0)=0, w(x,y,z,0)=0, ς(x,y,0)=0。

边界条件：海洋模式的边界包括垂向的海表面和海底，水平的固边界（岸界）和开边界（水界），根据本模式的简单假设，边界条件取为： 海表面边界条件： 运动学边界条件：w(x,y,0,t)=0动力学边界条件：ςτρτρςς==∂∂=∂∂==z at zv K z u K y z m x z m|,|（湖面）x τ，y τ分别为风应力矢量在x,y 方向的分量。

海底边界条件：运动学边界条件：w （x,y,-h,t ）=0，动力学边界条件：h z at zvK zuK by h z mbx h z m-==∂∂=∂∂-=-=τρτρ|,|（湖底）bx τ,by τ分别为底摩擦力矢量在x,y 方向的分量。

岸边界条件：0=t v ， 0=n v 。

t 表示岸边界的切向，n 表示法向。

2．方程离散及求解思路接下来对方程组进行离散，变量空间配置采用Arakawa C 格式，解释如下：取一个小立方体，在其中心求水位，在其前后左右四个面的中心求X ，Y 方向速度U ，V ，上下两个面中心求Z 方向速度W 。

平面图如下所示：其中r 表示水位，u 表示x 方向速度，v 表示y 方向速度。

整个模式数值计算方法的基本框架式欧拉前差，空间中央差，水位方程（连续方程）的求解采用隐式，整理成三对角方程组，矩阵形式后由MATLAB 求解。

模式中的动量方程分三步做时间积分。

第一步：fv zu u w yu u v xu u u tu u kj i k j i kj i kj i k j i kj i kj i k j i kj i nkj i n k j i +∆--∆--∆--=∆-++++,,,1,,,,,,1,,,,,,,1,,,,1,,12,1,,,,1,2,,1,,,,122yu u u A xu u u A kj i k j i k j i ykj i k j i k j i x∆+-+∆+-+-+-+fu zv v w yv v v xv v u tv v kj i k j i kj i kj i k j i kj i kj i k j i kj i n kj i n k j i +∆--∆--∆--=∆-++++,,,1,,,,,,1,,,,,,,1,,,,1,,12,1,,,,1,2,,1,,,,122yv v v A xv v v A kj i k j i k j i ykj i k j i k j i x∆+-+∆+-+-+-+在第一步中，仅考虑非线性项，科氏力项，水平涡动粘滞项，所有这些项均作显示处理。

第二步：xgt u u ji n j i n n ∆--=∆-++++,1,11112ςςygt v v ji n j i n n ∆--=∆-++++,11,1112ςς 在第二步中，仅考虑外重力波产生的水位梯度力，这个快过程作隐式处理，把它代入连续方程，用ADI 方法求解水位分布，消除了CFL 判据的严格限制，提高了计算效率。

第三步：21,,,,1,,11223z u u u K t u u k j i k j i k j i m n n ∆+-=∆--+++ 21,,,,1,,11223z v v v K t v v k j i k j i k j i m n n ∆+-=∆--+++注意：在上下边界处，zuK m∂∂用x τ代替，z v K m ∂∂用y τ代替。

在第三步中，仅考虑垂向涡动粘滞项，并作隐式处理，这样在垂向可大大提高分辨率。

由第三步可直接得出三维速度场的分布。

循环计算，可得到每一时间层的流场和水位。

下面，在方程的某些项作一个说明： 2.1科氏力项：由前面介绍知道，科氏参数f 为与纬度有关的一个变数，考虑到所取湖泊尺度不是很大，将科氏参数取为一个定值，即取定纬度为固定纬度，这就是通常的f 平面近似。

2.2下面介绍一下海表面风应力τ的计算（上面讨论假设是已知的，但在实际计算中必须给出实际表达式）。

a d a V W Cρτ=， ζτρ==∂∂z zV k mo处a ρ为空气密度，1.2×10-3kgm -3；d C 为大气对海面的拖曳系数；a V为风矢量，j V i U V a a a +=；W 为风速，22a aV U W +=。

关于C d 给出，有不同的表达式。

Large and Pond (1981)给出为C d = 1.2×10-3W<11 m/sC d =(0.49+0.065W)×10-3W ≥11 m/s在台风、热带气旋、飓风情况下，一般采用表达式（Miller, 1964; Smiter, 1980）：C d =(0.73+0.069W)×10-3在目前情况下，C d 应取随风速的函数，不取为常数。

风大浪高，可提升风对海水的作用2.3下面介绍海底摩擦力b τ的计算a db V W Cρτ=，底应力拖曳系数由近海底ab z 处的流速呈对数分布计算，]0025.0,)ln(/max[22z z k C ab d =， 其中k 为卡门常数，ab z 一般取为垂向分层底层厚度的一半，0z 为海底粗糙度，一般取为0.001－0.002m 。

2.4垂向湍流扩散系数变化较大，尤其是在河口浅海，它应是时间，空间的函数，本模型中，为了简单起见，还是取为了常数。