Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method

Aim:-

Simulation of a 1D supersonic nozzle flow simulation using

in this challenge you are going to simulation the isortopic flow through a 1D subsonic-supersonic nozzle .

you will derive both the conservation and non -conservation forms of the governingquatiorand solve them using the macrocosm's technique.

you need to determine the steady-state temperature distribution for the flow -field variable and investigate the difference between the two forms of governing quotations by comparing this solutions.

you need to show the following plots inside your report

1.Steady-State distribution of primitive variablrs inside the nozzle

2.Time -wise variation of the primitive variables

3.Variation of mass flow rate disteibution inside the nozzle at different tome steps during the time-marching process

4.Comparison if Normalized flow through the convergent nozzle as give in diagram below

INTRODUCTION TO THE GIVEN PROBLEM

We consider the steady isntropic flow theougth the convergent -divergent as given in diagram below .

We consider the steady isontropic flow throughr the convergent -divergent nozzle aas given in diagram below

Momentum equation

Energy equation

Calculation of Average

its basically the half of the sum of corrector and predictor value.

2)

3)

Density and Temperature at point 1 is held fixed.

The momentum equation

The energy equation :

The above are the non dimansional conservation form of the coninutiy ,Momentum and energy equation The quasi-one -dimansinal floe ,respectiveley.The equations for quasi one dimensional flow can be expressed in generic form. Let us define the elements of the solution vector U, flux vector F and the source term J as follows

From the above equations we can write that

These are the equations we wish to solve numerically by using the MacCormacks technique.

To obtain the primitve variables we must decode the elements U1,U2,U3 as follows.From the defination of U1,U2 and U3 given above we have

Expressing the flux and source term in form of solution vectors

Applying the boundry conditions:

Initial conditi

%inputs

n=31;

x = linspace(0,3,n);

dx = x(2)-x(1);

%calculate the initial profiles

rho = 1-0.3146*x;

t = 1-0.2314*x;

v =(0.1+1.09*x).*t.^0.5;

gamma = 1.4:

%calculate the are

a = 1+2.2*(x-1.5).^2;

pos = find(a==min(a))

%time step

nt = 1400;

cfl=0.5;

dt = min (cfl.*dx/sqrt(t)+v));

%outer time loop

for k = 1:nt

%preditor method

rho_old=rho;

v_old=v;

t-old=t;

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;

%continuty equation

drhdt_p(j) =-rho(j)*dvdx-rho(j)*(j(*dlogadx-v(j)*drhodx;

%momentum equation

dvdt_p(j)*dvdx-(1-gamma)*(ddtdx+(t(j)/rho(j)*drhodx);

%energy ewuation

dtdt_p(j)= -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);

%solution update

v(j)=v(j)+dvdt_p(j)*dt;

rho(j)=rho(j)-drhodt_p(j)*dt;

t(j)=t(j)+dtdt_p(j)*dt;

end

%corrector method

for j =2:n-1

dvdx =(v(j)-v(j-1))/dx;

drhodx=(rho(j)-rho(j-1))/dx;

dtdx=(t(j-1))/dx;

dlogadx=(log(a(j))-loga(j-1)))/dx;

%continuty equation

drhodt_c(j)=-rho(j)*dvdx-rho(j)*cv(j)*drhodx;

%momentum equation

dvdt_c(j) =-vdvdx-(1/gamma)*(dtdx+(t(j)/(j))*drhodx);