Conservative and non conservative form
non conservative
clear all
close all
clc
% inputs
n = 31
x = linspace(0,3,n);
dx = x(2) - x(1);
gama = 1.4;
% calculate initial profiles
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
% area
a = 1+2.2*(x-1.5).^2;
n_throat=find(a-1<0.001)
% time steps
nt = 50;
dt = 0.01;
% outer the loop
for k = 1:nt
rho_old = rho;
v_old = v;
t_old = t;
% predictor method
for j = 2:n-1;
%dhro_dt_p = -rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j)))/dx-(rho(j+1)-rho(j))/dx
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/gama)*(dtdx+(t(j)/rho(j))*drhodx);
% energy equation
dtdt_p(j) = -v(j)*dtdx-(gama-1)*t(j)*(dvdx+v(j)*dlogadx);
% 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
% corrector method
for j = 2:n-1;
%dhro_dt_p = -rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j)))/dx-(rho(j+1)-rho(j))/dx
dvdx = (v(j+1)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (t(j)-t(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/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/gama)*(dtdx+(t(j)/rho(j))*drhodx);
% energy equation
dtdt_c(j) = -v(j)*dtdx-(gama-1)*t(j)*(dvdx+v(j)*dlogadx);
end
% compute the average time derivative
drhodt = 0.5*(drhodt_p+drhodt_c)
dvdt = 0.5*(dvdt_p+dvdt_c)
dtdt = 0.5*(dtdt_p+dtdt_c)
% final solution update
for i = 2:n-1
rho(i) = rho_old(i)+drhodt(i)*dt
v(i) = v_old(i)+dvdt(i)*dt
t(i) = t_old(i)+dtdt(i)*dt
% apply the boundary condition
% 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) = 2*t(n-1)- t(n-2);
% calculating
% mach number
m = v./sqrt(t);
% mass flow rate
massflow = rho.*v.*a;
p = rho.*t;
% rho,t,p,m throat values
rho_throat(k) = rho(n_throat)
t_throat(k) = t(n_throat)
p_throat(k) = p(n_throat)
m_throat(k) = m(n_throat)
% check for time
hold on;
if k = 1
plot(x,massflow,\'--k\');
if k = 50;
plot(x,massflow,\'b\');
end
end
end
end
grid on
figure(2)
plot(x,rho)
legend(\'density ratio\')
ylabel(\'\\rho/\\rho_{o}\');
xlabel(\'x/l\');
grid on;
subplot(4,1,2)
plot(x,v,\'b\')
legend(\'velocity\')
xlabel(\'x/l\');
ylabel(\'t/t_{0}\');
grid on
subplot(4,1,3)
plot(x,t,\'b\')
legend(\'temparature\')
ylabel(\'t/t_{w}\');
xlabel(\'x/l\');
grid on;
figure(3)
plot(x,v,\'r\')
xlabel(\'non dimensional distance\')
ylabel(\'non dimensional velocity\')
legend(\'non conservative\')
figure(4)
plot(x,rho,\'p\')
xlabel(\'non dimensional density\')
ylabel(\'non dimensional distance\')
legend(\'non concervative\')
-----------------------------------------
conservative
clear all
close all
clc
% inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
rho = ones(1,n);
t = ones(1,n);
gamma = 1.4;
error = 1e9;
tol = 1e-5;
% area
a = 1+2.2*((x-1.5).^2);
n_throat = find(a-1<0.001)
for i = 1:n
if x(1) <= 0.5
rho(i) = 1;
t(i) = 1;
x(i)>0.5 & x(i) <= 1.5
rho(i) = 1-0.366*(x(i)-0.5);
t(i) = 1-0.67*(x(i)-0.5);
if x(i)>1.5&&x(i)<=3.5
rho(i) = 0.634-0.3879*(x(i)-1.5);
t(i) = 0.833-0.3507*(x(i)-1.5);
end
end
% velocity
u2 = 0.59; %for initial condition only
v = u2./(rho.*a);
p = rho(i).*t;
% initial condition
u1 = rho.*a;
u2 = (rho.*a).*v;
u3 = rho.*(t./(gamma-1))+(gamma/2*v.^2).*a;
k2 = 0;
% time marching
nt = 1000;
dt = 0.0267;
massflow_old = rho.*v.*a;
for k = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
% defining flux value
f1 = u2;
f2 = ((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(((gamma/2)*u2.^2)./u1));
f3 = (gamma*u2.*u3)./u1-(gamma*(gamma-1)/2)*((u2.^3)./(u1.^2));
% predictor value
for i = 2:n-1
j2 = (1/gamma)*rho(i)*t(i)*((a(i+1)-a(i))/dx);
% continuity equation
du1_dt_p(i) = -(f1(i+1)-f1(i))/dx;
% momentum equation
du2_dt_p(i) = -(f2(i+1)-f2(i))/dx+j2;
% energy equation
du3_dt_p(i) = -(f3(i+1)-f3(i))/dx;
%solution update
u1(i) = u1(i)+du1_dt_p(i)*dt;
u2(i) = u2(i)+du2_dt_p(i)*dt;
u3(i) = u3(i)+du3_dt_p(i)*dt;
% values
v = u2./u1;
rho(i) = u1(i)/a(i);
t(i) = (gamma-1)*((u3(i)/u1(i))-((gamma/2)*((u2(i)./u1(i)).^2)));
p = rho(i);
end
f1 = u2;
f2 = ((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(((gamma/2)*u2.^2)./u1));
f3 = (gamma*u2.*u3)./u1-(gamma*(gamma-1)/2)*((u2.^3)./(u1.^2));
% corrector step
for j = 2:n-1
j2 = (1/gamma)*rho(j)*t(j)*((a(j+1)-a(j))/dx);
% continuity equation
du1_dt_c(j) = -(f1(j)-f1(j-1))/dx;
% momentum equation
du2_dt_c(j) = -(f2(j)-f2(j-1))/dx+j2;
% energy equation
du3_dt_c(j) = -(f3(j)-f3(j-1))/dx;
end
% average solution update
du1_dt = 0.5*(du1_dt_p+du1_dt_c);
du2_dt = 0.5*(du2_dt_p+du2_dt_c);
du3_dt = 0.5*(du3_dt_p+du3_dt_c);
% final solution
for r = 2:n-1
u1(r) = u1_old(r)+du1_dt(r)*dt;
u2(r) = u2_old(r)+du2_dt(r)*dt;
u3(r) = u3_old(r)+du3_dt(r)*dt;
end
rho = u1./a;
v = u2./u1;;
t = (gamma-1)*(u3./u1-(gamma/2)*(v.^2));
% boundary conditiom
% inlet condition
u1(1) = a(1);
u2(1) = 2*u2(2)-u2(3);
v(1) = u2(1)/u1(1);
t(1) = 1;
u3(1)= u1(1)*(t(1)/(gamma-1)-(gamma/2)*(v(1)^2));
% outlet condition
u1(n) = 2*u1(n-1) - u1(n-2);
u2(n) = 2*u2(n-1) - u2(n-2);
u3(n) = 2*u3(n-1) - u3(n-3);
rho(n) = u1(n)/a(n);
v(n) = u2(n)/u1(n);
t(n) = (gamma-1)*((u3(n)/u2(n))-(gamma-1)*(v(n)^2));
p = rho(i).*t;
% calculation of other variable
m = v./sqrt(t);
% mass flow rate
massflow = rho.*v.*a;
p= rho.*t;
t_throatk(2) = t(i).*(n_throat);
p_throatk(2) = p(i).*(n_throat);
m_throatk(2) = m(i).*(n_throat);
% check for the convergence
error2 = max(abs(massflow-massflow_old));
if k2 == 100
plot(x,massflow_con,\'g\');
if k2 == 200
plot(x,massflow_con,\'r\');
if k2 == 700
plot(x,massflow_con,\'-.b\',\'linewidth\',1.25);
end
end
end
end
end
plot(x,rho,\'r\')
title(\'transient non-dimensional mass flow rates (non-conservative form)\');
xlabel(\'x/l\');
ylabel(\'\\rhova/\\rho{o}v_{o}a_{o}\');
legend(\'100\\deltat\',\'200\\deltat\',\'700\\deltat\');
axis([0 3 0.5 0.8]);
gri d on;
subplot(4,1,2)
plot(x,v,\'b\')
legend(\'velocity\')
xlabel(\'xl\');
ylabel(\'v/v{o}\');
grid on;
subpiot(4,1,3)
plot(x,t,\'b\')
legend(\'temparature\')
xlabel(\'x/l\');
ylabel(\'t/t_{o}\');
grid on;
figure(2)
plot(x,v,\'r\')
xlabel(\'non dimensionless distance\')
ylabel(\'non dimensional velocity\')
legend(\'conservative\')
figure(3)
subpolt(4,1,1);
plot(x,rho)
plot(x,rho)
legend(\'density ratio\');
xlabel(\'x/l\');
ylabel(\'\\rho\\/rho_{o}\');
grid on;
figure(4)
plot(x,rho,\'g\')
xlabel(\'non dimensionless distance\')
ylabel(\'non dimensional density\')
legend(\'conservative\')
-------------------------------------------------
the stability of the nonconservative from solver has been provided below;-
the flow rate is also high lighted at the inlet and an outlet portion.also we have found that all possible values at thne mach number =1 at the throat.
the simulation was run for a maximum of 900 time-steps with n = 31.
the tolerance was kept to 1e-5 where the change in rho,v and t had to fall under the tolerasnce the errror is considered pf 1e-9 for conservative form
we can say that conservative form solves in less number of iterations that non conservative form.however the net mass flow rate is in unstable conditions