transport_2d_upwind_sin.m

２Ｄ非连续性力学分析

close all
clear all

%*** declaration des variables

L=6
; % 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=500; % temps total 0<n<N

J=50; % nombre element en espace 0<j<J
dy=L/(J-1); % increment en espace
dt=0.1*dx/v; %increment de temps


%%%%%%%%%%%%%%%%%%%%***************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)=sin(1-x(i))+cos(y(j));
    end
end

for n=2:N
    for j=1:J
        phi(1,j,n)=sin(1-x(1)+v*n*dt)+cos(y(j)+v*n*dt);
    end
    for i=2:I
        phi(i,1,n)=sin(1-x(i)+v*n*dt)+cos(y(1)+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');
% hold on
% tmp=zeros(I,J);
% 
% mesh(tmp)
%aviobj = avifile('example.avi')
for n=1:10:N
    view([-J,-I,N])
    surfl(phi(:,:,n))
    %shading interp;
    shading faceted;
    colormap(pink)
    % axis xy
    pause
    %frame = getframe(gca);
    %aviobj = addframe(aviobj,frame);
end
%aviobj = close(aviobj);
%close all


for n=1:10:N
%     view([-J,-I,N])
%     surfl(phi(:,:,n))
    %shading interp;
%     shading faceted;
    colormap (hsv)
    [c,h] = contour(phi(:,:,n));
    h = clabel(c,h);
    pause
end