一维浅水模式 上文提到了浅水方程是由控制方程化简近似而来,形式相对简单且适合做数值算法的测试,其自身也具备实用价值。本文计划从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);
, l5 ]6 T) e- i7 B 上述的x可以理解为在整个区域等距的取了101个点,现在就需要根据初始条件给出这101个点的初始值。 hinit=10; %参考面水深10mfor i=1:N eta(i)=0; u(i)=0;endeta(49:53)=2; %480m-520m出的初始扰动h=hinit+eta; %参考面水深加波动等于水深
8 i I. @* i5 f, ^ 初始化网格后,需要对方程进行离散化。选择最简单的显式差分即可。 可以看出第二个式子的右端项上标是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;+ i- Z" S( b& Z2 P, h+ }) G6 U
本文所采用的交错网格,上文提到过,中央差分是二阶精度,比前向和后向差分的精度都高一阶,如果把 和的网格点交错开,则在计算的时候更容易实现中央差分。计算的时候已经在这么做了。对于 的计算,则更为复杂一些。下式看起来虽然复杂,但其实就是依据是正还是负,来选取和。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;end0 w: D2 z! ~4 g
这样就完成了一次迭代的计算,可以看出,在每次迭代时,un和etan所存储的是n+1时刻的计算结果,而u和eta所存储的是n时刻的计算结果。因此,在每轮迭代的最后,un和etan的值需要赋给u和eta,然后将此轮计算的结果写入到文件中。下一轮运算时,新结果将会覆盖旧结果。只需要不断迭代上述步骤,即可一步一步求解出所需时间长度的结果。 由于水面是连续的,在网格分辨率不够高时可能会出现求解出的 呈现锯齿状,因此可以在迭代求解时对做一个一阶Shapiro滤波器平滑。求解过程就变成如下所示。使用滤波器平滑之后,水面的起伏看起来会更加连续。本文所讲的是一维等深浅水方程的求解,形式相对简单。在遇到更复杂的地形时,求解难度会加大,需要考虑的问题会更多。在之后的文章中,将会探讨二维浅水方程的求解,之后会再考虑加上风应力和更复杂地形的情况。 版权声明本文创作的初衷是用于帮助数值模式的学习者。欢迎转载,转载请私信并注明作者和出处,请勿用于任何商业用途。 参考书目Ocean Modelling for Beginners. Jochen Kämpf. Springer Berlin Heidelberg, 2009. |