Projets numériques

Des projets numériques se rapportant au cours ont été conçus en 2017-2018 avec Victor DUPUIS, doctorant à cette époque. Victor a encadré ces projets et écrit les codes en Scilab.

MéthodeTD numérique 1 : Écoulement sur un obstacle

TD1.pdf [pdf]
1
%-------------------------------------------------------------------------%
2
%
3
% Cours ENSEEIHT "Hydraulique à surface libre"
4
% TD numerique n°1  Ecoulement sur un obstacle
5
% 2017-2018
6
%
7
%-------------------------------------------------------------------------%
8
9
clear all;
10
close all;
11
format long;
12
13
%---- Hauteur critique en m
14
hc = 0.05; 
15
16
%---- Maillage le long de x
17
xmin = -1;    % en m
18
xmax = 1;     % en m
19
nx = 201;     % nbr de points en x
20
x = linspace(xmin,xmax,nx);
21
dx = (xmax-xmin)/(nx-1); % Pas de maillage
22
23
%---- Bathymetrie du fond du cours d'eau 
24
25
%-- Obstacle de forme gaussienne
26
a = 0.03;     % hauteur maximale de l'obstacle en m
27
sigma = 0.2;  % largeur caractéristique de l'obstacle en m
28
Zf = @(x) 0 % A REMPLIR
29
% Calcul de la pente I = -dZf/ds
30
I = @(x) 0 % A REMPLIR
31
32
33
34
%---- Tracer bathymetrie et pente
35
figure(1),plot(x,Zf(x));
36
xlabel('x [m]') 
37
ylabel('Bathymétrie Z_f (m)')
38
39
figure(2),plot(x,I(x));
40
xlabel('x [m]') 
41
ylabel('Pente I')
42
43
44
%---- Definition de la condition initiale
45
h0 = 0.12; % en m
46
47
48
49
%---- Résolution EDO Euler
50
51
h_eu=zeros(nx,1);
52
53
% A REMPLIR
54
    
55
56
57
%---- Résolution EDO ode45
58
59
% A REMPLIR
60
    
61
    
62
%---- Résolution polynôme
63
64
h_poly=zeros(nx,3);
65
66
% A REMPLIR
67
    
68
69
%---- Résolution EDO 2D
70
71
% A REMPLIR
72
73
74
%---- Figures
75
76
%-- Hauteur d'eau
77
figure(3)
78
hold on
79
%plot(x,h_eu,'r-','Linewidth',3)
80
%plot(x,h_ode,'g-','Linewidth',3)
81
%plot(x,h_poly(:,1),'b-','Linewidth',3)
82
%plot(x,h_poly(:,2),'m-','Linewidth',3)
83
%plot(x_ode2D,h_ode2D,'c-','Linewidth',3)
84
plot(x,x*0+hc,'k--','Linewidth',3)
85
xlabel('x (m)')
86
ylabel('Hauteur d''eau h (m)')
87
legend('Euler','ode45','polynome 1/2','polynome 2/2','courbe param.')
88
89
%-- Cote d'eau
90
figure(4)
91
hold on
92
%plot(x,h_eu+Zf(x),'r-','Linewidth',3)
93
%plot(x,h_ode+Zf(x),'g-','Linewidth',3)
94
%plot(x,h_poly(:,1)+Zf(x),'b-','Linewidth',3)
95
%plot(x,h_poly(:,2)+Zf(x),'m-','Linewidth',3)
96
%plot(x_ode2D,h_ode2D++Zf(x_ode2D),'c-','Linewidth',3)
97
plot(x,Zf(x),'Color','k','Linestyle','-','Linewidth',3)
98
plot(x,Zf(x)+hc,'Color','k','Linestyle','--','Linewidth',3)
99
xlabel('x (m)')
100
ylabel('Cote d''eau z = Zf + h (m)')
101
legend('Euler','ode45','polynome 1/2','polynome 2/2','courbe param.')
102
103
MéthodeTD numérique 2 : Courbes de remous
TD2.pdf [pdf]
1
%-------------------------------------------------------------------------%
2
%
3
% Cours ENSEEIHT "Hydraulique à surface libre"
4
% TD numerique n°2  Courbes de remous
5
% 2017-2018
6
%
7
%-------------------------------------------------------------------------%
8
9
10
11
clear all;
12
close all;
13
format long;
14
15
%---- Cas fluvial (=1) ou cas torrentiel (=0)
16
fluvial=1;
17
18
%---- Parametres
19
Q = 150; % debit en m3/s
20
B = 50;  % largeur en m
21
hc = 0; % A REMPLIR  % hauteur critique en m
22
23
%---- pente I(x)
24
I = @(x) 0.001;
25
%---- rugosité Ks(x)
26
Ks = @(x) 30;
27
%---- hauteur normale hn(x)
28
hn = @(x) 0; % A REMPLIR
29
30
%---- condition initiale en m
31
h0=0; % A REMPLIR
32
33
%----- Maillage le long de l'abscisse curviligne
34
xmin = -10000;    % en m
35
xmax = 0;     % en m
36
nx = 2000;     % nbr de points en x
37
x = linspace(xmin,xmax,nx);
38
%----- Pas d'espace
39
dx = (xmax-xmin)/(nx-1);
40
41
42
43
%---- Résolution EDO Euler
44
45
% A REMPLIR
46
47
48
%---- Résolution EDO ode45
49
50
% A REMPLIR
51
    
52
    
53
54
%---- Figures
55
56
%---- Solution (hauteur d'eau)
57
figure(1)
58
hold on
59
%plot(x,h,'r-','Linewidth',3)
60
%plot(xode,hode,'g-','Linewidth',3)
61
plot(x,x*0+hc,'--k','Linewidth',3)
62
xlabel('x [m]')
63
ylabel('Hauteur d''eau h(x) [m]')
64
legend('Euler','ode45')
65
MéthodeTD numérique 3 : Propagation d'une crue
TD3.pdf [pdf]
1
%-------------------------------------------------------------------------%
2
%
3
% Cours ENSEEIHT "Hydraulique à surface libre"
4
% TD numerique n°3 : Propagation d'une crue 1/2
5
% 2017-2018
6
% 
7
%-------------------------------------------------------------------------%
8
9
clear all;
10
close all;
11
format long;
12
13
%---- Parametres
14
Ks = 35;                          % coefficient de Strickler
15
I = 0.002;                        % pente du fond
16
h_m = 1;                          % hauteur d'eau uniforme en m
17
U_m = 0; % A REMPLIR    	      % vitesse uniforme en m/s
18
a=5/3*U_m;                        % coefficient d'advection
19
L = 200*10^3;                     % Longueur domaine en m
20
T = L/a;                          % Temps total de propagation en s
21
22
%---- Discretisation spatio-temporelle
23
nx = 400;                         % Nombre de pas d'espace
24
nt = 1000;                          % Nombre de pas de temps
25
dx = L/nx;              	      % Pas d'espace en m
26
dt = T/nt;                        % Pas de temps en s                       
27
28
CFL = a*dt/dx;                    % Condition CFL
29
30
x = linspace(0,L,nx);             % Maillage spatial 
31
32
33
%---- Condition initiale sur la hauteur d'eau
34
h_a = 0.1; 
35
x_a = 20000;
36
sig = 4000;
37
h0 = 0;  % A REMPLIR
38
39
h_ftbs=zeros(nt,nx);
40
h_ctcs=zeros(nt,nx);
41
42
h_ftbs(1,:) = h0;
43
h_ctcs(1,:) = h0;
44
45
46
%---- Condition aux limites amont sur la hauteur d'eau
47
48
% A REMPLIR
49
50
51
52
%---- Résolution FTBS
53
54
% A REMPLIR
55
56
57
%---- Résolution CTCS ou Leapfrog
58
59
% A REMPLIR
60
61
62
63
64
65
figure(1)
66
plot(x,h0,'k-','Linewidth',3)
67
hold on
68
% plot(x,h_ftbs(round(nt/3),:),'r-','Linewidth',3)
69
% plot(x,h_ftbs(round(2*nt/3),:),'g-','Linewidth',3)
70
% plot(x,h_ctcs(round(nt/3),:),'r.','Linewidth',3)
71
% plot(x,h_ctcs(round(2*nt/3),:),'g.','Linewidth',3)
72
73
74
75
76
77
%-------------------------------------------------------------------------%
78
%
79
% Cours ENSEEIHT "Hydraulique à surface libre"
80
% TD numerique n°3 : Propagation d'une crue 2/2
81
% 2017-2018
82
% 
83
%-------------------------------------------------------------------------%
84
85
clear all;
86
close all;
87
format long;
88
89
%---- Parametres
90
Ks = 35;                          % coefficient de Strickler
91
I = 0.002;                        % pente du fond
92
h_m = 1;                          % hauteur d'eau uniforme en m
93
U_m = 0; % A REMPLIR    	      % vitesse uniforme en m/s
94
a=5/3*U_m;                        % coefficient d'advection
95
L = 200*10^3;                     % Longueur domaine en m
96
T = L/a;                          % Temps total de propagation en s
97
98
%---- Discretisation spatio-temporelle
99
nx = 400;                         % Nombre de pas d'espace
100
nt = 1000;                        % Nombre de pas de temps
101
dx = L/nx;              	      % Pas d'espace en m
102
dt = T/nt;                        % Pas de temps en s                   
103
104
CFL = a*dt/dx;                    % Condition CFL
105
106
x = linspace(0,L,nx);             % Maillage spatial 
107
t = linspace(0,T,nt);             % Maillage temporel 
108
109
%--- Condition initiale sur la hauteur d'eau
110
111
h_a = 1; % en m
112
x_a = 50000;
113
sig = 20000;
114
h0 = 0; % A REMPLIR
115
116
117
%---- Calcul des caractéristiques
118
119
x_car=zeros(nt,nx);
120
121
% A REMPLIR
122
123
%---- Tracé des caractéristiques dans le plan (t,x)
124
125
figure(1)
126
hold on
127
128
for j=1:10:nx
129
plot(x_car(:,j),t)    
130
end
131
132
133
%---- Calcul de la première coupure
134
135
t_c=0; % A REMPLIR
136
137
138
%---- Tracé de la solution à un temps t
139
140
figure(2)
141
hold on
142
plot(x_car(1,:),h0,'k-','Linewidth',3)
143
% plot(x_car(round(t_c/dt)/2,:),h0,'b-','Linewidth',3)
144
% plot(x_car(round(t_c/dt),:),h0,'g-','Linewidth',3)
145
% plot(x_car(round(2*t_c/dt),:),h0,'r-','Linewidth',3)
146
147
MéthodeTD numérique 4 : Ondes dans un port
TD4.pdf [pdf]
1
%-------------------------------------------------------------------------%
2
%
3
% Cours ENSEEIHT "Hydraulique à surface libre"
4
% TD numerique n°4 : Ondes dans un port 
5
% 2017-2018
6
% 
7
%-------------------------------------------------------------------------%
8
9
clear all;
10
close all;
11
format long;
12
13
%---- Parametres
14
L = 200;                     % Longueur domaine en m
15
T = 20;                          % Durée modélisée en s
16
g = 9.81;
17
18
%---- Discretisation spatio-temporelle
19
nx = 200;                         % Nombre de pas d'espace
20
nt = 4000000;                          % Nombre de pas de temps
21
dx = L/nx;              	      % Pas d'espace en m
22
dt = T/nt;                        % Pas de temps en s                       
23
24
x = linspace(0,L,nx);             % Maillage spatial 
25
26
27
%---- Condition initiale et aux limites sur la hauteur d'eau
28
29
h=zeros(nt,nx);
30
31
h_ref=15;
32
33
% gaussienne
34
h_a = 0.5; 
35
x_a = 100;
36
sig = 5;
37
h0_g = 0; % A REMPLIR
38
% sinusoïde
39
h_b = 0.5; 
40
lambda = 4/5*200;
41
h0_s =0;  % A REMPLIR
42
43
% A REMPLIR
44
45
%---- Condition initiale et aux limites sur la vitesse
46
47
u=zeros(nt,nx);
48
49
% A REMPLIR
50
51
%---- Résolution par différences finies
52
53
tic
54
55
% A REMPLIR
56
57
toc  % chronomètre
58
59
60
%--- Initialisation du graphique
61
figure; 
62
hold on;
63
graph = plot(x,h(1,:),'b','linewidth',2);
64
set(gca,'Linewidth',1.5,'Fontsize',16)
65
axis([0 L 14.4 15.6])
66
%axis([0 L -0.6 0.6])
67
set(gca,'Units','centimeters')
68
set(gcf,'Units','centimeters')
69
set(gca,'Position',[3,3,15,15])
70
set(gcf,'Position',[2,4,20,20])
71
xlabel('x (m)','Fontsize',16,'Fontname','Arial')
72
ylabel('h (m)','Fontsize',16,'Fontname','Arial')
73
set(gca,'nextplot','replacechildren');
74
grid on
75
76
for it = 1:3000:nt
77
    set(graph,'YData',h(it,:));  % Mise a jour du graphique
78
    movie_h(it) = getframe; 
79
    
80
end
81