Programation

Fermé
alkairaouni Messages postés 1 Date d'inscription mercredi 21 novembre 2012 Statut Membre Dernière intervention 21 novembre 2012 - 21 nov. 2012 à 10:44
salut,
j'ai un programme matlab de lax wendroff d'ordre 4 par la méthode de différence finie, je n'arrive pas à exécuter ce programme correctement.
et je n'arrive pas de tracer les figures de la vitesse radial sur la même figure.
**********************************
sujet:coalescence des deux bulles

**********************************
clear;
clc;
%paramétres d'entrée
N=121;
dr=0.05;
dt=0.00005;
h00=2;

%profil initial

for i=1:N
r(i)=(i-1)*dr;
h0(i)=h00+r(i)^2;
u0(i)=r(i)/(2*h0(i));
end
UHB=u0(N)*h0(N);
T_sim=0;

while(T_sim <7500*dt)

%=======================================================
SU(1)=-30*u0(1)/(12*dr^2);
SU(2)=(u0(2)+16*u0(1)-30*u0(2)+16*u0(3)-u0(4))/(12*dr^2);
FU(1)=(-u0(3)+8*u0(2)+8*u0(2)-u0(3))/(12*dr);
FU(2)=(-u0(2)-8*u0(1)+8*u0(3)-u0(4))/(12*dr);
for i=3:N-2
SU(i)=(-u0(i-2)+16*u0(i-1)-30*u0(i)+16*u0(i+1)-u0(i+2))/(12*dr^2);
FU(i)=(u0(i-2)-8*u0(i-1)+8*u0(i+1)-u0(i+2))/(12*dr);
end
FU(N-1)=(-u0(N-4)+6*u0(N-3)-18*u0(N-2)+10*u0(N-1)+3*u0(N))/(12*dr);
FU(N)=(3*u0(N-4)-16*u0(N-3)+36*u0(N-2)-48*u0(N-1)+25*u0(N))/(12*dr);
SU(N-1)=(-u0(N-4)+4*u0(N-3)+6*u0(N-2)-20*u0(N-1)+11*u0(N))/(12*dr^2);
SU(N)=(11*u0(N-4)-56*u0(N-3)+114*u0(N-2)-104*u0(N-1)+35*u0(N))/(12*dr^2);

%===============================================

%==================================================

FH(1)=0.;
SH(1)=(-h0(3)+16*h0(2)-30*h0(1)+16*h0(2)-h0(3))/(12*dr^2);
SH(2)=(-h0(2)+16*h0(1)-30*h0(2)+16*h0(3)-h0(4))/(12*dr^2);
FH(2)=(h0(2)-8*h0(1)+8*h0(3)-h0(4))/(12 *dr);

for i=3:N-2
SH(i)=(-h0(i-2)+16*h0(i-1)-30*h0(i)+16*h0(i+1)-h0(i+2))/(12*dr^2);
FH(i)=(h0(i-2)-8*h0(i-1)+8*h0(i+1)-h0(i+2))/(12*dr);

end
FH(N-1)=(-h0(N-4)+6*h0(N-3)-18*h0(N-2)+10*h0(N-1)+3*h0(N))/(12*dr);
FH(N)=(3*h0(N-4)-16*h0(N-3)+36*h0(N-2)-48*h0(N-1)+25*h0(N))/(12*dr);
SH(N-1)=(-h0(N-4)+4*h0(N-3)+6*h0(N-2)-20*h0(N-1)+11*h0(N))/(12*dr^2);
SH(N)=(11*h0(N-4)-56*h0(N-3)+114*h0(N-2)-104*h0(N-1)+35*h0(N))/(12*dr^2);
% plot(SH)
% return
%===============================================
F(1)=-SH(1);
for i=2:N
F(i)=(u0(i)*u0(i)/2)-0.5*(SH(i)+(1/r(i))*FH(i));
end
F(N)=(u0(N)*u0(N)/2)-2;
% plot(F)
% return
%=========================================================

%résolution de l'équation de Navier et Stokes(vitesse radialur(i))

u1(1)=0;

for i=2:N-1
u1(i)=u0(i)-0.5*(dt/dr)*(F(i+1)-F(i))+0.25*(dt*dt/dr^2)*(u0(i)*FU(i)*(u0(i+1)+u0(i))*(F(i+1)-F(i))-(u0(i)*FU(i)*(u0(i)+u0(i-1)))*(F(i)-F(i-1)));
end
u1(N)=UHB/h0(N);



%plot(r,u1)
%return
%==============================================================
SU(1)=-30*u1(1)/(12*dr^2);
SU(2)=(u1(2)+16*u1(1)-30*u1(2)+16*u1(3)-u1(4))/(12*dr^2);
FU(1)=(-u1(3)+8*u1(2)+8*u1(2)-u1(3))/(12*dr);
FU(2)=(-u1(2)-8*u1(1)+8*u1(3)-u1(4))/(12*dr);
for i=3:N-2
SU(i)=(-u1(i-2)+16*u1(i-1)-30*u1(i)+16*u1(i+1)-u1(i+2))/(12*dr^2);
FU(i)=(u1(i-2)-8*u1(i-1)+8*u1(i+1)-u1(i+2))/(12*dr);
end
FU(N-1)=(-u1(N-4)+6*u1(N-3)-18*u1(N-2)+10*u1(N-1)+3*u1(N))/(12*dr);
FU(N)=(3*u1(N-4)-16*u1(N-3)+36*u1(N-2)-48*u1(N-1)+25*u1(N))/(12*dr);
SU(N-1)=(-u1(N-4)+4*u1(N-3)+6*u1(N-2)-20*u0(N-1)+11*u1(N))/(12*dr^2);
SU(N)=(11*u1(N-4)-56*u1(N-3)+114*u1(N-2)-104*u1(N-1)+35*u1(N))/(12*dr^2);


% plot(SU)
% return

%============================================================
%resolution de l'équation de continuité(épaisseur du film h)

h1(1)=h0(1)-2*dt*h0(1)*FU(1);

for i=2:N-1
AA=h0(i)*u1(i)/r(i)+h0(i)*FU(i)+u1(i)*FH(i);
BB=2*u1(i)*FH(i)/r(i)+2*h0(i)*FU(i)/r(i)+2*FU(i)*FH(i)+h0(i)*SU(i)+u1(i)*SH(i);
h1(i)=h0(i)-dt*AA*(1-dt*FU(i)/2)+dt*dt*BB*u1(i)/2;
end

h1(N)=h0(N)-dt;
% plot(r,h1)
% return
%===============================================

h0=h1;
u0=u1;
%=========================================

T_sim=T_sim+ dt
end
figure(3);
plot(r,u1,'BLACK');
xlabel ('r');
ylabel('u1');
title('evolution u0 au cours de simulation ');
hold on
figure(4);
semilogy(r,h1,'black');
xlabel('r');
ylabel('h0');
title('evolution h0 au cours du simulation T_sim');
hold on
****************************************************
MERCI