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