## Simulation of a 1D Quasi Super-sonic nozzle flow simulation using Maccormack Method

Google drive link containing source code for both functions and all the plots

Observations and Conclusions:

1. Normalized Mass flow rate

We see that in case of conservative form of governing equations the mass flow rate along the nozzle is almost constant along the length of the nozzle. In contrast, the non conservation form results have sizeable variations and some oscillations at both inflow and outflow boundaries.

Also we see that the steady state mass flow results obtained with the conservation form are more close to analytical solution(0.579).

Since the mass flow rate obtained in conservative form is almost constant so we say that the mass flow is conserved as so the name conservative form of governing equations.

This behaviour of conservative form of governing equations can also be justified as we take mass flow rate as one of dependent variable i.e. U2=ρ'*A'*V' .

Non Conservative Case

Conservative Case

2. Mass Flow rate per unit area

The mass flow rate per unit area for both the cases are shown and it is in accordance with theoretical results i.e. the maximum mass flow rate per unit area occurs at the throat of the nozzle(x=1.5).

Non conservative

Conservative

3. Convergence Plots for various properties at the throat

Non-Conservative

Conservative

4. Variation of Pressure,density,velocity and temperature along the length of nozzle

We achieve identical plots for the properties for both the cases

Non Conservative Case

Conservative Case

5. Mach Number Variation

Similar to the previous plots we achieve ientical plots for the variation of mach number along the nozzle with mach number=1 at throat(x=1.5) of the nozzle.

Non-Conservative Case

Conservative Case

6. Minimum Number of Cycles for Convergence

The minimum number of cycles a simulation requires to converge depend upon following:

1) Courant Number(C)-which determines the time step size(dt)

dt=min(C*dx/(a+v))

where a=sqrt(T);      non dimensionlalized speed of sound

v=velocity;     non dimensionalized velocity along the nozzle

2) Error Tolerance/Residual

3) Spacing between two gridpoints(dx)/Number of Gridpoints(n) for a given domain

For the given simulation, considering the following criteria for the given parameters:

1) Courant Number=0.5

2) Error Tolerance=1e-06

3) Grid spacing(dx)=0.1;

No. of gridpoints(n)=31;

For the following the minimun number of cycles/minimum number of time steps required by the simulation to run were found out to be the following

Non-Conservative Case: 940 time steps

Conservative Case: 821 time steps

7. Simulation Time required

For the parameters given above for C,dx and error, the simulation time required by matlab to converge the simulation was measured by the inbuilt command in matlab (Run and Time) and the results for the two cases were:

1. Non-Conservative Case: 5.280 s

2. Conservative Case: 3.773 s

So we conclude by the results obtained that the simulation time required by conservative form of governing equations is less than that for non conservative case. The payoff for the same is it takes more time to code conservative form than the non conservative form as we have to decode the primitive variables (rho, Velovcity and Temperature) from the dependent variables ( U1, U2 and U3) after each time step.

8. Grid Dependence

As can be seen in the plots given in first point of mass flow rates that there are variations in the mass flow rate along the length of nozzle which is very less for conservative case.

Theoretically the mass flow rate should be constant and a flat line.

It was seen that as we refine the number of grid points we start getting more and more flat profile for the mass flow rate with a completely flat profile at about 301 gridpoints.

The following data was obtained for mass flow rate at the Throat from the grid dependence test:

 Time Steps Mass Flow Rate Gridpoints Conservative Form Non-Conservative Form 31 821 0.5853 940 0.5837 61 1680 0.5829 1716 0.5825 91 2311 0.5824 2533 0.5823 151 3715 0.5821 3819 0.5819 301 6575 0.5822 7399 0.5819 Analytical Mass Flow Rate 0.579

So we see that the results obtained from the simulation for 31 gridpoints are good enough if we can tolerate slight errors.

If we want more accurate results than we see that as we increase the number of gridpoints the Mass flow rate approaches more towards the analytical result and the simulation finally becomes grid independent at 151 gridpoints.

9. Source Code for the Non Conservative Form of Governing Equations

%Simulation of 1D Quasi Steady Supersonic Flow Through Nozzle
%Non-Conservative Form of governing equations
%MacCormack Method

clear all
close all
clc

%Inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
nt=1400;    %Time Steps
error_tolerance=1e-6;

%Initial Profiles
rho=1-0.3146*x;    %Density
t=1-0.2314*x;      %temperature
v=(0.1+1.09*x).*t.^0.5;   %Velocity

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

%time step size
C=0.5; %cfl
del_t=C*dx./(sqrt(t)+v);
dt=min(del_t);

rho_initial=rho;
v_initial=v;
t_initial=t;

drhodt_p=zeros(1,n);
dvdt_p=zeros(1,n);
dtdt_p=zeros(1,n);
drhodt_c=zeros(1,n);
dvdt_c=zeros(1,n);
dtdt_c=zeros(1,n);
error=9e9;

%Time Loop
for k=1:nt
%at beginning of time step
rho_old=rho;
t_old=t;
v_old=v;

M=v./sqrt(t);
p=rho.*t;

%For Plotting variables at throat
rho_plot(k)=rho(16);
t_plot(k)=t(16);
p_plot(k)=p(16);
M_plot(k)=M(16);

%Non-Conservative, Non- Dimensionalized form of governing equations

%Continuity Equation

%Momentum Equation
%dvdy=-v*dvdx-(1/gamma)*(dtdx+t*drhodx/rho)

%Energy Equation

for j=2:n-1
%Predictor Method
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;

%Forward Differencing space terms in continuity equation

%Forward Differencing space terms in momentum equation
dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+t(j)*drhodx/rho(j));

%Forward Differencing space terms in energy equation

%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

for j=2:n-1
%Corrector Method
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dtdx=(t(j)-t(j-1))/dx;

%Backward Differencing space terms in continuity equation

%Backward Differencing space terms in momentum equation
dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+t(j)*drhodx/rho(j));

%Backward Differencing space terms in energy equation
end

%computing average derivative
drhodt_av=0.5*(drhodt_p+drhodt_c);
dvdt_av=0.5*(dvdt_p+dvdt_c);
dtdt_av=0.5*(dtdt_p+dtdt_c);

%final calculation of values at next time step
rho=rho_old+drhodt_av*dt;
v=v_old+dvdt_av*dt;
t=t_old+dtdt_av*dt;

%Boundary conditions

%Inflow
v(1)=2*v(2)-v(3);

%Outflow
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);

error_rho=max(abs(rho-rho_old));
error_t=max(abs(t-t_old));
error_v=max(abs(v-v_old));
error=max([error_rho;error_t;error_v]);

if error

10. Source Code for the Conservative Form of Governing Equations

%Simulation of 1D Quasi Steady Supersonic Flow Through Nozzle
%Conservative Form of governing equations
%MacCormock Method

clear all
close all
clc

%Inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
nt=1400;
error_tolerance=1e-6;

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

%Density and Temperature

rho=zeros(1,n);
t=zeros(1,n);

%Initial Conditions
for i=1:n
rho(i)=1;
t(i)=1;
if x(i)==0.5
break;
end
end

for j=i:n
rho(j)=1-0.366*(x(j)-0.5);
t(j)=1-0.167*(x(j)-0.5);
if x(j)==1.5
break;
end
end

for k=j:n
rho(k)=0.634-0.3879*(x(k)-1.5);
t(k)=0.833-0.3507*(x(k)-1.5);
if x(k)==3.5
break;
end
end

%Velocity
v=((rho.*a)/0.59).^-1;
e=t;
%Initial Conditions for solution vectors
u1=rho.*a;
u2=rho.*a.*v;
u3=rho.*((e/(gamma-1))+(gamma/2)*v.^2).*a;

%timestep=min(C*dx/vel)
%for cfl number 0.5
dt=0.0267;
%time step size
% C=0.5; %cfl
% del_t=C*dx./(sqrt(t)+v);
% dt=min(del_t);

f1=zeros(1,n);
f2=zeros(1,n);
f3=zeros(1,n);
j2=zeros(1,n);
du1dt_p=zeros(1,n);
du2dt_p=zeros(1,n);
du3dt_p=zeros(1,n);
du1dt_c=zeros(1,n);
du2dt_c=zeros(1,n);
du3dt_c=zeros(1,n);

u1_initial=u1;
u2_initial=u2;
u3_initial=u3;

%Dependent Variables
%u1=rho*a
%u2=rho*a*v
%u3=rho*((e/(gamma-1))+(gamma/2)*v^2)*a

%e=t as non dimensionalized

%Time Loop
for l=1:nt
%Non dimensionalized pressure and mach number
p=rho.*t;
M=v./sqrt(t);

%for plotting variables at throat
rho_plot(l)=rho(16);
t_plot(l)=t(16);
p_plot(l)=p(16);
M_plot(l)=M(16);

%Working with dependent variables u1,u2 and u3 instead of primitive variables rho,t,v
u1_old=u1;
u2_old=u2;
u3_old=u3;

%Computing the flux and source terms in pure form
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f3=(gamma.*u2.*u3./u1)-((gamma*(gamma-1)*u2.^3)./(2*u1.^2));

f1_old=f1;
f2_old=f2;
f3_old=f3;

%Predictor Step
for m=2:n-1
%Modelling the governing equation
%du1dt=-df1dx
%du2dt=-df2dx+j2
%du3dt=-df3dx

%Forward Differncing
df1dx=(f1(m+1)-f1(m))/dx;
df2dx=(f2(m+1)-f2(m))/dx;
df3dx=(f3(m+1)-f3(m))/dx;

%Governing Equations
du1dt_p(m)=-df1dx;        %Continuity
du2dt_p(m)=-df2dx+j2;     %Momentum
du3dt_p(m)=-df3dx;        %Energy

%Solution Update
%Predicted Values of dependent variables/solution vectors
u1(m)=u1(m)+du1dt_p(m)*dt;
u2(m)=u2(m)+du2dt_p(m)*dt;
u3(m)=u3(m)+du3dt_p(m)*dt;
end

%Predicted Values of flux terms using predicted values of u1,u2 and u3
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f3=(gamma.*u2.*u3./u1)-((gamma*(gamma-1)*u2.^3)./(2*u1.^2));

%Corrector Step
for p=2:n-1
%Backward differencing for flux terms
df1dx=(f1(p)-f1(p-1))/dx;
df2dx=(f2(p)-f2(p-1))/dx;
df3dx=(f3(p)-f3(p-1))/dx;

%Governing Equations
du1dt_c(p)=-df1dx;        %Continuity
du2dt_c(p)=-df2dx+j2;     %Momentum
du3dt_c(p)=-df3dx;        %Energy
end

%Computing average derivatives
du1dt_av=0.5*(du1dt_p+du1dt_c);
du2dt_av=0.5*(du2dt_p+du2dt_c);
du3dt_av=0.5*(du3dt_p+du3dt_c);

%Final Corrected Values of dependent variables u1,u2 and u3 at next time step
u1=u1_old+(du1dt_av*dt);
u2=u2_old+(du2dt_av*dt);
u3=u3_old+(du3dt_av*dt);

%Applying Boundary Conditions
%Inflow Boundary
u1(1)=rho(1)*a(1);      %Fixed Value as rho and area at inflow is fixed
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));

%Outflow boundary
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);

%Finally obtaining the values of primitive variables by decoding u1,u2 and u3
rho=u1./a;
v=u2./u1;
t=(gamma-1)*((u3./u1)-(gamma/2)*(v.^2));

error_u1=max(abs(u1-u1_old));
error_u2=max(abs(u2-u2_old));
error_u3=max(abs(u3-u3_old));

error=max([error_u1; error_u2; error_u3]);

if error

### Flow over a backward facing step Utkarsh Garg · 2019-01-19 07:43:25

The backward facing step flow project was performed as a part of skill lync challenge. The simulation was performed on three base grid sizes: 1. 3e-4m 2. 2.5e-4m 3. 2e-4m The above grid sizes were chosen due to lack of computation power.   The Results a Read more

### Channel Flow using CONVERGE STUDIO Utkarsh Garg · 2019-01-04 18:07:51

The simulation was run by refining the grid size successively in three steps as follows: Grid Size 1. 2e-4 m 2. 1.5e-4 m 3. 1.2e-4 m   The post-processing results of the three cases are as follows: 1. 2e-4 m Velocity Magnitude Mass Flow Rate Spec Read more

### Centrifugal Pump design and Performance Analysis Utkarsh Garg · 2018-10-11 17:14:09

CENTRIFUGAL PUMP DESIGN AND PERFORMANCE ANALYSIS   1. Geometry: PUMP   IMPELLER      2. Performance Analysis of the Centrifugal Pump Performance Analysis of a Centrifugal Pump Velocity normal to face (Outlet Velocity) [m/ Read more

### Flowbench Simulation for Intake Port Utkarsh Garg · 2018-08-12 12:20:33

1. Geometry: All dimensions in mm.   2. Boundary Conditions: Total pressure inlet = 4 bar Static pressure outlet = 1 bar Real wall = all inner surfaces of the cylinder, valve, and intake.   3. Grid dependence study: Goal (Value) Design Point Read more

### Symmetry vs Wedge boundary conditions openFoam Utkarsh Garg · 2018-08-12 08:38:12

Project Report   Assumptions (same as in previous challenge): 1. Reynolds Number = 1000 (For Laminar Flow) 2. Temperature= 20 degree Celsius 3. Length of pipe = 1 m 4. Diameter of pipe = 20 mm 5. At 20 Degrees Celsius (Ref: Saturated Water tables from Cengel Read more

### Iterative Solvers for linear systems Utkarsh Garg · 2018-07-22 10:16:07

Google Drive Link for all sorce code files and functions files https://drive.google.com/open?id=1UE-_9nn7VMuXmUqr3XWFgwbxaaGaNpdX   Outputs: 1. Jacobi Method The program also showed the following output in the command window \"The Jacobi Solver could not co Read more

### Flow over NACA0008 for various Angles of attack Utkarsh Garg · 2018-06-25 18:41:04

Google drive link to all the contours and plot: https://drive.google.com/open?id=1fBI5Zcassvv2ut9o2Qsy7ktrIq3LQx10 The link contains all the plots for various angles of attack and the pressure and velocity distribution for all cases.   The airfoil used was NACA Read more