Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective: For a convergent and divergent nozzle
1. To derive non-conservative and conservative forms of continuity, momentum and energy equations.
2. To show solution technique using MacCormack's method (Discretization form).
3. Comparing conservative and non-conservative governing equations based on the final steady state results of Temperature,velocity and density as a function of space and time.
Literature survey:
The generalised continuity, momentum and energy equations (i.e navier stokes equations) in differential form are independent of area variation in space and time coordinates, as we consider infinitesimal small cubical structure where area does not change with respect to direction. In the present case, we are considering a convergent and divergent nozzle case where area does change with respect to space coordinates. This indicates that, the generalised equations which are developed previously cannot be used directly.
Figure-1 indicates control volume that is considered used for developing differential equation for generalised conditions.
Figure-2 indicates the case for which governing equations to be developed (convergent and divergent nozzle)
Sequentional steps towards governing equations:
Starting with Reynold's Transport Theorem:
`(DB_(sys))/(Dt)=(\partialB_(cv))/(\partialt)+oint\rho(vec(v).hat(n))bdA`--------(1)
- B-indicates extensive property. (this can be any extensive property namely mass, velocity and temperature)
- b- extensive property per unit mass. Also called intesive property of the system.
- LHS indicates extensive property rate of change within a system.
- RHS 1st term indicates extensive property change within the control volume. By definition this can be written as `(\partial oint \rhodV.b)/(\partialt)`
- RHS 2nd term indicates equivalent amount of property either entering or leaving the control volume.
- `v`indicates velocity
- `V`indicates volume
- `hatn`normal to the area
Finally equation-1 becomes:
`(DB_(sys))/(Dt)=(\partial oint \rhodV.b)/(\partialt)+oint\rho(vec(v).hat(n))bdA`--------(2)
Now lets take up our control volume where the governing equation must be applied as shown in below Figure (3).
Continuity equation:
Applying mass balance on this control volume,
indicates 'B=mass', thus indicates `b=m/m=1`.
Substituting this into reynolds transport theorem.
LHS indicates rate of change of mass of a system. As we know, in a system mass can neither be created nor destroyed. Thus, LHS becomes zero.
Finally, the equation jot down's to
`(\partial oint \rhodV)/(\partialt)+oint\rho(v.hatn)dA=0`.........(3)`
- `dV = 0.5*(A+A+dA) dx`(This is trapezoidal volume). After considering `dAdx` in the expansion, that leads to `dV=Adx`.
- The second term also indicates that, total mass that is entering/exiting from all the directions. Thus the integral can be replaced and written as, `-\rhovA+(\rho+d\rho)(v+dv)(A+dA)`. Further after little bit of mathematical calculation considering few multiples of infinitesimal terms to zero. That leads to,`d(\rhoVA)`
Thus, the resultant continuity equation converts into
`(d(\rhoAdx))/(dt)+d(\rhoAv)=0`......(4)
As dx element is a constant, thus equation 4 can be written as
`(d(\rhoA))/(dt)+(d(\rhoAv))/dx=0`........(5)`
Equation-5, indicates conservative form of continuity equation.
On expanding Equation-(5), that leads to non-conservative form of continuity equation
`(d(\rhoA))/dt+\rhoA(dv)/(dx)+\rhov(dA)/(dx)+vA(d\rho)/(dx)=0`...(6)
Equation-6 indicates non-conservative form of continuity equation.
Generally, for nozzle type problems many of the computational engineers prefer to write the governing equations in non-dimensional forms. So, the above governing equation for mass in non-dimensional form is written as in equation-(7),(8)
During non-dimensionalisation following dimensional terms are replaced by non-dimensional terms,
`T^'=T/T_0,\rho^'=\rho/\rho_0,x^'=x/L,a_0=sqrt(\gammaRT_0),v^'=v/a_0,t^'=ta_0/L,A^'=A/A^*,e^'=e/e_0,e_0=c_vT_0=RT_0/(\gamma-1)`
- `T_0 text( indicates temperature of the stagnant fluid)`
- `\rho_0 text( indicate density of the stagnant fluid)`
- `a_0 text( indicates acoustic wave speed in the stagnant fluid)`
- `L text( indicates Length of the nozzle)`
- `A^* text( indicates cross section area of the throat in the nozzle)`
- `e_0 text ( indicates internal energy at stagnant temperature)`
- `R text ( indicates gas constant)`
Substituting the above terms into equations (5) and (6) leads to equations(7) and (8) respectively. For my convenience each element after non-dimensionalisation is not represented in prime form in `T^' rather they are represented as `T`.
Non-dimensional conservative form of continuity equation,
`(\partial\rhoA)/(\partialt)+(\partial\rhoAv)/(\partialx)=0` ......(7)
Non-Dimensional non-conservative form of continuity equation,
`(\partial\rho)/(\partialt)=-\rho(\partialv)/(\partialx)-\rhov(\partiallnA)/(\partialx)-v(\partial\rho)/(\partialx)`........(8)
Momentum Equation in x-direction only:
As similar to continuity equation,
In the reynolds transport theorem (equation-2), if `B=momentum=mv, text(then ) b=(mv)/m=v` (here v indicates velocity in x-direction only)
Thus, after substituting 'B' and 'b' in equation-2.
- L.H.S represents rate of change of momentum of a system by lagrangian approach indicates summation of forces on a system. In our case, lets consider there are only pressure forces for our convenience.
- R.H.S 1st term reprsents rate of change of momentum within the control volume.
- R.H.S 2nd term indicates equivalent momentum that is either entering or leaving the system.
Finally, Momentum equation jot's down to
`\SigmaF=(\partial oint \rhovdV)/(\partialt)+oint\rhov(v.hatn)dA`.....(9)
- As previously done, integral on R.H.S 1st is removed and replaced `dV text( with ) Adx`
- R.H.S second term, momentum entering or leaving the control volume is written as follows, `oint\rhov(v.hatn)dA=-\rhov^2A+(\rho+d\rho)(v+dv)^2(A+dA)`.
This on further expansion, after considering few terms that are negligible ends up to `oint\rhov(v.hatn)dA=d(\rhov^2A)`
- As LHS in equation-9 represents summation of forces, it is represented in the form of Figure (4) as follows,
Summation of all the forces acting on the infinitesimal surface in the equation form is,`\SigmaF=pA+2*p((dA)/(2))-(p+dp)(A+dA)`. Which on further expansion leads to `\SigmaF=-Adp`.
Thus, now substituting all these terms into equation-9 leads to,
`(\partial(\rhovA))/(\partialt)+(\partial(\rhov^2A))/(\partialx)=-A(\partialp)/(\partialx)`........(10)
Equation-10 indicates conservative form of momentum equation.
Upon expansion of terms on LHS of equation-10 leads as follows,
`v*((\partial(\rhoA))/(\partialt)+(\partial(\rhoAv))/(\partialx))+\rhoA*((\partialv)/(\partialt)+v(\partialv)/(\partialx))=-A(\partialp)/(\partialx)`.....(11)
The summation of terms (i.e with in paranthesis) which are in product with velocity on LHS of equation indicate, continuity equation in coservative form which is equal to zero. Thus, equation-11 compresses down to
`\rhoA*((\partialv)/(\partialt)+v(\partialv)/(\partialx))=-A(\partialp)/(\partialx)`.......(12)
Equation-12 indicates momentum equation in x-direction in non-conservative form.
Substituting the non-dimensional variables in equation-10 and 12 leads to non-dimensional conservative and non-conservative forms of momentum equations in x-direction respectively. (For my convenience each element after non-dimensionalisation is not represented in prime form in T^' rather they are represented as `T`).
`(\partial(\rhovA))/(\partialt)+(\partial(\rhov^2A+(1/gamma)pA))/(\partialx)=1/(\gamma)p(\partialA)/(\partialx)`.......(13)
`(\partial v)/(\partialt)=-v(\partial v)/(\partialx)-(1/gamma)((\partialT)/(\partialx)+(T/\rho)((\partial \rho)/(\partial x)))`......(14)
Energy equation :
As similar to momentum equation derivation, starting with reynolds transport equation.
B=m(e+v^2/2)=internal energy+kinetic energy and b=B/m=e+v^2/2.
- In the final development of energy equation, LHS remains the same as equation-10 by just replacing `v` with `e+v^2/2`. (These steps are skipped, as the process is as similar to momentum equation)
- LHS of equation-2 in reynolds transport theorem indicates rate of change of energy in a system. It indicates summation of rate of work done on system and rate of conduction heat transfer. As there is only pressure work, `\Sigma dotW= pvA-(p+dp)(A+dA)(v+dv)`. Upon expansion leaving some negligible terms leads to `\Sigma dotW=-d(pAv)`.
Sustituting all the above points into equation form leads to,
`(\partial(\rho(e+v^2/2)A))/(\partialt)+(\partial(\rho(e+v^2/2)vA))/(\partialx)=-(\partial pAv)/(\partialx)`........(13)
Equation-13 is conservative form of energy equation.
multiplying equation-10 with `'v'`and substracting that equation in 13 leads to,
`(\partial(\rhoeA))/(\partialt)+(\partial(\rhoevA))/(\partialx)=-p(\partial Av)/(\partialx)`......(14)
Equation-14 indicates conservative form of internal energy.
Equation-14 on further expansion on LHS leads to,
`e((\partial(\rhoA))/(\partialt)+(\partial(\rhovA))/(\partialx))+\rhoA((\partiale)/(\partialt)+v(\partiale)/(\partialx))=-p(\partial Av)/(\partialx)`......(15)
The summation of terms (i.e with in paranthesis) which are in product with velocity on LHS of equation indicate, continuity equation in coservative form which is equal to zero. Thus, equation-15 compresses down to
`\rhoA((\partiale)/(\partialt)+v(\partiale)/(\partialx))=-pA(\partial v)/(\partialx)-pv(\partial A)/(\partialx)`.....(16)
Equation-16 indicates non-conservative of form of energy equation.
Writing internal energy and pressure in terms of temperature `e=c_vT` and `p=\rhoRT` ends up equation-16 to,
`\rhoc_v((\partialT)/(\partialt)+v(\partialT)/(\partialx))=-\rhoRT((\partial v)/(\partialx)-(v/A)(\partial A)/(\partialx))`........(17)
Upon non-dimensionalising the terms, equation-14 and 17. They are converted into conservative and non-conservative form of energy equation.(For my convenience each element after non-dimensionalisation is not represented in prime form in T^' rather they are represented as `T`).
Conservative form of energy equation in non-dimensional form,
`(\partial(\rho((e)/(gamma-1)+v^2/2)A))/(\partialt)+(\partial(\rho((e)/(gamma-1)+v^2/2)vA))/(\partialx)=-(\partial pAv)/(\partialx)`.....(18)
Non-conservative form of energy equation in non-dimensional form,
`(\partialT)/(\partialt)=-v(\partialT)/(\partialx)-(\gamma-1)T((\partial v)/(\partialx)-(v/A)(\partial A)/(\partialx))`..........(19)
Application of Macormack Method:
This is one such FDM explicit method used to solve transport equations to get an error of second order accurate in both space and time. The application procedure is explained for a very simple transport equation of the form shown in the equation-(20-22) for easy understanding.
The below equations are non-conservative form of continuity,momentum and energy equation in one-dimensional space,
`(\partial\rho)/(\partialt)=-(\rho(\partialu)/(\partialx)+u(\partial\rho)/(\partialx))`...................(20)
`(\partialu)/(\partialt)=-(u(\partialu)/(\partialx)+(1/\rho)(\partialp)/(\partialx))`......................(21)
`(\partiale)/(\partialt)=-((p/\rho)(\partialu)/(\partialx)+u(\partiale)/(\partialx))`......................(22)
Lets consider our nodal representation of space coordinates as shown in below form of Figure (5) and which is marching forward in time.
Procedure that is involved in estimating a variable in one time step is as shown below:
Predictor Step: All spatial coordinates on RHS are discretized forward in space and time differentials are calculated. As an example time derivative of density is calculated from continuity equation below.
`(((\partial\rho)/(\partialt))_i^{t))_p=-(\rho_i^t((u_(i+1)^t-u_(i)^t)/(\Deltax))+u_i^t((\rho_(i+1)^t-\rho_(i)^t)/(\Deltax))` ........(23)
The same procedure must be applied to all time dervatives (i.e for equation 21 and 22).
As all the flow variables are known on RHS of equation 23. Thus, partial derivative of density with respect to time is known.
Now predicted values of all the flow variables are estimated:
`(\rho_i^t)_p=\rho_i^t+(((\partial\rho)/(\partialt))_i^{t))_p\Deltat`.............(24)
As similar to density all the flow variables predicted values are estimated.
Corrector Step: Now as similar to predictor step, spatial coordinates must be discretized but now backward in space.
`(((\partial\rho)/(\partialt))_i^{t))_c=-(\rho_i^t((u_(i)^t-u_(i-1)^t)/(\Deltax))+u_i^t((\rho_(i)^t-\rho_(i-1)^t)/(\Deltax))`........(25)
Now average of time slopes are estimated and that are used for final estimation of flow variables in that particular time step.
`(((\partial\rho)/(\partialt))_i^{t))_a=0.5((((\partial\rho)/(\partialt))_i^{t))_p+(((\partial\rho)/(\partialt))_i^{t))_c)`.......(26)
`\rho_i^(t+\Deltat)=\rho_i^t+(((\partial\rho)/(\partialt))_i^{t))_a\Deltat`..................(27)
Thus, equation-27 ends up giving density value for the next time step. A similar procedure must be employed (i.e from equations 23-27) for all the flow variables to calculate the next time step value and the programme continues to run in a loop till the error converges to minimum value.
Applying Macormack method to convergent and divergent nozzle equation:
The continuity,momentum and energy equations are developed in conservative and non-conservative forms previously for convergent and divergent nozzle. Applying the Macormack method is just the same as discussed in the above section for non-conservative form of equations (i.e estimating the flow variables). But where as to estimate flow variables in a conservative form it is not just the same as the non-conservative form. It is expected that, all the variables which are in time derivative are considered to be of a single cauch (lets say it as U) and the remaining dependent variables as a function of spatial coordinates are also represented in the same new represented variable (i.e in terms of U).
For a better understanding, lets represent continuity, momentum and energy equations in terms of new variables so that conservation form still behaves to conservative even after discretized. (Note: Conversion is done for non-dimensionalised form of conservative equations)
- `U_1=\rhoA`
- `U_2=\rhoAv`
- `U_3=\rho(e/(\gamma-1)+(\gamma/2)(v^2))A`
- `F_1=\rhoAv text( which is ) F_1=U_2`
- `F_2=\rhoAv^2+1/(\gamma)(pA) text (which can be written as )(U_2^2/U_1)+(1/(\gamma-1))(U_3-(\gamma/2)(U_2^2/U_1))`
- `F_3=\rho(e/(\gamma-1)+(\gamma/2)(v^2))vA+pAv text ( which can be written as ) \gamma(U_2U_3)/U_1-(\gamma-1)/(\gamma)(U_2^3/U_1^2)`
- `J=(1/\gamma)(p(\partialA)/(\partialx)) text( which can be written as ) (\gamma-1)/(\gamma)(U_3-\gamma/2 (U_2^2/U_1))(\partiallnA)/(\partialx)`
Continuity Equation,
`(\partialU_1)/(\partialt)=(\partialF_1)/(\partialt)`
Momentum Equation,
`(\partialU_2)/(\partialt)=(\partialF_2)/(\partialt)+J`
Energy Equation,
`(\partialU_3)/(\partialt)=(\partialF_3)/(\partialt)`
A programmable code must be written in this form only (i.e in the new variable forms). If any of the flow variable is written in primitive variable form, It is found from my personal and previous authors experience that, program never converges (results in imaginary values).
Once after each time step, these new variables must be converted back into actual flow variables to understand actual behaviour of the fluid at every node.
Results and Discussions:
The Boundary and initial conditions that were used for conservative and non-conservative form of equations are written below,
Non-Conservative form:
Initial conditions Flow variables as a function of space coordinates and 'x' varies from 0 to 3,
- `\rho=1-0.3146x`
- `T=1-0.2314x`
- `v=(0.1-1.09x)T^(0.5)`
Boundary conditions,
- `\rho_1=1 \text( at 1st node at any time step)`
- `T_1=1 \text( at 1st node at any time step)`
- `v_1=2v_2-v_3`
- `v_n=2v_(n-1)-v_(n-2)`
- `T_n=2T_(n-1)-T_(n-2)`
- `\rho_n=2\rho_(n-1)-\rho_(n-2)`
Conservative form:
Initial conditions Flow variables as a function of space coordinates and 'x' varies from 0 to 3,
- `if 0 le x le 0.5, \rho=1,T=1,v=0.59/(\rhoA)`
- `if 0.5 le x le 1.5, \rho=1-0.366(x-0.5),T=1-0.167(x-0.5),v=0.59/(\rhoA)`
- `if 1.5 le x le 3, \rho=0.634-0.3879(x-1.5),T=0.833-0.3507(x-1.5),v=0.59/(\rhoA)`
Boundary conditions,
- `(U_1)_1=A_1`
- `(U_2)_1=2(U_2)_2-(U_2)_3`
- `(U_3)_1=(U_1)_1(T_1/(\gamma-1) +(gamma)/2(v_1^2))`, `T_1=1` is constant value which must be sustituted
- `(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)`
`\Deltat`value to be substituted in equation-24 and that is estimated as follows for both conservative and non-conservative forms,
`\Deltat=C(\Deltax)/(v_i+a_i), a_i \text ( is acoustic wave speed) =sqrt(T_i)`. `\Deltat`value is estimated at every node and minimum among all the node is considered as final time step length. In this challenge courant number is considered as 0.5.
In either of the forms total of 31 nodes were considered which are spatially equally placed.
An edittable MATLAB programme is made for the converging and diverging nozzle.
The first progamme is main file,
clc
clearvars
n=input('enter number of nodes');
%% Conservative form
[y_n,A_fn,rho_fn,v_fn,x_fn,T_fn]=non_con(n);
[y_c,A_fc,rho_fc,v_fc,x_fc,T_fc]=conse(n);
figure(11)
plot(x_fn,rho_fn.*(A_fn.*v_fn),'k',x_fc,rho_fc.*(A_fc.*v_fc),'r');
title('Steady state mass-flow rate at different nodes');
xlabel('spatial coordinate');
ylabel('mass flow rate');
legend('non-conservative form','conservative form');
The second file is a function file for non-conservative form,
%clc
%clearvars
function [y,A_f,rho_f,v_f,x_f,T_f]=non_con(n)
%% Initialization for the problem and defining the variables
x_f=linspace(0,3,n);
dx=x_f(2)-x_f(1);
A_f=1+2.2*(x_f-1.5).^2;
%% initial conditions for non-conservative form
rho_f=1-0.3146*x_f;
T_f=1-0.2314*x_f;
v_f=(0.1+1.09*x_f).*(T_f.^0.5);
p_f=rho_f.*T_f;
C=0.5;
gamma=1.4;
error= 9e9;
iter=1;
% plot(x_f,T_f);
%% Problem solving Macormack method
while(error>1e-5)
%calculation wave speed
ac=(T_f).^(0.5);
% assigning old values to another matrix
v_o=v_f;
T_o=T_f;
rho_o=rho_f;
%% Predictor steps
% primitive variable slopes wrt space
for i=2:n-1
dv_p(i)=(v_f(i+1)-v_f(i))/dx;
drho_p(i)=(rho_f(i+1)-rho_f(i))/dx;
dT_p(i)=(T_f(i+1)-T_f(i))/dx;
dA_p(i)=(log(A_f(i+1))-log(A_f(i)))/dx;
end
% Primitive variable slopes wrt time
for i=2:n-1
drho_dt(i)=-(rho_f(i)*dv_p(i))-(rho_f(i)*v_f(i)*dA_p(i))-(v_f(i)*drho_p(i));
dv_dt(i)= -(v_f(i)*dv_p(i))-((dT_p(i)+(T_f(i)*drho_p(i)/rho_f(i)))/gamma);
dT_dt(i)=-v_f(i)*dT_p(i)-(gamma-1)*T_f(i)*(dv_p(i)+v_f(i)*dA_p(i));
end
% Time step calculation
dt_f=dx*C*(v_f+ac).^(-1);
dt=min(dt_f);
% Predictor values for rho,Vel and temp.
for i=2:n-1
T_f(i)=T_f(i)+(dt*dT_dt(i));
v_f(i)=v_f(i)+(dt*dv_dt(i));
rho_f(i)=rho_f(i)+(dt*drho_dt(i));
end
%% Correction steps
% primitive variable slopes wrt space
for i=2:n-1
dv_p(i)=(v_f(i)-v_f(i-1))/dx;
drho_p(i)=(rho_f(i)-rho_f(i-1))/dx;
dT_p(i)=(T_f(i)-T_f(i-1))/dx;
dA_p(i)=(log(A_f(i))-log(A_f(i-1)))/dx;
end
% correction slopes
for i=2:n-1
drho_dtp(i)=-(rho_f(i)*dv_p(i))-(rho_f(i)*v_f(i)*dA_p(i))-(v_f(i)*drho_p(i));
dv_dtp(i)= -(v_f(i)*dv_p(i))-((dT_p(i)+(T_f(i)*drho_p(i)/rho_f(i)))/gamma);
dT_dtp(i)=-v_f(i)*dT_p(i)-(gamma-1)*T_f(i)*(dv_p(i)+v_f(i)*dA_p(i));
end
%% Average slopes
% average slope values
drho_dta=0.5*(drho_dt+drho_dtp);
dv_dta=0.5*(dv_dt+dv_dtp);
dT_dta=0.5*(dT_dt+dT_dtp);
%% Updation
% updating temperature, density and velocity
for i=2:n-1
rho_f(i)=rho_o(i)+dt*drho_dta(i);
v_f(i)=v_o(i)+dt*dv_dta(i);
T_f(i)=T_o(i)+dt*dT_dta(i);
end
% Boundary values updation
v_f(1)=2*v_f(2)-v_f(3);
v_f(n)=2*v_f(n-1)-v_f(n-2);
rho_f(n)=2*rho_f(n-1)-rho_f(n-2);
T_f(n)=2*T_f(n-1)-T_f(n-2);
% Pressure and Mach number updation
p_f=rho_f.*T_f;
M_f=v_f.*((T_f).^-0.5);
%% Error calculation
%calculating maximum error
error1=max(max(abs(T_f-T_o)));
error2=max(max(abs(v_f-v_o)));
error3=max(max(abs(rho_f-rho_o)));
error_ma=[error1,error2,error3];
error=max(error_ma);
%% Plot representation
figure(1)
plot(x_f,rho_f.*(A_f.*v_f));
title(sprintf('iteration number-%d,non-conservative',iter))
xlabel('spatial coordinates');
ylabel('mass flowrate (non-dimensional form)');
pause(0.00001);
rho_ft(iter)=rho_f(15);
T_ft(iter)=T_f(15);
v_ft(iter)=v_f(15);
p_ft(iter)=p_f(15);
M_ft(iter)=M_f(15);
%hold on
iter=iter+1;
end
% Plots
figure(2)
subplot(3,1,1)
plot(x_f,rho_f);
title('flow variables vs space coordinate after reaching steady state,non-conservative');
xlabel('spatial coordinate');
ylabel('density @ steady state');
subplot(3,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,v_f);
xlabel('spatial coordinate');
ylabel('velocity @ steady state');
subplot(3,1,3)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,T_f);
xlabel('spatial coordinate');
ylabel('Temperature @ steady state');
figure(3)
subplot(2,1,1)
plot(x_f,M_f);
title('flow variables vs space coordinate after reaching steady state,non-conservative')
xlabel('spatial coordinate');
ylabel('Mach No. @ steady state');
subplot(2,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,p_f);
xlabel('spatial coordinate');
ylabel('pressure @ steady state');
y_f=1:1:iter-1;
figure(4)
subplot(3,1,1)
plot(y_f,rho_ft);
title('flow variables vs iterations @ node-15,non-conservative');
xlabel('time coordinate as iteration number');
ylabel('density @ node-15');
subplot(3,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,v_ft);
xlabel('time coordinate as iteration number');
ylabel('velocity @ node-15');
subplot(3,1,3)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,T_ft);
xlabel('time coordinate as iteration number');
ylabel('Temperature @ node--15');
figure(5)
subplot(2,1,1)
plot(y_f,M_ft);
title('flow variables vs iterations @ node-15,non-conservative')
xlabel('time coordinate as iteration number');
ylabel('Mach No. @ node-15');
subplot(2,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,p_ft);
xlabel('time coordinate as iteration number');
ylabel('pressure @ node-15');
% figure(5)
% plot(x_f,rho_f.*(A_f.*v_f))
% title('mass flowrate after reaching steady state at various nodes');
% xlabel('space coordinate');
% ylabel('mass flowrate');
y=iter-1;
end
The Third file is a function file for conservative form,
%clc
%clearvars
function [y,A_f,rho_f,v_f,x_f,T_f]=conse(n)
% defining inputs
x_f=linspace(0,3,n);
dx=x_f(2)-x_f(1);
A_f=1+2.2*(x_f-1.5).^2;
for i=1:n
if i<=5
rho_f(i)=1;
T_f(i)=1;
elseif i>5 && i<=16
rho_f(i)=1-0.366*(x_f(i)-0.5);
T_f(i)=1-0.167*(x_f(i)-0.5);
elseif i>16
rho_f(i)=0.634-0.3879*(x_f(i)-1.5);
T_f(i)=0.833-0.3507*(x_f(i)-1.5);
end
v_f(i)=0.59/(rho_f(i)*A_f(i));
end
p_f=rho_f.*T_f;
C=0.5;
gamma=1.4;
error= 9e9;
iter=1;
U_1=rho_f.*A_f;
U_2=U_1.*v_f;
U_3=U_1.*((T_f/(gamma-1))+(gamma/2)*v_f.*v_f);
% Defining functions
iter=1;
while(error>1e-5)
%calculation wave speed
ac=(T_f).^(0.5);
% assigning old values to another matrix
U_1o=U_1;
U_2o=U_2;
U_3o=U_3;
F_1=U_2;
F_2=((U_2.^2).*U_1.^(-1))+((gamma-1)/gamma)*(U_3-((gamma*(U_2.^2).*U_1.^(-1))/2));
F_3=(gamma*(U_2.*U_3).*U_1.^(-1))-(gamma*(gamma-1)*((U_2.^3).*(U_1.^(-2)))/2);
K=(U_1*(gamma-1)/(gamma)).*((U_3.*(U_1.^(-1)))-(gamma/2)*(U_2.*(U_1.^(-1))).^2);
%% F1 F2 F3
for i=2:n-1
du1_dtp(i)=-(F_1(i+1)-F_1(i))/dx;
du2_dtp(i)=-((F_2(i+1)-F_2(i))/dx)+K(i)*(A_f(i+1)-A_f(i))/(dx*A_f(i));
du3_dtp(i)=-(F_3(i+1)-F_3(i))/dx;
end
% Time step calculation
dt_f=dx*C*(v_f+ac).^(-1);
dt=min(dt_f);
% Predictor values for rho,Vel and temp.
for i=2:n-1
U_1(i)=U_1o(i)+(dt*du1_dtp(i));
U_2(i)=U_2o(i)+(dt*du2_dtp(i));
U_3(i)=U_3o(i)+(dt*du3_dtp(i));
end
F_1=U_2;
F_2=((U_2.^2).*U_1.^(-1))+((gamma-1)/gamma)*(U_3-((gamma*(U_2.^2).*U_1.^(-1))/2));
F_3=(gamma*(U_2.*U_3).*U_1.^(-1))-(gamma*(gamma-1)*((U_2.^3).*(U_1.^(-2)))/2);
K=(U_1*(gamma-1)/(gamma)).*((U_3.*(U_1.^(-1)))-(gamma/2)*(U_2.*(U_1.^(-1))).^2);
for i=2:n-1
du1_dtc(i)=-(F_1(i)-F_1(i-1))/dx;
du2_dtc(i)=-((F_2(i)-F_2(i-1))/dx)+K(i)*(A_f(i)-A_f(i-1))/(dx*A_f(i));
du3_dtc(i)=-((F_3(i)-F_3(i-1))/dx);
end
% average slopes
du1_dta=0.5*(du1_dtp+du1_dtc);
du2_dta=0.5*(du2_dtp+du2_dtc);
du3_dta=0.5*(du3_dtp+du3_dtc);
% corrector values
for i=2:n-1
U_1(i)=U_1o(i)+(dt*du1_dta(i));
U_2(i)=U_2o(i)+(dt*du2_dta(i));
U_3(i)=U_3o(i)+(dt*du3_dta(i));
end
% Boundary condtions
U_1(1)=A_f(1);
U_2(1)=2*U_2(2)-U_2(3);
v_f(1) = U_2(1)/U_1(1);
U_3(1)=U_1(1)*((T_f(1)/(gamma-1))+(gamma/2)*v_f(1)^2);
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);
% representing back values into density,temp and velocity
rho_f=U_1.*(A_f).^(-1);
v_f=U_2.*(U_1).^(-1);
T_f=(gamma-1)*(U_3.*U_1.^(-1)-((gamma/2)*(v_f.^2)));
p_f=rho_f.*T_f;
M_f=v_f.*((T_f).^-0.5);
%calculating maximum error
error1=max(max(abs(U_1-U_1o)));
error2=max(max(abs(U_2-U_2o)));
error3=max(max(abs(U_3-U_3o)));
error_ma=[error1,error2,error3];
error=max(error_ma);
% plot(x_f,T_f);
% pause(0.3);
figure(6)
plot(x_f,rho_f.*(A_f.*v_f));
title(sprintf('iteration number-%d,conservative',iter))
xlabel('spatial coordinates');
ylabel('mass flowrate (non-dimensional form)');
pause(0.00001);
% Primitive variables are stored as a function of time
T_ft(iter)=T_f(15);
v_ft(iter)=v_f(15);
rho_ft(iter)=rho_f(15);
M_ft(iter)=M_f(15);
p_ft(iter)=p_f(15);
dtt(iter)=dt;
iter=iter+1;
end
% Plots
figure(7)
subplot(3,1,1)
plot(x_f,rho_f);
title('flow variables vs space coordinate after reaching steady state,conservative');
xlabel('spatial coordinate');
ylabel('density @ steady state');
subplot(3,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,v_f);
xlabel('spatial coordinate');
ylabel('velocity @ steady state');
subplot(3,1,3)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,T_f);
xlabel('spatial coordinate');
ylabel('Temperature @ steady state');
figure(8)
subplot(2,1,1)
plot(x_f,M_f);
title('flow variables vs space coordinate after reaching steady state,conservative')
xlabel('spatial coordinate');
ylabel('Mach No. @ steady state');
subplot(2,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(x_f,p_f);
xlabel('spatial coordinate');
ylabel('pressure @ steady state');
y_f=1:1:iter-1;
figure(9)
subplot(3,1,1)
plot(y_f,rho_ft);
title('flow variables vs iterations @ node-15,conservative');
xlabel('time coordinate as iteration number');
ylabel('density @ node-15');
subplot(3,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,v_ft);
xlabel('time coordinate as iteration number');
ylabel('velocity @ node-15');
subplot(3,1,3)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,T_ft);
xlabel('time coordinate as iteration number');
ylabel('Temperature @ node--15');
figure(10)
subplot(2,1,1)
plot(y_f,M_ft);
title('flow variables vs iterations @ node-15,conservative')
xlabel('time coordinate as iteration number');
ylabel('Mach No. @ node-15');
subplot(2,1,2)
%title(flow variable vs space coordinate after reaching steady state)
plot(y_f,p_ft);
xlabel('time coordinate as iteration number');
ylabel('pressure @ node-15');
% figure(12)
% plot(x_f,rho_f.*(A_f.*v_f))
% title('mass flowrate after reaching steady state at various nodes');
% xlabel('space coordinate');
% ylabel('mass flowrate');
y=iter-1;
end
Initial results were obtained for both conservative and non-conservative forms of governing equations when nozzle is equally divided into 30 elements (or 31 nodes).
Conservative form for Steady state results of flow variables are shown below :
Note: In the upcoming statements, if any Figure number is mentioned that indicates Results plots Screenshot number which is on top left of the screen shot.
Non-Conservative form for Steady state results of flow variables are shown below:
In Conservative form at node-15, time wise variation of primitive variables are plotted:
In Non-Conservative form at node-15, time wise variation of primitive variables are plotted below:
Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (conservative approach) are shown below, the title indicates the iteration number. It took 1085 iterations to reach steady in conservative form (this information can be drawn from Figure(9) and (10).
Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process (non-conservative approach) are shown below, the title indicates the iteration number. It took 731 iterations to achieve steady state (this information can be clearly drawn from Figure(4) and (5).
Comparison of Normalized mass flow rate distributions of both forms:
Observations from above results and literature survey:
| Conservative form | Non-conservative form |
| To achieve Steady results, it took 1085 iterations (more no.of time steps). (By considering figure (9,10) | To achieve Steady results, it took 731 iterations(less number of iterations).(By considering figure (4,5) |
| The above comparision also infers, error reduction rate between two consecutive time steps is low. | error reduction rate is high between two consecutive time steps. Here in the above case, error to converge it is considered as 1e-5 |
| All flow variables are in conservative form are in close agreement with Non-conservative form except the mass flow rate after reaching steady state, this indicates final steady state flow variables may or may not exactly match with analytical steady state results . The steady state mass flow rate error (this is a function of spatial coordinate) is found to have minimal error when compared with analytical steady state error. | All the final flow variable values which are found in non-conservative form are having almost the close values which are found conservative form, as already said there is minimal amount of mis-match with mass flow rate which is shown in figure 11. This infers, that conservative form is very essential to get a close agreement values with actual analytical state values of mass flowrates. |
| From literature, it was known that this can handle even if discontinuties exists with respect to flow variables. for example: shocks in the nozzle | This form cannot handle if any discontinuites or non-differentiable forms exists. |
Even grid independence test is also done, By increasing number of nodes to 61 in place of 31 nodes. It is found that, flow variables does not vary upto 3 decimals. That indicates, 31 nodes can be considered as optimum values for both the forms. In the below plots of various flow variables, for both 31 and 61 nodes they are in close agreement after reaching steady state. In correspondence to that,below plots and tables are made.
Note: As similar to 31 nodes case, flow variable plots are generated for 61 nodes too, but too make the report crunch enough those plots are not inputted here.
(For complete table, excel sheet has been attached)