海洋数值模型的理论及应用.ppt

《海洋数值模型的理论及应用.ppt》由会员分享，可在线阅读，更多相关《海洋数值模型的理论及应用.ppt（58页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第一节海洋数值模型的发展概况 观测，试验 理论分析 数值模型海洋的研究手段.海洋数值模型发展的客观要求 海水运动方程推导过程中的动力学假设，比如海水静力学假设，都是建立在大部分空间尺度大于10km，时间尺度大于惯性周期（2/f）的情况下的 海水运动方程不存在精确解：（1）偏微分方程本身是非线性的，（2）将一个实际问题具体化时所需要的环境场如海底地形、海岸线几何形状、表面驱动力等都无法给出解析函数表达式，只能在离散化的空间和时间区域取得 因此，求解海水运动方程只能通过近似求解的办法来确定，而这些数值解又可能引入明显的近似误差.海洋数值模型发展的历史 海洋模型到按其水平网格的离散方式以及所使用的垂

2、向坐标系的不同大致经历了如下几个发展阶段 最早出现并且还在使用的海洋模型是Bryan等人开发的基于原始方程的低阶精度的有限差分模型，它在水深方向采用z坐标系 目前常见的HOPS（Harvard Ocean Prediction System），MOM（GFDL Modular Ocean Model），POP(Parallel Ocean Program)，NCOM（NCAR Community Ocean Model）模型在某种程度上都可以认为是该模型改进版.海洋数值模型发展的历史 上个世纪70年代，sigma 坐标系开始应用于海洋模型在水深方向，比如目前被广泛使用的POM（Princeto

3、n Ocean Model）、ECOM（Estuarine Coastal and Ocean Model）、ROMS（Regional Ocean Modeling System）模型都属于这种类型的模型 与传统的 z 坐标系相比，sigma坐标可以更好地贴合海底地形变化.海洋数值模型发展的历史 自上世纪90年代以来，非结构化网格技术开始在海洋模型中得到应用，相应的也出现了基于非结构化网格离散方式的有限元模型，如SEOM 模型（Spectral Finite Element Ocean Model），和有限体积模型，如FVCOM 虽然与传统的结构化网格相比,非结构化网格可以更好地拟合陆地边界

4、，但是代码实现上的困难以及计算稳定性的问题使其迄今还没有得到非常广泛的应用.新一代的海洋数值模型 新一代的海洋模型广泛采用随地坐标系（terrain-following coordinates），进而促进了有关时间步长，对流项和压力梯度项等数值算法的改进 进入新世纪以来，下一代的海洋数值动力模型正在紧锣密鼓的研制中，代表性的是TOMS(Terrain-following Ocean Modeling System)，它融合了目前最先进的物理知识、数值方法和数据同化技术 http:/www.aos.princeton.edu/WWWPUBLIC/htdocs.pom/TOMS.htm.海洋模型的

5、商业软件 就海流的仿真建模研究而言，现阶段国外已形成了不少具有代表性的、具有强大前后处理功能及核心仿真建模技术的软件或系统（商业软件）如加拿大Dalhousie University 研究的DALCOAST 河口海岸预报系统；丹麦的DHI Water&Environment 机构开发的的MIKE 系列软件系统；荷兰的WL Delft Hydraulics 开发的的Delft 3D 软件.第二节基于POM 模型的一般介绍.水动力模型基本特点n 垂直方向的 坐标变换n 水平网格采用正交曲线坐标和ArakawaC差分格式n 水平时间差分采用显格式，而垂直差分为隐格式n 自由表面可以模拟水位变化n 垂

6、直和水平方向的混合扩散分别采用2.5阶的Mellor-Yamada湍流闭合模式和Smagorinski模式n 内外模态分别处理速度较慢的内重力波和速度较快的外重力波以提高整个模式计算效率n 包含了海水的热动力过程.Sigma 坐标变换：hs=-1z=0z=H(x,y)s=0 V(i,j+1)V(i,j)U(i+1,j)U(i,j)(i,j)模型采用（a）垂向sigma坐标系和（b）水平Arakawa C网格(1).控制方程连续方程海水运动方程(2)(3)(4).温度方程湍流动能方程盐度方程湍流混合长度方程(5)(6)(7)(8).其中，是基于sigma 坐标系的垂向流速矢量，它与笛卡儿坐标系下

7、的垂向流速W 之间的关系可表示为另外，方程中 另外，方程中 分别表示海水的实际密度，分别表示海水的实际密度，Boussinesq Boussinesq 近似密度以及密度扰动 近似密度以及密度扰动 其中，是实际水深(9).MY 湍流闭合模型（profq）在湍流动能及混合长度方程中在湍流动能及混合长度方程中 分别表示湍流混合系数，热扩散系数和湍流动能的垂直扩散系数；墙近似函数，其中，=0.4 是von Karman 常数 湍流动能和混合长度方程组由下列等式关系闭合 其中，是稳定函数(10).在2.5阶MY湍流闭合模型中，可表示为其中，。在不稳定层次 条件下，的上限值为0.023，对应的 值分别为2

8、.0145和2.4401；在稳定层次 条件下，它下限值为-0.28，对应的 值分别为 0.0470和 0.0461。是2.5阶M-Y湍流闭合模式参数，由实验测得(11).MY 湍流模型的适用条件M-Y湍流模型已被广泛应用于浅海潮汐，风生表面混合层的模拟以及底边界层的研究；该模型在混合较弱的层化流体中对湍流混合系数计算效果不佳，应用于河口的可行性也有待于探讨；虽然该模型无论是在物理上还是在数学上都存在着明显的缺陷，但是目前依然得到广泛应用。只有深入了解海水的微细结构以及海洋湍流结构和性质后，才有可能构建出更为完善的湍流闭合模型.海水运动方程的水平扩散项 其中水平切应力 是水平涡粘性系数，由Sma

9、gorinsky 模式计算(12)(13).Smagorinsky 水平扩散模式 其中参数C 为无量纲值，其取值为0.1或0.2，如果计算网格划分得足够细的话，C 可以取值为0 由Smagorinsky模式计算的水平涡粘性系数随着网格精度的改善和速度梯度的减小而减小(14).温度盐度的水平扩散项 其中 是水平热扩散系数，一般地，是一个小值，取0.1或0.2，甚至在某些情况下可以取为0(15)(16).垂向边界条件1.海面(17).2.海底其中，分别是湍流闭合模型参数和摩擦速度；分别是海面和海底的摩擦系数；(18).模型的外模态 控制近岸环流的海水运动方程包含了运动速度较快的外重力波和速度较慢的

10、内重力波，因此进行模式分裂可以大大提高模型的计算效率 模型的外模态就是计算自由水位和垂向平均的流速变化，需要较小的时间步长 内模态主要模拟流速矢量、温度、盐度等三维结构变化，对时间步长的要求相对宽松 采用模式分裂技巧可以有效地减少因外模态计算所消耗的计算时间，从而提高计算效率.垂向积分的海水运动方程其中：连续方程海水运动方程(19)(20)(21).水平扩散项水平数值耗散项(22)(23).模型代码中的一些符号 im,jm,kb,imm1,jmm1,kbm1,mode,isplit dte,dti,days,umol z,zz,dz,dzz aam2d,art,aru,arv,dum,dvm,

11、fsm dx,dy,h,el,d,ua,va,ut,vt,et,cor swrad wusurf,wvsurf,wubot,wvbot,wtsurf,wssurf,rad u,v,t,s,rho,l,q2 km,kh,kq,aam,aah rmean,tclim,sclim.数值积分流程图.外模态和内模态的相互作用.内模态的计算方法三维变量的计算（以T 为例）可以分解为垂直扩散项的计算和水平对流扩散项的计算，以温度方程为例方程写为求解过程分为两步，第一步计算水平对流扩散项，采用中心差分格式，由advt 子程序实现(24)(25)水平对流扩散项垂直扩散项.第二步计算垂向扩散项，数值方法采用“蛙跳

12、”格式，由proft 子程序实现由于“蛙跳”格式在奇数时间步时的解与在偶数时间步时的解会发生偏离，因此每一时间步的计算结束后还需对上述解进行平滑处理，即 为一个小值，一般取0.05。处理后，(26)(27).计算网格的安排二维外模态网格分布三维内模态网格分布.水平对流项的计算 虽然模型采用有限差分计算格式，但是对于每个网格对流项的计算都是按照有限容积的方法进行处理，即温度对流过程可表示为 速度的对流过程与温度相似，可写为 其中，是由z坐标转化为坐标后产生的弯曲项(28)(29).垂直扩散项的计算 垂直扩散项的计算公式（第k层，1kkb-1）可以改写为(30)(31).展开后合并同类项得到 其中

13、(32)(33).现在假定温度解的形式为 由(34)得到的 代入到(32)就可以得到其中 的值由(33)求得，的值由上一层的 值确定(34)(35).当k=1时，即表层海水，海水温度主要取决于海表面温度通量。公式（34）的解可以近似写为 这样，表层海水的温度只能根据式（31）求解，即 上式可以进一步表示为(36)(37).式（37）还可以进一步写为 再与温度的通解（34）比较 可以得到(34)(38)(39).当短波辐射通量 的计算公式 r，ad1，ad2 均为光辐射常数，根据不同水质取值如下(40)ntp 1 2 3 4 5Jerlov type I Ia Ib II IIIr 0.58 0

14、.62 0.67 0.70 0.78ad1(m)0.35 0.60 1.0 1.5 1.4ad2(m)23.0 20.0 17.0 14.0 7.9来源：Jerlov，1976；Paulson and Simpson，1977.当k=kb-1时，即底层海水，假定海底的热通量为0，根据前面的推导方式同样可以得到底层海水温度的计算公式 对于盐度方程垂向扩散项的计算与温度方程相同，只是不考虑太阳短波辐射这一项(41).时间步长的CFL 限制条件 由于模型水平方向采用显格式，因此时间步长的选取必须要满足CFL稳定条件。对于外模态，时间步长的限制条件为 其中，可能预见到的最大流速 对于内模态，时间步长的

15、限制条件较外模态的情形宽松很多，主要是速度较快的外重力波已经在外模态中考虑了。一般 取值3050即可(42).时间步长的其它限制条件 对于动量或标量还有其它一些时间限制条件 其中 或 以及 其中 分别为地球角速率和地理纬度(43)(44).侧开边界条件（bcond）陆地及岸线是由dum、dvm、fsm控制的，在陆地上这些变量的值设为0，有水的地方设为1 子程序bcond（idx），idx=1对应的水位边界条件；idx=2对应垂向平均流速边界条件；idx=3对应三维水平流速边界条件；idx=4对应温盐边界条件；idx=5对应垂向流速边界条件；idx=6对应湍流动能和湍流混合长度边界条件 模型对开

16、边界条件的要求很高，而开边界条件本身又具有很大的不确定性，因此有必要对模型的外模态和内模态的开边界分别进行处理.数值积分流程图.Formula Boundary Code Inflow condition EAST uaf(im,j)=2*bc(j)/(h(im,j)+elf(im,j)+h(imm1,j+elf(imm1,j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWEST uaf(2,j)=2*bc(j)/(h(1,j)+elf(1,j)+h(2,j)+elf(2,j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTH vaf(i,jm)=2*

17、bc(i)/(h(i,jm)+elf(i,jm)+h(i,jmm1)+elf(i,jmm1)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTH vaf(i,2)=2*bc(i)/(h(i,1)+elf(i,1)+h(i,2)+elf(i,2)elf(i,1)=elf(i,2)uaf(i,1)=set Elevation condition h=BC EAST elf(imm1,j)=bc(j)elf(im,j)=elf(imm1,j)cosmeticuaf(im,j)=uaf(imm1,j)vaf(im,j)=setWEST elf(2,j)=bc(j)uaf(2,

18、j)=uaf(3,j)vaf(1,j)=setNORTH elf(i,jmm1)=bc(i)elf(i,jm)=elf(i,jmm1)cosmeticvaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTH elf(i,2)=bc(i)vaf(i,2)=vaf(i,3)uaf(i,1)=set外模态开边界条件一.Formula Boundary CodeRadiation EAST uaf(im,j)=sqrt(grav/h(imm1,j)*el(imm1,j)+bc(j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWEST uaf(2,j)=-s

19、qrt(grav/h(2,j)*el(2,j)+bc(j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTH vaf(i,jm)=sqrt(grav/h(i,jmm1)*el(i,jmm1)+bc(i)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTH vaf(i,2)=-sqrt(grav/h(i,2)*el(i,2)+bc(i)elf(i,1)=elf(i,2)uaf(i,1)=setRadiation EAST gae=dte*sqrt(grav*h(im,j)/dx(im,j)uaf(im,j)=gae*ua(imm1,j)+(1.-gae

20、)*ua(im,j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWEST gae=dte*sqrt(grav*h(2,j)/dx(2,j)uaf(2,j)=gae*ua(3,j)+(1.-gae)*ua(2,j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTH gae=dte*sqrt(grav*h(i,jm)/dy(i,jm)vaf(i,jm)=gae*va(i,jmm1)+(1.-gae)*va(i,jm)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTH gae=dte*sqrt(grav*h(i,2)/dy(i,2

21、)vaf(i,2)=gae*va(i,3)+(1.-gae)*va(i,2)elf(i,1)=elf(i,2)uaf(i,1)=set外模态开边界条件二流速水位流速.Formula Boundary CodeRadiationEASTgae=dte*sqrt(grav*h(imm1,j)/dx(imm1,j)elf(imm1,j)=gae*el(imm2,j)+(1.-gae)*el(imm1,j)elf(im,j)=elf(imm1,j)uaf(im,j)=uaf(imm1,j)vaf(im,j)=setWEST gae=dte*sqrt(grav*h(2,j)/dx(1,j)elf(2,j

22、)=gae*el(3,j)+(1.-gae)*el(2,j)uaf(2,j)=uaf(3,j)vaf(2,jm)=setNORTH gae=dte*sqrt(grav*h(i,jmm1)/dy(i,jmm1)elf(i,jmm1)=gae*el(i,jmm2)+(1.-gae)*el(i,jmm1)elf(i,jm)=elf(i,jmm1)vaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTH gae=dte*sqrt(grav*h(i,2)/dy(i,2)elf(i,2)=gae*el(i,3)+(1.-gae)*el(i,2)vaf(i,2)=vaf(i,3)，ua

23、f(i,1)=set Cyclic EAST(I=IM)elf(im,j)=elf(3,j)uaf(im,j)=uaf(3,j)，vaf(im,j)=vaf(3,j)WEST(I=1)elf(1,j)=elf(imm2,j)，elf(2,j)=elf(imm1,j)uaf(2,j)=uaf(imm1,j)，vaf(2,j)=vaf(imm1,j)NORTH(J=JM)elf(i,jm)=elf(i,3)uaf(i,jm)=uaf(i,3)vaf(i,jm)=vaf(i,3)SOUTH(J=1)elf(i,1)=elf(i,jmm2)，elf(i,2)=elf(i,jmm1)uaf(i,2)=u

24、af(i,jmm1)vaf(i,2)=vaf(i,jmm1)外模态开边界条件三水位.Formula Boundary CodeInflow conditionEAST uf(im,j,k)=bc(j,k)vf(im,j,k)=setU=BC WEST uf(2,j,k)=bc(j,k)vf(1,j,k)=set NORTH vf(i,jm,k)=bc(i,k)uf(i,jm,k)=setSOUTH vf(i,2,k)=bc(i,k)uf(i,1,k)=setRadiation:EAST gai=sqrt(h(im,j)/hmax)uf(im,j,k)=gai*u(imm1,j,k)+(1.-g

25、ai)*u(im,j,k)vf(im,j,k)=setWEST gai=sqrt(h(2,j)/hmax)uf(2,j,k)=gai*u(3,j,k)+(1.-gai)*u(2,j,k)vf(1,j,k)=setNORTH gai=sqrt(h(i,jm)/hmax)vf(i,jm,k)=gai*v(i,jmm1,k)+(1.-gai)*v(i,jm,k)uf(i,jm,k)=setSOUTH gai=sqrt(h(i,2)/hmax)vf(i,2,k)=gai*v(i,3,k)+(1.-gai)*v(i,2,k)uf(i,1,k)=set内模态开边界条件（一）.Formula Boundar

26、y CodeUpstream advection on T or SEAST uf(im,j,k)=t(im,j,k)-dti/(dx(im,j)+dx(imm1,j)*(u(im,j,k)+abs(u(im,j,k)*(t(im,j,k)-t(imm1,j,k)+(u(im,j,k)-abs(u(im,j,k)*(tbe(j,k)-t(im,j,k)WEST uf(1,j,k)=t(1,j,k)-dti/(dx(1,j)+dx(2,j)*(u(1,j,k)+abs(u(1,j,k)*(t(1,j,k)-tbw(j,k)+(u(1,j,k)-abs(u(1,j,k)*(t(2,j,k)-t(1

27、,j,k)NORTH uf(i,jm,k)=t(i,jm,k)-dti/(dy(i,jm)+dy(i,jmm1)*(v(i,jm,k)+abs(v(i,jm,k)*(t(i,jm,k)-t(i,jmm1,k)+(v(i,jm,k)-abs(v(i,jm,k)*(tbn(i,k)-t(i,jm,k)SOUTH uf(i,1,k)=t(i,1,k)-dti/(dy(i,1)+dy(i,2)*(v(i,1,k)+abs(v(i,1,k)*(t(i,1,k)-tbs(i,k)+(v(i,1,k)-abs(v(i,1,k)*(t(i,2,k)-t(i,1,k)Cyclic Much the same a

28、s the external case except replace uaf with uf,etc.and t,s,q2 and q2l are handled similar to elf内模态开边界条件（二）.第三节POM 模型的实际应用渤海河流羽状流的数值模拟.渤海及入海河流.河流径流输入模型 将河流视为点源，而点源的模拟采用的是一个水平网格的面积 径流动量对动量方程的贡献可以忽略不计 受河流径流影响的二维连续方程.三维连续方程可表示为 内外模式相容.30 天时纯浮力羽流条件下渤海的盐度分布，等值线间隔为0.5 psu，最外层等值线为30 psu。（a）海水表层；（b）海水底层.30

29、天时潮致羽流条件下渤海的盐度分布，等值线间隔为0.5 psu，最外层等值线为30 psu。（a）海水表层；（b）海水底层.东风 南风风场加入10天后渤海表层海水流场.东风 南风风场加入10 天后渤海表层盐度的水平分布状况。等值线间隔为0.5 psu。.海河口附近站点（11745E，38 55N）盐度变化的时间序列海水表层 海水底层.7月份渤海表层盐度的水平分布状况，等值线间隔为1.0 psu。计算结果(左图）；实测结果（右图）.学习要点 了解基于sigma坐标变换的海水运动方程 了解POM模型的主要特征 熟悉模型运算的流程 掌握M-Y湍流闭合模式的基本特点 针对不同研究对象学会选择POM模式的开边界条件 河流模型与POM模型的耦合.