Buckley-Leverett方程Riemann问题的解。

clear all

%*****************************INPUT********************************************
geval = 2;  % = 1:     Initial condition: shock: example 9.4.4
    % not = 1: Initial condition: fan:   example 9.4.5
tend = 0.5;  % Final time
lambda = 0.75;  % lambda = dt/dx
J = 40;   % Number of nodes. 
%******************************************************************************

%*********************Uniform grid*********************
%                  #
%  o------o------o------o------o------o------o   #
%  1      2                                  J   #
%      x=-1                                      x=1  #
%******************************************************

dx = 2/(J-1); dt = lambda*dx;     % Stepsizes
ntime = floor(tend/dt); %ntime = 0;     % Number of time steps
y = [1:J]*dx - dx -1;      % Location of grid points

u = y;        % Preallocation of solution
if geval == 1, uleft =   1; uright = 0;     % Initial and boundary conditions 
else,           uleft = - 1; uright = 1;
end
for j = 1:J,
  if y(j) < -10*eps,    u(j) = uleft;
  elseif y(j) > 10*eps, u(j) = uright;
  else,      u(j) = 0.5*(uleft+uright);
  end 
end 
uinit = u;

for n = 1:ntime    %   Nonconservative first order upwind scheme (9.75)
  uold = u;
  j = 1; u(j) = uold(j) - 0.5*lambda*...
       ((uold(j) + abs(uold(j)))*  (uold(j) - uleft) + ...
        (uold(j) - abs(uold(j)))*(uold(j+1) - uold(j)));       
  for j = 2:J-1,  u(j) = uold(j) - 0.5*lambda*...
       ((uold(j) + abs(uold(j)))*  (uold(j) - uold(j-1)) + ...
        (uold(j) - abs(uold(j)))*(uold(j+1) - uold(j)));       
  end
  j = J; u(j) = uold(j) - 0.5*lambda*...
     ((uold(j) + abs(uold(j)))*(uold(j) - uold(j-1)) + ...
     (uold(j) - abs(uold(j)))* (uright - uold(j)));
end              
figure(1); clf, subplot(221), hold on, plot(y,u,'o')
title('Nonconservative scheme','fontsize',14)
 
u = uinit;
for n = 1:ntime    % Courant-Isaacson-Rees scheme (9.81)
  uold = u;
  j = 1;  s1 = sign(uold(j) + uold(j+1)); s2 = sign(uold(j) + uleft);
  u(j) = uold(j) - 0.25*lambda*((1+s1)*uold(j)^2 + ...
      (1-s1)*uold(j+1)^2 - (1+s2)*uleft^2 - (1-s2)*uold(j)^2);
  for j = 2:J-1
    s1   = sign(uold(j) + uold(j+1)); s2 = sign(uold(j) + uold(j-1));
    u(j) = uold(j) - 0.25*lambda*((1+s1)*uold(j)^2 + ...
       (1-s1)*uold(j+1)^2 - (1+s2)*uold(j-1)^2 - (1-s2)*uold(j)^2);
  end
  j = J; s1 = sign(uold(j) + uright); s2 = sign(uold(j) + uold(j-1));
  u(j) = uold(j) - 0.25*lambda*((1+s1)*uold(j)^2 + (1-s1)*uright^2 -...
      (1+s2)*uold(j-1)^2 - (1-s2)*uold(j)^2);
end              
subplot(222), hold on, plot(y,u,'o')
title('Courant-Isaacson-Rees scheme','fontsize',14)

u = uinit;
for n = 1:ntime    % Lax-Friedrichs scheme (9.80)
  uold = u;
  j = 1;      
  u(j) = 0.5*(uleft + uold(j+1)) - 0.25*lambda*(uold(j+1)^2 - uleft^2);
  for j = 2:J-1
    u(j) = 0.5*(uold(j-1)+uold(j+1)) - 0.25*lambda*(uold(j+1)^2 - uold(j-1)^2);
  end
  j = J;      
  u(j) = 0.5*(uold(j-1) + uright) -  0.25*lambda*(uright^2 - uold(j-1)^2);
end              
subplot(223), hold on, plot(y,u,'o')
title('Lax-Friedrichs scheme','fontsize',14)

u = uinit;  fplus = zeros(1,J); fminus = fplus;  
for n = 1:ntime    % Engquist-Osher scheme (9.82)
  uold = u;
  for j = 1:J
    if uold(j) > 0,   fplus(j) = 0.5*uold(j)^2; fminus(j) = 0;
    else,     fplus(j) = 0; fminus(j) = 0.5*uold(j)^2;
    end
  end
  if uleft > 0,     fplusleft = 0.5*uleft^2;
  else,       fplusleft = 0; 
  end
  if uright > 0,   fminusright = 0;
  else,       fminusright = 0.5*uright^2;
  end
  j = 1; 
  u(j) = uold(j) - lambda*(fplus(j) + fminus(j+1) - fplusleft - fminus(j));  
  for j = 2:J-1   
    u(j) = uold(j) - lambda*(fplus(j) + fminus(j+1) - fplus(j-1) - fminus(j));
  end
  j = J;
  u(j) = uold(j) - lambda*(fplus(j) + fminusright - fplus(j-1) - fminus(j));    
end              
subplot(224), hold on, plot(y,u,'o')
title('Engquist-Osher scheme','fontsize',14)

  % Exact solution
zz = -1:0.01:1; uexact = zeros(size(zz)); t = n*dt;
if geval == 1
  for j = 1:length(zz),
    if zz(j) < t/2,    uexact(j) = 1;
    else,    uexact(j) = 0;
    end
  end  
else
  for j = 1:length(zz),
    if zz(j) < - t,  uexact(j) = -1;
    elseif zz(j) > t,  uexact(j) = 1;
    else,    uexact(j) = zz(j)/t;
    end
  end  
end
for k = 1:4, subplot(2,2,k), plot(zz,uexact), end
%print solution -depsc;  % Figure is saved in solution.eps for printing

a = 0.5;  % Parameter in flux function on page 371
    % ............End of input.....................................
x = 0.0:0.05:1.0; den = x.*x + 0.5*(1-x).^2; f = (x.*x)./den;
figure(1), clf, plot(x,f,'-')

a = 0.5;  % Parameter in flux function on page 371
    % ............End of input.....................................

x = 0.0:0.05:1.0; den = x.*x + 0.5*(1-x).^2; 
f = (x.*x)./den;  % Buckley-Leverett flux function

x1 = 1/sqrt(3); 
f1 = (x1*x1)/(x1*x1 + 0.5*(1-x1)^2);  %Point where tangent meets curve

figure(1), clf, hold on, plot(x,f,'-')
xx = [0  1];  xx = x1*xx; yy = [0  1]; yy = f1*yy; plot(xx,yy,'-')
zz = [0  f1]; xx = [x1 x1];         plot(xx,zz,'--')
 

figure(1), clf, hold on
  % Left constant state
x = - 0.5:0.05:0; phi = ones(size(x)); plot(x,phi,'-')

  % Right constant state
x1 = 0.5 + 0.5*sqrt(3); x = x1 :0.05 :2; phi = zeros(size(x)); plot(x,phi,'-')

  % Expansion fan
phi = 0:0.05:1; phi = phi*(1 - 1/sqrt(3));  phi = 1/sqrt(3) + phi;
x = (phi - phi.^2)./((phi.^2 + 0.5*(1-phi).^2).^2); plot(x,phi,'-')
x = [x1   x1]; phi = [0   1/sqrt(3) ]; plot(x,phi,'-')