可压流的龙格-库塔时间步进求解法及其应用。

裴来凡

发布于 2022-05-28 15:37:28

3060

eulerstep_stag_com.m

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

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

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

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;