水利水能规划—水电站水库特征参数选择——C++程序

// 水利水能规划课设项目.cpp #include "stdafx.h"#include <fstream>#include <iostream>#include <iomanip>#include <cmath>const int  NumSeries = 20;//泄洪插值系列原始值数目double ZSeries[NumSeries], VSeries[NumSeries], qSeries[NumSeries];//插值原始系列值double ZVqChaZhi(double x[NumSeries], double y[NumSeries], double x0){//Z、V、q两两插值函数( 注意：返回q为单孔泄流流量)double y0 = -1;//返回为负值时，可以作为判定依据for(int i = 0; i < NumSeries - 1; i++){if(((x[i]>=x0)&&(x[i+1]<=x0))||((x[i]<=x0)&&(x[i+1]>=x0))) y0 = y[i] + (y[i+1] - y[i])/(x[i+1] - x[i])*(x0 -x[i]);}return y0;}#include "确定水库死水位.h"#include "设计洪水过程线.h"#include "调洪计算.h"#include "水能计算.h"#include "由工作容量求保证电能.h"void main(){using namespace std;ifstream infile;infile.open("infile_ZVqSeries.txt");for(int i = 0; i  < NumSeries; i++)//输入插值系列infile>>ZSeries[i]>>VSeries[i]>>qSeries[i];infile.close();//YuNiV();//计算淤泥体积WProcess();TiaoHong();//ChuLiJS();//NandERelation();}

//确定水库死水位.h//由历年月流量资料求年均流量系列，及淤泥体积const int StartYear = 1951,//Y = 50;//历年逐月入库流量统计年数double rho = 0.237,//多年平均含沙量(kg/m3)m_down = 0.9,//泥沙沉积率gamma = 1650,//泥沙干容重p_empty = 0.3,//孔隙率alpha = 0.15,//推移质与悬移质淤积量之比MonthQ[Y][12],//历年逐月入库流量统计资料AverageYearQ[Y] = {0},//历年逐月入库流量的全年平均ZS,//死水位YuNiVolume = 0;//设计运行年限淤泥累积体积,单位10^6(m3)void YuNiV(){using namespace std;ofstream outfile;outfile.open("outfile_AverageYearQ.txt");ifstream infile;infile.open("infile_MonthQ.txt");for(int i = 0; i < Y; i++){outfile<<setw(10)<<i+StartYear;for(int j = 0; j < 12; j++){infile>>MonthQ[i][j];YuNiVolume += MonthQ[i][j];AverageYearQ[i] += MonthQ[i][j];}AverageYearQ[i] /= 12;outfile<<setw(10)<<AverageYearQ[i]<<endl;}YuNiVolume *= ((1+alpha)*rho*2.63*m_down/(1-p_empty)/gamma);cout<<"水库淤积体积为："<<YuNiVolume<<endl;}

//设计洪水过程线.h//最小计时单位为3个小时，洪量单位统一(m3/s*3h)//W—洪量，T—3小时，J—极值，C—Classical，SJ—设计const int SumT = 52,//典型洪水时段数J = 3,//同频率法选择的时段NumJ[J] = {1, 8, 24},//同频率法选择的极值时段,WProT = 24;//设计洪水过程线总历时,即NumJ[J - 1]//NumJ[J]分别表示洪峰、最大一天、最大三天3个极值时段//三种设计标准,即下游防洪标准洪水、水库设计标准洪水、校核标准洪水double SJWPro[3][SumT],//设计洪水过程总历时由最长极值时段决定//SJWPro[3][NumJ[J - 1]]不行CTW[SumT],//典型洪水过程MaxSJW[3][J],//极值洪量的设计值,3是指三种标准MaxCW[J],//极值洪量的典型值RatioW[J],//不同极值的同频率倍比temp_W;//中间值int date[J];//记录典型极值时段，判断是否为包含关系void WProcess(){using namespace std;ifstream infile;infile.open("infile_ClassicalWProcess.txt");for(int i = 0; i < SumT; i++)infile>>CTW[i];infile.close();infile.open("infile_MaxSJW.txt");for(int k = 0; k<J; k++)for(int BZ = 0; BZ < 3; BZ++)infile>>MaxSJW[BZ][k];infile.close();//人为判断典型洪水选取是否合理cout<<"典型洪水极值："<<endl;for(int k = 0; k<J; k++){//求典型洪水的极值洪量 MaxCW[k] = 0;//初始化 for(int j = 0; j < SumT - NumJ[k] + 1; j++) { temp_W = 0;//初始化 for(int n = 0; n < NumJ[k]; n++) temp_W += CTW[j + n]; if(temp_W >MaxCW[k]) {MaxCW[k] = temp_W;//(m3/s*月)date[k] = j; } } if(k==0) cout<<"洪峰时间："<<3*date[k]; else cout<<"最大"<<NumJ[k]<<"个时段   时间范围："<<3*date[k]<<"至"<<3*(date[k]+NumJ[k] -1); cout<<setw(10)<<"洪量："<<MaxCW[k]<<endl;}cout<<"由以上时段包含关系判断典型洪水选取是否合适！"<<"合适请输入1，否则请关闭！"<<endl;cin>>temp_W;if(temp_W == 1){for(int BZ = 0; BZ < 3; BZ++){//三种标准for(int k = 0; k < J; k++){//求设计洪水过程if( k == 0||k==1) //洪峰和最大一天RatioW[k] = MaxSJW[BZ][k]/MaxCW[k];else RatioW[k] = (MaxSJW[BZ][k] - MaxSJW[BZ][k-1])/(MaxCW[k] - MaxCW[k-1]);for(int i = date[k]; i < date[k]+NumJ[k]; i++){if( k == 0) SJWPro[BZ][i] = CTW[i]*RatioW[k];else if(i < date[k-1]|| i>=date[k-1]+NumJ[k-1]) SJWPro[BZ][i] = CTW[i]*RatioW[k];}}}ofstream outfile;//输出设计洪水过程outfile.open("outfile_WProcess.txt");outfile<<"设计洪水过程："<<endl<<setw(10)<<"时间"<<setw(10)<<"典型洪水"<<setw(10)<<"下游防洪"<<setw(10)<<"水库设计"<<setw(10)<<"校核标准"<<endl;for(int i = date[J - 1]; i < date[J - 1] + WProT; i++){outfile<<setw(10)<<3*i<<setw(10)<<CTW[i];for(int BZ = 0; BZ < 3; BZ++)outfile<<setw(10)<<SJWPro[BZ][i];outfile<<endl;}outfile.close();cout<<"设计洪水过程结束,输入1继续！"<<endl;cin>>temp_W;//控制台暂停}}

//调洪计算.h//此步骤中，流量单位：m3/s//水库水量均换算成了10^6m3//用插值求水库泄水量q[i]时，一定要注意插值结果是单孔泄流//水量平衡算V[i]时一定要乘以流量转换洪量的算子const double  QToW = 3.0*60.0*60.0/pow(10.0, 6),//流量转换洪量的算子ZX = 199.3,//限制洪水位SafeQ = 10043,//下游防洪安全泄流流量epsilon_q = 0.1,//允许误差QFaDian = 1000,//发电引用流量QBottomHole = 1280;//泄洪底孔泄流流量,水位达到设计洪水位时启用const int  NumHole = 13;//泄洪道数目double q[SumT],//水库泄流流量V[SumT],//时段初始水库蓄水量Z[SumT],//水库水位///q_ZX,//限制水位对应的单孔泄流流量q_ZX,//限制水位对应的水库泄流流量ZFH,//防洪高水位ZSJ,//设计洪水位ZJH,//校核洪水位t_qMax, Q_qMax, q_qMax, V_qMax, Z_qMax,//最大泄水量对应的值q1, q2;//时段初、始单孔泄流流量int DestroyZFH = 0,//下游防洪安全还没破坏ShouWei = 0;//收尾了double qFunction(int BZ, int i,double QStart, double QEnd, double VStart, double qStart, int t){//求时段下泄流量过程//i——计算第i个泄流量//t 是时段二分割的次数,在求最大泄流时需要将时段细化,其他时段均为0int time = 0;using namespace std;double VEnd, qEnd = 0;//迭代结束判定条件q1 = qStart;//初始化while(qEnd == 0){VEnd =  VStart + pow(0.5, t)*QToW*(QStart + QEnd - qStart - q1)/2;q2 = ZVqChaZhi(VSeries, qSeries, VEnd);//由水量平衡求V[i]if(q2 < 0)//由V[i]插值求单孔泄流{cout<<"库存水量"<<VEnd<<endl//1455.94超了插值域限<<"失败！！！求BZ = "<<BZ<<"的q["<<3*i<<"]时插值超过了域限!";cin>>q2;}q2 = QFaDian + NumHole*q2;//由单孔泄流求水库泄流if(Z[i - 1]>ZSJ) q2 +=1280;//求校核洪水时打开泄洪底孔if(fabs(q2 - q1) < epsilon_q) qEnd = q1;else{q1 = (q1 + q2)/2;time++;if(time == 2000) cout<<"Time" ;}}return qEnd;}void qMax(int BZ, int i, double QStart, double QEnd, double VStart, double qStart){//求最大泄洪量using namespace std;double qEnd, deltaT = 0;int t = 1;//t 是时段二分割的次数,在求最大泄流时需要将时段细化,其他时段均为0    qEnd = qFunction(BZ, i, QStart, (QStart+QEnd)/2, VStart, qStart, t);     if(DestroyZFH ==0&&qEnd > SafeQ) qEnd = SafeQ;while(fabs(qEnd - (QStart+QEnd)/2) >= epsilon_q){if(t>100){//分割了100次还找不到,说明程序或数据有问题！cout<<"未找到BZ = "<<BZ<<"的最大流量！"<<endl<<"时段初来流"<<SJWPro[BZ][i - 1]<<endl<<"时段初泄流"<<q[i - 1]<<endl<<"时段初水位"<<Z[i - 1]<<endl<<"时段末来流"<<SJWPro[BZ][i]<<endl<<"时段末泄流"<<q[i]<<endl<<"时段末水位"<<Z[i]<<endl<<"得到的最大泄流量近似值为："<<qEnd<<endl<<"对应的入库流量为："<<(QStart+QEnd)/2<<endl; cin>>t;}if(qEnd> (QStart+QEnd)/2) QEnd = (QStart+QEnd)/2;else{deltaT += pow(0.5, t);VStart += pow(0.5, t)*QToW*(QStart + (QStart+QEnd)/2 - qEnd - qStart); qStart = qEnd; QStart =  (QStart+QEnd)/2;}t++;qEnd = qFunction(BZ, i, QStart, (QStart+QEnd)/2, VStart, qStart, t);if(DestroyZFH ==0&&qEnd > SafeQ) qEnd = SafeQ;}deltaT += pow(0.5, t);VStart += pow(0.5, t)*QToW*(QStart + (QStart+QEnd)/2 - qEnd - qStart);t_qMax = 3*(i - 1+deltaT);Q_qMax = (QStart+QEnd)/2;q_qMax = qEnd;V_qMax = VStart;Z_qMax = ZVqChaZhi(VSeries, ZSeries, V_qMax);if(BZ==0) {ZFH = Z_qMax; cout<<"防洪高水位ZFH ="<<ZFH<<endl;}if(BZ==1) {ZSJ = Z_qMax; cout<<"设计洪水位ZSJ ="<<ZSJ<<endl;}if(BZ==2) {ZJH = Z_qMax;cout<<"校核洪水位ZJH ="<<ZJH<<endl;}}void TiaoHong(){//求调洪过程using namespace std;system("cls");ifstream infile;infile.open("infile_XiuYunWProcess.txt");for(int i = date[J - 1]; i < date[J - 1] + WProT; i++){//输入修匀后的设计洪水过程for(int BZ = 0; BZ < 3; BZ++)infile>>SJWPro[BZ][i];}q_ZX =QFaDian + NumHole*ZVqChaZhi(ZSeries, qSeries, ZX); ZFH = ZSJ = 100000;//初始大化防洪高水位和设计洪水位,方便后面判断infile.close();ofstream outfile ;outfile.open("outfile_TiaoHong.txt");outfile<<"限制水位："<<ZX<<endl//199.3<<"限制水位对应的泄流流量："<<q_ZX<<endl//6502.9<<"对应的库容："<<ZVqChaZhi(ZSeries, VSeries, ZX)<<endl;for(int BZ =0; BZ < 3; BZ++){//默认初始状态为限制水位对应的泄流流量,和库存水量outfile<<endl<<endl<<endl;//(3)if(BZ==0) outfile<<"防洪高水位过程：";if(BZ==1) outfile<<"设计洪水位过程：";if(BZ==2) outfile<<"校核水位过程：";outfile<<endl<<setw(10)<<"时间"<<setw(10)<<"入库洪量" <<setw(10)<<"下泄流量"<<setw(10)<<"水库蓄水"<<setw(10)<<"水库水位"<<endl;SJWPro[BZ][date[J - 1] - 1] = q[date[J - 1] - 1] = q_ZX;V[date[J - 1] - 1] = ZVqChaZhi(ZSeries, VSeries, ZX);int qMaxGet = 0;//没有过最大泄流量的时刻DestroyZFH = ShouWei = 0;//下游防洪安全还没破坏for(int i = date[J - 1]; i < date[J - 1] + WProT; i++){if((SJWPro[BZ][i]<q_ZX&&qMaxGet == 0)||ShouWei == 1)//i = 8~13时满足此条件 SJWPro[0][13] = 5628.24{//无须进行调洪计算q[i] = SJWPro[BZ][i];//V[i] = ZVqChaZhi(ZSeries, VSeries, ZX);Z[i] = ZX;}else{//开始调洪计算q[i] = qFunction(BZ, i, SJWPro[BZ][i-1],SJWPro[BZ][i], V[i-1], q[i-1], 0);V[i] =  V[i - 1] + QToW*(SJWPro[BZ][i - 1] + SJWPro[BZ][i] - q[i - 1] - q[i])/2;Z[i] = ZVqChaZhi(VSeries, ZSeries, V[i]);if(Z[i]>ZFH) DestroyZFH =1;if(DestroyZFH ==0&&q[i] > SafeQ) {//此时人为控制泄流，水位泄流关系破坏q[i] = SafeQ;V[i] =  V[i - 1] + QToW*(SJWPro[BZ][i - 1] + SJWPro[BZ][i] - q[i - 1] - q[i])/2;Z[i] = ZVqChaZhi(VSeries, ZSeries, V[i]);}if(qMaxGet == 1&&q[i]<q_ZX) {//q[i] = SJWPro[BZ][i];//此次洪水已过,且水位放至限制洪水位了//V[i] = V[i - 1];//Z[i] = Z[i - 1];Z[i] = ZX;V[i] = ZVqChaZhi(ZSeries, VSeries, ZX);q[i] = SJWPro[BZ][i - 1] + SJWPro[BZ][i]  - q[i - 1]- (V[i] - V[i - 1])*2.0/QToW;ShouWei = 1;//V[i] =  V[i - 1] + QToW*(SJWPro[BZ][i - 1] + SJWPro[BZ][i] - q[i - 1] - q[i])/2;}if(q[i]>=SJWPro[BZ][i]&&qMaxGet == 0)//要求泄流最大值了{qMax(BZ, i, SJWPro[BZ][i -1], SJWPro[BZ][i], V[i-1], q[i -1]);qMaxGet = 1;//已经过了最大泄流量的时刻}}outfile<<setw(10)<<3*i<<setw(10)<<SJWPro[BZ][i]   <<setw(10)<<q[i]<<setw(10)<<V[i]<<setw(10)<<Z[i]<<endl;}outfile<<"最大泄流时各值如下："<<endl<<setw(10)<<t_qMax<<setw(10)<<Q_qMax<<setw(10)<< q_qMax<<setw(10)<< V_qMax<<setw(10)<< Z_qMax;}outfile.close();cout<<"调洪结束，请关闭！"<<endl;cin>>temp_W;//控制台暂停}

//水能计算.h//单位：水量10 ^6m3，流量m3/s，出力 万kw//事先确定好供水蓄水关系const double N_XZ = 68.4;//装机容量(若求无限装机情况,此值改成极大值)100000;//const int HLowSeries = 14,//下游水位流量插值系列数据HUpSeries = 60;//上游水位库容插值系列数据const double kFactor = 8.2,//出力系数CostH = 0.5,//水头损失VDead = 185.68;//死库容,    double deltaH,//净水头QLoss_XZ = N_XZComeQ[12], LossQ[12],//入库、出库流量ThrowQ[12],//弃水VAbs[12],//绝对库存水量 ChuLi[12],//出力HLow_h[HLowSeries],HLow_Q[HLowSeries],HUp_h[HUpSeries],HUp_V[HUpSeries];double HLowChaZhi(double x[HLowSeries], double y[HLowSeries], double x0){//下游水位流量插值double y0 = -1;//返回为负值时，可以作为判定依据for(int i = 0; i < HLowSeries - 1; i++){if(((x[i]>=x0)&&(x[i+1]<=x0))||((x[i]<=x0)&&(x[i+1]>=x0))) y0 = y[i] + (y[i+1] - y[i])/(x[i+1] - x[i])*(x0 -x[i]);}return y0;}double HUpChaZhi(double x[HUpSeries], double y[HUpSeries], double x0){//上游水位库容插值,库容单位10^6m3double y0 = -1;//返回为负值时，可以作为判定依据for(int i = 0; i < HUpSeries - 1; i++){if(((x[i]>=x0)&&(x[i+1]<=x0))||((x[i]<=x0)&&(x[i+1]>=x0))) y0 = y[i] + (y[i+1] - y[i])/(x[i+1] - x[i])*(x0 -x[i]);}return y0;}void ChuLiJS(){//出力计算using namespace std;ifstream infile;infile.open("infile_HLowSeries.txt");for(int i = 0; i < HLowSeries; i++)infile>>HLow_h[i]>>HLow_Q[i];infile.close();infile.open("infile_HUpSeries.txt");for(int i = 0; i < HUpSeries; i++)infile>>HUp_h[i]>>HUp_V[i];infile.close();infile.open("infile_ComeLossQ.txt");for(int i = 0; i < 12; i++)infile>>ComeQ[i]>> LossQ[i];infile.close();ofstream outfile;outfile.open("outfile_ChuLi.txt");outfile<<setw(10)<<"出力"<<setw(10)<<"弃水"<<endl;for(int i = 0; i < 12; i++){if(i == 0) {VAbs[i] = VDead + (ComeQ[i] - LossQ[i])*2.63;deltaH = HUpChaZhi(HUp_V, HUp_h, (VAbs[i]+VDead)/2.0);} else {VAbs[i] = VAbs[i - 1] + (ComeQ[i] - LossQ[i])*2.63; deltaH= HUpChaZhi(HUp_V, HUp_h, (VAbs[i-1]+VAbs[i])/2.0);}deltaH -= (HLowChaZhi(HLow_Q, HLow_h, LossQ[i]) + CostH);ChuLi[i] =deltaH*kFactor* LossQ[i]/10000.0;//单位 万kwThrowQ[i] = 0; if(ChuLi[i] > N_XZ) { ChuLi[i] = N_XZ;ThrowQ[i] = LossQ[i] - 10000.0*N_XZ/kFactor/deltaH; }outfile<<setw(10)<<ChuLi[i]<<setw(10)<<ThrowQ[i]<<endl;}}

//由工作容量求保证电能.h//本课设中供水期为冬季的//单位  负荷:万kw  工作容量: 万kw   电能: 亿kw*hconst int NumN = 20;//假定的工作容量数目const double MaxWorkLoad[3] = {260, 255, 250};//12、1、2三个月对应的最大负荷，降序排序double N_g[NumN],//假定的工作容量系列(对应供水期保证出力)E_g[NumN],//假定的工作容量对应的供水期保证电能DayWorkLoad[24],//冬季日负荷占对应月的最大负荷比率LowLine;void NandERelation(){//工作容量与供水期保证电能的关系using namespace std;ifstream infile;infile.open("infile_DayWorkLoad.txt");for(int i = 0; i < 24; i++)infile>>DayWorkLoad[i];infile.close();ofstream outfile;outfile.open("outfile_NandERelation.txt");cout<<"输入假设的工作容量下限"<<endl;cin>>N_g[0];// 假设的最小工作容量cout<<"输入假设的工作容量上限"<<endl;cin>>N_g[NumN  - 1];// 假设的最大工作容量for(int i = 0; i < NumN; i++){N_g[i] = N_g[0] + (double)(N_g[NumN  - 1] - N_g[0])/(NumN - 1)*i;E_g[i] = 0;for(int Month = 0; Month < 3; Month++){LowLine =MaxWorkLoad[0] - N_g[i];for(int Hour = 0; Hour < 24; Hour++)if((DayWorkLoad[Hour]/100.0*MaxWorkLoad[Month])>LowLine)E_g[i] += ((DayWorkLoad[Hour]/100.0*MaxWorkLoad[Month]) - LowLine);}E_g[i] *= (30.4/10000.0);outfile<<setw(10)<<N_g[i]<<setw(10)<<E_g[i]<<endl;}outfile.close();}