1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define NL 100
int main(int argc, char *argv[])
{
char *cheminProg = argv[0];
float u[2][NL]; //la vitesse
float h[2][NL]; //la hauteur d'eau
float Q[2][NL]; //le debit
float JC[2][NL]; //...
float tempsloc; //l'instant dans lequel on est
float dx,dt,g;
float B,m,hn,LT,Jf,n;//caractéristiques du canal
float S,P,RH,U0,Q0;//caractéristiques d'écoulement uniforme
float TP,TPD,TTH,Qmax,QE; //caractéristique de l'hydrogramme de crue
float BETA,GAMA;
float tempsfinal; //l'instant jusqu'auquel on dois calculer
FILE *vitesseFile, *hauteurFile, *debitFile;
int i, z; // iterateurs
//TODO: données consernants le canal
printf(" LARGEUR DU FOND DU CANAL ? B (m) = ");
scanf("%lf",&B);
printf(" PENTE DES BERGES ? m (-) = ");
scanf("%lf",&m);
printf(" PENTE DU FOND DU CANAL ? Jf (-) = ");
scanf("%lf",&Jf);
printf(" COEFFICIENT DE MANNING ? n (m^-1/3 s) = ");
scanf("%lf",&n);
printf(" PROFONDEUR POUR ECOULEMENT UNIFORME ? hn (m) = ");
scanf("%lf",&hn);
printf(" LONGUEUR TOTALE DU CANAL ? LT (m) = ");
scanf("%lf",<);
// calcul du debit pour écoulement uniforme
S = hn*B;
P = B+(2*hn) ;
RH = S/P;
U0 = pow(RH,(2.0/3.0))*sqrt(Jf)/n;
Q0 = U0*S;
//TODO: initialisation
for(i=0;i<NL;i++){
h[0][i]=hn;
u[0][i]=U0;
Q[0][i]=Q0;
}
// ********
//ici on calcule les pas de temps d'espace
dx=LT/NL;
dt=0.1;
// parametres consernant le bord d'entree amont
printf(" TEMPS DE MONTEE (s)? = ");
scanf("%lf",&TP);
printf(" DEBIT DE POINTE (m3/s)? = ");
scanf("%lf",&Qmax);
printf(" TEMPS DE DESCENTE (s)? = ");
scanf("%lf",&TPD);
TTH = TP + TPD;
printf("quelle est la duree des calcul (s)? =");
scanf("%lf",&tempsfinal);
tempsloc=0.;
while(tempsloc<tempsfinal){
tempsloc=tempsloc+dt;
i=0;
// condition de bord en entrée
if(tempsloc<=TP){
QE = Q0 + tempsloc * (Qmax - Q0) / TP;
}else if(tempsloc<TP+TPD){
QE = Qmax - (tempsloc - TP) * (Qmax - Q0) / TPD;
}else{
QE = Q0;
}
h[1][i]=h[0][i]-(dt/dx)*(u[0][i]*(h[0][i]-h[0][i+1])+h[0][i]*(u[0][i]-u[0][i+1]));
S=h[1][i]*B;
Q[1][i]=QE;
u[1][i]=Q[1][i]/S;
// ....
// le milieu
for(i=1;i<NL;i++){
h[1][i]=h[0][i]-dt/(2.0*dx)*(u[0][i]*(h[0][i+1]-h[0][i-1])+h[0][i]*(u[0][i+1]-u[0][i-1]));
BETA=u[0][i]-dt/(2.0*dx)*(u[0][i]*(u[0][i+1]-u[0][i-1])+g*(h[0][i+1]-h[0][i-1]))+dt*g*Jf;
RH = (B*h[1][i]) / (B + (2 * h[1][i]));
GAMA = pow(RH,(4.0/3.0))/(n*n*g*dt);
u[1][i] = 0.5*(sqrt((GAMA*GAMA)+(4.0*GAMA*BETA))-GAMA);
S=B*h[1][i];
Q[1][i]=u[1][i]*S;
printf("Q[%d]= %f\n",i,Q[1][i]);
system("pause");
}
i=NL;
// condition de bord en sortie
h[1][i]=h[0][i]-(dt/dx)*(u[0][i]*(h[0][i]-h[0][i-1])+h[0][i]*(u[0][i]-u[0][i-1]));
S=B*h[1][i];
RH=(B*h[1][i])/(B+(2*h[1][i]));
u[1][i] = pow(RH,(2./3.)) * sqrt(Jf) / n;
Q[1][i]=S*u[1][i];
// ....
//TODO procédure de stockage du pas n+1
vitesseFile = fopen(strcat(cheminProg,"\\vitesse.txt"),"w");
hauteurFile = fopen(strcat(cheminProg,"\\hauteur.txt"),"w");
debitFile = fopen(strcat(cheminProg,"\\debit.txt"),"w");
fprintf(vitesseFile, "#----Commencement d'un bloc----#\n\n");
fprintf(hauteurFile, "#----Commencement d'un bloc----#\n\n");
fprintf(debitFile, "#----Commencement d'un bloc----#\n\n");
for(i = 0; i<=2; i++)
{
fprintf(vitesseFile, "%d:", i);
fprintf(hauteurFile, "%d:", i);
fprintf(debitFile, "%d:", i);
for(z = 0; z<=NL; z++)
{
fprintf(vitesseFile, " u[%d][%d]:%8f", i, z, u[i][z]);
fprintf(hauteurFile, " h[%d][%d]:%8f", i, z, u[i][z]);
fprintf(debitFile, " Q[%d][%d]:%8f", i, z, u[i][z]);
}
fprintf(vitesseFile, "Fin de %d", i);
}
fprintf(vitesseFile, "\n#---Fin du bloc----#\n\n\n");
fprintf(hauteurFile, "\n#----Fin du bloc----#\n\n\n");
fprintf(debitFile, "\n#----Fin du bloc----#\n\n\n");
//TODO permutaion
for (i = 0;i <= NL;i++){
h[0][i] = h[1][i];
u[0][i] = u[1][i];
Q[0][i] = Q[1][i];
}
// ******
}
//fin !!!
return 0;
} |
Partager