我们知道,由质点的机械振动引起的波动被称为机械波。根据质点的振动方向与传播方向之间的关系,可以把机械波分为横波和纵波。横波,是指振动方向与传播方向垂直的一类波动。纵波是指振动方向与传播方向相同的一类波动。生活中有很多常见的波动现象,那么水波究竟是种什么类型的波动呢?
当我们向水中抛一个石子时,以石子为中心,会出现一圈一圈的涟漪,看起来水波的振动方向与传播方向相同,似乎水波是纵波;在传播过程中,水面又会上下起伏,看起来水波的振动方向与传播方向是垂直的,似乎水波又是横波,一种波动怎么可能即是纵波又是横波呢?
对于生活常见的水波,也很早就有科学家开始研究了,并给出了令人较为满意的解释。其实,水波本质上是受重力场和表面张力控制的一种复杂的波动。按波长与水体的平衡深度之间关系可分为重力短波型水波和重力长波型水波。重力短波型水波是指波长远小于水体平衡深度的一类水波,重力长波型水波是指波长远大于水体平衡深度的一类水波。接下来,我们将简要地介绍一下重力短波型和重力长波型水波的波动方程,并在文中给出两列重力长波型水波干涉过程的演示视频,文末附代码。
重力短波型水波
首先,我们不加推导地给出重力短波型水波运动的控制方程:
其中,φ是水波的速度势,g是重力加速度,式(1)是速度势的拉普拉斯方程。对于重力短波型水波,波长非常小,故振幅很小,所以式(2)中可以用z = 0时的导数值代替z = ζ 的导数,因此式(2)可修正为:
对于公式(3)这类波动方程,其通解可以表示为:
小结
综上可知,重力短波型水波是一种被限制在一个波长内振动的正常色散(频散)的表面波。
重力长波型水波
重力长波型水波的波动方程,适用于波长远远大于水深情况下的波动,其波动方程为:
模型1
水体的平衡深度为10 m,xoy平面的网格间距为10 m,x方向长度为3200 m,y方向长度1600 m。以模型的左上角为原点,两个波源的坐标分别是(1000 m, 800 m)和(2000 m, 800 m),波源的信号是雷克子波,波源的详细参数请查看代码中的具体细节。
模型2
水体的平衡深度为40 m,其他参数与模型1相同。
与模型1相比,水深变深,波的传播速度加快,因此两列波相遇,并发生干涉现象。
模型3
水体的平衡深度在x方向上发生变化,左侧为10 m,右侧为40 m,其他参数与模型1相同。
由于,水深在x方向发生变化,左侧与右侧水波的传播速度不同,在深度交界处存在速度差异,会发生“反射”现象。小结
重力长波本质上是一种传播速度受控于重力和水体平衡深度的非色散(频散)的横波,其性质与声波相似。需要注意的是在真实世界中,重力短波与重力长波都不会单独存在,水波是二者的复合波,因此表现出非常复杂的波动现象。
代码
主程序,TestInterference.m。
%% 重力长波型水波干涉演示程序%+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+%% 请将程序所在文件夹设为Matlab的当前路径,然后运行 %% 水波按波长与水体的平衡深度的关系可分为重力短波和重力长波,重力短波是指 %% 波长远小于水体平衡深度的水波,重力长波是指波长远大于水体平衡深度的水波 %% 。%% 此处通过数值模拟演示两重力长波类型的水波的干涉,并录制干涉过程的视频。%% 重力长波的波动方程的求解方法是交错网格高阶有限差分,需要指出的是在波 %% 动方程数值求解过程中,并没有加入吸收边界条件。需要在程序中加入吸收边 %% 界条件和需要重力短波求解程序的同学,可以通过以下方式联系我。%% %% 微信:Groundfrog %% 邮箱:wallace2012y@outlook.com; 523631898@qq.com %% 微信公众号:地震成像与模拟 %% %% 参考资料:%% 朗道, 栗弗席兹, 李植. 理论物理学教程:流体动力学[M]. 高等教育出版社, %% 2013. %% %% 作者:Yingwei Yan %% 日期:2020.06.13 %%+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+%% 例1,平衡深度为10 m的水体中两列重力长波型水波的干涉
space_x = 10; % xoy 平面网格的间距
h0 = 10*ones(161,321); % 水体的平衡深度,为 10 m% 扰动源采取雷克子波
fm = 0.05;delay_T = 20; dt = 0.25;nt=240;sid='ricker';maxamp=1; % 扰动源参数
s = calcswavelet(fm,delay_T,dt,nt,sid,maxamp); % 扰动源
sx = [101 201];sy = [81 81]; % 扰动源位置索引
g = 9.8; % 重力加速度
filename = 'record1.avi'; % 视频文件名
u1 = CalcLongGravityWave(h0,s,sx,sy,dt,space_x,g,filename);% 例2,平衡深度为40 m的水体中两列重力长波型水波的干涉
space_x = 10; % xoy 平面网格的间距
h0 = 40*ones(161,321); % 水体的平衡深度,为 40 m
fm = 0.05;delay_T = 20; dt =0.25;nt=220;sid='ricker';maxamp=1;
s = calcswavelet(fm,delay_T,dt,nt,sid,maxamp);
sx = [101 201];sy = [81 81];
g = 9.8;
filename = 'record2.avi';
u2 = CalcLongGravityWave(h0,s,sx,sy,dt,space_x,g,filename);% 例3,水体平衡深度在xoy平面发生变化时的两列重力长波型水波的干涉
space_x = 10; % xoy 平面网格的间距
h0 = 40*ones(161,321); % 水体的平衡深度在x = 1600 m处发生变化,左侧10 m,右侧40 m
h0(1:161,1:180) = 10;
fm = 0.05;delay_T = 10; dt =0.25;nt=200;sid='ricker';maxamp=1;
s = calcswavelet(fm,delay_T,dt,nt,sid,maxamp);
sx = [101 201];sy = [81 81];
g = 9.8;
filename = 'record3.avi';
u3 = CalcLongGravityWave(h0,s,sx,sy,dt,space_x,g,filename);
重力长波型水波的位移计算函数,CalcLongGravityWave.m。
function u = CalcLongGravityWave(h0,s,sx,sy,dt,space_x,g,filename)% Summary of this function goes here.% u = CalcGravityLongWave(h0,s,sx,sy,dt,space_x,g,filename)% Detailed explanation goes here.% The function is for solving the wave phenomenon of long gravity wave .% The function is generally adopted to analyse the wave phenomenon of water.%% IN% h0: the balanceable depth (m) of water, a 2D matrix in xoy coordinate% x-direction is controlled by column, y-direction is controlled% by row.% s: the perturbation source.% sx: the index of the perturbation source in x-direction.% sy: the index of the perturbation source in y-direction.% dt: the sampling interval (s) in time domain.% space_x: the spacing (m) of the grid.% g: the gravitational acceleration.% filename: the filename of the output video.%% OUT% u: a 3D matrix, the height of water wave of all particle in time domain.%% References:% 朗道, 栗弗席兹, 李植. 理论物理学教程:流体动力学[M]. 高等教育出版社, 2013.%% Author(s): Yan Yingwei% Copyright: 2020-2025% Revision: 1.0 Date: 6/13/2020%% Academy of Opto-Electronics, China Electronic Technology Group Corporation (AOE CETC)% 输入参数个数判断if nargin==6
g = 9.8;
filename = 'record.avi';elseif nargin==7
filename = 'record.avi';end
[M,N] = size(h0); % 获得平衡深度h0的行列维度
[X,Y] = meshgrid(0:space_x:(N-1)*space_x,0:space_x:(M-1)*space_x); % xoy平面网格坐标
P = length(s); % 扰动源的长度
u = zeros(M,N,P+2); % 时间域xoy平面的水波高度i = 2:M-1;j = 2:N-1;
dtx_square = (dt*dt)/(space_x*space_x);
videof = VideoWriter(filename); % 写入视频文件
videof.FrameRate = 8; % 设置视频帧数为8
open(videof); % 打开视频文件,等待写入for k=2:P+1
u(sy,sx,k) = s(k-1);
u(i,j,k+1) = dtx_square*g*h0(i,j).*((u(i+1,j,k)-2*u(i,j,k)+u(i-1,j,k))+...
(u(i,j+1,k)-2*u(i,j,k)+u(i,j-1,k)))+2*u(i,j,k)-u(i,j,k-1);
pcolor(X(i,j),Y(i,j),u(i,j,k));
colormap(gray);% colorbar;
shading interp; %用细致添加描影法渲染图形
axis image; % 沿每个坐标区使用相同的数据单位长度,并使坐标区框紧密围绕数据。
axis off; % 去除坐标轴
frame = getframe(gcf); % 获得帧
writeVideo(videof,frame); % 将视频帧写入end
close(videof); % 关闭写入end
PS:完整的代码例子已经整理到我的百度企业云盘中,感兴趣的同学可以从以下链接下载(有效期7天):
https://eyun.baidu.com/s/3i6dOg8P 需要指出的是,在重力长波型水波的计算程序中没有加入防止边界反射的吸收边界条件,想知道如何加入吸收边界条件的同学可以联系我哦~~~。 参考资料 朗道, 栗弗席兹, 李植. 理论物理学教程:流体动力学[M]. 高等教育出版社, 2013.