Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Before a simulation is run, it is important to ascertain the governing differential equation for the phenomena. Since it has to be run for one dimensional orientation, the G.D.E(Governing Differential Equation) will be transformed into one dimensional form.
So as it is known that there are three form of conservation equations that govern fluid flow:-
1)Continuity Equation
`del/(delt)(ubrace(intintint)_("Volume")rho.dnu)+ubrace(intint)_("Surface")rhoV.dS=0`
2)Momentum Equation(for inviscid flow)
`del/(delt)(ubrace(intintint)_("Volume")(rhou).dnu)+ubrace(intint)_("Surface")(rhouV).dS=-ubrace(intint)_("Surface")(p dS)_x`
3)Energy Equation
`del/(delt)(ubrace(intintint)_("Volume")rho(e+(V^2/2)).dnu)+ubrace(intint)_("Surface")rho(e+V^2/2)V.dS=-ubrace(intint)_("Surface")(pV).dS`
Where `nu`=volume
`rho`=density
V=Velocity vector
u=x-component of velocity
S=surface
p=pressure
e=Internal energy
All the above equations are simplified and first derived in a form of governing differential equation which is called Conservation form
1)Continuity Equation
`(del(rho.A))/(delt)+(del(rho.A.V))/(delx)=0`
2)Momentum Equation(for inviscid flow)
`(del(rho.V.A))/(delt)+(del(rho.V^2.A))/(delx)=-A(delp)/(delx)`
3)Energy Equation
`(del(rho.(e+V^2/2).A))/(delt)+(del(rho.(e+V^2/2).V.A))/(delx)=-A(del(p.A.V))/(delx)`
These are equations in their purest form.
But the terms which are of physical significance are p,rho,T and V which are called primitive variables.
So these conservative form of equation are further simplified to obtain differential terms with primitive variables. The G.D.E obtained in this form is called non-conservative form of differential equations.
1)Continuity Equation
`(del(rho.A))/(delt)+rho.A(delV)/(delx)+rho.V(delA)/(delx)+V.A(delrho)/(delx)=0`
2)Momentum Equation(for inviscid flow)
`rho(delV)/(delt)+rhoV(delV)/(delx)=-(delp)/(delx)`
By further simplifying `p=rho.R.T` and eliminating pressure term`rho(delV)/(delt)+rhoV(delV)/(delx)=-R(rho(delT)/(delx)+T(delrho)/(delx))`
3)Energy Equation
`rho(dele)/(delt)+rhoV(dele)/(delx)=-p(delV)/(delx)-pV(del(ln A))/(delx)`
By eliminating p like in momentum equation and replacing e as `e=c_v.T`
From further simplification is
`rhoc_v(delT)/(delt)+rhoVc_v(delT)/(delx)=-rhoRT((delV)/(delx)+V(del(ln A))/(delx))`
Now the derivations for all the above derived equations can be located in John Anderson Textbook.
Now there are a subtle differences when numerical solution is obtained from non-conservation form of equation and when it is obtained from conservation form of the equation.
Now first it will be solved using non-conservation form of the equation.
Non-conservational form
Now the above G.D.E are in terms of dimensional terms. Though from here it can be proceeded onwards with the simulation but it is important they are converted into non-dimensional terms.It is done as follows:-
`T'=T/T_0`
`rho'=rho/rho_0`...where `T_0`and `rho_0`are temperature and density values respectively of reservoir i.e first point in space.
`x'=x/L`where L=length of nozzle
let `a_0`be the speed of sound in the reservoir such that `a_0=sqrt(gammaRT_0)`
`V'=V/a_0`
`t'=t/(L/a_0)`
`A'=A/A^(ST)`
`A^(ST)`=sonic throat area
These dimensionless terms are basically ratio as defined above. This is done so that no matter what the physical specifications are, after finding the solution , the values of these ratios will give the values of original or dimensional variables, once physical specifications like reservoir density or Temperature,length of the nozzle are defined. But at any point if these one of these specifications are changed, the solution will be still obtained without running the program again by just changing the substitution in non-dimensional variable.
But first these non-conservation form must be converted from dimensional to non-dimensional terms.
So the set of equations that are obtained are:
1)from continuity equation
`(delrho')/(delt')=-rho'(delV')/(delx')-rho'V'(del(ln A'))/(delx')-V'(delrho')/(delx')`
2)from momentum equation
`(delV')/(delt')=-V'(delV')/(delx')-(1/gamma)((delT')/(delx')+(T')/(rho')(delrho)/(delx'))`
3)from energy equation
`(delT')/(delt')=-V'(delT')/(delx')-(gamma-1)T'((delV')/(delx')+V'(delln A')/(delx'))`
Now since Differential equations have been simplified to represent the primitive variables in derivatives of time, numerical simulation will commence.
The method that will be applied is explicit method. Maccormak method is a type of explicit method which will be used to simulate the nozzle flow.
Basic structure of Maccormak method is:-
In a time loop:-
1)First the value of variable is stored in another variable for eg. pressure is stored in "pressure initial"
2)Then time derivative of a variable such as density or pressure is calculated by forward differencing of all the space derivative terms in their formula.This is called predictor step.
3)Then the variables are updated using the time derivatives calculated in predictor step.
4)Then these updated variables are used to calculate the time derivatives of the variables but this time the space derivatives are backward differenced. This is called corrector step.
5)An average of time derivatives in predictor step and corrector step, is calculated.
6)Now the updated values of the variables which was done in between predictor and corrector step is neglected. And the updation is done on the variable ,which were stored at the starting point of time loop, using the average time derivative.
7)This final updation is the solution of the time loop and is repeated for the finite number of time loops untill a steady state reaches.
Steady state time:-
It is important to find the time at which the solution reaches steady state. But it has converged or not is just a matter of perspective. So the criteria of convergence may be that two consecutive solutions have an error less than predetermined tolerance and if that happens then it can be concluded that convergence has reached. But it just so happens that after few hundred time steps the increment in the values of solutions of two consecutive time step is so miniscule that the solution seems to have become stagnant but over a course of few hundred time steps it was found that change was there,though negligible when measured over consecutive time step, but significant over a hundred time step. So that time has to be found out when change in the values measured between solutions 10 time steps apart, the change is still less than the tolerance. And it is then when it can be assumed that solution has slowed down to the point that it can be concluded that the solution has converged.
First there are two inputs that have to be determined to find the time of convergence
1)K stability:- This value determines how far apart the solutions must be in terms of time steps for error to be calculated.
2)Tolerance:- Its the value which determines how accurate the solution must be when it stabilizes
Methodology to be incorporated in the code:-
Lets just say K_stability=100 and tolerance is 1e-4. Lets just say at time step or K=1700,
error=abs(Solution(1700)-Solution(1700-K_stability))
error=abs(Solution(1700)-Solution(1600))
And this error falls below the tolerance, then at K-k_s (1600) the system had reached convergence which was checked only when system reached at k=1700.
So time at k=1600 must be noted down as time of convergence.
Time step size:-
Adaptive time steping strategy is applied to get at the steady state faster than keeping time step constant which also ensures that the solution does not blow up due to large size of dt.
Now it is known that
Couriant number=`C=((a+V)*dt)/dx` Where C ranges from 0 to 1 and a is speed of sound and V is velocity of fluid and dt and dx are increments in time and space respectively.
But since non-dimensional terms are being utilized, this is transformed to represent non-dimensional terms.
given as
`C=(sqrt(T')*V')*(dt')/(dx')`
And further calculating the biggest time step that can be utilised without causing the solution to blow up:
`dt=(C.dx')/(sqrt(T')+V')`
Now all the strategy for time step size and simulation is converted in the form of the code:-
Note:-The initial condition has been designed to be closer to actual solution so that steady state is achieved faster. This initial condition has been given in the John D Anderson Textbook.
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
p=rho.*T;
%area
A=1+2.2*(x-1.5).^2;
%Time loop
nt=5000;
t=0;%to find the net time
%to find the time at which system stabilizes
tol=1e-4;
time=0;%to find time st which system stabilizes
tn=0;%to find time step number st which system stabilizes
C=0.5;%couriant Number
%outer time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
t=t+dt;
%storing old solution vector values
rho_old=rho;
v_old=v;
T_old=T;
p_old=p;
%predictor step
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTdx=(T(j+1)-T(j))/dx;
d_logA_dx=(log(A(j+1))-log(A(j)))/dx;
%Continuity Equation
drhodt_pred(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-v(j)*drhodx;
%momentum equation
dvdt_pred(j)=-v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
%energy equation
dTdt_pred(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*d_logA_dx);
%solution update
rho(j)=rho(j)+drhodt_pred(j)*dt;
v(j)=v(j)+dvdt_pred(j)*dt;
T(j)=T(j)+dTdt_pred(j)*dt;
endfor
%Corrextor method
for j=2:n-1
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dTdx=(T(j)-T(j-1))/dx;
d_logA_dx=(log(A(j))-log(A(j-1)))/dx;
%Continuity Equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-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+v(j)*d_logA_dx);
endfor
% compute average derivatives
drhodt=0.5*(drhodt_c+drhodt_pred);
dvdt=0.5*(dvdt_c+dvdt_pred);
dTdt=0.5*(dTdt_c+dTdt_pred);
%Final Solution
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;
endfor
%apply 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);
%pressure calculation
p=rho.*T;
%Steady state coding to find time at which solution stabilizes
a(k,1:n)=v;
b(k,1:n)=rho;
c(k,1:n)=T;
d(k,1:n)=p;
k_s=10;%time step difference to find error
if k>k_s
error_v=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_rho=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_T=max(abs(c(k,1:n)-c(k-k_s,1:n)));
error_p=max(abs(d(k,1:n)-d(k-k_s,1:n)));
if time==0
if tol>=error_v
if tol>=error_rho
if tol>=error_T
if tol>=error_p
time=t;
tn=k;
end
end
end
end
end
end
%steady state coding end
endfor
%Calculations for graphic display start
ymax=sqrt(max(A)/pi);
ny=1001;
nx=1001;
v_cont=-0.1*ones(ny,nx);
T_cont=-0.1*ones(ny,nx);
p_cont=-0.1*ones(ny,nx);
rho_cont=-0.1*ones(ny,nx);
y=linspace(-ymax,ymax,ny);
dy=y(2)-y(1);
x_new=linspace(0,3,nx);
for i=1:nx
i_n=floor((n-1)*(i-1)/(nx-1))+1;
if nx>i
v_x(i)=v(i_n)+((v(i_n+1)-v(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
A_x(i)=A(i_n)+((A(i_n+1)-A(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
else
v_x(i)=v(i_n);
A_x(i)=A(i_n);
end
end
for i=1:nx
b=sqrt(4*A_x(i)/pi);
ny_start=round(((-b/2)+ymax)/dy)+1;
ny_end=length(y)-ny_start+1;
for j=ny_start:ny_end
v_cont(j,i)=v_x(i);
endfor
endfor
%Calculations for graphic display end
%plot results
figure(1)
contourf(x_new,y,v_cont)
colorbar
colormap(hot)
xlabel('Length of the Convergence Divergence Nozzle')
ylabel('Breadth')
text_time=sprintf('Steady state time=%d s',time);
title({'Velocity profile contour in the approximate shape of nozzle',text_time})
figure(2)
plot(x,v,'color','r');
hold on
plot(x,T,'color','b');
hold on
plot(x,rho,'color','g');
hold on
plot(x,p,'color','m');
hold off
xlabel('x')
ylabel('output ratios')
legend('Velocity','Temperature','Density','Pressure')
title_text1=sprintf('Time=%d',t);
title_text2=sprintf('Steady state time=%d and No. of time steps=%d',time,tn);
title({title_text1,title_text2});
axis([0 3 -0.5 2.5]);In this code the end result has been displayed in the form of linear graph as well as a contour plot for velocity profile taking the shape of the nozzle itself.The two figures are as follows:-
figure 1
figure 2
As it can been seen steady state time has been found out such that error between values of primitive variables spaced between 10 time steps spaced between the tolerance of 1e-4.
Conservative form
During derivation of governing equation resulted into a form of equations which were termed as conservative form of equations. These equations consisted of derivative terms which were a combination of primitive variables in their prime form.They are as follows:-
1)Continuity Equation
`(del(rho.A))/(delt)+(del(rho.A.V))/(delx)=0`
2)Momentum Equation(for inviscid flow)
`(del(rho.V.A))/(delt)+(del(rho.V^2.A))/(delx)=-A(delp)/(delx)`
3)Energy Equation
`(del(rho.(e+V^2/2).A))/(delt)+(del(rho.(e+V^2/2).V.A))/(delx)=-A(del(p.A.V))/(delx)`
But they are in dimensional form and has to be converted into non-dimensional form.
1)Continuity Equation
`(del(rho'.A'))/(delt')+(del(rho'.A'.V'))/(delx')=0`
2)Momentum Equation(for inviscid flow)
`(del(rho'.V'.A'))/(delt')+(del(rho'.A'.(V')^2+(1/gamma)p'.A'))/(delx')=(1/gamma).p'.(delA')/(delx')`
3)Energy Equation
`(del(rho'.((e')/(gamma-1)+(gamma/2).(V')^2).A'))/(delt')+(del(rho'.((e')/(gamma-1)+(gamma/2).(V')^2).V'.A'+p'A'V'))/(delx')=0`
Now definitions for all non-dimensional form have been explained but e' was not explained.
since `e=C_p.T`
and `e_0=C_p.T_0`
Therefore
`e'=(C_p.T)/(C_p.T_0)=T/T_0=T'`
Now if it can be observed all these equations are of the form:-
`(delU)/(delt')+(delF)/(delx')=J`
where U=Solutions vector
F=the flux vector
J=source term
So by conforming this nomenclature to the non-dimensional form of governing differential equation in conservation form
1)Continuity Equation
`U_1=rho'.A'`
`F_1=rho'.A'.V'`
`J1=0`
The equation
`(delU_1)/(delt')=-(delF_1)/(delx')`
2)Momentum Equation
`U_2=rho'.A'.V'`
`F_2=rho'.A'.(V')^2+(1/gamma)p'.A'`
`J_2=(1/gamma).p'.(delA')/(delx')`
The equation
`(delU_2)/(delt')=-(delF_2)/(delx')+J2`
3)Energy Equation
`U_3=rho'.((e')/(gamma-1)+(gamma/2).(V')^2).A'`
`F_3=rho'.((e')/(gamma-1)+(gamma/2).(V')^2).V'.A'+p'A'V'`
`J_3=0`
The equation
`(delU_3)/(delt')=-(delF_3)/(delx')`
Now it has been found whenever values of flux vector is obtained from values of primitive variables which were inturn calculated from solution vector,the solution blew up.Therefore it becomes important to directly derive flux vector variable and source values from solution vector terms.Though getting values from primitive variable or from solution vector directly is theoretically the same but when incorporated in numerical simulation methodology, values calculated from primitive variables causes error to creep in.
So flux vector values have to obtained directly from flux vector.
`F_1=U_2`
`F_2=(U_2)^2/U_1+((gamma-1)/gamma)*(U_3-(gamma/2).(U_2)^2/U_1)`
`F_3=gamma.(U_2.U_3)/(U_1)-((gamma(gamma-1))/2).(U_2)^3/(U_1)^2`
`J_2=((gamma-1)/gamma)*(U_3-(gamma/2).(U_2)^2/U_1)(del(ln A'))/(delx')`
These formulas are incorporated in the coding for numerical simulation.
Here as well adaptive time steping is applied to reach faster to the solution.
The initial condition are again designed to achieve steady state faster and has been obtained from John D Anderson Textbook.
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
x1=0.5;
x2=1.5;
n1=((x1-x(1))/dx)+1;
n2=((x2-x(1))/dx)+1;
for i=1:n1
rho(i)=1;
T(i)=1;
end
for i=n1:n2
rho(i)=1-0.366*(x(i)-0.5);
T(i)=1-0.167*(x(i)-0.5);
end
for i=n2:n
rho(i)=0.634-0.3879*(x(i)-1.5);
T(i)=0.833-0.3507*(x(i)-1.5);
end
%area
A=1+2.2*(x-1.5).^2;
%velocity input
v=0.59./(rho.*A);
v_init=v;
%Time loop
nt=5000;
C=0.5;
t=0;
%to find the time at which system stabilizes
tol=1e-3;
time=0;
%defining the solution vector
U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*A.*((T./(gamma-1))+(0.5*gamma.*v.^2));
% Time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
t=t+dt;
%storing old solution vector values
U1_old=U1;
U2_old=U2;
U3_old=U3;
%defining flux vector
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));
%defining source vector for predictor step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(U3(j)-(gamma/2)*((U2(j)^2)/U1(j)))*((log(A(j+1))-log(A(j)))/dx);
endfor
%predictor step
for j=2:n-1
dU1dt_pred(j)=-1*(F1(j+1)-F1(j))/dx;
dU2dt_pred(j)=(-1*(F2(j+1)-F2(j))/dx)+J2(j);
dU3dt_pred(j)=-1*(F3(j+1)-F3(j))/dx;
%Solution Update
U1(j)=U1(j)+dU1dt_pred(j)*dt;
U2(j)=U2(j)+dU2dt_pred(j)*dt;
U3(j)=U3(j)+dU3dt_pred(j)*dt;
endfor
%Updation of flux vector and source terms
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));
%defining source vector for corrector step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(U3(j)-(gamma/2)*((U2(j)^2)/U1(j)))*((log(A(j))-log(A(j-1)))/dx);
endfor
%Corrector step
for j=2:n-1
dU1dt_cor(j)=-1*(F1(j)-F1(j-1))/dx;
dU2dt_cor(j)=(-1*(F2(j)-F2(j-1))/dx)+J2(j);
dU3dt_cor(j)=-1*(F3(j)-F3(j-1))/dx;
endfor
%average derivatives
dU1dt=0.5*(dU1dt_pred+dU1dt_cor);
dU2dt=0.5*(dU2dt_pred+dU2dt_cor);
dU3dt=0.5*(dU3dt_pred+dU3dt_cor);
%Final solution
for i=2:n-1
U1(i)=U1_old(i)+dU1dt(i)*dt;
U2(i)=U2_old(i)+dU2dt(i)*dt;
U3(i)=U3_old(i)+dU3dt(i)*dt;
endfor
%Primitive variable
rho=U1./A;
v=U2./U1;
T=(gamma-1)*((U3./U1)-(0.5*gamma.*(v.^2)));
p=rho.*T;
%Boundary Condition
%inlet
%U1(1) is fixed
%U2(1)=rho(1)*A(1)*v(1);
U2(1)=2*U2(2)-U2(3);
U3(1)= rho(1)*A(1)*((T(1)/(gamma-1))+(0.5*gamma*(v(1)^2)));
%Outlet
U1(n)=2*U1(n-1)-U1(n-2);
%U1(n)=rho(n)*A(n);
U2(n)=2*U2(n-1)-U2(n-2);
%U2(n)=rho(n)*A(n)*v(n);
U3(n)=2*U3(n-1)-U3(n-2);
%U3(n)= rho(n)*A(n)*((T(n)/(gamma-1))+(0.5*gamma*(v(n)^2)));
%to find time at which solution stabilizes
%Steady state coding to find time at which solution stabilizes
a(k,1:n)=U1;
b(k,1:n)=U2;
c(k,1:n)=U3;
k_s=10;%time step difference to find error
if k>k_s
error_U1=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_U2=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_U3=max(abs(c(k,1:n)-c(k-k_s,1:n)));
if time==0
if tol>=error_U1
if tol>=error_U2
if tol>=error_U3
time=t;
tn=k;
end
end
end
end
end
%steady state coding end
endfor
%Primitive variable
rho=U1./A;
v=U2./U1;
T=(gamma-1)*((U3./U1)-(0.5*gamma.*(v.^2)));
p=rho.*T;
%Calculations for graphic display start
ymax=sqrt(max(A)/pi);
ny=1001;
nx=1001;
v_cont=-0.1*ones(ny,nx);
T_cont=-0.1*ones(ny,nx);
p_cont=-0.1*ones(ny,nx);
rho_cont=-0.1*ones(ny,nx);
y=linspace(-ymax,ymax,ny);
dy=y(2)-y(1);
x_new=linspace(0,3,nx);
for i=1:nx
i_n=floor((n-1)*(i-1)/(nx-1))+1;
if nx>i
v_x(i)=v(i_n)+((v(i_n+1)-v(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
A_x(i)=A(i_n)+((A(i_n+1)-A(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
else
v_x(i)=v(i_n);
A_x(i)=A(i_n);
end
end
for i=1:nx
b=sqrt(4*A_x(i)/pi);
ny_start=round(((-b/2)+ymax)/dy)+1;
ny_end=length(y)-ny_start+1;
for j=ny_start:ny_end
v_cont(j,i)=v_x(i);
endfor
endfor
%Calculations for graphic display end
%plot results
figure(1)
contourf(x_new,y,v_cont)
colorbar
colormap(hot)
xlabel('Length of the Convergence Divergence Nozzle')
ylabel('Breadth')
text_time=sprintf('Steady state time=%d s',time);
title({'Velocity profile contour in the approximate shape of nozzle',text_time})
figure(2)
plot(x,v,'color','r');
hold on
plot(x,T,'color','b');
hold on
plot(x,rho,'color','g');
hold on
plot(x,p,'color','m');
hold off
xlabel('x')
ylabel('output ratios')
legend('Velocity','Temperature','Density','Pressure')
title_text1=sprintf('Time=%d',t);
title_text2=sprintf('Steady state time=%d and No. of time steps=%d',time,tn);
title({title_text1,title_text2});
axis([0 3 -0.5 2.5]);Now the plot is
And the graph
Now it is quite evident that the simulation can be done using both the form of governing equation:Non-Conservative and Conservative. But it is important to perform quantitative analysis on these two forms to establish which method is better.
To perform such a analysis, it is important that a function be set up for both methodologies so that it becomes convenient to employ any methodology with ease.
Function for Non-conservative form:-
function [v,rho,T,tim,time,tn]=non_conservative_form_solver(v_init,rho_init,T_init,A,x,nt,C,tol,k_stability)
%inputs
n=length(x);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho=rho_init;
T=T_init;
v=v_init;
p=rho.*T;
%Time loop
tim=0;%to find the net time
%to find the time at which system stabilizes
time=0;%to find time st which system stabilizes
tn=0;%to find time step number st which system stabilizes
%outer time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
tim=tim+dt;
t(k)=tim;
%storing old solution vector values
rho_old=rho;
v_old=v;
T_old=T;
p_old=p;
%predictor step
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTdx=(T(j+1)-T(j))/dx;
d_logA_dx=(log(A(j+1))-log(A(j)))/dx;
%Continuity Equation
drhodt_pred(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-v(j)*drhodx;
%momentum equation
dvdt_pred(j)=-v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
%energy equation
dTdt_pred(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*d_logA_dx);
%solution update
rho(j)=rho(j)+drhodt_pred(j)*dt;
v(j)=v(j)+dvdt_pred(j)*dt;
T(j)=T(j)+dTdt_pred(j)*dt;
endfor
%Corrextor method
for j=2:n-1
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dTdx=(T(j)-T(j-1))/dx;
d_logA_dx=(log(A(j))-log(A(j-1)))/dx;
%Continuity Equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-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+v(j)*d_logA_dx);
endfor
% compute average derivatives
drhodt=0.5*(drhodt_c+drhodt_pred);
dvdt=0.5*(dvdt_c+dvdt_pred);
dTdt=0.5*(dTdt_c+dTdt_pred);
%Final Solution
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;
endfor
%apply 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);
%pressure calculation
p=rho.*T;
%Steady state coding to find time at which solution stabilizes
a(k,1:n)=v;
b(k,1:n)=rho;
c(k,1:n)=T;
d(k,1:n)=p;
k_s=k_stability;%time step difference to find error
if k>k_s
error_v=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_rho=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_T=max(abs(c(k,1:n)-c(k-k_s,1:n)));
error_p=max(abs(d(k,1:n)-d(k-k_s,1:n)));
if time==0
if tol>=error_v
if tol>=error_rho
if tol>=error_T
if tol>=error_p
time=t(k-k_s);
tn=k-k_s;
end
end
end
end
end
end
%steady state coding end
endfor
endfunctionFunction for conservative form
function [U1,U2,U3,tim,time,tn]=conservative_form_solver(v_init,rho_init,T_init,A,x,nt,C,tol,k_stability)
%inputs
n=length(x);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho=rho_init;
T=T_init;
v=v_init;
%Time loop
tim=0;
%to find the time at which system stabilizes
time=0;
tn=0;
%defining the solution vector
U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*A.*((T./(gamma-1))+(0.5*gamma.*v.^2));
% Time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
tim=tim+dt;
t(k)=tim;
%storing old solution vector values
U1_old=U1;
U2_old=U2;
U3_old=U3;
%defining flux vector
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));
%defining source vector for predictor step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(1/A(j))*(U3(j)-((gamma/2)*((U2(j)^2)/U1(j))))*((A(j+1)-A(j))/dx);
endfor
%predictor step
for j=2:n-1
dU1dt_pred(j)=-1*(F1(j+1)-F1(j))/dx;
dU2dt_pred(j)=(-1*(F2(j+1)-F2(j))/dx)+J2(j);
dU3dt_pred(j)=-1*(F3(j+1)-F3(j))/dx;
%Solution Update
U1(j)=U1(j)+dU1dt_pred(j)*dt;
U2(j)=U2(j)+dU2dt_pred(j)*dt;
U3(j)=U3(j)+dU3dt_pred(j)*dt;
endfor
%Updation of flux vector and source terms
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));
%defining source vector for corrector step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(1/A(j))*(U3(j)-((gamma/2)*((U2(j)^2)/U1(j))))*((A(j)-A(j-1))/dx);
endfor
%Corrector step
for j=2:n-1
dU1dt_cor(j)=-1*(F1(j)-F1(j-1))/dx;
dU2dt_cor(j)=(-1*(F2(j)-F2(j-1))/dx)+J2(j);
dU3dt_cor(j)=-1*(F3(j)-F3(j-1))/dx;
endfor
%average derivatives
dU1dt=0.5*(dU1dt_pred+dU1dt_cor);
dU2dt=0.5*(dU2dt_pred+dU2dt_cor);
dU3dt=0.5*(dU3dt_pred+dU3dt_cor);
%Final solution
for i=2:n-1
U1(i)=U1_old(i)+dU1dt(i)*dt;
U2(i)=U2_old(i)+dU2dt(i)*dt;
U3(i)=U3_old(i)+dU3dt(i)*dt;
endfor
%Boundary Condition
%inlet
%U1(1) is fixed
%U2(1)=rho(1)*A(1)*v(1);
U2(1)=2*U2(2)-U2(3);
v(1)=U2(1)/U1(1);
U3(1)= rho(1)*A(1)*((T(1)/(gamma-1))+(0.5*gamma*(v(1)^2)));
%Outlet
U1(n)=2*U1(n-1)-U1(n-2);
%U1(n)=rho(n)*A(n);
U2(n)=2*U2(n-1)-U2(n-2);
%U2(n)=rho(n)*A(n)*v(n);
U3(n)=2*U3(n-1)-U3(n-2);
%U3(n)= rho(n)*A(n)*((T(n)/(gamma-1))+(0.5*gamma*(v(n)^2)));
%to find time at which solution stabilizes
%Steady state coding to find time at which solution stabilizes
%Primitive variable
rho=U1./A;
v=U2./U1;
T=(gamma-1)*((U3./U1)-(0.5*gamma.*(v.^2)));
p=rho.*T;
a(k,1:n)=v;
b(k,1:n)=rho;
c(k,1:n)=T;
d(k,1:n)=p;
k_s=k_stability;%time step difference to find error
if k>k_s
error_v=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_rho=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_T=max(abs(c(k,1:n)-c(k-k_s,1:n)));
error_p=max(abs(d(k,1:n)-d(k-k_s,1:n)));
if time==0
if tol>=error_v
if tol>=error_rho
if tol>=error_T
if tol>=error_p
time=t(k-k_s);
tn=k-k_s;
end
end
end
end
end
end
%steady state coding end
endfor
endfunctionNow a function to get required inputs to plot a contour plot to display output(i.e velocity,pressure,temperature,pressure or mass flow rate) contour in the approximate shape of the nozzle.
function [x_new,y,output_cont]=Graphic_contour_plot_function(A,output)
n=length(A);
x=linspace(0,3,n);
%Calculations for graphic display start
ymax=sqrt(max(A)/pi);
ny=1001;
nx=1001;
output_cont=-0.1*ones(ny,nx);
y=linspace(-ymax,ymax,ny);
dy=y(2)-y(1);
x_new=linspace(0,3,nx);
for i=1:nx
i_n=floor((n-1)*(i-1)/(nx-1))+1;
if nx>i
output_x(i)=output(i_n)+((output(i_n+1)-output(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
A_x(i)=A(i_n)+((A(i_n+1)-A(i_n))*((x_new(i)-x(i_n))/(x(i_n+1)-x(i_n))));
else
output_x(i)=output(i_n);
A_x(i)=A(i_n);
end
end
for i=1:nx
b=sqrt(4*A_x(i)/pi);
ny_start=round(((-b/2)+ymax)/dy)+1;
ny_end=length(y)-ny_start+1;
for j=ny_start:ny_end
output_cont(j,i)=output_x(i);
endfor
endfor
endfunctionQuantitative Analysis:-
1)Grid Independence test
This test is done to test whether solution values change if the number of grid points are varied. It is to test whether the current coding is sensitive to size of increment in space.
This will be done in two manner:-
Values will be compared at x=1.5 for velocity,Temperature and density at n=31,61 and 101.
Note:-the values will be checked for nt=5000 time steps.
Test for Non-Conservative Form:-
The coding:
clear all
close all
clc
%Input for the analysis
x_val=1.5;
number_grid_points=[31 61 101];
%Time loop
nt=5000;
C=0.5;
%to find the time at which system stabilizes
tol=1e-4;
k_stability=10;
for k=1:length(number_grid_points)
n=number_grid_points(k);
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
%calculate initial profile
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;
[v_out(k,1:n),rho_out(k,1:n),T_out(k,1:n),t(k),time(k),tn(k)]=non_conservative_form_solver(v,rho,T,A,x,nt,C,tol,k_stability);
num(k)=round(((x_val-x(1))/dx)+1);
Velocity_comparison(k)=v_out(k,num(k));
Density_comparison(k)=rho_out(k,num(k));
Temperature_comparison(k)=T_out(k,num(k));
Machnumber_comparison(k)=Velocity_comparison(k)/sqrt(Temperature_comparison(k));
endfor
number_grid_points
Velocity_comparison
Density_comparison
Temperature_comparison
Machnumber_comparisonAnd the output in command window:-
As it can be seen as number of grid points increase, the value of primitive variables at the throat remains the same value more or less. The accuracy increases as numbr of grid points increase. So failing of grid test would have implied large changes in the values of the variables. But the values remain approximately the same. So even at n=31, the methodology is independent of grid discretization. And the accuracy of these values can be affirmed by checking the values of Mach number at throat values. Since Mach number analytically is 1 at the throat and numerically it can be seen that the mach number is very close to 1 at the throat.This confirms the accuracy of the numerical solution to analytical values.
Test for Conservative test:-
clear all
close all
clc
%Input for the analysis
x_val=1.5;
number_grid_points=[31 61 101];
%Time loop
nt=5000;
C=0.5;
%to find the time at which system stabilizes
tol=1e-4;
k_stability=10;
for k=1:length(number_grid_points)
n=number_grid_points(k);
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
x1=0.5;
x2=1.5;
n1=round(((x1-x(1))/dx)+1);
n2=round(((x2-x(1))/dx)+1);
for i=1:n1
rho(i)=1;
T(i)=1;
end
for i=n1:n2
rho(i)=1-0.366*(x(i)-0.5);
T(i)=1-0.167*(x(i)-0.5);
end
for i=n2:n
rho(i)=0.634-0.3879*(x(i)-1.5);
T(i)=0.833-0.3507*(x(i)-1.5);
end
%area
A=1+2.2*(x-1.5).^2;
%velocity input
v=0.59./(rho.*A);
[U1,U2,U3,t(k),time(k),tn(k)]=conservative_form_solver(v,rho,T,A,x,nt,C,tol,k_stability);
%Primitive variable
rho_out(k,1:n)=U1./A;
v_out(k,1:n)=U2./U1;
T_out(k,1:n)=(gamma-1)*((U3./U1)-(0.5*gamma.*(v.^2)));
p_out(k,1:n)=rho_out(k,1:n).*T_out(k,1:n);
num(k)=round(((x_val-x(1))/dx)+1);
Velocity_comparison(k)=v_out(k,num(k));
Density_comparison(k)=rho_out(k,num(k));
Temperature_comparison(k)=T_out(k,num(k));
Pressure_comparison(k)=p_out(k,num(k));
Machnumber_comparison(k)=Velocity_comparison(k)/sqrt(Temperature_comparison(k));
endfor
number_grid_points
Velocity_comparison
Density_comparison
Temperature_comparison
Pressure_comparison
Machnumber_comparisonAnd the input in the command window
Now it can be observed that at n=31, the values of the variables become lax in terms of precision. Compared to values for n=61 and n=101, values at n=31 are considerably away. Even the mach number whose analytical value is known at the throat, it can be seen that at n=31 the mach number deviates to a value of 0.99145 when analytical value is 1 at the throat. Moreover Mach number at n=61 and n=101, are 1.00471 and 1.00522 respectively which is much more accurate than mach number at n=31. Now since at n-31, the solution seems to be a little away in terms of precision, but is still close when comes to failing grid dependency test. So n=31 will not give bad results but it is recommended to use n=61 or above in future.
Note:-
In the above test Mach number has been used. But it is confused with velocity ratio profile due to similar formula.
Now velocity ratio formula:
`v'=v/a_0=v/sqrt(gamma.R.T_0)` ....where T_0 is temperature of the reservoir.
Mach number formula:
`M=v/a=v/sqrt(gamma.R.T)` .....where T is the temperature at the point in space at which velocity is measured.
But since v',T' and rho' is calculated in the numerical simulation. Therefore Mach number should also be as a result of formula of these ratios if possible.
Now since
`M=v/sqrt(gamma.R.T)=(v/sqrt(gamma.R.T_0)).(sqrt(gamma.R.T_0)/sqrt(gamma.R.T))`
`M=(v')/sqrt(T')`
And this is the formula applied in the above coding to get value of Mach number.
2) Accuracy and Stability Test
Now since its difficult to get an analytical solution using the analytical equation but it is possible to check the solutions developed numerically and analytically at sonic throat of the nozzle. Since it is known that at throat Mach number is equal to 1 analytically, therefore values of other varibles can be found out analyticall using the formula:-
`p/p_0=(1+((gamma-1)/2)*M^2)^(-gamma/(gamma-1))`
`rho/rho_0=(1+((gamma-1)/2)*M^2)^(-1/(gamma-1))`
`T/T_0=(1+((gamma-1)/2)*M^2)^(-1)`
These formulas will be converted into a function
Function for analytical pressure
function [p]=analytical_p(M)
gamma=1.4;
term_a=1+((gamma-1)/2)*(M)^2;
p=term_a^(-gamma/(gamma-1));
endfunctionFunction for analytical density
function [rho]=analytical_rho(M)
gamma=1.4;
term_a=1+((gamma-1)/2)*(M)^2;
rho=term_a^(-1/(gamma-1));
endfunctionFunction for analytical Temperature
function [T]=analytical_T(M)
gamma=1.4;
term_a=1+((gamma-1)/2)*(M)^2;
T=term_a^(-1);
endfunctionNow a code will be written which will check the value of variables at the throat against time and will be plotted for both conservation as well as non-conservation form. All the initial values,k_stability,tolerance,x,nt and C will be same for both the forms.At the same time comparison will be done the analytical value which the system should reach in steady state.
This will tell not only the accuracy of the solution but also the time it takes to reach convergence for both the form of codes.
First a function is written which will take initial values,k_stability,tolerance,x,nt and C as inputs and Solution values at x=1.5 for each time step and time and time step at which solution converges as output.
function for non_conservation form
function [M_t,T_t,rho_t,p_t,tim,time,tn]=Stability_Non_Conservation_function(n,nt,v_init,rho_init,T_init,A,x_val,C,tol,k_stability)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho=rho_init;
T=T_init;
v=v_init;
p=rho.*T;
%Time loop
tim=0;
t=0;%to find the net time
time=0;%to find time st which system stabilizes
tn=0;%to find time step number st which system stabilizes
%outer time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
tim=tim+dt;
t(k)=tim;
%storing old solution vector values
rho_old=rho;
v_old=v;
T_old=T;
p_old=p;
%predictor step
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTdx=(T(j+1)-T(j))/dx;
d_logA_dx=(log(A(j+1))-log(A(j)))/dx;
%Continuity Equation
drhodt_pred(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-v(j)*drhodx;
%momentum equation
dvdt_pred(j)=-v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
%energy equation
dTdt_pred(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*d_logA_dx);
%solution update
rho(j)=rho(j)+drhodt_pred(j)*dt;
v(j)=v(j)+dvdt_pred(j)*dt;
T(j)=T(j)+dTdt_pred(j)*dt;
endfor
%Corrextor method
for j=2:n-1
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dTdx=(T(j)-T(j-1))/dx;
d_logA_dx=(log(A(j))-log(A(j-1)))/dx;
%Continuity Equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*d_logA_dx-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+v(j)*d_logA_dx);
endfor
% compute average derivatives
drhodt=0.5*(drhodt_c+drhodt_pred);
dvdt=0.5*(dvdt_c+dvdt_pred);
dTdt=0.5*(dTdt_c+dTdt_pred);
%Final Solution
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;
endfor
%apply 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);
%pressure and Mach number calculation
p=rho.*T;
M=v./sqrt(T);
%Steady state coding to find time at which solution stabilizes
a(k,1:n)=v;
b(k,1:n)=rho;
c(k,1:n)=T;
d(k,1:n)=p;
k_s=k_stability;%time step difference to find error
if k>k_s
error_v=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_rho=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_T=max(abs(c(k,1:n)-c(k-k_s,1:n)));
error_p=max(abs(d(k,1:n)-d(k-k_s,1:n)));
if time==0
if tol>=error_v
if tol>=error_rho
if tol>=error_T
if tol>=error_p
time=t(k-k_s);
tn=k-k_s;
end
end
end
end
end
end
%over time values analysis
%over time values analysis
num=round(((x_val-x(1))/dx)+1);
T_t(k)=T(num);
rho_t(k)=rho(num);
p_t(k)=p(num);
M_t(k)=M(num);
%end
endfor
endfunctionAnd the function for conservation form
function [M_t,T_t,rho_t,p_t,tim,time,tn]=Stability_Conservation_function(n,nt,v_init,rho_init,T_init,A,x_val,C,tol,k_stability)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho=rho_init;
T=T_init;
v=v_init;
p=rho.*T;
%Time loop
tim=0;
t=0;%to find the net time
time=0;%to find time st which system stabilizes
tn=0;%to find time step number st which system stabilizes
%defining the solution vector
U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*A.*((T./(gamma-1))+(0.5*gamma.*(v.^2)));
% Time loop
for k=1:nt
%dt change as stated in anderson textbook
delta_t=C.*(dx./((T.^0.5)+v));
dt=min(delta_t);
%calculating time
tim=tim+dt;
t(k)=tim;
U1_old=U1;
U2_old=U2;
U3_old=U3;
%defining flux vector
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)));
%defining source vector for predictor step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(1/A(j))*(U3(j)-((gamma/2)*((U2(j)^2)/U1(j))))*((A(j+1)-A(j))/dx);
endfor
%predictor step
for j=2:n-1
dU1dt_pred(j)=-1*(F1(j+1)-F1(j))/dx;
dU2dt_pred(j)=(-1*(F2(j+1)-F2(j))/dx)+J2(j);
dU3dt_pred(j)=-1*(F3(j+1)-F3(j))/dx;
%Solution Update
U1(j)=U1(j)+dU1dt_pred(j)*dt;
U2(j)=U2(j)+dU2dt_pred(j)*dt;
U3(j)=U3(j)+dU3dt_pred(j)*dt;
endfor
%Updation of flux vector and source terms
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));
%defining source vector for corrector step
for j=2:n-1
J2(j)=((gamma-1)/gamma)*(1/A(j))*(U3(j)-((gamma/2)*((U2(j)^2)/U1(j))))*((A(j)-A(j-1))/dx);
endfor
%Corrector step
for j=2:n-1
dU1dt_cor(j)=-1*(F1(j)-F1(j-1))/dx;
dU2dt_cor(j)=(-1*(F2(j)-F2(j-1))/dx)+J2(j);
dU3dt_cor(j)=-1*(F3(j)-F3(j-1))/dx;
endfor
%average derivatives
for j=2:n-1
dU1dt(j)=0.5*(dU1dt_pred(j)+dU1dt_cor(j));
dU2dt(j)=0.5*(dU2dt_pred(j)+dU2dt_cor(j));
dU3dt(j)=0.5*(dU3dt_pred(j)+dU3dt_cor(j));
endfor
%Final solution
for i=2:n-1
U1(i)=U1_old(i)+dU1dt(i)*dt;
U2(i)=U2_old(i)+dU2dt(i)*dt;
U3(i)=U3_old(i)+dU3dt(i)*dt;
endfor
%Boundary Condition
%inlet
%U1(1) is fixed
U2(1)=2*U2(2)-U2(3);
v(1)=U2(1)/U1(1);
U3(1)= rho(1)*A(1)*((T(1)/(gamma-1))+(0.5*gamma*(v(1)^2)));
%Outlet
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-2);
%to find time at which solution stabilizes
%Primitive variable
rho=U1./A;
v=U2./U1;
T=(gamma-1)*((U3./U1)-(0.5*gamma.*(v.^2)));
p=rho.*T;
M=v./sqrt(T);
%Steady state coding to find time at which solution stabilizes
a(k,1:n)=v;
b(k,1:n)=rho;
c(k,1:n)=T;
d(k,1:n)=p;
k_s=k_stability;%time step difference to find error
if k>k_s
error_v=max(abs(a(k,1:n)-a(k-k_s,1:n)));
error_rho=max(abs(b(k,1:n)-b(k-k_s,1:n)));
error_T=max(abs(c(k,1:n)-c(k-k_s,1:n)));
error_p=max(abs(d(k,1:n)-d(k-k_s,1:n)));
if time==0
if tol>=error_v
if tol>=error_rho
if tol>=error_T
if tol>=error_p
time=t(k-k_s);
tn=k-k_s;
end
end
end
end
end
end
%over time values analysis
num=((x_val-x(1))/dx)+1;
T_t(k)=T(num);
rho_t(k)=rho(num);
p_t(k)=p(num);
M_t(k)=M(num);
%end
endfor
endfunctionAnd then a main code is written which will have same inputs for both conservation and non-conservation form and the outputs of both are analysed. The mach number,temperature,pressure,density and mass flow rate is analysed.
clear all
close all
clc
%inputs
n=61;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate initial profile
rho_init=1-0.3146*x;
T_init=1-0.2314*x;
v_init=(0.1+1.09*x).*T_init.^0.5;
%area
A=1+2.2*(x-1.5).^2;
%Time loop
nt=5000;
%to find the time at which system stabilizes
x_val=1.5;
tol=1e-4;
k_stability=500;
C=0.5;%couriant Number
[M_t_non_con,T_t_non_con,rho_t_non_con,p_t_non_con,tim_non_con,time_non_con,tn_non_con]=Stability_Non_Conservation_function(n,nt,v_init,rho_init,T_init,A,x_val,C,tol,k_stability);
[M_t_con,T_t_con,rho_t_con,p_t_con,tim_con,time_con,tn_con]=Stability_Conservation_function(n,nt,v_init,rho_init,T_init,A,x_val,C,tol,k_stability);
t_num=1:1:nt;
M_anal=1;
T_anal=analytical_T(M_anal);
rho_anal=analytical_rho(M_anal);
p_anal=analytical_p(M_anal);
%plot for Mach number
figure(1)
plot(t_num,M_t_non_con,'r')
hold on
plot(t_num,M_t_con,'b')
hold on
plot([t_num(1) t_num(end)],[M_anal M_anal],'--','color','m')
hold on
plot([tn_non_con tn_non_con],[0 M_anal],'--','color','r');
hold on
plot([tn_con tn_con],[0 M_anal],'--','color','b');
xlabel('time steps')
ylabel('Mach number')
title_text_non_con=sprintf('time at which non conservation form stabilises = %d s',time_non_con);
title_text_con=sprintf('time at which conservation form stabilises = %d s',time_con);
title({'Mach number vs time steps',title_text_non_con,title_text_con})
text_legend_non_con=sprintf('time step of convergence(non con)=%d',tn_non_con);
text_legend_con=sprintf('time step of convergence(con)=%d',tn_con);
legend('Non conservation form','Conservation form','Analytical Value',text_legend_non_con,text_legend_con)
inc=0.01
axis([t_num(1) t_num(end) M_anal-inc M_anal+inc])
%plot for Temperature
figure(2)
plot(t_num,T_t_non_con,'r')
hold on
plot(t_num,T_t_con,'b')
hold on
plot([t_num(1) t_num(end)],[T_anal T_anal],'--','color','m')
hold on
plot([tn_non_con tn_non_con],[0 T_anal],'--','color','r');
hold on
plot([tn_con tn_con],[0 T_anal],'--','color','b');
xlabel('time steps')
ylabel('Temperature')
title_text_non_con=sprintf('time at which non conservation form stabilises = %d s',time_non_con);
title_text_con=sprintf('time at which conservation form stabilises = %d s',time_con);
title({'Temperature vs time steps',title_text_non_con,title_text_con})
text_legend_non_con=sprintf('time step of convergence(non con)=%d',tn_non_con);
text_legend_con=sprintf('time step of convergence(con)=%d',tn_con);
legend('Non conservation form','Conservation form','Analytical Value',text_legend_non_con,text_legend_con)
axis([t_num(1) t_num(end) T_anal-inc T_anal+inc])
%plot for Density
figure(3)
plot(t_num,rho_t_non_con,'r')
hold on
plot(t_num,rho_t_con,'b')
hold on
plot([t_num(1) t_num(end)],[rho_anal rho_anal],'--','color','m')
hold on
plot([tn_non_con tn_non_con],[0 rho_anal],'--','color','r');
hold on
plot([tn_con tn_con],[0 rho_anal],'--','color','b');
xlabel('time steps')
ylabel('Density')
title_text_non_con=sprintf('time at which non conservation form stabilises = %d s',time_non_con);
title_text_con=sprintf('time at which conservation form stabilises = %d s',time_con);
title({'Density vs time steps',title_text_non_con,title_text_con})
text_legend_non_con=sprintf('time step of convergence(non con)=%d',tn_non_con);
text_legend_con=sprintf('time step of convergence(con)=%d',tn_con);
legend('Non conservation form','Conservation form','Analytical Value',text_legend_non_con,text_legend_con)
axis([t_num(1) t_num(end) rho_anal-inc rho_anal+inc])
%plot for Pressure
figure(4)
plot(t_num,p_t_non_con,'r')
hold on
plot(t_num,p_t_con,'b')
hold on
plot([t_num(1) t_num(end)],[p_anal p_anal],'--','color','m')
hold on
plot([tn_non_con tn_non_con],[0 p_anal],'--','color','r');
hold on
plot([tn_con tn_con],[0 p_anal],'--','color','b');
xlabel('time steps')
ylabel('Pressure')
title_text_non_con=sprintf('time at which non conservation form stabilises = %d s',time_non_con);
title_text_con=sprintf('time at which conservation form stabilises = %d s',time_con);
title({'Pressure vs time steps',title_text_non_con,title_text_con})
text_legend_non_con=sprintf('time step of convergence(non con)=%d',tn_non_con);
text_legend_con=sprintf('time step of convergence(con)=%d',tn_con);
legend('Non conservation form','Conservation form','Analytical Value',text_legend_non_con,text_legend_con)
axis([t_num(1) t_num(end) p_anal-inc p_anal+inc])
%Mass flow rate
%since U2=rho.A.v=rho.A.(M*sqrt(T))
%A=1 since at throat A'=1
U2_t_non_con=rho_t_non_con.*(M_t_non_con.*sqrt(T_t_non_con));
U2_t_con=rho_t_con.*(M_t_con.*sqrt(T_t_con));
U2_anal=rho_anal.*(M_anal*sqrt(T_anal));
%plot for mass flow rate
figure(5)
plot(t_num,U2_t_non_con,'r')
hold on
plot(t_num,U2_t_con,'b')
hold on
plot([t_num(1) t_num(end)],[U2_anal U2_anal],'--','color','m')
hold on
plot([tn_non_con tn_non_con],[0 U2_anal],'--','color','r');
hold on
plot([tn_con tn_con],[0 U2_anal],'--','color','b');
xlabel('time steps')
ylabel('Mass flow rate')
title_text_non_con=sprintf('time at which non conservation form stabilises = %d s',time_non_con);
title_text_con=sprintf('time at which conservation form stabilises = %d s',time_con);
title({'Mass flow rate vs time steps',title_text_non_con,title_text_con})
text_legend_non_con=sprintf('time step of convergence(non con)=%d',tn_non_con);
text_legend_con=sprintf('time step of convergence(con)=%d',tn_con);
legend('Non conservation form','Conservation form','Analytical Value',text_legend_non_con,text_legend_con)
axis([t_num(1) t_num(end) U2_anal-inc U2_anal+inc])
%Stable values
sprintf('Analytical values Non conservation form Conservation form')
M_values=[M_anal M_t_non_con(end) M_t_con(end)]
T_values=[T_anal T_t_non_con(end) T_t_con(end)]
rho_values=[rho_anal rho_t_non_con(end) rho_t_con(end)]
p_values=[p_anal p_t_non_con(end) p_t_con(end)]
U2_values=[U2_anal U2_t_non_con(end) U2_t_con(end)]And the respective plots:
And the respective table of convergence values
First it is important to comment about the accuracy of the two forms. If the values of both the forms are compared at x=1.5 for n=61, both are fairly accurate. But it seems that non-conservation form is more accurate than conservation form as is represented by the table above and the curves as well. In 4 out of 5 curves the curve for non-conservation form is closer to the true value curve(analytical value curve) than conservation form curve. This means that Non-conservation form is more accurate than conservation form.
Now The speed is important to be discussed
The time step at which convergence happens for non-conservation form is 1605 and for conservation form is 1611.
That form is faster which acheives convergence faster. Now since non-conservation form achieves convergence faster, therefore non-conservation form is faster. But approximately time for convergence is same for both non-convergence form and conservation form. There is only difference for few time steps. It is important to understand that the transient state must be approximately the same for the two forms of G.D.E and if there is gross error occuring then there must be error in the coding. But the curves for bothe formes are approximately the same so the coding is right.
Which method to choose:-
If it is left upto me,then I would choose conservation form of equation. The reason is that the with conservation form of equations shock capturing can also be done but with non-conservation form shock capturing cannot be done. Though simulation using non-conservation form has proven to be faster and more accurate than conservation form of equations but the difference in accuracy and speed is small and hence insignificant. The simulation using conservation form covers the spectrum from subsonic-supersonic to shock capturing to pure subsonic flow but non -conservation can simulate either subsonic-supersonic or pure subsonic flow.
Hence it is recommended to use conservation form of equations for simulation since it covers the entire spectrum of possibility of flow contour in a nozzle.