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

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

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

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


# I0 J" U" ~, x                               
登录/注册后可看大图

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


! w$ |4 w$ K  W: @1 M& O3 H                               
登录/注册后可看大图

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

0 z0 r9 W5 {" l/ J2 w6 [- |  G+ @
                               
登录/注册后可看大图

% T' \( b5 r# l# C                               
登录/注册后可看大图
总是相差半个网格距离。在二维的情况,速度是
0 E& i, F: L$ y
                               
登录/注册后可看大图

! e# d+ L" d- j/ H/ D                               
登录/注册后可看大图
两个变量,因此我们需要考虑速度

9 ]2 J9 n& o) F0 ^6 D                               
登录/注册后可看大图

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


+ j+ |3 M$ R; J5 m% @9 L                               
登录/注册后可看大图
- S7 `, _$ d  L- e( |* |" j
                               
登录/注册后可看大图
的导数,因此为了让对
3 q% g/ W& \6 U. b
                               
登录/注册后可看大图
的中央差分落在和

- W- \# E) |  U. I3 G                               
登录/注册后可看大图
同样的位置上,
. d1 q* [2 b$ n$ ~, A# n7 P
                               
登录/注册后可看大图
也需要和

' @) y" y3 B( W: ~                               
登录/注册后可看大图
间隔半个网格距。这就是大名鼎鼎的Arakawa C网格。在使用交错网格的时候,我们要明确一个准则,就是要确保计算时各个变量都要在相同的空间位置上。比如要算
. e3 j0 Z. R' f6 m4 C  s
                               
登录/注册后可看大图
在t+1时刻的值,就需要
3 c- ~0 |) L# c1 h$ ~- h
                               
登录/注册后可看大图
两侧的
3 |; r' i! a9 k! ~: g
                               
登录/注册后可看大图
做中央差分,这样确保等号左边的
. ^" g, l5 `1 \! E
                               
登录/注册后可看大图
和等号右边的导数是相同的空间位置。这一点在之后讲解更复杂的模式中会体现的更明显。例如POM模式,会有更多的参变量,为了保证处于同一空间位置上需要做平均运算。而本文的后半部分,在考虑底摩擦的浅水方程中也涉及到了。


& f* o! m( \- |8 W" B


4 }7 L4 \6 ^# I. u! L: T% q                               
登录/注册后可看大图

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


+ P+ @: m3 ]+ a. G2 ?4 u                               
登录/注册后可看大图
是水平方向的流速,而东边界(图中右边)的

2 y3 h: {3 }) d' c                               
登录/注册后可看大图
是在陆地边界上,所以要使

+ [" B' c/ d( D; t- E0 C5 ]                               
登录/注册后可看大图
时变为0,不能让水质点穿出边界,而
+ }' `" ?* P3 p6 a6 F/ Z# d' V( E( j4 u$ W
                               
登录/注册后可看大图
或者

& j( g8 @# U- A' Q$ K3 U                               
登录/注册后可看大图
时则不需要给予约束,这样就不用考虑西边界,因为西边界没有
, ^4 m* J+ s- F4 w9 B' z
                               
登录/注册后可看大图
的格点。对北边界的

9 q* s, V, h  y# R                               
登录/注册后可看大图
处理方式一样,在此不做赘述。


2 N3 j: R( P9 H


$ e' D& M* n6 C                               
登录/注册后可看大图

un(j,k)=0;uu=u(j,k);duu=du(j,k);if(wet(j,k)==1)2 b8 u/ e, `- U; h# }5 I3 N
    if((wet(j,k+1)==1)||(duu>0)): P5 ~* J& P/ Z
        un(j,k)=uu+duu;: l; {; Z3 D8 [, z( Y/ n
    endelse  p( N& U' P- k
    if((wet(j,k+1)==1)&&(duu<0))
" a" }9 |6 k, i0 u, _* Y8 X        un(j,k)=uu+duu;6 _6 k9 E) d: V/ |# G$ @
    endend

  对

$ Y. \& i$ g9 b+ T0 {5 o  [: K% t& \' q
                               
登录/注册后可看大图
求解的代码实现过程如上,其中if就是在判断边界和速度方向,把不在边界和在边界但是
4 P& B  y1 L( V: j+ L' S  @. J
                               
登录/注册后可看大图
不大于0的值赋计算结果,使得在边界且大于0的点赋为0。

0 e( b0 I# Q6 \7 M- x                               
登录/注册后可看大图
的计算和
1 _+ m8 |8 V1 [8 m
                               
登录/注册后可看大图

的计算同理。

9 d) b0 S* `5 ^) W

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

" K7 r; z& {& n4 y9 u6 @
                               
登录/注册后可看大图

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


2 F2 n: P& @. ?. P                               
登录/注册后可看大图

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


: t4 a7 m, p4 |                               
登录/注册后可看大图

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


1 N8 U# v3 v0 N7 ]                               
登录/注册后可看大图

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


# v6 o# L. J4 C7 A+ l% ~* A                               
登录/注册后可看大图

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

) X5 d' K. o2 B1 z# D; v
                               
登录/注册后可看大图

3 o; B  N1 l$ R) A# D: ^
                               
登录/注册后可看大图

,具体表达式如下。


( Y, Q9 S& M% S9 X. W$ F% J2 e


0 R4 |5 n' g; S! X                               
登录/注册后可看大图

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

! ]/ z: S" |$ q7 G
                               
登录/注册后可看大图

7 {$ L$ L6 C$ s* l                               
登录/注册后可看大图
这种变量。这就是前文已经提到的那个问题,Arakawa C网格的

6 F" p( W; e- t8 v, \3 {6 g# K                               
登录/注册后可看大图
+ ~$ c* p4 M0 N# u9 U
                               
登录/注册后可看大图
在网格的不同位置处,但是做运算只能在相同网格处。因此
7 D2 j, K, _- Q+ K9 L  ]& S
                               
登录/注册后可看大图

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


0 t  u, ?: I0 S; N7 h, c0 Q( r# p2 T                               
登录/注册后可看大图
,而

+ k8 ?* A& f& r" S                               
登录/注册后可看大图
的意思是u网格所在位置的速度

7 s) O% K3 l& t& @$ h                               
登录/注册后可看大图
。以
1 M7 D( v" M. J) Z' K0 M  G% v
                               
登录/注册后可看大图
为例,如下图所示,可由周围的四个
; x6 M# o, q. K- O- h3 Y  r. ]
                               
登录/注册后可看大图
平均求得。

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


% w  E) ^  K: V3 u1 |" [) d- ?8 x7 i                               
登录/注册后可看大图

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


1 s0 [! Y1 W) n3 J/ Q                               
登录/注册后可看大图

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

for k=1:N5 ?1 G9 m6 l! ^6 s* c) ]
    h(k)=hzero(k) + eta(k);: n7 K! P+ V. @) \5 g  h' X* ~
    wet(k)=1;& \# X6 K' Y3 _6 p# x
    if(h(k)<hmin)" y8 c6 Z+ H0 L
        wet(k)=0;
& e8 m) O) q0 N    end7 ?7 i3 @1 e# T2 X- w
    u(k)=un(k);end

版权声明

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

参考书目

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

; K# ]/ l+ e: R0 c1 a
回复

举报 使用道具

相关帖子

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