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

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

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

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

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

3 b2 E2 E# o& n/ Y( y$ B
                               
登录/注册后可看大图

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

7 f4 ^1 y" y! i7 m% g. ^9 l
                               
登录/注册后可看大图

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

+ Y2 P+ W! C6 C' ^  E- m  l+ {- x
                               
登录/注册后可看大图
9 z6 J6 u( k3 Z' w
                               
登录/注册后可看大图
总是相差半个网格距离。在二维的情况,速度是

& T% D: ~: Z( G; A1 F3 A                               
登录/注册后可看大图
3 b3 t0 h! N3 @+ H( m+ I7 a) H
                               
登录/注册后可看大图
两个变量,因此我们需要考虑速度
+ U8 w: `3 k6 K+ h. L4 C! A
                               
登录/注册后可看大图

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

% |) X/ G4 b6 z' B( t) t& u
                               
登录/注册后可看大图
& w4 ?0 B( S" T2 \+ f2 s+ ~' J1 E& ^+ Q
                               
登录/注册后可看大图
的导数,因此为了让对

0 ?6 h! p0 r9 _: K9 Z# g1 ^2 U                               
登录/注册后可看大图
的中央差分落在和

2 n! _: y, K4 p: W                               
登录/注册后可看大图
同样的位置上,
( ~9 Q: J5 D% j- |. n
                               
登录/注册后可看大图
也需要和
4 ^3 [/ p4 S* E5 R3 J
                               
登录/注册后可看大图
间隔半个网格距。这就是大名鼎鼎的Arakawa C网格。在使用交错网格的时候,我们要明确一个准则,就是要确保计算时各个变量都要在相同的空间位置上。比如要算

5 V! k, k9 m4 c5 w                               
登录/注册后可看大图
在t+1时刻的值,就需要

8 Q$ t/ q5 Q& l7 h0 O# L" l                               
登录/注册后可看大图
两侧的

+ g( ?( Z( G; r. q3 M                               
登录/注册后可看大图
做中央差分,这样确保等号左边的
0 j' R8 G; ?9 p3 y3 u
                               
登录/注册后可看大图
和等号右边的导数是相同的空间位置。这一点在之后讲解更复杂的模式中会体现的更明显。例如POM模式,会有更多的参变量,为了保证处于同一空间位置上需要做平均运算。而本文的后半部分,在考虑底摩擦的浅水方程中也涉及到了。


5 w5 {8 k7 }5 u4 D


. `5 z/ t  i  c                               
登录/注册后可看大图

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


& u! u- n& f& ?                               
登录/注册后可看大图
是水平方向的流速,而东边界(图中右边)的
% Q6 ^1 ^1 Y! R4 o7 W: Y8 k
                               
登录/注册后可看大图
是在陆地边界上,所以要使

% \1 J* _5 C, F                               
登录/注册后可看大图
时变为0,不能让水质点穿出边界,而

- W; p: d( U! D) Y0 t9 h                               
登录/注册后可看大图
或者

: \1 s6 S0 c, ~- [0 T8 ~                               
登录/注册后可看大图
时则不需要给予约束,这样就不用考虑西边界,因为西边界没有
2 Z$ i6 o' u. M' e
                               
登录/注册后可看大图
的格点。对北边界的

9 E9 @; G3 @8 z( S                               
登录/注册后可看大图
处理方式一样,在此不做赘述。


/ n$ E/ U8 x  p+ K  F; m


4 c7 f9 e3 _) y- ?" ^                               
登录/注册后可看大图

un(j,k)=0;uu=u(j,k);duu=du(j,k);if(wet(j,k)==1)
6 d& B) l% q" I& x9 W* K    if((wet(j,k+1)==1)||(duu>0))5 M5 H5 c/ H7 ^! \5 a3 R
        un(j,k)=uu+duu;
; v: N7 t& H% ]9 B. m- z) e0 b    endelse# ^8 l6 C. D; v" z2 }9 b0 r; h4 w
    if((wet(j,k+1)==1)&&(duu<0))
' x; U3 r% w. L, q4 D! I: S0 q( u        un(j,k)=uu+duu;. `  p0 {- n" c& h9 }6 Z/ r
    endend

  对


9 ]' Q5 l* J6 ]% ^% l                               
登录/注册后可看大图
求解的代码实现过程如上,其中if就是在判断边界和速度方向,把不在边界和在边界但是
2 x" z# L. x, V1 g7 m2 ~
                               
登录/注册后可看大图
不大于0的值赋计算结果,使得在边界且大于0的点赋为0。
! J% [+ ~% m" s6 |0 f
                               
登录/注册后可看大图
的计算和
; V* B- {' B+ u$ G+ d: F
                               
登录/注册后可看大图

的计算同理。

: U9 S5 [% T5 V- q# |: D

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


/ G. r( V; f. b- v; O5 V/ r* P                               
登录/注册后可看大图

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


& s* j- ]0 M/ m. ]1 d4 f+ _9 `7 ?                               
登录/注册后可看大图

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

- w& t2 M) ^$ \
                               
登录/注册后可看大图

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

# h9 ~8 F* J- \+ S8 W
                               
登录/注册后可看大图

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

! B7 n/ t/ o% Y! p$ l, C" _
                               
登录/注册后可看大图

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

% l" X2 t4 M! N3 x% h- W" @
                               
登录/注册后可看大图


& H5 Z: h4 p, F* E( i8 U* c2 h% ~                               
登录/注册后可看大图

,具体表达式如下。

6 D4 }$ Z5 o& v$ X/ ^2 t/ a


% M  p2 b9 `7 S                               
登录/注册后可看大图

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

3 ?! j) S; D9 q2 c! y/ n9 L7 d- E
                               
登录/注册后可看大图
4 Z8 ^6 x7 |6 P
                               
登录/注册后可看大图
这种变量。这就是前文已经提到的那个问题,Arakawa C网格的

& m2 i- r" e7 a& Q                               
登录/注册后可看大图
1 Q/ _! T. l5 [( e; M# U$ ?; q
                               
登录/注册后可看大图
在网格的不同位置处,但是做运算只能在相同网格处。因此

. w$ k" u6 _( G; ~1 u/ x                               
登录/注册后可看大图

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

9 K. q- f8 M5 X" E. K. R
                               
登录/注册后可看大图
,而
: t8 ]+ @; K! `- b0 j8 W
                               
登录/注册后可看大图
的意思是u网格所在位置的速度
* h2 Q# i+ ~) w  {
                               
登录/注册后可看大图
。以

, C( S; k# d% g0 M6 y+ ?# h! P! w                               
登录/注册后可看大图
为例,如下图所示,可由周围的四个

$ Q  @* }- ]  k: z                               
登录/注册后可看大图
平均求得。

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

; l" D3 M" V2 b9 \) F, p, J  b
                               
登录/注册后可看大图

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


2 B$ d5 M) k+ E3 I3 r) f                               
登录/注册后可看大图

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

for k=1:N) W2 R" ~6 g5 }3 P& E- k
    h(k)=hzero(k) + eta(k);. B! k9 F, @" w7 E6 V# M9 ^4 a% {
    wet(k)=1;1 i# p2 i7 Z( j5 `- {
    if(h(k)<hmin)
( H( F4 Z& i+ x, V- U        wet(k)=0;  |' s) U6 O( t3 E6 ~2 I1 l# {
    end; r- N4 z' V! K. a. _7 z
    u(k)=un(k);end

版权声明

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

参考书目

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

! S0 o) D# k% X! _$ r
回复

举报 使用道具

相关帖子

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