alkairaouni
Messages postés1Date d'inscriptionmercredi 21 novembre 2012StatutMembreDernière intervention21 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
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
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