//有限差分时间域二阶、空间域四阶二维声波方程加吸收边界常速模拟研究 (F=30 C=2500 时好) 震源不可放在边界
//
// u(i,k,n+1)=(p*p/12)*{16[u(i+1,k,n)+u(i-1,k,n)+u(i,k+1,n)+u(i,k-1,n)]-[u(i-2,k,n)+u(i+2,k,n)+u(i,k+2,n)+u(i,k-2,n)]}
// +(2-5*p*p)*u(i,k,n)-u(i,k,n-1)+C*C*DT*DT*S(m,n,k)
//
// p=C*DT/h
//
// 左右边界条件
// for(k=0;k<=NZ-1;k++)
// { u2[1][k]=u1[1][k]+u1[2][k]-u0[2][k]+C*DT/DX*((u1[2][k]-u1[1][k])-(u0[3][k]-u0[2][k]));
// u2[0][k]=u1[0][k]+u1[1][k]-u0[1][k]+C*DT/DX*((u1[1][k]-u1[0][k])-(u0[2][k]-u0[1][k]));
// u2[NZ-2][k]=u1[NZ-2][k]+u1[NZ-3][k]-u0[NZ-3][k]+C*DT/DX*((u1[NZ-3][k]-u1[NZ-2][k])-(u0[NZ-4][k]-u0[NZ-3][k]));
// u2[NZ-1][k]=u1[NZ-1][k]+u1[NZ-2][k]-u0[NZ-2][k]+C*DT/DX*((u1[NZ-2][k]-u1[NZ-1][k])-(u0[NZ-3][k]-u0[NZ-2][k]));
// }
// 上下边界条件
// for(i=0;i<=NX-1;i++)
// { u2[i][1]=u1[i][1]+u1[i][2]-u0[i][2]+C*DT/DZ*((u1[i][2]-u1[i][1])-(u0[i][3]-u0[i][2]));
// u2[i][0]=u1[i][0]+u1[i][1]-u0[i][1]+C*DT/DZ*((u1[i][1]-u1[i][0])-(u0[i][2]-u0[i][1]));
// u2[i][NZ-2]=u1[i][NZ-2]+u1[i][NZ-3]-u0[i][NZ-3]+C*DT/DZ*((u1[i][NZ-3]-u1[i][NZ-2])-(u0[i][NZ-4]-u0[i][NZ-3]));
// u2[i][NZ-1]=u1[i][NZ-1]+u1[i][NZ-2]-u0[i][NZ-2]+C*DT/DZ*((u1[i][NZ-2]-u1[i][NZ-1])-(u0[i][NZ-3]-u0[i][NZ-2]));
// }
//
#include <STDIO.H>
#include <MATH.H>
#include <IOSTREAM.H>
#include <IOMANIP.H> //输入输出格式控制
#include <FSTREAM.H>
#define F 30.0 //雷克子波主频，单位:hz
#define N0 2 //震源起爆时间
#define NZ 51 //纵向最大采样点号
#define NX 51 //横向最大采样点号
#define DT 0.002 //时间采样间隔，单位：s
#define NT 100 //时间最大样点号
#define DZ 9 //纵向采样间隔
#define DX 9 //横向采样间隔
#define XS 25 //震源位置横坐标样点号
#define ZS 25 //震源位置纵坐标样点号
#define C 2500 //波传播速度，单位：m/s
#define T1 80 //输出t=T1*DT时的切片
double u0[NX][NZ],u1[NX][NZ],u2[NX][NZ];
double wavelet(double f,int t0,int t,double dt)
{double wav,a,b;
b=dt*(t-t0);
a=pow(3.1415,2)*pow(f,2)*pow(b,2)*(-1);
wav=(1+2*a)*exp(a);
return wav;
}
double shot(int a,int b,int e)
{double m,n,mn,s[NT];
if(a==XS&&b==ZS)
{ m=pow(C,2)*pow(DT,2);
n=DX*DZ;
s[e]=wavelet(F,N0,e,DT);
mn=m*s[e]/n;
}
else mn=0.0;
return mn;
}
void printu2()
{
int i,k;
ofstream outfile("sjxs.out");
outfile<<NZ<<ENDL; if(n="=0)" else 第三步：边界条件处理及11点差分计算波场延拓 } u1[i][k]="0.0;" for(k="0;k<=NZ-1;k++)" for(i="0;i<=NX-1;i++)" { u0[i][k]="0.0;" ZS)给出一个脉冲震源S(t)，其它点波场P="0；" (XS, 时刻dt时，在炮点S 0)="0）" j, 第二步：赋初值;初始时刻的全波场值均为零，P(i, for(n="0;n<=NT-1;n++)" b="pow(C,2)*pow(rz,2);" a="pow(C,2)*pow(rx,2);" rz="DT/DZ;" rx="DT/DX;" 第一步：速度文件的载入及相关整理（替换） rx,rz,a,b,w; double i,k,n,g; int main() void 调用成员函数close outfile.close(); success!?<<?\n?; file the cout<<?write outfile<<setw(20)<<u2[i][k]<<endl; outfile<<setiosflags(ios::left)<<setiosflags(ios::scientific); outfile<<DX<<endl; outfile<<NX<<endl;>=2)
{ g=n-1;
for(i=2;i<=NX-3;i++)
for(k=2;k<=NZ-3;k++)
{ w=(2-2.5*a-2.5*b)*u1[i][k]-u0[i][k];
if(i==XS&&k==ZS)
u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w+shot(i,k,g);
else
u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w;
}
u2[0][0]=a/12*(16*u1[1][0]-u1[2][0])+b/12*(16*u1[0][1]-u1[0][2])+(2-2.5*a-2.5*b)*u1[0][0]-u0[0][0];
u2[1][0]=a/12*(16*(u1[0][0]+u1[2][0])-u1[3][0])+b/12*(16*u1[1][1]-u1[1][2])+(2-2.5*a-2.5*b)*u1[1][0]-u0[1][0];
u2[0][1]=a/12*(16*u1[1][1]-u1[2][1])+b/12*(16*(u1[0][0]+u1[0][2])-u1[0][3])+(2-2.5*a-2.5*b)*u1[0][1]-u0[0][1];
u2[1][1]=a/12*(16*(u1[0][1]+u1[2][1])-u1[3][1])+b/12*(16*(u1[1][0]+u1[1][2])-u1[1][3])+(2-2.5*a-2.5*b)*u1[1][1]-u0[1][1];
u2[NX-1][0]=a/12*(16*u1[NX-2][0]-u1[NX-3][0])+b/12*(16*u1[NX-1][1]-u1[NX-1][2])+(2-2.5*a-2.5*b)*u1[NX-1][0]-u0[NX-1][0];
u2[NX-2][0]=a/12*(16*(u1[NX-3][0]+u1[NX-1][0])-u1[NX-4][0])+b/12*(16*u1[NX-2][1]-u1[NX-2][2])+(2-2.5*a-2.5*b)*u1[NX-2][0]-u0[NX-2][0];
u2[NX-1][1]=a/12*(16*u1[NX-2][1]-u1[NX-3][1])+b/12*(16*(u1[NX-1][0]+u1[NX-1][2])-u1[NX-1][3])+(2-2.5*a-2.5*b)*u1[NX-1][1]-u0[NX-1][1];
u2[NX-2][1]=a/12*(16*(u1[NX-3][1]+u1[NX-1][1])-u1[NX-4][1])+b/12*(16*(u1[NX-2][0]+u1[NX-2][2])-u1[NX-2][3])+(2-2.5*a-2.5*b)*u1[NX-2][1]-u0[NX-2][1];
u2[0][NZ-1]=a/12*(16*u1[1][NZ-1]-u1[2][NZ-1])+b/12*(16*u1[0][NZ-2]-u1[0][NZ-3])+(2-2.5*a-2.5*b)*u1[0][NZ-1]-u0[0][NZ-1];
u2[0][NZ-2]=a/12*(16*u1[1][NZ-2]-u1[2][NZ-2])+b/12*(16*(u1[0][NZ-3]+u1[0][NZ-1])-u1[0][NZ-4])+(2-2.5*a-2.5*b)*u1[0][NZ-2]-u0[0][NZ-2];
u2[1][NZ-1]=a/12*(16*(u1[0][NZ-1]+u1[2][NZ-1])-u1[3][NZ-1])+b/12*(16*u1[1][NZ-2]-u1[1][NZ-3])+(2-2.5*a-2.5*b)*u1[1][NZ-1]-u0[1][NZ-1];
u2[1][NZ-2]=a/12*(16*(u1[0][NZ-2]+u1[2][NZ-2])-u1[3][NZ-2])+b/12*(16*(u1[1][NZ-3]+u1[1][NZ-1])-u1[1][NZ-4])+(2-2.5*a-2.5*b)*u1[1][NZ-2]-u0[1][NZ-2];
u2[NX-1][NZ-1]=a/12*(16*u1[NX-2][NZ-1]-u1[NX-3][NZ-1])+b/12*(16*u1[NX-1][NZ-2]-u1[NX-1][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-1]-u0[NX-1][NZ-1];
u2[NX-2][NZ-1]=a/12*(16*(u1[NX-3][NZ-1]+u1[NX-1][NZ-1])-u1[NX-4][NZ-1])+b/12*(16*u1[NX-2][NZ-2]-u1[NX-2][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-1]-u0[NX-2][NZ-1];
u2[NX-1][NZ-2]=a/12*(16*u1[NX-2][NZ-2]-u1[NX-3][NZ-2])+b/12*(16*(u1[NX-1][NZ-3]+u1[NX-1][NZ-1])-u1[NX-1][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-2]-u0[NX-1][NZ-2];
u2[NX-2][NZ-2]=a/12*(16*(u1[NX-3][NZ-2]+u1[NX-1][NZ-2])-u1[NX-4][NZ-2])+b/12*(16*(u1[NX-2][NZ-3]+u1[NX-2][NZ-1])-u1[NX-2][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-2]-u0[NX-2][NZ-2];
for(i=2;i<=NX-3;i++)
u2[i][0]=a/12*(16*(u1[i-1][0]+u1[i+1][0])-(u1[i-2][0]+u1[i+2][0]))+b/12*(16*u1[i][1]-u1[i][2])+(2-2.5*a-2.5*b)*u1[i][0]-u0[i][0];
for(i=2;i<=NX-3;i++)
u2[i][NZ-1]=a/12*(16*(u1[i-1][NZ-1]+u1[i+1][NZ-1])-(u1[i-2][NZ-1]+u1[i+2][NZ-1]))+b/12*(16*u1[i][NZ-2]-u1[i][NZ-3])+(2-2.5*a-2.5*b)*u1[i][NZ-1]-u0[i][NZ-1];
for(k=2;k<=NZ-3;k++)
u2[0][k]=a/12*(16*u1[1][k]-u1[2][k])+b/12*(16*(u1[0][k-1]+u1[0][k+1])-(u1[0][k-2]+u1[0][k+2]))+(2-2.5*a-2.5*b)*u1[0][k]-u0[0][k];
for(k=2;k<=NZ-3;k++)
u2[NX-1][k]=a/12*(16*u1[NX-2][k]-u1[NX-3][k])+b/12*(16*(u1[NX-1][k-1]+u1[NX-1][k+1])-(u1[NX-1][k-2]+u1[NX-1][k+2]))+(2-2.5*a-2.5*b)*u1[NX-1][k]-u0[NX-1][k];
for(i=2;i<=NX-3;i++)
u2[i][1]=a/12*(16*(u1[i-1][1]+u1[i+1][1])-(u1[i-2][1]+u1[i+2][1]))+b/12*(16*(u1[i][0]+u1[i][2])-u1[i][3])+(2-2.5*a-2.5*b)*u1[i][1]-u0[i][1];
for(i=2;i<=NX-3;i++)
u2[i][NZ-2]=a/12*(16*(u1[i-1][NZ-2]+u1[i+1][NZ-2])-(u1[i-2][NZ-2]+u1[i+2][NZ-2]))+b/12*(16*(u1[i][NZ-3]+u1[i][N
评论0