1D Shock-Tube (MacChromac Model)

Solve 1D shock equation with MacChromac Model

MATLAB Script: 

%% This code solves the shock tube problem

close all

clear all

%Properties @ T= 298 Kelvin

Cp = 1.0049*1000; % specific heat at constant pressure

Cv = 0.7178*1000; % specific heat at constant volume

gamma = 1.4;

prandtl = 0.707;

R = Cp-Cv;

mu0 = 1.846*10^-5;  %initial dynamic viscosity

%Given Initial Condition

P1 = 10^4; %Pa

P4 = 10^5; %Pa

Pr = 3.031; %pressure ratio, P2/P1

initial_rho1 = 0.125;    %kg/m^3

initial_rho4 = 1;        %kg/m^3 

Lx = 10;   %meters

Nx = 101;  %grid pts

Cx = 0.3;  %artificial damping coefficient 

dx = Lx/(Nx-1);

xarray = 0:dx:10;

x0 = 5;         %diaphram start location (change)

n0 = x0/dx+1;   %node of diaphram start

tn = 6.1/1000;  %total time (change) 

%grid gen

%initialize

P = zeros(Nx,1);

P(1:n0) = P4;

P(n0+1:end) = P1;   %initial pressure values

u = zeros(Nx,1);    %initial velocity

tau_xx = zeros(Nx,1);

qx = zeros(Nx,1);

rho = zeros(Nx,1);

rho(1:n0) = initial_rho4;

rho(n0+1:end) = initial_rho1;   %initial density values

T=P./(R.*rho); %intial temp. Kelvin

e = T*Cv;      %initial energy v

%U, E and S matrix size Nx by 3

U = [rho, rho.*u, rho.*(e+u.^2/2)];

U1 = U;     %one time step before

U2 = U1;    %two time steps before

E = [rho.*u, rho.*u.*u+P-tau_xx, (U(:,3)+P).*u-u.*tau_xx+qx];

E1 = E;

E2 = E1;

S = zeros(Nx,3);

dUdt_star = zeros(Nx,3);

dUdt_bar = zeros(Nx,3);

dTdx = zeros(Nx,1);

dudx = zeros(Nx,1);

%finding dt

mu = mu0*(T/298).^(3/2).*(298+110)./(T+110);

k = mu.*Cp./prandtl;

nu_eff = 4/3*mu./rho;

big_Gamma = k./(rho.*Cv);

nu_prime = max(nu_eff, big_Gamma);

dt_cfl = 1./(abs(u)./dx+sqrt(gamma*R.*T)./dx+2*nu_prime*(1./dx./dx));

dt = 0.7*min(dt_cfl); %time step (will change)

t=0;

count = 0;

count_t = 0;

time = [10 20 30 40 50 60]

for tt = 1:length(time)

    count_t = count_t+1;

for w = 1:1:time

    count = count +1;

% while t <=tn

    %MacCormack Solver

    %step 1 (predictor step)

    for i=2:Nx-1

        ratio = Cx.*abs(P(i+1)-2*P(i)+P(i-1))./(P(i+1)+2*P(i)+P(i-1));

        S(i,:) = [ratio ratio ratio].*(U(i+1,:)-2*U(i,:)+U(i-1,:));  %artificial damping

        dUdt_star(i,:) = -xdf2(E, i);    %forward diff (will error at end of matrix)

    end

    %inlet/outlet BC

    dUdt_star(1,:) = -xdf2(E, 1);

    dUdt_star(Nx,:) = dUdt_star(Nx-1,:);  %for now just say dUdx = 0

    %intermediate values

    U_bar = U + dUdt_star.*dt+S;

    rho = U_bar(:,1);

    u = U_bar(:,2)./rho;

    e = U_bar(:,3)./rho-u.*u/2;

    T = e./Cv;

    P = rho.*R.*T;

    %stress

    mu = mu0*(T/298).^(3/2).*(298+110)./(T+110);

    lambda = -2/3.*mu;

    k = mu.*Cp./prandtl;

    for i=2:Nx-1

        dTdx(i) = (dTdx(i+1)-dTdx(i-1))/(2*dx);

        dudx(i) = (dudx(i+1)-dudx(i-1))/(2*dx);

    end

    dudx(1) = dudx(2);

    dudx(Nx) = dudx(Nx-1);

    dTdx(1) = dTdx(2);

    dTdx(Nx) = dTdx(Nx-1);

    qx = -k.*dTdx;

    tau_xx = lambda.*dudx + 2*mu.*dudx;

    E_bar = [rho.*u, rho.*u.*u+P-tau_xx, (U_bar(:,3)+P).*u-u.*tau_xx+qx];

    %Step 2 (corrector step)

    for i=2:Nx-1

        ratio = Cx.*abs(P(i+1)-2*P(i)+P(i-1))./(P(i+1)+2*P(i)+P(i-1));

        S(i,:) = [ratio ratio ratio].*(U_bar(i+1,:)-2*U_bar(i,:)+U_bar(i-1,:));  %artificial damping

        dUdt_bar(i,:) = -xdb2(E_bar, i);    %backward diff (will error at end of matrix)

    end

    %inlet/outlet BC

    dUdt_bar(1,:) = dUdt_bar(2,:);

    dUdt_bar(Nx,:) = -xdb2(E_bar, Nx);

    %final values

    dUdt_avg = 0.5*(dUdt_star + dUdt_bar);

    U = U + dUdt_avg.*dt + S;

    rho = U(:,1);

    u(:,count) = U(:,2)./rho;

    e = U(:,3)./rho-u(:,count).*u(:,count)/2;

   %e=U(:,3)./rho;

    T = e./Cv;

    P(:,count) = rho.*R.*T;

    %find stresses

    mu = real(mu0*(T/298).^(3/2).*(298+110)./(T+110));

    lambda = real(-2/3.*mu);

    k = real(mu.*Cp./prandtl);

    for i=2:Nx-1

        dTdx(i) = (dTdx(i+1)-dTdx(i-1))/(2*dx);

        dudx(i) = (dudx(i+1)-dudx(i-1))/(2*dx);

    end

    dudx(1) = dudx(2);

    dudx(Nx) = dudx(Nx-1);

    dTdx(1) = dTdx(2);

    dTdx(Nx) = dTdx(Nx-1);

    qx = -k.*dTdx;

    tau_xx = lambda.*dudx + 2*mu.*dudx;

    E = [rho.*u(:,count), rho.*u(:,count).*u(:,count)+P(:,count)-tau_xx, (U(:,3)+P(:,count)).*u(:,count)-u(:,count).*tau_xx+qx];

    %find dt

    nu_eff = 4/3*mu./rho;

    big_Gamma = k./(rho.*Cv);

    nu_prime = max(nu_eff, big_Gamma);

    dt_cfl = (1./(abs(u(:,count))./dx+sqrt(gamma*R.*T)./dx+2*nu_prime*(1./dx./dx)));

    dt_cfl = real(dt_cfl);

    dt = 0.7*min(dt_cfl); %time step

    dt = real(dt);

    t = (t + dt);

    t = real(t);

    count = count +1;

    if t>=0.0061

        t=0.0061;

    end

    P_history(:,count_t) = P(:,end);

    t_history(:,count_t) = t;

    u_history(:,count_t) = u(:,end);

    rho_history(:,count_t) = rho;

    T_history(:,count_t) = T;

end

end

%% Comparison with analytical calculation

 load('tr.mat');

% Pressure

figure()

subplot(2,2,1)

plot(x,P_x,'b','linewidth',2.5)

hold on

plot(xarray,P,'--r','linewidth',2.5)

set(gca,'fontsize',14)

xlabel('X-Coordinate [m]')

ylabel('Pressure [Pa]')

legend('Analytic sol^n','Numeric sol^n')

title(sprintf('Pressure vs Postion'),'fontsize',14);

ylim([0 1.2e5])

grid on

% Velocity

subplot(2,2,2)

plot(x,U_x,'b','linewidth',2.5)

hold on

plot(xarray,u,'--r','linewidth',2.5)

set(gca,'fontsize',14)

xlabel('X-Coordinate [m]')

ylabel('Velocity [m/s]')

legend('Analytic sol^n','Numeric sol^n','location','northwest')

title(sprintf('Velocity vs Position'),'fontsize',14);

grid on

ylim([0 400])

% Density

subplot(2,2,3)

plot(x,rho_x,'b','linewidth',2.5)

hold on

plot(xarray,rho,'--r','linewidth',2.5)

set(gca,'fontsize',14)

xlabel('X-Coordinate [m]')

ylabel('Density [kg/m^3]')

legend('Analytic sol^n','Numeric sol^n')

title(sprintf('Density vs Position'),'fontsize',14);

ylim([0 1.2])

grid on

% Temperature

subplot(2,2,4)

plot(x,T_x,'b','linewidth',2.5)

hold on

plot(xarray,T,'--r','linewidth',2.5)

set(gca,'fontsize',14)

xlabel('X-Coordinate [m]')

ylabel('Temperature [k]')

legend('Analytic sol^n','Numeric sol^n','location','northwest')

title(sprintf('Temperature vs Position'),'fontsize',14);

ylim([200 500])

grid on