cpp1.cpp

油藏数值模拟的cchengxu

/*定流压条件下的计算过程*/


#include <stdio.h>
#include <math.h>
/********初始参数**********/
#define K 0.1
#define p 0.2
#define m 0.5
#define C 0.0005/*单位：1/atm*/
#define h 500./*厘米*/
#define Nx 12
#define Ny 12
#define dx 5000.
#define dy 5000.
#define dt 864000
#define Pwf 120.
int tmax=180*86400;
int u=tmax/dt;
#define Rw 10.
double vp=p*dx*dy*h;
/******主函数****/
void main()
{int k,i,j,n,d;
double qv[Nx*Ny],b[Nx*Ny],W[18];
double T[Nx][Ny],a[Nx*Ny][Nx*Ny],Q[18][Nx*Ny];
double l,ll,s;
double P[18][Nx*Ny];

	/****定义传导系数******/
for(n=1;n<=u;n++)
{
for(i=0;i<=Nx;i++)
	for(j=0;j<=Ny;j++)
	{
		if (j==0||j==Ny||i==0||i==Nx)
		T[i][j]=0;
		else T[i][j]=K*h/m;
		//printf("T[%d][%d]=%f\n",i,j,T[i][j]);
	}
/*定义系数矩阵*/
for(i=0;i<=Nx*Ny;i++)
	for(j=0;j<=Nx*Ny;j++)
	{if (j==0||j==Ny||i==0||i==Nx)
	a[i][j]=0;
	}
for(j=1;j<=Ny-1;j++)
	for(i=1;i<=Nx-1;i++)
	{ 
	k=i+(j-1)*(Nx-1);
	for(d=0;d<=(Nx-1)*(Ny-1);d++)
		if((d!=k+1)&&(d!=k-1)&&(d!=k+Nx-1)&&(k!=k-Nx+1))
			a[k][d]=0;
		else	 
		{
		if(T[i][j]*T[i-1][j]==0) a[k][k-1]=0;
		else
		a[k][k-1]=2*T[i-1][j]*T[i][j]/(T[i][j]+T[i-1][j]);
		if(T[i][j]*T[i+1][j]==0) a[k][k+1]=0;
		else 
		a[k][k+1]=2*T[i+1][j]*T[i][j]/(T[i][j]+T[i+1][j]);
		if(T[i][j-1]*T[i][j]!=0) 
		a[k][k-Nx+1]=2*T[i][j-1]*T[i][j]/(T[i][j]+T[i][j-1]); 
		else
		a[k][k-Nx+1]=0;
		if(T[i][j+1]*T[i][j]==0) a[k][k+Nx-1]=0;
		else	
		a[k][k+Nx-1]=2*T[i][j]*T[i][j+1]/(T[i][j]+T[i][j+1]); 
		}
	a[k][k]=-(a[k][k+1]+a[k][k-1]+a[k][k-Nx+1]+a[k][k+Nx-1]+1.0*vp*C/dt);
	
//printf("a[%d][k-1]=%f,a[k][k+1]=%f,a[k][k-Nx]=%f,a[k][k+Nx]=%f,a[k][k]=%f\n",k,a[k][k-1],a[k][k+1],a[k][k-Nx+1],a[k][k+Nx-1],a[k][k]);
}
	/******输出矩阵A*****/
	/*for(j=1;j<=(Nx-1)*(Ny-1);j++)
		for(i=1;i<=(Nx-1)*(Ny-1);i++)
			{if(a[j][i]!=0)
			{printf("a[%d][%d]=%-10.5f",j,i,a[j][i]);
			if(i%((Ny-1)*(Nx-1))==0) printf("\n");}
			}*/


	/**********计算结果矩阵**********/
	for(j=1;j<=Ny-1;j++)
		for(i=1;i<=Nx-1;i++)
		{
			k=i+(j-1)*(Nx-1);P[0][k]=150.;
			if(i==Nx/2&&j==Ny/2)/*取中心网格*/
			
{Q[n][k]=-2*3.1415926*T[i][j]*(P[n-1][k]-Pwf)/(log(0.208*dx/Rw)-0.75)/100;
			
qv[k]=Q[n][k];}/*若为定流压条件Q[n]=2*3.1415926*T[i][j]*(P[n][k]-Pwf)/(log((0.208*dx/Rw)-0.75)*/
			else qv[k]=0;
			b[k]=-(qv[k]+vp*C*P[n-1][k]/dt);
		//printf("qv[%d]=%f,b[%d]=%f\n",k,qv[k],k,b[k]);
			

		}
	/***************矩阵求解**********************/
				/***主元消去*****/
	for(i=1;i<=(Nx-1)*(Ny-1);i++)
	{	l=a[i][i];
		b[i]=b[i]/l;
		for(j=1;j<=(Nx-1)*(Ny-1);j++)
			a[i][j]=a[i][j]/l;
		for(k=i+1;k<=(Nx-1)*(Ny-1);k++)
			{ll=a[k][i];
			for(j=1;j<=(Nx-1)*(Ny-1);++j)
			a[k][j]=a[k][j]-a[i][j]*ll;
			/*printf("a[%d][%d]=%f,b[%d]=%f\n",k,j,a[k][j],k,b[k]);*/
			b[k]=b[k]-b[i]*ll;//printf("b[%d]=%f\n",k,b[k]);
			}
	}		/***主元消去*****/
	/****输出此时的矩阵****/				
/*for(i=1;i<=(Nx-1)*(Ny-1);i++)
{
	for(j=1;j<=(Nx-1)*(Ny-1);j++)
	{
	printf("a[%d][%d]=%10.8f,",i,j,a[i][j]);
	}
	printf("b[%d]=%10.4f,",i,b[i]);
}*/
			/*******回归过程*****/
	for(i=(Nx-1)*(Ny-1);i>=1;i--)
		{
		s=0;
			for(j=i+1;j<=(Nx-1)*(Ny-1);j++)
			{s=s+a[i][j]*P[n][j];}
		P[n][i]=b[i]-s;
		printf("P[%d][%d]=%f,\n",n,i,P[n][i]);
		}			
}
				/*输出压力值*/
for(n=0;n<=u;n++)
{s=0;
for(i=1;i<=(Nx-1)*(Ny-1);i++)
s=s+P[n][i];
W[n]=s/121;

//printf("%10f\n",W[n]);  }
}
}