111.cpp

一维油水两相模拟计算程序

#include<iostream.h>
#include<iomanip.h>
#include<fstream.h>
#include<math.h>
//void zhuigan(int n,double*a, double*b, double*c, double*d, double*x);
double defkro(double s);
double defkrw(double s);
double *p,*ss,*lmd,*a,*b,*c,*d;
int N=30,M;//N为空间节点个数
double xiangshen[13][3]=
{
		{0.2,0,0.85},
		{0.25,0.03,0.75},
		{0.3,0.06,0.62},
		{0.35,0.1,0.49},
		{0.4,0.14,0.31},
		{0.45,0.17,0.19},
		{0.5,0.27,0.14},
		{0.55,0.35,0.1},
		{0.6,0.42,0.07},
		{0.65,0.52,0.05},
		{0.7,0.65,0.03},
		{0.75,0.79,0.01},
		{0.8,0.9,0}
};
	//各点的渗透率sw,油相渗透率kro,水相渗透率krw,相渗曲线的含水率S,相渗曲线的相对渗透率KRO,KRW;
void main()
{
	int i;
	double uo=8,uw=0.5,K=3.5,f=0.28,swc=0.2,A=500*10000,L=30000,Qv=-0.10*1000000/86400,Tmax=365*86400,dt=5*86400,dx=L/N;
	M=sizeof(xiangshen)/sizeof *xiangshen;//M=13;
    p=new double[N];
    ss=new double[N];
    lmd=new double[N];
	a=new double[N];
    b=new double[N];
    c=new double[N];
    d=new double[N];
	for(i=0;i<N;i++)
	{
		ss[i]=swc;
	}
	ofstream fp3;
	fp3.open("一维.xls",ios::out);
	for (int j = 1; j < Tmax/dt+1; j++) //时间            
    {   
        fp3<<"____________________________"<<endl<<"时间:天"<<j*dt/86400<<endl;;
		fp3<<setw(10)<<"位置/m"<<setw(10)<<"压力/atm"<<setw(10)<<"含水饱和度"<<endl;
	    //zhuigan(N,a, b, c, d, p);
		for(i=0;i<N;i++)
		{
		    lmd[i]=K*(defkro(ss[i])/uo+defkrw(ss[i])/uw);
		}
	    a[0]=1;   c[0]=-1; d[0]=dx*Qv/A/lmd[0];
	    a[N-1]=-1;b[N-1]=1;d[N-1]=dx*Qv/A/lmd[N-2];
	    for(i=1;i<N-1;i++)
		{
			a[i]=-lmd[i-1]-lmd[i];
		    b[i]=lmd[i-1];
		    c[i]=lmd[i];
		    d[i]=0;
		}
		p[N-1]=150;//流出端压力
		p[N-2]=(d[N-1]-a[N-1]*p[N-1])/b[N-1];
		for(i=N-2;i>=2;i--)
		{
			p[i-1]=(d[i]-c[i]*p[i+1]-a[i]*p[i])/b[i];
		}

        p[0]=(d[0]-c[0]*p[1])/a[0];

		ss[0]=ss[0]+dt/f/dx*(Qv/A-K*defkrw(ss[0])/uw*(p[0]-p[1])/dx);
		ss[N-1]=ss[N-1]+dt/f/dx*(K*defkrw(ss[N-2])/uw/dx*(p[N-2]-p[N-1])-Qv/A*(K*defkrw(ss[N-1])/uw/(lmd[N-2])));//*(p[N-2]-p[N-1])/dx-Qv/A*K*K*defkrw(ss[N-2])/uw/lmd[N-2]);
		for(i=1;i<N-1;i++)
		{
			ss[i]=ss[i]+dt/f/dx/dx*(K*defkrw(ss[i])/uw*(p[i+1]-p[i])-K*defkrw(ss[i-1])/uw*(p[i]-p[i-1]));
		}
		for(i=0;i<N;i++)
		{
			fp3<<setw(10)<<(i)*dx/100<<setw(10)<<p[i]<<setw(10)<<ss[i]<<endl;
		}
    }   
	fp3.close();
	cout<<"计算结束！"<<endl;
}
double defkrw(double s)
{
	int i=0;
	for(i=0;i<M;i++)
	{
		if(s<xiangshen[i][0]) break;
	}
	if(i==M) return xiangshen[M-1][1];
	else if(i==0) return xiangshen[0][1];
	else
	return xiangshen[i][1]+(xiangshen[i][1]-xiangshen[i-1][1])/(xiangshen[i][0]-xiangshen[i-1][0])*(s-xiangshen[i][0]);
}
double defkro(double s)
{
	int i=0;
	for(i=0;i<M;i++)
	{
		if(s<xiangshen[i][0]) break;
	}
	if(i==M) return xiangshen[M-1][2];
	else if(i==0) return xiangshen[0][2];
	else
	return xiangshen[i][2]+(xiangshen[i][2]-xiangshen[i-1][2])/(xiangshen[i][0]-xiangshen[i-1][0])*(s-xiangshen[i][0]);
}

/*void zhuigan(int n,double*a, double*b, double*c, double*d, double*x)
{
	//求解AX=d;
	int i=0;
    for(i=1;i<n;i++)
	{
		a[i]=a[i]-c[i-1]/a[i-1]*b[i];
		d[i]=d[i]-d[i-1]/a[i-1]*b[i];
	}
	x[n-1]=d[n-1]/a[n-1];
    for(i=n-2;i>=0;i--)
	{
		x[i]=(d[i]-c[i]*x[i+1])/a[i];
	}
 }*/