Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Main object of this project is to solve the 1D super-sonic nozzle flow using numerical method (Maccormack method) for both conservative and non-conservative form.
We consider the steady, isentropic flow through a convergent-divergent nozzle as sketched in the following figure
The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by P0 and T0, respectively. The cross sectional area of the reservoir is large (theoretically, A`rightarrow``infty`), and hence the velocity is very small(V`rightarrow`0). Thus, p0 and T0 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 pe, Te, Ve, and Me, respectively. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat and supersonic at the divergent section.
Mathematical model of the 1D supersonic nozzle flow in non-conservative form:
For the quasi one dimensional subsonic-supersonic isentropic nozzle flow, the non conservation form of the governing equation can be expressed as:
Continuity equation
`Afrac{partialrho}{partialt}+Vfrac{partial(rhoA)}{partialx}=-rhoAfrac{partialV}{partialx}`
Momentum equation
`rhofrac{partialV}{partialt}+rhoVfrac{partialV}{partialx}=-frac{partialp}{partialx}`
Energy equation
`rhofrac{partiale}{partialt}+rhoVfrac{partiale}{partialx}=-pfrac{partialV}{partialx}-pVfrac{partial(lnA)}{partialx}`
Generally, for nozzle flows, the flow field variables are expressed as non-dimensionalize variables.This is done by dividing all the dimensional flow variables by their respective reservoir values.
Therefore, the non-dimensional flow variables are
`Tprime=frac{T}{T_0}, rhoprime=frac{rho}{rho_0}, xprime=frac{x}{L}, a_0=sqrtgammaRT_0, Vprime=frac{V}{a_0}, tprime=frac{t}{frac{L}{a_0}}, Aprime=frac{A}{A^ast}`
Let us now express the non-conservative form in non-dimensional form
Continuity equation(non-dimensionaize)
`frac{partialrhoprime}{partialtprime}=-rhoprimefrac{partialVprime}{partialx}-rhoprimeVprimefrac[partial(lnAprime)}{partialxprime}-Vprimefrac{partialpprime}{partialxprime}`
Momentum equation(non-dimensionalize)
`frac{partialVprime}{partialtprime}=-Vprimefrac{partialVprime}{partialxprime}-(1/gamma)(frac{partialTprime}{partialxprime}+frac{Tprime}{pprime}frac{partialpprime}{partialx}`
Energy equation(non-dimensionalize)
`frac{partialTprime}{partialtprime}=-Vprimefrac{partialTprime}{partialxprime}-(gamma-1)Tprime[frac{partialVprime}{partialxprime}+Vprimefrac{partial(lnAprime)}{partialxprime}]`
Predictor step
First,we consider the predictor step. we set up the spatial derivatives as forward differences. Also, to reduce the complexity of the notation, we will drop the use of the prime to denote a dimensionless variable. in what follows, all variables are the nondimensionalize variables.
now we discretized the above equations by using Maccormack method
`(frac{partialrho}{partialt})_(i)^t=-rho_(i)^tfrac{V_(i+1)^t-V_(i)^t}{Deltax}-rho_(i)^tV_i^tfrac{lnA_(i+1)-lnA_(i)}{Deltax}-V_i^tfrac{rho_(i+1)^t-rho_i^t}{Deltax}`
`(frac{partialV}{partialt})_(i)^t=-V_(i)^tfrac{V_(i+1)^t-V_(i)^t}{Deltax}-1/gamma(frac{T_(i+1)^(t)-T_(i)^(t)}{Deltax}-frac{T_(i)^t}{rho_(i)^t}frac{rho_(i+1)^t-rho_i^t}{Deltax})`
`(frac{partialT}{partialt})_(i)^t=-V_(i)^tfrac{T_(i+1)^t-T_(i)^t}{Deltax}-(gamma-1)T_(i)^t(frac{V_(i+1)^(t)+V_(i)^(t)}{Deltax}+V_(i)^tfrac{lnA_(i+1)-lnA_i}{Deltax})`
We obtain predicted values of `rho, V`and T, denoted by barred quantities, from
`barrho_(i)^(t+Deltat)=rho_(i)^t+(frac{partialrho}{partialt})_(i)^tDeltat`
`barV_(i)^(t+Deltat)=V_(i)^t+(frac{partialV}{partialt})_(i)^tDeltat`
`barT_(i)^(t+Deltat)=T_(i)^t+(frac{partialT}{partialt})_(i)^tDeltat`
in above equation `rho_i^t, V_i^t` and `T_i^t`are known values at time t. Numbers for the time derivative in above Eqs are supplied directly to `(frac{partialrho}{partialt})_(i)^t`, `(frac{partialV}{partialt})_(i)^t``, `(frac{partialT}{partialt})_(i)^t``
Corrector step
In corrector step, we use the spatial derivatives with backward differences, using the predicted(barred) quantities.
`bar(frac{partialrho}{partialt})_(i)^(t+Deltat)=-rho_(i)^(t+Deltat)frac{barV_(i)^(t+Deltat)-barV_(i-1)^(t+Deltat)}{Deltax}-barrho_(i)^(t+Deltat)barV_i^(t+Deltat)frac{lnA_(i)-lnA_(i-1)}{Deltax}-barV_i^(t+Deltat)frac{barrho_(i)^(t+Deltat)-barrho_(i-1)^(t+Deltat)}{Deltax}`
`(frac{partialbarV}{partialt})_(i)^(t+Deltat)=-barV_(i)^(t+Deltat)frac{barV_(i)^(t+Deltat)-barV_(i-1)^(t+Deltat)}{Deltax}-1/gamma(frac{barT_(i)^(t+Deltat)-barT_(i-1)^(t+Deltat)}{Deltax}-frac{barT_(i)^(t+Deltat)}{barrho_(i)^(t+Deltat)}frac{barrho_(i)^(t+Deltat)-barrho_(i-1)^(t+Deltat)}{Deltax})`
`(frac{partialbarT}{partialt})_(i)^(t+Deltat)=-barV_(i)^(t+Deltat)frac{barT_(i)^(t+Deltat)-barT_(i-1)^(t+Deltat)}{Deltax}-(gamma-1)barT_(i)^(t+Deltat)(frac{barV_(i)^(t+Deltat)+barV_(i-1)^(t+Deltat)}{Deltax}+barV_(i)^(t+Deltat)frac{lnA_(i)-lnA_(i-1)}{Deltax})`
The average time derivatives are given by
`(frac{partialrho}{partialt})_(av)=0.5[(frac{partialrho}{partialt})_(i)^(t)+bar(frac{partialrho}{partialt})_(i)^(t+Deltat)]`
`(frac{partialV}{partialt})_(av)=0.5[(frac{partialV}{partialt})_(i)^(t)+(bar(frac{partialV}{partialt}))_(i)^(t+Deltat)]`
`(frac{partialT}{partialt})_(av)=0.5[(frac{partialT}{partialt})_(i)^(t)+(bar(frac{partialT}{partialt}))_(i)^(t+Deltat)]`
we have for the corrected values of the flow-field variables at time `t+Deltat`
`rho_(i)^(t+Deltat)=rho_(i)^t+(frac{partialrho}{partialt})_(av)Deltat`
`V_(i)^(t+Deltat)=V_(i)^t+(frac{partialV}{partialt})_(av)Deltat`
`T_(i)^(t+Deltat)=T_(i)^t+(frac{partialT}{partialt})_(av)Deltat`
Calculation of time step
In this case, we used the time step value
`Deltat=Cfrac{Deltax}{sqrtT+V}`
Nozzle shape and initial conditions
The nozzle shap, A=A(x), is specified and held constant, independent of time. we choose the parabolic area distribution given by
`A=1+2.2(x-1.5)^2, 0lexle3`
The initial conditions must be chosen as close to the steady state as possible through the knowledge and experience. otherwise solution is going to unstable. in this problem rho and T decreases and V increases as the flow expands tgrough the nozzle. For simplicity, let us assume linear variations of the flow field variables, as a function of x. we assume the following values at time t=0.
`rho=1-0.3146x`
`T=1-0.2314x`
`V=(0.1+1.09x)T^(1/2)`
Boundary Conditions
Inlet:
At the subsonic inlet, we keep rho and T fixed and allow V to float, we calculate this value by using linear extrapolation
V(1)=2*V(2)-V(3)
oulet
At the supersoinc outlet, all the three flow field variables are allowed to float
`V_N=2V_(N-1)-V_(N-2)`
`rho_N=2rho_(N-1)-rho_(N-2)`
`T_N=2T_(N-1)-T_(N-2)`
Mathematical model of the 1D supersonic nozzle flow in conservative form:
For the quasi one dimensional subsonic-supersonic isentropic nozzle flow, the conservation form of governing equation can be expressed as;
Continuity equation
`frac{partial( rhoA)}{partialt}+frac{partial( rhoAV)}{partialx}=0`
Momentum equation
`frac{partial( rhoAV)}{partialt}+frac{partial( rhoAV^2)}{partialx}=-Afrac{partialp}{partialx}`
Energy equation
`frac{partial[rho(e+V^2/2)A]}{partialt}+frac{[rhoAV^2+(1/gamma)pA]}{partialx}=(1/gamma)pfrac{partialA}{partialx}`
Expressing the above conservation form equation in non-dimensional variables, we get
Continuity equation
`frac{partial( rhoprimeAprime)}{partialt}+frac{partial( rhoprimeAprimeVprime)}{partialx}=0`
Momentum equation
`frac{partial( rhoprimeAprimeVprime)}{partialtprime}+frac{partial( rhoprimeAprimeVprime^2)+(1/gamma)pprimeAprime}{partialxprime}=1/gammapprimefrac{partialAprime}{partialxprime}`
Energy equation
`frac{partial[rhoprime(frac{eprime}{gamma-1}+frac{gamma}{2})Aprime]}{partialtprime}+frac{partial[rhoprime(frac{eprime}{gamma-1}+frac{gamma}{2}Vprime^2)VprimeAprime+pprimeAprimeVprime]}{partialxprime}=0`
Let us define the elements of the solutions vector U, the flux vector F and the source term J as follows
`U1=rhoprimeAprime`
`U2=rhoprimeAprimeVprime`
`U3=rhoprime(frac{eprime}{gamma-1}+frac{gamma}{2}Vprime^2)Aprime`
`F1=rhoprimeAprimeVprime`
`F2=rhoprimeAprimeVprime^2+frac{1}{gamma}pprimeAprime`
`F3=rhoprime(frac{eprime}{gamma-1}+frac{gamma}{2}Vprime^2)VprimeAprime+pprimeAprimeVprime`
`J_2=frac{1}{gamma}pprimefrac{partialAprime}{partialxprime}`
put the value of U1,U2,U3,F1,F2,F3 and J2 into the non-dimensionalize form of continuity,momentum and energy equation, we obtain the following forms
Continuity equation
`frac{partialU_1}{partialt} = -frac{partialF_1}{partialx}`
Momentum equation
`frac{partialU_2}{partialt} = -frac{partialF_2}{partialx}+J_2`
Energy equation
`frac{partialU_3}{partialt} = -frac{partialF_3}{partialx}`
while writing the code we will solve for the solution vector and then immediately calculate the primitive variables using the following relations
`rhoprime=frac{U_1}{Aprime}`
`Vprime=frac{U_2}{U1}`
`Tprime=eprime=(gamma-1)(frac{U_3}{U_1}-frac{gamma}{2}Vprime^2)`
`pprime = rhoprimeTprime`
it is the best choice to express the F1,F2 and F3 in terms of the solution vector U1,U2 and U3 and avoid the use of the primitive variables, which is given below
`F_1=U_2`
`F_2=frac{U_2^2}{U1}+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}`
`J_2=frac{gamma-1}{gamma}(U_3-frac{gamma}{2}frac{U_2^2}{U_1})frac{partial(lnAprime)}{partialxprime}`
Discretizated the conservation form of governing equation using maccormack method.
Predictor Method
First consider the predictor step.We set up the spatial derivatives as forward differences.
`(frac{partialU_1}{partialtprime})_(i )^t= -frac{F_(1(i+1))^t-F_(1(i))^t}{Deltaxprime}`
`(frac{partialU_2}{partialtprime})_(i )^t= -frac{J_(2(i+1))^t-J_(2(i))^t}{Deltaxprime} -frac{F_(2(i+1))^t-F_(2(i))^t}{Deltaxprime}`
`(frac{partialU_3}{partialtprime})_(i )^t= -frac{F_(3(i+1))^t-F_(3(i))^t}{Deltaxprime`
These time derivative values are used to calculate the values of U1,U2 and U3 at `t+Deltat`time step.
`barU_(1(i))^(t+Deltat)=(frac{partialU_1}{partialtprime})_(i )^tDelta t`
`barU_(2(i))^(t+Deltat)=(frac{partialU_2}{partialtprime})_(i )^tDelta t`
`barU_(3(i))^(t+Deltat)=(frac{partialU_3}{partialtprime})_(i )^tDelta t`
Corrector step
In Corrector step, we setup the spatial derivative as a backward difference
`(barfrac{partialU_1}{partialtprime})_(i )^(t+Deltat)= -frac{barF_(1(i))^(t+Deltat)-barF_(1(i-1))^(t+Deltat}}{Deltaxprime}`
`(barfrac{partialU_2}{partialtprime})_(i )^(t+Deltat)= -frac{barJ_(2(i))^(t+Deltat)-barJ_(2(i-1))^(t+Deltat)}{Deltaxprime} -frac{barF_(2(i))^(t+Deltat)-barF_(2(i-1))^(t+Deltat)}{Deltaxprime}`
`(barfrac{partialU_3}{partialtprime})_(i )^(t+Deltat)= -frac{barF_(3(i+1))^(t+Deltat)-barF_(3(i))^(t+Deltat}}{Deltaxprime}`
The average value of the predictor and corrector step derivatives is calculated
`(frac{partialU_1}{partialtprime})_(av)=frac{1}{2}[(frac{partialU_1}{partialtprime})_(i )^t+(barfrac{partialU_1}{partialtprime})_(i )^(t+Deltat)]`
`(frac{partialU_2}{partialtprime})_(av)=frac{1}{2}[(frac{partialU_2}{partialtprime})_(i )^t+(barfrac{partialU_2}{partialtprime})_(i )^(t+Deltat)`
`(frac{partialU_3}{partialtprime})_(av)=frac{1}{2}[(frac{partialU_3}{partialtprime})_(i )^t+(barfrac{partialU_3}{partialtprime})_(i )^(t+Deltat)]`
The corrected values of U1,U2 and U3 are computed using the average of derivatives of values
`U_(1(i))^(t+Deltat)=(frac{partialU_1}{partialtprime})_(av)Delta tprime`
`U_(2(i))^(t+Deltat)=(frac{partialU_2}{partialtprime})_(av)Delta tprime`
`U_(3(i))^(t+Deltat)=(frac{partialU_3}{partialtprime})_(av)Delta tprime`
Nozzle Shape and initial conditions
The nozzle shape and area is the same as mention in above. The area as a function of x is given as
`A=1+2.2(x-1.5)^2`
The initial conditions for U1,U2 and U3 were synthesized by assuming the following variations of `rhoprime`and `Tprime`:
`for 0 lexprimele0.5{rhoprime = 1, Tprime=1}`
`for 0.5 lexprimele1.5{rhoprime = 1.0-0.366(xprime-0.5), Tprime=1.0-0.167(xprime-0.5)}`
`for 1.5 lexprimele3.5{rhoprime = 0.634-0.3879(xprime-0.5), Tprime=0.833-0.3507(xprime-0.5)}`
The initial conditions for the variation of `Vprime`is synthesized by taking the constant mass flow rate through the nozzle and calculated `Vprime`as
`Vprime=frac{U_2}{rhoprimeAprime}=frac{0.59}{rhoprimeAprime}`
But the dependent variables being solved by the governing equatios are U1,U2,U3. We need initial conditions for them. The same initial conditions can be applied to U1,U2 and U3 using the relations
`U_1=rhoprimeAprime`
`U_2=rhoprimeVprimeAprime=0.59`
`U_3=rhoprime(frac{eprime}{gamma-1}+frac{gamma}{2}Vprime^2)Aprime`
Boundary Conditions
At the subsonic inflow boundary two properties are held fixed and one is allowed to float and at the supersonic outflow boundary all properties are allowed to float.
Inlet:
We hold `rhoprime`and `Tprime`fixed at the inflow boundary both equal to 1.0, and allow`Vprime` to float. By holding `rhoprime`fixed, then `U_1`at grid point i=1 fixed, independent of time, via `U_1=rhoprimeAprime`.that is,
```U_(1(i=1))=(rhoprimeAprime)_(i=1)=Aprime_(i=1)=`fixed value.
The floating value of `Vprime`at the inflow boundary is calculated at the end of each time step by linearly extrapolating `U_2`from the known values at the internal grid points i=2 and 3, that is,
`U_(2(i=1))=2U_(2(i=2))-U_(2(i=3))`
and then obtaining `Vprime=frac{U_2}{U_1}`at i=1. Since `Vprime`floats at the inflow boundary, so does the value of `U_3`, which is given by
`U_3=U_1(frac{Tprime}{gamma-1}+frac{gamma}{2}Vprime^2)`
Outlet:
The flow properties at the downstream, supersonic outflow boundary are obtained by linear extrapolation from the two adjacent internal points. If N denotes the grid point at the outflow boundary, then
`(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)`
Matlab Codes
% Simulation of 1D supersoic nozzle flow using Macmormack method
clear all
close all
clc
%global V,rho,T,m,M,P,Vc,rhoc,Tc,mc,Mc,Pc,n,nt
n=31; % grid points
x = linspace(0,3,n);
dx = x(2)-x(1);
%area
A = 1+2.2*(x-1.5).^2;
gamma = 1.4;
%Time step
nt=1200;
C = 0.5; %courant number
throat = ((n-1)/2)+1; %discuss about this
% Non-conservative form function
tic
[V,rho,T,m,M,P,V_t,rho_t,T_t,m_t,M_t,P_t]= non_consv(n,nt,throat);
non_consv_form = toc;
% conservative form function
tic
[Vc,rhoc,Tc,mc,Mc,Pc,V_tc,rho_tc,T_tc,m_tc,M_tc,P_tc]= consv_func_updated(n,nt,throat);
consv_form = toc;
% plots for conservative vs non-conservative form
% Mass flow rate
figure(3)
plot(x,m,\'color\',\'c\',\'markersize\',10)
hold on
plot(x,mc,\'color\',\'r\',\'linewidth\',1)
% %Exact analytical solution rho*A*V=0.579
xlabel(\'length of nozzle(x)\');
ylabel(\'mass flow rate(m)\');
text1 = sprintf(\'conservation analysis time=%d\',non_consv_form);
text2 = sprintf(\'non-conservation analysis time=%d\',consv_form);
title({\'comparison of mass flow rate consv vs non consv\';text1;text2});
legend(\'Non-conservative form\',\'conservative form\',\'location\',\'northwest\');
grid on
%plotting different parameters like density, velocity,temperature,mach
%number and pressure.variation in non-dimensional flow parameter along x-axis
%Mach no
figure(4)
subplot(5,1,1)
plot(x,M,\'r\',\'linewidth\',1.5);
hold on
plot(x,Mc,\'*\',\'color\',\'r\',\'markersize\',10);
grid on;
xlabel(\'length(x)\');
ylabel(\'Mach Number\')
legend(\'Non conservative Form\',\'Conservative Form\',\'location\',\'northwest\');
title(\'Mach Number profile\')
% %Density
subplot(5,1,2)
plot(x,rho,\'r\',\'linewidth\',1.5);
hold on
plot(x,rhoc,\'*\',\'color\',\'r\',\'markersize\',10);
grid on;
xlabel(\'length(x)\');
ylabel(\'Density\')
legend(\'Non conservative Form\', \'Conservative Form\',\'location\',\'northwest\');
title(\'change in density in x-axis\')
% %Pressure
subplot(5,1,3)
plot(x,P,\'r\',\'linewidth\',1.5);
hold on
plot(x,Pc,\'*\',\'color\',\'r\',\'markersize\',10);
grid on;
xlabel(\'length(x)\');
ylabel(\'Pressure\')
legend(\'Non conservative Form\', \'Conservative Form\',\'location\',\'northwest\');
title(\'Pressure profile\')
% %Temperature
subplot(5,1,4)
plot(x,T,\'r\',\'linewidth\',1.5);
hold on
plot(x,Tc,\'*\',\'color\',\'r\',\'markersize\',10);
grid on;
xlabel(\'length(x)\');
ylabel(\'Temperature\')
legend(\'Non conservative Form\', \'Conservative Form\',\'location\',\'northwest\');
title(\'change in Temperature\')
% %Velocity
subplot(5,1,5)
plot(x,V,\'color\',\'g\',\'linewidth\',1.5);
hold on
plot(x,Vc,\'*\',\'color\',\'r\',\'markersize\',10);
grid on;
xlabel(\'length(x)\');
ylabel(\'Velocity\')
legend(\'Non conservative Form\', \'Conservative Form\',\'location\',\'northwest\');
title(\'study of velocity along x-axis\')
%Variation of flow properties at throat
figure(5)
subplot(4,1,1);
plot(1:nt,rho_tc);
hold on
plot(1:nt,rho_t,\'color\',\'r\');
grid on;
ylabel(\'Density\');
legend(\'conservative Form\', \'Non_Conservative Form\');
title(\'study of density variation at each time step at throat\');
subplot(4,1,2);
plot(1:nt,T_tc);
hold on
plot(1:nt,T_t);
grid on
ylabel(\'Temperature\')
legend(\'conservative Form\', \'Non_Conservative Form\');
title(\'study of Temperature variation at each time step at throat\');
subplot(4,1,3);
plot(1:nt,P_tc);
hold on
plot(1:nt,P_t);
grid on
ylabel(\'Pressure\')
legend(\'conservative Form\', \'Non_Conservative Form\');
title(\'study of Pressure variation at each time step at throat\');
subplot(4,1,4);
plot(1:nt,M_tc);
hold on
plot(1:nt,M_t);
grid on
xlabel(\'Time step\')
ylabel(\'Mach number\')
legend(\'conservative Form\', \'Non_Conservative Form\');
title(\'study of Mach no variation at each time step at throat\');Non_Conservative_code
%Creating mesh and assigning points
%global V,rho,T,m,M,P,Vc,rhoc,Tc,mc,Mc,Pc,n,nt
n=31; % grid points
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
C = 0.5;
%initial profiles
rho = 1-0.3146*x;
T = 1-0.2314*x;
V = (0.1+1.09*x).*T.^0.5;
A = 1+2.2*(x-1.5).^2;
%Time step calculated using courant number for optimal value
dt = min(C.*dx./(sqrt(T)+V));
%outer time loop
for k=1:nt
% copying old values
rho_old = rho;
V_old = V;
T_old = T;
%predictor method
for i = 2:n-1 %space loop
%drho_dt_p= -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;
dvdx = (V(i+1)-V(i))/dx;
drhodx= (rho(i+1)-rho(i))/dx;
dtdx = (T(i+1)-T(i))/dx;
dlogadx = (log(A(i+1))-log(A(i)))/dx;
%continuity equation
drhodt_p(i)= -rho(i)*dvdx-rho(i)*V(i)*dlogadx-V(i)*drhodx;
%Momentum equation
dvdt_p(i)= -V(i)*dvdx-(1/gamma)*(dtdx+(T(i)/rho(i))*drhodx);
%Energy equation
dtdt_p(i) = -V(i)*dtdx-(gamma-1)*T(i)*(dvdx+V(i)*dlogadx);
%solution update
rho(i)= rho(i)+drhodt_p(i)*dt;
V(i) = V(i)+dvdt_p(i)*dt;
T(i) = T(i)+dtdt_p(i)*dt;
end
%Corrector method
for i = 2:n-1
%drho_dt_p= -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;
dvdx = (V(i)-V(i-1))/dx;
drhodx= (rho(i)-rho(i-1))/dx;
dtdx = (T(i)-T(i-1))/dx;
dlogadx = (log(A(i))-log(A(i-1)))/dx;
%continuity equation
drhodt_c(i)= -rho(i)*dvdx-rho(i)*V(i)*dlogadx-V(i)*drhodx;
%Momentum equation
dvdt_c(i)= -V(i)*dvdx-(1/gamma)*(dtdx+(T(i)/rho(i))*drhodx);
%Energy equation
dtdt_c(i) = -V(i)*dtdx-(gamma-1)*T(i)*(dvdx+V(i)*dlogadx);
end
%compute the average time derivative
drhodt = 0.5*(drhodt_p+drhodt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dtdt = 0.5*(dtdt_p+dtdt_c);
%Final solution update
for i = 2:n-1
rho(i)= rho_old(i)+drhodt(i)*dt;
V(i) = V_old(i)+dvdt(i)*dt;
T(i) = T_old(i)+dtdt(i)*dt;
end
%Applying the boundary conditions
%inlet
V(1) = 2*V(2)-V(3);
%outlet
V(n) = 2*V(n-1)-V(n-2);
T(n) = 2*T(n-1)-T(n-2);
rho(n) = 2*rho(n-1)-rho(n-2);
%Non-dimens pressure
P = rho.*T;
%Mach number
M = V./sqrt(T);
%Mass flow rate
m = rho.*V.*A;
% %throat values
rho_t(k)= rho(throat);
T_t(k) = T(throat);
P_t(k) = P(throat);
M_t(k) = M(throat);
V_t(k) = V(throat);
m_t(k) = m(throat);
%plotting mass flow rate at different time step
figure(1)
if k==1
plot(x,m,\'*\',\'color\',\'b\',\'markersize\',10)
hold on
grid on
elseif k==50
plot(x,m,\'*\',\'color\',\'y\',\'linewidth\',1.5)
elseif k==100
plot(x,m,\'*\',\'color\',\'c\',\'linewidth\',1.5)
elseif k==200
plot(x,m,\'*\',\'color\',\'g\',\'linewidth\',1.5)
elseif k==600
plot(x,m,\'*\',\'color\',\'r\',\'linewidth\',1.5)
elseif k==800
plot(x,m,\'*\',\'color\',\'k\',\'linewidth\',1.5)
elseif k==1000
plot(x,m,\'--\',\'color\',\'bla\',\'linewidth\',1.5)
hold on
%labeling the plot
xlabel(\'Length of Nozzle(x)\');
ylabel(\'Mass flow rate(m)\');
title({\'Non dimensional mass flow rate at different Time steps\';\'(Non conservative)\'});
legend({\'1 Time steps\',\'50 Time step\',\'100 Time steps\',\'200 Time steps\',\'600 Time steps\',\'800 Time steps\',\'1000 Time steps\'});
end
end
endConservative_code
%Function for solving conservative form
function[Vc,rhoc,Tc,mc,Mc,Pc,V_tc,rho_tc,T_tc,m_tc,M_tc,P_tc] = consv_func_updated(n,nt,throat)
%global V,rho,T,m,M,P,Vc,rhoc,Tc,mc,Mc,Pc,n,nt
%Creating mesh and assign inputs
n=31; % grid points
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
C=0.5;
%Calculating the initial profiles
for i=1:n
if (x(i)>=0 && x(i)<=0.5)
rhoc(i) = 1;
Tc(i) = 1;
elseif (x(i)>=0.5 && x(i)<=1.5)
rhoc(i) = 1-0.366*(x(i)-0.5);
Tc(i) = 1-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3.5)
rhoc(i) = 0.634-0.3879*(x(i)-1.5);
Tc(i) = 0.833-0.3507*(x(i)-1.5);
end
end
% Calculating the nondimensional numbers
A = 1+2.2*(x-1.5).^2;
Vc = 0.59./(rhoc.*A); %velocity
%Pc = rhoc.*Tc;
%Solution vector parameters
U1 = rhoc.*A;
U2 = rhoc.*A.*Vc;
U3 = rhoc.*A.*((Tc./(gamma-1))+((gamma/2)*(Vc.^2)));
%outer loop
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
for k = 1:nt
%Time step
dt = min(C.*dx./(sqrt(Tc)+Vc));
%copying values to old
U1_old = U1;
U2_old = U2;
U3_old = U3;
%Flux terms
% F1 = U2;
%F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
%F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
for i= 2:n-1
F1(i) = U2(i);
F2(i) = ((U2(i).^2)./U1(i)) + (((gamma-1)/gamma)*(U3(i) - ((gamma*0.5*(U2(i).^2))./U1(i))));
F3(i) = ((gamma*U2(i).*U3(i))./U1(i)) - ((