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 |
Partager