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

https://drive.google.com/open?id=1q-2do7AV4YP2U5ilL54n37U3CXuNoeNK

 

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
   %drhodt=-rho*(dvdx)-rho*v*(dlogadx)-v*(drhodx)
   
   %Momentum Equation
   %dvdy=-v*dvdx-(1/gamma)*(dtdx+t*drhodx/rho)
   
   %Energy Equation
   %dtdt=-v*dtdx-(gamma-1)*t(dvdx+v*dlogadx)
   
   for j=2:n-1 
       %Predictor Method
      dvdx=(v(j+1)-v(j))/dx;
      dlogadx=(log(a(j+1))-log(a(j)))/dx;
      drhodx=(rho(j+1)-rho(j))/dx;
      dtdx=(t(j+1)-t(j))/dx;
      
      %Forward Differencing space terms in continuity equation
      drhodt_p(j)=-rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
      
      %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
      dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
      
      
      %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;
      dlogadx=(log(a(j))-log(a(j-1)))/dx;
      drhodx=(rho(j)-rho(j-1))/dx;
      dtdx=(t(j)-t(j-1))/dx;
      
      %Backward Differencing space terms in continuity equation
      drhodt_c(j)=-rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;
      
      %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
      dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);    
   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;
       dadx=(a(m+1)-a(m))/dx;
       j2=rho(m)*t(m)*dadx/gamma;
       
       %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;
        dadx=(a(p)-a(p-1))/dx;
        j2=rho(p)*t(p)*dadx/gamma;
        
        %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

Projects by Utkarsh Garg

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

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   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

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

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

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

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

Hagen Poiseuille s Flow
Utkarsh Garg · 2018-06-15 20:28:52

Google drive link to all the blockMeshDict file for various angles: https://drive.google.com/open?id=1GuYMGjgkJHC8hovIfkqA1sp4k1OdNR1t Google drive link to all the plots and source code: https://drive.google.com/open?id=1ntyenSwIed9JrNxDXwMKVZIeMk0-xNrk   Assu Read more


Loading...

The End