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
148
149
150
151
152
153
154
155
156
157
158
159
| clear all
close all
clc
if callback 'pushbutton1';
[filename, pathname] = uigetfile('*.dat');
carteperm = [pathname, filename];
entree= fopen (carteperm, 'r');
a = fgetl (entree);
sprintf('a');
aHndl =findobj(gcbf,'Tag','titre');
set(aHndl,'String',a);
toto = fscanf (entree, '%c', 6);
nlign = fscanf (entree, '%d',3);
toto = fscanf (entree, '%c', 6);
ncol = fscanf (entree, '%d', 3);
toto = fscanf (entree, '%c', 4);
dx = fscanf (entree, '%d', 1);
toto = fscanf (entree, '%c', 5);
dy = fscanf (entree, '%d', 1);
toto = fscanf (entree, '%c', 5);
dt = fscanf (entree, '%d', 1);
toto = fscanf (entree, '%c', 10);
hwater = fscanf (entree, '%d', 3);
toto = fscanf (entree, '%c', 9);
niveausol= fscanf (entree, '%d', 3);
Kin = fscanf (entree,'%f %f', [nlign,ncol]);
fclose (entree);
K= 10.^Kin;
h = ones(nlign, ncol)*hwater;
Lx = ncol*(dx);
Ly = nlign*(dy);
j=h(1, ncol);
x = 0 : dx : ncol-1; %position de chaques colonnes
y = 0 : dy : nlign-1;
phi=1; % porosité
%conditions initiales
h(1,:) = 0;
end
for t=dt:dt:8000 %le tps va de 10s à 20000s par pas de 10s
%conditions aux limites
h(:,1) = h(:,2);
h(:,ncol)= h(:,ncol-1);
%calculs surfaces
Sx=(h(:, 1:ncol-1)+ h(:, 2:ncol))/2*dy; %pour toutes les lignes, on calcule la moyenne des surfaces de la premiere à l'avant dernière colonne et de la 2eme à la dernière
Sy=(h(1:nlign-1, :)+ h(2:nlign,:))/2*dx; %pour toutes les colonnes, on calcule la moyenne des surfaces de la premiere à l'avant dernière ligne et de la 2eme à la dernière
%calculs débits
Qx= (-K(:, 1:ncol-1)) .* ( h(:, 2:ncol)-h(:, 1:ncol-1) )./dx .* Sx; %-K(:, 1:ncol-1)
Qy= (-K(1:nlign-1,:)) .* ( h(2:nlign ,:)-h(1:nlign-1,:) )./dy .*Sy; %(2:nlign ,:)-K(
dhx=(Qx(:,1:ncol-2)-Qx(:,2:ncol-1))/dx/dy*dt;
dhy=(Qy(1:nlign-2,:)-Qy(2:nlign-1,:))/dx/dy*dt;
h(2:nlign-1,2:ncol-1)= h(2:nlign-1,2:ncol-1)+dhx(2:nlign-1,:)+ dhy(:,2:ncol-1);
subplot (2,1,1); surf (x,y,h, 'EdgeColor', 'none'); % vue de surface
shading interp
axis([0 Lx 0 Ly 0 j])
subplot (2,1,2); surf (x,y,h), view(-90,0);% vue en coupe
axis([0 Lx 0 Ly 0 j])
drawnow
end
q=input('entrez coordonnées du forage en x : ');
p = input('entrez coordonnées du forage en y (en m) : ');
z = input('donnez une profondeur de forage (en m) : ');
%conversion des coordonnées du forage m -> pixels
p = (p./dy); % !!!!!!!! Attention à améliorer pour éviter d'avoir un résultat non entier
q = (q./dx);
r = j-z; %hauteur du fond du forage
phi = 0.6; % porosité
Qm = input('définir un débit : ');
for t = dt:dt:15000 %le tps va de 10s à 20000s par pas de 10s
h(:,1) = h(:,2);
h(:,ncol)= h(:,ncol-1);
%variations de hauteur d'eau ds forage
hp = Qm.*dt./(phi*dx*dy); % hauteur d'eau potentiellement pompée avec un débit de pompage Qm à chaque itération
he = h(p, q) - r; %hauteur d'eau disponible ds forage
deltah = min (hp, he) ;
h(p, q) = h(p, q) - deltah;
Qforage(t) = deltah.*phi.*dx.*dy./dt;
%calculs surfaces
Sx=(h(:, 1:ncol-1)+ h(:, 2:ncol))/2*dy; %pour toutes les lignes, on calcule la moyenne des surfaces de la premiere à l'avant dernière colonne et de la 2eme à la dernière
Sy=(h(1:nlign-1, :)+ h(2:nlign,:))/2*dx; %pour toutes les colonnes, on calcule la moyenne des surfaces de la premiere à l'avant dernière ligne et de la 2eme à la dernière
%calculs débits
Qx= (-K(:, 1:ncol-1)) .* ( h(:, 2:ncol)-h(:, 1:ncol-1) )./dx .* Sx; %-K(:, 1:ncol-1)
Qy= (-K(1:nlign-1,:)) .* ( h(2:nlign ,:)-h(1:nlign-1,:) )./dy .*Sy; %(2:nlign ,:)-K(
dhx=(Qx(:,1:ncol-2)-Qx(:,2:ncol-1))/dx/dy*dt;
dhy=(Qy(1:nlign-2,:)-Qy(2:nlign-1,:))/dx/dy*dt;
h(2:nlign-1,2:ncol-1)= h(2:nlign-1,2:ncol-1)+dhx(2:nlign-1,:)+ dhy(:,2:ncol-1);
cla %efface ce qui avait sur la figure
hold on
subplot (3,1,1); surf (x,y,h); % vue de surface
shading interp
axis([0 Lx 0 Ly 0 j])
subplot (3,1,2); surf (x,y,h, 'EdgeColor', 'none'), view(-90,0);% vue en coupe
shading interp
axis([0 Lx 0 Ly 0 j])
drawnow
end
subplot(3,1,3);
plot(t,Qforage,'k');
axis([0 t 0 Qm])
drawnow |
