收藏本站 劰载中...网站公告 | 吾爱海洋论坛交流QQ群:835383472

数值海洋与大气模式(三):从0到1实现浅水模式(下) ...

[复制链接]
二维浅水模式

  上文讨论了浅水方程在一维情况下的求解,本文探讨浅水方程在二维情况下的求解,并考虑加上底摩擦和风应力的情况。

  首先考虑最简单的二维浅水方程,形式如下。

& y. o% ~* J# B! ~" e0 S! U
                               
登录/注册后可看大图

  还是同样的过程,对其做离散化,得到下式。


. V/ r% }$ s/ w- ^                               
登录/注册后可看大图

  上文提到,为了方便实现中央差分来取得二阶精度,选择了交错网格,使得

2 n/ O/ r3 L1 e& Q% F# k7 {
                               
登录/注册后可看大图
# n6 y/ {9 M! A7 \: s: _: V* U7 h
                               
登录/注册后可看大图
总是相差半个网格距离。在二维的情况,速度是

1 e! E# ^' I' L7 s4 d                               
登录/注册后可看大图
! t( `$ d1 I* h( L: q
                               
登录/注册后可看大图
两个变量,因此我们需要考虑速度

' @4 B; @* m6 S+ J, {! C% Q# @                               
登录/注册后可看大图

放在网格的哪个位置。从上式的第二个式子可以看出,式子右边是

" e' A" T0 Z+ P3 z) q
                               
登录/注册后可看大图

  \: V) x6 F* U! _' Y1 i                               
登录/注册后可看大图
的导数,因此为了让对

* n2 \5 \" [& h                               
登录/注册后可看大图
的中央差分落在和
5 i+ e' H* ?( u0 X9 \. _: z2 o6 i) q; M
                               
登录/注册后可看大图
同样的位置上,

  P" [  y& G9 N1 W8 e1 i) W9 Z                               
登录/注册后可看大图
也需要和
$ _) f; X  n  L
                               
登录/注册后可看大图
间隔半个网格距。这就是大名鼎鼎的Arakawa C网格。在使用交错网格的时候,我们要明确一个准则,就是要确保计算时各个变量都要在相同的空间位置上。比如要算

) M, {$ e  Y7 n                               
登录/注册后可看大图
在t+1时刻的值,就需要

  i5 e* W, A* ]: N                               
登录/注册后可看大图
两侧的
  @' K$ m3 a# G, R3 x
                               
登录/注册后可看大图
做中央差分,这样确保等号左边的
4 k4 O. F% E. k7 ]* j0 q
                               
登录/注册后可看大图
和等号右边的导数是相同的空间位置。这一点在之后讲解更复杂的模式中会体现的更明显。例如POM模式,会有更多的参变量,为了保证处于同一空间位置上需要做平均运算。而本文的后半部分,在考虑底摩擦的浅水方程中也涉及到了。


2 [+ O; l, {' `# x& ^/ K4 h# d; |


& N2 E$ `7 f5 d: R% N                               
登录/注册后可看大图

  使用了Arakawa C网格后,可以神奇的发现,在这种情况下,只需要讨论两个边界的速度即可,如下图所示。考虑到


5 K2 G5 a* t- X9 r1 I0 m' d                               
登录/注册后可看大图
是水平方向的流速,而东边界(图中右边)的
" N! ]: ^6 K3 R7 c; a* [
                               
登录/注册后可看大图
是在陆地边界上,所以要使

1 E# S% x# O7 _# n3 v9 @+ }                               
登录/注册后可看大图
时变为0,不能让水质点穿出边界,而

' x$ f6 `$ C6 Y/ t                               
登录/注册后可看大图
或者
9 I5 c! g" W( n0 t% j
                               
登录/注册后可看大图
时则不需要给予约束,这样就不用考虑西边界,因为西边界没有

7 x* h# }& ~/ H! \- E                               
登录/注册后可看大图
的格点。对北边界的
" I! b9 I$ c/ |; P$ o  }
                               
登录/注册后可看大图
处理方式一样,在此不做赘述。

, j! H: u: A) M

3 v4 i8 d4 n9 g5 I5 s+ l$ ]) G( N& Y5 i
                               
登录/注册后可看大图

un(j,k)=0;uu=u(j,k);duu=du(j,k);if(wet(j,k)==1)3 z* \; Z% W6 k* n
    if((wet(j,k+1)==1)||(duu>0))  i/ g+ {1 }5 {# U' V3 r2 ^
        un(j,k)=uu+duu;: Q3 Y# c. K$ I" c, [  A! ]
    endelse
$ Y& s7 C; h* S    if((wet(j,k+1)==1)&&(duu<0))- e! c2 I3 ], ^& K
        un(j,k)=uu+duu;
- v4 d; U& S9 X) d    endend

  对


  q" J) q# V( t4 p3 _% `: F                               
登录/注册后可看大图
求解的代码实现过程如上,其中if就是在判断边界和速度方向,把不在边界和在边界但是
7 u& ^" N% @+ G# L5 k
                               
登录/注册后可看大图
不大于0的值赋计算结果,使得在边界且大于0的点赋为0。

2 R6 O. x# g% W0 b                               
登录/注册后可看大图
的计算和
$ U( O* i+ a& B7 T$ N1 Z$ v
                               
登录/注册后可看大图

的计算同理。


5 M( m" k, o& {( g0 ?) a

  接下来,开始考虑更复杂的形式,引入风应力和底摩擦。方程组的形式立马变得更加复杂起来。


8 A7 C: y, u2 Q" v4 [3 j* d                               
登录/注册后可看大图

  底摩擦和风应力这两项由上式可以看出,是相减的关系,因此可以在计算时,分别看成两项来对待。如果只考虑底摩擦,不考虑其他项的话,方程是如下形式。


. R3 c5 ~# m: I2 @                               
登录/注册后可看大图

  对其进行半隐式差分得到如下形式。显式差分应用于底摩擦会出现数值稳定性问题。


  ]0 r! o$ b7 L9 V                               
登录/注册后可看大图

  考虑完了底摩擦,还需要考虑风应力问题。在原方程去掉底摩擦项后再进行差分得到下式,其中下标j和k分别代表x和y方向。


7 M# K2 Z8 [4 e; ~; j                               
登录/注册后可看大图

  将其求完后回代到推的底摩擦差分方程中即可。

' I9 Q3 P) M# c/ N2 C
                               
登录/注册后可看大图

  为了表达式清晰简洁,引入了


/ N- }3 b% c6 |3 W6 X7 o4 {                               
登录/注册后可看大图


$ y, t- U% [3 a                               
登录/注册后可看大图

,具体表达式如下。

$ s, U( u) ?! Q# W: H3 T+ M


2 _7 V) H2 x& g& [. w2 ~* J# I7 C                               
登录/注册后可看大图

  值得注意的是,上式中出现了


. ~( i/ q- F, |  F2 @: H0 X4 N                               
登录/注册后可看大图
7 C1 H6 ]8 b: B3 |  d6 X
                               
登录/注册后可看大图
这种变量。这就是前文已经提到的那个问题,Arakawa C网格的

( L$ l' n0 [6 k' ?' J                               
登录/注册后可看大图

% @. f5 l6 w6 A# U5 g, J! ~                               
登录/注册后可看大图
在网格的不同位置处,但是做运算只能在相同网格处。因此
  S+ v/ L( @, g1 @# ]4 r, F$ r# K8 z
                               
登录/注册后可看大图

的意思是v网格所在位置的速度


# L: y+ J( I4 y) V, x                               
登录/注册后可看大图
,而
* f( [/ {' O) M; v1 E/ p
                               
登录/注册后可看大图
的意思是u网格所在位置的速度
& M, ]: c* }5 S. t
                               
登录/注册后可看大图
。以

9 }/ k3 `( i6 h+ m3 e                               
登录/注册后可看大图
为例,如下图所示,可由周围的四个

! f! K. }$ S( y2 O8 n+ }, C                               
登录/注册后可看大图
平均求得。

u_v = 0.25*(v(j+1,k)+v(j,k)+v(j,k-1)+v(j+1,k-1))


, V: }* d* t' P4 s, J, s" `. e                               
登录/注册后可看大图

  这样就实现了二维考虑底摩擦和风应力的浅水方程求解。之前说过,干网格和湿网格分别对应网格中的陆地和水,我们在计算时只需要计算流体部分,而不计算陆地的部分。之前我们提到的案例中,陆地边界都在网格的四周。倘若需要模拟的场景是大洋中的一个小岛,即陆地出现在网格内部时,会有什么不同?倘若大洋中的小岛海拔高度有限,在模拟台风时,小岛被水淹没时,干网格和湿网格该怎么转换?当考虑了干湿网格的转换之后,一个具备基本功能的浅水模式就完成了。

1 s4 g/ e9 C1 x5 j! \
                               
登录/注册后可看大图

  以一维的情况为例,从图中可以看出来,这个扰动显然会在某些时刻漫过这个岛屿。为了解决这个问题,我们只需要在每轮时间步迭代的最后,加上干湿网格的更新即可。需要定义一个比较小的hmin,意思是当陆地上的积水高度超过hmin就意味着被淹没,下次计算时就把这个网格点当做水面,即把wet(k)这一点赋值为1。

for k=1:N
( A+ t* G0 w$ r    h(k)=hzero(k) + eta(k);
% @& y0 G. F5 D    wet(k)=1;9 \; `% d3 V: y) s, x
    if(h(k)<hmin)  w0 l: p, C& N9 p* r9 |
        wet(k)=0;7 C& W/ d1 x3 u& e) D* `  l  S
    end
; _) `/ R, Y: \" F* o: F    u(k)=un(k);end

版权声明

  本文创作的初衷是用于帮助数值模式的学习者。欢迎转载,转载请私信并注明作者和出处,请勿用于任何商业用途。

参考书目

Ocean Modelling for Beginners. Jochen Kämpf. Springer Berlin Heidelberg, 2009.

0 f- N2 `. L, L9 T) F/ L
回复

举报 使用道具

相关帖子

全部回帖
暂无回帖,快来参与回复吧
懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
jnsbmj
活跃在2022-11-6
快速回复 返回顶部 返回列表