transport_2d_upwind_plan.m

２Ｄ非连续性力学分析

%*******************
%  transport 1D
%*******************

clear all
close all

%*** declaration des variables

L=50; % longueur de la poutre
v=0.02; % vitesse constante

I=50; % nombre element en espace 0<i<I
dx=L/(I-1); % increment en espace

N=200; % temps total 0<n<N

J=50; % nombre element en espace 0<j<J
dy=L/(J-1); % increment en espace

% dt=0.9*sqrt(dx*dx+dy*dy)/v; %increment de temps
dt=0.1*dx/v; %increment de te

%***************initialisation****
phi=zeros(I,J,N);

for i=1:I
    x(i)=(i-1)*dx;
end
for j=1:J
    y(j)=(j-1)*dy;
end

for i=1:I
    for j=1:J
        phi(i,j,1)=1-x(i)-y(j);
    end
end
%***** boucle sur le temps

for n=2:N
    for j=1:J
        phi(1,j,n)=1-x(1)-y(j)+2*v*n*dt;
    end
    for i=2:I
        phi(i,1,n)=1-x(i)-y(1)+2*v*n*dt;
        
    end
end



%boucle sur temps n>2
for n=2:N
    
    for i=2:I
        for j=2:J
        phi(i,j,n)=-(v*(dt/dx))*(phi(i,j,n-1)-phi(i-1,j,n-1))-(v*(dt/dy))*(phi(i,j,n-1)-phi(i,j-1,n-1))+phi(i,j,n-1);
        end
    end
end

 title('Transport de la surface initiale')
 xlabel('I');
 ylabel('J');
 zlabel('Phi');
for n=1:5:N
hold on
view([J,-I,N])
tmp=zeros(I,J);
surfl(phi(:,:,n))
mesh(tmp)
shading interp;
colormap(pink);
pause
end
% close


%  for n=1:N
%  [c,h] = contour(phi(:,:,n));
% h = clabel(c,h);
% pause(0.05)
% end
% close