利用有限差分方法合成地震记录
1、程序的编写。
program seismic
implicit none
real::lamda1,mu1,r1,dt,dx,dz,vp1,vs1,freq,lamda2,mu2,r2,vp2,vs2
real::c11,c12,c21,c22,c31,c32,d1,d2,d3,pi
real,dimension(101)::x,z
real,dimension(101,101)::uz1,uz2,uz3,u3,ux1,ux2,ux3
integer::i,j,k,n,c
character*10::text,nfile1,nfile2,tp
write(*,*)'input ytpe(mid-reflect,len-direct)'
read(*,'(a)') tp
if(tp.eq.'mid') then
print *, 'dx,dz'
read(*,*) dx,dz
else if(tp.eq.'len') then
print *, 'dx,dz'
read(*,*) dx,dz
end if
vp1=1500.0;vs1=1100.0
dt=0.00025
freq=10.0
r1=300.0
lamda1=(vp1**2-2*vs1**2)*r1
mu1=vs1**2*r1
pi=atan(1.0)*4.00
c11=0.5*(lamda1+2.0*mu1)/r1
c21=0.5*(lamda1+mu1)/r1
c31=mu1/(2.0*r1)
d1=(dt/dx)**2
d2=dt**2/(4.0*dx*dz)
d3=(dt/dz)**2
n=500
ux1=0.0;uz1=0.0
ux2=0.0;uz2=0.0
ux3=0.0;uz3=0.0
do i=2,100
do j=2,100
ux2(i,j)=ux1(i,j)+c11*d1*(ux1(i+1,j)+ux1(i-1,j)-2*ux1(i,j)) &
+c21*d2*((uz1(i+1,j+1)-uz1(i-1,j+1))-(uz1(i+1,j-1)-uz1(i-1,j-1))) &
+c31*d3*(ux1(i,j+1)+ux1(i,j-1)-2*ux1(i,j))
uz2(i,j)=uz1(i,j)+c11*d1*(uz1(i,j+1)+uz1(i,j-1)-2*uz1(i,j)) &
+c21*d2*((ux1(i+1,j+1)-ux1(i-1,j+1))-(ux1(i+1,j-1)-ux1(i-1,j-1))) &
+c31*d3*(uz1(i+1,j)+uz1(i-1,j)-2*uz1(i,j))
end do
end do
ux2(51,51)=0.0
uz2(51,51)=0.0
ux2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
ux2(52,50)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(52,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
ux2(50,52)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(50,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
ux2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
ux2(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
ux2(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
uz2(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
do k=2,n
do i=2,100
do j=2,100
ux3(i,j)=-ux1(i,j)+2*ux2(i,j)+2*c11*d1*(ux2(i+1,j)+ux2(i-1,j)-2*ux2(i,j)) &
+2*c21*d2*((uz2(i+1,j+1)-uz2(i-1,j+1))-(uz2(i+1,j-1)-uz2(i-1,j-1)))&
+2*c31*d3*(ux2(i,j+1)+ux2(i,j-1)-2*ux2(i,j))
uz3(i,j)=-uz1(i,j)+2*uz2(i,j)+2*c11*d1*(uz2(i,j+1)+uz2(i,j-1)-2*uz2(i,j)) &
+2*c21*d2*((ux2(i+1,j+1)-ux2(i-1,j+1))-(ux2(i+1,j-1)-ux2(i-1,j-1)))&
+2*c31*d3*(uz2(i+1,j)+uz2(i-1,j)-2*uz2(i,j))
u3(i,j)=sqrt(ux3(i,j)**2+uz3(i,j)**2)
end do
end do
ux3(51,51)=0.0
uz3(51,51)=0.0
ux3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
ux3(52,50)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(52,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
ux3(50,52)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(50,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
ux3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
ux3(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
ux3(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
uz3(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)
write(text,'(i3)') k
c=mod(k,10)
if(c .eq. 0) then
nfile1=text(1:3)//'ux.dat'
open(88,file=nfile1)
do i=1,101
do j=1,101
x(i)=(i-1)*dx
z(j)=(j-1)*dz
write(88,'(1x,i3,f10.4,3x,f10.4,3x,f16.4,3x,f16.4,3x,f16.4)')k,x(i),z(j),ux3(i,j),&
uz3(i,j),u3(i,j)
end do
end do
close(88)
end if
do i=2,100
do j=2,100
ux1(i,j)=ux2(i,j)
ux2(i,j)=ux3(i,j)
uz1(i,j)=uz2(i,j)
uz2(i,j)=uz3(i,j)
end do
end do
end do
stop
end program seismic
2、实例采用网格数为100*100，x方向和z方向的间隔dx=dz=1，时间间隔dt=0.00025s， 模拟地质体的纵波速度vp=1500.0m/s，横波速度为vs=1100.0m/s，震源的频率为freq=10.0 Hz，地质体的密度为r=300.0.震源放在模型区域的中心位置，其振动为爆炸震源。
通过程序计算结果如下：
当T=10*dt时间后，由于处于爆炸的初期，波刚刚向前传播了10*dt*vp的距离。如下图一所示。
3、当T=100*dt时，波则传播到100*dt*vp=100*0.00025*1500=37.5m，则上界达到100-37.5=12.5。下界达到50+37.5=87.5m。左边界到达12.5m，右边界为87.5m如下图所示。
4、当T=200*dt时，波往外传播200*dt*vp=500*0.00025*1500=75m，则波发生反射，如下图三所
