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
| %----------------------------------------------------------------------------
% Script de calcul par differences finies 2D, stationnaire.
% E. Lefrancois, 12/2008
%
% Resolution de l'equation de la chaleur :
%
% kd (T,xx + T,yy) + f0 = 0
%
% sur un domaine rectangulaire de longueur L (selon x) et hauteur H (selon y)
%
% Conditions aux limites :
% T(x,0)= 10 C
% T(O,y)=T(L,y)=T(x,H)=100
%
% Maillage : i=1,nx selon x et j=1,ny selon y
%
% Variables classees avec sens de lecture horizontal : T(1,1); T(2,1)...
%
%
% Forme discrete : ligne (i,j) indice : (j-1)nx+i
%
% A1*T(i,j-1)+A2*T(i-1,j)+A3*T(i,j)+A4*T(i+1,j)+A5*T(i,j+1)+f(i,j)=0
%
% avec : A1=kd/dy^2, A2=kd/dx^2, A3=-2(kd/dy^2+kd/dy^2), A4=kd/dx^2, A5=kd/dy^2
%
% "5"
% Les indices associes sont : (i,j+1)
% (i,j-1) => N1 = i+(j-2)*nx |
% (i-1,j) => N2 = i-1+(j-1)*nx |
% (i, j) => N3 = i+(j-1)*nx (i-1,j)-----(i,j)-----(i+1,j)
% (i+1,j) => N4 = i+1+(j-1)*nx "2" |"3" "4"
% (i,j+1) => N5 = i+j*nx |
% (i,j-1)
% "1"
%----------------------------------------------------------------------------
clear all % nettoyage de la memoire
close all % fermeture des fenetres graphiques
%----- parametres physiques
L = 10; % longueur du domaine m
H = 10; % hauteur du domaine m
kd=2; % coeff de conductivite w/°c-m
f0=30; % production w/m3
%----- parametres numeriques
nx=25; %input('Entrer le nombre de noeuds : ');
ny=25; %input('Entrer le nombre de noeuds : ');
dx = L / (nx - 1); dy = H / (ny - 1); % pas de discretisation
nnt=nx*ny; % nombre de noeuds total
k=zeros(nnt,nnt); % initialisation de la matrice
f=zeros(nnt,1); % initialisation du second membre
%
%----- Generation du maillage
for j=1:ny
for i=1:nx
X(i,j)=(i-1)*dx;
Y(i,j)=(j-1)*dy;
end
end
%
%----- Coefficients de la matrice
Ax=kd./dx^2; Ay=kd/dy^2; fij=-f0;
%
for j=2:ny-1
for i=2:nx-1
% Calcul de l'indice du noeud (i,j)
N1 = i+(j-2)*nx;
N2 = i-1+(j-1)*nx;
N3 = i+(j-1)*nx; % Noeud (i,j) !!
N4 = i+1+(j-1)*nx;
N5 = i+j*nx;
%
% On remplit la ligne "N3" du systeme
k(N3,[N1 N2 N3 N4 N5])=[Ay Ax -2*(Ax+Ay) Ax Ay];
f(N3)=fij;
end
end
%----- conditions aux limites de Dirichlet
% T=10 en (x,O)
liste=1:nx;
for m=1:nx
k(liste(m),liste(m))=1; f(liste)=10;
end
% T=100 en (x,H)
liste=(ny-1)*nx+[1:nx];
for m=1:nx
k(liste(m),liste(m))=1; f(liste)=100;
end
% T=100 en (0,y)
liste=[1:nx:nnt];
for m=1:ny
k(liste(m),liste(m))=1; f(liste)=100;
end
% T=100 en (L,y)
liste=[nx:nx:nnt];
for m=1:ny
k(liste(m),liste(m))=1; f(liste)=100;
end
%----- resolution
T = k\f; % calcul et affichage de la solution
%----- affichage
m=0;
for j=1:ny
for i=1:nx
m=m+1;
matT(i,j)=T(m);
end
end
surf(X,Y,matT)
xlabel('X (m)')
ylabel('Y (m)')
zlabel('Temperature (C)') |
Comprendre un programme Matlab?
Répondre avec citation
Envoyé par Caro-Line
Partager