📄 test.cpp
字号:
// test.cpp : 定义控制台应用程序的入口点。
//
#include <errno.h> //读文件用的
#include "stdafx.h"
#include<stdio.h>
#include<math.h>
#define N1 150 //断面数
#define N2 20 //断面的分层数
#define N3 50 //计算次数
#define G 9.8f
int IZ;
float A1[N1],B1[N1],C1[N1],D1[N1],E1[N1],A2[N1];
float B2[N1],DKDZJ[N1];
float C2[N1],D2[N1],E2[N1],H1[N1],H2[N1],H3[N1];
float DKDZS[N1][N2],DBDZS[N1][N2],DBDZJ[N1];
float TT;//后加上去的
int JS,JE,IK[2];//断数的起始个数
//*********************************共享去了*****************************************
float Z0[N1],Q0[N1]; //初始水深和流量
float DZ[N1],DQ[N1]; //水深和流量的增量
float BCUP[N2][2],BCDN[N2][2]; //上下游边界条件
float FJ[N1],GJ[N1];
//**********************************************************************************
//计算水面积
//*****************************************************************
// 各段面分层水深 AS对应过水面积 Z0初始水深 过水面积(要求的量) 断面起始 结束 每个断面的层数
void SufaceArea( float A[][N2], float B[][N2], float C[N1], float D[N1], int JS, int JE, int IZJ[N1])
{
//int N1=150,N2=20,N3=50;
// float A[N1][N2],
// float B[N1][N2],C[N1],D[N1],IZJ[N1];
int LevelNum,j;
float InitWaterHigh;
for(j=JS;j<=JE;j++) //求每个断面的过水面积，断面循环
{
LevelNum=IZJ[j]; //每个断面的层数
InitWaterHigh=C[j]; //断面初始水位高度
for(int i=1;i<LevelNum;i++) //求每层的过水面积
{
if(InitWaterHigh<=A[j][i]) //当前水位低于该层水位时，计算该层过水面积
break;
}
D[j]=B[j][i-1]+(B[j][i]-B[j][i-1])*(InitWaterHigh-A[j][i-1])/(A[j][i]-A[j][i-1]);//求过水面积
//计算理论参见PPT
}
}
//计算水位流量
//*****************************************************************
// 上游边界条件 要求的量 实际分层数
void WaterHigh( float A[][2], float A1, float* A2, int IA)
{
//float A[N2][2];
for(int i=1;i<IA;i++) //每一层的流量
{
if(A1<A[i][0]) //如果A1<上游边界流量条件
break;
}
*A2=A[i-1][1]+(A[i][1]-A[i-1][1])*(A1-A[i-1][0])/(A[i][0]-A[i-1][0]);//计算流量
//流量=上游入水流量+（本段流量-上段流量）*（当前-上）/（本断流量-）
}
//输入边界
//*****************************************************************
void BOUND(float DATA[][2], float A1, int IA,int J,int IGG,int KD)
{
float ZT=0.0f;
float QT,QZD,HIQ,DFQ,QIQ;
// J=J-1;
if(IGG==1)
{
WaterHigh(DATA,A1,&ZT,IA);
DZ[J]=ZT-Z0[J];
if(KD==1)
{
FJ[J]=1000000.0;
GJ[J]=-FJ[J]*DZ[J];
}
else
DQ[J]=FJ[J]*DZ[J]+GJ[J];
}
if(IGG==2)
{
WaterHigh(DATA,A1,&QT,IA);
DQ[J]=QT-Q0[J];
if(KD==1)
{
FJ[J]=0.0;
GJ[J]=DQ[J];
}
else
DZ[J]=(DQ[J]-GJ[J])/FJ[J];
}
if(IGG==3)
{
WaterHigh(DATA,Z0[J],&QZD,IA);
HIQ=Z0[J]+0.1f;
WaterHigh(DATA,HIQ,&QIQ,IA);
DFQ=(QIQ-QZD)/0.1f;
if(KD==1)
{
FJ[J]=DFQ;
GJ[J]=QZD-Q0[J];
}
else
{
DZ[J]=(QZD-GJ[J]-Q0[J])/(FJ[J]-DFQ);
DQ[J]=FJ[J]*DZ[J]+GJ[J];
}
}
}
int _tmain(int argc, _TCHAR* argv[])
{
float ZN[N3][N1],QN[N3][N1]; //计算的水深和流量
float T,temp,DZS,DBS; //
// AS面积, BS面宽， PS湿周 各断面分层水深 流量模数
float AS[N1][N2],BS[N1][N2], PS[N1][N2], ZS[N1][N2], RKS[N1][N2];
//对应ZJ[]的过水面面积，水面宽度，湿周和流量模数，用于计算流量模数，断面坐标值
float AJ[N1],BJ[N1],RKJ[N1],RNJ[N1],XL[N1];
int IZJ[N1]; //各断面对应ZS[][]的分层数
float ZETA,DT;
int NX;
FILE * fp1;
if(NULL==(fp1=fopen("E:\\WangStu\\Water Simulation\\N当前成果代码\\test\\test\\input.txt","rb+")))
{
printf("Open file 1 error:\n");
return 1;
}
printf("Open file(s) successfull!\n");
printf("Get data from file 1 : \n");
fscanf(fp1,"%f %d %f %f", &TT, &NX, &DT, &ZETA);//时间间隔，段数，时间间隔，代表权值,应该>0且<1 一般取0.5
printf("\n时间长度(HOURS)\t河道段数\t时间间隔(HOURS)\t代表权值\n%f\t%f\t%f\t%f\t\n",TT,NX,DT,ZETA);
float NT,RS; //后加上去的,未声明的标识符
int IBU,IBD; //后加上去的,未声明的标识符
NT=int(TT/DT)+1.001f;
JS=0;
JE=NX;
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f", &XL[0],&XL[1],&XL[2],&XL[3],&XL[4],&XL[5],&XL[6],&XL[7],&XL[8],&XL[9],&XL[10]); //READ(1,*) (XL(J),J=JS,JE);//各断面坐标值
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f", &RNJ[0],&RNJ[1],&RNJ[2],&RNJ[3],&RNJ[4],&RNJ[5],&RNJ[6],&RNJ[7],&RNJ[8],&RNJ[9],&RNJ[10]); //READ(1,*) (RNJ(J),J=JS,JE);
fscanf(fp1,"%d %d %d %d %d %d %d %d %d %d %d", &IZJ[0],&IZJ[1],&IZJ[2],&IZJ[3],&IZJ[4],&IZJ[5],&IZJ[6],&IZJ[7],&IZJ[8],&IZJ[9],&IZJ[10]); //READ(1,*) (IZJ(J),J=JS,JE);各断面对应zs分层数
int i;
printf("\n各断面分层水深值，水面宽度 一层 二层\n");
for(i=JS;i<JE;i++)//段循环
{
IZ=IZJ[i];//各段的分层数
fscanf(fp1,"%f %f %f %f",&ZS[i][0],&BS[i][0],&ZS[i][1],&BS[i][1]);//各断面分层水深值，水面宽度
printf("%d:\t%f\t%f\t%f\t%f\t\n",i+1,ZS[i][0],BS[i][0],ZS[i][1],BS[i][1]);
}
fscanf(fp1,"%d %d", &IK[0], &IK[1]);//READ(1,*) (IK(I),I=1,2);//边界条件指示数，IK[1]上游边界条件IK[2]下游边界条件，IK[*]=1给定水位时间过程Z(t) =2给定流量时间过程Q(t) =3给定流量水位关系Q(Z)
fscanf(fp1,"%d %d", &IBU, &IBD); //READ(1,*) IBU,IBD;上下游分层数
printf("\n边界条件类型值\n[%d]\t\t[%d]\t\t\n上下游条件的个数\n[%d]\t\t[%d]\t\t\n",IK[0],IK[1],IBU,IBD);
printf("\n读取 上 游各层的边界条件（水位时间过程） 一层\t二层\n");
for(i=0;i<IBU;i++)//读取 上 游各层的边界条件//以列表形式给出时，数据的个数
{
fscanf(fp1,"%f %f", &BCUP[i][0],&BCUP[i][1] );//READ(1,*) ((BCUP(J,I),I=1,2),J=1,IBU); //上游边界条件
printf("[%f]\t[%f]\t\n",BCUP[i][0],BCUP[i][1]);
}
printf("\n读取 下 游各层的边界条件（水位时间过程） 一层\t二层\n");
for(i=0;i<IBD;i++)//读取 下 游各层的边界条件
{
fscanf(fp1,"%f %f", &BCDN[i][0],&BCDN[i][1] );//下游边界条件
printf("[%f]\t[%f]\t\n",BCDN[i][0],BCDN[i][1]);
}
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f",&Q0[0], &Q0[1], &Q0[2], &Q0[3], &Q0[4], &Q0[5], &Q0[6], &Q0[7], &Q0[8], &Q0[9], &Q0[10]);//READ(1,*) (Q0(J),J=JS,JE); //各断面初识流量
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f",&Z0[0], &Z0[1], &Z0[2], &Z0[3], &Z0[4], &Z0[5], &Z0[6], &Z0[7], &Z0[8], &Z0[9], &Z0[10]);//READ(1,*) (Z0(J),J=JS,JE); //各断面初识水深
printf("\n各断面初始流量\n%f\t%f\t%f\t%f\t%f\t%f\t\n",Q0[0], Q0[1], Q0[2], Q0[3], Q0[4], Q0[5]);
printf("\n各断面初始水深\n%f\t%f\t%f\t%f\t%f\t%f\t\n",Z0[0], Z0[1], Z0[2], Z0[3], Z0[4], Z0[5]);
printf("--------------------文件读取结束-----------------\n");
fclose(fp1);
int j;
//未知变量赋初值
for(j=JS;j<JE;j++)//DO 30 J=JS,JE
{
AS[j][0]=0.0; //过水面积
PS[j][0]=0.0; //湿周
RKS[j][0]=0.0; //流量模数
DQ[j]=0.0; //水深增量
DZ[j]=0.0; //流量增量
}
//计算过水面面积，水面宽，流量模数 dbdz dkdz 与z的关系
for(j=JS;j<JE;j++)
{
IZ=IZJ[j];
for(int i=1;i<IZ;i++)//DO 20 i=2,IZ
{
DZS=ZS[j][i]-ZS[j][i-1]; //各断分层水深
DBS=(BS[j][i]-BS[j][i-1])/2.0f; //水面宽
AS[j][i]=AS[j][i-1]+(BS[j][i-1]+BS[j][i])*DZS/2.0f; //计算面积=下面的面积+（上底+下底）*高/2
PS[j][i]=PS[j][i-1]+2.0f*sqrtf(DBS*DBS+DZS*DZS); //计算湿周=可以理解为水与土壤的接触长度
RS=AS[j][i]/PS[j][i]; //面积/湿周
RKS[j][i]=AS[j][i]*powf(float(RS),float(2.0/3.0))/RNJ[j];//流量模数
DKDZS[j][i]=(RKS[j][i]-RKS[j][i-1])/DZS;
DBDZS[j][i]=(BS[j][i]-BS[j][i-1])/DZS;
}
DKDZS[j][0]=DKDZS[j][1];
DBDZS[j][0]=DBDZS[j][1];
}
//计算系数 A,B,C,D,E
float DX,CR1,BJM,AJM,DQJ,DZJ; //后加上去的,未声明的标识符
int JEM1,JSP1;
JEM1=JE-1;//断面数-1
JSP1=JS+1;//断面数起始+1
NT=(TT/DT)+1;//TT为读进来的值．DT为读进来的值 时间间隔，TT可以假设是时间总长度
for(int N=0;N<NT;N++)//模拟次数循环
{
T=float(N+1)*DT;
//Z0初始水深，AJ=对应ZJ过水断面积，JS=0，JE=最大断面数，IZJ=各断面对应ZS的分层数
SufaceArea(ZS,AS,Z0,AJ,JS,JE,IZJ); //ZS=各段面分层水深值，AS=对应ZS过水断面积,求对应ZJ过水断面积AJ
SufaceArea(ZS,BS,Z0,BJ,JS,JE,IZJ); //BS=对应ZS过水面宽度，BJ=对应ZJ水面宽度,求对应ZJ过水面宽度BJ
SufaceArea(ZS,RKS,Z0,RKJ,JS,JE,IZJ); //RKS=对应ZS流量模数，RKJ=对应ZJ流量模数,求对应ZJ流量模数RKJ
SufaceArea(ZS,DKDZS,Z0,DKDZJ,JS,JE,IZJ);//DKDZS=各断对应不同水深的dk/dz值，DKDZJ=（dk/dz）f值,求对应ZJ（dk/dz）f值DKDZJ
SufaceArea(ZS,DBDZS,Z0,DBDZJ,JS,JE,IZJ);//DBDZS=各断对应不同水深的db/dz值，DKDZJ=（db/dz）f值,求对应ZJ（db/dz）f值DBDZJ
for(int o=0;o<5;o++)
// printf("##### AJ=%.2f,BJ=%.2f,RKJ=%.2f,DKDZJ=%.2f,DBDZJ=%.2f\n",AJ[o],BJ[o],RKJ[o],DKDZJ[o],DBDZJ[o]);
for(j=JS;j<JEM1;j++)//断数-1个循环
{
DX = XL[j+1]-XL[j];//各断面坐标值，这里理解为段河道的长度
CR1 = ZETA*DT/DX;//ZETA读进来的
BJM = BJ[j]+BJ[j+1];//水面宽度
AJM = AJ[j]+AJ[j+1];//水断面面积
DQJ = Q0[j+1]-Q0[j];
DZJ = Z0[j+1]-Z0[j];
//计算各个参数或系数
A1[j] = -4.0f*CR1/BJM;
B1[j] = 1.0f-4.0f*CR1*ZETA*DQJ*DBDZJ[j]/BJM/BJM;
C1[j] = 4.0f/BJM*CR1;
D1[j] = 1.0f-4.0f*CR1*DQJ*DBDZJ[j+1]/BJM/BJM;
E1[j] = -4*CR1*DQJ/BJM;
A2[j] = 1.0f-4.0f*CR1*Q0[j]/AJ[j]+2.0f*CR1*DX*G*AJ[j]*fabsf(Q0[j])/RKJ[j]/RKJ[j];
B2[j] = CR1*(2.0f*Q0[j]*Q0[j]*BJ[j]/AJ[j]/AJ[j]-G*AJM+G*DZJ*BJ[j]);
temp = G*CR1*DX*Q0[j]*fabsf(Q0[j])/RKJ[j]/RKJ[j]*(BJ[j]-2*AJ[j]*DKDZJ[j]/RKJ[j]);
B2[j] = B2[j]+temp;
C2[j] = 1.0f+4.0f*CR1*Q0[j+1]/AJ[j+1]+2.0f*G*CR1*DX*AJ[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1];
temp = -2.0f*Q0[j+1]*Q0[j+1]*BJ[j+1]/AJ[j+1]/AJ[j+1]+G*AJM+G*DZJ*BJ[j+1];
D2[j] = CR1*temp;
temp = BJ[j+1]-2.0f*AJ[j+1]*DKDZJ[j+1]/RKJ[j+1];
D2[j] = D2[j]+G*CR1*DX*Q0[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1]*temp;
E2[j] = DT/DX*(-2.0f*Q0[j+1]*Q0[j+1]/AJ[j+1]+2.0f*Q0[j]*Q0[j]/AJ[j]-G*AJM*DZJ);
E2[j]=E2[j]-G*DT*(AJ[j+1]*Q0[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1]+Q0[j]*fabsf(Q0[j])*AJ[j]/RKJ[j]/RKJ[j]);
}
float AFB1,AFB2,AHC2;
int JR,IGG;
//输入边界条件 上游
IGG=IK[0];
BOUND(BCUP,T,IBU,JS,IGG,1);
//THE FIRST SWEEP
for(j=JS;j<JEM1;j++)
{
AFB1=A1[j]*FJ[j]+B1[j];
H1[j]=-C1[j]/AFB1;
H2[j]=-D1[j]/AFB1;
H3[j]=(E1[j]-A1[j]*GJ[j])/AFB1;
AFB2=A2[j]*FJ[j]+B2[j];
AHC2=AFB2*H1[j]+C2[j];
FJ[j+1]=-(AFB2*H2[j]+D2[j])/AHC2;
GJ[j+1]=(E2[j]-AFB2*H3[j]-A2[j]*GJ[j])/AHC2;
}
//输入边界条件 下游
IGG=IK[1];
BOUND(BCDN,T,IBD,JE-1,IGG,2);
//THE SECOND SWEEP
for(j=JS;j<JEM1;j++)
{
JR=JE-j-2;
DZ[JR]=H1[JR]*DQ[JR+1]+H2[JR]*DZ[JR+1]+H3[JR]; //水深的增量
DQ[JR]=FJ[JR]*DZ[JR]+GJ[JR]; //流量的增量
}
//JE个断数循环，重置参数
for(j=JS;j<JE;j++) //RESET INITIAL CONDITION
{
Q0[j]=Q0[j]+DQ[j]; //初识流量+流量增量
Z0[j]=Z0[j]+DZ[j]; //初始水深+水深增量
QN[N][j]=Q0[j]; //计算的流量
ZN[N][j]=Z0[j]; //计算的水深
}
}
//输出 计算次数，时间点，水深，流量
for(N=0;N<NT;N++)
{
float TIME=float(N+1)*DT;
printf("\n%d\t%f\t\n\t\t",N+1,TIME);
for(i=0;i<JE;i++)
{
printf("断面%d ",i);
}
printf("\n计算的流量\t");
for(i=0;i<JE;i++)
{
printf("%.3f ",QN[N][i]);
}
printf("\n计算的水深\t");
for(i=0;i<JE;i++)
{
printf("%.3f ",ZN[N][i]);
}
printf("\n");
}
return 0;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -