一维波动方程传播模拟

一维波动传播模拟

  • 引言
    各向异性介质之中机械波的传播过程是复杂的,为了简化展示说明,以下我们先对一维的均匀介质中地震波的传播过程进行理论模拟。
    波在均匀介质中的传播,其波动方程可以表示为:
    file
    u表示位移。
  • 有限差分来近似求解上述波动方程
    右边部分表示为:
    file
    左边部分表示为:
    file
    下一时刻的位移表示为:
    file
  • matlab模拟
  • 波场快照
    file
    假定左端(0km)为自由边界条件,右端(100km)处为固定边界条件(位移=0),一下为模拟出1s-30s的地震波传播图像。
  • 模拟程序
    clc;clear all;
    filename='wave1D.gif';
    dx=1;dt=0.1;tlen=3;c=4;  %初始化变量,tlen为震源持续时间,c为波传播的速度
    u1=zeros(101,1);u2=u1;u3=u1; 
    t=0;
    jj=0;   
    while (t<=50)    %模拟的最长时间为33秒
    jj=jj+1;
      for ii=2:100   
          rhs=c^2*(u2(ii+1)-2*u2(ii)+u2(ii-1))/dx^2;  %以当前时刻求方程的右边部分
          u3(ii)=dt^2*rhs+2*u2(ii)-u1(ii);  
      end
      u3(1)=u3(2);    
      u3(101)=0.0;    %右边为固定边界条件
      if(t<=tlen)
          u3(51)=(sin(pi*t/tlen)).^2;   %地震震源时间函数(3-1-8)
    end
    u1=u2;  %将当前时刻的位移变为前一个时刻的位移
    u2=u3;   %将后一个时刻的位移变为当前时刻的位移
    plot(u2);   
    axis([0, 101 -1.2 1.2]);  
     xlabel('x/km')  
     ylabel('波的振幅')  
     f = getframe(gcf);  
     imind=frame2im(f);  %将frame格式变为图像(image)的格式,imind为图像文件
    [imind,cm] = rgb2ind(imind,256);   
    %将RGB图像文件imind转换为编号图像
    if jj==1    
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.05);  
      else
          imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1);
      end
     t=t+dt;   %时间延长
       end
  • 参考自万永革著地震学导论
  • 2019/4/8修改,补end结尾