Bonjour,

Je dois modeliser un transfer de chaleur sur un cible ayant la forme d'un disque.

J'ai a la fois de la conduction partout sur mon disque, et aussi de la radiation sur le pourtour de mon disque.

J'utilise l'equation du transfert de charler pour modeliser ma conduction :

dU/dt = k * Nabla * U

ou Nabla est le laplacien que j'ai mis en coordonne polaire. Mon probleme c'est que puisque je suppose que ma chaleur s'en va uniquement par radiation, et non par conduction, j'ai beaucoup de misere a imposer cette condition dans mon probleme.

Je resoud le probleme numeriquement, avec la methode des differences fini et des condition contour periodique en theta et sans transfert de chaleur sur le pourtout de mon disque. Malgre tout, je constate numeriquement que je perd de la chaleur a travers mon processus de conduction ... ou est l'erreur ????

Si ca peut aider quelqu'un voici le code matlab :


Code:
function GAUSS (fx,fy,phi,wx,wy,ax,ay,per)
% wx : half width of the guassian in X direction
% wy : half width of the guassian in Y direction
% fx : frequency of the oscillation in X direction  (in Hz)
% fy : frequency of the oscillation in Y direction  (in Hz)
% phi : phase difference between of oscillation in x and oscillation in Y
%       (in degree)
% ax : amplitude of the oscillation in x
% ay : amplitude of the oscillation in y
% per : period to do the average  (in second)

nelm=50; % space resolution
temp=10000; % time resolution
L=0.013*1e6; % length of the target
ka=57.5/1e6; % thermal conductivity of the target
e=0.92; % radiation coefficent
sr=5.6704e-8/1e12; % radiation constant
ray=0.019*1e6/2; % radius of the target
r=linspace(ray/nelm,ray,nelm);
teta=linspace(0,2*pi,nelm);
[TETA,R]=meshgrid(teta,r);
x=r'*cos(teta);
y=r'*sin(teta);
dr=ray/nelm;
dteta=2*pi/nelm;
dx=2*ray/nelm;
% rr=repmat(r',[1,nelm]);
T=linspace(0,per,temp);
dt=per/temp;
ave = zeros(nelm,nelm);
gaussm=ave;
dz=L/525;

norm=1.7264e10/(pi*ray^2/nelm^2);      
alpha=0.000024582*1e12;

cenx=ax*sin(2*pi*fx*T);
ceny=ay*sin(2*pi*fy*T+phi/360*2*pi);

Qcond=ave;
Qradl=ave;

mask = (R==ray);


for j=1:temp
    gaussm=0.4*exp(-4*log(2)*(1/wx)^2*(x-cenx(j)).^2).*exp(-4*log(2)*(1/wy)^2*(y-ceny(j)).^2);

    if j>1

        Qcond(2:nelm-1,1)=alpha.*((1./R(2:nelm-1,1)).*(ave(3:nelm,1)-ave(1:nelm-2,1))./(2.*dr)+(ave(3:nelm,1)-2.*ave(2:nelm-1,1)+ave(1:nelm-2,1))/(dr^2)+(1./R(2:nelm-1,1).^2).*(ave(2:nelm-1,2)-2*ave(2:nelm-1,1)+ave(2:nelm-1,nelm))./(dteta.^2));

            Qcond(2:nelm-1,2:nelm-1)=alpha.*((1./R(2:nelm-1,2:nelm-1)).*(ave(3:nelm,2:nelm-1)-ave(1:nelm-2,2:nelm-1))./(2.*dr)+(ave(3:nelm,2:nelm-1)-2.*ave(2:nelm-1,2:nelm-1)+ave(1:nelm-2,2:nelm-1))/(dr^2)+(1./R(2:nelm-1,2:nelm-1).^2).*(ave(2:nelm-1,3:nelm)-2*ave(2:nelm-1,2:nelm-1)+ave(2:nelm-1,1:nelm-2))./(dteta.^2));

        Qcond(1,2:nelm-1)=alpha.*((1./R(1,2:nelm-1)).*(ave(2,2:nelm-1)-ave(1,2:nelm-1))./(dr)+(ave(3,2:nelm-1)-2.*ave(2,2:nelm-1)+ave(1,2:nelm-1))/(2*dr^2)+(1./R(1,2:nelm-1).^2).*(ave(1,3:nelm)-2*ave(1,2:nelm-1)+ave(1,1:nelm-2))./(dteta.^2));
        Qcond(1,1)=alpha.*((1./R(1,1)).*(ave(2,1)-ave(1,1))./(dr)+(ave(3,1)-2.*ave(2,1)+ave(1,1))/(2*dr^2)+(1./R(1,1).^2).*(ave(1,2)-2*ave(1,1)+ave(1,nelm))./(dteta.^2));
        Qcond(1,nelm)=alpha.*((1./R(1,nelm)).*(ave(2,nelm)-ave(1,nelm))./(dr)+(ave(3,nelm)-2.*ave(2,nelm)+ave(1,nelm))/(2*dr^2)+(1./R(1,nelm).^2).*(ave(1,1)-2*ave(1,nelm)+ave(1,nelm-1))./(dteta.^2));

        Qcond(nelm,2:nelm-1)=alpha.*((2*ave(nelm-1,2:nelm-1)-2.*ave(nelm,2:nelm-1))/(dr^2)+(1./R(nelm,2:nelm-1).^2).*(ave(nelm,3:nelm)-2*ave(nelm,2:nelm-1)+ave(nelm,1:nelm-2))./(dteta.^2));
        Qcond(nelm,1)=alpha.*((2*ave(nelm-1,1)-2.*ave(nelm,1))/(dr^2)+(1./R(nelm,1).^2).*(ave(nelm,2)-2*ave(nelm,1)+ave(nelm,nelm))./(dteta.^2));
        Qcond(nelm,nelm)=alpha.*((2*ave(nelm-1,nelm)-2.*ave(nelm,nelm))/(dr^2)+(1./R(nelm,nelm).^2).*(ave(nelm,1)-2*ave(nelm,nelm)+ave(nelm,nelm-1))./(dteta.^2));

        Qcond(2:nelm-1,nelm)=alpha.*((1./R(2:nelm-1,nelm)).*(ave(3:nelm,nelm)-ave(1:nelm-2,nelm))./(2.*dr)+(ave(3:nelm,nelm)-2.*ave(2:nelm-1,nelm)+ave(1:nelm-2,nelm))/(dr^2)+(1./R(2:nelm-1,nelm).^2).*(ave(2:nelm-1,1)-2*ave(2:nelm-1,nelm)+ave(2:nelm-1,nelm-1))./(dteta.^2));

        Qradl = mask.*e.*sr.*R.*dteta.*dz.*((norm*ave).^4);
        Qdiss = Qcond-Qradl;
        ave = ave + dt*Qdiss;

    end
    ave = ave + dt*gaussm;
end


% ave normalisation :
ave=ave.*norm;
aveT=sum(sum(ave))/nelm^2


figure
mesh(x,y,Qdiss)
xlabel('X')
ylabel('Y')
zlabel('Z')