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

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

上文提到了浅水方程是由控制方程化简近似而来,形式相对简单且适合做数值算法的测试,其自身也具备实用价值。本文计划从0到1,一步一步实现出一个完整的浅水模式。为了便于说明代码实现的细节,本文考虑分成上下两篇,本篇从最简单的一维浅水方程入手。选择的算例参考了《Ocean Modelling for beginners》中的练习题。考虑到MATLAB的界面对于初学者相对友好,本文选择使用MATLAB进行讲解和说明。

假定波动周期小于惯性周期,即科氏力可忽略,忽略非线性项和摩擦力项,控制方程变成了如下形式。其中


. }8 e7 ^; t! {: ^' {                               
登录/注册后可看大图
可以理解为水面波动与参考面的距离,在水面平静时

2 i5 z* r2 _8 ^1 K7 \8 \7 {7 o- Z8 h                               
登录/注册后可看大图

为0。

- V& d* q0 m+ F, z7 n& C% r
                               
登录/注册后可看大图
为水质点的水平速度。

' w- @) H' J, C) i" i
                               
登录/注册后可看大图

该场景如下图所示,是一个水平长度为1000m,两端封闭,水深为10m。倘若初始时刻在480m到520m处给一个扰动,使得那一段水面突然升高2m,计算扰动之后整个水面会发生怎样的运动。

. O  Z+ A4 [. r& p. c5 Y! q
                               
登录/注册后可看大图

空间步长


4 N9 k3 v/ u& }1 z+ p) \! V+ g                               
登录/注册后可看大图
,为了满足CFL条件,时间步长
9 D4 r' V  p( t) Z
                               
登录/注册后可看大图
。确定好了空间和时间步长后,我们就得到了计算所需的网格。此处之所以不对t做linspace,而只取时间间隔是因为在该格式下每次计算只需要存储n和n+1时刻的

) P! z; ~$ q) K) p" H6 P; m& D$ J                               
登录/注册后可看大图

! [2 o. }% a; I; b0 i6 B# h                               
登录/注册后可看大图
,在得到相应时刻的结果后,将内存中存储的结果输出到文件中,以节省内存开支。

L=1000;dx=10;dt=0.1;N=L/dx+1;x=linspace(0,L,N);/ n0 |( A/ k' O. Y0 _9 E

上述的x可以理解为在整个区域等距的取了101个点,现在就需要根据初始条件给出这101个点的初始值。

hinit=10; %参考面水深10mfor i=1:N     eta(i)=0;     u(i)=0;endeta(49:53)=2; %480m-520m出的初始扰动h=hinit+eta; %参考面水深加波动等于水深
4 J/ `. e2 W; S, ]" a, |4 {. @) {

初始化网格后,需要对方程进行离散化。选择最简单的显式差分即可。

2 L: q, {% I% O9 f
                               
登录/注册后可看大图

可以看出第二个式子的右端项上标是n+1时刻,这是因为每次迭代计算的时候,先算


1 l4 s0 r: X% l  u, ?4 @& B$ A                               
登录/注册后可看大图
而后算

/ r2 c* E; R2 C( I( _' P! l8 O+ N: ^                               
登录/注册后可看大图
,因此在计算第二个式子的时候已经有了n+1时刻的速度

" R+ n, _: t% E# o# [  z1 p                               
登录/注册后可看大图
了。因此右端项虽然是n+1时刻,但依然是显式差分的格式。代码中wet用于表示网格是干网格还是湿网格。为了表示网格区域内是否为陆地,定义了wet这一变量。在本案例中,只有网格两端是封闭的,因此两端wet(1)和wet(N)的值是0,其余的值是1。

for k=1:N     if(wet(k)==1)         if(wet(k+1)==1)             un(k)=u(k)-dt*g*(eta(k+1)-eta(k))/dx;         else             un(k)=0;         end     endend% 上面的代码便于理解,为简洁可写成如下形式,计算效率更高% un=u-dt*g(eta(2:end)-eta(1:end-1))/dx% un(1)=0;un(end-1:end)=0;% B2 a( F* K, [/ W7 r% |; I. V, Y

本文所采用的交错网格,上文提到过,中央差分是二阶精度,比前向和后向差分的精度都高一阶,如果把

3 A$ F$ |" B4 @7 t6 t, C
                               
登录/注册后可看大图
, v+ W4 B" [/ v; r  u; h2 G
                               
登录/注册后可看大图
的网格点交错开,则在计算的时候更容易实现中央差分。计算

$ Q/ ^, H, f; r                               
登录/注册后可看大图
的时候已经在这么做了。


0 [7 X* J+ f/ S4 O4 X                               
登录/注册后可看大图

对于

! |: O5 K' S7 @/ {% N5 X( E# G
                               
登录/注册后可看大图
的计算,则更为复杂一些。下式看起来虽然复杂,但其实就是依据

8 a7 o: r2 U% l- v+ W/ E                               
登录/注册后可看大图
是正还是负,来选取
0 ~. e) p6 E. d4 @
                               
登录/注册后可看大图
( j1 ^' l* d7 S8 s$ |$ A* N1 }
                               
登录/注册后可看大图


. e$ s* Y$ Q. A' O; g. Z                               
登录/注册后可看大图

for k=2:N-1     he1=0.5*(un(k)+abs(un(k)))*h(k);     he2=0.5*(un(k)-abs(un(k)))*h(k+1);     hue=he1+he2;     hw1=0.5*(un(k-1)+abs(un(k-1)))*h(k-1);     hw2=0.5*(un(k-1)-abs(un(k-1)))*h(k);     huw=hw1+hw2;     etan(k)=eta(k)-dt*(hue-huw)/dx;end
$ T$ N6 ]1 p7 [# g% |# ]8 y3 H; i( ]

这样就完成了一次迭代的计算,可以看出,在每次迭代时,un和etan所存储的是n+1时刻的计算结果,而u和eta所存储的是n时刻的计算结果。因此,在每轮迭代的最后,un和etan的值需要赋给u和eta,然后将此轮计算的结果写入到文件中。下一轮运算时,新结果将会覆盖旧结果。只需要不断迭代上述步骤,即可一步一步求解出所需时间长度的结果。

由于水面是连续的,在网格分辨率不够高时可能会出现求解出的


, h( n* _( H; o- g                               
登录/注册后可看大图
呈现锯齿状,因此可以在迭代求解时对

* b7 x" }  M3 I+ g# q0 c2 f                               
登录/注册后可看大图
做一个一阶Shapiro滤波器平滑。求解过程就变成如下所示。使用滤波器平滑之后,水面的起伏看起来会更加连续。

; Z2 H  H6 i3 o0 e  w
                               
登录/注册后可看大图

本文所讲的是一维等深浅水方程的求解,形式相对简单。在遇到更复杂的地形时,求解难度会加大,需要考虑的问题会更多。在之后的文章中,将会探讨二维浅水方程的求解,之后会再考虑加上风应力和更复杂地形的情况。

版权声明

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

参考书目

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

回复

举报 使用道具

相关帖子

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