Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Introduction
Let us consider a quasi 1-D flow through convergent divergent nozzle which is initially at subsonic speed and exiting at supersonic speed. The flow is along the direction of x axis making it one dimensional flow. All the parameters and conditions vary along the length of the nozzle that is along x axis. Eventhough the area of cross section of nozzle changes as the function of distance along the nozzle lenght x, which means that in reality the flow field is of two dimensional, it is assumed that flow properties varies only along x and assuming uniform flow properties across any cross section, such a flow is called Quasi-one-dimensional flow.
Inlet parameters :
The temperature and pressure at inlet is denoted by P0 and T0 respectively, it is considered that thee flow is incoming from a reservoir and the reservoir conditions are taken as the paramenters for inlet conditons. The area of cross section is very large and it is considered as A -> ∞, this implies that the velocity in the reservoir is very small and is considered as V-> 0 . Since the velocity tends to zero the pressure and temperature is in their stagnation conditon which means that P0 and T0 is total pressure and total temperature respectively.
The flow from inlet (reservoir) to the throat is through a convergent nozzle section and the velocity in this section is subsonic in nature.
Outlet parameters :
The flow exits the the nozzle through the divergent section where it isentropically expands to supersonic speed. The exit pressue , temperature , velocity and mach number is denoted by Pe ,Te ,Ve and Me respectively.
Throat conditions:
At the throat area the flow has a mach number = 1. The area of cross-section of the nozzle is minimum at the throat of the nozzle. The sonic flow here means that the velocity at this location is equal to the local velocity of sound at that location. The velocity, pressure, temperature and area at the throat is denoted by V* , P* , T* and A* respectively .
Understanding the concept
We are here to observe and deduce what happens when a one dimensional quasi flow through a convergent divergent nozzle from subsonic to supersonic at the exit.
To see how the primitive variables such as pressure, density, temperature(energy), velocity and the mach number of the flow varies along the length of the nozzle.
In order to do that, the physical aspects of the flow is described in the form of equations and all the factors that affect the flow is taken into consideration.
The governing continuity, momentum and energy equations for this quasi one dimensional, steady state isentropic flow can be expressed as :
Continuity : ρ1 V1 A1=ρ2 V1 A1
Momentum : `p 1 A 1 + ρ 1 V 1^ 2 A 1 + ∫ p d A = p 2 A 2 + ρ 2 V 2 ^2 A 2`
Energy : `h1 + (V1^2)/2 = h2 + (V2^2)/2`
The above governing equations can be derived into two forms of equations:
>Non-conservation form
>Conservation form
which can be further solved to obtain the final solution for all the primitive variables (pressure, velocity, mach number, temperature, density, momentum) at all the nodes along the length of the convergent divergent nozzle using Macormack method.
Governing equations of Non conservation form -
- Continuity: `(del rho')/(delt') = -rho' (del V')/(delx') - rho' V' (del ln A')/(delx') - V (del rho')/(delx')`
- Momentum : `(del V')/ (del t') = -V' (del V')/ (del x') - 1/gamma ((del T')/ (del x') + (T')/(rho')(del rho')/ (del x') )`
- Energy : `(delT')/(delt') = -V' (delT')/(delx') - (gamma-1)T'[(delV')/(delx') + V' (dellog A')/(delx')]`
Governing equations of Conservation form -
- Continuity : `(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
- Momentum : `(del(rho'A'V'))/(delt') + (del(rho'A'V') +(1/gamma)p'A')/(delx') = 1/gamma p' (delA')/(delx')`
- Energy : `(del[rho'((e')/(gamma-1) + gamma/2 V')A'])/(delt') + (del[rho'((e')/(gamma-1) + gamma/2 V')V'A' + p'A'V'])/(delx') =0`
Code :
Main program code to call the functons of conservative and non conservative forms solver and to plot the comparison between the two and the analytical solution with respect to mass flow rate along the length of the nozzle.
The length of the nozzle, number of nodes and the total number of time steps to be solved is assigned here and those values are sent to the function.
All the primitives and their values at each node inside the nozzle is obtained for further comparison.
%main program
clear all
close all
clc
%inputs
%grid points
n = 31;
%length of the nozzle
L=3;
%total time steps
nt = 1400;
%Allocating grid points along length of the nozzle
x = linspace(0,L,n);
%Function call for non conserrvative form
[rhonc, Tnc, vnc, Mnc, mfnc] = non_conservative(n,L,nt);
%Function call for conserrvative form
[rhoc, Tc, vc, Mc, mfc] = conservative(n,L,nt);
%rho * a * v = constant
analytical_result = 0.579;
figure (9)
plot (x,mfnc)
hold on
plot (x, mfc)
hold on
plot (x,analytical_result,'.','color','black')
legend('Non-conservation form' , 'Conservation form','analytical result');
xlabel('x/L' );
ylabel('Mass flow rate');
Non - conservation form
The governing equations are converted into suitable partial diffrential equation so that the time marching soluton on quasi 1-D flow can be exectued using these equations.
Governing equations of Non conservation form in PDE form is given by -
Continuity: `(del rho')/(delt') = -rho' (del V')/(delx') - rho' V' (del ln A')/(delx') - V (del rho')/(delx')`
Momentum : `(del V')/ (del t') = -V' (del V')/ (del x') - 1/gamma ((del T')/ (del x') + (T')/(rho')(del rho')/ (del x') )`
Energy : `(delT')/(delt') = -V' (delT')/(delx') - (gamma-1)T'[(delV')/(delx') + V' (dellog A')/(delx')]`
where
`rho` is the density of the flow
A is the area of cross section of the nozzle at that node
V is the velocity of the flow
T is the temperature as the form of energy
`delt` time step
`gamma` is the ratio of of heat capacity
Non dimensional form of quantities are used in the governing equations since the values of these lies between 0 and 1 which makes it easier to present the result, hence this is preferred in nozzle flow problems.
The non dimensional forms are given by :
`T' = T/(To)`
`rho' = rho/(rhoo)`
`x' = x/L` where x' is the dimensionless component of length L
`V' = V/(ao)` where ao is the speed of sound in the reservoir
`t' = t/(L/(ao))` where t' is the dimensionless time
`A' = A/(Aast)` dimensionless area is given by the ration of area at that point to area at throat of thee nozzle.
initial conditions:
We know that, from start to end of the flow, density and temperature decreases and velocity of the flow increases. This relation is put into equation as shown below and the variation of these quantities are considered to be linear and changes as a function of x . Hence at t=0, we get the flow field variables as :
`rho = 1 - 0.3146*x`
`T = 1 -0.2314*x;`
`v = (0.1 + 1.09*x)*T^0.5;`
The area of cross section of the nozzle is such that it converges into a throat and diverges from throat to exit. Hence the area of cross section decreases at first and reaches it's minimum at throat and increases from there. Considering the throat to be at the center of the nozzle along the length of it, we get the equation as :
`a = 1 + 2.2 * (x - 1.5)^2;`
Time step is calculated considering that the system is stable and taking cfl number
`dt = cfl*(dx/(sqrtT+v))`
The time step is calculated for at every node and the minimum of the obtained result is used since all the flow-field variables converge to the same physical time.
After obtaining the respecctive governing equations and setting up the initial conditions, finite-differrence expressions pertaining to MacCormack's method is applied.
The solution is converged for the determined time length and provided total number of time steps, following steps are iterated until the solution is converged
i here defines the node at that time step
>Predictor step :
This step is basically the forward differencing of the governing equation
`drhodt(i)= -(rho(i) * (v(i+1) - v(i))/dx) - (rho(i) * v(i) * (log(a(i+1)) - log(a(i)))/dx) - (v(i) * (rho(i+1)- rho (i))/dx)`
`dvdt(i) = -(v(i)*(v(i+1) - v(i))/dx) - 1/gamma*((T(i+1) - T(i))/dx)+ ((T(i))/(rho(i))*(rho(i+1)- rho (i))/dx))`
`dTdt(i) = -(v(i)*(T(i+1) - T(i))/dx) - ((gamma -1)*T(i) * ((v(i+1) - v(i))/dx)+ (v(i) * (log(a(i+1)) - log(a(i)))/dx)`
>Corrector step :
This step is basically the backward differencing of the governing equation
`drhodt(i)= -(rho(i) * (v(i) - v(i-1))/dx) - (rho(i) * v(i) * (log(a(i)) - log(a(i-1)))/dx) - (v(i) * (rho(i)- rho (i-1))/dx)`
`dvdt(i) = -(v(i)*(v(i) - v(i-1))/dx) - 1/gamma*((T(i) - T(i-1))/dx)+ ((T(i))/(rho(i))*(rho(i)- rho (i-1))/dx))`
`dTdt(i) = -(v(i)*(T(i) - T(i-1))/dx) - ((gamma -1)*T(i) * ((v(i) - v(i-1))/dx)+ (v(i) * (log(a(i)) - log(a(i-1)))/dx)`
>Averaging :
Averaging is done by taking the average of the result from predictor step and corrector step
>Final step
In this step the averaged result is added with the initial value of the quantity at the start of that iteration
Boundary conditions
At inlet
Velocity = V(2) - V(3)
At outlet:
velocity: `v(n) = 2*v(n-1) - v(n-2)`
Temperature :`T(n) = 2*T(n-1) - T(n-2)`
Density : `rho(n) = `2*rho(n-1) - rho(n-2)`
Mach number : `M = v / (sqrtT)`
Pressure : `p = rho *T`
Non dimensional mass flow rate = `rho* a* v`
n -> determines the value of that primitive at nth node.
Non conservation form function program code:
function [rho,T,v,M,ndmf] = non_conservative(n,L,nt)
%inputs
x = linspace(0,L,n);
%distance between the nodes
dx = x(2) - x(1);
%ratio of specific heat
gamma = 1.4;
%Courant-Friedrichs-Levy nummber for stable flow;
cfl = 0.49 ;
%initial conditions at t = 0
%density
rho = 1 - 0.3146*x;
%Temperature
T = 1 -0.2314*x;
%velocity
v = (0.1 + 1.09*x).*T.^0.5;
%Area of cross section along the length of the nozzle
a = 1 + 2.2 * (x - 1.5).^2;
%calculation of 1 time step
for i=1:n
dt(i) = cfl*(dx/((T(i).^0.5)+v(i)));
end
dt = min(dt);
%outer time loop
for k = 1:nt
ts(k) = k;
%Saving the initial values
rho_old = rho;
T_old = T;
v_old = v;
%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;
dlogadx = (log(a(j+1)) - log(a(j)))/dx;
%continuity equation
drhodt_p(j) = -(rho(j) * dvdx) - (rho(j) * v(j) * dlogadx) - (v(j) * drhodx);
%momentum equation
dvdt_p(j) = -(v(j)*dvdx) - ((1/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
v(j) = v(j) + dvdt_p(j)*dt;
T(j) = T(j) + dTdt_p(j)*dt;
rho(j) = rho(j) + drhodt_p(j)*dt;
end
%Corrector step
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;
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 average
drhodt = 0.5 * (drhodt_p + drhodt_c);
dvdt = 0.5 * (dvdt_p + dvdt_c);
dTdt = 0.5 * (dTdt_p + dTdt_c);
%Final solution update
for i = 2:n-1
rho(i) = rho_old(i) + drhodt(i) * dt ;
v(i) = v_old(i) + dvdt(i) * dt ;
T(i) = T_old(i) + dTdt(i) * dt ;
end
%Applying boundary conditions
%inlet
v(1) = 2*v(2) - v(3);
%outlet
%velocity
v(n) = 2*v(n-1) - v(n-2);
%Temperature
T(n) = 2*T(n-1) - T(n-2);
%density
rho(n) = 2*rho(n-1) - rho(n-2);
%Mach number
M = v ./ (T.^0.5);
%Pressure
p = rho .* T;
%Non dimensional mass flow rate
ndmf = rho.* a.* v ;
%Mach number and pressure at the throat
dtM(k) = M(16);
dtp(k) = p(16);
figure(1)
hold on
if k == 50
title ('Non Conservation form');
plot(x,ndmf,'y');
elseif k== 100
plot (x,ndmf,'r')
elseif k==150
plot(x,ndmf,'g')
elseif k==200
plot(x,ndmf,'b')
elseif k==350
plot(x,ndmf,'c')
elseif k==700
plot(x,ndmf,'black')
legend('dt 50','dt100','dt150','dt 200','dt 350','dt 700');
xlabel('Non dimensional distance(x) through nozzle');
ylabel('Non dimensional mass flow');
end
end
%pressure change at a node with respect to time
figure(2)
hold on
plot(ts,dtp);
title('Non conservation form');
xlabel('Number of Time step');
ylabel('Pressure');
%Mach number change at a node with respect to time
figure (3)
hold on
plot (ts,dtM);
title('Non conservation form');
xlabel('Number of Time step');
ylabel('Mach number');
figure(4)
subplot(4,1,1);
plot(x,p,'y');
title ('Non conservation form');
xlabel('Distace(x) along the nozzle');
ylabel('Pressure ratio');
subplot(4,1,2);
plot (x,T,'r');
xlabel('Distace(x) along the nozzle');
ylabel('Temperature ratio');
subplot(4,1,3);
plot(x,rho,'b');
ylabel('Density ratio');
xlabel('Distace(x) along the nozzle');
subplot(4,1,4);
plot(x,M,'g');
ylabel('Mach number');
xlabel('Distace(x) along the nozzle');
end
Plot :
> The plot below shows the mass flow distribution through the nozzle at different times during the time-marching process from the results from the non- conservation form of the governing equation.
>The plot below depicts the variation of primitive variables along the length of the nozzle for isentropic steady state sub sonic to supersonic solution for non-conservation form of governing equation.
It is observed that the mach number increases and is at supersonic speed at the exit of the nozzle, while all other variables such as pressure density and temperature ration decreases along the length of the nozzle.
>Below we can see how the pressure and the mach number in non conservation form of governing equation converges to final solution as the number of iteration increases. In both the cases we can see that the solution is converged at around 500-600 iterations.
Conservation form
The governing equations are converted into suitable partial diffrential equation so that the time marching soluton on quasi 1-D flow can be exectued using these equations.
Governing equations of Conservation form -
- Continuity : `(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
- Momentum : `(del(rho'A'V'))/(delt') + (del(rho'A'V') +(1/gamma)p'A')/(delx') = 1/gamma p' (delA')/(delx')`
- Energy : `(del[rho'((e')/(gamma-1) + gamma/2 V')A'])/(delt') + (del[rho'((e')/(gamma-1) + gamma/2 V')V'A' + p'A'V'])/(delx') =0`
where
`rho` is the density of the flow
A is the area of cross section of the nozzle at that node
V is the velocity of the flow
T is the temperature as the form of energy 'e'
`delt` time step
`gamma` is the ratio of of heat capacity
Non dimensional form of quantities are used in the governing equations since the values of these lies between 0 and 1 which makes it easier to present the result, hence this is preferred in nozzle flow problems.
The non dimensional forms are given by :
`T' = T/(To)`
`rho' = rho/(rhoo)`
`x' = x/L` where x' is the dimensionless component of length L
`V' = V/(ao)` where ao is the speed of sound in the reservoir
`t' = t/(L/(ao))` where t' is the dimensionless time
`A' = A/(Aast)` dimensionless area is given by the ration of area at that point to area at throat of thee nozzle.
Initial conditions
We know that, from start to end of the flow, density and temperature decreases and velocity of the flow increases. This relation is put into equation as shown below and the variation of these quantities are considered to be linear and changes as a function of x . Hence at t=0, we get the flow field variables as :
The area of cross section of the nozzle is such that it converges into a throat and diverges from throat to exit. Hence the area of cross section decreases at first and reaches it's minimum at throat and increases from there. Considering the throat to be at the center of the nozzle along the length of it, we get the equation as :
`a = 1 + 2.2 * (x - 1.5)^2;`
Time step is calculated considering that the system is stable and taking cfl number
`dt = cfl*(dx/(sqrtT+v))`
The time step is calculated for at every node and the minimum of the obtained result is used since all the flow-field variables converge to the same physical time.
governing equations are non dimensional form of the continuity , momentum and energy equations for quasi 1-D flow.
`(del U)/(del t) + (del F)/(del x) + (del G)/(dely) + (del H) / (delz) = J`
The above equation is generic form of governing equation in 3-D flow.
since we are working on 1-D flow we will consider only x component and in order to express governing equations in generic form, let us consider defining in the form of solution vectors U, the flux vector F and the source term J :
`U1 = rho'*A'`
`U2 = rho'*A'*V'`
`U3 = rho'* ((T')/(gamma-1)) + ((gamma*(V')^2)/2))*A'`
`F1 = rho'*A'*V'`
`F2 = rho'*A'(V')^2 + 1/gamma*p'A'`
`F3 = rho' *((T')/(gamma-1) + gamma/2 * (V')^2) V'*A' + p'*A'*V'`
`J2 = 1/gamma * p' * (delA')/(delx')`
After obtaining the respecctive governing equations and setting up the initial conditions, finite-differrence expressions pertaining to MacCormack's method is applied.
The solution is converged for the determined time length and provided total number of time steps, following steps are iterated until the solution is converged
i here defines the node at that time step
>Predictor step
`(dU_1)/dt(i) = -(F1(i+1) - F1(i))/dx`
`(dU_2)/dt(i) = -(F2(i+1) - F2(i))/dx + J2`
`(dU_3)/dt(i) = -(F3(i+1) - F3(i))/dx`
where J2 is given by `(1/gamma) * rho(i) *T(i)*((A(i+1)-A(i))/dx)`
>Corrector step
`(dU_1)/dt(i) =-(F1(i) - F1(i-1))/dx`
`(dU_1=2)/dt(i) =-(F2(i) - F2(i-1))/dx + J2`
`(dU_3)/dt(i) =-(F3(i) - F3(i-1))/dx`
where J2 is given by `(1/gamma) * rho(i) *T(i)*((a(i)-a(i-1)))/dx`
>Averaging :
Averaging is done by taking the average of the result from predictor step and corrector step
>Final step
In this step the averaged result is added with the initial value of the quantity at the start of that iteration
Boundary conditions
At inlet
`U1(1) = a(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^2(1))`
At 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)`
`rho=(U1)/A`
`v=(U2)/(U1)`
`T = (gamma -1)*(U3)/(U1)-(1/2*gamma)*v^2`
Mach number : `M = v / (sqrtT)`
Pressure : `p = rho *T`
Non dimensional mass flow rate = `rho* a* v`
n -> determines the value of that primitive at nth node.
Conservation form function program code:
function [rho,T,v,M,ndmf] = conservative (n,L,nt)
%inputs
x = linspace(0,L,n);
%distance between the nodes
dx = x(2) - x(1);
%ratio of specific heat
gamma = 1.4;
%Courant-Friedrichs-Levy nummber for stable flow;
cfl = 0.5 ;
%initial conditions at t = 0
%Area of cross section along the length of the nozzle
a = 1 + 2.2 * (x - 1.5).^2;
%Density and temperature for dependent variabels U1 U2 and U3 at t=0
for i = 1:n
if x(i) < 0.5
rho(i) = 1;
T(i) = 1;
elseif 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);
elseif 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
v(i)= 0.59/(rho(i) * a(i));
end
%solutions of vector U
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.* ((T./(gamma-1)) + ((gamma.*v.^2)/2)).*a;
dt = zeros(1,length(a));
%calculation of 1 time step
for i=1:n
dt(i) = cfl*(dx/(sqrt(T(i))+v(i)));
end
dt = min(dt);
%outer time loop
for k = 1:nt
ts(k) = k;
%saving the initial values of dependent variables
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1 = U2;
F2 = U2.^2./U1 + ((gamma-1)/gamma).*(U3-(gamma/2)*U2.^2./U1);
F3 = (gamma.*U2.*U3./U1) - (0.5* gamma*(gamma-1).*U2.^3./U1.^2);
%predictor step
for i = 2:n-1
%Source term J
J2 = (1/gamma) * rho(i) *T(i)*((a(i+1)-a(i))/dx);
dU1dtp(i) = -(F1(i+1) - F1(i))/dx;
dU2dtp(i) = -((F2(i+1)-F2(i))/dx)+J2;
dU3dtp(i) = -(F3(i+1) - F3(i))/dx;
%solution update
U1(i) = U1_old(i) + dU1dtp(i)*dt;
U2(i) = U2_old(i) + dU2dtp(i)*dt;
U3(i) = U3_old(i) + dU3dtp(i)*dt;
end
%updating the variables before corrector step
F1 = U2;
F2 = U2.^2./U1 + ((gamma-1)/gamma).*(U3-(gamma/2)*U2.^2./U1);
F3 = (gamma.*U2.*U3./U1) - (0.5* gamma*(gamma-1).*U2.^3./U1.^2);
rho = U1./a;
T = (gamma-1).*(U3./U1-(gamma/2).*(U2./U1).^2);
%corrector step
for i = 2:n-1
%Source term J
J2 = (1/gamma) * rho(i) *T(i)*((a(i)-a(i-1))/dx);
dU1dtc(i) =-(F1(i) - F1(i-1))/dx;
dU2dtc(i) =-(F2(i) - F2(i-1))/dx+J2;
dU3dtc(i) =-(F3(i) - F3(i-1))/dx;
end
%computing average
du1dt = 0.5*(dU1dtp+dU1dtc);
du2dt = 0.5*(dU2dtp+dU2dtc);
du3dt = 0.5*(dU3dtp+dU3dtc);
%solution update
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 ;
end
%Applying boundary conditions
%at inlet
U1(1) = a(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);
rho=U1./a;
v=U2./U1;
T = (gamma -1).*(U3./U1-(gamma*0.5).*v.^2);
M = v./(T.^(0.5));
ndmf = rho.*a.*v;
p = rho .* T;
%Mach number and pressure at the throat
dtM(k) = M(16);
dtp(k) = p(16);
figure(5)
hold on
if k == 50
title ('Conservation form');
plot(x,ndmf,'y');
elseif k== 100
plot (x,ndmf,'r')
elseif k==150
plot(x,ndmf,'g')
elseif k==200
plot(x,ndmf,'b')
elseif k==350
plot(x,ndmf,'c')
elseif k==700
plot(x,ndmf,'black')
legend('dt 50','dt100','dt150','dt 200','dt 350','dt 700');
xlabel('Non dimensional distance(x) through nozzle');
ylabel('Non dimensional mass flow');
end
end
%pressure change at a node with respect to time
figure(6)
hold on
plot(ts,dtp);
title('Conservation form');
xlabel('Number of Time step');
ylabel('Pressure');
%Mach number change at a node with respect to time
figure (7)
hold on
plot (ts,dtM);
title('Conservation form');
xlabel('Number of Time step');
ylabel('Mach number');
figure(8)
subplot(4,1,1);
plot(x,p,'y');
title ('Conservation form');
xlabel('Distace(x) along the nozzle');
ylabel('Pressure ratio');
subplot(4,1,2);
plot (x,T,'r');
xlabel('Distace(x) along the nozzle');
ylabel('Temperature ratio');
subplot(4,1,3);
plot(x,rho,'b');
ylabel('Density ratio');
xlabel('Distace(x) along the nozzle');
subplot(4,1,4);
plot(x,M,'g');
ylabel('Mach number');
xlabel('Distace(x) along the nozzle');
end
Plots of the results:
> The plot below shows the mass flow distribution through the nozzle at different times during the time-marching process from the results from the conservation form of the governing equation.
>The plot below depicts the variation of primitive variables along the length of the nozzle for isentropic steady state sub sonic to supersonic solution for conservation form of governing equation.
It is observed that the mach number increases and is at supersonic speed at the exit of the nozzle, while all other variables such as pressure density and temperature ration decreases along the length of the nozzle.
>Below we can see how the pressure and the mach number in conservation form of governing equation converges to final solution as the number of iteration increases. In both the cases we can see that the solution is converged at around 700-800 iterations.
Comparison
The plot shows the comparison of the steady state mass flow variations obtained with the non conservation and conservation forms of the governing flow equations.
The result obtained from the conservation form is much closer to the exact analytical solution of `rho' A' V' = 0.579`
Grid independence test:
Grid convergence is the term used to describe the improvement of results by using successively smaller cell sizes for the calculations. A calculation should typically approach the correct answer as the mesh becomes finer, hence the term grid convergence.
Here we know that the mach number at the throat of the the nozzle is 1, so as we converge the solution using MacCormack method the solution should be very close to 1 for any number of grid points.
Also we will see how the density changes at the throat by changing the number of grid points.
Here is a code to obtain the values for the pressure, density , temperature and mach number for different number of grid points. The function returns the values of the variables at the throat of the nozzle for grid points of 31,51,101,151,201,251,301 and 1400 time steps.
%main program
clear all
close all
clc
%inputs
%grid points
n = [31,51,71,101,151,201,251,301] ;
%length of the nozzle
L=3;
%total time steps
nt = 1400;
for i = 1: length(n)
%Function call for non conserrvative form
[rhonc(i),pnc(i), Tnc(i), Mnc(i)] = non_conservative(n(i),L,nt);
%Function call for conserrvative form
[rhoc(i),pc(i), Tc(i), Mc(i)] = conservative(n(i),L,nt);
end
figure (1)
plot(n,rhonc);
hold on
plot(n,rhoc);
legend('Non conservative','Conservative')
xlabel('Number of grid points')
ylabel('Pressure ratio at the throat')
figure(2)
plot(n,Mnc);
hold on
plot (n,Mc);
legend('Non conservative','Conservative')
xlabel('Number of grid points')
ylabel('Mach number at the throat')
From the above plots we can see that the result from non conservation form is changing considerbly compared to conservative form. Hence we can conclude that the non conservation form is more dependent on the number of grid points comparitively.
The error with change of gridpoints is less for conservation form for finer mesh.