1dtwophasesimulator.cpp

1D TWOPHASE SImulator.cpp //一维油水两相驱油的数值模拟 //说明：使用达西单位制

//一维油水两相驱油的数值模拟
//说明：使用达西单位制
#include<iostream.h>
#include<fstream.h>
#include<stdlib.h>
const int SIZE = 100;
const int N = 40;   //网格数
//线性差值求相对渗透率
double Interp(double S, double S1, double S2, double K1, double K2)
{
	return K2-(K2-K1)*(S2-S)/(S2-S1);
}
//追赶法解三对角方程组
void Thomas(double x[N+1], int n, double a[N+1], double d[N+1], 
			double c[N+1], double b[N+1])
{
	int k;
	double pp[N+1], q[N+1], y[N+1];
	if(d[1]==0)
	{
		cout<<"求解方程组失败1!"<<endl;
		abort();
	}

	pp[1] = d[1];
	q[1] = c[1]/d[1];
	for(k=2; k<n; k++)
	{
		pp[k] = d[k]-a[k]*q[k-1];
		if(pp[k] == 0)
		{
			cout<<"求解方程组失败2!"<<endl;
			abort();
		}
		q[k] = c[k]/pp[k];
	}
	pp[n] = d[n]-a[n]*q[n-1];
	y[1] = b[1]/pp[1];
	for(k=2; k<=n-1; k++)
	{
		y[k] = (b[k]-a[k]*y[k-1])/pp[k];
	}
//第n个结点压力直接赋为0
	x[n] = 0;
	for(k=n-1; k>=1; k--)
	{
		x[k] = y[k]-q[k]*x[k+1];
	}
}
void main()
{
	int i, j;
	char buf[SIZE];
	fstream OutputFile, InputFile;
	InputFile.open("data.dat", ios::in);
	if(!InputFile)
	{
		cout<<"打开数据文件错误!"<<endl;
		abort();
	}
	InputFile.getline(buf, SIZE, '#');
	//读取粘度、绝对渗透率、空隙度、岩心长度、网格长度
	double MuO,MuW,K,Phi,L,Dx;
	InputFile>>MuO>>MuW>>K>>Phi>>L>>Dx;
	InputFile.getline(buf, SIZE, '#');
	InputFile.getline(buf, SIZE, '#');
	//读取相渗曲线
	double Sw[14], KrW[14], KrO[14];
	for(i=1; i<=13; i++)
	{
		InputFile>>Sw[i]>>KrW[i]>>KrO[i];
	    InputFile.getline(buf, SIZE, '#');
	}
	InputFile.getline(buf, SIZE, '#');
	//读取最大模拟时间、时间步长、产量、束缚水饱和度、初始压力、截面积
	int Tmax, Dt;
	double Qv, Swc, Pi, A;
	InputFile>>Tmax>>Dt>>Qv>>Swc>>Pi>>A;
	InputFile.close();
	double SwR[N+1];
	for(i=1; i<=N; i++)
	{
		SwR[i] = Swc;
	}
	double P[N+1];
	double CKrW, CKrO;   //差值计算得到的相对渗透率
	double LamdW[N+1], LamdO[N+1], Lamd[N+1];//流度 
	int T;//模拟时间
	T = 0;
	double a[N+1], d[N+1], c[N+1], b[N+1];
	OutputFile.open("Output.dat", ios::out);
	if(!OutputFile)
	{
		cout<<"打开输出文件错误!"<<endl;
		abort();
	}
	OutputFile<<"一维油水两相水驱油数值模拟结果"<<endl;
	while(SwR[N]<=0.2)
	{
		T = T+Dt;
		for(i=1; i<=N; i++)
		{		
			for(j=1; j<=13; j++)
			{
				if(SwR[i] > Sw[13])
				{
					cout<<"含水饱和度超出!"<<endl;
					abort();
				}	
				if(SwR[i] < Sw[j])
				{
					CKrW = Interp(SwR[i], Sw[j-1], Sw[j],
						KrW[j-1], KrW[j]);
					CKrO = Interp(SwR[i], Sw[j-1], Sw[j],
						KrO[j-1], KrO[j]);
					break;
				}
			}
			LamdW[i] = K*CKrW/MuW;
			LamdO[i] = K*CKrO/MuO;
			Lamd[i] = LamdW[i] + LamdO[i];
		}
		//求三对角方程组的系数
		for(i=2; i<=N-1; i++)
		{
			a[i] = Lamd[i-1];           //下次对角线
			d[i] = -(Lamd[i-1]+Lamd[i]);//主对角线
			c[i] = Lamd[i];             //上次对角线
			b[i] = 0.0;                   //常数项
		}
		d[1] = 1.0;
		c[1] = -1.0;
		a[N] = 1.0;
		d[N] = -1.0;
		b[1] = (Dx*Qv)/(A*Lamd[1]);
		b[N] = (Dx*Qv)/(A*Lamd[N-1]);
		//求解压力
		Thomas(P, N, a, d, c, b);
		//求含水饱和度
		for(i=2; i<=N-1; i++)
		{
			SwR[i] = SwR[i]+Dt*(LamdW[i]*(P[i+1]-P[i])
				-LamdW[i-1]*(P[i]-P[i-1]))/(Phi*Dx*Dx);
		}
		SwR[1] = SwR[1]+Dt*(Qv/A-LamdW[1]*(P[1]-P[2])/Dx)/(Phi*Dx);
		SwR[N] = SwR[N]+Dt*(LamdW[N-1]*(P[N-1]-P[N])/Dx-Qv*
			(LamdW[N-1]/Lamd[N-1])/A)/(Phi*Dx);
		//输出结果
		while(T==100 || T==200 || T==300 || T==400 || T==500)
		{
			//压力
			OutputFile<<"T = "<<T<<" s时的压力(单位：0.1MPa)随岩芯长度的分布："<<endl;
			for(i=1; i<=N; i++)
			{
				OutputFile<<P[i]<<' '<<endl;	
			}
			//饱和度
			OutputFile<<"T = "<<T<<" s时的含水饱和度随岩芯长度的分布："<<endl;
			for(i=1; i<=N; i++)
			{
				OutputFile<<SwR[i]<<' '<<endl;
			}
			break;
		}
	}
	OutputFile<<"水的突破时间为："<<T<<" s"<<endl;
	OutputFile.close();
}