close all
clear all
clc
%initial conditions
n=31
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.01+1.09*x).*t.^0.5;
a=1+2.2*(x-1.5).^2;
%time steps
nt=5;
dt=0.001;
%outer time loop
for kt=1:nt
%continuity equation
for j=2:n-1
%predictor method
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%continuity equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%momentum equation
dvdt_p(j)=v(j)*dvdx-(1/gamma)*dtdx+(t(j)/rho(j))*drhodx;
%energy Equation
dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*dvdx;
%solution update
rho(j)=rho(j)+drhodt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
end
for j=2:n-1
%corrector method
dvdx=(v(j-1)-v(j))/dx;
drhodx=(rho(j-1)-rho(j))/dx;
dtdx=(t(j-1)-t(j))/dx;
dlogadx=(log(a(j-1))-log(a(j)))/dx;
%continuity equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%momentum equation
dvdt_c(j)=v(j)*dvdx-(1/gamma)*dtdx+(t(j)/rho(j))*drhodx;
%energy Equation
dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*dvdx;
end
%computing the value
dvdt=0.5*(dvdt_p+dvdt_c);
drhodt=0.5*(drhodt_p+drhodt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
%final Values
for i=2:n-1
v(i)= v(i)+dvdt(i)*dt;
rho(i)=rho(i)+drhodt(i)*dt;
t(i)=t(i)+dtdt(i)*dt;
end
%applying the boundary conditions
%inlet
v(1)=2*v(2)-v(3);
%outlet
v(n)=2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
t(n)=t(n-1)-t(n-2);
plot(x,v)
end