【实验题目】声波谐振管实验
【实验仪器】
【实验记录与分析】
（长度保留到毫米，所有图标明图题、各轴代表的量和单位、实验数据点、拟合直线及其方程和相关系数）
实验室温度：T＝24℃,谐振管直径：L0=3.2cm，声速理论值：v0=346.1m/s（保留一位小数）
一、 最小谐振频率的测定
计算开管最小谐振频率的理论值f0= Hz
所测闭管的管长L＝ cm，闭管最小谐振频率的理论值f0’= Hz
分别作开管和闭管的f与n的曲线，用直线拟合求出各自的最小谐振频率的实验平均值
图1 图2
最小谐振频率的实验结果：
开管：f实=179.47 Hz，|(f实-f0)/f0|=4.02 %;闭管：f实’= Hz，|(f实’-f0’)/f0’|= %
二、 测量驻波波形
1、 开管驻波波形（压强波）（从0开始每隔2cm测一点）
开管驻波波形图：
图3
2、闭管驻波波形（压强波）（从0开始每隔2cm测一点）
闭管驻波波形图：
图4
三、管长和谐振模式
以n为横坐标，管长为纵坐标作图，并进行直线拟合，求声速：
图5
四、（选做）管中的声速 1、测量闭管中的声速（麦克风位置接近零点）
作L与Δt曲线，直线拟合，求声速。
麦克风距离波源1.5cm，波反射回来经过的路程s(cm)=2L(cm)-3(cm),以此为横坐标作图,则声速v=1/k:
图6
闭管中的声速：v= 341.3 m/s，相对误差：|(v-v0)/v0|= 1.39 %
2、 测量开管中的声速（自行设计）
图7
【总结与讨论】
本实验总体来讲较为成功。实验一中，闭管测得基频与经验公式的结果相差很大，可能是实验过程中仪器不太稳定，而且采用了较长的管长，导致测量管长的时候不太准确。实验二、三、四都是探究性质的实验，因此没有做误差分析。但是根据常识，实验二中的波形图大体正确，实验三四测出的声速也与24摄氏度时的声速（346.1m/s）相差不多。
成绩（满分30分）： 指导教师签名： 日期：
第二篇：声波方程数值模拟实验报告
声波方程数值模拟实验报告
实验要求：
1、应用声波方程作为正演模拟的波动方程；
2、将所提供震源函数离散后绘图；
3、给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图形来；
4、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；
5、对于大模型，定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。
实验目的：
1. 通过本次作业，加深对波动方程的理解，明白波动方程所代表的物理意义。
2. 通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规律。
3. 通过模拟水平层状速度模型，体会地震波在两种介质分界面的传播规律，并能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。
4. 通过模拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波形变化趋势。
需要的已知条件包括：
1）震源函数
2）地层速度（波速）
3）边界条件
2．弹性波方程：
声波方程的有限差分法数值模拟
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
（4-1）
是介质在点（x , z）处的纵波速度，
为描述速度位或者压力的波场，
为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。
设x、z方向的网格间隔长度为
，
为时间采样步长，则有：
（i为正整数）
(j为正整数)
(n为正整数)
表示在(i,j)点，k时刻的波场值。
将
在(i,j)点k时刻用Taylor展式展开：
(4-2)
将
在(i,j)点k时刻用Taylor展式展开：
(4-3)
将上两式相加，略去高阶小量，整理得(i,j)点k时刻的二阶时间微商为：
(4-4)
对于空间微分，采用四阶精度差分格式，（以X方向为例）即将
、
、
分别在(i,j)点k时刻展开到四阶小量，消除四阶小量并解出二阶微分得：
(4-5)
同理可得：
(4-6)
这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得：
}
(4-7)
式中
为介质速度的空间离散值，是空间离散步长，为时间离散步长，
为震源函数，关于一般使用一个理论的雷克型子波代替，即：
（4-8）
上式中，
为时间， 
为中心频率，一般取为20-40HZ，
为控制频带宽度的参数，一般取3-5。在实际计算过程中，需把此震源函数离散，参与波场计算。
确定震源位置。
稳定性条件：
(4-8)
这里
表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。
频散关系式： 
(4-9)
式中
为最小速度，
为Nyquist频率。一般取震源子波中的主频的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分
取8，而对于空间为四阶差分的情况则取4方能有效减少频散。
二．实验步骤
1、 应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；
2、 将所提供震源函数离散后绘图；
震源函数为雷克子波，离散绘图如下：
fm=30;
r=3;
t=0.002;
for n=1:200
w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);
end
figure(1),plot(w);
3、 对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻。
由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系
程序，图形如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.141593
#define FM 30
#define R 3
#define KN 200
#define XN 101
#define ZN 101
#define DH 5
#define DT 0.001
void main()
{
FILE *fp;
int i,j,k,m,n;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0
float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
for(k=0;k<KN;k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0
for(j=0;j<ZN;j++)
{if(i==50&&j==50)
f[i][j]=1;
else f[i][j]=0;
}
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
{ u1[i][j]=0.0;
u2[i][j]=0.0;
u3[i][j]=0.0;
u4[i][j]=0.0;
v[i][j]=2000; //速度相同表示同一介质
}
for(k=0;k<KN;k++)
{for(i=2;i<XN-2;i++)
for(j=2;j<ZN-2;j++)
{uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
{u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if(k==100)
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照,中间点
}
if((fp=fopen("wavefront.dat","w"))!=NULL)
{fprintf(fp,"%d\n",XN);
fprintf(fp,"%d\n",ZN);
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
fprintf(fp,"%f\n",u4[i][j]);
fclose(fp);
}
}
第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；
程序图形如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.141593
#define FM 30
#define R 3
#define KN 200
#define XN 101
#define ZN 101
#define DH 5
#define DT 0.001
void main()
{
FILE *fp;
int i,j,k,m,n;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0
float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
for(k=0;k<KN;k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0
for(j=0;j<ZN;j++)
{if(i==50&&j==50)
f[i][j]=1;
else f[i][j]=0;
}
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
{ u1[i][j]=0.0;
u2[i][j]=0.0;
u3[i][j]=0.0;
u4[i][j]=0.0;
v[i][j]=2000; //速度相同表示同一介质
}
for(k=0;k<KN;k++)
{for(i=2;i<XN-2;i++)
for(j=2;j<ZN-2;j++)
{uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
{u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if(k==160) //只需改动K值，即显示边界反射
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照，边界
}
if((fp=fopen("wavefront.dat","w"))!=NULL)
{fprintf(fp,"%d\n",XN);
fprintf(fp,"%d\n",ZN);
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
fprintf(fp,"%f\n",u4[i][j]);
fclose(fp);
}
}
4、 对于大模型，定义为水平层状速度模型；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律，指出哪是反射波，哪是透射波；
这时取小模型的常量，为减少频散，速度v至少为2400
程序图形如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.141593
#define FM 30
#define R 3
#define KN 200
#define XN 200
#define ZN 100
#define DH 5
#define DT 0.001
void main()
{
FILE *fp;
int i,j,k,m,n;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0
float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
for(k=0;k<KN;k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
{ u1[i][j]=0.0;
u2[i][j]=0.0;
u3[i][j]=0.0;
u4[i][j]=0.0;
f[i][j]=0.0;
if (j<=30)
v[i][j]=2400; //第一层速度为2400
else
v[i][j]=3000; //第二层速度为3000
}
f[100][10]=1; //定义f函数，在100，10为1
for(k=0;k<KN;k++)
{for(i=2;i<XN-2;i++)
for(j=2;j<ZN-2;j++)
{uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
u5[i][k]=u1[i][10]; //地震记录
}
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
{u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if(k==100)
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照
}
if((fp=fopen("wavefront.dat","w"))!=NULL)
{fprintf(fp,"%d\n",XN);
fprintf(fp,"%d\n",ZN);
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
fprintf(fp,"%f\n",u4[i][j]);
fclose(fp);
}
}
二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波？
这时取小模型的常量，为减少频散，速度v至少为2400
程序图像如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.141593
#define FM 30
#define R 3
#define KN 200
#define XN 200
#define ZN 200
#define DH 5
#define DT 0.001
void main()
{
FILE *fp;
int i,j,k,m,n;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],u5[XN][ZN]; //不能直接初值为0
float f[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
for(k=0;k<KN;k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++)
for(j=0;j<ZN;j++)
{ u1[i][j]=0.0;
u2[i][j]=0.0;
u3[i][j]=0.0;
u4[i][j]=0.0;
f[i][j]=0.0;
if (j<=30)
v[i][j]=2400; //第一层速度为2400
else
v[i][j]=3100; //第二层速度为3100
}
f[100][10]=1; //定义f函数，在100，10为1
for(k=0;k<KN;k++)
{for(i=2;i<XN-2;i++)
for(j=2;j<ZN-2;j++)
{uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];
uu2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];
u3[i][j]=2*u2[i][j]-u1[i][j]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
u5[i][k]=u1[i][10]; //地震记录
}
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
{u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if(k==100)
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照
}
if((fp=fopen("sei_record.dat","w"))!=NULL)
{fprintf(fp,"%d\n",XN/2);
fprintf(fp,"%d\n",KN);
for(i=0;i<XN;i=i+2)
for(j=0;j<KN;j++)
fprintf(fp,"%f\n",u5[i][j]);
fclose(fp);
}
}
从上向下看的两条直线为与震源同深度点的各点所有时刻的波场值，在向下与两直线紧邻的双曲线为第二层界面反射的记录。
三． 通过实验所发现的问题和认识
(1) 在用C语言写程序时，很多小问题不注意比如将i=i+2写成i+2,在m，n循环中出现i，j等等，导致无限循环，数据过大，不能出现dat文件等问题。经过这次试验对C语言和matlab的使用有了更加深刻的认识。
（2）程序成图时，图形不合适，调试数据时出现The data file length error! 按照， 对数据进行多次修改，得到较合理的图形。
（3）通过这次试验对地下介质中地震波的传播规律和声波方程有了深刻的认识。