|
一维浅水模式 上文提到了浅水方程是由控制方程化简近似而来，形式相对简单且适合做数值算法的测试，其自身也具备实用价值。本文计划从0到1，一步一步实现出一个完整的浅水模式。为了便于说明代码实现的细节，本文考虑分成上下两篇，本篇从最简单的一维浅水方程入手。选择的算例参考了《Ocean Modelling for beginners》中的练习题。考虑到MATLAB的界面对于初学者相对友好，本文选择使用MATLAB进行讲解和说明。 假定波动周期小于惯性周期，即科氏力可忽略，忽略非线性项和摩擦力项，控制方程变成了如下形式。其中 可以理解为水面波动与参考面的距离，在水面平静时为0。 为水质点的水平速度。该场景如下图所示，是一个水平长度为1000m，两端封闭，水深为10m。倘若初始时刻在480m到520m处给一个扰动，使得那一段水面突然升高2m，计算扰动之后整个水面会发生怎样的运动。 空间步长 ，为了满足CFL条件，时间步长。确定好了空间和时间步长后，我们就得到了计算所需的网格。此处之所以不对t做linspace，而只取时间间隔是因为在该格式下每次计算只需要存储n和n+1时刻的和，在得到相应时刻的结果后，将内存中存储的结果输出到文件中，以节省内存开支。L=1000;dx=10;dt=0.1;N=L/dx+1;x=linspace(0,L,N);
- b ^* _8 O0 J 上述的x可以理解为在整个区域等距的取了101个点，现在就需要根据初始条件给出这101个点的初始值。 hinit=10; %参考面水深10mfor i=1:N eta(i)=0; u(i)=0;endeta(49:53)=2; %480m-520m出的初始扰动h=hinit+eta; %参考面水深加波动等于水深& u. k* q0 ?- d" Y/ e
初始化网格后，需要对方程进行离散化。选择最简单的显式差分即可。 可以看出第二个式子的右端项上标是n+1时刻，这是因为每次迭代计算的时候，先算 而后算，因此在计算第二个式子的时候已经有了n+1时刻的速度了。因此右端项虽然是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;6 X/ p8 C. u" C; C' @
本文所采用的交错网格，上文提到过，中央差分是二阶精度，比前向和后向差分的精度都高一阶，如果把 和的网格点交错开，则在计算的时候更容易实现中央差分。计算的时候已经在这么做了。对于 的计算，则更为复杂一些。下式看起来虽然复杂，但其实就是依据是正还是负，来选取和。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* M0 S: _9 R: u4 u# m
这样就完成了一次迭代的计算，可以看出，在每次迭代时，un和etan所存储的是n+1时刻的计算结果，而u和eta所存储的是n时刻的计算结果。因此，在每轮迭代的最后，un和etan的值需要赋给u和eta，然后将此轮计算的结果写入到文件中。下一轮运算时，新结果将会覆盖旧结果。只需要不断迭代上述步骤，即可一步一步求解出所需时间长度的结果。 由于水面是连续的，在网格分辨率不够高时可能会出现求解出的 呈现锯齿状，因此可以在迭代求解时对做一个一阶Shapiro滤波器平滑。求解过程就变成如下所示。使用滤波器平滑之后，水面的起伏看起来会更加连续。本文所讲的是一维等深浅水方程的求解，形式相对简单。在遇到更复杂的地形时，求解难度会加大，需要考虑的问题会更多。在之后的文章中，将会探讨二维浅水方程的求解，之后会再考虑加上风应力和更复杂地形的情况。 版权声明本文创作的初衷是用于帮助数值模式的学习者。欢迎转载，转载请私信并注明作者和出处，请勿用于任何商业用途。 参考书目Ocean Modelling for Beginners. Jochen Kämpf. Springer Berlin Heidelberg, 2009. | 
. @6 q1 x7 n) \; V4 I
% Y8 N) z4 F( j( i
0 Q7 S; W9 O% B; O) n
; s/ Q0 i+ m. X7 P E8 c/ g