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.