下游防洪标准演算过程.cpp

调洪演算程序

/*
  M--库容曲线插值节点数
  M--泄流曲线节点数
  K1--计算时段总数
  Z0--起调水位
  S[0]--起调水位相应库容
  S[K]--时段K的水库库容
  O[K]--时段K下泄流量
  Q[K]--对段K入库流量（设计洪水过程）
  Z，V--库容曲线的水位，库容数组
  V，q--泄流曲线的库容，流量数组
  QE--防洪控制点的安全泄流量。
  Y（K）--闸门全开时的泄流能力
  VF--防洪库容
  ZF--防洪高水位
*/
#include<fstream.h>
#include<math.h>
#define K1 45
#define M  25
#define N  35
void main()
{
	int i,j,k,l;
	double Z0,V0,QE,DV,SMAX,ZF,VF,A,C,D;
	double  S[K1],O[K1],Q[K1],Z[K1],Y[K1];
	double Z1[M],V[M],q[N],Z2[N];
	ifstream ZV("D:\\水位库容关系.txt");
	for(i=0;i<M;i++)
		ZV>>Z1[i]>>V[i];
	ZV.close;
	ifstream Zq("D:\\泄流能力曲线.txt");
	for(i=0;i<N;i++)
		Zq>>Z2[i]>>q[i];
	Zq.close;
	ifstream Qin("D:\\防洪标准洪水.txt");
	for(i=0;i<K1;i++)
		Qin>>Q[i];
	Qin.close;
	Z0=98.0;
	QE=18000;
	for(i=0;i<M-1;i++)
		if(Z0>=Z1[i]&&Z0<Z1[i+1])
			V0=V[i]+(Z0-Z1[i])*(V[i+1]-V[i])/(Z1[i+1]-Z1[i]);
	for(j=0;j<N-1;j++)
	    if(Z0>=Z2[j]&&Z0<Z2[j+1])
			Y[0]=q[j]+(Z0-Z2[j])*(q[j+1]-q[j])/(Z2[j+1]-Z2[j]);
    S[0]=V0;
	Z[0]=Z0;
	O[0]=Q[0];
		if(Q[0]>Y[0])
			O[0]=Y[0];
		else
			O[0]=Q[0];
		for(k=0;k<K1-1;k++)
		{
			if(Q[k+1]<=QE)
			{
				if(Q[k+1]<=Y[k])
				{
					O[k+1]=Q[k+1];
					Y[k+1]=Y[k];
					S[k+1]=S[k];
				}
			    else
				{
					    
						O[k+1]=O[k];
			            l=0;
						while(l<200)
						{
							S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
							C=S[k+1];
							{
								for(i=0;i<M-1;i++)
									if(C>=V[i]&&C<V[i+1])
							           A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
								if(C>V[M-1])
						            A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
					            for(j=0;j<N-1;j++)
						            if(A>=Z2[j]&&A<Z2[j+1])
					                   D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
									if(A>Z2[N-1])   
									   D=q[N-1];
							}
							if(fabs(O[k+1]-D)>0.0001)
							{
								O[k+1]=(O[k+1]+D)/2;
				            	l++;
							}
				            else
							{
				               	O[k+1]=D;
					            break;
							}					
						}	
						Y[k+1]=D;
				}
			}
			if(Q[k+1]>QE) 
			{
				
				DV=(Q[k]+Q[k+1]-O[k]-QE)*3600*6/200000000;
				C=S[k]+DV;
				for(i=0;i<M-1;i++)
					if(C>=V[i]&&C<V[i+1])
						A=Z1[i]+(C-V[i])*(Z1[i+1]-Z1[i])/(V[i+1]-V[i]);
                    if(C>V[M-1])
					    A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
				for(j=0;j<N-1;j++)
					if(A>=Z2[j]&&A<Z2[j+1])
						D=q[j]+(A-Z2[j])*(q[j+1]-q[j])/(Z2[j+1]-Z2[j]);
					if(A>Z2[N-1])   
						D=q[N-1];
				if(D>QE)
				{
					O[k+1]=QE;
					Y[k+1]=D;
					S[k+1]=C;
				}
				else
				{
					    double O1;
						O[k+1]=O[k];
			            l=0;
						while(l<200)
						{
							S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
							C=S[k+1];
							{
								for(i=0;i<M-1;i++)
									if(C>=V[i]&&C<V[i+1])
							           A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
								if(C>V[M-1])
						            A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
					            for(j=0;j<N-1;j++)
						            if(A>=Z2[j]&&A<Z2[j+1])
					                   D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
								if(A>Z2[N-1])   
									D=q[N-1];
							    if(D>QE) 
									O1=QE;
								else 
									O1=D;
							}
							if(fabs(O[k+1]-O1)>0.0001)
							{
								O[k+1]=(O[k+1]+O1)/2;
					            l++;								
							}
				            else
					           	break;
						}
						Y[k+1]=D;					    
				}
			}
		}
		SMAX=0.0;
	    for(k=0;k<K1;k++)
		{
			if(SMAX<S[k])
			{
				SMAX=S[k];
				A=k;
			}
		}
		for(k=A;k<K1-1;k++)
		{
			if(S[k]>V0)
			{
		       double O1;
	           O[k+1]=O[k];
               l=0;
		       while(l<200)
			   {
			         S[k+1]=(Q[k]+Q[k+1]-O[k]-O[k+1])*6*3600/200000000+S[k];
			         C=S[k+1];
					 {
				        for(i=0;i<M-1;i++)
				        if(C>=V[i]&&C<V[i+1])
		                   A=Z1[i]+(Z1[i+1]-Z1[i])*(C-V[i])/(V[i+1]-V[i]);
				        if(C>V[M-1])
				           A=Z1[M-1]+(C-V[M-1])*(Z1[M-1]-Z1[M-2])/(V[M-1]-V[M-2]);
				        for(j=0;j<N-1;j++)
				        if(A>=Z2[j]&&A<Z2[j+1])
				           D=q[j]+(q[j+1]-q[j])*(A-Z2[j])/(Z2[j+1]-Z2[j]);
				        if(A>Z2[N-1])   
				           D=q[N-1];
				        if(D>QE) 
				           O1=QE;
				        else 
				           O1=D;
					 }
			         if(fabs(O[k+1]-O1)>0.0001)
					 {
				        O[k+1]=(O[k+1]+O1)/2;
			        	l++;								
					 }
			         else
					 {
				        O[k+1]=(Q[k]+O1)*0.4;
			            break;
					 }
			   }
		       Y[k+1]=D;
			}
            else
			{
				O[k+1]=Q[k+1];
				S[k+1]=S[k];
				Y[k+1]=Y[k];
			}
		}
        for(k=0;k<K1;k++)
		{
			for(i=0;i<M-1;i++)
				if(S[k]>=V[i]&&S[k]<V[i+1])
					Z[k]=Z1[i]+(Z1[i+1]-Z1[i])*(S[k]-V[i])/(V[i+1]-V[i]);
			if(S[k]>V[M-1])
				Z[k]=Z1[M-1]+(Z1[M-1]-Z1[M-2])*(S[k]-V[M-1])/(V[M-1]-V[M-2]);
		}
        for(i=0;i<M-1;i++)
		if(SMAX>=V[i]&&SMAX<V[i+1])
			ZF=Z1[i]+(Z1[i+1]-Z1[i])*(SMAX-V[i])/(V[i+1]-V[i]);
		VF=SMAX-S[0];
        ofstream file1("D:\\防洪标准调洪结果.txt");
		file1<<"来水过程"<<"\t"<<"下泄能力"<<"\t"<<"下泄过程"<<"\t"<<"水库蓄水量"<<"\t"<<"水库水位"<<endl;
			for(k=0;k<K1;k++)
				file1<<Q[k]<<"\t"<<Y[k]<<"\t"<<O[k]<<"\t"<<S[k]<<"\t"<<Z[k]<<endl;
        file1<<"防洪高水位水位"<<ZF<<endl;
		file1<<"防洪库容"<<VF<<endl;
			file1.close();
}