Simulation of a 1D Super-sonic nozzle flow
In this project we will observe quasi 1D flow through a convergent - divergent nozzle. We will be refering to the case set up for this problem from Chapter 7 of Computational Fluid Dynamics - John D Anderson.
Physical description of problem
- The flow at the inlet comes from a reservior where the pressure and temperature is `p_0` and `T_0`
- The area at the inlet is very large (`A -> Inf`) so accordingly velocity becomes very low (`v -> 0`). However, `v` can not be equal to zero as there will be no mass flow rate in that case.Thus, `p_0` and `T_0` are stagnation values.
- At the inlet the flow is subsonic i.e `M <1`, as the flow expands isentropically ,the flow becomes sonic at throat i.e `M =1` and supersonic at the outlet `M>1`.
- We assume that at given cross section the flow properties are uniform at that section , even though the area is function of `xy` we make the assumption that the flow properties vary only along `x` direction. Such a flow is termed as quasi-one dimensional flow.
Governing Equations
The governing equations are driven by considering a small control volume. If the control volume is fixed in space while the fluid around it is in motion, then the resulting equations would take a conservative form, whereas when a control volume of fixed mass moving in fluid is considered , then non-conservative form of equations are obtained. Mathematically, it can be proved that both the equations are same.
Equations in Non-Conservative Form:
Continuity Equation : `frac{delrho}{delt} = -rhofrac{delV}{delx} - rhoVfrac{del(lnA)}{delx} - Vfrac{delrho}{delx}` `---(1)`
Momentum Equation: `frac{delV}{dt} = -Vfrac{delV}{delx} - frac{1}{gamma}(frac{delT}{delx} + frac{T}{rho}frac{delp}{delx})` `---(2)`
Energy Equation: `frac{delT}{delt} = -Vfrac{delT}{delx} -(gamma - 1)T[frac{delV}{delx} + Vfrac{del(lnA)}{delx}]` `---(3)`
Equations in Conservative Form:
Continuity Equation : `frac{del(rhoA)}{delt} + frac{del(rhoAV)}{delx} = 0` `---(4)`
Momentum Equation: `frac{del(rhoAV)}{delt} + frac{del[rhoAV^2 + (1/gamma)pA]}{delx} = frac{1}{gamma}p(frac{delA}{delx}) ` `---(5)`
Energy Equation: `frac{del[rho(frac{e}{gamma-1} + frac{gamma}{2}V^2)A]}{delt} + frac{del[rho(frac{e}{gamma-1} + frac{gamma}{2}V^2)AV + rhoAV]}{delx} = 0` `---(6)`
where all variables `rho, A, V, T, p, e` are in dimensionless form.
Density, `rho = frac{rho}{rho_0}` Area, `A = frac{A}{A^**}` Length of nozzle, `x = frac{x}{L}`
Velocity, `V = frac{V}{a_0}` Time, `t = frac{t}{(L/a_0)}` Temperature, `T = frac{T}{T_0}` Speed of sound, `a_0 = sqrt(gammaRT_0)`
The conservative form of equation contains divergence terms . Hence they can be moulded in in a generic form described as follows :
`frac{delU}{delt} + frac{delF}{delx} + frac{delG}{dely} + frac{delH}{delz} = J`
where `F,G,H` are called as flux terms, `J` represents source term and `U` represents solution vector.
Since the equations `(4), (5), (6)` of conservative form can not be discretized directly , let the equations be described in form of solution vector `U`, flux terms `F` and source terms `J` .
`U_1 = rhoA`
`U_2 = rhoAV`
`U_3 = rho(frac{e}{gamma-1} + frac{gamma}{2}V^2)A`
`F_1 = rhoAV`
`F_2 = rhoAV^2 + frac{1}{gamma}rhoA`
`F_3 = rho(frac{e}{gamma-1} + frac{gamma}{2}V^2)AV + rhoAV`
`J_2 = frac{1}{gamma}rhofrac{delA}{delx}`
Substituting above values in eqs `(4), (5), (6)` we get
`frac{delU_1}{delt} = - frac{delF_1}{delx}` `---(7)`
`frac{delU_2}{delt} = - frac{delF_2}{delx} + J_2` `---(8)`
`frac{delU_3}{delt} = - frac{delF_3}{delx}` `---(9)`
Further simplifying above terms we get,
`F_1 = U_2`
`F_2 = frac{U_2^2}{U_1} + frac{gamma - 1}{gamma}(U_3 - frac{gamma}{2}frac{U_2^2}{U_1})`
`F_3 = gammafrac{U_2U_3}{U_1} - frac{gamma(gamma - 1)}{2}frac{U_2^3}{U_1^2}`
`J_2 = frac{1}{gamma}rhoTfrac{delA}{delx} = frac{gamma - 1}{gamma}(U_3 - frac{gamma}{2}frac{U_2^2}{U_1})frac{del(lnA)}{delx}`
For conservative form equations `(7), (8), (9)` will be discretized using above terms, while for non-conservative form equations `(1), (2), (3)` will be discretized.
Solver : MacCormack\'s Technique
This technique is a variant of Lax-Wendroff approach. This method is also an explicit finite difference technique , accurate in second order in both space and time.
In the predictor step the space terms will be discretized forward in space , while in corrector step the terms will be discretized backward in space. At the end of each time step iteration average of these values will be used to update the solution.
Calculation of time Step
As we know that at larger time step values solution of explicit difference scheme tend to blow up , hence a constraint has to be put on the value of time step that will be used.
`Deltat = Cfrac{delx}{a + V}`
For linear hyperbolic system , the Courant Number should follow `Cle1`. However, our current system under consideration is a system of non -linear equations. Hence , the relation `Cle1` will be considered as general guidance system.
At a given grid point `i` and at given time step `t` value of time step can be calculated as ,
`(Deltat)_i^t = Cfrac{delx}{a_i^t + V_i^t}`
At an adjacent grid point the value of time step will be
`(Deltat)_(i+1)^t = Cfrac{delx}{a_(i+1)^t + V_(i+1)^t}`
Now to perform time marching solution , we have two options
1. To calculate `(Deltat)` at each grid point `i` at each step. However, the shortcoming of this method is that at `t +Deltat` it becomes artificial time wrap for flow field variables at given grid point. The value of flow field variables will correspond to non physical time. Hence , such a method will not replicate an actual physical time marching solution and hence can not be used for an accurate solution for unsteady flow.
2. To calculate `(Deltat)` at all grid points from `i=1,2,3...N` and calculate the minimum value among them. In this way the value calculated at `t +Deltat` will correspond to actual physical time. The shortcoming of this method will be that it will require more time steps to approach steady state.
For calculations of time step we will be using approach (2)
Initial Conditions
For Non - Conservative Form we will consider following profiles to set up initial conditions :
`rho = 1 - 0.3146x`
`T = 1 - 0.23146x`
`V = (0.1 + 1.09x)T^(1/2)`
For Conservative Form :
for `0 le x le 0.5` `{(rho = 1.0),(T = 1.0):}`
for `0.5 le x le 1.5` `{(rho = 1.0 - 0.366(x - 0.5)),(T = 1.0 - 0.167(x - 0.5)):}`
for `0 le x le 0.5` `{(rho = 0.634 - 0.3789(x - 1.5)),(T = 0.833 - 0.3507(x - 1.5)):}`
A constant mass flow rate through nozzle is assumed for calculating initial velocity.Hence, `U_2 = rhoAV = 0.59`
Hence `V = frac{U_2}{rhoA} = frac{0.59}{rhoA}`
Boundary Conditions
For Non-Conservative Form
The unsteady, inviscid flow is governed by hyperbolic equations and therefore for 1D flow there exists two characertistic lines through any point in `xt` plane. This is shown in figure below. Point 1 refers to the inlet point , where the flow is subsonic. Point N refers to outlet of nozzle where the flow is supersonic. The direction of stream line flow is from left to right i.e. from point `1` to point `N`.
At the subsonic point, since the local velocity is less velocity of sound hence, the left running characteristic at point 1 travels upstream against the low velocity sub sonic flow which is moving slowly from left to right. Hence, the resultant velocity of left running characteristic will be `a_1 - V_1`.This left running characteristic is moving out of the domain as the domain under consideration is from point `1` to point `N`.
However, the right running characteristic is moving to the right at i.e into the domain. This is due to the fact that the fluid element at point `1` is already moving into the domain. Also the right running characteristic is moving to the right at velocity of `a_1 + V_1`
The choice of boundary conditions will be decided using the above method of characteristics. The method of characteristics tells us that at a boundary where one characteristic flows into the domain then value of one dependent flow field variable must be specified at that boundary and if one characteristic flows out of the domain then value of one dependent flow field variable must be allowed to float at that boundary.
Also as noted above , stream line also enters the flow of domain, hence another flow field variable must be specified at the inlet boundary.
So at the inlet two values must be specified while one value of flow field variable must be allowed to float.
Applying similar logic to the super sonic point at the outlet, we observe that the left running characteristic is propagating towards the left. However, since the fluid element itself is super sonic the left running characteristic gets swept downstream relative to nozzle with a relative velocity of `V_N - a_N`. The right running characteristic also flows downstream with a relative velocity of \'V_N + a_N\'. As both the characteristics are moving out of the domain at the outlet, hence no flow field variables need to be specified at the outlet boundary. All dependent flow field variables will be allowed to float.
The variable chosen to float at the sub sonic inlet is Velocity `V`. Hence the value at point `1` will be calculated by linear extrapolation.
`V_1 = 2V_2 - V_3`. Other flow field variables will be fixed and hence will be independent with respect to time.
At super sonic outflow all flow field variables are allowed to float, hence by applying linear extrapolation we get following relations.
`V_N = 2V_(N-1) - V_(N-2)`
`rho_N = 2rho_(N-1) - rho_(N-2)`
`T_N = 2T_(N-1) - T_(N-2)`
For Conservative Form
Applying similar method of characteristic approach for solving conservative form of equations , we get following relations:
Inlet : Two dependent flow field variables to be specified at boundary while one is allowed to float.
`U_(1(i=1)) = rhoA`
`U_(2(i=1)) = 2U_(2(i=2)) - U_(2(i=3))`
`U_(3(i=1)) = U_1(frac{T}{gamma - 1} + frac{gamma}{2}V^2)`
Outlet : All flow field variables are allowed to float.
`(U_1)_N = 2(U_1)_(N-1) - (U_1)_(N-2)`
`(U_2)_N = 2(U_2)_(N-1) - (U_2)_(N-2)`
`(U_3)_N = 2(U_3)_(N-1) - (U_3)_(N-2)`
Code
Main Driver script:
close all;
clear all;
clc;
grid_points = [31 61];
solver_type = [\'1\' \'2\'];
% Solver 1 Conservative Form
% Solver 2 Non Conservative Form
numerical_values = [];
analytical_values_at_throat = [0.634 0.833 0.528 1.000];
throat_points = [16 31];
for g=1:length(grid_points)
throat_point = throat_points(g);
n = grid_points(g);
for solver=1:length(solver_type)
x=linspace(0,3,n);
dx = x(2) - x(1);
gamma_p = 1.4;
L = 1;
C = 0.5;
nt = 1400;
if solver_type(solver) == \'1\'
[elapsed_time_conservative_form, A_conservative,rho_conservative, v_conservative, T_conservative, M_conservative, p_conservative] = function_conservative_form(n, x, dx, gamma_p, L, C, throat_point, nt);
elseif solver_type(solver) == \'2\'
[elapsed_time_non_conservative_form, A_non_conservative, rho_non_conservative, v_non_conservative, T_non_conservative, M_non_conservative, p_non_conservative] = function_non_conservative_form(n, x, dx, gamma_p, L, C, throat_point, nt);
end
end
disp(sprintf(\"Grid Points %d\", n));
disp(sprintf(\"Elapsed Time for Conservative Form %d\", elapsed_time_conservative_form));
disp(sprintf(\"Elapsed Time for Non-Conservative Form %d\", elapsed_time_non_conservative_form));
vals = [rho_conservative(throat_point), T_conservative(throat_point), p_conservative(throat_point), M_conservative(throat_point)];
numerical_values(g, 1:length(vals)) = vals ;
f3 = figure(3, \'name\', \'Conservative vs Non-Conservative\');
plot(x, rho_conservative .* A_conservative .* v_conservative);
hold on;
plot(x, rho_non_conservative .* A_non_conservative .* v_non_conservative);
xlabel(\"Length of domain\");
ylabel(\"Mass flow rate\")
title(\"Mass flow rate after steady state\");
legend(\'ConservativeForm\', \'Non-ConservativeForm\');
saveas(f3, sprintf(\'Conservative_vs_NonConservativeForm_%d_GridPoints\',n), \'jpg\');
f4 = figure(4, \'name\', \"Qualitative Aspects of Quasi 1D nozzle flow\");
title(\'Qualitative Aspects : Conservative flow vs Non-Conservative flow\');
subplot(4, 1, 1)
plot(x, M_conservative)
hold on;
plot(x, M_non_conservative, \"*\")
xlabel(\'x\')
ylabel(\'Mach Number\')
legend(\'Conservative Form\', \'Non-Conservative Form\');
subplot(4, 1, 2)
plot(x, p_conservative)
hold on;
plot(x, p_non_conservative, \"*\")
xlabel(\'x\')
ylabel(\'Pressure ratio\')
legend(\'Conservative Form\', \'Non-Conservative Form\');
subplot(4, 1, 3)
plot(x, rho_conservative)
hold on;
plot(x, rho_non_conservative, \"*\")
xlabel(\'x\')
ylabel(\'Density ratio\')
legend(\'Conservative Form\', \'Non-Conservative Form\');
subplot(4, 1, 4)
plot(x, T_conservative)
hold on;
plot(x, T_non_conservative, \"*\")
xlabel(\'x\')
ylabel(\'Temperature ratio\')
legend(\'Conservative Form\', \'Non-Conservative Form\');
saveas(f4, sprintf(\'QualitativeAspects_Conservative_vs_NonConservativeForm_%d_GridPoints\',n), \'jpg\');
end
numerical_values(end+1, 1:length(analytical_values_at_throat)) = analytical_values_at_throat;
numerical_values = num2str(numerical_values);
rowNames = [\'31Points\';\'61Points\'; \'ExactSolution\'];
output = horzcat(rowNames, numerical_values);
disp(\"Comparison of throat values for varying grid points\");
fprintf(\"\\tRho\\tT\\tp\\tM\\n\");
disp(output);
Function for Conservative Form:
function [elapsed_time, A, rho, v, T, M, p] = function_conservative_form(n, x, dx, gamma_p, L, C, throat_point, nt)
close all;
%Density and Temperature Profile
for i=1:n
if (x(i) >= 0) && (x(i) <=0.5)
rho(i) = 1;
T(i) = 1;
elseif (x(i) >= 0.5) && (x(i) <=1.5)
rho(i) = 1 - 0.366*(x(i) - 0.5);
T(i) = 1 - 0.167*(x(i) - 0.5);
elseif (x(i) >= 1.5 && x(i)<=3.5)
rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
T(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
A = 1 + 2.2*(x - 1.5).^2;
p = rho .* T;
%Flux Variables
U1 = rho.* A;
U2 = 0.59*ones(1, n);
v = U2 ./ U1;
U3 = U1 .* (T./(gamma_p - 1) + 0.5 * gamma_p * (v .^ 2));
f1 = figure(1, \'name\', \'Mass flow rate variation\');
plot(x, rho .* A .*v)
hold on;
for i=1:n
dt_array(i) = min(C * (dx ./(sqrt(T(i)) + v(i)) ));
end
dt = min(abs(dt_array));
tic;
%Outer Loop
for k=1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
%Predictor Step
for j=2:n-1
J2 = (1/gamma_p) * rho(j) * T(j) * (A(j+1) - A(j))/dx;
F1_discrete = U2(j+1) - U2(j);
F2_term_first = ((U2(j+1)) ^ 2 / U1(j+1)) + ((gamma_p - 1)/gamma_p) * (U3(j+1) - 0.5*gamma_p*( ((U2(j+1))^2)/U1(j+1) ));
F2_term_second = ((U2(j)) ^ 2 / U1(j)) + ((gamma_p - 1)/gamma_p) * (U3(j) - 0.5*gamma_p*( ((U2(j))^2)/U1(j) ));
F2_discrete = F2_term_first - F2_term_second;
F3_term_first = gamma_p * ((U2(j+1) * U3(j+1))/U1(j+1)) - 0.5*gamma_p*(gamma_p-1)* ((U2(j+1))^3 / ((U1(j+1)))^ 2);
F3_term_second = gamma_p * ((U2(j) * U3(j))/U1(j)) - 0.5*gamma_p*(gamma_p-1)* ((U2(j))^3 / ((U1(j)))^ 2) ;
F3_discrete = F3_term_first - F3_term_second;
du1_dt_p(j) = -(F1_discrete)/dx;
du2_dt_p(j) = -(F2_discrete)/dx + J2;
du3_dt_p(j) = -(F3_discrete)/dx;
U1(j) = U1(j) + du1_dt_p(j) * dt;
U2(j) = U2(j) + du2_dt_p(j) * dt;
U3(j) = U3(j) + du3_dt_p(j) * dt;
rho(j) = U1(j)/A(j);
v(j) = U2(j)/U1(j);
% disp(v(j));
T(j) = (gamma_p - 1)*(U3(j)/U1(j) - 0.5 * gamma_p*v(j)**2);
end
% [(1:n)\' x\' U1\' U2\' U3\' rho\' T\']
%Corrector Step
for j=2:n-1
J2(j) = (1/gamma_p)*rho(j)*T(j)*(A(j) - A(j-1))/dx;
F1_discrete = U2(j) - U2(j-1);
F2_term_first = ((U2(j)) ^ 2 / U1(j)) + ((gamma_p - 1)/gamma_p) * (U3(j) - 0.5*gamma_p*( ((U2(j))^2)/U1(j) ));
F2_term_second = ((U2(j-1)) ^ 2 / U1(j-1)) + ((gamma_p - 1)/gamma_p) * (U3(j-1) - 0.5*gamma_p*( ((U2(j-1))^2)/U1(j-1) ));
F2_discrete = F2_term_first - F2_term_second;
F3_term_first = gamma_p * ((U2(j) * U3(j))/U1(j)) - 0.5*gamma_p*(gamma_p-1)* ((U2(j))^3 / ((U1(j))) ^ 2);
F3_term_second = gamma_p * ((U2(j-1) * U3(j-1))/U1(j-1)) - 0.5*gamma_p*(gamma_p-1)* ((U2(j-1))^3 / ((U1(j-1))) ^ 2);
F3_discrete = F3_term_first - F3_term_second;
du1_dt_c(j) = -(F1_discrete)/dx;
du2_dt_c(j) = -(F2_discrete)/dx + J2(j);
du3_dt_c(j) = -(F3_discrete)/dx;
end
% tt = [du1_dt_c\' du2_dt_c\' du3_dt_c\']
%Compute the average time
du1_dt_avg = 0.5 * (du1_dt_c + du1_dt_p);
du2_dt_avg = 0.5 * (du2_dt_c + du2_dt_p);
du3_dt_avg = 0.5 * (du3_dt_c + du3_dt_p);
%Final Solution Update
for j=2:n-1
U1(j) = U1_old(j) + du1_dt_avg(j) * dt;
U2(j) = U2_old(j) + du2_dt_avg(j) * dt;
U3(j) = U3_old(j) + du3_dt_avg(j) * dt;
end
%Apply Boundary Conditions
%Inlet
U2(1) = 2*U2(2) - U2(3);
%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);
%Finding Values of Primitive variables
for j=2:n-1
rho(j) = U1(j)/A(j);
v(j) = U2(j)/U1(j);
T(j) = (gamma_p - 1)*((U3(j)/U1(j)) - 0.5 * gamma_p*(v(j)^2));
end
p = rho.*T;
M = v./sqrt(T);
throat_rho(k) = rho(throat_point);
throat_v(k) = v(throat_point);
throat_t(k) = T(throat_point);
throat_p(k) = p(throat_point);
throat_mach(k) = M(throat_point);
f1 = figure(1,\'name\', \'Mass flow rate variation\') ;
if k == 50
plot(x, rho .* A .*v)
hold on;
elseif k == 100
plot(x, rho .* A .*v)
hold on;
elseif k == 150
plot(x, rho .* A .*v)
hold on;
elseif k == 200
plot(x, rho .* A .*v)
hold on;
elseif k == 700
plot(x, rho .* A .*v)
xlabel(\'Length of domain\');
ylabel(\'Non Dimensional Mass flow rate\')
title(\'Variation of mass flow rate at various time steps\')
end
end
legend(\'0\\Deltat\', \'50\\Deltat\', \'100\\Deltat\', \'150\\Deltat\', \'200\\Deltat\', \'700\\Deltat\')
elapsed_time = toc;
% [x\' A\' rho\' v\' T\' p\' M\']
f2 = figure(2, \'name\', \'Flow field variables variation\');
title(\"Flow field variables variation at throat\");
subplot(4, 1, 1)
plot(1:nt, throat_rho)
xlabel(\"Time steps\");
ylabel(\"Density\");
title(\'Density\');
subplot(4, 1, 2)
plot(1:nt, throat_t)
xlabel(\"Time steps\");
ylabel(\"Temperature\");
title(\'Temperature\');
subplot(4, 1, 3)
plot(1:nt, throat_p)
xlabel(\"Time steps\");
ylabel(\"Pressure\");
title(\'Pressure\');
subplot(4, 1, 4)
plot(1:nt, throat_mach)
xlabel(\"Time steps\");
ylabel(\"Mach Number\");
title(\'Mach Number\');
saveas(f1, sprintf(\'ConservativeForm_MassFlowRate_%d_GridPoints\', n), \'jpg\');
saveas(f2, sprintf(\'ConservativeForm_FlowVariables_%d_GridPoints\', n), \'jpg\');
endFunction for Non-Conservative Form
function [elapsed_time, A, rho, v, T, M, p] = function_non_conservative_form(n, x, dx, gamma_p, L, C, throat_point, nt)
close all;
%initial conditions
rho = 1 - 0.3146 * x;
T = 1 - 0.2314*x;
v = (0.1 + 1.09*x).* T.^0.5 ;
A = 1 + 2.2*(x - 1.5).^2;
dt = min(C * (dx ./ (sqrt(T) + v) ));
f1 = figure(1, \'name\', \'Mass Flow rate Variation\');
plot(x, rho .* A .*v)
hold on;
% Outer time loop
tic;
for k=1:nt
rho_old = rho;
v_old = v;
T_old = T;
% Predictor Method
for j=2:n-1
dvdx_p = (v(j+1) - v(j))/dx;
drhodx_p = (rho(j+1) - rho(j))/dx;
dtdx_p = (T(j+1) - T(j))/dx;
dloga_dx_p = (log(A(j+1)) - log(A(j)))/dx;
%continuity equation
drho_dt_p(j) = -rho(j) * dvdx_p - rho(j) * v(j) * dloga_dx_p - v(j) * drhodx_p;
%momentum equation
dvdt_p(j) = - v(j) * dvdx_p - (1 / gamma_p) * ( dtdx_p + (T(j)/ rho(j)) * drhodx_p );
%energy equation
dTdt_p(j) = -v(j) * dtdx_p - (gamma_p -1 ) * T(j) * (dvdx_p + v(j) * dloga_dx_p);
% Solution update
rho(j) = rho(j) + drho_dt_p(j) * dt;
v(j) = v(j) + dvdt_p(j) * dt;
T(j) = T(j) + dTdt_p(j) * dt;
end
% Corrector Method
for j=2:n-1
dvdx_c = (v(j) - v(j-1))/dx;
drhodx_c = (rho(j) - rho(j - 1))/dx;
dtdx_c = (T(j) - T(j -1))/dx;
dloga_dx_c = (log(A(j)) - log(A(j - 1)))/dx;
%continuity equation
drho_dt_c(j) = -rho(j) * dvdx_c - rho(j) * v(j) * dloga_dx_c - v(j) * drhodx_c;
%momentum equation
dvdt_c(j) = - v(j) * dvdx_c - (1 / gamma_p) * ( dtdx_c + (T(j)/ rho(j)) * drhodx_c );
%energy equation
dTdt_c(j) = -v(j) * dtdx_c - (gamma_p -1 ) * T(j) * (dvdx_c + v(j) * dloga_dx_c);
end
% Compute Average Values
drho_dt_avg = 0.5*(drho_dt_p + drho_dt_c);
dv_dt_avg = 0.5*(dvdt_p + dvdt_c);
dT_dt_avg = 0.5*(dTdt_p + dTdt_c);
% Final Solution Update
for j=2:n-1
rho(j) = rho_old(j) + drho_dt_avg(j) * dt;
v(j) = v_old(j) + dv_dt_avg(j) * dt;
T(j) = T_old(j) + dT_dt_avg(j) * dt;
end
p = rho .* T;
% Apply the boundary conditions
%Inlet
v(1) = 2*v(2) - v(3);
%Outlet
rho(n) = 2 *rho(n-1) - rho(n-2);
v(n) = 2 *v(n-1) - v(n-2);
T(n) = 2 *T(n-1) - T(n-2);
p(n) = rho(n) * T(n);
M = v ./ sqrt(T);
throat_rho(k) = rho(throat_point);
throat_v(k) = v(throat_point);
throat_t(k) = T(throat_point);
throat_p(k) = p(throat_point);
throat_mach(k) = M(throat_point);
throat_drho_dt_avg(k) = abs(drho_dt_avg(throat_point));
throat_dv_dt_avg(k) = abs(dv_dt_avg(throat_point));
f1 = figure(1, \'name\', \'Mass flow rate variation\');
if k == 50
plot(x, rho .* A .*v)
hold on;
elseif k == 100
plot(x, rho .* A .*v)
hold on;
elseif k == 150
plot(x, rho .* A .*v)
hold on;
elseif k == 200
plot(x, rho .* A .*v)
hold on;
elseif k == 700
plot(x, rho .* A .*v)
xlabel(\'Length of domain\');
ylabel(\'Non Dimensional Mass flow rate\')
title(\'Variation of mass flow rate at various time steps\')
end
end
legend(\'0\\Deltat\', \'50\\Deltat\', \'100\\Deltat\', \'150\\Deltat\', \'200\\Deltat\', \'700\\Deltat\')
elapsed_time = toc;
% [x\' A\' rho\' v\' T\' p\' M\']
f2 = figure(2, \'name\', \'Flow field variables variation\');
title(\"Flow field variables variation at throat\");
subplot(4, 1, 1)
plot(1:nt, throat_rho)
xlabel(\"Time steps\");
ylabel(\"Density\");
title(\'Density\');
subplot(4, 1, 2)
plot(1:nt, throat_t)
xlabel(\"Time steps\");
ylabel(\"Temperature\");
title(\'Temperature\');
subplot(4, 1, 3)
plot(1:nt, throat_p)
xlabel(\"Time steps\");
ylabel(\"Pressure\");
title(\'Pressure\');
subplot(4, 1, 4)
plot(1:nt, throat_mach)
xlabel(\"Time steps\");
ylabel(\"Mach Number\");
title(\'Mach Number\');
saveas(f1, sprintf(\'NonConservativeForm_MassFlowRate_%d_GridPoints\', n), \'jpg\');
saveas(f2, sprintf(\'NonConservativeForm_FlowVariables_%d_GridPoints\', n), \'jpg\');
endComparison of Mass flow rate between Conservative form and Non Conservative Form:
Conservative and Non-Conservative Mass Flow rate with Grid Points = 31
Conservative and Non-Conservative Mass Flow rate with Grid Points = 61
* From above figures it can be seen that in all the cases the mass flow rate fluctuates at initial time steps(`0Deltat, 50Deltat`). However, after 200 time steps it slowly starts to settle down and attains steady value after 700 time steps.
* The mass flow rate obtained using conservative form of equations is more milder in the sense that lesser flucutations of values are observed. This is due to the fact that we had assumed a constant mass flow in our initial conditions.
Comparison between Conservative and Non-Conservative mass flow rate after steady state :
For grid points = 31 For Grid Points = 61
However, if we compare the mass flow rate of both the approaches after 1400 steps we would observe that the value obtained by conservative form lies closer to the exact analytical solution.
The conservative form involves more computation as it involves conversion of primitive variables to flux terms and vice-versa.Hence, the time required to attain steady state is more in case of conservative form as compared to non-conservative form.This can be seen from the above screenshot of output. The non-conservative form is more faster than conservative form.
Attaining Convergence or Steady State:
To check when convergence has been attained we can either iterate our code upto a tolerance value of error or we can visually inspect for a steady state value of flow field variables.
Variation of Flow field variable for 31 Grid Points (Conservative and Non-Conservative Form)
Variation of Flow field variable for 61 Grid Points (Conservative and Non-Conservative Form)
From above graphs we can observe after about 700 time step iterations the profiles are in form of straight lines. Hence, the values of flow field variables doesn\'t change or there is no significant change . So it can be assert that steady state has been achieved.
Grid dependence test:
Let us assume we are using `N` grid points to find the solution of a particular case. The output obtained using `N` grid points is marked as output 1.
If we change the grid points to `2N` we are reducing the values of `Deltax` and possibly `Deltay` also if its a 2D case. If the solution is not grid indepedent then the output 2 obtained using `2N` grid points can vary from output 1 significantly or may not be stable and hence can not be accepted.
To solve this situation the number of grid points should be increased till a solution which is independent of grid points is obtained.
For a solution to be acceptable, the output obtained should not be a function of number of grid points i.e. the solution must be independent of grid points, and whatever, difference obtained is within acceptable tolerance limits. Considering 1-2% of deviation as acceptable we observe that in present case the solution is grid indepedent at 31 grid points.