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
| %% initialisation
clear all;
clc;
%% paramètres numériques
L1=10; % longuer du modèl [cm]
L2=10; % largeur du modèle [cm]
TOL = 1e-6; % erreur max
T0=37; % Tempérture ambiante [°C]
nt=5000; % Nombre d'itérations
A=0.9
%% propriétés physiques
K=0.25; % conductivité thermique [Wm-1 C-1]
rho = 4000; % density [kg m-3]
Cp = 450; % capacité thermique [Jkg-1C-1]
D=(K/rho*Cp); % dufisité thermique [m2/s]
C=(1/rho*Cp); % coefficient
%% Mesh Generation
dx=0.1; % Pas du maillage selon la direction X
dy=dx; % Pas du maillage selon la direction Y
[x,y]=meshgrid(0:dx:L1,0:dy:L2); % Maillage
m = L1/dx + 1; % number of rows in mesh
n = L2/dy + 1; % number of columns in mesh
%% boundary conditions
T=T0*ones(n,m); % initialisation a t=0
T(1,:)=T0; % Temperature a y=0 (bottom side)
T(:,m)=T0; % Temperature a x=m (right side)
T(n,:)=T0; % Temperature a y=n (Top side)
%% cacactéristique du laser
%% cacactéristique du laser
p=100; % power [W]
R=0.00005; % le rayon du faisceau laser [m]
r=0.00005; % La distance radiale [m]
Qa=((2*A*p)/(pi*(R^2)))...
*(exp((-2*(r^2))/(R^2))) % Intensité pic du laser [W/m2]
Q(1:n,1:m)=0;
Q(round(n/2),round(m/2))=Qa;
%% Resolution
dt = dx^2/4;
sx=D*dt/dx^2;
sy=D*dt/dy^2;
gamma=C*dt;
tic
for k=1:nt
Told = T;
for j = 2:m-1
for i = 2:n-1
T(i,j) = sx*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j)) ...
+sy*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1)) ...
+ Told(i,j)+(gamma*Q(i,j));
end
end
if (mod(k,10)==0);% plot tous les 5 pas du temps
i=1:n;j=1:m;
imagesc(i,j,T);
colormap(jet);
colorbar;
xlabel('X');
ylabel('Y');
drawnow;
end
end
toc |
Partager