Numerical analysis of Quasi 1D subsonic-supersonic flow Through convergent-divergent nozzle using MacCormack's Technique
OBJECTIVE:
1. To solve Quasi-1D isentropic subsonic supersonic flow throught convergent divergent nozzle using Maccormack's Technique in Matlab.
2.To solve using both conservative and non-conservative form of governing equations and compare their results.
3. To determine mininmum number of cycles required in order for convergence.
Computational Flow domain:
Nozzle: A nozzle is often a pipe or tube of varying cross sectional area, and it can be used to direct or modify the flow of a fluid (liquid or gas). Nozzles are frequently used to control the rate of flow, speed, direction, mass, shape, and/or the pressure of the stream that emerges from them. In a nozzle, the velocity of fluid increases at the expense of its pressure energy.
We are considering unsteady, isentropic flow through a convergent- divergent nozzle as shown in figure below.
The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by `p_o` and `T_o` respectively.The cross- sectional area of the reservoir is large(theoretically, `A rightarrow oo`), and hence the velocity is very small (`V rightarrow 0`). Thus,`p_o` and `T_o` are the stagnation values, or total pressure and total temperature,respectively. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity and mach number are denoted by `p_e`, `T_e`,`V_e` and `M_e` respectively. The flow is locally subsonic in the convergent section of the nozzle,sonic at throat area(minimum area) and supersonic at the divergent section. The sonic flow(M=1) at the throat means that the local velocity at this location is equal to the local speed of sound (At throat velocity, `V^* = a^*`). The area of the sonic throat is denoted by `A^*`
Assumption: We assume Quasi-one-dimensional flow,i.e. we assume that at a give section, where the cross sectional area is A, the flow properties are uniform across that section.Hence, although the area of the nozzle changes as a function of distance along the nozzle,x, and therefore in reality the flow field is two-dimensional, we make the assumption that the flow properties vary only with x, this is tantamount to assuming uniform flow properties across any given cross section.
MacCormack's Technique:
The MacCormack method is an explicit finite-difference technique which is second order accurate in both space and time. This method is a popular technique as it provides satisfactory results for many fluid flow applications.MacCormack method is a variant of the Lax- Wendroff but is much simpler in its application because there is no need to evaluate the second order time derivatives with was the case for Lax-Wendroff method..MacCormack method is a two-step predictor- corrector sequence is used with forward differences on the predictor and with rearward differences on the corrector is a second order accurate method.
As we assume that the flow field at each grid point in time t is known, the flow field variables at each grid point at time `t+ Delta t` can be calculated using maccormack method in a predictor-corrector sequence.
In maccormack method,this can be obtained from,
`rho_(i,j)^(t+Deltat) = rho_(i,j)^t + (frac{del rho}{del t})_(av)* Deltat`
where, `(frac{delrho}{delt})_(av)= 0.5* [(frac{delrho}{delt})_(predicted)+ (frac{delrho}{delt})_(corrected)]`
Predicted step:
All the flow variables at time t are known, and we determine a predicted value of density, `(barrho)^(t+Deltat)` from the first two terms of Taylor series as follows,
`(barrho)^(t+Deltat) = rho_(i,j)^t + (frac{delrho}{delt})_(i,j)^t * Deltat`
here, `(frac{delrho}{delt})_(i,j)^t` can be obtained from the continuity equation and replacing the spatial derivatives on the right side with forward differences.
Corrected step:
In the corrected step, we first obtain a predicted value of time derivative at time `t+Deltat`,`(frac{bar(delrho)}{delt})_(i,j)^(t+Deltat)` (also known as corrected time derivative), by substituting the predicted values of density (`(barrho)_(i,j)^(t+Deltat)`),velocity (`(barv)_(i,j)^(t+Deltat)`) into the right side of the continuity equation, also replacing the spactial derivatives with rearward differences.
This sequence is repected at all grid points to obtain the density throught the flow field at time `t+Deltat`.To calculate velocity(v), temperature(t) and energy(e),the same technique is used and utilizing the momemtum and energy equations to obtain the average time derivatives via the predictor-corrector sequence, using forward differences on the predictor and rearward differences on the corrector.
1. Non- conservative form
Governing Equations:
Continuity equation:
`frac{delrho}{delt}= -rho*frac{delv}{delx}-rho*v*frac{del(ln A)}{dx}-v*frac{delrho}{delx}`
Momemtum equation:
`frac{delv}{delt}=-v*frac{delv}{delx}-frac{1}{gamma}*(frac{delT}{delx}+frac{T}{rho}*frac{delrho}{delx})`
Energy Equation:
`frac{delT}{delt}=-v*frac{delT}{delx}-(gamma-1)*T*(frac{delv}{delx}+v*frac{del(lnA)}{delx})`
For our solution of the quasi one dimensional nozzle flow of a calorically perfect gas, temperature(T) is a fundamental variable, therefore it is convinient to deal with it as the primary dependent variable in the energy equation.For a calorically perfect gas, energy, `e=c_v*T`.
Here, Density(`rho`), velocity (v), Temperature(T), nozzle distance(x), pressure (P), crossectional area (A), time(t) are expressed in terms of Non-dimensional variables.
Boundary conditions:
Since the above governing equations are hyperbolic in nature,therefore there are two real characteristic lines.The method of characteristics tells us that at a boundary where one characteristic propagates into the domian, then the value of one dependent flow field variable must be specified at that boundary, and if one characteristic line propagates out of the domian, then the value of another dependent flow field variable must be allowed to float at the boundary(i.e. it must be calculated in steps of time as a function of timewise solution of flow field).
At subsonic inflow boundary:
Streamline moving into the domian at this boundary stipulates that the value of a second flow feild variable must be specified.Therefore at subsonic inflow boundary, we must stipulate the values of two dependent flow field, whereas the value of one other variable must be allowed to float.
`rho_1=1,T_1=1` (fixed, independent of time)
`v_1=2v_2-v_3`
At supersonic outflow boundary:
We have both characteristics propagating out of the domian, so does the streamline at point N. Therefore,there are no flow field variables which require their values to be stipulated at the supersonic outflow boundary.Therefore all variables must be allowed to float at this boundary.
`v_n= 2v_(n-1)-v_(n-2)`
`rho_n= 2rho_(n-1)- rho_(n-2)`
`T_n= 2T_(n-1)-T_(n-2)`
Nozzle shape and initial conditions:
For nozzle shape, we specify a parabolic area distribution
`A= 1+2.2*(x-1.5)^2 , 0<=x<=3`
In present problem, we know that density and temperature decrease and velocity increases as the flow expands through the nozzle. Hence we choose initial conditions that qualitatively behave in same fashion.
`rho= 1-0.3146x`
`T= 1-0.2314x`
`v= (0.1+1.09x)T^(1/2)`
2. Conservative Form:
Governing Equation:
Continuity Equation:
`frac{del(rhoA)}{delt}= -frac{del(rhoAv)}{delx}`
where, we assume
solution vector, `U_1= rhoA`
flux vector, `F_1= rhoAv`
Momemtum Equation:
`frac{del(rhoAv)}{delt}+ frac{del(rhoAv^2+(frac{1}{gamma}*PA))}{delx}= frac{1}{gamma}*P*frac{delA}{delx}`
where, we assume
solution vector, `U_2= rhoAv`
flux vector, `F_2= rhoAv^2+ frac{1}{gamma}PA`
source term, `J_2= frac{1}{gamma}Pfrac{delA}{delx}`
Energy Equation:
`frac{del(rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)A)}{delt}+frac{del(rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)vA+PAv)}{delx}=0`
where, we assume
solution vector, `U_3= rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)A`
flux vector, `F_3=rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)vA+PAv`
Therefore, above governing equations can be rewritten as
`frac{delU_1}{delt}=-frac{delF_1}{delx}`
`frac{delU_2}{delt}= -frac{delF_2}{delx}+J_2`
`frac{delU_3}{delt}= -frac{delF_3}{delx`
In Conservative form of governing equations, the dependent variables( the variables for which we directly obtain numbers) are not the primitive variables.Therefore to obtain the primitive variables from above governing equations,
`rho= frac{U_1}{A}, v= frac{U_2}{U_1}, P= rhoT`
`T=e=(gamma-1)(frac{U_3}{U_1}-frac{gamma}{2}v^2)`
Pure Form of Flux terms: The flux terms are represented in terms of solution vectors instead to primitive variables(`rho,v,T`) in order to obtain a stable solution.Therefore,
`F_1= U_2`
`F_2= frac{U_2^2}{U_1}+ frac{gamma-1}{gamma}(U_3-frac{gamma}{2}frac{U_2^2}{U_1})`
`F_3= gammafrac{U_2U_3}{U_1}-frac{gamma(gamma-1)}{2}frac{U_2^3}{U_1^2}`
Boundary Conditions:
At subsonic inflow boundary:
`U_(2(i=1))= 2U_(2(i=2))-U_(2(i=3))`
At supersonic outflow boundary:
`(U_1)_N= 2(U_1)_(N-1)-(U_1)_(N-2)`
`(U_2)_N= 2(U_2)_(N-1)-(U_2)_(N-2)`
`(U_3)_N= 2(U_3)_(N-1)-(U_3)_(N-2)`
Initial Conditions:
`rho=1,T=1,` For `0<=x<=0.5`
`rho= 1-0.366(x-0.5), T=1-0.167(x-0.5)` For `0.5<=x<=1.5`
`rho=0.634-0.3879(x-1.5), T=0.833-0.3507(x-1.5)` For `1.5<=x<=3.5`
`U_2` is physically the local mass flow rate, for initial conditions let us assum a constant mass flow rate, Therefore
`v= frac{U_2}{rhoA}= frac{0.59}{rhoA}`
Calculation of time step:
According to the Courant- Friedrichs-Lowry(CFL) criterion,
`(Deltat)_i^t = Cfrac{Deltax}{a_i^t+v_i^t`
where, `v` is the local flow velocity at a point in the flow and `a` is the local speed of sound
For a stable solution, `C<=1` , we assume C=0.5.
we consider the minimum time step calculated at each grid point,i.e.,
`Deltat= min(Deltat_1^t,Deltat_2^t,....,Deltat_N^t)`
Convergence criteria:
The convergence criteria used to determine minimum number of time steps required to attain the convergence is,
mass flow rate= `rho vA`= constant.
Tolerance: 1e-5.
Main Program:
% nozzle flow
close all
clear all
clc
Num_of_grid=[31 61];
for i=1:length(Num_of_grid)
if Num_of_grid(i)==31
n_31=Num_of_grid(i);
x_31= linspace(0,3,n_31);
dx_31= x_31(2)-x_31(1);
nt_31= 1400;
tic
[rho_non_31,T_non_31,v_non_31,dt_31,M_non_31,P_non_31,count_non_31,mfr_non_31,drhodt_av_th_31,dvdt_av_th_31,dTdt_av_th_31]=nonconservation_form(n_31,nt_31,x_31,dx_31);
time_non_31= toc;
tic
[rho_con_31,T_con_31,v_con_31,dt_con_31,M_con_31,P_con_31,count_con_31,mfr_con_31,dU1dt_av_th_31,dU2dt_av_th_31,dU3dt_av_th_31]= conservation_form(n_31,nt_31,x_31,dx_31);
time_con_31= toc;
else
if Num_of_grid(i)==61
n_61=Num_of_grid(i);
x_61= linspace(0,3,n_61);
dx_61= x_61(2)-x_61(1);
nt_61=1400;
tic
[rho_non_61,T_non_61,v_non_61,dt_61,M_non_61,P_non_61,count_non_61,mfr_non_61,drhodt_av_th_61,dvdt_av_th_61,dTdt_av_th_61]=nonconservation_form(n_61,nt_61,x_61,dx_61);
time_non_61=toc;
tic
[rho_con_61,T_con_61,v_con_61,dt_con_61,M_con_61,P_con_61,count_con_61,mfr_con_61,dU1dt_av_th_61,dU2dt_av_th_61,dU3dt_av_th_61]= conservation_form(n_61,nt_61,x_61,dx_61);
time_con_61= toc;
end
end
end
figure(1)
subplot(4,1,1)
plot(x_31,M_non_31,'r',x_61,M_non_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('non- dimensional mach no')
title('Non dimensional mach num variation of Non-conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,2)
plot(x_31,P_non_31,'r',x_61,P_non_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Pressure ratio')
title('Non dimensional pressure ratio variation of Non-conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,3)
plot(x_31,rho_non_31,'r',x_61,rho_non_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Density ratio')
title('Non dimensional Density ratio variation of Non-conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,4)
plot(x_31,T_non_31,'r',x_61,T_non_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Temperature ratio')
title('Non dimensional Temperature ratio variation of Non-conservation Form')
legend('31 grids','61 grids','location','best')
figure(7)
plot(x_31,mfr_con_31,'color','r','linewidth',2)
hold on
plot(x_31,mfr_non_31,'color','b','linewidth',2)
plot([0 3],[0.579 0.579],'color','m','linewidth',2)
xlabel('non dimensional length')
ylabel('non dimensional mass flow rate')
title('mass flow rate comparision of conservative form and non conservative form')
legend('conservative form','non conservative form','exact analytical solution')
figure(4)
subplot(4,1,1)
plot(x_31,M_con_31,'r',x_61,M_con_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('non- dimensional mach no')
title('Non dimensional mach num variation of conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,2)
plot(x_31,P_con_31,'r',x_61,P_con_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Pressure ratio')
title('Non dimensional pressure ratio variation of conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,3)
plot(x_31,rho_con_31,'r',x_61,rho_con_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Density ratio')
title('Non dimensional Density ratio variation of conservation Form')
legend('31 grids','61 grids','location','best')
subplot(4,1,4)
plot(x_31,T_con_31,'r',x_61,T_con_61,'b*');
hold on
xlabel('Non- dimensional length')
ylabel('Temperature ratio')
title('Non dimensional Temperature ratio variation of conservation Form')
legend('31 grids','61 grids','location','best')
1.nonconservation_form is a function used to solve the nozzle flow using non conservative form of governing equations and conservation_form is a function used to solve the nozzle flow using conservative form of governing equations.
Function nonconservation_form:
function [rho,T,v,dt,M,P,count,mfr,drhodt_av_th,dvdt_av_th,dTdt_av_th]= nonconservation_form(n,nt,x,dx)
%non- conservation
%initial condition
A=1+2.2.*(x-1.5).^2;
rho= 1-0.3146.*x;
T= 1-0.2314.*x;
v= (0.1+1.09.*x).*T.^0.5;
gamma=1.4;
C=0.5; %courant number
for i=1:length(x)
if x(i)==1.5
th=i;
end
end
Tol=1e-5;
count=0;
mfr= rho.*v.*A;
for k= 1:nt
%update values
rho_old= rho;
T_old= T;
v_old=v;
mfr_old=mfr;
%to find dt with CFL condition
for j= 1:n
d_t(j)= C*(dx/(T(j)^0.5+v(j)));
end
dt(k)= min(abs(d_t));
for j= 2:n-1
%predictor step
dvdx_p(j)= (v(j+1)-v(j))/dx;
dAdx_p(j)= (log(A(j+1))-log(A(j)))/dx;
drhodx_p(j)= (rho(j+1)-rho(j))/dx;
dTdx_p(j)= (T(j+1)-T(j))/dx;
%continuity eqn
drhodt_p(j)= -rho(j)*dvdx_p(j)-rho(j)*v(j)*dAdx_p(j)-v(j)*drhodx_p(j);
%momemtum eqn
dvdt_p(j)= -v(j)*dvdx_p(j)-(1/gamma)*(dTdx_p(j)+(T(j)/rho(j))*drhodx_p(j));
%energy eqn
dTdt_p(j)=-v(j)*dTdx_p(j)-(gamma-1)*T(j)*(dvdx_p(j)+v(j)*dAdx_p(j));
%predicted values
rho(j)=rho(j)+drhodt_p(j)*dt(k);
v(j)=v(j)+dvdt_p(j)*dt(k);
T(j)=T(j)+dTdt_p(j)*dt(k);
end
for i= 2:n-1
%corrected step
dvdx_c(i)= (v(i)-v(i-1))/dx;
dAdx_c(i)= (log(A(i))-log(A(i-1)))/dx;
drhodx_c(i)= (rho(i)-rho(i-1))/dx;
dTdx_c(i)= (T(i)-T(i-1))/dx;
%continuity eqn
drhodt_c(i)= -rho(i)* dvdx_c(i)-rho(i)*v(i)*dAdx_c(i)-v(i)*drhodx_c(i);
%momemtum eqn
dvdt_c(i)= -v(i)* dvdx_c(i)-(1/gamma)*(dTdx_c(i)+(T(i)/rho(i))*drhodx_c(i));
%energy eqn
dTdt_c(i)=-v(i)*dTdx_c(i)-(gamma-1)*T(i)*( dvdx_c(i)+v(i)*dAdx_c(i));
end
%to find avg derivative
for j=2:n-1
drhodt_av(j)= 0.5*(drhodt_p(j)+drhodt_c(j));
dvdt_av(j)= 0.5*(dvdt_p(j)+dvdt_c(j));
dTdt_av(j)= 0.5*( dTdt_p(j)+dTdt_c(j));
end
% final value
for j=2:n-1
rho(j)= rho_old(j)+drhodt_av(j)*dt(k);
v(j)= v_old(j)+ dvdt_av(j)*dt(k);
T(j)= T_old(j)+ dTdt_av(j)*dt(k);
end
%boundary conditions
% at subsonic inlet
v(1)= 2*v(2)-v(3);
%at supersonic 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);
%mass flow rate
mfr= rho.*v.*A;
if n==31
figure(5)
hold on
if k==1
plot(x,mfr,'linewidth',2)
else
if k==50
plot(x,mfr,'linewidth',2)
else
if k==100
plot(x,mfr,'linewidth',2)
else
if k==200
plot(x,mfr,'linewidth',2)
else
if k==300
plot(x,mfr,'linewidth',2)
else
if k==700
plot(x,mfr,'linewidth',2)
xlabel('Non dimensional distance through nozzle(x)')
ylabel('Non dimensional mass flow rate')
title('Non dimensional mass flow rate variation through nozzle at different no of timesteps')
legend('1 timestep','50 time steps','100 time steps','200 time steps','300 time steps','700 time steps')
end
end
end
end
end
end
end
%number of timesteps taken to converge
error_mfr=max(mfr-mfr_old);
if error_mfr>=Tol
count = count +1;
end
%residual
drhodt_av_th(k)=abs(drhodt_av(th));
dvdt_av_th(k)=abs(dvdt_av(th));
dTdt_av_th(k)=abs(dTdt_av(th));
M= v./(T.^0.5); %mach no
P= rho.*T;
rho_th(k)=rho(th);
T_th(k)= T(th);
M_th(k)= M(th);
P_th(k)= rho(th)*T(th);
end
%Analytical solution
rho_an= 0.634*ones(1,1400);
T_an= 0.833*ones(1,1400);
p_an=0.528*ones(1,1400);
M_an= ones(1,1400);
M= v./(T.^0.5); %mach no
P= rho.*T;
if n==31
figure(3)
plot(linspace(1,nt,nt),drhodt_av_th,'color','r','linewidth',2)
hold on
plot(linspace(1,nt,nt),dvdt_av_th,'color','b','linewidth',2)
plot(linspace(1,nt,nt),dTdt_av_th,'color','k','linewidth',2)
xlabel('Number of time steps')
ylabel('residual')
title('Timewise variation of the absolute values of time derivatives of non dimensional density,velocity and temperature at the nozzle throat')
legend('density time derivative','velocity time derivative','temperature time derivative')
end
if n==31
figure(8)
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_th,'r',linspace(1,nt,nt),rho_an,'k')
ylabel('density ratio')
title('Timewise variation of non dimensional density,velocity,Mach no and temperature at the nozzle throat')
legend('numerical solution','analytical solution','location','best')
hold on
subplot(4,1,2)
plot(linspace(1,nt,nt),T_th,'b',linspace(1,nt,nt),T_an,'k')
ylabel('temperature ratio')
legend('numerical solution','analytical solution','location','best')
subplot(4,1,3)
plot(linspace(1,nt,nt),M_th,'m',linspace(1,nt,nt),M_an,'k')
ylabel('non-dimensional mach no')
legend('numerical solution','analytical solution','location','best')
subplot(4,1,4)
plot(linspace(1,nt,nt),P_th,'g',linspace(1,nt,nt),p_an,'k')
xlabel('Number of time steps')
ylabel('pressure ratio')
legend('numerical solution','analytical solution','location','best')
end
end
Function conservation_form:
function [rho,T,v,dt,M,P,count,mfr,dU1dt_av_th,dU2dt_av_th,dU3dt_av_th]= conservation_form(n,nt,x,dx)
%conservation
gamma=1.4;
C=0.5;
for i=1:length(x)
if x(i)==1.5
th=i;
end
end
Tol=1e-5;
count=0;
%initial condition
A=1+2.2.*(x-1.5).^2;
for i= 1:n
if x(i)>=0 && x(i)<=0.5
rho(i)=1;
T(i)=1;
else
if x(i)>= 0.5 && x(i)<=1.5
rho(i)= 1-0.366*(x(i)-0.5);
T(i)= 1-0.167*(x(i)-0.5);
else
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
end
end
v= 0.59./(rho.*A);
e=T;
mfr= rho.*v.*A;
%solution vectors
U1= rho.*A;
U2= rho.*A.*v;
U3= rho.*(T/(gamma-1)+(gamma/2)*v.^2).*A;
for i=1:nt
U1_old= U1; %update values
U2_old= U2;
U3_old= U3;
rho_old=rho;
v_old=v;
T_old=T;
mfr_old=mfr;
% to determine dt with CFL condition
for j= 1:n
d_t(j)= C*(dx/(T(j)^0.5+v(j)));
end
dt(i)= min(abs(d_t));
%flux vectors pure form
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);
for j= 2:n-1
%source term
adx(j)= (A(j+1)-A(j))/dx;
J2(j)= (1/gamma)*(rho(j)*T(j)*adx(j));
end
%predicted step
for j= 2:n-1
dU1dt_p(j)= -(F1(j+1)-F1(j))/dx;
dU2dt_p(j)= -(F2(j+1)-F2(j))/dx + J2(j);
dU3dt_p(j)= -(F3(j+1)-F3(j))/dx;
U1(j)=U1(j)+ dU1dt_p(j)*dt(i);
U2(j)=U2(j)+ dU2dt_p(j)*dt(i);
U3(j)=U3(j)+ dU3dt_p(j)*dt(i);
end
%update predicted values
rho = U1./A;
v= U2./U1;
T= (gamma-1).*(((U3./U1)-(gamma/2).*v.^2));
e=T;
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);
for j=2:n-1
adx(j)= (A(j)-A(j-1))/dx;
J2(j)= (1/gamma)*(rho(j)*T(j)*adx(j));
end
%corrected step
for j= 2:n-1
dU1dt_c(j)= -(F1(j)-F1(j-1))/dx;
dU2dt_c(j)= -(F2(j)-F2(j-1))/dx + J2(j);
dU3dt_c(j)= -(F3(j)-F3(j-1))/dx;
end
%avg time derivative
dU1dt_av= 0.5*(dU1dt_p+dU1dt_c);
dU2dt_av= 0.5*(dU2dt_p+dU2dt_c);
dU3dt_av= 0.5*(dU3dt_p+dU3dt_c);
for j=2:n-1
U1(j)= U1_old(j)+ dU1dt_av(j)*dt(i);
U2(j)= U2_old(j)+ dU2dt_av(j)*dt(i);
U3(j)= U3_old(j)+ dU3dt_av(j)*dt(i);
end
%boundary condition
U2(1)= 2*U2(2)-U2(3);
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);
%final update of primitive terms
rho = U1./A;
v= U2./U1;
T= (gamma-1)*((U3./U1)-(gamma/2)*v.^2);
e=T;
%number of timesteps taken to converge
mfr= rho.*v.*A;
error_mfr=max(mfr-mfr_old);
if error_mfr>=Tol
count = count +1;
end
%mass flow rate at different time steps
if n==31
figure(6)
hold on
if i==1
plot(x,mfr,'linewidth',2)
else
if i==50
plot(x,mfr,'linewidth',2)
else
if i==100
plot(x,mfr,'linewidth',2)
else
if i==200
plot(x,mfr,'linewidth',2)
else
if i==300
plot(x,mfr,'linewidth',2)
else
if i==700
plot(x,mfr,'linewidth',2)
xlabel('Non dimensional distance through nozzle(x)')
ylabel('Non dimensional mass flow rate')
title('Non dimensional mass flow rate variation through nozzle at different no of timesteps')
legend('1 timestep','50 time steps','100 time steps','200 time steps','300 time steps','700 time steps')
end
end
end
end
end
end
end
%residual
dU1dt_av_th(i)=abs(dU1dt_av(th));
dU2dt_av_th(i)=abs(dU2dt_av(th));
dU3dt_av_th(i)=abs(dU3dt_av(th));
M= v./(T.^0.5); %mach no
P= rho.*T;
rho_th(i)=rho(th); %throat values
T_th(i)= T(th);
M_th(i)= M(th);
P_th(i)= rho(th)*T(th);
end
%Analytical solution
rho_an= 0.634*ones(1,1400);
T_an= 0.833*ones(1,1400);
p_an=0.528*ones(1,1400);
M_an= ones(1,1400);
M= v./(T.^0.5); %mach no
P= rho.*T;
if n==31
figure(2)
plot(linspace(1,nt,nt),dU1dt_av_th,'color','r','linewidth',2)
hold on
plot(linspace(1,nt,nt),dU2dt_av_th,'color','b','linewidth',2)
plot(linspace(1,nt,nt),dU3dt_av_th,'color','k','linewidth',2)
xlabel('Number of time steps')
ylabel('residual')
title('Timewise variation of the absolute values of time derivatives of U1,U2 and U3 at the nozzle throat')
legend('U1 time derivative','U2 time derivative','U3 time derivative')
end
if n==31
figure(9)
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_th,'r',linspace(1,nt,nt),rho_an,'k')
ylabel('density ratio')
title('Timewise variation of non dimensional density,velocity,Mach no and temperature at the nozzle throat')
legend('numerical solution','analytical solution','location','best')
hold on
subplot(4,1,2)
plot(linspace(1,nt,nt),T_th,'b',linspace(1,nt,nt),T_an,'k')
ylabel('temperature ratio')
legend('numerical solution','analytical solution','location','best')
subplot(4,1,3)
plot(linspace(1,nt,nt),M_th,'m',linspace(1,nt,nt),M_an,'k')
ylabel('non-dimensional mach no')
legend('numerical solution','analytical solution','location','best')
subplot(4,1,4)
plot(linspace(1,nt,nt),P_th,'g',linspace(1,nt,nt),p_an,'k')
xlabel('Number of time steps')
ylabel('pressure ratio')
legend('numerical solution','analytical solution','location','best')
end
endGrid Dependense Test:
Definition: If we increase the of grid points(or decrease the grid size) of our computational domain, we will obtain a more accurate solution. If we continue to increase the number of grid points until we reach a solution which is no longer sensitive to the number of grid points, Then we have achieved grid independence of our solution.
Grid dependense test have been performed for our both non conservative form and conservative form of governing equations with 31 and 61 number of grid points. Non- dimensional varaibles such as pressure ratio,density ratio, temperature ratio and Mach number along the nozzle flow have been plot in order to check grid independence.
Varaition of Non- Dimensional Mach number, pressure ratio, density ratio and temperature ratio along the nozzle of non-conservative form of governing equations:
Varaition of Non- Dimensional Mach number, pressure ratio, density ratio and temperature ratio along the nozzle of conservative form of governing equations:
1. For both conservative and non-conservative form of governing equations, If we increase the number of grid points from 31 to 61 number of grid points, There is negligible change in variation of non dimensional pressure ratio,density ratio , temperature ratio and mach number along the nozzle flow.Therefore we can consider our solution with 31 number of grid points to be grid indenpendent.
2. By observing the variation of mach number along the nozzle flow for both conservative and non conservative form, mach number increases from inlet and almost equals to 1 at throat of convergence divergence nozzle(at x=1.5) and further increases proceeding towards outet. In other words, the flow remains subsonic(mach number less than 1) at the convergent part of the nozzle and flow becomes supersonic(greater at 1) at divergent part of the nozzle.
3. By observing the variation of non- dimesional pressure ratio, density ratio and temperature ratio along the nozzle flow for both conservative and non conservative form, As the velocity increases from flow through nozzle, Pressure ratio, density ratio and temperature ratio decreases as the fluid(calorically perfect gas) passes throught the nozzle. The functionality of the nozzle can be observed from above plots.
2. Variation of mass flow rate along the nozzle flow at different time steps:
2.1 Non- conservative form:
2.2 Conservative Form:
For Non-conservative Form, mass flow rate variation at timestep 1 is behaving sinelike variations due to the initial conditions and parabolic variation of area.At 100 timestep, mass flow rate distribution has changed considerably from initial conditions due to the transient variations of flow field variables.However after 200 timesteps, mass flow rate distribution is beginning to settle down and after 700 timesteps, the mass flow rate is almost constant value of 0.587 throughout the nozzle.
For conservative Form, mass flow rate variation at timestep 1 is more flatter than non-conservative form because the initial conditions of conservative form were more close to the exact value.At 50 timesteps, mass flow rate distribution deviates away from the initial conditions.After 100 timesteps, the mass flow rate distribution begins to settle down and at 700 timesteps, mass flow rate distribution is almost constant value(around 0.585).
3. Timewise variation of Non-dimensional Mach no, pressure ratio,density ratio and temperature at nozzle throat:
3.1 Non-conservative Form:
3.2 Conservative Form:
For both conservative and non-conservative form, Non-dimensional Mach number and other primitive variables begins to settle down after 400 timesteps. After 700 timesteps, Primitive variables and Mach number is almost constant.
Howerever, For non-conservative form the primitive variables( pressure ratio, density ratio,temperature ratio) and Mach number is more close to exact analytical value.And conservative form is not close to exact analytical value comparetively.(As shown below)
At Throat,
Considering Primitive variables and Mach number, Non- conservation form is more close to the exact analytical value.
4.Residual Plots:
4.1 Timewise variation of absolute values of time derivatives of non-dimensional density, velocity and temperature at nozzle throat(Non-conservation form)
4.2 Timewise variation of absolute values of time derivatives of U1,U2 and U3(source terms) at the nozzle throat(conservation Form)
Residuals are the absolute values of the time derivatives of non dimensional pressure, density and temperature(in case of non-conservative) or source terms(in case of conservative form).For both conservative and non-conservative forms, the residuals ossilate at early time steps but for non-conservation the residuals settles down at around 600 timesteps where in case conservative form , the residual settles down only after 800 timesteps. After 1400 timesteps, residuals of non conservative form drops in the range of 1e-08 but the residuals of conservative form drops only upto 1e-05.Therefore, the non conservation form leads to smaller residuals.The amount by which the residuals decay is often used as an index of quality of the numerical algorithmn.
5. Steady state mass flow rate comparison of Non-conservation form and conservation form:
Comparing the steady state mass flow rate of non-conservation form and conservation form, mass fow rate distribution of conservative form is more horizontal and is more closer to the exact analytic value of mass flow rate. This is due to mass flow rate is one of the dependent variables in the conservative form of governing equations(U2),therefore mass flow rate is a primary results form these equations.In constrast, the dependent variables in non conservative form of govening equations are the primitive variables(pressure,density and temperature) and mass flow rate is obtained only as an secondary result.
6. Comparison of number of time steps taken for solution to converge:
6.1 Non-conservative Form
For 31 grid points, Number of timesteps taken for solution to converge= 670. , Simulation time= 3.1880 sec
For 61 grid points, Number of time taken for solution to converge= 1207., Simulation time= 0.0302sec ( lesser simulation time because some graphs were not plot for 61 grid points since 61 grid points solution were only used for grid dependence test).
6.2 Conservative Form
For 31 grid points, Number of timesteps taken for solution to converge= 932. , Simulation time= 2.9080 sec
For 61 grid points, Number of time taken for solution to converge= 941., Simulation time= 0.0464sec ( lesser simulation time because some graphs were not plot for 61 grid points since 61 grid points solution were only used for grid dependence test).
Therefore, at 31 grid points, Non-conservative form takes the least number of timesteps taken to converge and at 61 grid points, conservative form takes the least number of timesteps taken to converge.
Conclusion:
1. In non- conservative form of governing equations, the primitive variables(non dimesional pressure, density and temperature) are the dependent varaibles, Therefore in this case primitive varaibles are more close to the exact analytical solution.
2. In conservative form of governing equations, mass flow rate is one of the dependent variables of governing equations and primitive variables are obtained as a secondary result, Therefore in this case mass flow rate is more close to the exact analytical solution.
3. After 1400 timesteps, residuals of non conservative form drops in the range of 1e-08 but the residuals of conservative form drops only upto 1e-05.Therefore, the non conservation form leads to smaller residuals.The amount by which the residuals decay is often used as an index of quality of the numerical algorithmn.
4.At 31 grid points, Non-conservative form takes the least number of timesteps taken to converge and at 61 grid points, conservative form takes the least number of timesteps taken to converge.
5. By observing the variation of non- dimesional pressure ratio, density ratio and temperature ratio along the nozzle flow for both conservative and non conservative form, As the velocity increases from flow through nozzle, Pressure ratio, density ratio and temperature ratio decreases as the fluid(calorically perfect gas) passes throught the nozzle.