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;
