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
| !Exercice sur le coup de bélier
!Ce programme modélise un réservoir par lequel s'échappe un débit ini Q0,
!dans un tuyaux régulé par une vanne à son extrémité.
!Pour ce premier modèle nous ne prenons pas en comtpe les pertes de charges, ainsi que les fluctuations.
!du réservoir.
program Vanne_Reservoir
!Déclaration des variables
implicit none
real, dimension(11) :: Hp, Qp, H, Q
real :: a, Din, Pi, B, K, g, L, E, ep, epsilon, rho, S, D
real :: H0, Q0, Cp, Cm, Hr
real :: t, Tmax, dt, th, Tc, mu, Tau, dtv
logical :: fini
integer :: Ns, N, i
! Données initiales
D = 100.*1e-3 !diametre de la section de la conduite
ep = 5.*1e-3 !épaisseur de la conduite
E = 2.2*1e11 !module de young du matériau de la conduite
epsilon = 2.*1e9 !déformation de la conduite
rho = 1.e3 !masse volumique de l'eau
g = 9.81 !gravité terrestre
Pi = 3.14159
Q0 = 0.015 !en m^3/s
H0 = 50 !en m.CE
Tc = 2.1 !tau = (1 - t/tc)^Em
N = 5
! Calcul des grandeurs initiales
a = 1/sqrt(rho*(1/epsilon + Din/(E*ep))) !Célérité de l'onde
S = Pi*D**2/4 !Section
NS = N + 1
B = a/(g*S)
dt = L/(a*float(N))
t = 0
! Trouver les conditions d'écoulement stationnaire
! et sauvegarder les variables initiales
do i = 1, NS
H(i) = H0
Q(i) = Q0
end do
write(*, FMT='a, L, Din=',2F8.1,2F8.4/' Hr, H0, Q0=',2F8.2,F8.3/2'Tc, &
& E=',3F8.3/' G, Tmax,dt,B=',F8.3,F8.1,2F8.3/' N =', &
& 32I4//' Heads and discharges along the pipe'//' Time X/L= 40. .2 .4 .6 .8 1. TAU') &
a, L, Din, Hr, H(NS), Q0, Tc, E, g, Tmax, dt, B, N
write(6, FMT=1H0,F7.3,5H 'H=',6F8.2,F7.3/10X,3H 'Q=',6F8.3) t, (H(i),i=1,NS),Tau, (Q(i), i=1, NS)
t = t + dt
fini = .false.
do while (t.GT.Tmax)
! Calcul des noeuds intérieur
do i = 2, N
Cp = H(i-1) + B*Q(i-1)
Cm = H(i+1) - B*Q(i+1)
Hp(i) = .5*(Cp + Cm)
Qp(i) = (Cp - Hp(i))/B
end do
!Condition limite
Hp(1) = H0
Cm = H(2) - B*Q(2)
Qp(1) = (Hp(1)-Cm)/B
Cp = H(N) + B*Q(N) !valve at end of pipe
if(t.LT.Tc) then
dtv = Tc/10
i = ifix(t/dtv)+1
Th = (t - (i-1)*dtv)/dtv
Hp(NS) = Cp - B*Qp(NS)
else
Qp(NS) = 0
Hp(NS) = Cp
end if
Do i = 1, NS
Q(i) = Qp(i)
H(i) = Hp(i)
end do
end do
end program Vanne_Reservoir |
