Skip to content

Instantly share code, notes, and snippets.

@beeflamian
Created December 27, 2013 04:12
Show Gist options
  • Save beeflamian/8142480 to your computer and use it in GitHub Desktop.
Save beeflamian/8142480 to your computer and use it in GitHub Desktop.
%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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment