1D Unsteady Compressible Navier

In the current project a Shock wave tube or a Riemann tube problem is being studied.This is aparticularly useful case as it simulates a shock wave, expansion wave and discontinuity.From a numerical point of view, this problem constitutes, since the exact solution is known, aninevitable and difficult test case for any numerical method dealing with discontinuous solutions.Finally, there is a practical interest, since this model is used to describe real shock tubeexperimental devices.

The physical description and the method used to calculate the numericalsolution is discussed in the Introduction part of the report. The report concludes by comparing theanalytical and numerical solution.

The fundamental idea of the shock tube is the following: consider a long one-dimensional (1D)tube, closed at its ends and divided into two equal regions by a thin diaphragm (Fig 1). Each regionis filled with the same gas, but with different thermodynamic parameters (pressure, density, andtemperature). The region with the highest pressure is called the driver section of the tube, while thelow-pressure part is the working section. The gas being initially at rest, the sudden breakdown ofthe diaphragm generates a high-speed flow, which propagates in the working section.

Solution Method:

The MacCormack method is a widely used finite difference scheme that works well for solvingnonlinear differential equations suited for marching solutions. The MacCormack is derived from theLax-Wendroff technique. MacCormack method is a simplified form of the Lax-Wendroff in a waythat the second order derivatives are replaced by a representative mean of first order derivatives attime t and t + Δt.

By introducing the average values, the lengthy algebra of calculating values forsecond order derivatives is circumvented still preserving the second order of accuracy for thesolutions. In this work, the MacCormack method has been implemented in order to solve the shocktube problem.

General Idea:

Results discussion:

[I] Analyzing results for different spatial stepsNumerical calculation using MacCormack method is calculated for different number of spatialsteps, and analyzed its effect on solution accuracy. Time step is calculated based on the stabilitycriteria required for each time iteration to prevent solution divergence. Hence for stability,numerical simulation is time adaptive.

Results obtained, suggests that increasing the number of spatial steps, i.e. reducing step size,improves the accuracy of the solution.

These results were expected because, reducing the the stepsize, will reduce the truncation error i.e. error due to higher order terms. If dx is large, then thetruncation error will be large hence, the solution will become less accurate.But, apart from accuracy, one has to take into account of - stability of solution.

To dampen thegrowth of error with each iteration, scheme must satisfy stability criteria, so that solution proceedsin a stable manner, otherwise error will progressively become larger.Hence for satisfying stability criteria as discussed in previous section, time step size decreases aswell, as we decrease the step size. Hence lower the spatial step size, more are the number of timeiterations, and hence it takes more time to complete the numerical simulations.

Comparing the effect of order of numerical scheme: 4rth order and 2nd order

Results:

After examining results for 2nd and 4rth order, it was observed that 4rth order was more accurateand has less wiggles in the numerical solution as compared to 2nd order, which is quite visible intemperature and velocity plot.

[III] Effect of artificial viscosity on solution’s accuracy and stability (by varying Cx)Comparisons are carried out spatial steps n=101 and 4rth order scheme:

As discussed in previous section, artificial viscosity, is added to the numerical solution to make itmore stable. So as observed in the plots, for the lowest value for Cx, numerical solution haswiggles, i.e. tending towards instability, and if further the value of Cx is reduced, solution becomesunstable.Now, though adding artificial viscosity, makes the solution stable but it decreases the accuracy ofsolution which is quite observable in velocity and temperature plot for Cx-0.2 &0.3.

After several experimentation of Cx values, value of 0.1 was chosen as the optimum for stabilityand accuracy of solution.

Conclusion:

Based on the numerical analysis carried out using Maccormack method, following conclusions weremade:

1. It was observed, as the grid size or spatial step size is decreased, solution becomes moreaccurate. This is because, the truncation error for lesser grid size is much lower as whencompared to higher grid sizes. Hence, smaller the grid size, lower is the truncation error,higher is the accuracy.

2. As, the grid size was decreased in the numerical scheme, computation time increased. Thisis because, with lower grid size, time step should be proportionately reduced to satisfy thestability criteria, unless the error will grow with each iteration, making the solution unstable.Hence, lower the grid size, lesser is the time step size, hence computation time increases.

3. Compared to a 2nd order scheme, 4rth order scheme is found to be more accurate. This wasexpected because, the truncation error is much lower in the 4rth order scheme, makingsolution more accurate and closer to the analytical solution

4. Artificial viscosity was varied using Cx. It was observed that lower the value of Cx, lower isthe stability of solution, and higher the value of Cx, higher is the stability but lower is theaccuracy. So, basically, value of Cx, is calculated experimentally to achieve the optimalstate of solution which is stable as well as close to the analytical solution.

Matlab Code:

%%%%%%%%%%%%%%%%% Analytical solution at  time instance-0.06.%%%%%%%%%%%%%%

 clear all;close all;clc
 
 ta = 0.006;  % Analytical time instance
 
%  Parameters:
 cp  = 1.005e3;
 cv  = 0.718e3;
 gamma   = cp/cv;
 alpha   = (gamma+1)/(gamma-1);
 Pr  = 0.707;
 R   = 286.9;
 x0  = 5;
 PR = 3.031;
 
 P4_a      = 1e5;
 P1_a    = 1e4;
 rho4_a   = 1;
 rho1_a   = 0.125;
 T1_a      = P1_a/rho1_a/R;
 T4_a    = P4_a/rho4_a/R;
 a1_a     = sqrt(gamma * R * T1_a);
 a4_a     = sqrt(gamma * R * T4_a);
 u1_a     = 0;
 u4_a     = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%calculate analytical results%%%%%%%%%%%%%%%%%%%%%
 Na   = 1000;%grid steps for analytical
 xx  	= linspace(0,10,Na);
 Pa    = zeros(size(xx));
 Ta    = zeros(size(xx));
 ua    = zeros(size(xx));
 rhoa    = zeros(size(xx));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Region 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 P2_a      = P1_a*PR;
 rho2_a    = rho1_a*(1+alpha*PR)/(alpha+PR);
 T2_a      = P2_a/rho2_a/R;
 u2_a     = u1_a + a1_a*(PR-1)/(1+alpha*PR)^(1/2)/sqrt(gamma*(gamma-1)/2);
 
%%%%%%%%%%%%%%%%%%%% Contact discontinuity location <-- u2%%%%%%%%%%%%%%%%
 xcont   = x0 + u2_a*ta;
 
%%%%%%%%%%%%%%%%%%%% Shock wave location <-- u2%%%%%%%%%%%%%%%%%%%%%%%%%%%
 W_a      = u1_a + a1_a*(PR-1)*a1_a/gamma/(u2_a-u1_a);
 xshock  = x0 + W_a*ta;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Region 3:%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 P3_a = P2_a;
 u3_a = u2_a;
 rho3_a    = (P3_a/P4_a*rho4_a^gamma)^(1/gamma);
 a3_a     = a4_a - (gamma-1)/2*u3_a;
 T3_a     = a3_a^2 / gamma / R;
 
%%%%%%%%%%%%%%%%%%%%%%%Expansion wave location <-- u3%%%%%%%%%%%%%%%%%%%%%
 xl      = x0 - a4_a*ta;        % left boundary 
 xr      = x0 + (u3_a-a3_a)*ta;   % right boundary

%%%%%%%%%%%%%%%%%%%%%%% Expansion wave region:%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 xexp 	= xx(xx>=xl & xx<xr);
 uexp 	= 2/(gamma+1)*(a4_a+(xexp-x0)/ta);
 Texp  	= T4_a*(1-(gamma-1)/2*uexp/a4_a).^2;
 rhoexp  = rho4_a*(1-(gamma-1)/2*uexp/a4_a).^(2/(gamma-1));
 Pexp    = rhoexp.*Texp * R;
 
 Pa(xx>=xshock) = P1_a;    % Region 1
 Ta(xx>=xshock) = T1_a;
 ua(xx>=xshock) = u1_a;
 rhoa(xx>=xshock) = rho1_a;
 Pa(xx>=xcont & xx<xshock) = P2_a; % Region 2
 Ta(xx>=xcont & xx<xshock) = T2_a;
 ua(xx>=xcont & xx<xshock) = u2_a;
 rhoa(xx>=xcont & xx<xshock) = rho2_a;
 Pa(xx>=xr & xx<xcont) = P3_a;     % Region 3
 Ta(xx>=xr & xx<xcont) = T3_a;
 ua(xx>=xr & xx<xcont) = u3_a;
 rhoa(xx>=xr & xx<xcont) = rho3_a;
 Pa(xx<xl) = P4_a;                 % Region 4
 Ta(xx<xl) = T4_a;
 ua(xx<xl) = u4_a;
 rhoa(xx<xl) = rho4_a;
 Pa(xx>=xl & xx<xr)     = Pexp;  % Expansion wave region
 ua(xx>=xl & xx<xr)     = uexp;
 Ta(xx>=xl & xx<xr)     = Texp;
 rhoa(xx>=xl & xx<xr)     = rhoexp;
 
 aa  = sqrt(gamma * R * Ta);
 Ma  = ua./aa;
 sa  = Pa./rhoa.^gamma;
 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Numerical Solution-Maccormack scheme%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%given values
R=286.9;
gamma=1.4;
Pr=0.707;
gamma_Pr=gamma/Pr;
cp=1.005e3;
cv=0.718e3;
Tf=298;
muf=1.846e-5;
R_cfl=0.7;

%Grid
n=101;
x=linspace(0,10,n);
dx=x(2)-x(1);

%Setting array 
U=zeros(3,n);
F=zeros(3,n);
J=zeros(3,n);

%Initial conditions
%setting for pressure, temperature, and density
P=ones(1,n);
rho=ones(1,n);
T=ones(1,n);
for i=1:((n-1)/2)+1 %Pressure,Density,Temperatute setting<=x0
    P(i)=1e5;
    rho(i)=1;
    T(i)=P(i)/(rho(i)*R);
end
for j=((n-1)/2)+2:n  %Pressure,Density,Temperatute setting>=x0
    P(j)=1e4;
    rho(j)=0.125;
    T(j)=P(j)/(rho(j)*R);
end

%BC conditions
u=zeros(1,n);
T(1)=2*T(2)-T(3); %Temperature
T(n)=2*T(n-1)-T(n-2);

P(1)=2*P(2)-P(3);   %Pressure
P(n)=2*P(n-1)-P(n-2);

e=cv*T; %Internal energy

[Uini,Fini,Jini]=Code(rho,u,P,T,n,e,dx); %setting initial values of Uini,Fini,Jini,Sini

%[rho1,u1,P1,T1,e1]=Decode(Uini,n); %To Check Decode Function

%time step-calculating time step based on given dx for stability
% mu=muf*(Tf+110)/((Tf)^1.5)*((T).^1.5)./(T+110);
% visc=mu./rho;
% visc_d=max((4/3)*visc,visc.*gamma_Pr);
% sonic=sqrt(gamma*R*T);
% dt_cfl=((abs(u)./dx)+(sonic./(dx)^2)+((2*visc_d)./(dx)^2)).^-1;
% dt=R_cfl*min(dt_cfl); %dt optimized  for given dx
  

nt=1e5; %time steps
dt=1e-4;
%  for Initial Predictor values
Up=Uini; 
Fp=Fini;
Jp=Jini;
%Sp=Sini;
U_bar=zeros(3,n);

%%%%%%%%%%%%%%%%%%%%%Outer loop- time loop%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=1:nt     
    
    %%%%%%%%%%%%%%%%%%%%%%%%%Predictor loop%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%% 2nd order in space%%%%%%%%%%%%%%%%%%%%%
    U_bar(:,1)=Up(:,1)+dt*(Jp(:,1)-(Fp(:,2)-Fp(:,1))/dx);%Sp(:,1);
    for j=2:n-1
        U_bar(:,j)=Up(:,j)+dt*(Jp(:,j)-(Fp(:,j+1)-Fp(:,j))/dx);%+Sp(:,j);
    end
    U_bar(:,n)=Up(:,n)+dt*(Jp(:,n)-(Fp(:,n)-Fp(:,n-1))/dx);%+Sp(:,n);
    
    %%%%%%%%%%%%%%%%%%%%%%%%Fourth order in space%%%%%%%%%%%%%%%%%%%%%
%     U_bar(:,1)=Up(:,1)+dt*(Jp(:,1)-(Fp(:,2)-Fp(:,1))/dx)+Sp(:,1);
%     U_bar(:,2)=Up(:,2)+dt*(Jp(:,2)-(Fp(:,3)-Fp(:,2))/dx)+Sp(:,2);
%     for j=3:n-1
%         U_bar(:,j)=Up(:,j)+dt*(Jp(:,j)-(7*Fp(:,j)-8*Fp(:,j-1)+Fp(:,j-2))/(6*dx))+Sp(:,j);
%     end
%     U_bar(:,n)=Up(:,n)+dt*(Jp(:,n)-(Fp(:,n)-Fp(:,n-1))/dx)+Sp(:,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%predictor end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%calculating primitive variables%%%%%%%%%%%%%%%%%%%
     [rhop,up,Pp,Tp,ep]=Decode(U_bar,n); 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
 
    %%%%%%%%%%%%%%%setting boundary conditions for predictor steps%%%%%%
    up(1)=0; %predictor velocity
    up(n)=0;
    
    Tp(1)=2*Tp(2)-Tp(3); %predictor temperature
    Tp(n)=2*Tp(n-1)-Tp(n-2);
    
    ep(1)=cv*Tp(1);
    ep(n)=cv*Tp(n);
    
    Pp(1)=2*Pp(2)-Pp(3); %predictor pressure
    Pp(n)=2*Pp(n-1)-Pp(n-2);
    
    rhop(1)= Pp(1)/(R*Tp(1));%predictor density
    rhop(n)=Pp(n)/(R*Tp(n));
    
    [U_bar,F_bar,J_bar]=Code(rhop,up,Pp,Tp,n,ep,dx);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%Corrector step %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%2nd order%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    U_new(:,1)=0.5*(Up(:,1)+U_bar(:,1)+dt*(J_bar(:,1)-(F_bar(:,2)-F_bar(:,1))/dx));%+S_bar(:,1);
    for j=2:n-1
    U_new(:,j)=0.5*(Up(:,j)+U_bar(:,j)+dt*(J_bar(:,j)-(F_bar(:,j)-F_bar(:,j-1))/dx));%+S_bar(:,j);
    end
    U_new(:,n)=0.5*(Up(:,n)+U_bar(:,n)+dt*(J_bar(:,n)-(F_bar(:,n)-F_bar(:,n-1))/dx));%+S_bar(:,n);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%Fourth order%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     U_new(:,1)=0.5*(Up(:,1)+U_bar(:,1)+dt*(J_bar(:,1)-(F_bar(:,2)-F_bar(:,1))/dx))+S_bar(:,1);
%     for j=2:n-2
%     U_new(:,j)=0.5*(Up(:,j)+U_bar(:,j)+dt*(J_bar(:,j)-(-7*F_bar(:,j)+8*F_bar(:,j+1)-F_bar(:,j+2))/(6*dx)))+S_bar(:,j);
%     end
%     U_new(:,n-1)=0.5*(Up(:,n-1)+U_bar(:,n-1)+dt*(J_bar(:,n-1)-(F_bar(:,n-1)-F_bar(:,n-2))/dx))+S_bar(:,n-1);
%     U_new(:,n)=0.5*(Up(:,n)+U_bar(:,n)+dt*(J_bar(:,n)-(F_bar(:,n)-F_bar(:,n-1))/dx))+S_bar(:,n);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    %%%%%%%%%%%%%%%%%%%%%%%%%primitive variables  calculated%%%%%%%%%%%%%%
    [rho_new,u_new,P_new,T_new,e_new]=Decode(U_new,n); 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%setting boundary conditions for corrector steps%%%%%%%%
    u_new(1)=0; 
    u_new(n)=0;
    
    T_new(1)=2*T_new(2)-T_new(3); 
    T_new(n)=2*T_new(n-1)-T_new(n-2);
    
    e_new(1)=cv*T_new(1);
    e_new(n)=cv*T_new(n);
    
    P_new(1)=2*P_new(2)-P_new(3); 
    P_new(n)=2*P_new(n-1)-P_new(n-2);
    
    rho_new(1)= P_new(1)/(R*T_new(1));
    rho_new(n)=P_new(n)/(R*T_new(n));
   [Up,Fp,Jp]=Code(rho_new,u_new,P_new,T_new,n,e_new,dx);
   plot(x,rho_new,'o');
   grid on
   pause(0.0001)
   axis([0 10 0 1])
   if (dt*k>0.006)
       break;
   end
   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%Adaptive grid%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   mu=muf*(Tf+110)/((Tf)^1.5)*((T_new).^1.5)./(T_new+110);
%   visc=mu./rho_new;
%   visc_d=max((4/3)*visc,visc.*gamma_Pr);
%   sonic=sqrt(gamma*R*T_new);
%   dt_cfl=((abs(u_new)./dx)+(sonic./(dx)^2)+((2*visc_d)./(dx)^2)).^-1;
%   dt=R_cfl*min(dt_cfl); %dt optimized 
%   dt_n(k)=dt; %array of-dt at different time steps
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
   subplot(221)
   plot(x,rho_new,'o');
   hold on
   plot(xx,rhoa,'-');
   xlabel('x')
   ylabel('Density')
   
   subplot(222)
   plot(x,P_new,'o');
   hold on
   plot(xx,Pa,'-');
   xlabel('x')
   ylabel('Pressure')
   
   subplot(223)
   plot(x,T_new,'o');
   hold on
   plot(xx,Ta,'-');
   xlabel('x')
   ylabel('Temperature')
   
   subplot(224)
   plot(x,u_new,'o');
   hold on
   plot(xx,ua,'-');
   xlabel('x')
   ylabel('velocity')
   grid on


%%%%%%%%%%%%%%%%%%Functions for coding and decoding%%%%%%%%%%%%%%%%%%%%%%%
% coding values for U,F,J
function [U,F,J]=Code(rho,u,P,T,n,e,dx)
%given values
Pr=0.707;
cp=1.005e3;
cv=0.718e3;
Tf=298;
muf=1.846e-5;

Cx=0.1; %Cx- chosen as per suggested in Project

%Setting Array
U=zeros(3,n);
F=zeros(3,n);
J=zeros(3,n);
S=zeros(3,n);
du_dx=zeros(1,n);
dT_dx=zeros(1,n);
mu=zeros(1,n);
dmu_dx=zeros(1,n);
St_xx=zeros(1,n);
K=zeros(1,n);
d2u_dx2=zeros(1,n);
d2T_dx2=zeros(1,n);

%Constant calculated for Mu
A=muf*(Tf+110)/((Tf)^1.5);

%Boundary conditions setting for differentials-velocity,Temperature
du_dx(1)=(-3*u(1)+4*u(2)-u(3))/(2*dx); %du/dx at BC=1
du_dx(n)=(3*u(n)-4*u(n-1)+u(n-2))/(2*dx);%du/dx at BC=n

d2u_dx2(1)=(-u(4)+4*u(3)-5*u(2)+2*u(1))/(dx)^2;%d2u/dx2 at BC=1
d2u_dx2(n)=(2*u(n)-5*u(n-1)+4*u(n-2)-u(n-3))/(dx)^2;%d2u/dx2 at BC=n

dT_dx(1)=0; %dT/dx at BC=1
dT_dx(n)=0;%dT/dx at BC=n

d2T_dx2(1)=0; %d2T/dx2 at BC=1
d2T_dx2(n)=0;%d2T/dx2 at BC=n

%setting the values of differentials at inner nodes-Velcity,Temperature
for i=2:n-1
    du_dx(i)=(u(i+1)-u(i-1))/(2*dx);
    dT_dx(i)=(T(i+1)-T(i-1))/(2*dx);
    
    d2u_dx2(i)=(u(i+1)-2*u(i)+u(i-1))/(dx)^2;
    d2T_dx2(i)=(T(i+1)-2*T(i)+T(i-1))/(dx)^2;
end

for i=1:n
    
mu(i)=A*((T(i))^1.5)/(T(i)+110); %mu calculated
K(i)=mu(i)*(cp/Pr); %Thermal conductivity
dmu_dx(i)=A*((-(T(i)^1.5)/(T(i)+110)^2)+1.5*((T(i)^0.5)/(T(i)+110)))*dT_dx(i); %rate of change of mu
St_xx(i)=(du_dx(i)*mu(i))*(4/3); %viscous stress
dq_dx(i)=(K(i)*d2T_dx2(i)+dT_dx(i)*dmu_dx(i)*(cp/Pr)); %rate of transfer of heat

%Calculation of U,F,J using primitive variables
U(1,i)=rho(i); 
U(2,i)=rho(i).*u(i);
U(3,i)=rho(i).*e(i);

F(1,i)=U(2,i);
F(2,i)=rho(i).*(u(i).^2)+P(i);
F(3,i)=U(3,i).*u(i);

J(1,i)=0;
J(2,i)=0;%(4/3)*(d2u_dx2(i)*mu(i)+du_dx(i)*dmu_dx(i));
J(3,i)=-P(i)*du_dx(i);%+dq_dx(i);%+St_xx(i)*du_dx(i);
end
%calculation of artifical viscous term
% for i=2:n-1
% S(1,i)=Cx*(abs(P(i+1)-2*P(i)+P(i-1))/(P(i+1)+2*P(i)+P(i-1)))*(U(1,i+1)-2*U(1,i)+U(1,i-1));
% S(2,i)=Cx*(abs(P(i+1)-2*P(i)+P(i-1))/(P(i+1)+2*P(i)+P(i-1)))*(U(2,i+1)-2*U(2,i)+U(2,i-1));
% S(3,i)=Cx*(abs(P(i+1)-2*P(i)+P(i-1))/(P(i+1)+2*P(i)+P(i-1)))*(U(3,i+1)-2*U(3,i)+U(3,i-1));
% end
end

%%%%%%%%%%%%%%%%%%%%%%%Decoding values for u,P,p,e,T%%%%%%%%%%%%%%%%%%%%%%
function [rho,u,P,T,e]=Decode(U,n)
%given values
cp=1.005e3;
cv=0.718e3;
R=286.9;
%setting array
rho=zeros(1,n);
u=zeros(1,n);
e=zeros(1,n);
T=zeros(1,n);
P=zeros(1,n);
%primitive variables calculation
for i= 1:n
 rho(i)= U(1,i);
 u(i)=U(2,i)./U(1,i);
 e(i)=U(3,i)./U(1,i);
 T(i)=e(i)/cv;
 P(i)=rho(i)*R*T(i);
end
end