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

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

Schematic representation of convergent divergent nozzle.

In the above described convergent divergent nozzle, while simulating for the flow parameters, an assumptions has been made that the flow property will only vary with x. Thus the uniform flow property along cross section is assumed. The assumption made here, makes the simulation, quasi-one-Dimensional flow simulation.

 

 

The physical phenomenon of the stated problem is governed by the following equations:-

The equations are presented in 2 forms, Conservative and Non-Conservative Form.

Conservative Form

Continuity Equation:-

`(del(rho A))/(del t) + (del(rho A V))/(del x) = 0`

Momentum Equation:-

`(del(rho A V))/(del t) + (del(rho A V^2))/(del x) = -A (del p)/(del x)`

Energy Equation:-

`(del[rho (e + V^2/2) A])/(del t) + (del[rho (e + V^2/2) VA])/(del x) = - (del(p A V))/(del x)`

 

Non Conservation Form 

Continuity Equation:-

`A (del rho)/(del t) + V (del(rho A))/(del x) = -rho A (del V)/(del x)`

Momentum Equation:-

`rho (del V)/(del t) + rho V(del V)/(del x) = -(del p)/(del x)`

Energy Equation:-

`rho (del e)/(del t) + rho V(del e)/(del x) = - p(del V)/(del x) - pV(del(ln A))/(del x)`

 

 

 

 

In order to non dimensionalise the variables, we will be using reservoir values as the base. 

Non-Dimensional values will be represented with ' superscript

Velocity

`V' = V/a_0`

Temperature

`T' = T/T_0`

Density

`rho' = rho/rho^0`

Length

`x' = x/L`

time

`t' = {t}/{{L}/{a_0}}`

Pressure

`P' = P/P^0`

Area

`A' = A/A^0`

Where, 

`A^0` is the throat area.

`a_0` is `sqrt(gammaRT_0)` represents speed of sound in the reservoir.

`P^0,rho^0,T^0` is the reservoir pressure, density and temperature.

Maccormack Method

It is a 2 step method for discretization of hyperbolic partial differential equations. The method is second order accurate in time and space domain. The first step called predictor step is used to find so called predicted value for the next time step. Then, in corrector step, the predicted value for next time step is corrected using backward differencing in space domain.

Governing Equations:-

1. Conservation Form using dimentionless variables:-

Continuity Equation:-

`{del(rho'A')}/{delt'} + {del(rho'A'V')}/{dx'} = 0`

Momentum Equation:-

`{del(rho'A'V')}/{dt'} + {del[rho'A'V^('2) + ({1}/{gamma})p'A']}/{delx'} = {1}/{gamma}p'{delA'}/{delx'}`

Energy Equation :-

`{del[rho'({e'}/{gamma-1}+{gamma}/{2}V^('2))A']}/{delt'} +{del[rho'({e'}/{gamma - 1} +{gamma}/{2}V^('2))V'A' + p'A'V']}/{delx'} = 0`

 

Dependent variables are not the primitive variables while solving conservation form equations, So they need to be converted into solution vectors for solving and later they need to be decoded into primitive variables after solving.

 

Defining Source(J), Flux(F) and Solution Vectors(U) for the solution of above equation we get,

`U_1 = rho'A'`

`U_2 = rho'A'V'`

`U_3 = rho'({e'}/{gamma -1} + {gamma}/{2}V^('2))A'`

`F_1 = rho'A'V'`

`F_2 = rho'A'V^('2) + {1}/{gamma}p'A'`

`F_3 = rho'({e'}/{gamma - 1} + {gamma}/{2}V^('2))V'A' +p'A'V'`

`J_2 = {1}/{gamma}p'(delA')`

 

time deravatives of solution vectors:-

`(del U_1)/(del t') = -(del F_1)/(del x')`

`(del U_2)/(del t') = J2 - (del F_2)/(del x')`

`(del U_3)/(del t') = -(del F_3)/(del x')`

 

The premitive variables can be decoded back from solution vectors using below expressions:-

`rho' = {U'}/{A'}`

`V' = {U_2}/{U_1}`

`T' = e' = (gamma - 1)({U_3}/{U_1} - {gamma}/{2}V^('2))`

`p' = rho'T'`

 

Transforming the governing equation entirely in the form of solution vecors

`F_1 = U_2`

`F_2 = U_2^2/U_1 + (gamma-1)/gamma(U_3 - gamma/2 U_2^2/U_1)`

`F_3 = gamma (U_2 U_3)/U_1 - (gamma(gamma-1))/2 U_2^3/U_1^2`

`J_2 = (gamma-1)/gamma(U_3 - gamma/2 U_2^2/U_1)(del(ln A'))/(del x')`

 

Predictor Step :-

Calculating time deravative of the solution vectors:-

`((del U_1)/(del t'))_i^t = -(F_(1(i+1))^t - F_(1(i))^t)/(Delta x')`

`((del U_2)/(del t'))_i^t = (J_(2(i+1))^t-J_(2(i))^t)/(Delta x') - (F_(2(i+1))^t - F_(2(i))^t)/(Delta x')`

`((del U_3)/(del t'))_i^t = -(F_(3(i+1))^t - F_(3(i))^t)/(Delta x')`

 

Claculating the predictor values

`bar(U)_(1(i))^(t+Delta t) = ((del U_1)/(del t'))_i^t .Delta t'`

`bar(U)_(2(i))^(t+Delta t) = ((del U_2)/(del t'))_i^t .Delta t'`

`bar(U)_(3(i))^(t+Delta t) = ((del U_3)/(del t'))_i^t .Delta t'`

 

Corrector Step:-

Calculating time deravitive:-

 

`(bar(del U_1)/(del t'))_i^(t+Delta t) = -(bar(F)_(1(i))^(t+Delta t) - bar(F)_(1(i-1))^(t+Delta t))/(Delta x')`

`(bar(del U_2)/(del t'))_i^(t+Delta t) = (bar(J)_(2(i))^(t+Delta t)-bar(J)_(2(i-1))^(t+Delta t))/(Delta x') - (bar(F)_(2(i))^(t+Delta t) - bar(F)_(2(i-1))^(t+Delta t))/(Delta x')`

`(bar(del U_3)/(del t'))_i^(t+Delta t) = -(bar(F)_(3(i))^(t+Delta t) - bar(F)_(3(i-1))^(t+Delta t))/(Delta x')`

 

Claculating average deravative:-

`((del U_1)/(del t'))_(av)=1/2[((del U_1)/(del t'))_i^t+(bar(del U_1)/(del t'))_i^(t+Delta t)]`

`((del U_2)/(del t'))_(av)=1/2[((del U_2)/(del t'))_i^t+(bar(del U_2)/(del t'))_i^(t+Delta t)]`

`((del U_3)/(del t'))_(av)=1/2[((del U_3)/(del t'))_i^t+(bar(del U_3)/(del t'))_i^(t+Delta t)]`

 

Calculating corrected values:-

`U_(1(i))^(t+Delta t) = ((del U_1)/(del t'))_(av) .Delta t'`

`U_(2(i))^(t+Delta t) = ((del U_2)/(del t'))_(av).Delta t'`

`U_(3(i))^(t+Delta t) = ((del U_3)/(del t'))_(av).Delta t'`

 

After the above step the solution vectors are decoded into non dimensional primitive values.

 

2. Non Conservation Form 

Continuity Equation:-

`{delrho'}/{delt'} = -rho'{delV'}/{delx'} - rho'V'{del(lnA')}/{delx'} - V'{delrho'}/{delx'}`

 

Momentum Equation:-

`{delV'}/{delt'} = -V'{delV'}/{delx'} - {1}/{gamma}({delT'}/{delx'} + {T'}/{rho'}{delrho'}/{delx'})`

 

Energy Equation:-

`{delT'}/{delt'} = -V'{delT'}/{delx'} - (gamma-1)T'[{delV'}/{delx'} + V'{del(lnA')}/{delx'}]`

 

 

Maccormak Strategy for Non-Conservation part:-

 

Calculating time deravatives for current time step, to be used to calculate values for predictor step.

`({delrho}/{delt})_i^t = -rho_i^t[{V_(t+1)^t - V_i^t}/{Deltax}] - rho_i^tV_i^t[{ln(A_(i+1)) - ln(A_i)}/{Deltax}] - V_i^t[{rho_(i+1)^t - rho_i^t}/{Deltax}]`

 

`((del V)/(del t))_i^t = -V_i^t[(V_(i+1)^t - V_i^t)/(Delta x)] - 1/gamma[(T_(i+1)^t - T_i^t)/(Delta x) + T_i^t/rho_i^t( (rho_(i+1)^t - rho_i^t)/(Delta x))]`

 

`((del T)/(del x))_i^t = -V_i^t[(T_(i+1)^t - T_i^t)/(Delta x)] - (gamma-1)T_i^t [(V_(i+1)^t-V_i^t)/(Delta x)+ V_i^t((ln(A_(i+1)) - ln(A_i))/(Delta x))]`

 

Computing predictor values for rho, velocity and temperature.

`(bar rho)_i^(t+Delta t) = rho_i^t + ((del rho)/(del t))_i^t Delta t`

`(bar V)_i^(t+Delta t) = V_i^t + ((del V)/(del t))_i^t Delta t`

`(bar T)_i^(t+Delta t) = T_i^t + ((del T)/(del t))_i^t Delta t`

Calculating deravatives for corrector steps:-

`(bar((del rho)/(del t)))_i^(t+Delta t) = -(bar rho)_i^(t+Delta t)[((bar V)_i^(t+ Delta t) - (bar V)_(i-1)^(t+Delta t))/(Delta x)] - (bar rho)_i^(t+Delta t)(bar V)_i^(t+Delta t)[(ln(A_i) - ln(A_(i-1)))/(Delta x)] - (bar V)_i^(t+Delta t)[((bar rho)_i^(t+Delta t) - (bar rho)_(i-1)^(t+Delta t))/(Delta x)]`

 

`(bar((del V)/(del t)))_i^(t+ Delta t) = -(bar V)_i^(t+Delta t)[((bar V)_i^(t+Delta t) - (bar V)_(i-1)^(t+Delta t))/(Delta x)] - 1/gamma[((bar T)_i^(t+Delta t) - (bar T)_(i-1)^(t+Delta t))/(Delta x) + (bar T)_i^(t+Delta t)/(bar rho)_i^(t+Delta t)( ((bar rho)_i^(t+Delta t) - (bar rho)_(i-1)^(t+Delta t))/(Delta x))]`

 

`(bar((del T)/(del x)))_i^(t+Delta t) = -(bar V)_i^(t+Delta t)[((bar T)_i^(t+Delta t) - (bar T)_(i-1)^(t+Delta t))/(Delta x)] - (gamma-1)T_i^(t+Delta t) [((bar V)_i^(t+Delta t)- (bar V)_(i-1)^(t+Delta t))/(Delta x)+ (bar V)_i^(t+Delta t)((ln(A_i) - ln(A_(i-1)))/(Delta x))]`

 

Calculating average of deravatives:-

`((del rho)/(del t))_(av) = 1/2[((del rho)/(del t))_i^t + (bar((del rho)/(del t)))_i^(t+Delta t)]`

`((del V)/(del t))_(av) = 1/2[((del V)/(del t))_i^t + (bar((del V)/(del t)))_i^(t+Delta t)]`

`((del T)/(del t))_(av) = 1/2[((del T)/(del t))_i^t + (bar((del T)/(del t)))_i^(t+Delta t)]`

 

Calculating the final values for next time step:-

`rho_i^(t+Delta t) = rho_i^t + ((del rho)/(del t))_(av) Delta t`

`V_i^(t+Delta t) = V_i^t + ((del V)/(del t))_(av) Delta t`

`T_i^(t+Delta t) = T_i^t + ((del T)/(del t))_(av) Delta t`

`p_i^(t+Delta t) = rho_i^(t+Delta t) .T_i^(t+Delta t)`

 

The Time Step Calculation of both the scheme has been discribed below:-

In order for the solution to be stable the effect of courant no. as been taken into consideration.

`Deltat = C{Deltax}/{a+V}`

where,

C is the courant No. As in the literatue, for stable analysis, C is taken to be 0.5.

`Deltat`for each forward step is calculated for each grid point at once using the previous time step values of a, v. Then the minimum `Deltat` is choosen for the present time step calculation. This has been imlemented in the code using a seperate function called updateDt.

 The Total time steps taken as 1400. The value has been derived from the below refrenced literature, which suggests, it to be way beyond the steady state.

The initial and boundary Condition has been used as prescribed in the Anderson, P. J. D., Anderson, J. D. (1995). Computational Fluid Dynamics. Colombia: McGraw-Hill Education.

 

The program developed to study the above case is as follows:-

 The program developed is organized in classes and a main driver script.

It contains, class initialCondition for storing the initial class for the simulation, quasiOnedim for storing the simulation functions and parameters, Result class to store the result of the simulation.

 The quasiOnedim class:-

classdef quasiOnedim
    %Class implementing Maccormack's method for quasi-one-Dim convergent-divergent nozzle 
    %   Implementing the class requires following parameters:-
    %   n = No. of Node Points
    %   L = Length of the domain
    %   gamma = Specific Heat
    %   initObj = initialCondition object
    %   C = Courant No.
    %   nt = Total time steps
    
    properties
        n
        L
        % x-axis grid
        x
        % x-axis grid size
        dx
        % calculated time step
        dt
        % Area of the given profile
        A
        gamma
        % Values of the non dimensional variables at time t = 0
        % Non Dimentional density
        rho
        % Non - Dimentional Velocity
        v
        % Non - Dimentional temperature
        T
        
        C
        nt

        
        
    end
    
    methods
        function obj = quasiOnedim(n,L,gamma,initObj,C,nt)
            %Initializes the declared properties            
            obj.n = n;
            obj.L = L;
            obj.x = linspace(0,L,n);
            obj.dx = obj.x(2) - obj.x(1);
            obj.A = 1 + 2.2*(obj.x - 1.5).^2;
            obj.gamma = gamma;
            obj.nt = nt;
            
            obj.rho = initObj.rho;
            obj.v = initObj.v;
            obj.T = initObj.T;
            
            obj.C = C;       
        end
        
        function obj = updateDt(obj)
            %Update the time step using scheme discussed in the theory
            %section
            a = sqrt(obj.T);
            delt = (obj.C * obj.dx) ./ (a + obj.v);
            dt_temp = min(delt);
            obj.dt = dt_temp;
        end
        
        function [obj,result] = conservativeCD(obj)
            % Implementing Conservative Governing Equation using
            % Maccormak's Method.
            
            % Declaring the result variable for saving transient result
            transient_result = ones(obj.nt,obj.n,7);
            dt_simulation = ones(obj.nt,1);
            % Initializing result object
            result = Result(1,1,1,1,1,1,1);
            % Solution vectors declaration.
            u1 = obj.rho .* obj.A;
            u2 = obj.rho .* obj.A .* obj.v;
            u3 = obj.rho .* (obj.T/(obj.gamma-1) + (obj.gamma/2) * (obj.v.^2)) .* obj.A;
            % Time loop
                for z = 1:obj.nt 
                    % updating time step size for stability
                    obj = updateDt(obj);
                    % Storing the initial values as a previous time step
                    % data
                    u1Prev = u1;
                    u2Prev = u2;
                    u3Prev = u3;
                    
                    % Declaring the Flux Vectors
                    f1 = u2;
                    f2 = (u2.^2 ./ u1) + ((obj.gamma-1)/obj.gamma) * (u3 - (obj.gamma/2) * (u2.^2 ./ u1));
                    f3 = (obj.gamma * ((u2 .* u3) ./ u1)) - (obj.gamma*(obj.gamma-1)/2) * (u2.^3 ./ u1.^2);
                    
                    
                    % Predictor Step claculation as discussed int the theory
                    for i = 2: obj.n-1
                        % Deravative of the flux term 
                        df1Dx = (f1(i+1) - f1(i)) / obj.dx;
                        df2Dx = (f2(i+1) - f2(i)) / obj.dx;
                        df3Dx = (f3(i+1) - f3(i)) / obj.dx;
                        daDx = (obj.A(i+1) - obj.A(i)) / obj.dx;
                        
                        % Source Term  calculation
                        J2(i) = (1/obj.gamma) * obj.rho(i) * obj.T(i) * daDx;
                        
                        % Deravative of Solution vectors
                        du1DtPred(i) = -df1Dx;
                        du2DtPred(i) = -df2Dx + J2(i);
                        du3DtPred(i) = -df3Dx;
                        
                        
                        % Updating the predictor step values
                        u1(i) = u1(i) + du1DtPred(i) * obj.dt;
                        u2(i) = u2(i) + du2DtPred(i) * obj.dt;
                        u3(i) = u3(i) + du3DtPred(i) * obj.dt;
                        
                    end
                    
                    
                    % Update Flux Vectors
                    
                    f1 = u2;
                    f2 = (u2.^2 ./ u1) + ((obj.gamma-1)/obj.gamma) * (u3 - (obj.gamma/2) * (u2.^2 ./ u1));
                    f3 = (obj.gamma * ((u2 .* u3) ./ u1)) - (obj.gamma*(obj.gamma-1)/2) * (u2.^3 ./ u1.^2);
                    
                    
                    % Decoding the variables 
                    
                    obj.rho = u1 ./ obj.A;
                    obj.v = u2 ./ u1;
                    obj.T = (obj.gamma - 1) * (u3 ./ (obj.rho .* obj.A) - (obj.gamma/2) * obj.v.^2);
                    
                    % Corrector step Calculation
                    
                    for i = 2: obj.n-1
                        % Deravative of updated flux 
                        df1Dx = (f1(i) - f1(i-1)) / obj.dx;
                        df2Dx = (f2(i) - f2(i-1)) / obj.dx;
                        df3Dx = (f3(i) - f3(i-1)) / obj.dx;
                        daDx = (obj.A(i) - obj.A(i-1)) / obj.dx;
                        
                        
                        J2(i) = (1/obj.gamma) * obj.rho(i) * obj.T(i) * daDx;
                        du1DtCorrec(i) = -df1Dx;
                        du2DtCorrec(i) = -df2Dx + J2(i);
                        du3DtCorrec(i) = -df3Dx;
                        
                    end
                    
                    % Calculating the Average Derivative
                    
                    du1Dt = 0.5 * (du1DtPred + du1DtCorrec);
                    du2Dt = 0.5 * (du2DtPred + du2DtCorrec);
                    du3Dt = 0.5 * (du3DtPred + du3DtCorrec);
                    
                    
                    % Updating the solution vectors
                    
                    for j = 2: obj.n-1
                        
                        u1(j) = u1Prev(j) + du1Dt(j) * obj.dt;
                        u2(j) = u2Prev(j) + du2Dt(j) * obj.dt;
                        u3(j) = u3Prev(j) + du3Dt(j) * obj.dt;
                        
                    end
                    
                    % Applying Boundary Conditions
                    
                    % Inlet BC
                    obj.rho(1) = 1;
                    obj.T(1) = 1;
                    
                    u1(1) = obj.rho(1) * obj.A(1);
                    u2(1) = 2 * u2(2) - u2(3);
                    u3(1) = u1(1) * (obj.T(1)/(obj.gamma-1) + (obj.gamma/2) * ((u2(1) / u1(1))^2));
                    
                    % Outlet BC
                    u1(obj.n) = 2 * u1(obj.n-1) - u1(obj.n-2);
                    u2(obj.n) = 2 * u2(obj.n-1) - u2(obj.n-2);
                    u3(obj.n) = 2 * u3(obj.n-1) - u3(obj.n-2);
                    
                    % Decoding the properties using solution vectors
                    obj.rho = u1 ./ obj.A;
                    obj.v = u2 ./ u1;
                    obj.T = (obj.gamma - 1) * (u3 ./ (obj.rho .* obj.A) - (obj.gamma/2) * obj.v.^2);
                    P = obj.rho .* obj.T;
                    M = obj.v ./ sqrt(obj.T);
                    massFlow = obj.rho .* obj.v .* obj.A;

                    %Saving time series data
                    transient_result(z,:,1) = obj.rho;
                    transient_result(z,:,2) = obj.T;
                    transient_result(z,:,3) = obj.v;
                    transient_result(z,:,4) = P;
                    transient_result(z,:,5) = M;
                    transient_result(z,:,6) = massFlow;
                    dt_simulation(z) = obj.dt;

                end
 
            % Updating the declared Result object
            result.rho = transient_result(:,:,1);
            result.T = transient_result(:,:,2);
            result.v = transient_result(:,:,3);
            result.P = transient_result(:,:,4);
            result.M = transient_result(:,:,5);
            result.massFlow = transient_result(:,:,6);
            result.dt = dt_simulation;
            
            
        end
               
        function [obj,result] = nonConservativeCD(obj)
            % Implements nonconservative form for the simulation
            
            % Initializing the result variable for storing the transient
            % results
            transient_result = ones(obj.nt,obj.n,6);
            dt_simulation = ones(obj.nt,1);
            % Initializing the results object
            result = Result(1,1,1,1,1,1,1);
            for z = 1:obj.nt
               % Storing the inital values to be used as previous time step
               % data
                rhoOld = obj.rho;
                vOld = obj.v;
                TOld = obj.T;
                
                % Time step calculation based on Stability
                obj = updateDt(obj);
                
                % Predictor Step claculation as discussed int the theory
                
                for i = 2: obj.n-1
                    % Calculating the deravitives
                    dRhoDx = (obj.rho(i+1) - obj.rho(i)) / obj.dx;
                    dvDx = (obj.v(i+1) - obj.v(i)) / obj.dx;
                    dtDx = (obj.T(i+1) - obj.T(i)) / obj.dx;
                    dlogaDx = (log(obj.A(i+1)) - log(obj.A(i))) / obj.dx;
                    
                    % Continuity, Momentum and Energy Equation
                    dRhoDtPred(i) = -(obj.rho(i) * dvDx) - (obj.rho(i) * obj.v(i) * dlogaDx) - (obj.v(i) * dRhoDx);
                    dvDtPred(i) = -(obj.v(i) * dvDx) - (1/obj.gamma) * (dtDx + (obj.T(i)/obj.rho(i)) * dRhoDx);
                    dTdtPred(i) = -(obj.v(i) * dtDx) - (obj.gamma - 1) * obj.T(i) * (dvDx + obj.v(i) * dlogaDx);
                    
                    % Updating the predicted values
                    obj.rho(i) = obj.rho(i) + dRhoDtPred(i) * obj.dt;
                    obj.v(i) = obj.v(i) + dvDtPred(i) * obj.dt;
                    obj.T(i) = obj.T(i) + dTdtPred(i) * obj.dt;
                    
                end
                
                % Calculation of Corrector Step
                
                for i = 2: obj.n-1
                    % Coalculating the deravitives
                    dRhoDx = (obj.rho(i) - obj.rho(i-1)) / obj.dx;
                    dvDx = (obj.v(i) - obj.v(i-1)) / obj.dx;
                    dtDx = (obj.T(i) - obj.T(i-1)) / obj.dx;
                    dlogaDx = (log(obj.A(i)) - log(obj.A(i-1))) / obj.dx;
                    
                    % Continuity, Momentum and Energy Equation calculation
                    dRhoDtCorrect(i) = -(obj.rho(i) * dvDx) - (obj.rho(i) * obj.v(i) * dlogaDx) - (obj.v(i) * dRhoDx);
                    dvDtCorrect(i) = -(obj.v(i) * dvDx) - (1/obj.gamma) * (dtDx + (obj.T(i)/obj.rho(i)) * dRhoDx);
                    dTDtCorrect(i) = -(obj.v(i) * dtDx) - (obj.gamma - 1) * obj.T(i) * (dvDx + obj.v(i) * dlogaDx);
                    
                end
                
                % Calculating the Average Derivative
                
                dRhoDtAvg = 0.5 * (dRhoDtPred + dRhoDtCorrect);
                dvDtAvg = 0.5 * (dvDtPred + dvDtCorrect);
                dTDtAvg = 0.5 * (dTdtPred + dTDtCorrect);
                
                % Updating the final values in the variables
                
                for j = 2: obj.n-1
                    
                    obj.rho(j) = rhoOld(j) + dRhoDtAvg(j) * obj.dt;
                    obj.v(j) = vOld(j) + dvDtAvg(j) * obj.dt;
                    obj.T(j) = TOld(j) + dTDtAvg(j) * obj.dt;
                    
                end
                
                % Apply=ing Boundary Conditions
                
                % Inlet BC
                obj.v(1) = 2 * obj.v(2) - obj.v(3);
                % Outlet BC
                obj.rho(obj.n) = 2 * obj.rho(obj.n-1) - obj.rho(obj.n-2);
                obj.v(obj.n) = 2 * obj.v(obj.n-1) - obj.v(obj.n-2);
                obj.T(obj.n) = 2 * obj.T(obj.n-1) - obj.T(obj.n-2);
                
                P = obj.rho .* obj.T;
                M = obj.v ./ sqrt(obj.T);
                massFlow = obj.rho .* obj.v .* obj.A;
        
                % Updating the transient values
                transient_result(z,:,1) = obj.rho;
                transient_result(z,:,2) = obj.T;
                transient_result(z,:,3) = obj.v;
                transient_result(z,:,4) = P;
                transient_result(z,:,5) = M;
                transient_result(z,:,6) = massFlow;
                dt_simulation(z) = obj.dt;


            end
            % Updating the result object
            result.rho = transient_result(:,:,1);
            result.T = transient_result(:,:,2);
            result.v = transient_result(:,:,3);
            result.P = transient_result(:,:,4);
            result.M = transient_result(:,:,5);
            result.massFlow = transient_result(:,:,6);
            result.dt = dt_simulation;
        end   
    end
end

the initialCondition class

classdef initialCondition
    %The class stores the initial condition for the numerical analysis
    %   It has following parameters:
    %   rho = non-dim- density
    %   T = non-dim Temperature
    %   v = non-dim Velocity
    
    properties
        rho
        T
        v
    end
    
    methods
        function obj = initialCondition(rho,T,v)
            %Implements the initial condition and initializes the object
            
            obj.rho = rho;
            obj.T = T;
            obj.v = v;
        end
        
    end
end

The Result class

classdef Result
    %Defines the Result class for storing time-series simulation data
    % It has following parameters:-
    % rho - non-dimensional-density
    % T = non-dimensiona-Temperature
    % v = non-dimensiona-velocity
    % P = non-dimensiona-Pressure
    % M = Mach No.
    % massFlow = non dimensional mass flow
    % dt = time step size
    
    properties
        rho
        T
        v
        P
        M
        massFlow
        dt
    end
    
    methods
        function obj = Result(rho,T,v,P,M,massFlow,dt)
            %Initializes the result variable
            
            obj.rho = rho;
            obj.T = T;
            obj.v = v;
            obj.P = P;
            obj.M = M;
            obj.massFlow = massFlow;
            obj.dt = dt;
        end
        
    end
end

The main driver script

close all
clear all
clc

n = 31;
gamma = 1.4;
C = 0.5;
totalTimeSteps = 1400;
%% Domain declaration
L = 3;
x = linspace(0,L,n);
A = 1 + 2.2*(x - 1.5).^2;
%% initial conditions
% Initial Temperature and velocity values in domain at time 0.
% Conservative 
Rho = ones(1,n);
T = ones(1,n);
for m = 1:n
    if x(m) >= 0 && x(m) < 0.5
        Rho(m) = 1;
        T(m) = 1;
    end
    
    if x(m) >= 0.5 && x(m) < 1.5
        Rho(m) = 1 - 0.366 * (x(m) - 0.5);
        T(m) = 1 - 0.167 * (x(m) - 0.5);
    end
    
    if x(m) >= 1.5 && x(m) <= 3
        Rho(m) = 0.634 - 0.3879 * (x(m) - 1.5);
        T(m) = 0.833 - 0.3507 * (x(m) - 1.5);
    end
        
end
v = 0.59 ./ (Rho .* A);

%% Simulation Setup
% Conservative
conservative_initialCond = initialCondition(Rho,T,v);
conservative = quasiOnedim(n,L,gamma,conservative_initialCond,C,totalTimeSteps);
[conservative,result_conservative] = conservativeCD(conservative);

% Non - Conservative
Rho = 1 - 0.3146*x;
T = 1 - 0.2314*x;
v = (0.1 + 1.09*x) .* T.^0.5;
nonconservative_initialCond = initialCondition(Rho,T,v);
nonconservative = quasiOnedim(n,L,gamma,nonconservative_initialCond,C,totalTimeSteps);
[nonconservative,result_nonconservative] = nonConservativeCD(nonconservative);


%% Plotting

close all
%% Plotting
timeStepAxis = linspace(0,totalTimeSteps,totalTimeSteps);
throat_index = find(A==min(A));
% 1. Steady-state distribution of primitive variables inside the nozzle
figure(1)
subplot(1,2,1)
plot(x,result_conservative.rho(totalTimeSteps,:), 'linewidth',2, 'color', 'b')
hold on
plot(x,result_conservative.v(totalTimeSteps,:), 'linewidth',2, 'color', 'g')
hold on
plot(x,result_conservative.T(totalTimeSteps,:), 'linewidth',2, 'color', 'r')
hold on
plot(x,result_conservative.P(totalTimeSteps,:), 'linewidth',2, 'color', 'c')
xlabel("X Domain")
ylabel("Non Dimension Magnitude")
title("Conservative Form")
legend("rho", 'velocity','Temperature','Pressure','location','best')


subplot(1,2,2)
plot(x,result_nonconservative.rho(totalTimeSteps,:), 'linewidth', 2, 'color', 'b')
hold on
plot(x,result_nonconservative.v(totalTimeSteps,:), 'linewidth',2, 'color', 'g')
hold on
plot(x,result_nonconservative.T(totalTimeSteps,:), 'linewidth',2, 'color', 'r')
hold on
plot(x,result_nonconservative.P(totalTimeSteps,:), 'linewidth',2, 'color', 'c')
xlabel("X Domain")
ylabel("Non Dimensional Magnitude")
title("Non - Conservative Form")
sgtitle("Steady-State Distribution of Non Dimensional Primitive Variables")
legend("rho", 'velocity','Temperature','Pressure','location','best')

% 2. Time-wise variation of the primitive variables
figure(2)
subplot(1,2,1)
plot(timeStepAxis,result_conservative.rho(:,throat_index), 'linewidth',2, 'color', 'b')
hold on
plot(timeStepAxis,result_conservative.v(:,throat_index), 'linewidth',2, 'color', 'g')
hold on
plot(timeStepAxis,result_conservative.T(:,throat_index), 'linewidth',2, 'color', 'r')
hold on
plot(timeStepAxis,result_conservative.P(:,throat_index), 'linewidth',2, 'color', 'c')
hold on
plot(timeStepAxis,result_conservative.M(:,throat_index), 'linewidth',2, 'color', 'm')
xlabel("Time Steps")
ylabel("Non Dimensional Magnitude")
title("Conservative Form")
legend("rho","velocity","Temperature","Pressure","Mach No.")
xlim([0,totalTimeSteps])
ylim([0, 2])


subplot(1,2,2)
plot(timeStepAxis,result_nonconservative.rho(:,throat_index), 'linewidth',2, 'color', 'b')
hold on
plot(timeStepAxis,result_nonconservative.v(:,throat_index), 'linewidth',2, 'color', 'g')
hold on
plot(timeStepAxis,result_nonconservative.T(:,throat_index), 'linewidth',2, 'color', 'r')
hold on
plot(timeStepAxis,result_nonconservative.P(:,throat_index), 'linewidth',2, 'color', 'c')
hold on
plot(timeStepAxis,result_nonconservative.M(:,throat_index), 'linewidth',2, 'color', 'm')
xlabel("Time Steps")
ylabel("Nin Dimensional Magnitude")
title("Non - Conservative Form")
sgtitle("Time-wise variation of the Non Dimensional variables at throat section")
legend("rho","velocity","Temperature","Pressure","Mach No.",'location','best')
xlim([0,totalTimeSteps])
ylim([0, 2])


%3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

figure(3)
subplot(1,2,1)
plot(x,result_conservative.massFlow(1,:), 'linewidth',2, 'color', 'b')
hold on
plot(x,result_conservative.massFlow(50,:), 'linewidth',2, 'color', 'g')
hold on
plot(x,result_conservative.massFlow(100,:), 'linewidth',2, 'color', 'r')
hold on
plot(x,result_conservative.massFlow(150,:), 'linewidth',2, 'color', 'c')
hold on
plot(x,result_conservative.massFlow(200,:), 'linewidth',2, 'color', 'm')
hold on
plot(x,result_conservative.massFlow(totalTimeSteps,:), 'linewidth',2, 'color', 'y')
xlabel("X-Domain")
ylabel("Non Dimensional Mass Flow Rate")
title("Conservative Form")
legend("1dt", "50dt", "100dt", "150dt", "200dt", "steady-state", 'location', 'best')
ylim([0,2])

subplot(1,2,2)
plot(x,result_nonconservative.massFlow(1,:), 'linewidth',2, 'color', 'b')
hold on
plot(x,result_nonconservative.massFlow(50,:), 'linewidth',2, 'color', 'g')
hold on
plot(x,result_nonconservative.massFlow(100,:), 'linewidth',2, 'color', 'r')
hold on
plot(x,result_nonconservative.massFlow(150,:), 'linewidth',2, 'color', 'c')
hold on
plot(x,result_nonconservative.massFlow(200,:), 'linewidth',2, 'color', 'm')
hold on
plot(x,result_nonconservative.massFlow(totalTimeSteps,:), 'linewidth',2, 'color', 'y')
title("Non-Conservative Form")
xlabel("X-Domain")
ylabel("Non Dimensional Mass Flow Rate")
legend("1dt", "50dt", "100dt", "150dt", "200dt", "steady-state", 'location', 'best')
ylim([0,2])
sgtitle("Variation of Non Dimensional Mass flow rate distribution at different Time-Steps")


% 4. Comparison of Normalized mass flow rate distributions of both forms
figure(4)
plot(x,result_conservative.massFlow(totalTimeSteps,:), 'linewidth',2, 'color', 'r')
hold on
plot(x,result_nonconservative.massFlow(totalTimeSteps,:), 'linewidth',2, 'color', 'c')

xlabel("X-Domain")
ylabel("Non Dimensional Mass Flow Rate")
legend('Conservative Form', 'Non Conservative Form')
title('Comparison of Normalized mass flow rate distributions')

%% Grid Independence
% For grid independence Lets run the simulation for n = 61 points for
% conservation form
n = 61;
x = linspace(0,L,n);
A = 1 + 2.2*(x - 1.5).^2;
Rho = ones(1,n);
T = ones(1,n);
for m = 1:n
    if x(m) >= 0 && x(m) < 0.5
        Rho(m) = 1;
        T(m) = 1;
    end
    
    if x(m) >= 0.5 && x(m) < 1.5
        Rho(m) = 1 - 0.366 * (x(m) - 0.5);
        T(m) = 1 - 0.167 * (x(m) - 0.5);
    end
    
    if x(m) >= 1.5 && x(m) <= 3
        Rho(m) = 0.634 - 0.3879 * (x(m) - 1.5);
        T(m) = 0.833 - 0.3507 * (x(m) - 1.5);
    end
        
end
v = 0.59 ./ (Rho .* A);
conservative_initialCond_grid_i = initialCondition(Rho,T,v);
conservative_grid_i = quasiOnedim(n,L,gamma,conservative_initialCond_grid_i,C,totalTimeSteps);
[conservative_grid_i,result_conservative_grid_i] = conservativeCD(conservative_grid_i);




% throat index when grid size is changed
throat_index_grid_i = find(A==min(A));

diff_rho = max(max(abs(result_conservative.rho(totalTimeSteps,throat_index)...
    - result_conservative_grid_i.rho(totalTimeSteps,throat_index_grid_i))));
diff_M = max(max(abs(result_conservative.M(totalTimeSteps,throat_index)...
    - result_conservative_grid_i.M(totalTimeSteps,throat_index_grid_i))));
diff_P = max(max(abs(result_conservative.P(totalTimeSteps,throat_index)...
    - result_conservative_grid_i.P(totalTimeSteps,throat_index_grid_i))));
diff_T = max(max(abs(result_conservative.T(totalTimeSteps,throat_index)...
    - result_conservative_grid_i.T(totalTimeSteps,throat_index_grid_i))));

percent_change_rho = (diff_rho/result_conservative.rho(totalTimeSteps,throat_index))*100
percent_change_M = (diff_rho/result_conservative.M(totalTimeSteps,throat_index))*100
percent_change_P = (diff_rho/result_conservative.P(totalTimeSteps,throat_index))*100
percent_change_T = (diff_rho/result_conservative.T(totalTimeSteps,throat_index))*100

 

 

 

 Findings from above numerical analysis are as follows:-

1. Steady-state distribution of primitive variables inside the nozzle

The above plot represents the steady state variation of different parameters along the domain. Visually the conservation and non conservation form results appear similar. As described below from the table where, the parameter valeus at the throat is compared, It can be seen that the non conservation form is more near to analytical soultion when the soultion has reached steady state.

 

2. Time-wise variation of the primitive variables

In the above plot where, transient variation in variables is presented. It can be seen that, in conservative form there is less variation in the inital time steps, when comared with non conservative form. The initial condition is more improved in the conservative form. After around 400 time steps the solution is seen to be stable visualy.

3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

Transient mass flow rate suffers from similar trend. Less variation is seen in the conservative form as compared with non conservative form. 

 

4. Comparison of Normalized mass flow rate distributions of both forms

Above plot represents the magnified view of steady form of both conservative as well as non conservative form of governing equations. It can be seen that there is less variation in conservative form as compared to the non conservative form.

 

Comparision of Steady state results, between conservation and Non Conservation form with Changing Grid Size

	No. of Grid Point	non dim `rho`	non dim `T`	non dim `P`	M
Exact Solution		0.634	0.833	0.528	1.0
Non Conservation form	31	0.6387	0.8365	0.5342	0.9994
Conservation Form	31	0.6507	0.8401	0.5466	0.9812
Conservation	61	0.6386	0.8353	0.5334	1.0001

 

Also Change in percentage change in values of parameters between no. of grid point 31 and 61 is too low.

for Non Dimensional Density its 1.8516%

for Mach No. its 1.2278%

for Non Dimensional pressure its 2.2039%

for Non Dimensional Temprature its 1.4340%

 

With this much less variation when grid size is almost doubled, we can say that we have achieved grid independence.