范阿尔巴达限制器图包括asic激波管关系式等。

figure(1), clf, hold on

x = 0.0:0.01:4.0; y = (x.*x + x)./(1 + x.*x);  % Graph of van Albada limiter
plot(x,y)

x = 0.0:0.01:1.1; y = 0.5*(1+sqrt(2))*x; plot(x,y)
x = 0.0:0.01:4.0; y = 0.5*(1+sqrt(2))*ones(size(x)); plot(x,y)
y = ones(size(x)); plot(x,y)
x = 0.0:0.01:1.3; plot(x,x)
y = 0.0:0.01:1; x = ones(size(y)); plot(x,y,'--')

function[zL,zR] = extrap(z, limtype);  % Slope-limited extrapolation
          % Function called: limiter 
J = length(z);
zL(1) = z(1); 
for j = 2:J-1
  b = (z(j+1) - z(j));   a = (z(j) - z(j-1));
  zL(j)   = z(j) + 0.5*limiter(a, b, limtype)*(z(j+1)-z(j));
  zR(j-1) = z(j) + 0.5*limiter(b, a, limtype)*(z(j-1)-z(j));
end
zR(J-1) = z(J);

function y = limiter(a,b,limtype)
% Limiter function as defined in Sect. 4.8 and as used in Sect. 10.8
if  b==0,  y=0; else
  if limtype == 2            % van Albada
    y = (a^2 + a*b)/(a^2 + b^2);
  elseif limtype == 1    % Minmod
    y = max(0, min(a/b,1));
  else        % No limiting, no MUSCL
    y = 0; 
  end   
end

clear all
global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
  uright  tend  epsi  lambda    % lambda = dt/dx

    % .....................Input............................
gamma = 1.4;   % Ratio of specific heats
J = 48;    % Number of grid cells
limtype = 2;  % Enter 1 for minmod or 2 for van Albada or something else
    %  if MUSCL not wanted

gammab = 1/(gamma - 1); gam1 = gamma-1; gamgam = gamma^gamma;
problem_specification  

h = 1/J;          % Cell size  
dt = lambda*h;        % Time step
n = floor(tend/dt);      % Number of time-steps

%     Definition of grid numbering 
%       x=0               x=1
% grid   |---o---|---o---|---o---  ...  --|---o---|
%            1   1   2   2   3           J-1  J 

xcenter = h*[1:J] - h/2;    % Location of cell centers

press = zeros(size(xcenter));    % Preallocation of pressure, 
rhoold = press; uold = press;    %       density, velocity,
rhonew = press; mnew = press;    %  momentum, 
totenew = press;      %  total energy

for j = 1:length(xcenter)    % Initial conditions
  if xcenter(j) < 0.5, press(j) = pleft; rhoold(j) = rholeft;  uold(j) = uleft;
  else,          press(j) = pright; rhoold(j) = rhoright; uold(j) = uright;
  end
end

  % toten = rho*(total energy)
totenold = rhoold.*(0.5*uold.*uold + gammab*press./rhoold);
totenleft = totenold(1); totenright = totenold(J);
mold = rhoold.*uold;      % m is momentum
c = sqrt(gamma*press./rhoold);    % Sound speed
mach = uold./c;
enthalpy =  0.5*uold.*uold + gammab*c.^2;

    % Preallocation of auxiliary variables
machplus = mach; machminus = mach;
presplus = mach; presminus = mach;

    % Preallocation of Liou-Steffen fluxes
flux1 = zeros(J-1,1); flux2 = flux1; flux3 = flux1;

    % Preallocation of extrapolated states
Z1L = zeros(J-1,1); Z2L = Z1L; Z3L = Z1L;
Z1R = Z1L; Z2R = Z1L; Z3R = Z1L;
U1L = Z1L; U1R = Z1L; U2L = Z1L; U2R= Z1L;
U3R = Z1L; U3L = Z1L;

flopcount = flops;      % Operations counter
t = 0;
for i = 1:n,  t = t + dt;
  rhostar = rhoold; mstar = mold; totenstar = totenold;
  rkalpha = 0.25;  listeeulerstep
  rkalpha = 1/3;   listeeulerstep
  rkalpha = 0.5;   listeeulerstep
  rkalpha = 1;     listeeulerstep
  
  rhonew = rhostar; mnew = mstar; totenew = totenstar;
  
  uold = mnew./rhonew; press = gam1*(totenew - 0.5*mnew.*uold);
  rhoold = rhonew; totenold = totenew; mold = mnew;
  c = sqrt(gamma*press./rhoold); mach = uold./c;
end

flopcount = flops - flopcount;
entropy = log(press./rhoold.^gamma);

figure(1), clf
subplot(2,3,1),hold on,title('DENSITY','fontsize',14),plot(xcenter,rhonew,'o')
subplot(2,3,2),hold on,title('VELOCITY','fontsize',14),plot(xcenter,uold,'o')
subplot(2,3,3),hold on,title('PRESSURE','fontsize',14),plot(xcenter,press,'o')
subplot(2,3,4),hold on,title('MACHNUMBER','fontsize',14),plot(xcenter,mach,'o')
subplot(2,3,5),hold on,title('ENTROPY','fontsize',14),plot(xcenter,entropy,'o')
subplot(2,3,6), axis('off'), hold on, title('Liou-Steffen scheme','fontsize',14)
text(0,0.9,['lambda = ', num2str(lambda),'  t = ',num2str(n*dt)])
text(0,0.75,' Primitive extrapolation'),  text(0,0.6,'SHK Runge-Kutta')
if limtype == 0,    s1 = 'No MUSCL';
elseif limtype == 1,    s1 = 'MUSCL, minmod limiter';
else,        s1 = 'MUSCL, van Albada limiter';
end
text(0,0.45,s1)
text(0,0.3,['flopcount = ',num2str(flopcount)])

Riemann    % Plot exact solution