Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
AIM:
To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle and to derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique and also to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
OBJECTIVE:
To show the plots of the following;
1. Steady-state distribution of primitive variables inside the nozzle.
2. Time-wise variation of the primitive variables.
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
4. Comparison of Normalized mass flow rate distributions of both forms.
CONVERGENT-DIVERGENT NOZZLE:
Schematic for subsonic-supersonic isentropic nozzle flow
QUASI ONE-DIMENSIONAL FLOW:
Quasi one-dimensional flow means that the cross-section is allowed to vary smoothly but flow quantities are changing in one direction only that is in X-direction.In this study we will simulate quasi one-dimensional flow in nozzle. We are assuming flow is steady and isentropic. The fluid is flowing through a convergent-divergent nozzle as shown in above figure.
EXPLANATION:
The flow at the nozzle comes from a reservoir where the pressure and temperature are denoted by `P_0`and `T_0`. The cross-sectional area of the reservoir is large(`A=oo`) and hence the velocity is very small (`V=0`).Thus `P_0`and `T_0` 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`. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat section and supersonic at the divergent section.The sonic flow is `(M=1)`at the throat that is equal to the speed of the sound . `P^ast,T^ast,V^ast (V^ ast=a^ast) and A^ ast`are the values denoted at the throat section.
Steps to be followed to solve the 1D supersonic nozzle flow simulation using Maccormack's method:
1. The governing equations should be in Non-dimensionalized form suitable for the time marching solutions of quasi-one dimensional flow.(Both conservative and Non-conservative part)
2. The finite difference equations derived by the Maccormack's technique as applied to this problem will be set up.
3. Other details such as calculation of timestep and applying of initial and boundary conditions , Nozzle shape for the numerical solution will be formulated.
STEP 1:NON-DIMENSIONALIZED GOVERNING EQUATIONS
Non-conservative form:
1. Continuity Equation
2. Momentum Equation
3. Energy Equation:
Conservative form:
1. Continuity Equation
2. Momentum Equation
3. Energy Equation
STEP 2: FINITE DIFFERENCE EXPRESSIONS USING MACCORMACK'S EXPLICIT TECHNIQUE FOR THE NUMERICAL SOLUTION
To implement a finite-difference solution, we divide the x-axis along the nozzle into a number of discrete grid points.
Grid point distribution along the nozzle
REFERENCE: CFD BY JOHN D ANDERSON
The first grid point labeled as point 1 is assumed to be in the reservoir. The points are evenly distributed along the X-axis, with `Delta x`denoting the spacing between grid points. The last point namely at that the nozzle exit is denoted by N; we have a total number of N grid points distributed along the axis. Point 'i' is simply an arbitrary grid point, with points 'i-1'and 'i+1'as the adjacent points. MacCormack's technique is a predictor-corrector method. In the time marching approach, remember that we know that the flow-field variables at time t, and we use the finite difference equations to solve explicitly for the variables at time `t+Deltat`.
Non-Conservative form
Predictor Step:-
First of all, we consider the Predictor step and then we set up the spatial derivatives as forward differences.
Continuity Equation:
Momentum Equation:
Energy Equation:
Using the above values we will obtain the predicted values of `rho,V and T`.
The above values are First order accurate only. So we have to do the corrector step to derive the second order accurate.
Corrector Step:-
In the corrector step we first obtained the predicted values of the time derivative at time `t+Deltat,((del rho)/(delt))_(i,j) ^(t+ trianglet)`By substituting the predicted values of `rho,V and T` into the right side of the continuity equation, replacing the spatial derivatives with backward or rearward differences.
Similarly Momentum and Energy equations will be,
Using above equations,the average time derivatives are;
Then we can substitute the values into the below formula to find the corrected values of the flow field variables at time `t+Deltat`i.e..next time steps
Conservative form
The equations for quasi one-dimensional flow can be expressed in generic form.Let us define the elements of the solutions vector U, the flux vector F and the source term J as follows;
With these elements, we can rewrite the non-dimensionalized continuity, momentum and energy equations in the format as follows,
To obtain the primitive variables `(rho,V,T,p,etc.)`,we must decode the elements `U_1,U_2,and U_3`as follows. From the definitions of `U_1,U_2, and U_3`given above,we have
In above Energy equation, we know that `e'=T'`, or
Therefore, After we obtain the solution vectors `U_1,U_2, and U_3`at each time step from the numerical solution. We can immediately calculate the corresponding primitive variables at each time step `rho',V',T', and p'`from above eqns.
Pure form of flux terms and source term
The flux vectors should be in terms of solution vectors as follows to avoid the numerical difficulties occasionally arise when the primitive variables are used to construct `F_1,F_2, and F_3`.
Predictor step:-
Here, we set up the spatial derivatives as forward differences.
`((delU_1)/(del t'))_i^(t')=-((F_1)_(i+1)-(F_1)_i)/(Delta x')`
`((delU_2)/(del t'))_i^(t')=-((F_2)_(i+1)-(F_2)_i)/(Delta x')+J_2`
`((delU_3)/(del t'))_i^(t')=-((F_3)_(i+1)-(F_3)_i)/(Delta x')`
Then we proceed to find the predicted values of `U_1,U_2 and U_3` denoted by the barred quantities. For example , we are going to calculate predicted values at grid point i=16;
From the above equations , we can decode the primitive variables of the predicted values. For example at node i=16,
The above Equation for the predicted `rho'and T'`will be needed on the corrector step. Now update the flux vector `F_1,F_2,F_3`with the new density, temperature and velocity.
Corrector step:-
The predicted time derivatives of `U_1,U_2 and U_3`are obtained above and correcter time derivatives can be find using rearward differences for the x derivatives.
`(bar(delU_1)/(del t'))_i^(t'+Delta t')=-((barF_1)_(i)-(barF_1)_(i-1))/(Delta x')`
`(bar(delU_2)/(del t'))_i^(t'+Delta t')=-((barF_2)_(i)-(barF_2)_(i-1))/(Delta x')+(1/gammabarrho' barT' (delA')/(del x'))`
`(bar(delU_3)/(del t'))_i^(t'+Delta t')=-((barF_3)_(i)-(barF_3)_(i-1))/(Delta x')`
Then the average time derivatives are formed as follows;
So, the final corrected values of `U_1,U_2,U_3`will become,
Finally, the corrected values of the primitive variables are obtained by decoding `U_1,U_2 and U_3`obtained above.
Step 3: Setting Time step calculation,initial and boundary conditions
Time step can be calculated using the following formula:-
Where,
C= courant number,
a= speed of sound,
V= velocity,
`Deltax`=Distance between two consecutive points.
Non-conservative form
Boundary conditions:-
Study of boundary conditions: Subsonic inflow and supersonic outflow
At the subsonic inflow boundary, we must stipulate the values of two dependent flow-field variables, whereas the value of one other variable must be allowed to float.
But at the Supersonic outflow boundary,all variables must be allowed to float at this boundary.
At Subsonic inflow :-
Since point 1 is viewed as essentially the reservoir, we stipulate the density and temperature at point 1 to be their respective stagnation values `rho_o andT_o` respectively . These are held fixed, independent of time.
At Supersonic outflow:-
Nozzle shape:-
Initial conditions:-
Conservative form
Nozzle shape:-
Initial conditions:-
Boundary condition at inflow:-
Boundary condition at outflow:-
MATLAB CODE AND EXPLANATION:-
Main Code:
clear all
close all
clc
%inputs
n=input('Number of grid points:');%number of nodes
x=linspace(0,3,n); %grids
dx=x(2)-x(1); %grid size
gamma=1.4; %adiabatic constant
nt=input('Number of time steps:');%number of time steps
% area
a=1+(2.2*(x-1.5).^2);
%courant number
c=0.5;
% calling the function
call=input('Type 1 for Non-conservative form and 2 for conservative form: ');
if call==1
%Non-conservative form
[rho,v,t,P,M,m,rho_t,v_t,t_t,P_t,M_t,m_t]=non_conservative_form(n,nt,x,dx,gamma,a,c);
%Plots for Non-conservative form
%(Mach number,Presssure,temperature,Density,Mass flow rate and Velocity Vs distance at 1400 time steps)
figure(1)
plot(x,M)
title_text=sprintf('Non-Conservative form \n Mach number vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Mach number')
figure(2)
plot(x,P)
title_text=sprintf('Non-Conservative form \n Pressure vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Pressure')
figure(3)
plot(x,t)
title_text=sprintf('Non-Conservative form \n Temperature vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Temperature')
figure(4)
plot(x,rho)
title_text=sprintf('Non-Conservative form \n Density vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Density')
figure(5)
plot(x,m)
title_text=sprintf('Non-Conservative form \n Mass flow rate vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Mass flow rate')
figure(6)
plot(x,v)
title_text=sprintf('Non-Conservative form \n Velocity vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Velocity')
%Plotting time wise variation of variables
%(Mach number,Presssure,temperature,Density,Mass flow rate and Velocity Vs Number of time steps)
nt=1:1400;
figure(7)
plot(nt,M_t)
title_text=sprintf('Non-Conservative form \n Mach number vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Mach number')
figure(8)
plot(nt,P_t)
title_text=sprintf('Non-Conservative form \n Pressure vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Pressure')
figure(9)
plot(nt,t_t)
title_text=sprintf('Non-Conservative form \n Temperature vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Temperature')
figure(10)
plot(nt,rho_t)
title_text=sprintf('Non-Conservative form \n Density vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Density')
figure(11)
plot(nt,m_t)
title_text=sprintf('Non-Conservative form \n Mass flow rate vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Mass flow rate')
figure(12)
plot(nt,v_t)
title_text=sprintf('Non-Conservative form \n Velocity vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Velocity')
end
if call==2
%Conservative form
[rho_c,v_c,t_c,P_c,M_c,m_c,rho_t_c,v_t_c,t_t_c,P_t_c,M_t_c,m_t_c]=conservative_form(n,nt,x,dx,gamma,a,c);
%Plots for conservative form
%(Mach number,Presssure,temperature,Density,Mass flow rate and Velocity Vs distance at 1400 time steps)
figure(1)
plot(x,M_c)
title_text=sprintf('Conservative form \n Mach number vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Mach number')
figure(2)
plot(x,P_c)
title_text=sprintf('Conservative form \n Pressure vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Pressure')
figure(3)
plot(x,t_c)
title_text=sprintf('Conservative form \n Temperature vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Temperature')
figure(4)
plot(x,rho_c)
title_text=sprintf('Conservative form \n Density vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Density')
figure(5)
plot(x,m_c)
title_text=sprintf('Conservative form \n Mass flow rate vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Mass flow rate')
figure(6)
plot(x,v_c)
title_text=sprintf('Conservative form \n Velocity vs Distance (No. of timestep=%d)',nt);
title(title_text)
xlabel('Distance')
ylabel('Velocity')
%Plotting time wise variation of variables
%(Mach number,Presssure,temperature,Density,Mass flow rate and Velocity Vs Number of time steps)
nt=1:1400;
figure(7)
plot(nt,M_t_c)
title_text=sprintf('Conservative form \n Mach number vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Mach number')
figure(8)
plot(nt,P_t_c)
title_text=sprintf('Conservative form \n Pressure vs Number of time steps ');
title(title_text)
xlabel('Number of time steps')
ylabel('Pressure')
figure(9)
plot(nt,t_t_c)
title_text=sprintf('Conservative form \n Temperature vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Temperature')
figure(10)
plot(nt,rho_t_c)
title_text=sprintf('Conservative form \n Density vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Density')
figure(11)
plot(nt,m_t_c)
title_text=sprintf('Conservative form \n Mass flow rate vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Mass flow rate')
figure(12)
plot(nt,v_t_c)
title_text=sprintf('Conservative form \n Velocity vs Number of time steps');
title(title_text)
xlabel('Number of time steps')
ylabel('Velocity')
end
%Mass flow rate distribution inside the nozzle
%Analytical form
m_an=linspace(0.579,0.579,n);
figure(14)
plot(x,m)
hold on
plot(x,m_c)
hold on
plot(x,m_an)
title_text=sprintf('Mass flow rate distribution inside the nozzle');
title(title_text)
xlabel('Distance')
ylabel('Mass flow rate')
legend('Non-Conservative form','Conservative form','Analytical form')
Explanation:
1. Here, I have given code to input the calling of function that is function 1 or 2
2. Used 'if' function to call the non conservative or conservative form function.
3. Plotted the graphs for steady state variation of primitive variables , time wise variation of primitive variables and Comparison of Normalized mass flow rate distributions of both forms with Analytical solution of mass flow rate.
Function code for Non-conservative form:
function[rho,v,t,P,M,m,rho_t,v_t,t_t,P_t,M_t,m_t]=non_conservative_form(n,nt,x,dx,gamma,a,c)
%calculating initial profiles in terms of x
%density
rho=1-0.3146*x;
%temperature or internal energy
t=1-0.2314*x;
%velocity
v=(0.1+1.09*x).*t.^0.5;
%time step calculation
dt=min(c*(dx./(sqrt(t)+v)));
%time marching loop
for k=1:nt
rho_old=rho;
v_old=v;
t_old=t;
%Predictor step
for j=2:n-1
drhodx=(rho(j+1)-rho(j))/dx;
dvdx=(v(j+1)-v(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%drho_dt_p(j)=-rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
%Continuity equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%Momentum equation
dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
%Energy equation
dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+(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 step
for j=2:n-1
drhodx=(rho(j)-rho(j-1))/dx;
dvdx=(v(j)-v(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/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
%Energy equation
dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+(v(j)*dlogadx));
end
%computing the average time derivatives
drhodt_avg=0.5*(drhodt_p+drhodt_c);
dvdt_avg=0.5*(dvdt_p+dvdt_c);
dtdt_avg=0.5*(dtdt_p+dtdt_c);
%solution updates
for i=2:n-1
rho(i)=rho_old(i)+drhodt_avg(i)*dt;
v(i)=v_old(i)+dvdt_avg(i)*dt;
t(i)=t_old(i)+dtdt_avg(i)*dt;
end
%Applying the boundary conditions
%inlet
v(1)=2*v(2)-v(3);
%outlet
v(n)=2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
t(n)=2*t(n-1)-t(n-2);
%properties
%Mach number
M=(v./sqrt(t));
%pressure
P=rho.*t;
%mass flow rate
m=rho.*v.*a;
%Throat datas
throat=find(a==1);
rho_t(k)=rho(throat);
v_t(k)=v(throat);
t_t(k)=t(throat);
P_t(k)=P(throat);
M_t(k)=M(throat);
m_t(k)=m(throat);
%Using if else condition to capture the variation of non dimensional mass flow rate at specific time steps and plotting it
if k==50
%plotting mass flow rate at time step 50
figure(13)
plot(x,m,'linewidth',2)
hold on
elseif k==100
%plotting mass flow rate at time step 100
plot(x,m,'linewidth',2)
elseif k==150
%plotting mass flow rate at time step 150
plot(x,m,'linewidth',2)
elseif k==200
%plotting mass flow rate at time step 200
plot(x,m,'linewidth',2)
elseif k==700
%plotting mass flow rate at time step 700
plot(x,m,'linewidth',2)
title('Plot of mass flow rate at different time steps in Non-conservative form')
xlabel('Length')
ylabel('Mass flow rate')
legend('50th time step','100th time step','150th time step','200th time step','700th time step')
end
end
endExplanation:
1. Initial conditions , time step calculation and Boundary conditions for inlet and outlet are inputed here.
2. Then Using Maccormack's technique ,we have found the predicted and corrected values of primitive variables.
3. From the new solutions , By applying the boundary conditions we calculate the value for the next timestep.
4. Then the 'for' loop goes on till the given value of "2 to n-1"
Function code for Conservative form:
function[rho_c,v_c,t_c,P_c,M_c,m_c,rho_t_c,v_t_c,t_t_c,P_t_c,M_t_c,m_t_c]=conservative_form(n,nt,x,dx,gamma,a,c)
%Calculating initial profile
rho_initial=ones(1,n);
t_initial=ones(1,n);
for i=1:n
if(x(i)>=0&&x(i)<=0.5)
rho_initial(i)=1;
t_initial(i)=1;
elseif (x(i)>=0.5&&x(i)<=1.5)
rho_initial(i)=1-0.366*(x(i)-0.5);
t_initial(i)=1-0.167*(x(i)-0.5);
elseif (x(i)>=1.5&&x(i)<=3)
rho_initial(i)=0.634-0.3879*(x(i)-1.5);
t_initial(i)=0.833-0.3507*(x(i)-1.5);
end
end
%velocity profile
v=0.59./(rho_initial.*a);
%initial solution vectors
u1=rho_initial.*a;
u2=rho_initial.*a.*v;
u3=u1.*((t_initial./(gamma-1))+((gamma/2)*(v.^2)));
P=rho_initial.*t_initial;
%time step calculation
dt=min(c*(dx./(sqrt(t_initial)+v)));
rho=rho_initial;
t=t_initial;
%time marching loop
for k=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
%flux terms
f1=u2;
term=((gamma-1)/gamma)*(u3-(gamma/2)*((u2.^2)./u1));
f2=((u2.^2)./u1)+term;
f3=((gamma.*u2.*u3)./u1)-((gamma*((gamma-1)/2))*((u2.^3)./(u1.^2)));
%predictor step
for j=2:n-1
%source term
%j2(j)=(rho(j)*t(j)*(a(j+1)-a(j)))/dx;
dadx=(a(j+1)-a(j))/dx;
j2(j)=(rho(j)*t(j)*dadx)/gamma;
%continuity equation
du1dt_p(j)=-((f1(j+1)-f1(j))/dx);
%momentum equation
du2dt_p(j)=-((f2(j+1)-f2(j))/dx)+j2(j);
%energy equation
du3dt_p(j)=-((f3(j+1)-f3(j))/dx);
%solution update
u1(j)=u1_old(j)+du1dt_p(j)*dt;
u2(j)=u2_old(j)+du2dt_p(j)*dt;
u3(j)=u3_old(j)+du3dt_p(j)*dt;
end
%solution update by solving primitive variables
rho=u1./a;
v=u2./u1;
t=(gamma-1)*((u3./u1)-((gamma/2)*(v.^2)));
%updating flux terms f1,f2,f3
%flux terms in terms of solution vectors
f1=u2;
term=((gamma-1)/gamma)*(u3-(gamma/2)*((u2.^2)./u1));
f2=((u2.^2)./u1)+term;
f3=((gamma.*u2.*u3)./u1)-(gamma*((gamma-1)/2)*(u2.^3)./(u1.^2));
%corrector step % spatial derivatives in backward differencing
for j=2:n-1
%source term
%j2(j)=(rho(j)*t(j)*(a(j)-a(j-1)))/dx;
dadx=(a(j)-a(j-1))/dx;
j2(j)=(rho(j)*t(j)*dadx)/gamma;
%continuity equation
du1dt_c(j)=-((f1(j)-f1(j-1))/dx);
%momentum equation
du2dt_c(j)=-((f2(j)-f2(j-1))/dx)+j2(j);
%energy equation
du3dt_c(j)=-((f3(j)-f3(j-1))/dx);
end
%Average time derivatives
du1dt_avg=0.5*(du1dt_p+du1dt_c);
du2dt_avg=0.5*(du2dt_p+du2dt_c);
du3dt_avg=0.5*(du3dt_p+du3dt_c);
%final solution or next timesteps update
for i=2:n-1
u1(i)=u1_old(i)+(du1dt_avg(i)*dt);
u2(i)=u2_old(i)+(du2dt_avg(i)*dt);
u3(i)=u3_old(i)+(du3dt_avg(i)*dt);
end
%Primitive variables are calculated from updated solution vector
%'_c'is the conservative form
rho_c=u1./a;
v_c=u2./u1;
t_c=(gamma-1)*((u3./u1)-((gamma/2)*(v.^2)));
%Applying the boundary conditions
%inlet
u1(1)=rho_c(1)*a(1);% where rho_c(1)=1 at node 1
u2(1)=2*u2(2)-u2(3);
v(1)=u2(1)/u1(1);
u3(1)=u1(1)*(t(1)/(gamma-1)+(gamma/2)*(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);
%properties
%Mach number
M_c=(v_c./(t_c.^0.5));
%pressure
P_c=rho_c.*t_c;
%Mass flow rate
m_c=rho_c.*a.*v_c;
%Throat datas
throat=find(a==1);
rho_t_c(k)=rho_c(throat);
v_t_c(k)=v_c(throat);
t_t_c(k)=t_c(throat);
P_t_c(k)=P_c(throat);
M_t_c(k)=M_c(throat);
m_t_c(k)=m_c(throat);
%Using if else condition to capture the variation of non dimensional mass flow rate at specific time steps and plotting it
figure(13)
hold on
if k==50
%plotting mass flow rate at time step 50
plot(x,m_c,'linewidth',2)
hold on
elseif k==100
%plotting mass flow rate at time step 100
plot(x,m_c,'linewidth',2)
elseif k==150
%plotting mass flow rate at time step 150
plot(x,m_c,'linewidth',2)
elseif k==200
%plotting mass flow rate at time step 200
plot(x,m_c,'linewidth',2)
elseif k==700
%plotting mass flow rate at time step 700
plot(x,m_c,'linewidth',2)
title('Plot of mass flow rate at different time steps in conservative form')
xlabel('Length')
ylabel('Mass flow rate')
legend('50th time step','100th time step','150th time step','200th time step','700th time step')
end
end
end
Explanation:
1. The conservative function is also the same as the non-conservative function.
2. Initial conditions , time step calculation and Boundary conditions for inlet and outlet are inputed here.
3. Then Using Maccormack's technique ,we have found the predicted and corrected values of primitive variables.
4. From the new solutions , By applying the boundary conditions we calculate the value for the next timestep.
5. Then the 'for' loop goes on till the given value of "2 to n-1"
OUTPUT:
Non-Conservative form:
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
Conservative form:
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms
Grid independence test:
| Cases/variables | Density | Pressure | Temperature | Mach number |
| Exact Analytical solution | 0.634 | 0.528 | 0.833 | 1 |
| Non conservative form numerical results(Nodes 31) | 0.638 | 0.534 | 0.836 | 0.999 |
| Conservative form numerical results(Nodes 31) | 0.650 | 0.546 | 0.840 | 0.981 |
| Conservative form numerical results(Nodes 61) | 0.638 | 0.533 | 0.835 | 0.999 |
From the above result values of Density, Temperature , Pressure and Mach number of Analytical, conservative and Non-conservative form,we can see that the Non conservative form results are really closer to the analytical form results and also we can see that for the conservative form results (Nodes 61) , the steady state numerical results are slightly closer to the analytical values of exact solution. For doubling the nodes for conservative form, it gives us closer values to exact solution.
Observations:
1. From the plotting of graphs, we can see that the the mach number increasing from the inlet to the outlet.wkt, At inlet the mach number is less than 1 i.e..Subsonic.At throat M=1, sonic and at outlet, M>1 i.e..Supersonic.At the outlet, the velocity increases in the divergent section so, the temperature and pressure decreases.As the pressure decreases the density also decreases.
2. Observing the time-wise variation of primitive variables, at first the largest changes takes place in the starting timesteps after which the steady state values are obtained at the final stages as later the time derivatives rapidly grows smaller over a finite number of timesteps.
3. While comparing the steady state variation of the mass flow rate obtained from the Non-conservation form with the conservative form , we can see that the conservation form is much more satisfactory than the non conservative form.The conservative form gives a distribution that is much closer to being a constant. But the non conservative form results have a sizeable variation with some huge oscillations at both the inflow and outflow boundaries. Also, the conservative form is much closer to the exact solution.
Conclusion:
1. Mass flow rate distribution is acceptable in conservative form solutions.
2. Non-conservative form is easy for programming convenience because they are solved in terms of primitive variables in the governing equations and their values are obtained directly from the Maccormack's technique. But in conservation form, the equations are expressed in terms of flux variables,first we have to solve solution vector and decoding it to find the primitive variables and updating the flux terms to go into the next spatial step. which makes things complicated.
3. For Accuracy , the Non-conservative form is better due to the direct calculation of primitive variables while in conservative form increased number of steps will lead to more errors.
4. In Non-conservative form ,Behavior of the residuals are good as compared to the conservation form. The residuals decays faster in non conservative form when compared to the conservative form.The amount by which the residual decay is often used as an index 'Quality' of the numerical algorithm.So, the non conservative form is better.
5. The simulation time for the non conservative form is shorter than the conservative because of the primitive variables converge to a numerical steady-state faster when using the non-conservative form.