各向异性VTI介质地震波场模拟，速度-应力弹性波方程交错网格有限差分法，PML吸收边界（C++和MATLAB两种语言的代码都有）

已于 2024-07-26 10:44:08 修改
阅读量488 收藏 7
点赞数 5
分类专栏： 地球物理 计算物理 各向异性理论方法 文章标签： c++ c语言 matlab 算法
于 2023-10-18 21:11:46 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42262783/article/details/133914854
版权 
参考各向同性的程序，进一步修改为VTI介质，不懂的随时交流。

这里是C/C++写的：

/*************************************************************************************/
/***********Elastic Velocity-Stress Finite-difference Modeling with PML***************/
/*************High Order (Highest: sixteen order) Finite-difference*******************/
/***************************Written By Zhang Jianming,2022.10.10**********************/
/***************************************CopyRight*************************************/
#include<stdio.h>   
#include<math.h>   
#include<iostream>   
#include<fstream>   
#include<iomanip>   
using namespace std;

#define PI 3.1415926   
#define dt 0.0005   
#define dx 8 
#define dz 8   
#define NX 401   
#define NZ 401   
#define NT 401   
#define N 1 
#define pml 25

//float Txx[NX][NZ], Tzz[NX][NZ], Txz[NX - 1][NZ - 1], Vx[NX - 1][NZ], Vz[NX][NZ - 1], L[NX][NZ], M[NX][NZ], e[NX][NZ];
//float Txx_x[NX][NZ], Tzz_x[NX][NZ], Txz_x[NX - 1][NZ - 1], Vx_x[NX - 1][NZ], Vz_x[NX][NZ - 1];
//float Txx_z[NX][NZ], Tzz_z[NX][NZ], Txz_z[NX - 1][NZ - 1], Vx_z[NX - 1][NZ], Vz_z[NX][NZ - 1];
//float C11[NX][NZ], C33[NX][NZ], C44[NX][NZ], C66[NX][NZ], C13[NX][NZ], Rou[NX][NZ];

//float Txx[NX][NZ],Txx_x[NX][NZ],Txx_z[NX][NZ];
//float Tzz[NX][NZ],Tzz_x[NX][NZ],Tzz_z[NX][NZ];
//float Txz[NX][NZ],Txz_x[NX][NZ],Txz_z[NX][NZ];
//float Vx[NX][NZ],Vx_x[NX][NZ],Vx_z[NX][NZ];
//float Vz[NX][NZ],Vz_x[NX][NZ],Vz_z[NX][NZ]; 
//float L[NX][NZ], M[NX][NZ], e[NX][NZ];
//float C11[NX][NZ], C33[NX][NZ], C44[NX][NZ], C66[NX][NZ], C13[NX][NZ], Rou[NX][NZ];

int main()
{
	float **Txx=new float*[NX],**Txx_x=new float*[NX],**Txx_z=new float*[NX];
	float **Tzz=new float*[NX],**Tzz_x=new float*[NX],**Tzz_z=new float*[NX];
	float **Txz=new float*[NX],**Txz_x=new float*[NX],**Txz_z=new float*[NX];
	float **Vx=new float*[NX],**Vx_x=new float*[NX],**Vx_z=new float*[NX];
	//A wierd question: icpc compeler broke down using the below Vz, why?
	float **Vz=new float*[NX],**Vz_x=new float*[NX],**Vz_z=new float*[NX];
	
	float **L=new float*[NX],**M=new float*[NX],**e=new float*[NX];
	float **C11=new float*[NX],**C33=new float*[NX],**C44=new float*[NX],**C66=new float*[NX],**C13=new float*[NX],**Rou=new float*[NX];
	for(int i=0;i<NX;i++)
	{
		Txx[i]=new float[NZ];Txx_x[i]=new float[NZ];Txx_z[i]=new float[NZ];
		Tzz[i]=new float[NZ];Tzz_x[i]=new float[NZ];Tzz_z[i]=new float[NZ];
		Txz[i]=new float[NZ];Txz_x[i]=new float[NZ];Txz_z[i]=new float[NZ];
		Vx[i]=new float[NZ];Vx_x[i]=new float[NZ];Vx_z[i]=new float[NZ];
		Vz[i]=new float[NZ];Vz_x[i]=new float[NZ];Vz_z[i]=new float[NZ];
		L[i]=new float[NZ];M[i]=new float[NZ];e[i]=new float[NZ];
		C11[i]=new float[NZ];C33[i]=new float[NZ];C44[i]=new float[NZ];C66[i]=new float[NZ];C13[i]=new float[NZ];Rou[i]=new float[NZ];
		for(int j=0;j<NZ;j++)
		{
			Txx[i][j]=0;Txx_x[i][j]=0;Txx_z[i][j]=0;
			Tzz[i][j]=0;Tzz_x[i][j]=0;Tzz_z[i][j]=0;
			Txz[i][j]=0;Txz_x[i][j]=0;Txz_z[i][j]=0;
			Vx[i][j]=0;Vx_x[i][j]=0;Vx_z[i][j]=0;
			Vz[i][j]=0;Vz_x[i][j]=0;Vz_z[i][j]=0;
			L[i][j]=0;M[i][j]=0;e[i][j]=0;
			C11[i][j]=0;C33[i][j]=0;C44[i][j]=0;C66[i][j]=0;C13[i][j]=0;Rou[i][j]=0;
		}
	}

	int i,j,k,f0=30;float t0=1.2/f0;
	FILE *fp1,*fp2,*fp3,*fp4,*fp5;
	float ddx[NX][NZ];float ddz[NX][NZ];
	float R=0.001,Vmax=2500;
	int x,z,plx,plz;
	for(i=0;i<NX;i++)
		for(j=0;j<NZ;j++)
		{
			Txx[i][j]=0;Txx_x[i][j]=0;Txx_z[i][j]=0;
			Tzz[i][j]=0;Tzz_x[i][j]=0;Tzz_z[i][j]=0;
			Txz[i][j]=0;Txz_x[i][j]=0;Txz_z[i][j]=0;
			Vx[i][j]=0;Vx_x[i][j]=0;Vx_z[i][j]=0;
			Vz[i][j]=0;Vz_x[i][j]=0;Vz_z[i][j]=0;
			L[i][j]=1.19*pow(10, 10);M[i][j]=5.4*pow(10, 9);e[i][j]=2.0*pow(10, 3);

			//if(j<241)
			if(j<55551)
			{
				//Iso
				//C11[i][j]=L[i][j]+2*M[i][j];
				//C33[i][j]=C11[i][j];
				//C44[i][j]=M[i][j];
				//C66[i][j]=M[i][j];
				//C13[i][j]=L[i][j];
				//Rou[i][j]=e[i][j];
				//Ice 
				//C11[i][j]=13.33*pow(10, 9);
				//C33[i][j]=14.28*pow(10, 9);
				//C44[i][j]=3.26*pow(10, 9);
				//C66[i][j]=3.65*pow(10, 9);
				//C13[i][j]=5.08*pow(10, 9);
				//Rou[i][j]=0.91*pow(10, 3);
				//Quartz 
				//C11[i][j]=116.6*pow(10, 9);
				//C33[i][j]=110.40*pow(10, 9);
				//C44[i][j]=36.06*pow(10, 9);
				//C66[i][j]=49.95*pow(10, 9);
				//C13[i][j]=32.8*pow(10, 9);
				//Rou[i][j]=2.65*pow(10, 3);
				//Mantle 
				C11[i][j]=212.8*pow(10, 9);
				C33[i][j]=249.18*pow(10, 9);
				C44[i][j]=70.47*pow(10, 9);
				C66[i][j]=66.72*pow(10, 9);
				C13[i][j]=76.24*pow(10, 9);
				Rou[i][j]=3.22*pow(10, 3);
			}
			else
			{
				C11[i][j]=116.6*pow(10, 9);
				C33[i][j]=110.40*pow(10, 9);
				C44[i][j]=36.06*pow(10, 9);
				C66[i][j]=49.95*pow(10, 9);
				C13[i][j]=32.8*pow(10, 9);
				Rou[i][j]=2.65*pow(10, 3);
			}
		}

	//PML
	plx=pml*dx;plz=pml*dz;
	for(i=0;i<NX;i++)
		for(j=0;j<NZ;j++)
		{
			if(i>=0&&i<pml&&j>=0&&j<pml)
			{
				x=pml-i;z=pml-j;
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}

			else if(i>=0&&i<pml&&j>NZ-pml&&j<NZ)
			{
				x=pml-i;z=j-(NZ-pml);
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}
			else if(i>NX-pml&&i<NX&&j>=0&&j<pml)
			{
				x=i-(NX-pml);z=pml-j;
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}
			else if(i>NX-pml&&i<NX&&j>NZ-pml&&j<NZ)
			{
				x=i-(NX-pml);z=j-(NZ-pml);
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}
			else if(i>=pml&&i<=NX-pml&&j>=0&&j<pml)
			{
				x=i-pml;z=pml-j;
				ddx[i][j]=0;
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}
			else if(i>=pml&&i<=NX-pml&&j>NZ-pml&&j<NZ)
			{
				x=i-pml;z=j-(NZ-pml);
				ddx[i][j]=0;
				ddz[i][j]=-log(R)*3*Vmax*z*z/(2*plz*plz);
			}
			else if(i>=0&&i<pml&&j>=pml&&j<=NZ-pml)
			{
				x=pml-i;z=j-pml;
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=0;
			}
			else if(i>NX-pml&&i<NX&&j>=pml&&j<=NZ-pml)
			{
				x=i-(NX-pml);z=j-pml;
				ddx[i][j]=-log(R)*3*Vmax*x*x/(2*plx*plx);
				ddz[i][j]=0;
			}
			else if(i>=pml&&i<=NX-pml&&j>=pml&&j<=NZ-pml)
			{
				x=i-pml;z=j-pml;
				ddx[i][j]=0;
				ddz[i][j]=0;
			}
		}

	for(k=0;k<NT;k++)
	{
		for(i=N;i<NX-N;i++)
			for(j=N;j<NZ-N;j++)
			{
				//Txx[i][j]=Txx[i][j]+dt*((L[i][j]+2*M[i][j])*(Vx[i][j]-Vx[i-1][j])/dx+(Vz[i][j]-Vz[i][j-1])/dz*L[i][j]);
				//Tzz[i][j]=Tzz[i][j]+dt*((L[i][j]+2*M[i][j])*(Vz[i][j]-Vz[i][j-1])/dz+(Vx[i][j]-Vx[i-1][j])/dx*L[i][j]);
				Txx_x[i][j]=(1-0.5*dt*ddx[i][j])/(1+0.5*dt*ddx[i][j])*Txx_x[i][j]+dt*(C11[i][j])/(1+0.5*dt*ddx[i][j])*(Vx[i][j]-Vx[i-1][j])/dx;
				Txx_z[i][j]=(1-0.5*dt*ddz[i][j])/(1+0.5*dt*ddz[i][j])*Txx_z[i][j]+dt*C13[i][j]/(1+0.5*dt*ddz[i][j])*(Vz[i][j]-Vz[i][j-1])/dz;
				Tzz_x[i][j]=(1-0.5*dt*ddx[i][j])/(1+0.5*dt*ddx[i][j])*Tzz_x[i][j]+dt*C13[i][j]/(1+0.5*dt*ddx[i][j])*(Vx[i][j]-Vx[i-1][j])/dx;
				Tzz_z[i][j]=(1-0.5*dt*ddz[i][j])/(1+0.5*dt*ddz[i][j])*Tzz_z[i][j]+dt*(C33[i][j])/(1+0.5*dt*ddz[i][j])*(Vz[i][j]-Vz[i][j-1])/dz;
				Txx[i][j]=Txx_x[i][j]+Txx_z[i][j];
				Tzz[i][j]=Tzz_x[i][j]+Tzz_z[i][j];
				if(i==NX/2&&j==NZ/2)
				{
					float tmp=pow(PI*f0*(k*dt-t0),2);
					//Ricker
					Txx_x[i][j]+=0.5*exp(-1.0*tmp)*(1-2.0*tmp);
					Txx_z[i][j]+=0.5*exp(-1.0*tmp)*(1-2.0*tmp);
					Tzz_x[i][j]+=0.5*exp(-1.0*tmp)*(1-2.0*tmp);
					Tzz_z[i][j]+=0.5*exp(-1.0*tmp)*(1-2.0*tmp);
					//Gauss
					//Tzz[i][j]+=-1.0*PI*PI*f0*f0*(k*dt-t0)*exp(-1.0*tmp)*(3.0-2.0*tmp);
					//Txx[i][j]+=-1.0*PI*PI*f0*f0*(k*dt-t0)*exp(-1.0*tmp)*(3.0-2.0*tmp);
				}
			}

		for(i=N;i<NX-N-1;i++)
			for(j=N;j<NZ-N-1;j++)
			{
				//Txz[i][j]=Txz[i][j]+dt*M[i][j]*((Vz[i+1][j]-Vz[i][j])/dx+(Vx[i][j+1]-Vx[i][j])/dz);
				Txz_x[i][j]=(1-0.5*dt*ddx[i][j])/(1+0.5*dt*ddx[i][j])*Txz_x[i][j]+dt*C44[i][j]/(1+0.5*dt*ddx[i][j])*(Vz[i+1][j]-Vz[i][j])/dx;
				Txz_z[i][j]=(1-0.5*dt*ddz[i][j])/(1+0.5*dt*ddz[i][j])*Txz_z[i][j]+dt*C66[i][j]/(1+0.5*dt*ddz[i][j])*(Vx[i][j+1]-Vx[i][j])/dz;
				Txz[i][j]=Txz_x[i][j]+Txz_z[i][j];
			}

		for(i=N;i<NX-N-1;i++)
			for(j=N;j<NZ-N;j++)
			{	
				//Vx[i][j]=Vx[i][j]+dt/e[i][j]*((Txx[i+1][j]-Txx[i][j])/dx+(Txz[i][j]-Txz[i][j-1])/dz);
				Vx_x[i][j]=(1-0.5*dt*ddx[i][j])/(1+0.5*dt*ddx[i][j])*Vx_x[i][j]+dt/Rou[i][j]/(1+0.5*dt*ddx[i][j])*(Txx[i+1][j]-Txx[i][j])/dx;
				Vx_z[i][j]=(1-0.5*dt*ddz[i][j])/(1+0.5*dt*ddz[i][j])*Vx_z[i][j]+dt/Rou[i][j]/(1+0.5*dt*ddz[i][j])*(Txz[i][j]-Txz[i][j-1])/dz;
				Vx[i][j]=Vx_x[i][j]+Vx_z[i][j];
			}


		for(i=N;i<NX-N;i++)
			for(j=N;j<NZ-N-1;j++)
			{
				//Vz[i][j]=Vz[i][j]+dt/e[i][j]*((Tzz[i][j+1]-Tzz[i][j])/dz+(Txz[i][j]-Txz[i-1][j])/dx);
				Vz_x[i][j]=(1-0.5*dt*ddx[i][j])/(1+0.5*dt*ddx[i][j])*Vz_x[i][j]+dt/Rou[i][j]/(1+0.5*dt*ddx[i][j])*(Txz[i][j]-Txz[i-1][j])/dx;
				Vz_z[i][j]=(1-0.5*dt*ddz[i][j])/(1+0.5*dt*ddz[i][j])*Vz_z[i][j]+dt/Rou[i][j]/(1+0.5*dt*ddz[i][j])*(Tzz[i][j+1]-Tzz[i][j])/dz;
				Vz[i][j]=Vz_x[i][j]+Vz_z[i][j];
			}
	}
	fp1=fopen("Txx.dat", "wb+");
	for (i = 0; i < NX; i++)
		for (j = 0; j < NZ; j++)
			fwrite(&Txx[i][j], sizeof(float), 1, fp1);
	fp2=fopen("Tzz.dat", "wb+");
	for (i = 0; i < NX; i++)
		for (j = 0; j < NZ; j++)
		{
			fwrite(&Tzz[i][j], sizeof(float), 1, fp2);
		}
	fp3=fopen("Txz.dat", "wb+");
	for (i = 0; i < NX-1; i++)
		for (j = 0; j < NZ-1; j++)
			fwrite(&Txz[i][j], sizeof(float), 1, fp3);
	fp4=fopen("Vx.dat", "wb+");
	for (i = 0; i < NX-1; i++)
		for (j = 0; j < NZ; j++)
			fwrite(&Vx[i][j], sizeof(float), 1, fp4);
	fp5=fopen("Vz.dat", "wb+");
	for (i = 0; i < NX; i++)
		for (j = 0; j < NZ-1; j++)
			fwrite(&Vz[i][j], sizeof(float), 1, fp5);
}

这里是MATLAB写的：

clear;clc
dt=0.002;
dx=20;
dz=20;
NX=801;
NZ=801;
NT=801;
N=1;
pml=20;
Uxx=zeros(NX,NZ);
Uzz=zeros(NX,NZ);
Uxz=zeros(NX-1,NZ-1);
Vx=zeros(NX-1,NZ);
Vz=zeros(NX,NZ-1);
L=zeros(NX,NZ);
M=zeros(NX,NZ);
e=zeros(NX,NZ);

Uxx_x=zeros(NX,NZ);
Uzz_x=zeros(NX,NZ);
Uxz_x=zeros(NX-1,NZ-1);
Vx_x=zeros(NX-1,NZ);
Vz_x=zeros(NX,NZ-1);

Uxx_z=zeros(NX,NZ);
Uzz_z=zeros(NX,NZ);
Uxz_z=zeros(NX-1,NZ-1);
Vx_z=zeros(NX-1,NZ);
Vz_z=zeros(NX,NZ-1);

C11=zeros(NX,NZ);
C33=zeros(NX,NZ);
C44=zeros(NX,NZ);
C66=zeros(NX,NZ);
C13=zeros(NX,NZ);
Den=zeros(NX,NZ);

f0=10;
t0=1.2/f0;
t=-t0:dt:(NT-1)*dt-t0;
%FILE *fp1,*fp2,*fp3,*fp4,*fp5;
ddx=zeros(NX,NZ);
ddz=zeros(NX,NZ);
R=0.000001;
Vmax=2500;
for i=1:NX
    for j=1:NZ
        Uxx(i,j)=0;
        Uxx_x(i,j)=0;
        Uxx_z(i,j)=0;
        Uzz(i,j)=0;
        Uzz_x(i,j)=0;
        Uzz_z(i,j)=0;
        L(i,j)=1.19*(10^10);
        M(i,j)=5.4*(10^9);
        e(i,j)=2.0*(10^3);
        if j<=241/401*NX
            C11(i,j)=13.33*(10^9);
            C33(i,j)=14.28*(10^9);
            C44(i,j)=3.26*(10^9);
            C66(i,j)=3.65*(10^9);
            C13(i,j)=5.08*(10^9);
            Den(i,j)=0.91*(10^3);
        else
            C11(i,j)=116.6*(10^9);
            C33(i,j)=110.40*(10^9);
            C44(i,j)=36.06*(10^9);
            C66(i,j)=49.95*(10^9);
            C13(i,j)=32.8*(10^9);
            Den(i,j)=2.65*(10^3);
        end
    end
end

for i=1:NX-1
    for j=1:NZ
        Vx(i,j)=0;
        Vx_x(i,j)=0;
        Vx_z(i,j)=0;
    end
end

for i=1:NX
    for j=1:NZ-1
        Vz(i,j)=0;
        Vz_x(i,j)=0;
        Vz_z(i,j)=0;
    end
end

for i=1:NX-1
    for j=1:NZ-1
        Uxz(i,j)=0;
        Uxz_x(i,j)=0;
        Uxz_z(i,j)=0;
    end
end

%PML

plx=pml*dx;
plz=pml*dz;

for i=1:NX
    for j=1:NZ
        if i>=1&&i<=pml&&j>=1&&j<=pml
            x=pml-i;z=pml-j;
            ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
            ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
        else
            if i>=1&&i<=pml&&j>=NZ-pml&&j<=NZ
                x=pml-i;z=j-(NZ-pml);
                ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
                ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
            else
                if i>=NX-pml&&i<=NX&&j>=1&&j<=pml
                    x=i-(NX-pml);z=pml-j;
                    ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
                    ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
                else
                    if i>=NX-pml&&i<=NX&&j>=NZ-pml&&j<=NZ
                        x=i-(NX-pml);z=j-(NZ-pml);
                        ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
                        ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
                    else
                        if i>=pml+1&&i<NX-pml&&j>=1&&j<=pml
                            x=i-pml;z=pml-j;
                            ddx(i,j)=0;
                            ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
                        else
                            if i>=pml+1&&i<NX-pml&&j>=NZ-pml&&j<=NZ
                                x=i-pml;z=j-(NZ-pml);
                                ddx(i,j)=0;
                                ddz(i,j)=-log(R)*3*Vmax*z*z/(2*plz*plz);
                            else
                                if i>=1&&i<=pml&&j>=pml+1&&j<NZ-pml
                                    x=pml-i;z=j-pml;
                                    ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
                                    ddz(i,j)=0;
                                else
                                    if i>=NX-pml&&i<=NX&&j>=pml+1&&j<NZ-pml
                                        x=i-(NX-pml);z=j-pml;
                                        ddx(i,j)=-log(R)*3*Vmax*x*x/(2*plx*plx);
                                        ddz(i,j)=0;
                                    else
                                        if i>=pml+1&&i<NX-pml&&j>=pml+1&&j<NZ-pml
                                            x=i-pml;z=j-pml;
                                            ddx(i,j)=0;
                                            ddz(i,j)=0;
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end


figure('units','normalized','position',[0,0,0.5,0.75])%设置成图框的大小
for k=1:NT
    for i=N+1:NX-N
        for j=N+1:NZ-N
            
            Uxx_x(i,j)=(1-0.5*dt*ddx(i,j))/(1+0.5*dt*ddx(i,j))*Uxx_x(i,j)+dt*(C11(i,j))/(1+0.5*dt*ddx(i,j))*(Vx(i,j)-Vx(i-1,j))/dx;
            Uxx_z(i,j)=(1-0.5*dt*ddz(i,j))/(1+0.5*dt*ddz(i,j))*Uxx_z(i,j)+dt*C13(i,j)/(1+0.5*dt*ddz(i,j))*(Vz(i,j)-Vz(i,j-1))/dz;
            Uzz_x(i,j)=(1-0.5*dt*ddx(i,j))/(1+0.5*dt*ddx(i,j))*Uzz_x(i,j)+dt*C13(i,j)/(1+0.5*dt*ddx(i,j))*(Vx(i,j)-Vx(i-1,j))/dx;
            Uzz_z(i,j)=(1-0.5*dt*ddz(i,j))/(1+0.5*dt*ddz(i,j))*Uzz_z(i,j)+dt*(C33(i,j))/(1+0.5*dt*ddz(i,j))*(Vz(i,j)-Vz(i,j-1))/dz;
            Uxx(i,j)=Uxx_x(i,j)+Uxx_z(i,j);
            Uzz(i,j)=Uzz_x(i,j)+Uzz_z(i,j);
            if i==(NX+1)/2&&j==(NZ+1)/2
                Uxx(i,j)=Uxx(i,j)+1000*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
                Uxx_x(i,j)=Uxx_x(i,j)+500*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
                Uxx_z(i,j)=Uxx_z(i,j)+500*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
                Uzz(i,j)=Uzz(i,j)+1000*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
                Uzz_x(i,j)=Uzz_x(i,j)+500*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
                Uzz_z(i,j)=Uzz_z(i,j)+500*(1-2*pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0))*exp(-pi*f0*(k*dt-t0)*pi*f0*(k*dt-t0));
            end
        end
    end
    for i=N+1:NX-N-1
        for j=N+1:NZ-N-1
            Uxz_x(i,j)=(1-0.5*dt*ddx(i,j))/(1+0.5*dt*ddx(i,j))*Uxz_x(i,j)+dt*C44(i,j)/(1+0.5*dt*ddx(i,j))*(Vz(i+1,j)-Vz(i,j))/dx;
            Uxz_z(i,j)=(1-0.5*dt*ddz(i,j))/(1+0.5*dt*ddz(i,j))*Uxz_z(i,j)+dt*C66(i,j)/(1+0.5*dt*ddz(i,j))*(Vx(i,j+1)-Vx(i,j))/dz;
            Uxz(i,j)=Uxz_x(i,j)+Uxz_z(i,j);
        end
    end
    for i=N+1:NX-N-1
        for j=N+1:NZ-N
            Vx_x(i,j)=(1-0.5*dt*ddx(i,j))/(1+0.5*dt*ddx(i,j))*Vx_x(i,j)+dt/Den(i,j)/(1+0.5*dt*ddx(i,j))*(Uxx(i+1,j)-Uxx(i,j))/dx;
            Vx_z(i,j)=(1-0.5*dt*ddz(i,j))/(1+0.5*dt*ddz(i,j))*Vx_z(i,j)+dt/Den(i,j)/(1+0.5*dt*ddz(i,j))*(Uxz(i,j)-Uxz(i,j-1))/dz;
            Vx(i,j)=Vx_x(i,j)+Vx_z(i,j);
        end
    end
    for i=N+1:NX-N
        for j=N+1:NZ-N-1
            Vz_x(i,j)=(1-0.5*dt*ddx(i,j))/(1+0.5*dt*ddx(i,j))*Vz_x(i,j)+dt/Den(i,j)/(1+0.5*dt*ddx(i,j))*(Uxz(i,j)-Uxz(i-1,j))/dx;
            Vz_z(i,j)=(1-0.5*dt*ddz(i,j))/(1+0.5*dt*ddz(i,j))*Vz_z(i,j)+dt/Den(i,j)/(1+0.5*dt*ddz(i,j))*(Uzz(i,j+1)-Uzz(i,j))/dz;
            Vz(i,j)=Vz_x(i,j)+Vz_z(i,j);
        end
    end



    %Uxx=Uxx+Uzz;
    mx=0.2*max(max(Uxx));
    mn=0.2*min(min(Uxx));
    if(mod(k+5,5)==0)
        Uxx=Uxx';
        imagesc(NZ,NX,Uxx)
        colorbar;
        colormap(flipud(gray))
        axis image           %纵、横轴采用等长刻度，且坐标框紧贴数据范围
        caxis([mn mx])
        set(gcf,'color','w');
        string=sprintf('%s%f%s','t=',k*dt ,'s');
        title(string,'fontsize',12)
        pause(0.0000000000001)
    end
end

张建明0704
关注
专栏目录
基于Hudson理论的VTI介质相速度模拟
05-15
天然裂缝储层的勘探和开发需要对地震波通过裂缝岩石的特征进行研究解释。由于其裂缝发育状况的复杂性,等效介质理论被广泛应用于裂缝岩石特征解释中。基于Hudson理论模型参数表征VTI介质P波、SV波和SH波的相速度。通过改变裂隙密度的大小,研究不同裂隙密度条件下P波、SV波和SH波的相速度随相位角的变化状况。
各向异性TTI介质地震波场模拟，速度-应力弹性波方程交错网格高阶有限差分法，PML吸收边界（C++和MATLAB两种语言的代码都有）
weixin_42262783的博客
10-18 931
逐渐改进，从各向同性到VTI各向异性，再到TTI介质。模拟效果可以在走时计算的文章里看到，有对比。不懂的欢迎随时交流。
FDTD方法中的完美匹配层（PML）研究（Matlab代码实现）
weixin_46039719的博客
05-10 858
1]尤双双,李孝坤.基于PML的FDTD法正演模拟[J].中国锰业,2018,36(06):191-194.DOI:10.14101/j.cnki.issn.1002-4336.2018.06.045.[3]李志. 基于FDTD的高阶PML算法及其在集成电路仿真中的应用[D].天津工业大学,2021.DOI:10.27357/d.cnki.gtgyu.2021.000394.但是在截断时，我们面临着在其边界上的反射问题。假设一个波向外传播，它最终会到达允许空间的边缘，这是由矩阵在程序中的尺寸决定的。
matlab绘制vti群速度,声波测井中的相速度与群速度讨论.PDF
weixin_39695672的博客
03-16 797
声波测井中的相速度与群速度讨论2012 年3 月 APPLIED GEOPHYSICS (应用地球物理) 第9 卷，第1 期声波测井中的相速度与群速度讨论//Studies on phase and group velocities from acoustic logging, APPLIED GEOPHYSICS, 2012,9(1), P. 108 ...
VTI知识点回顾（第二次修改）
bwdvzzlf415492的博客
07-03 455
VTI =Virtual Tunnel Interface 优势：支持IP组播、动态路由协议。VTI比GRE OVER IPSEC节省了GRE头部的4个字节，能承载更多的数据。 限制：只支持IP单播、组播。而GRE支持更多协议。不支持IPsec Stateful failover。可通过动态路由协议实现冗余。 VTI是IPSEC隧道技术，因此它的transform set需要调整为tunnel。（...
PDMS二次开发（一）——PML类型程序类型与概念
qq_41810415的博客
02-23 3820
PDMS二次开发需要.net 有自带的PML语言和C#.net一般通常泛指的是C#语言模型数据借助.NET的接口可以转换成数据库中的关系型数据，将固化在PDMS软件中的工程数据，转移到更广阔的应用领域。.NET二次开发的两种方式。一种是使用.NET插件；另一种是使用PML.NET。两种方式的机制都可以使用程序被动态加载到使用的模块中。我采用的是纯.NET开发的方式，熟悉PML的人可能更喜欢PML.NET的方式。编辑器的话 PML 可以用。
二维VTI介质弹性波高阶有限差分正演模拟+pml边界条件
11-17
基于matlab编写的，VTI介质弹性波方程的高阶交错网格有限差分正演模拟+pml吸收边界条件。带有注释说明，添加了pml边界代码，很适合基础新手参考学习，波场模拟的结果以动画的方式展示。同时也可以进行波场快照的输出...
VTI numerical stimulation.zip_VTI介质_matlab_数值模拟_波场_波场快照
07-15
本资料包"VTI numerical stimulation.zip"包含了一套基于Matlab的数值模拟程序，专门用于模拟VTI（各向异性纵波速度）介质中的地震波场。VTI介质是指在不同方向上具有不同速度特性的地下岩石层，这种特性在实际的...
【物理应用】旋转交错网格弹性波正演模拟实验附matlab代码.zip
04-13
1.版本：matlab2014/2019a，内含运行结果，不会运行可私信 2.领域：智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，更多内容可点击博主头像 3.内容：标题所示，对于介绍可点击主页搜索博客 4.适合人群：本科，硕士等教研学习使用 5.博客介绍：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可si信
各向同性介质弹性波高阶交错网格有限差分数值模拟
11-17
基于matlab，各向同性介质弹性波方程的高阶交错网格有限差分正演模拟+pml吸收边界条件。带有注释说明，添加了pml边界代码，很适合基础新手参考学习，波场模拟的结果以动画的方式展示。同时也可以进行波场快照的输出。注释里面还带有地震记录的代码，可以自行演示。
PML二维声波方程有限差分matlab程序
11-22
带PML条件的二维声波方程有限差分正演模拟程序，采用matlab编写。
5种pml的边界条件
09-30
FDTD计算的边界条件，比较实用，经过改编就可以自己编写FDTD的算法
PML吸收边界三维FDTD算法MATLAB代码
10-28
PML吸收边界三维FDTD算法MATLAB代码
新建文本文档 (2)_波场快照_matlab_地震有限差分_声波正演_
09-29
这篇文档主要探讨了如何利用MATLAB编程环境以及GPU加速技术来实现二维VTI（各向异性速度梯度）介质中的一阶拟声波方程的有限差分地震正演模拟。 一、波场快照 波场快照是地震模拟过程中的一个重要概念，它指的是在...
针对波动计算复杂性的吸收边界条件(PML 用于一般波动方程)（Matlab代码实现）
最新发布
weixin_46039719的博客
04-23 699
然而，即使在存在电势的情况下，如果你的初始波函数足够扩展，它最终会演变成由计算域本身的模式（通常是一个正方形）主导的一堆振荡模式。也就是说，如果你试图探索特定标量电势的模式——而不是你的计算机的模式——你需要找到一种方法，在波函数到达计算域边缘之前衰减它。对于一些电势，比如谐振子，这通常是不必要的，但对于其他电势来说是必要的。这些条件不是用来解决麦克斯韦方程，而是用于标准波动方程，比如用于电势的波动方程。哦，是的，麦克斯韦方程无穷有用，但拜托，就算是一桶激光也永远不会像需要一个波动方程是复值那样有趣。
FDTD脚本设置PML层数的代码
qq_46007076的博客
04-11 447
【代码】FDTD脚本代码。
二维声波方程的有限差分法数值模拟
Neverlevsun的博客
08-11 6131
二维声波方程的有限差分法数值模拟 文章目录二维声波方程的有限差分法数值模拟一、实现效果二、代码分享 一、实现效果 二、代码分享 close all clear clc % 此程序是有限差分法实现声波方程数值模拟 %% 参数设置 delta_t = 0.001; % s delta_s = 10; % 空间差分：delta_s = delta_x = delta_y (m) nx = 800; ny = 800; nt = 1000; fmain = 12.5; %loop:按照10000s为一次大
二维弹性波应力-位移方程
04-02
二维弹性波应力-位移方程描述了在二维弹性介质中的应力和位移之间的关系。在直角坐标系下，它可以表示为： $$\begin{aligned} \frac{\partial^2 \sigma_{xx}}{\partial x^2} + \frac{\partial^2 \sigma_{xy}}{\partial x \partial y} &= \rho \frac{\partial^2 u_x}{\partial t^2}\\ \frac{\partial^2 \sigma_{yy}}{\partial y^2} + \frac{\partial^2 \sigma_{xy}}{\partial x \partial y} &= \rho \frac{\partial^2 u_y}{\partial t^2}\\ \frac{\partial^2 \sigma_{xy}}{\partial x^2} + \frac{\partial^2 \sigma_{yy}}{\partial x \partial y} &= \rho \frac{\partial^2 u_x}{\partial t^2}\\ \frac{\partial^2 \sigma_{xy}}{\partial y^2} + \frac{\partial^2 \sigma_{xx}}{\partial x \partial y} &= \rho \frac{\partial^2 u_y}{\partial t^2}\\ \end{aligned}$$ 其中，$\sigma_{xx}$，$\sigma_{yy}$ 和 $\sigma_{xy}$ 是应力张量的分量，$u_x$ 和 $u_y$ 是位移向量的分量，$\rho$ 是介质的密度，$t$ 是时间。这些方程描述了应力和位移之间的相互作用，即介质中的应力会导致位移的变化，而位移的变化也会引起应力的变化。这些方程对于研究地震波、声波等波动现象非常重要。
评论 4
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值