北部湾三维潮流数值模拟.pdf

第 19 卷 第 2 期 1997-03海 洋 学 报 ACTA OCEANOLOGICA SINICAVol. 19, No. 2 M arch, 1997北部湾三维潮流数值模拟夏华永 殷忠斌 郭芝兰 陈 明 剑( 广西海洋监测预报中心, 北海) (广西科委, 南宁)摘 要 模型以经 Sigma 坐标变换后具有自由表面的三维非线性 Navier-Stokes 方程为基本方程, 以分裂算子法剖分动量方程, 用全隐差分格式求解连续方程得到自由表面, 最 后计算完整的速度场. 此外, 采用紊流封闭模型求解垂直方向紊动粘滞系数, 准确获得了摩擦影响层的潮流结构. 应用本模型计算了北部湾的潮波运动, 重现了 K1、 M2分潮潮波系统并揭示了潮流的空间结构特征.关键词 三维模型 紊流封闭模型 北部湾潮流本文于 1996-01-02 收到,修改稿于 1996-08-02 收到.1)Wang J.L A Myasak and R G Ingram.A three-dimensional numerical simulation of Hudson Bay summer circulation.J.p. o. , 1994 ( in press) .2) Wang J and M Ikeda. On stability of finite difference schemes for inertial oscillations in ocean general circulation models.Submitted to International J. Numo M ethods in Fluids, July, 1994.前 言北部湾是南海的主要海湾. 由于潮流是海洋资源开发、 海洋工程发展及海洋环境管理中必须考虑的动力因素, 因此研究北部湾潮流的空间结构具有重要意义. 迄今, 已有一些北部湾的 二维潮波研究[ 1, 2], 但三维潮波研究尚不多见[ 3], 且尚未深入探讨三维潮波的空间结构特征. 在其他海域的三维潮波研究中, 窦振兴等[ 4]采用 Sigma 坐标系下的三维模型, 克服了固定分层模型[ 5]在浅水部分垂直方向分辨率很低的缺点, 在应用中获得了良好的效果, 但在数值方法上,以往 Sigma 坐标系下的三维模型[ 4, 6]都采用过程分裂法, 先用二维模型求解快过程——表面重力波——获得一个自由表面, 再用三维模型求解慢过程——内重力波——获得三维流场. 过 程分裂法有较好的计算精度, 但计算量较大, 程序也较复杂. 本文采用分裂算子法1, 2), 不须要将三维潮流中的快过程与慢过程分开, 可直接求解自由表面及速度场.应用本模型对北部湾 K1、 M2分潮做了 6 层模式的模拟, 潮波系统的模拟结果与观测值符 合较好. 并获得了潮流空间结构特征.1 基本方程组1. 1 潮流动力学方程组北部湾是一个较小的海湾, 海平面可用笛卡尔直角坐标系表明, 垂直方向采用 ? 坐标变 换, ?= ( z- ?) / ( H + ?) , H ( x , y) 为平均海平面到海底的距离, ?( x, y, t) 为海面升降, 在 z = ? 时, ?= 0, 在 z= - H 时, ?= - 1, 在 ?坐标系下, 正压环流的控制方程组可写为 ?? ?t+?uD ?x+?uD ? y+?w ??= 0,( 1)?uD ?t+?u2D ?x+?uvD ?y+?uw ??- f vD+ gD?? ?x=? ??KM D?u ??+? ?x2AMD?u ?x+? ?yAMD?u ?y+? v ?x,( 2)?uD ?t+?uvD ?x+?v2D ?y+?vw ??+ f uD+ gD?? ?y=? ??KM D?v ??+? ?t2AMD?v ?y+? ?xAMD?u ?y+?v ?x,( 3)式中, D 为水深, D= H + ?; f 为科氏力系数, AM为水平方向紊动粘滞系数; KM分垂直方向紊 动粘滞系数, ?为坐标变换后产生的一个垂向速度, 实际上, 它垂直于 Sigma 分层, 三维笛卡尔 直角坐标系中的垂向速度 ? 与 Sigma 坐标系中的速度场有如下关系:?= w - u ??D ?x+? ? ?x- v ?? D ?y+?? ?y-??D ?t+?? ?t.( 4)1. 2 紊流封闭方程许多三维模型都将垂直紊动粘滞系数 KM取为常量, 实际上垂直紊动粘滞系数是随水深 变化的. 叶安乐[ 7]在探讨潮流椭圆长轴方向随深度变化的特征时发现, 最大流速方向对 KM的 取值是敏感的, 不同的 KM值选用可得取截然不同的结论; 叶安乐[ 8]也发现潮流最大流速发生 时刻随深度的变化率也受 KM值的大小影响; 沈育疆等[ 9]在模拟中将 KM取为常量得出了一 个过厚的摩擦影响层, 因此将 KM取为常量不能完全真实的反映潮流的铅直结构. 本文采用近 年来应用效果很好的紊流封闭模型[ 10]计算垂计紊动粘滞系数 ?q2D ?t+?uq2D ?x+?vq2D ?y+??q2 ??=? ??Kq D?q2??+2KM D?u2??+?v ??2+2g ?0KH?? ??-2Dq2 A1+? ?xDAH?q2 ?x+? ?yDAH?q2 ?y,( 5)?q2lD ?t+?uq2lD ?x+?vq2D ?y+?wq2l ??=? ??Kq D?q2l ??+E1lKM D?u ??2 +?v ??2+lE1g ?0KH?? ? ?-Dq3 B1W-+? ? xDAH? q2l ?x+? ?yDAH?q2l ?y,( 6)式中, q2/ s 为紊流动能; l 为紊流宏观尺度( turbulence macroscale) ; ?0为参考密度, ?0取为22海洋学报 19 卷1 000kg/ m3; ?为海水密度, 对于正压环流海水密度为常量; W-为面壁近似函数( wall proximi- ty function) ,W-= 1+ E2(l kL)2,( 7)( L)- 1= ( ?- z)- 1+ ( H + z)- 1,( 8)AH为水平方向的物质扩散系数, 垂向紊动粘滞系数 KM和垂向扩散系数 KH、 Kq分别由下列 公式确定KM= lqSM,( 9)KH= lqSH,( 10)Kq= lqSq,( 11)SM、 SH、 Sq为稳定性函数( stability function) , SM、 SH、 Sq由下列方程组求解:GM=l2 q2D?u ??2 +?v ??21/ 2 ,( 12)GH=l2q2Dg ?0?? ? ?,( 13)SM( 6A1A2GM) + SH( 1- 2A2B2GH- 12A2A2GH) = A2,( 14)SM( 1 + 6A21GM- 9A1A2GH) - SH( 12A21GH+ 9A1A2GH)= A1( 1 - 3C1) ,( 15) Sq= 0. 20,( 16)在正压环流中, 海水密度为常数, 在式( 13) 中, GH= 0. A1、 A2、 B1、 B2、 C1、 E1、 E2为经验常数, 其 值由试验所得,( A1, A2, B1, B2, C1) = ( 0. 92, 0. 74, 16. 6, 10. 1, 0. 08) ,( E1, E2) = ( 1. 8, 1. 33) .1. 3 边界条件动力学边界条件为: 在自由表面上满足?0KM D?u ??,?v ??= ( ?ax, ?ay, ?→0,( 17)对于纯天文潮, ?a= 0; 在近海底处满足 KM D?u ??,?v ??= Cz( u2+ v2)1/ 2( u, v) , ?→- 1,( 18)式( 18) 中,Cz= maxk2 ln2( z/ z 0), 0. 002 5 ,( 19)k 是卡门常数, k= 0. 4; z0是海底粗糙度, 在潮流模拟中 z0可取为 0. 002～0. 01m; z 是离海底 最近网格与海底的距离. 垂直边界条件满足 ?( 0) = ?( - 1) = 0.( 20)232期 夏华永等: 北部湾三维潮流数值模拟海岸边界条件上满足 uN( x, y, ?, t) = 0,( 21)N 为岸线外法线方向. 在开边界上, 给定潮位值?= ? cos2? Tt- ? .( 22)2 计算模式为表达离散方程的方便起见, 先定义如下差分和求和算子:f ( x, y, ?, t)x=fx+?x 2, y, ?, t + fx-?x 2, y, ?, t2,( 23)?xf ( x, y, ?, t) =fx +?x 2, y, ?, t - fx-?x 2, y, ?, t?x,( 24)图 1 变量在空间交错网格中的位置f ( x, y, ?, t)xy= f , ( x, y, ?, t)xy = f , ( x, y, ?, t)yx.( 25) 变量分布在交错网格上, 其相对位置如图 1 所示. 在对方程进行离散时, 先将动量方程( 2) 、 ( 3) 分裂成如下形式:Dnun+ 1/ 2- ( Du)n?t+?u2D?x+?uvD ?y+?uw ? ?- f vD=? ??KM D?u ??+? ?x2AMD?u ?x+? ? yAMD?u ?y+?v ?x,( 26)( Du)n+ 1- Dnun+ 1/ 2?t= - Gd?? ? x,( 27)Dnun+ 1/ 2- ( Du)n+?t+? u vD ? x+?v2D y+?vw ??+ f uD=?u ??KM D?v ??+? ?y2AMD?v ?y+? ?xAMD?u ?y+?v ?x,( 28)( Dv)n+ 1- Dnvn+ 1/ 2?t= - gD?? ?y.( 29)在 u、 v 所处的格点位置处, 分别将方程( 26) 、 ( 27) 、 ( 28) 、 ( 29) 、 离散成如下代数方程:Dnun+ 1/ 2- ( Du)n?t+ ?x( Dxxuux) + ?y( Dxyvuy) + ? ?( ?xu?) - f vyDx= ??KMxDx??( u)n+ 1/ 2+ ?x2 AMDx?x( u)+ ?y{AMDy[ ?y( u) + ?x( v) ] },( 30)( Du)n+ 1- Dnun+ 1/ 2?t= - gDx?x( ?)n+ 1, ( 31)24海洋学报 19 卷Dnvn+ 1/ 2- ( Dv)n?t+ ?x( Dxy u vx) + ?y( Dyy v vy) + ??( ?yv?) + f uxDy= ??KMyDy??( v)n+ 1/ 2+ ?y[ 2AMDy?y( v) ] + ?x{AMDx[ ?y( u) ] + ?x( v) ] },( 32)=( Dv)n+ 1- Dnvn+ 1/ 2?t= - G Dy?y( ?)n+ 1. ( 33)在方程式( 30) 、 ( 32) 中, 水平扩散项中的流量通量 Du、 Dv以及科氏力中的速度 u、 v 取 n 时刻的值或迭代计算所得的 n+ 1 时刻的值, 其他没有时间指标的变量全部为 n 时刻的值. 连 续方程( 1) 和紊流模型( 5) 、 ( 6) 离散成如下形式:?n+ 1- ?n?t+ ?x( Dxu)n+ 1+ ?y( Dyv)n+ 1+ ??( ?)n+ 1= 0,( 34)( q2D)n+ 1- ( q2D)n?t+ ?x( u?Dxq2x+ ?y( v?Dyq2y) + ??( ??q2?)= ??Kq D??q2n+ 1 +2KM D[ ( ??ux)2+ ( ??uy)2] +2g ?0KH???-2Dq3 A1+ ?x( DxA H?? xq2) + ?y( DyAH?? yq2) ,( 35)( q2lD)n+ 1- ( q2lD)n?t+ ?x( u?Dxq2lx) + ?y( v?Dyq2ly) + ??( ??q2l?)= ??Kq D??q2ln。