|
二维浅水模式 上文讨论了浅水方程在一维情况下的求解，本文探讨浅水方程在二维情况下的求解，并考虑加上底摩擦和风应力的情况。 首先考虑最简单的二维浅水方程，形式如下。 还是同样的过程，对其做离散化，得到下式。 上文提到，为了方便实现中央差分来取得二阶精度，选择了交错网格，使得 和总是相差半个网格距离。在二维的情况，速度是和两个变量，因此我们需要考虑速度 放在网格的哪个位置。从上式的第二个式子可以看出，式子右边是 对的导数，因此为了让对的中央差分落在和同样的位置上，也需要和间隔半个网格距。这就是大名鼎鼎的Arakawa C网格。在使用交错网格的时候，我们要明确一个准则，就是要确保计算时各个变量都要在相同的空间位置上。比如要算在t+1时刻的值，就需要两侧的做中央差分，这样确保等号左边的和等号右边的导数是相同的空间位置。这一点在之后讲解更复杂的模式中会体现的更明显。例如POM模式，会有更多的参变量，为了保证处于同一空间位置上需要做平均运算。而本文的后半部分，在考虑底摩擦的浅水方程中也涉及到了。
& i# c! v6 F/ \6 Q' ^ 使用了Arakawa C网格后，可以神奇的发现，在这种情况下，只需要讨论两个边界的速度即可，如下图所示。考虑到 是水平方向的流速，而东边界（图中右边）的是在陆地边界上，所以要使 1 k& d9 l5 _0 N9 o, J/ e1 O; y- F( B: [
登录/注册后可看大图
6 d- z) @& i) w' F8 p3 J8 y7 |un(j,k)=0;uu=u(j,k);duu=du(j,k);if(wet(j,k)==1)
) Y9 d. p% a8 p- R$ g if((wet(j,k+1)==1)||(duu>0))
6 |0 F3 b# t# K6 Z un(j,k)=uu+duu;
6 K+ n0 z. ]5 E7 w endelse
8 o# O4 x3 e( k& R if((wet(j,k+1)==1)&&(duu<0))
' c0 g3 U, E- e7 X ?/ W2 } x# m un(j,k)=uu+duu;- F: R" l2 E/ Z. P) n1 O" L& m5 e$ \
endend 对 求解的代码实现过程如上，其中if就是在判断边界和速度方向，把不在边界和在边界但是不大于0的值赋计算结果，使得在边界且大于0的点赋为0。的计算和 的计算同理。
8 t& H; k2 B, M4 s1 ?9 Z! ^; s 接下来，开始考虑更复杂的形式，引入风应力和底摩擦。方程组的形式立马变得更加复杂起来。 底摩擦和风应力这两项由上式可以看出，是相减的关系，因此可以在计算时，分别看成两项来对待。如果只考虑底摩擦，不考虑其他项的话，方程是如下形式。 对其进行半隐式差分得到如下形式。显式差分应用于底摩擦会出现数值稳定性问题。 考虑完了底摩擦，还需要考虑风应力问题。在原方程去掉底摩擦项后再进行差分得到下式，其中下标j和k分别代表x和y方向。 将其求完后回代到推的底摩擦差分方程中即可。 为了表达式清晰简洁，引入了 和 ，具体表达式如下。 " a% N4 b1 @. F. [6 D7 s
值得注意的是，上式中出现了 和这种变量。这就是前文已经提到的那个问题，Arakawa C网格的和在网格的不同位置处，但是做运算只能在相同网格处。因此 的意思是v网格所在位置的速度 ，而的意思是u网格所在位置的速度。以为例，如下图所示，可由周围的四个平均求得。u_v = 0.25*(v(j+1,k)+v(j,k)+v(j,k-1)+v(j+1,k-1)) 这样就实现了二维考虑底摩擦和风应力的浅水方程求解。之前说过，干网格和湿网格分别对应网格中的陆地和水，我们在计算时只需要计算流体部分，而不计算陆地的部分。之前我们提到的案例中，陆地边界都在网格的四周。倘若需要模拟的场景是大洋中的一个小岛，即陆地出现在网格内部时，会有什么不同？倘若大洋中的小岛海拔高度有限，在模拟台风时，小岛被水淹没时，干网格和湿网格该怎么转换？当考虑了干湿网格的转换之后，一个具备基本功能的浅水模式就完成了。 以一维的情况为例，从图中可以看出来，这个扰动显然会在某些时刻漫过这个岛屿。为了解决这个问题，我们只需要在每轮时间步迭代的最后，加上干湿网格的更新即可。需要定义一个比较小的hmin，意思是当陆地上的积水高度超过hmin就意味着被淹没，下次计算时就把这个网格点当做水面，即把wet(k)这一点赋值为1。 for k=1:N0 X. J* g& p5 x4 q* `" J/ [
h(k)=hzero(k) + eta(k);/ }) h: C) p( B# \* @0 D8 x
wet(k)=1;
5 l( W' H5 e5 j7 j( l+ l) P if(h(k)<hmin)6 g2 m9 n1 X# @, F
wet(k)=0;
$ E+ a Q# O$ C# L8 C end
# j1 n& P3 X% e, D" | u(k)=un(k);end 版权声明 本文创作的初衷是用于帮助数值模式的学习者。欢迎转载，转载请私信并注明作者和出处，请勿用于任何商业用途。 参考书目Ocean Modelling for Beginners. Jochen Kämpf. Springer Berlin Heidelberg, 2009.
. t# D$ @: g( j( m1 e+ x. R | 
$ T: w0 K- b2 M( `! N7 z
( j9 p! s% L5 q$ x( t
* a; z, {9 z9 W# A- N3 i3 I
! j) Z4 v8 n3 H; w- M* @' |
y& u4 t' V$ o' G4 @
1 Z% w9 D" v9 O/ n8 h. w" }
1 C( R+ s5 W* s/ a) E5 z2 h$ R* U
- `* G' H. l' M, X: ?