Created
December 27, 2013 04:12
-
-
Save beeflamian/8142480 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%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