M.Pilotti: Classwork of Environmental Hydraulics

M.Pilotti: Classwork of Environmental Hydraulics

Classwork 1

A computer code for computing basic hydraulic quantities for free surface flow in a rectangular cross section

Why classwork 1 ?

In this classwork we propose students to make acquaintance with a programming tool that will reveal itself invaluable in the calculations of some problems of great cultural and technical relevance.

To this purpose it is not always possible to use a spreadsheet and in some cases it seems more appropriate to write your own code. On the other hand, this effort has a great educative relevance because by implementing an algorithm students have the unique opportunity of understanding how a computer code is structured. This will give them a considerable insight also on the limitations and problems of commercial codes that they will used during their future activities. In this sense, the use of specific languages ​​for numerical computation (eg, Matlab) appears here not fully suitable to our purposes. On the other hand, Matlab will be here extensively used for graphics production.

We are aware that this goal is ambitious because of the unfamiliarity of most students with the basics of programming. We will attempt to get through a phased approach leading step by step to the implementation of simple routines that calculate geometric quantities and hydraulic fundamentals. It is believed that this effort will enrich the spirit of criticism of the future practitioner and will gratifies the future researcher with a tool that will prove an invaluable tool for your reserch and professional activities.

The programming language that will be used is PASCAL. The choice has been done because of its extreme semplicity and clarity, because many students already had previous experiences with it and, finally, because it is suitable also for scientific computations without significative limitations. In addition, a freeware versione of Borland Pascal 7 is available in the web. This software, altough suitable only for non-graphical application, perfectly suits our needs.

In order to import the lines of code contained in this practical work into the compiler, save this file as a text file and then read it into the compiler, following your instructor’s indications.

Final targets of Classwork 1

• Learning to write an algorithm in a programming language
• Understanding that it is possible to extend greatly our capabilities by joining the depth of theory to the potentials of numerical computation.
• Start writing a preliminary bricks of a larger code for free surface flow computations, so understanding that the latin motto divide et impera is a true also in the programming acticvity.

The content of classwork 1

• Basic structure of a program in PASCAL: in this example, for clarity’s sake, reserved words ar written in upper cases whilst variables are written in lower cases. Mind: this is not a rule of this language that in itself in unsensitive to the type of cases !

PROGRAM prova;
CONST
g = 9.8065;
namefile = ‘c:getta.via’;
VAR
I,j : INTEGER;
arch : TEXT;
value : REAL;
FUNCTION cubo(x:REAL):REAL;
BEGIN
Cubo := x*x*x;
END;
PROCEDURE aprifile_in_scrittura;
BEGIN
ASSIGN(arch,namefile);
REWRITE(arch);
END;
PROCEDURE chiudifile;
BEGIN
CLOSE(arch);
END;
BEGIN
WRITELN(‘This program makes useless things…’);
aprifile_in_scrittura;
FOR i := 1 TO 10 DO
BEGIN
WRITELN(arch,i:4,’ ‘,cubo(i/2):7:3);
WRITELN(i:4,’ ‘,cubo(i/2):7:3);
END;
WRITELN(arch,’ciclo WHILE’);
WHILE i > 0 DO
BEGIN
WRITELN(arch,i:4,’ ‘,cubo(i/2):7:3);
WRITELN(i:4,’ ‘,cubo(i/2):7:3);
i := i –1;
END;
WRITELN(arch,’ciclo REPEAT UNTIL’);
i := 5;
REPEAT
WRITELN(arch,i:4,’ ‘,cubo(i/2):7:3);
WRITELN(i:4,’ ‘,cubo(i/2):7:3);
i := i –1;
UNTIL i < 0;
chiudifile;
END.

• Preliminary statements to be placed at the beginning of the code

CONST
g = 9.8065;
ni = 1.1e-06; {viscosità cinematica acqua}
row = 998.2;
n_max_trial = 100;
FUNCTION pt(a,potenza:real):REAL;
{calcola a elevato ad una potenza}
BEGIN
IF a = 0 THEN pt := 0 ELSE
BEGIN
IF a > 0 THEN pt := exp(potenza*ln(a))
ELSE writeln('base negativa ed esponente reale');
END;
END;

• Area, wetted perimeter and hydraulic radius as a function of water depth Y: A(Y), P(Y), R(Y)

FUNCTION Area(h,Base: REAL):REAL;
{area per sezione rettangolare}
BEGIN
Area := (Base)*h;
END;
FUNCTION Perimetro(h,Base: REAL):REAL;
{perimetro idraulico per sezione rettangolare}
BEGIN
Perimetro := Base+2*h;
END;
FUNCTION R_idraulico(h,Base: REAL):REAL;
{raggio idraulico per sezione rettangolare}
BEGIN
R_idraulico := Area(h,Base)/ Perimetro (h,Base);
END;

• Specific energy as a function of Y: E(Y) and specific discharge for constant E

FUNCTION E(alfa,h,Q,base:real):REAL;
{calcolo di E(h) per sezioni rettangolare}
BEGIN
E := h +alfa*sqr(Q)/(2*g*sqr(Area(h,base)));
END;
n.b.: sqr indica nel linguaggio pascal l'operazione di elevazione di elevazione alla 2° potenza; sqrt indica invece la radice quadrata. Si tratta di funzioni predefinite.
FUNCTION q_spec(alfa,E,h:real):REAL;
{calcolo di q(h)=Q/B(h) ad E assegnata}
BEGIN
q_spec := h*sqrt( (2*g*(E-h))/alfa );
END;

• Specific Force as a function of water depth Y: (Y).

FUNCTION SpintaM(beta,h,Q,base,row:real):REAL;
{calcolo di M(h) flusso di quantità di moto per sezioni rettangolare}
BEGIN
SpintaM := +beta*row*sqr(Q)/Area(h,base);
END;
FUNCTION SpintaP(h,Q,base, row:real):REAL;
{calcolo di (h) ,risultante della forza di pressione, per sezioni rettangolare }
BEGIN
SpintaP := +1/2*row*g*h*Area(h,base);
END;
FUNCTION SpintaTot(beta,h,Q,base,row:real):REAL;
BEGIN
SpintaTot := SpintaM(beta,h,Q,base,row)+ SpintaP(h,Q,base,row);
END;

• Critical depth

FUNCTION k_critica(Q,base, alfa:real):real;
{valutazione della k critica data Q e la base della sezione rettangolare}
BEGIN
k_critica := pt(alfa*sqr(Q/base)/g,1/3)
END;

• Application of Chezy’s equation: compute Q given h, ks, j, and the cross section base

FUNCTION Q_Chezy(h,ks,j,base:real):REAL;
{calcolo della portata data Ks di Gauckler-Strickler, la pendenza j,la base e h}
VAR
A,R : real;
BEGIN
A := Area(h,Base);
R := R_idraulico(h,Base);
Q_Chezy := A*ks*pt(R,1/6)*sqrt(R*j);
END;

• Computing the head loss j

FUNCTION J_darcy( Ks,h,Q,base:real):REAL;
BEGIN
J_darcy := sqr(Q/ ks*Area(h,base)*pt(R_idraulico (h,base),2/3));
END;

• We write a simple code to compute and write onto disk E(h), (Y) and Q(Y) in uniform flow.

program Esercitazione0;
{scrittura su disco di alcune quantità relative alla sezione rettangolare}
const
dest = 'c:\atemp\';
nameout = 'E_rettangolare.out';
Q = 120.0;
base = 7.5;
alfa = 1.;
beta = 1.;
ks = 50;
pendenza = 0.001;
var
arch : text;
i : integer;
h,deltah,M,P,energia,spinta_totale,Q : real;
BEGIN
assign(arch,dest+nameout);
rewrite(arch);
writeln(arch,' [m] E [m] M[N] S_totale[N] Q[mc/s]');
writeln('Scrittura dei dati su file :',dest+nameout);
h := 0.0;
deltah := 0.05;
FOR i := 1 TO 100 DO
BEGIN
h := h +deltah;
energia := E(alfa,h,Q,base);
M := SpintaM(beta,h,Q,base,row);
P := SpintaP(h,Q,base,row);
spinta_totale := M +P;
Q := Q_Chezy(h,ks,pendenza,base);
writeln(arch,h:5:2,' ',energia:7:2,' ',M:7:2,' ',P:7:2,' ',spinta_totale :7:2, ,' ',Q :7:2);
END;
close(arch);
writeln('programma terminato: premere un tasto per uscire');
readln;
END.

M.Pilotti: Classwork of Environmental Hydraulics

M.Pilotti: Classwork of Environmental Hydraulics

(da Open Channel Hydraulics, T.Sturm. McGraw Hill)