gistfile1.txt

%Solution of quasi-one-dimensional isentropic flow through a
%diverging-converging nozzle, by the explicit Mac Cormack finite difference method.

%Nozzle geometry
dx=0.1;
x=0:dx:3;
A=(1+2.2*(x.^2-3*x+2.25))';
cfl=0.5; %Courant number
gama=1.4; %Specific heat ratio
N=length(x);%Number of grid points
tt=1400; %Number of time steps

%Initial values
rho=(1-0.3146*x)';
T=(1-0.2314*x)';
for i=1:N
    V(i)=((0.1+x(i)*1.09)*(T(i)^0.5));
    deltat(i)=(cfl*(dx/((T(i)^0.5+V(i)))))';
end
dt=min(deltat); %Time step difference

%%Mac Cormack method%%

for j=1:tt
    
    %Prediction step%
    
    for k=2:N-1
    
rho_der(k)=(1/dx)*(-rho(k)*(V(k+1)-V(k))-rho(k)*V(k)*(log(A(k+1))...
    -log(A(k)))-V(k)*(rho(k+1)-rho(k)));

V_der(k)=(1/dx)*(-V(k)*(V(k+1)-V(k))-(1/gama)*((T(k+1)-T(k))+...
    (T(k)/rho(k))*(rho(k+1)-rho(k))));

T_der(k)=(1/dx)*(-V(k)*(T(k+1)-T(k))-(gama-1)*T(k)*((V(k+1)-V(k))...
    +(V(k)*(log(A(k+1))-log(A(k))))));

rho_bar(k)=rho(k)+rho_der(k)*dt;
V_bar(k)=V(k)+V_der(k)*dt;
T_bar(k)=T(k)+T_der(k)*dt;

    end

%Intermediate boundary conditions
V_bar(1)=2*V_bar(2)-V_bar(3);
rho_bar(1)=rho(1);
T_bar(1)=T(1);

    %Correction step%

    for m=2:N-1

rho_der_bar(m)=(1/dx)*(-rho_bar(m)*(V_bar(m)-V_bar(m-1))-rho_bar(m)*V_bar(m)...
    *(log(A(m))-log(A(m-1)))-V_bar(m)*(rho_bar(m)-rho_bar(m-1)));

V_der_bar(m)=(1/dx)*(-V_bar(m)*(V_bar(m)-V_bar(m-1))-(1/gama)*...
    ((T_bar(m)-T_bar(m-1))+(T_bar(m)/rho_bar(m))*(rho_bar(m)-rho_bar(m-1))));

T_der_bar(m)=(1/dx)*(-V_bar(m)*(T_bar(m)-T_bar(m-1))-(gama-1)*T_bar(m)*((V_bar(m)-V_bar(m-1))...
    +(V_bar(m)*(log(A(m))-log(A(m-1))))));

rho(m)=rho(m)+dt*0.5*(rho_der(m)+rho_der_bar(m));
V(m)=V(m)+dt*0.5*(V_der(m)+V_der_bar(m));
T(m)=T(m)+dt*0.5*(T_der(m)+T_der_bar(m));

    end

%Boundary conditions%
V(1)=2*V(2)-V(3);
V(31)=2*V(30)-V(29);
rho(31)=2*rho(30)-rho(29);
T(31)=2*T(30)-T(29);

Q(:,1)=rho;
Q(:,2)=V;
Q(:,3)=T;

j
Q

end