eulerstep_stag_com.m
代码语言：javascript复制totenthalpy = etotstar pstar; totenthL = extrap(totenthalpy, limtype);
rhoL = extrap(rhostar, limtype);
etotstar(1) = etotold(1) - rkalpha*lambda*(ustar(2)*totenthL(1)...
- uleft*(gamma*pleft/(gamma-1) 0.5*rholeft*uleft^2));
rhostar(1) = rhoold(1) - rkalpha*lambda*(ustar(2)*rhoL(1) - uleft*rholeft);
for j = 2:J-1
etotstar(j) = etotold(j) - rkalpha*lambda*(ustar(j 1)*totenthL(j)...
- ustar(j)*totenthL(j-1));
rhostar(j) = rhoold(j) - rkalpha*lambda*(ustar(j 1)*rhoL(j)...
- ustar(j)*rhoL(j-1));
end
um = ustar.*mstar; umL = extrap(um, limtype);
ustar(1) = uleft; mstar(1) = rholeft*uleft;
for j = 2:J-1
mstar(j) = mold(j) - rkalpha*lambda*(umL(j) - umL(j-1))...
- rkalpha*lambda*(pstar(j) - pstar(j-1));
ustar(j) = 2*mstar(j)/(rhostar(j) rhostar(j-1));
end
mstar(J) = mold(J) - rkalpha*lambda*(umL(J) - umL(J-1))...
- rkalpha*lambda*(pright - pstar(J-1));
ustar(J) = 2*mstar(J)/(rhostar(J-1) rhoright);
for j = 1:J-1
pstar(j) = (gamma-1)*(etotstar(j) - 0.25*(mstar(j)*ustar(j) ...
mstar(j 1)*ustar(j 1)));
end
eulerstep_stag_unif.m
代码语言：javascript复制allam = rkalpha*lambda;
rhoL = extrap(rhostar, limtype);
rhostar(1) = rhoold(1) - allam*(ustar(2)*rhoL(1) - uleft*rholeft);
for j = 2:J-1
rhostar(j) = rhoold(j) - allam*(ustar(j 1)*rhoL(j) - ustar(j)*rhoL(j-1));
end
um = mstar.*ustar; umL = extrap(um, limtype);
mstar(1) = rholeft*uleft;
for j = 2:J-1
mstar(j) = mold(j) - allam*(umL(j) - umL(j-1) pold(j) - pold(j-1));
end
mstar(J) = mold(J) - allam*(um(J) - umL(J-1) pright - pold(J-1));
rhowall(1) = rholeft; rhowall(J) = rhostar(J-1);
rhoL = extrap(rhostar, limtype);
for j = 2:(J-1)
rhowall(j) = 0.5*(rhostar(j-1) rhostar(j));
end
ustar = mstar./rhowall;
extrap.m
代码语言：javascript复制function zL = extrap(z, limtype); % Slope-limited extrapolation (MUSCL)
% Staggered scheme
% Function called: limiter
% Numbering for pressure nodes:
% 1 1 2 2 J J
% |-o-|-o-|-o-|-o-|-o-|-o-|-o-|-o-|
% Numbering for momentum nodes:
% 1 1 2 2 J-1 J
% |-o-|-o-|-o-|-o-|-o-|-o-|-o-|---o---|
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));
end
zL(J) = z(J);
f.m
代码语言：javascript复制function y = f(p34) % f is the shocktube relation (10.51)
global PRL CRL MACHLEFT gamma
wortel = sqrt(2*gamma*(gamma-1 (gamma 1)*p34));
yy = (gamma-1)*CRL*(p34-1)/wortel;
yy = (1 MACHLEFT*(gamma-1)/2-yy)^(2*gamma/(gamma-1));
y = yy/p34 - PRL;
limiter.m
代码语言：javascript复制function y = limiter(a,b,limtype) % Called by extrap
% Theory in Section 4.8 of:
% P. Wesseling: Principles of Computational Fluid Dynamics
% Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html
if abs(b) < 1/10^9, y=0;
else
if limtype == 3 % superbee
r = a/b; y = max(0, max(min(2*r,1), min(r,2)));
elseif limtype == 2 % van Albada
y = (a^2 a*b)/(a^2 b^2);
elseif limtype == 1 % Minmod (Roe(1986))
y = max(0, min(a/b,1));
elseif limtype == 4 % PL limiter
r = a/b;
if r<0, y = 0;
elseif r<0.4, y = 2*r;
elseif r<4, y = (2 r)/3;
else, y = 2;
end
elseif limtype == 5 % Chakravarty/Osher limiter (Kr"oner (1997) p. 109)
r = a/b; y = max(0, min(r,1.1));
else % no MUSCL: first order upwind
y = 0;
end
end
