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

1ADI方法数值离散海洋原始方程组并在理想湖泊风生环流中的应用数学系薛鹏飞B00111624指导导师：朱建荣摘要：随着海洋科学的发展，对海洋现象的研究越来越倚重于定量和预测。计算机速度，容量和计算方法的飞速发展，使得海洋数值模式的应用越来越广泛，海洋数值模式在定量和预测海洋动力过程的研究和应用中已起着不可替代的作用。本文首先应用ADI方法对海洋原始方程组数值离散，再在理想湖泊的假设前提下，结合物理海洋学，有限差分法，利用MATLAB数学软件，进行三维数值模拟，以期对海洋数值模式做出最基础的实现，并从动力机制上分析，验证其模拟结果的近似准确性。Abstract:Withthescientificdevelopmentofocean,thestudyonmarinephenomenonreliesonmoreandmoretherationandpredicts.Thedevelopmentatfullspeedofcapacityandcomputingtechnologyandthespeedofthecomputer,makingtheapplicationofthemarinenumericalmodelmoreandmoreextensive,andmarinenumericalmodelwithpredictmarinepowerresearchofcourseandalreadyplayanirreplaceableroleofusinginration.ThistextuseADItodispersedequationgroupatfirst,onthepremiseofassumptionoftheideallake,combiningphysicaloceangraphy,finitedifference,MATLAB,carryonthreedimensionnumericalsimulation,expecttomakethemostbasicrealizationtothemarinenumericalmodel,analyse,proveitssimulationaccuracyofresultfrommotiveforcemechanism关键词：海洋控制方程组，f平面近似，三维数值模拟，理想湖泊，半动量格式，ADI方法，MATLABKeywords:Controltheequationgroupintheocean,f-levelapproximate,threedimensionvaluesimulation,ideallake,halfamomentumform,ADImethod,MATLAB一、海洋运动控制方程组海洋运动可以通过求解一组数学方程来描写。这些方程包括（1）运动方程，（2）连续方程以及（3）热量和盐量守恒方程。可利用质量守恒定律和牛顿第二定律导出这些方程。因本文为了进行简单的三维数值模拟，设计了一个理想湖泊，形状为一个长方体，四周都是固边界，底形无起伏，湖泊里的水盐度为0，纯净，无泥沙等杂质，为均质不可压流体，即密度为常数，所以只需要考虑动力学因素。控制方程组由动量方程，连续方程组成：)()()(1zuKzyuAyxuAxxPfvzuwyuvxuutumyx)()()(1zvKzyvAyxvAxyPfuzvwyvvxvutvmyx2gzP0zwyvxu由于密度取为常数,再由连续方程,利用边界条件，可以得到如下一组方程：)()()(zuKzyuAyxuAxxgfvzuwyuvxuutumyx)()()(zvKzyvAyxvAxygfuzvwyvvxvutvmyx0tvdzyudzxhh0|zhzhzxvdzyudzxw其中wvu,,分别为x,y,z方向的速度，myxKAA,,,分别为水位波动，水平湍流系数和垂向湍流系数。上述海洋方程组为时间，空间的变量，作为数学物理的适定问题，还必须给出初始条件和边界条件。初始条件：因海洋的动力（流场和水位）过程调整较快，初值取为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动力学边界条件：zatzvKzuKyzmxzm|,|（湖面）x，y分别为风应力矢量在x,y方向的分量。海底边界条件：运动学边界条件：w（x,y,-h,t）=0，动力学边界条件：3hzatzvKzuKbyhzmbxhzm|,|（湖底）bx,by分别为底摩擦力矢量在x,y方向的分量。岸边界条件：0tv，0nv。t表示岸边界的切向，n表示法向。2．方程离散及求解思路接下来对方程组进行离散，变量空间配置采用ArakawaC格式，解释如下：取一个小立方体，在其中心求水位，在其前后左右四个面的中心求X，Y方向速度U，V，上下两个面中心求Z方向速度W。平面图如下所示：其中r表示水位，u表示x方向速度，v表示y方向速度。整个模式数值计算方法的基本框架式欧拉前差，空间中央差，水位方程（连续方程）的求解采用隐式，整理成三对角方程组，矩阵形式后由MATLAB求解。模式中的动量方程分三步做时间积分。第一步：fvzuuwyuuvxuuutuukjikjikjikjikjikjikjikjikjinkjinkji,,,1,,,,,,1,,,,,,,1,,,,1,,12,1,,,,1,2,,1,,,,122yuuuAxuuuAkjikjikjiykjikjikjixfuzvvwyvvvxvvutvvkjikjikjikjikjikjikjikjikjinkjinkji,,,1,,,,,,1,,,,,,,1,,,,1,,12,1,,,,1,2,,1,,,,122yvvvAxvvvAkjikjikjiykjikjikjix在第一步中，仅考虑非线性项，科氏力项，水平涡动粘滞项，所有这些项均作显示处理。第二步：xgtuujinjinn,1,111124ygtvvjinjinn,11,1112在第二步中，仅考虑外重力波产生的水位梯度力，这个快过程作隐式处理，把它代入连续方程，用ADI方法求解水位分布，消除了CFL判据的严格限制，提高了计算效率。第三步：21,,,,1,,11223zuuuKtuukjikjikjimnn21,,,,1,,11223zvvvKtvvkjikjikjimnn注意：在上下边界处，zuKm用x代替，zvKm用y代替。在第三步中，仅考虑垂向涡动粘滞项，并作隐式处理，这样在垂向可大大提高分辨率。由第三步可直接得出三维速度场的分布。循环计算，可得到每一时间层的流场和水位。下面，在方程的某些项作一个说明：2.1科氏力项：由前面介绍知道，科氏参数f为与纬度有关的一个变数，考虑到所取湖泊尺度不是很大，将科氏参数取为一个定值，即取定纬度为固定纬度，这就是通常的f平面近似。2.2下面介绍一下海表面风应力的计算（上面讨论假设是已知的，但在实际计算中必须给出实际表达式）。adaVWC，zzVkmo处a为空气密度，1.2×10-3kgm-3；dC为大气对海面的拖曳系数；aV为风矢量，jViUVaaa；W为风速，22aaVUW。关于Cd给出，有不同的表达式。LargeandPond(1981)给出为Cd=1.2×10-3W11m/sCd=(0.49+0.065W)×10-3W≥11m/s在台风、热带气旋、飓风情况下，一般采用表达式（Miller,1964;Smiter,1980）：Cd=(0.73+0.069W)×10-3在目前情况下，Cd应取随风速的函数，不取为常数。风大浪高，可提升风对海水的作用2.3下面介绍海底摩擦力b的计算adbVWC，5底应力拖曳系数由近海底abz处的流速呈对数分布计算，]0025.0,)ln(/max[202zzkCabd，其中k为卡门常数，abz一般取为垂向分层底层厚度的一半，0z为海底粗糙度，一般取为0.001－0.002m。2.4垂向湍流扩散系数变化较大，尤其是在河口浅海，它应是时间，空间的函数，本模型中，为了简单起见，还是取为了常数。下面我们对ADI（AlternatingDirectionImplicit）格式的稳定性作一简单的分析：交叉方向隐式格式是对x,y方向交叉使用隐式格式（同时也交叉使用了显示格式），使得求解过程简化。现以下面方程为例，在概念，方法上作一介绍。)(22222yuxutu写成差分形式：nnunuDuDtuu2211212/，（2.1）122112112/nnnnuDuDtuu（2.2）其中2,1DD分别是222xu，222xu的差分算子。ADI思想是在一个时间步长t内，分成两步。第一步在x方向隐式，y方向显式求解（2.1）式，第二步在第一步的基础上，x方向显式，y方向隐式求解（2.2）式。这就是交叉方向隐式的意思，把一个二维线性代数方程组问题转化为一维线性代数方程组问题，求解有多种方法，如追赶法。再考虑其稳定性：)(ykxujinneBU（2.3）将（2.3）式带入（2.1）式，)]2cos2()2cos2([222222121yyaBxxaBtBBnnnn(2.4)记，2sin22221xuxtaa，2sin22222xyta，改写式（2.4）6212121nnnnBBBB)1()1(2121nnBB122111nnBB，将式（2.3）代入式（2.2），重复以上步骤可得2121111nnBB，这样，加辐因子AF＝)1)(1()1)(1(12121nnBB由此我们可以得出对任意时间步长t，|AF|1,即ADI方法是无条件绝对稳定的。所以ADI方法无论是在计算精度（证明从略），计算量，计算稳定性都是相当令人满意的，在海洋模式中已得到了广泛的应用。在模式的第一步中，我们使用了显式求解，而且所用到的方程组是复杂的非线性偏微分方程，对相应线性化方程的稳定性分析只能给出非线性差分方程计算稳定性的必要条件，即使满足这个条件，也可能会产生虚假的能量转移，扩大，甚至产生计算湍流，而造成计算不稳定。因非线性作用产生的不稳定称为非线性不稳定。非线性不稳定常常是突变的，在满足线性稳定性判据的条件下，非线性差分方程可以稳定的计算一段时间，随后计算结果突然迅速增长而产生不稳定，这种不稳定不能用缩短时间步长来克服。克服非线性不稳定的方法：（1）进行空间和时间平滑，滤去短波分量。（2）构造守恒的差分格式，如能量，动量，涡度拟能守恒。所以，为了提高第一步中（仅考虑非线性项，科氏力项，水平涡动粘滞项，所有这些项均作显示处理）的稳定性，采用了著名的半动量格式。该数值差分格式为：yyyxxxvuv（2.48）它即为处理非线性项的一，二次守恒格式。一次守恒对正负数值变化很大是，总的效应可以体现不出来，而二次守恒对每点都不能变化很大，因此二次守恒更能保证计算的稳