transport_galerkin1d_upwind.m

１Ｄ非连续性力学分析

%******************************
%    Galerkin dicontinu
%*****************************


clear all
close all


%*** declaration des variables

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


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

N=200; % temps total 0<n<N
p=1;
dt=p*0.1*dx/v; %increment de temps
%dt=30;

%*******************************

Me=[(13*dx/35) (9*dx/70);(9*dx/70) (13*dx/35)];
Qe=[(-0.5) (-0.5);(0.5) (0.5)];
CLe=[-1 0;0 1];
%%%%%%%%Assemblage%%%%%%%%%%%%%%%%%%%%a la fin%%%%%%%%%%%%%%%%%
M(1:(I+1),1:(I+1))=zeros;
for i=1:I
    M(i:(i+1),i:(i+1))=M(i:(i+1),i:(i+1))+Me;
end

Q(1:(I+1),1:(I+1))=zeros;
for i=1:I
    Q(i:(i+1),i:(i+1))=Q(i:(i+1),i:(i+1))+Qe;
end

CL(1:(I+1),1:(I+1))=zeros;
for i=1:I
    CL(i:(i+1),i:(i+1))=CL(i:(i+1),i:(i+1))+CLe;
end

invM=inv(M);
U=invM*Q;
V=invM*CL;

%initialisation

phi=sparse(zeros(I,N));
for i=1:(I+1)
    x(i)=(i-1)*dx;
    phi(i,1)=1-x(i);
end

% r閟olution%
for n=1:N
        phi(:,n+1)=v*dt*(U-V)*phi(:,n)+phi(:,n);
end

for n=1:N
    for i=1:I+1
        phiex(i,n)=1-x(i)+v*n*dt;
    end
end

%recherche du front stabilite et diffusion
for n=1:N
    for i=1:I+1
        if (phi(i,n)<0)
            k=i;
            a=(phi(k-1,n)-phi(k,n))/(x(k-1)-x(k));
            b=phi(k,n)-a*x(k);
            xk(n)=-(b/a);
            x_exacte(n)=1+v*n*dt;
            %stabilit