close all
clear all
clc
%Inputs
L = 3;%length in m
n = 31;%Number of nodes
x = linspace(0,L,n);
g =1.4;%Gamma
R = 278;%Specific gas constant
dx1 = abs(x(1)-x(2));
[d,n1]=grid(L,g,R);%User-defined function keeps decreasing dx until error<1e-2
%returns the element size as d
dx = d;
n=n1;
x = linspace(0,L,n);
dt(L,n,x,dx,g,R,100);%User-defined function to determine mass flow rate
%returns the plot of mass flow rate at different time
%steps for conservative flow based on which asuitable
%dt value can be chosen, but for 100 time steps only
dtn(L,n,x,dx,g,R,100);%User-defined function to determine mass flow rate
%returns the plot of mass flow rate at different time
%steps for non-conservative flow based on which asuitable
%dt value can be chosen, but for 100 time steps only
cons(L,n,x,dx,g,R);%User-defined function to determine conservative flow
%returns the plot of Pressure, Temperature, density,
%mach-number, Velocity along the nozzle length
ncons(L,n,x,dx,g,R);%User-defined function to determine non-conservative flow
%returns the plot of Pressure, Temperature, density,
%mach-number, Velocity along the nozzle length
function [d,n]= grid(L,g,R)%grid convergence test
clear all
close all
clc
%inputs
n=3;
err = 1;
e = 1;
dx=1;
while dx>0.01
L = 3;
n = n;
x = linspace(0,L,n);
g =1.4;
R = 278;
dx = abs(x(1)-x(2));
r0 =1.225;
T0 = 288;
a0 = sqrt(T0);
r = 1 - 0.3146 *x;
T = 1 - 0.2314 *x;
V = (0.1+1.09*x).*(T.^0.5);
ai = sqrt(T);
p = r.*R.*T;
A = 1+2.2*(x-1.5).^2;
C = 0.5;
dt = min(C*dx./(ai+V));
rold =r;
told = T;
pold = p;
Vold = V;
Rinit = r;
tinit = T;
pinit = p;
Vinit = V;
r= Rinit;
T = tinit;
p =pinit;
V = Vinit;
rn=r;
rold =r;
told = T;
pold = p;
Vold = V;
for j = 2: n-1
t1(j) = -r(j)*(V(j+1)-V(j))/dx;
t2(j) = -r(j)*V(j)*(log(A(j+1))-log(A(j)))/dx;
t3(j) = -V(j)*(r(j+1)-r(j))/dx;
drdt_p(j) = t1(j)+t2(j)+t3(j);
ru(j) = r(j)+(drdt_p(j)*dt);
end
for k=2:n-1
t7(k) = -r(k)*(V(k)-V(k-1))/dx;
t8(k) = -r(k)*V(k)*(log(A(k))-log(A(k-1)))/dx;
t9(k) = -V(k)*(r(k)-r(k-1))/dx;
drdt_c(k) = t7(k)+t8(k)+t9(k);
drdt_a(k) = 0.5*(drdt_p(k)+drdt_c(k));
rnew(k) = r(k)+drdt_a(k)*dt;
end
r=rnew;
r(1) = 1;
r(n) = 2*r(n-1)-r(n-2);
rnn(e) = r(round(n/2));
nn(e) = n;
e=e+1;
ro=rn;
rn=r;
err = max(abs(rn-ro));
if err<1e-2
d=dx;
dx=1e-10;
nx= n;
end
n=n+1;
end
figure
plot(nn,rnn)
title('Number of grid v/s density');
xlabel('Number of grids')
ylabel('Density')
end
function [m] = dt(L,n,x,dx,g,R,t)%Conservative mass flow rate for different time steps
%inputs
L = 3;
n = n;
x = linspace(0,L,n);
g =1.4;
R = 278;
dx = dx;
r0 =1.225;
T0 = 288;
a0 = sqrt(T0);
r = 1 - 0.3146 *x;
T = 1 - 0.2314 *x;
V = (0.1+1.09*x).*(T.^0.5);
ai = sqrt(T);
p = r.*R.*T;
A = 1+2.2*(x-1.5).^2;
C = 0.0001;
dt= [0.001 0.0005 0.00025 0.000125 0.0000625];
ndt = max(size(dt));
tz = t;
rold =r;
told = T;
pold = p;
Vold = V;
Rinit = r;
tinit = T;
pinit = p;
Vinit = V;
for b=1:5
r= Rinit;
T = tinit;
p =pinit;
V = Vinit;
for i = 1:tz
rold =r;
told = T;
pold = p;
Vold = V;
for j = 2:n-1
drdt_p(j) = -((r(j+1)-r(j))*(log(A(j+1))-log(A(j)))*(V(j+1)-V(j)))/dx;
t1(j) = -(V(j)/r(j))*drdt_p(j);
t2(j) = 1/(g*r(j))*((r(j+1)-r(j))*(T(j+1)-T(j)));
dvdt_p(j) = t1(j)+t2(j);
t3(j) = (g-1)*((A(j+1)-A(j))*(V(j+1)-V(j)))/dx;
t4(j) = ((r(j+1)-r(j))*(T(j+1)-T(j))*(V(j+1)-V(j))*(A(j+1)-A(j)))/dx;
dtdt_p(j) = -drdt_p(j)+t3(j)+t4(j);
ru(j) = r(j)+drdt_p(j)*dt(b);
vu(j) = V(j)+dvdt_p(j)*dt(b);
tu(j) = T(j)+dtdt_p(j)*dt(b);
end
for k =2:n-1
drdt_c(k) = -((ru(k)-ru(k-1))*(log(A(k))-log(A(k-1)))*(vu(k)-vu(k-1)))/dx;
drdt_a(k) = 0.5*(drdt_p(k)+drdt_c(k));
t5(k) = -(vu(k)/ru(k))*drdt_a(k);
t6(k) = 1/(g*ru(k))*((ru(k)-ru(k-1))*(tu(k)-tu(k-1)));
dvdt_c(k) = t5(k)+t6(k);
dvdt_a(k) = 0.5*(dvdt_p(k)+dvdt_c(k));
t7(k) = (g-1)*((A(k)-A(k-1))*(vu(k)-vu(k-1)))/dx;
t4(k) = ((ru(k)-ru(k-1))*(tu(k)-tu(k-1))*(vu(k)-vu(k-1))*(A(k)-A(k-1)))/dx;
dtdt_c(k) = -drdt_c(k)+t3(k)+t4(k);
dtdt_a(k) = 0.5 * (dtdt_c(k)+dtdt_p(k));
rnew(k) = rold(k)+(drdt_a(k) * dt(b));
vnew(k) = Vold(k)+(dvdt_a(k) * dt(b));
tnew(k) = told(k)+(dtdt_a(k) * dt(b));
end
r = rnew;
r(1) = 1;
r(n) = 2*r(n-1) - r(n-2);
V = vnew;
V(1) = 2*V(1) - V(2);
V(n) = 2*V(n-1) - V(n-2);
T =tnew;
T(1) = 1;
T(n) = 2*T(n-1) - T(n-2);
end
p = r(1:30).*R.*T(1:30);
M = (V(1:30).*a0)./ai(1:30);
m = r.*V.*A; %Mass flow rate
o=find(x==1.5);
drdt_am = drdt_a;
dvdt_am = dvdt_a;
dtdt_am=dtdt_a;
pt = p;
if b==1
m1=m;
else if b==2
m2=m;
else if b==3
m3=m;
else if b==4
m4=m;
else if b==5
m5=m;
end
end
end
end
end
end
figure
plot(x,m1,x,m2,x,m3,x,m4,x,m5)
title('Mass flow rate at different time steps for conservative flow');
xlabel('Nozzle length')
ylabel('Mass flow rate')
legend('dt=0.001','dt=0.0005','dt=0.00025','dt=0.000125','dt=0.0000625');
end
function [m] = dtn(L,n,x,dx,g,R,t)%Non-conservative mass flow rate for different time steps
%inputs
L = 3;
n = n;
x = linspace(0,L,n);
g =1.4;
R = 278;
dx = dx;
r0 =1.225;
T0 = 288;
a0 = sqrt(T0);
r = 1 - 0.3146 *x;
T = 1 - 0.2314 *x;
V = (0.1+1.09*x).*(T.^0.5);
ai = sqrt(T);
p = r.*R.*T;
A = 1+2.2*(x-1.5).^2;
C = 0.0001;
dt= [0.001 0.0005 0.00025 0.000125 0.0000625];
ndt = max(size(dt));
tz = t;
rold =r;
told = T;
pold = p;
Vold = V;
Rinit = r;
tinit = T;
pinit = p;
Vinit = V;
for b=1:5
r= Rinit;
T = tinit;
p =pinit;
V = Vinit;
for i=1:tz
rold =r;
told = T;
pold = p;
Vold = V;
for j = 2: n-1
t1(j) = -r(j)*(V(j+1)-V(j))/dx;
t2(j) = -r(j)*V(j)*(log(A(j+1))-log(A(j)))/dx;
t3(j) = -V(j)*(r(j+1)-r(j))/dx;
drdt_p(j) = t1(j)+t2(j)+t3(j);
t4(j) =-V(j)*(V(j+1)-V(j))/dx;
t5(j) = -(1/g)*(T(j+1)-T(j))/dx;
t6(j) = -(T(j)/(g*r(j)))*(r(j+1)-r(j))/dx;
dvdt_p(j) = t4(j)+t5(j)+t6(j);
t7(j) = -V(j)*(T(j+1)-T(j))/dx;
t8(j) = -(g-1)*T(j)*(V(j+1)-V(j))/dx;
t9(j) = -(g-1)*T(j)*V(j)*(log(A(j+1))-log(A(j)))/dx;
dtdt_p(j) = t7(j)+t8(j)+t9(j);
ru(j) = r(j)+(drdt_p(j)*dt(b));
vu(j) = V(j)+dvdt_p(j)*dt(b);
tu(j) = T(j)+dtdt_p(j)*dt(b);
end
for k=2:n-1
t7(k) = -r(k)*(V(k)-V(k-1))/dx;
t8(k) = -r(k)*V(k)*(log(A(k))-log(A(k-1)))/dx;
t9(k) = -V(k)*(r(k)-r(k-1))/dx;
drdt_c(k) = t7(k)+t8(k)+t9(k);
drdt_a(k) = 0.5*(drdt_p(k)+drdt_c(k));
t10(k) =-V(k)*(V(k)-V(k-1))/dx;
t11(k) = -(1/g)*(T(k)-T(k-1))/dx;
t12(k) = -(T(k)/(g*r(k)))*(r(k)-r(k-1))/dx;
dvdt_c(k) = t10(k)+t11(k)+t12(k);
dvdt_a(k) = 0.5*(dvdt_p(k)+dvdt_c(k));
t13(k) = -V(k)*(T(k-1)-T(k))/dx;
t14(k) = -(g-1)*T(k)*(V(k)-V(k-1))/dx;
t15(k) = -(g-1)*T(k)*V(k)*(log(A(k))-log(A(k-1)))/dx;
dtdt_c(k) = t13(k)+t14(k)+t15(k);
dtdt_a(k) = 0.5*(dtdt_p(k)+dtdt_c(k));
rnew(k) = r(k)+drdt_a(k)*dt(b);
vnew(k) = V(k)+dvdt_a(k)*dt(b);
tnew(k) = T(k)+dtdt_a(k)*dt(b);
end
r=rnew;
r(1) = 1;
r(n) = 2*r(n-1)-r(n-2);
V = vnew;
V(1) = 2*V(1)-V(2);
V(n) = 2*V(n-1)-V(n-2);
T =tnew;
T(1) = 1;
T(n) = 2*T(n-1)-T(n-2);
end
P = r.*R.*T;
M = (V*a0)./ai;
m = r.*V.*A;%non-conservative mass flowrate
o=find(x==1.5);
drdt_am = drdt_a;
dvdt_am = dvdt_a;
dtdt_am=dtdt_a;
pt = p;
if b==1
m1=m;
else if b==2
m2=m;
else if b==3
m3=m;
else if b==4
m4=m;
else if b==5
m5=m;
end
end
end
end
end
end
figure
plot(x,m1,x,m2,x,m3,x,m4,x,m5)
title('Mass flow rate at different time steps for non-conservative flow');
xlabel('Nozzle length')
ylabel('Mass flow rate')
legend('dt=0.001','dt=0.0005','dt=0.00025','dt=0.000125','dt=0.0000625');
end
function [r,T,V,p,M] = cons(L,n,x,dx,g,R)%Conservative flow
r0 =1.225;
t0 = 288;
a0 = sqrt(g*r0*t0);
r = 1 - 0.3146 *x;
T = 1 - 0.2314 *x;
V = (0.1+1.09*x).*(T.^0.5);
ai = sqrt(T);
p = r.*R.*T;
A = 1+2.2*(x-1.5).^2;
C = 0.00125;
dt = min(C*dx./(ai+V));
rold =r;
told = T;
pold = p;
Vold = V;
Rinit = r;
tinit = T;
pinit= p;
Vinit= V;
for i = 1:1500
rold =r;
told = T;
pold = p;
Vold = V;
for j = 2:n-1
drdt_p(j) = -((r(j+1)-r(j))*(log(A(j+1))-log(A(j)))*(V(j+1)-V(j)))/dx;
t1(j) = -(V(j)/r(j))*drdt_p(j);
t2(j) = 1/(g*r(j))*((r(j+1)-r(j))*(T(j+1)-T(j)));
dvdt_p(j) = t1(j)+t2(j);
t3(j) = (g-1)*((A(j+1)-A(j))*(V(j+1)-V(j)))/dx;
t4(j) = ((r(j+1)-r(j))*(T(j+1)-T(j))*(V(j+1)-V(j))*(A(j+1)-A(j)))/dx;
dtdt_p(j) = -drdt_p(j)+t3(j)+t4(j);
ru(j) = r(j)+drdt_p(j)*dt;
vu(j) = V(j)+dvdt_p(j)*dt;
tu(j) = T(j)+dtdt_p(j)*dt;
end
for k =2:n-1
drdt_c(k) = -((ru(k)-ru(k-1))*(log(A(k))-log(A(k-1)))*(vu(k)-vu(k-1)))/dx;
drdt_a(k) = 0.5*(drdt_p(k)+drdt_c(k));
t5(k) = -(vu(k)/ru(k))*drdt_a(k);
t6(k) = 1/(g*ru(k))*((ru(k)-ru(k-1))*(tu(k)-tu(k-1)));
dvdt_c(k) = t5(k)+t6(k);
dvdt_a(k) = 0.5*(dvdt_p(k)+dvdt_c(k));
t7(k) = (g-1)*((A(k)-A(k-1))*(vu(k)-vu(k-1)))/dx;
t4(k) = ((ru(k)-ru(k-1))*(tu(k)-tu(k-1))*(vu(k)-vu(k-1))*(A(k)-A(k-1)))/dx;
dtdt_c(k) = -drdt_c(k)+t3(k)+t4(k);
dtdt_a(k) = 0.5 * (dtdt_c(k)+dtdt_p(k));
rnew(k) = rold(k)+(drdt_a(k) * dt);
vnew(k) = Vold(k)+(dvdt_a(k) * dt);
tnew(k) = told(k)+(dtdt_a(k) * dt);
end
r = rnew;
r(1) = 1;
r(n) = 2*r(n-1) - r(n-2);
V = vnew;
V(1) = 2*V(1) - V(2);
V(n) = 2*V(n-1) - V(n-2);
T =tnew;
T(1) = 1;
T(n) = 2*T(n-1) - T(n-2);
rmax(i) = max(r);
rt(i) = i;
end
p = r.*R.*T;
M = (V*a0)./ai;
m = r.*V.*A; %Mass flow rate
figure
subplot(18,1,[1,2,3])
plot(x,r)
title({' ';'';'Conservative density v/s nozzle length'})
ylabel('Density(r/r0)');
subplot(18,1,[4,5,6])
plot(x,V)
title({' ';'';'Conservative Velocity v/s nozzle length'});
ylabel('Velocity(V/a0)');
subplot(18,1,[7,8,9])
plot(x,T)
title({' ';'';'Conservative Temperature v/s nozzle length'});
ylabel('Temperature(T/T0)');
subplot(18,1,[10,11,12])
plot(x,p)
title({' ';'';'Conservative Pressure v/s nozzle length'});
ylabel('Pressure(P/P0)');
subplot(18,1,[13,14,15])
plot(x,M)
title({' ';'';'Conservative Mach number v/s nozzle length'});
ylabel('Mach number (V/a0)');
subplot(18,1,[16,17,18])
plot(x,m)
title({' ';'';'Conservative mass flowrate v/s nozzle length'});
xlabel('Nozzle length');
ylabel('Mass flowrate m');
end
function [r,T,V,p,M] = ncons(L,n,x,dx,g,R)%Non-conservative flow
r0 =1.225;
T0 = 288;
a0 = sqrt(T0);
r = 1 - 0.3146 *x;
T = 1 - 0.2314 *x;
V = (0.1+1.09*x).*(T.^0.5);
ai = sqrt(T);
p = r.*R.*T;
A = 1+2.2*(x-1.5).^2;
C = 0.0001;
dt = min(C*dx./(ai+V));
rold =r;
told = T;
pold = p;
Vold = V;
Rinit = r;
tinit = T;
pinit= p;
Vinit= V;
for i=1:1500
rold =r;
told = T;
pold = p;
Vold = V;
for j = 2: n-1
t1(j) = -r(j)*(V(j+1)-V(j))/dx;
t2(j) = -r(j)*V(j)*(log(A(j+1))-log(A(j)))/dx;
t3(j) = -V(j)*(r(j+1)-r(j))/dx;
drdt_p(j) = t1(j)+t2(j)+t3(j);
t4(j) =-V(j)*(V(j+1)-V(j))/dx;
t5(j) = -(1/g)*(T(j+1)-T(j))/dx;
t6(j) = -(T(j)/(g*r(j)))*(r(j+1)-r(j))/dx;
dvdt_p(j) = t4(j)+t5(j)+t6(j);
t7(j) = -V(j)*(T(j+1)-T(j))/dx;
t8(j) = -(g-1)*T(j)*(V(j+1)-V(j))/dx;
t9(j) = -(g-1)*T(j)*V(j)*(log(A(j+1))-log(A(j)))/dx;
dtdt_p(j) = t7(j)+t8(j)+t9(j);
ru(j) = r(j)+(drdt_p(j)*dt);
vu(j) = V(j)+dvdt_p(j)*dt;
tu(j) = T(j)+dtdt_p(j)*dt;
end
for k=2:n-1
t7(k) = -r(k)*(V(k)-V(k-1))/dx;
t8(k) = -r(k)*V(k)*(log(A(k))-log(A(k-1)))/dx;
t9(k) = -V(k)*(r(k)-r(k-1))/dx;
drdt_c(k) = t7(k)+t8(k)+t9(k);
drdt_a(k) = 0.5*(drdt_p(k)+drdt_c(k));
t10(k) =-V(k)*(V(k)-V(k-1))/dx;
t11(k) = -(1/g)*(T(k)-T(k-1))/dx;
t12(k) = -(T(k)/(g*r(k)))*(r(k)-r(k-1))/dx;
dvdt_c(k) = t10(k)+t11(k)+t12(k);
dvdt_a(k) = 0.5*(dvdt_p(k)+dvdt_c(k));
t13(k) = -V(k)*(T(k-1)-T(k))/dx;
t14(k) = -(g-1)*T(k)*(V(k)-V(k-1))/dx;
t15(k) = -(g-1)*T(k)*V(k)*(log(A(k))-log(A(k-1)))/dx;
dtdt_c(k) = t13(k)+t14(k)+t15(k);
dtdt_a(k) = 0.5*(dtdt_p(k)+dtdt_c(k));
rnew(k) = r(k)+drdt_a(k)*dt;
vnew(k) = V(k)+dvdt_a(k)*dt;
tnew(k) = T(k)+dtdt_a(k)*dt;
end
r=rnew;
r(1) = 1;
r(n) = 2*r(n-1)-r(n-2);
V = vnew;
V(1) = 2*V(1)-V(2);
V(n) = 2*V(n-1)-V(n-2);
T =tnew;
T(1) = 1;
T(n) = 2*T(n-1)-T(n-2);
end
P = r.*R.*T;
M = (V*a0)./ai;
m = r.*V.*A;%non-conservative mass flowrate
figure
subplot(18,1,[1,2,3])
plot(x,r)
title({'';'Non-Conservative density v/s nozzle length'})
ylabel('Density(r/r0)');
subplot(18,1,[4,5,6])
plot(x,V)
title({'';'Non-Conservative Velocity v/s nozzle length'});
ylabel('Velocity(V/a0)');
subplot(18,1,[7,8,9])
plot(x,T)
title({'';'Non-Conservative Temperature v/s nozzle length'});
ylabel('Temperature(T/T0)');
subplot(18,1,[10,11,12])
plot(x,P)
title({'';'Non-Conservative Pressure v/s nozzle length'});
ylabel('Pressure(P/P0)');
subplot(18,1,[13,14,15])
plot(x,M)
title({'';'Non-Conservative Mach number v/s nozzle length'});
xlabel('Nozzle length');
ylabel('Mach number (V/a0)');
subplot(18,1,[16,17,18])
plot(x,m)
title({'';'Non-Conservative mass flowrate v/s nozzle length'});
xlabel('Nozzle length');
ylabel('Mass flowrate m');
end