Skill-Lync Launch Pad – Your Gateway to a Core Engineering Job by 2025! Only 1̶0̶0̶ +50 Seats Available.

00D 14H 45M 50S

Menu

Executive Programs

Workshops

Projects

Blogs

Careers

Placements

Student Reviews

For Business

More

Academic Training

Informative Articles

Find Jobs

We are Hiring!

All Courses

Choose a category

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

All Courses

All Courses

logo

CHOOSE A CATEGORY

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

Top Job Leading Courses

Automotive

CFD

FEA

Design

MBD

Med Tech

Courses by Software

Design

Solver

Automation

Vehicle Dynamics

CFD Solver

Preprocessor

Courses by Semester

First Year

Second Year

Third Year

Fourth Year

Courses by Domain

Automotive

CFD

Design

FEA

Tool-focused Courses

Design

Solver

Automation

Preprocessor

CFD Solver

Vehicle Dynamics

Machine learning

Machine Learning and AI

POPULAR COURSES

coursePost Graduate Program in Hybrid Electric Vehicle Design and Analysis
coursePost Graduate Program in Computational Fluid Dynamics
coursePost Graduate Program in CAD
coursePost Graduate Program in CAE
coursePost Graduate Program in Manufacturing Design
coursePost Graduate Program in Computational Design and Pre-processing
coursePost Graduate Program in Complete Passenger Car Design & Product Development
Executive Programs
Workshops
For Business

Success Stories

Placements

Student Reviews

More

Projects

Blogs

Academic Training

Find Jobs

Informative Articles

We're Hiring!

phone+91 9342691281Log in
  1. Home/
  2. Amol Patel/
  3. Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

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

Simulation of 1D supersonic nozzel flow using Macormack method Objectives:1. Steady-state distribution of premitive variable inside the nozzel.2. Time-wise distribution of premitive variables.3. Variation of mass flow rate distribution inside the nozzel at different time steps during the time marching process.4. Comparison…

    • Amol Patel

      updated on 03 Jul 2021

    Simulation of 1D supersonic nozzel flow using Macormack method

    Objectives:
    1. Steady-state distribution of premitive variable inside the nozzel.
    2. Time-wise distribution of premitive variables.
    3. Variation of mass flow rate distribution inside the nozzel at different time steps during the time marching process.
    4. Comparison of normalised mass flow rate distribution of both forms.

     

    We will be simulating the isentropic flow through a Quasi 1D subsonic-supersonic nozzel. For this simulation we will be using both forms of governing equations, namely non-conservative and conservative forms. We will first go through the idea behid each form and then we will see the code to simulate the quasi 1D nozzel flow.

    The nozzel model has a convergent section as the inlet and the divergent section as the outlet, the inlet is attached to a resiorvior so that we can assume that teh velocity of flow at the inlet is approaching zero and the area tends to infinity the pressure and temperature  at the inlet are the total pressure and total temperature for the resiorvior. The mach number for the inlet is less than one that is the flow is subsonic. The outlet of the nozzel is supersonic so the mach number will be greater than one. Also at the throat of the nozzel we assume that the flow has a mach number value equal to unity.

    Consider the following image for the quasi 1D supersonic nozzel flow

     

    Non-Conservative Form: 

    Now we will see the non conservative form of governing equation that are associated with this flow.

    Continuity equation:

    ∂ρ∂t=−ρ∂v∂x−ρv∂lna∂x−v∂ρ∂x∂ρ∂t=-ρ∂v∂x-ρv∂lna∂x-v∂ρ∂x

    Momentum equation:

    ∂v∂t=−v∂v∂x−1γ(∂T∂x+tTρ∂ρ∂x)∂v∂t=-v∂v∂x-1γ(∂T∂x+tTρ∂ρ∂x)

    Enegry equation:

    ∂T∂t=−v∂T∂x−(γ−1)T(∂v∂x+v∂(lna)∂x)∂T∂t=-v∂T∂x-(γ-1)T(∂v∂x+v∂(lna)∂x)

    here all the variable are in nondimensionalized forms.

    We willbe using macormack method for our numerical simulations

    the time step is taken as: dt=min⎛⎜⎝CFL⋅dx(√T)+V⎞⎟⎠dt=min(CFL⋅dx(T)+V)

    Now we will see the function code for the non conservative form:

    function [rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n, nt,x, gamma, CFL, a, dx, throat)

%% function for Quasi 1D nozzel flow using non-conservative form 

% calculate initial profiles for non conservative form
rho_n = 1 - 0.3146*x;             % density
t_n = 1- 0.2314*x;                % temperature
v_n = (0.1 + 1.09*x).*t_n.^0.5;   % velocity



% time loop
for k = 1:nt
    
    % minimum time step value among all the nodes approach should be
    % employed.
    dt = min(CFL*dx./((t_n.^0.5)+v_n));
    
    
    rho_old = rho_n;
    v_old = v_n;
    t_old = t_n;
    
    %% predictor method
    
    for j = 2:n-1
        
        % spatial derivativs
        dvdx = (v_n(j+1) - v_n(j))/dx;
        drhodx = (rho_n(j+1) - rho_n(j))/dx;
        dlogadx = (log(a(j+1)) - log(a(j)))/dx;
        dtdx = (t_n(j+1) - t_n(j))/dx;
        
        % predictive step
        % continuity equation
        drho_dt_p(j) = - rho_n(j)*dvdx - rho_n(j)*v_n(j)*dlogadx - v_n(j)*drhodx;
        % momentum equation
        dv_dt_p(j) = - v_n(j)*dvdx - (1/gamma)*(dtdx + (t_n(j)/rho_n(j)*drhodx));
        % energy equation
        dt_dt_p(j) = -v_n(j)*dtdx - (gamma -1)*t_n(j)*(dvdx + v_n(j)*dlogadx);
        
        % solution update
        v_n(j) = v_n(j) + dv_dt_p(j) * dt;
        rho_n(j) = rho_n(j) + drho_dt_p(j) * dt;
        t_n(j) = t_n(j) + dt_dt_p(j) * dt;
        
    end
    
    %% corrector method
    
    for j = 2:n-1
        
        % spatial derivatives
        dvdx = (v_n(j) - v_n(j-1))/dx;
        drhodx = (rho_n(j) - rho_n(j-1))/dx;
        dlogadx = (log(a(j)) - log(a(j-1)))/dx;
        dtdx = (t_n(j) - t_n(j-1))/dx;
        
        % corrector step
        % continuity equation
        drho_dt_c(j) = - rho_n(j)*dvdx - rho_n(j)*v_n(j)*dlogadx - v_n(j)*drhodx;
        % momentum equation
        dv_dt_c(j) = - v_n(j)*dvdx - (1/gamma)*(dtdx + (t_n(j)/rho_n(j)*drhodx));
        % energy equation
        dt_dt_c(j) = -v_n(j)*dtdx - (gamma -1)*t_n(j)*(dvdx + v_n(j)*dlogadx);
        
    end
    
    %% compute the average time derivative
    
    dv_dt = 0.5*(dv_dt_p + dv_dt_c);
    drho_dt = 0.5*(drho_dt_p + drho_dt_c);
    dt_dt = 0.5*(dt_dt_p + dt_dt_c);
    
    
    %% final solution update
    for i = 2:n-1
        v_n(i) = v_old(i) + dv_dt(i) * dt;
        rho_n(i) = rho_old(i) + drho_dt(i) * dt;
        t_n(i) = t_old(i) + dt_dt(i) * dt;
    end
    
    %% apply boundary conditions
    % inlet
    v_n(1) = 2* v_n(2) - v_n(3);
    
    % outlet
    v_n(n) = 2*v_n(n-1) - v_n(n-2);
    rho_n(n) = 2*rho_n(n-1) - rho_n(n-2);
    t_n(n) = 2*t_n(n-1) - t_n(n-2);
    
    % other properties 
    p_n = rho_n.*t_n;           % pressure 
    M_n = v_n./(t_n.^0.5);      % Mach number 
    mfr_n=rho_n.*v_n.*a;        % mass flow rate
    
    % calculating the values for throat
    rho_th_n(k) = rho_n(throat);
    v_th_n(k) = v_n(throat);
    t_th_n(k) = t_n(throat);
    p_th_n(k) = p_n(throat);
    M_th_n(k) = M_n(throat);
    mfr_th_n(k) = mfr_n(throat);
    
    
    %% plotting the mass flow rate for various time steps
    figure(1)
    hold on 
    title(' distribution of mass flow rate over the length of nozzel at various time steps')
    if k == 1                                      % at time step = 1 
        plot(x,mfr_n,'--','linewidth',2)
    elseif k == 10                                 % at time step = 10
        plot(x,mfr_n,'b','linewidth',2)
    elseif k == 50                                 % st time step = 50
        plot(x,mfr_n,'r','linewidth',2)
    elseif k == 100                                % at time step = 100
        plot(x,mfr_n,'y','linewidth',2)
    elseif k == 200                                % at time step = 200
        plot(x,mfr_n,'c','linewidth',2)
    elseif k == 700                                % at time step = 700
        plot(x,mfr_n,'g','linewidth',2)
        legend('1','10','50','100','200','700')
        xlabel('nondimensional distance through nozzel')
        ylabel('mass flow rate')
    end
    
    
end


    %% plottting the timewise distribution of the variables at the throat
    figure(2)
    subplot(5,1,1)
    plot(rho_th_n)
    ylabel('density')
    title('timewise distribution of premitive variables at the throat')
    subplot(5,1,2)
    plot(v_th_n)
    ylabel('velocity')
    subplot(5,1,3)
    plot(t_th_n)
    ylabel('temperature')
    subplot(5,1,4)
    plot(p_th_n)
    ylabel('pressure')
    subplot(5,1,5)
    plot(M_th_n)
    ylabel('mach number')
    xlabel('number of time steps')
    
    %% ploting the distribution of variables over the length of the nozzel
    figure(3)
    subplot(5,1,1)
    plot(x,v_n)
    ylabel('velocity')
    title(' distribution of premitive variables over the length of nozzel')
    hold on
    subplot(5,1,2)
    plot(x,t_n)
    ylabel('temperature')
    hold on
    subplot(5,1,3)
    plot(x,rho_n)
    ylabel('density')
    hold on
    subplot(5,1,4)
    plot(x,p_n)
    ylabel('pressure')
    hold on
    subplot(5,1,5)
    plot(x,M_n)
    ylabel('Mach number')
    xlabel('X')
    

end

     

    Constervative Form:

    We will now see the conservative form of governing equations associated with the flow

    Continuity equation:

    ∂ρA∂t+∂ρAv∂x=0∂ρA∂t+∂ρAv∂x=0

    Momentum equation:

    ∂ρAv∂t+∂ρAv2+1γpA∂x=1γp∂A∂x∂ρAv∂t+∂ρAv2+1γpA∂x=1γp∂A∂x

    Energy equation: 

    ∂(ρ(eγ−1+γ2v2)A)∂t+∂(ρ(eγ−1+γ2v2)Av+pAv)∂x=0∂(ρ(eγ-1+γ2v2)A)∂t+∂(ρ(eγ-1+γ2v2)Av+pAv)∂x=0 

    all the above equations are in the form of non dimensioinalised variables.

    We will be converting our above equation in the form of solution vector , flux vectors and source terms

    solution vectors: 

    U1=ρAU1=ρA

    U2=ρAvU2=ρAv

    U3=ρ(eγ−1+γ2v2)AU3=ρ(eγ-1+γ2v2)A

    flux vectors:

    F1=ρAvF1=ρAv

    F2=ρAv2+1γpAF2=ρAv2+1γpA

    F3=ρ(eγ−1+γ2v2)vA+pAvF3=ρ(eγ-1+γ2v2)vA+pAv

    source terms:

    J2=1γp∂A∂xJ2=1γp∂A∂x

    now the after substituting the values for the solution vectors , flux vectors and source terms our governing equation will look like the following

    Continuity equation:

    ∂U1∂t=−∂F1∂x∂U1∂t=-∂F1∂x

    Momentum equation:

    ∂U2∂t=−∂F2∂x+J2∂U2∂t=-∂F2∂x+J2

    Energy equation:

    ∂U3∂t=−∂F3∂x∂U3∂t=-∂F3∂x

     

    we will use this above form of governing equation for our macormack method in our code 

    Now we will see the function code for our conservative form:

    function [ rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n, nt,x, gamma, CFL , a , dx, throat)

%% Function for the quasi 1D nozzel flow using conservative form


% calculate initial profile for conservative form
n_x_1 = ceil(0.5*(n-1)/3.5 + 1); % number of node at x = 0.5
n_x_2 = ceil(1.5*(n-1)/3.5 + 1); % number of node at x = 1.5
% density
rho_c(1:n_x_1) = 1.0;
rho_c(n_x_1:n_x_2) = 1.0 - 0.366*(x(n_x_1:n_x_2)-0.5);
rho_c(n_x_2:n) = 0.634 - 0.3879*(x(n_x_2:n)-1.5);
% temperature
t_c(1:n_x_1) = 1.0;
t_c(n_x_1:n_x_2) = 1.0 - 0.167*(x(n_x_1:n_x_2)-0.5);
t_c(n_x_2:n) = 0.833 - 0.3507*(x(n_x_2:n)-1.5);
% velocity
v_c = 0.59./(rho_c.*a);


% converting the premitive variables into solution variables
% solution vectors
U1 = rho_c.*a;                                        % continuity
U2 = rho_c.*a.*v_c;                                   % momentum
U3 = rho_c.*(t_c./(gamma-1) + gamma/2*v_c.^2).*a;     % energy

% time loop
for k = 1:nt
    
    % minimum time step value among all the nodes approach should be
    % employed. 
    dt = min(CFL*dx./(sqrt(t_c(2:n-1)) + v_c(2:n-1)));
        
    % saving the solution vectors
    U1_old = U1; 
    U2_old = U2;
    U3_old = U3;

    % flux vectors
    F1 = U2;
    F2 = U2.^2./U1 + (gamma -1)/gamma * (U3 - gamma/2 * U2.^2./U1);
    F3 = gamma*U2.*U3./U1 - gamma*(gamma-1)/2* U2.^3./U1.^2;

        
        
    %% perdictor method
    for i = 2:n-1
        
        % source term
        %J2 = (gamma -1)/gamma * (U3 - gamma/2 * U2.^2./U1) * dadx;
        J2_p(i) = (1/gamma)*rho_c(i)*t_c(i)*((a(i+1) - a(i))/dx);
        
        % predictor step
        % forward diferencing
        dU1dt_p(i) = - (F1(i+1) - F1(i))/dx;             % continuity
        dU2dt_p(i) = - (F2(i+1) - F2(i))/dx + J2_p(i);   % momentum
        dU3dt_p(i) = - (F3(i+1) - F3(i))/dx;             % energy
        
        % updating the solution variables 
        U1(i) = U1_old(i) + (dU1dt_p(i)*dt);
        U2(i) = U2_old(i) + (dU2dt_p(i)*dt);
        U3(i) = U3_old(i) + (dU3dt_p(i)*dt);
    end
    
    %% decoding the primitive variables
    rho_c_p = U1./a;
    v_c_p = U2./U1;
    t_c_p = (gamma -1)*(U3./U1 - (gamma/2) * ((v_c_p).^2));
    
    F1_c = U2;
    F2_c = U2.^2./U1 + (gamma -1)/gamma * (U3 - gamma/2 * U2.^2./U1);
    F3_c = gamma*U2.*U3./U1 - gamma*(gamma-1)/2* U2.^3./U1.^2;
    
    %% corrector method
    for i = 2:n-1
         
        % source terms
        J2_c(i) = (1/gamma)*rho_c_p(i)*t_c_p(i)*((a(i)-a(i-1))/dx);
        % corrector step
        % backward differencing
        dU1dt_c(i) = - (F1_c(i) - F1_c(i-1))/dx;                % continutiy
        dU2dt_c(i) = - (F2_c(i) - F2_c(i-1))/dx + J2_c(i);      % momentum 
        dU3dt_c(i) = - (F3_c(i) - F3_c(i-1))/dx;                % energy
        
        
    end
    
    % average time derivatives
    dU1dt_avg = 0.5*(dU1dt_p+ dU1dt_c);
    dU2dt_avg = 0.5*(dU2dt_p+ dU2dt_c);
    dU3dt_avg = 0.5*(dU3dt_p+ dU3dt_c);
    
    for i = 2:n-1
        % final updated  values 
        U1(i) = U1_old(i) + dU1dt_avg(i)*dt;      % continuity 
        U2(i) = U2_old(i) + dU2dt_avg(i)*dt;      % momentum
        U3(i) = U3_old(i) + dU3dt_avg(i)*dt;      % energy
        
    end
    
    %% apply boundary conditions
    % inlet
    U1(1)=rho_c_p(1).*a(1);
    U2(1) =2*U2(2) - U2(3);
    v_c_p(1) = U2(1)/U1(1);
    U3(1)= U1(1)*(t_c_p(1)/(gamma-1)+(gamma/2)*v_c_p(1).^2);
    
    % outlet
    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);
    
    %% final updated primitive variables 
    rho_c = U1./a ;
    v_c = U2./U1;
    t_c = (gamma - 1)* (U3/U1 - (gamma/2)*v_c.^2);
    % other premitive variables    
    p_c = rho_c.*t_c;            % pressure
    M_c = v_c./(t_c.^0.5);       % Mach number 
    mfr_c =rho_c.*v_c.*a;        % mass flow rate
    
    % values at throat
    rho_th_c(k) = rho_c(throat);
    v_th_c(k) = v_c(throat);
    t_th_c(k) = t_c(throat);
    M_th_c(k) = M_c(throat);
    p_th_c(k) = p_c(throat);
    mfr_th_c(k) = mfr_c(throat);
    
    %% plotting the mass flow rate for various time steps
    figure(4)
    hold on 
    title(' distribution of mass flow rate over the length of nozzel at various time steps')
    if k == 1
        plot(x,mfr_c,'--','linewidth',2)
    elseif k == 10
        plot(x,mfr_c,'b','linewidth',2)
    elseif k == 50
        plot(x,mfr_c,'r','linewidth',2)
    elseif k == 100
        plot(x,mfr_c,'y','linewidth',2)
    elseif k == 200
        plot(x,mfr_c,'c','linewidth',2)
    elseif k == 700
        plot(x,mfr_c,'g','linewidth',2)
        legend('1','10','50','100','200','700')
        xlabel('nondimensional distance through nozzel')
        ylabel('mass flow rate')
    end
    
    
    
end

    %% plottting the timewise distribution of the variables at the throat
    figure(5)
    subplot(5,1,1)
    plot(rho_th_c)
    ylabel('density')
    title('timewise distribution of primitive variables at the throat')
    subplot(5,1,2)
    plot(v_th_c)
    ylabel('velocity')
    subplot(5,1,3)
    plot(t_th_c)
    ylabel('temperature')
    subplot(5,1,4)
    plot(p_th_c)
    ylabel('pressure')
    subplot(5,1,5)
    plot(M_th_c)
    ylabel('mach number')
    xlabel('number of time steps')
    
    %% ploting the distribution of variables over the length of the nozzel
    figure(6)
    subplot(5,1,1)
    plot(x,v_c)
    ylabel('velocity')
    title(' distribution of premitive variables over the length of nozzel')
    hold on
    subplot(5,1,2)
    plot(x,rho_c)
    ylabel('density')
    hold on
    subplot(5,1,3)
    plot(x,t_c)
    ylabel('temperature')
    hold on
    subplot(5,1,4)
    plot(x,p_c)
    ylabel('pressure')
    hold on
    subplot(5,1,5)
    plot(x,M_c)
    ylabel('mach number')
    xlabel('X')
    

end

    Now using both the functions for Non-consesrvative and conservative forms in our main program we will see the results for a 1D nozzel geometry with 31 nodes.

     

    clear all
close all
clc

%inputs
n = 31; % number of nodes
x = linspace(0,3,n); 
dx = x(2)-x(1); % element length 
gamma = 1.4; 

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

%CFL criteria
CFL = 0.5;

% throat of the nozzel
throat = find(a==1);

% number of time steps
nt = 1400;

%% using the functions for non-conservative and conservative form of governing equations     
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);


%% plotting the mass flow rate comparison at the steady state for both forms     
figure(7)
plot(x,mfr_c,'linewidth',2)
hold on
plot(x,mfr_n,'linewidth',2)
legend('conservative','non-conservative')
xlabel('non dimensional length of nozzel')
ylabel('non dimensional mass flow rate')
title('comparison of normalised mass-flow-rate distributions for both forms')

    OUTPUTS:

    The results for the main program is given below:

    Mass flow rate at various time steps for non conservative form: 

    Time wise distribution of variables at the throat for non conservative form: 

    here we can see athat after 800 timesteps the values of the variables does not change.

    Distribution of properties across the length of the nozzel for non conservative form:

    Mass flow rate at various time steps for conservative form:

     

    Time wise distribution of variables at the throat for conservative form: 

    here we can see athat after 800 timesteps the values of the variables does not change.

    Distribution of properties across the length of the nozzel for conservative form:

    Also we will see the comparison of mass flow rate for both forms

    Comparison of mass flow rates:

     

    From the above plot we can see that the mass flow rate for the conservative form remains almost constant throughout the length of the nozzzel where as in the non conservative form the mass flow rate fluctuates more. in the Non-conservative forms the residuals and the round off error are very less as compared to the conservative forms. This is due to the conversion from the primitive variables to solution vectors, and back then to primitive variables in the conservative form.

    Grid Independence test:

    Now we will perform the grid independence test for both forms:

    we will now use 61 nodes in our simulation

    the code in the main program will have 61 as the number of nodes and all the other things will remian same for the code.

    similarly ,we will also use 91 grid point to perform the test.

    the code for grid indepndence test is given below:

    clear all
close all
clc

%inputs
for n = 31; % number of nodes
x = linspace(0,3,n); 
dx = x(2)-x(1); % element length 
gamma = 1.4; 

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

%CFL criteria
CFL = 0.5;

% throat of the nozzel
throat = find(a==1);

% number of time steps
nt = 1400;

%% using the functions for non-conservative and conservative form of governing equations     
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);

rho_th_n_gr(1) = rho_th_n(1400)
v_th_n_gr(1) = v_th_n(1400)
t_th_n_gr(1) = t_th_n(1400)

rho_th_c_gr(1) = rho_th_c(1400)
v_th_c_gr(1) = v_th_c(1400)
t_th_c_gr(1) = t_th_c(1400)

end


%inputs
for n = 61; % number of nodes
x = linspace(0,3,n); 
dx = x(2)-x(1); % element length 
gamma = 1.4; 

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

%CFL criteria
CFL = 0.5;

% throat of the nozzel
throat = find(a==1);

% number of time steps
nt = 1400;

%% using the functions for non-conservative and conservative form of governing equations     
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);

rho_th_n_gr(2) = rho_th_n(1400)
v_th_n_gr(2) = v_th_n(1400)
t_th_n_gr(2) = t_th_n(1400)

rho_th_c_gr(2) = rho_th_c(1400)
v_th_c_gr(2) = v_th_c(1400)
t_th_c_gr(2) = t_th_c(1400)

end


%inputs
for n = 91; % number of nodes
x = linspace(0,3,n); 
dx = x(2)-x(1); % element length 
gamma = 1.4; 

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

%CFL criteria
CFL = 0.5;

% throat of the nozzel
throat = find(a==1);

% number of time steps
nt = 1400;

%% using the functions for non-conservative and conservative form of governing equations     
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);

rho_th_n_gr(3) = rho_th_n(1400)
v_th_n_gr(3) = v_th_n(1400)
t_th_n_gr(3) = t_th_n(1400)

rho_th_c_gr(3) = rho_th_c(1400)
v_th_c_gr(3) = v_th_c(1400)
t_th_c_gr(3) = t_th_c(1400)

end

%% plotting the mass flow rate comparison at the steady state for both forms     
figure(7)
plot(x,mfr_c,'linewidth',2)
hold on
plot(x,mfr_n,'linewidth',2)
legend('conservative','non-conservative')
xlabel('non dimensional length of nozzel')
ylabel('non dimensional mass flow rate')
title('comparison of normalised mass-flow-rate distributions for both forms')

    From the simulation we see that the result

    we will first compare the results of the non conservative from

    now the results of the conservative form:

     

     from the results we can see that there is not much change in the values of the variables so our numerical simulation is grid independent and the grid is good for teh results.

    Also we can see that there is less change in the values of the non conservational form as compared to the comservational form. 

    Leave a comment

    Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.

    Please  login to add a comment

    Other comments...

    No comments yet!
    Be the first to add a comment

    Read more Projects by Amol Patel (31)

    Week 6: Conjugate Heat Transfer Simulation

    Objective:

    AIM-  To simulate a Conjugate Heat Transfer flow through a pipe, while the inlet Reynolds number should be 7000. To run the grid independance test on 3 grids and show that the outlet temperature converges to a particular value To observe the effect of various supercycle stage interval on the total simulation time.…

    calendar

    09 Nov 2022 06:55 AM IST

      Read more

      Week 7: Shock tube simulation project

      Objective:

      AIM - To set up a transient shock tube simulation Plot the pressure and temperature history in the entire domain Plot the cell count as a function of time    SHOCK TUBE- The shock tube is an instrument used to replicate and direct blast waves at a sensor or a model in order to simulate actual explosions…

      calendar

      07 Nov 2022 09:18 PM IST

        Read more

        Week 5: Prandtl Meyer Shock problem

        Objective:

        AIM - 1. To understand what is a shock wave. 2. To understand the what are the different boundary conditions used for a shock flow problems. 3. To understand the effect of SGS parameter on shock location. 4. To simulate Prandalt Meyer Shcok Wave.   OBJECTIVE -  Que 1. What is Shock Wave? A shock wave or shock,…

        calendar

        01 Nov 2022 06:36 PM IST

          Read more

          Week 4.2: Project - Transient simulation of flow over a throttle body

          Objective:

          AIM - Transient simulation of flow over a throttle body.   OBJECTIVE -  Setup and run transient state simulation for flow over a throttle body. Post process the results and show pressure and velocity contours. Show the mesh (i.e surface with edges) Show the plots for pressure, velocity, mass flow rate and total…

          calendar

          12 Feb 2022 07:08 AM IST

            Read more

            Schedule a counselling session

            Please enter your name
            Please enter a valid email
            Please enter a valid number

            Related Courses

            coursecard

            Design loads considered on bridges

            Recently launched

            10 Hours of Content

            coursecard

            Design of Steel Superstructure in Bridges

            Recently launched

            16 Hours of Content

            coursecard

            Design for Manufacturability (DFM)

            Recently launched

            11 Hours of Content

            coursecard

            CATIA for Medical Product Design

            Recently launched

            5 Hours of Content

            coursecardcoursetype

            Accelerated Career Program in Embedded Systems (On-Campus) Courseware Partner: IT-ITes SSC nasscom

            Recently launched

            0 Hours of Content

            Schedule a counselling session

            Please enter your name
            Please enter a valid email
            Please enter a valid number

            logo

            Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.

            https://d27yxarlh48w6q.cloudfront.net/web/v1/images/facebook.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/insta.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/twitter.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/youtube.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/linkedin.svg

            Our Company

            News & EventsBlogCareersGrievance RedressalSkill-Lync ReviewsTermsPrivacy PolicyBecome an Affiliate
            map
            EpowerX Learning Technologies Pvt Ltd.
            4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
            mail
            info@skill-lync.com
            mail
            ITgrievance@skill-lync.com

            Top Individual Courses

            Computational Combustion Using Python and CanteraIntroduction to Physical Modeling using SimscapeIntroduction to Structural Analysis using ANSYS WorkbenchIntroduction to Structural Analysis using ANSYS Workbench

            Top PG Programs

            Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in Electric Vehicle Design & Development

            Skill-Lync Plus

            Executive Program in Electric Vehicle Embedded SoftwareExecutive Program in Electric Vehicle DesignExecutive Program in Cybersecurity

            Trending Blogs

            Heat Transfer Principles in Energy-Efficient Refrigerators and Air Conditioners Advanced Modeling and Result Visualization in Simscape Exploring Simulink and Library Browser in Simscape Advanced Simulink Tools and Libraries in SimscapeExploring Simulink Basics in Simscape

            © 2025 Skill-Lync Inc. All Rights Reserved.

                        Do You Want To Showcase Your Technical Skills?
                        Sign-Up for our projects.