Simulating a quasi-one-directional subsonic-supersonic nozzle flow using MacCormack s technique and drawing comparison between numerical solutions obtained from conservation and non-conservation forms of governing equations
Objectives
1. To write code to solve the Quasi 1D subsonic- supersonic nozzle flow equations in conservation and non- conservation forms of the governing equations.
2. To perform the grid dependency test.
3. To determine the minimum number of iterations required for convergence and to find out which form of governing equations results in faster convergence.
4. To compare the normalized mass flow rate between the conservation and non- conservation forms and arrive at suitable conclusions.
Problem Definition
This problem has been sourced from the book titled \'COMPUTER FLUID DYNAMICS- The Basics with Applications\' written by John D Anderson Jr in Chapter 7 and Section-7.3 & 7.5. The numerical solution is computed and compared with the solutions mentioned in the book.
We consider a steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir with a large cross-sectional area ( A tends to infinity), and hence the velocity is very small ( V tends to 0). The flow expands isentropically to supersonic speed at the nozzle exit. It is subsonic in the convergent section, sonic at the throat and supersonic in the divergent section. The flow properties are assumed to vary along only the x-direction although the area is two-dimensional. Such flow is called quasi-one-dimensional flow.
Governing equations in Non-Conservation form
Governing equations can be generally expressed in two forms, the non-conservation form and the conservation form. For the quasi 1D subsonic-supersonic isentropic flow nozzle, the non-conservation form can be expressed as;
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)`
Generally, for nozzle flows, the flow field variables are expressed as non-dimensional variables. This is done by dividing all the dimensional flow variables by their respective reservoir values.
Therefore, the non-dimensional flow variables are;
`x' = x/L ;" " rho' = rho/rho^o;" " T' = T/T^o;" " V' = V/a_o;" " a_o = sqrt(gamma R T_o);" " t'= t/(L/a_o) and A'= A/(A^**)`
where
`rho^o,T^o,a_o` are the density, temperature and speed of sound at the reservoir
`A^**` is the throat area, and
'L' is the nozzle length.
Let us now express the above non-conservation form using the non-dimensional variables.
Continuity equation (non dimensional)
`(del rho\')/(del t\') = -rho\' (del V\')/(del x\') - rho\' V\' (del (ln A\'))/(del x\') - V\'(del rho\')/(del x\')`
Momentum equation (non dimensional)
`(del V\')/(del t\') = -V\'(del V\')/(del x\') - 1/gamma((del T\')/(del x\') + (T\')/(rho\') (del rho\')/(del x\'))`
Energy equation (non dimensional)
`(del T\')/(del t\') = -V\'(del T\')/(del x\') - (gamma-1)T\'[(del V\')/(del x\') + V\' (del(ln A\'))/(del x\')]`
The non-dimensional variables x\',rho\',T\',V\' and t\' will be represented without the prime symbol on further usage for convenience and dimensional variables will not be used hereon further in the equations for non-conservation form.
Discretization
The discretization of the governing equation in non-conservation form is done by using the MacCormack technique. The MacCormack method is an explicit finite difference method that is 2nd order accurate in both time and space. Usually, to employ a time-marching solution that is 2nd order accurate, it is necessary to calculate the second order derivatives. The MacCormack method overcomes that by using a predictor-corrector set-up that can produce 2nd order accurate results without having to compute the 2nd order derivatives, thus making the process a whole lot simpler.
The predictor step computes the first order time derivatives using the forward difference for spatial derivatives by substituting previous state values of the flow variables. It then updates the flow field variables. These updated values are then used to compute corrector step time derivatives using the backward difference approach. Finally, the corrected solutions for the flow variables are obtained by using the average value of time derivatives computed in the predictor and the corrector steps. These steps are repeated until a steady state is reached or until predefined time steps.
Let us discretize the non-conservation form equations using MacCormack\'s approach.
Predictor step
* Computing time derivatives using forward difference scheme for space derivatives by substituting known values of flow variables at time t
`((del rho)/(del t))_i^t = -rho_i^t[(V_(i+1)^t - V_i^t)/(Delta x)] - rho_i^t V_i^t[(ln(A_(i+1)) - ln(A_i))/(Delta x)] - V_i^t[(rho_(i+1)^t - rho_i^t)/(Delta x)]`
`((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 the predicted values of flow field variables at time (t+dt) using the time derivatives at time t.
`(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`
where dt is the time step size satisfying the CFL condition (with C = 0.5)
`Delta t = C((Delta x)/(sqrt(T)+V))`
The time step size is calculated for all nodes and the minimum of all time step size values is takes as the value for \'dt\'.
Corrector step
* Computing time derivative values at time (t+dt) with the predicted values of flow variables.
`(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 the average of time derivative values obtained from foward difference and backward difference
`((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)]`
* Compute the corrected values of flow field variables using the average value of time derivatives.
`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)`
Nozzle shape and Initial conditions
The nozzle is a convergent-divergent nozzle of non-dimensional length 3. The nozzle area is minimum at the throat and increases as we progress on either side. So we define the area as a function of x,
`A = 1 + 2.22(x - 1.5)^2 | 0<=x<=3`
The initial conditions must be chosen as close to the steady state as possible through knowledge and experience. Otherwise, it can cause the solution to go unstable. We know that density and temperature increase and the velocity decreases as the flow expands inside the nozzle. Therefore, we choose the initial conditions that quantitatively describe the progressive increase and decrease.
At t=0
`rho = 1 - 0.3146 x`
`T = 1 - 0.2314 x`
`V = (0.1 + 1.09 x)T^(1/2)`
Boundary conditions
Inlet
At the subsonic inlet, we keep rho and T fixed and allow V to float i.e., it must be calculated as a function of timewise solution by linear extrapolation of the inner nodes.
V(1) = 2*V(2) - V(3)
Outlet
At the supersonic outlet, all three flow field variables are allowed to float.
rho(n) = 2*rho(n-1) - rho(n-1)
V(n) = 2*V(n-1) - V(n-2)
T(n) = 2*T(n-1) - T(n-2)
With this information, a matlab code can be written to simulate how the flow field variables vary across the nozzle and with time.
Governing equations in the Conservation from
For the quasi 1D subsonic-supersonic isentropic flow nozzle, the conservation form of the governing equations can be expressed as;
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)`
Expressing the above conservation form equation using non-dimensional variables, we get;
Continuity equation (non-dimensional)
`(del(rho\' A\'))/(del t\') + (del(rho\' A\' V\'))/(del x\') = 0`
Momentum equation (non-dimensional)
`(del(rho\' A\' V\'))/(del t\') + (del[rho\' A\' V\'^2 + (1/gamma)p\'A\'])/(del x\') = 1/gamma p\' (del A\')/(del x\')`
Energy equation (non-dimensional)
`(del[rho\'((e\')/(gamma-1) + gamma/2 V\'^2)A\'])/(del t\') + (del[rho\'((e\')/(gamma-1) + gamma/2 V\'^2)V\'A\' + rho\'A\'V\'])/(del x\')`
Let U be the solutions vector, F be the flux vector and J be the source term, such that
`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\' (del A\')/(del x\')`
Substituting these terms in the conesrvation form equations, we get;
`(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\')`
While writing the code, we\'ll solve for the solution vectors and then consecutively calculate the primitive variables using the following relations;
`rho\' = (U\')/(A\')`
`V\' = U_2/U_1`
`T\' = e\' = (gamma-1)(U_3/U_1-gamma/2 V\'^2)`
`p\' = rho\' T\'`
Now that we\'re solving for the solution vectors, the RHS consisting of flux terms and source terms must also be expressed in terms of U1, U2 and U3 to obtain the pure form of the flux terms. Otherwise, the solution may blow up after just few time steps.
`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\')`
Discretization
Discretization of governing equations in the conservation form will be done using the MacCormack method.
Predictor step
* The flux terms and the source terms are calculated at nodes i and i+1 using previous state values of U1, U2 and U3 and forward difference is applied to compute the value of the time derivatives.
`((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\')`
*These time derivative values are used to calculate the predicted values of U1, U2 and U3.
`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
* New set of flux and source term values are calculated using the predicted values of U1, U2 and U3 and backward difference is applied to compute the time derivative values for the corrector step.
`(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\')`
* The average value of the predictor and corrector step derivatives is calculated.
`((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)]`
*The corrected values of U1, U2 and U3 are computed using the average of derivate 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\'`
Nozzle shape and Initial conditions
The nozzle shape and area is the same as above. The area as a function of x is given as;
`A = 1 + 2.22(x - 1.5)^2`
Since the conservation form governing equations are more sensitive to initial conditions, it is necessary to make them more realistic.
`0<=x\'<=0.5 {(rho\',=,1),(T\',=,1):}`
`0.5<=x\'<=1.5 {(rho\',=,1-0.366(x\'-0.5)),(T\',=,1-0.167(x\'-0.5)):}`
`1.5<=x\'<=3 {(rho\',=,0.634-0.3879(x\'-1.5)),(T\',=,0.833-0.3507(x\'-1.5)):}`
The variation of V\' for initial conditions can be synthesized by assuming a constant mass flow (U2) for initial conditions only. U2 = 0.59 is chose for initial conditions.
`V\' = U_2/(rho\' A\')`
But the dependent variables being solved by the governing equations are U1, U2 and U3. We need initial conditions for them. The same initial conditions can be applied to U1, U2 and U3 using the relations;
`U_1 = rho\' A\'`
`U_2 = rho\' A\' V\' = 0.59`
`U_3 = rho\'((e\')/(gamma-1)+gamma/2 V\'^2)A\'`
Boundary Conditions
Inlet
At the sub-sonic inlet, rho\' and T\' are fixed and V\' is allowed to float.
In terms of the dependent variables, U1 will be fixed, U2 will float and as U3 is a function of V\', it will also float.
`U_1(1) = rho\'(1) *A\'(1)`
`U_2(1) = 2 U_2(2) - U_2(3)`
`V\'(1) = (U_2(1))/(rho\'(1)*A\'(1))`
`U_3(1) = U_1(1)((T\'(1))/(gamma-1) + gamma/2 V\'(1)^2)`
Outlet
At super-sonic outlet, all the dependent variables are allowed to float.
U1(1) = 2*U1(n-1) - U1(n-2)
U2(1) = 2*U2(n-1) - U2(n-2)
U3(1) = 2*U3(n-1) - U3(n-2)
With all the above information, let us write a matlab code to solve quasi 1d nozzle flow problem using both non-conservation and conservation forms.
Matlab Script
Main program
% Main script- Simulation of quasi 1D subsonic-supersonic isentropic flow
% nozzle in non-conservation and conservation forms
clear all
close all
clc
L = 3;
n = 31;
x = linspace(0,3,n);
nt = 1400;
C = 0.5;
tol = 1e-4;
m_analytical = 0.579*ones(n);
% solving governing equations in the non-conservation form
[rho_nc,T_nc,p_nc,M_nc,m_nc,n_nc] = non_conservation_form(L,n,nt,C,tol);
% solving giverning equations in the consevtion form
[rho_c,T_c,p_c,M_c,m_c,n_c] = conservation_form(L,n,nt,C,tol);
% Comparing the normalized mass flow rate between conservation &
% non-conservation forms
figure(7)
plot(x,m_nc,x,m_c,x,m_analytical,\'g\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_oV_oA*\');
title(sprintf([\'Comparison of steady state mass flow variations obtained with\'...
\'\\nconservation and non-conservation forms\']));
legend(\'Non-conservation form\',\'Conservation form\',\'Analytical solution\');
% comparing the convergence rate
figure(8)
c = categorical({\'Non-conservation form\',\'Conservation form\'});
b = bar(c, [n_nc n_c]);
b.FaceColor = \'flat\';
b.CData(2,:) = [1 1 0];
ylabel(\'No. of iterations\');
title(sprintf([\'Numer of iterations taken to converge to a steady-state\'...
\'\\n for a tolerance of %g\'], tol));
% Computing execution times accurately using timeit()
myfun_nc = @()non_conservation_form(L,n,nt,C,tol);
time_nc = timeit(myfun_nc,6);
myfun_c = @()conservation_form(L,n,nt,C,tol);
time_c = timeit(myfun_c,6);
% displaying the execution times on the output window
fprintf(\'Time taken to execute non-conservation form = %0.3gs\', time_nc);
fprintf(\'\\nTime taken to execute conservation form = %0.3gs\', time_c);
Function to solve non-conservation form
% function to compute solution to quasi 1D subsonic-supersonic isentropic
% flow nozzle in non-conservation form
function [rho,T,p,M,m,n_iter]= non_conservation_form(L,n,nt,C,tol)
% initialization
x = linspace(0,L,n);
dx = L/(n-1);
A = 1 + 2.2*(x-1.5).^2;
throat = (n+1)/2; % assuming equal no of grid points on either side of throat
nt_x = linspace(1,nt,nt);
gamma = 1.4;
error = 1;
n_iter = 1;
m_old = 0;
% Initial conditions
rho = 1 - 0.3146*x;
T = 1 - 0.2314*x;
v = (0.1 + 1.09.*x).*T.^(1/2);
% analytical steady-state mass flow rate solution
m_an = 0.579*ones(n);
% copy of initial values
rho_initial = rho;
T_initial = T;
v_initial = v;
m_initial = rho.*v.*A;
% loop to calculate minimum value of dt among all nodes
for k = 1:n
dt_x(k) = C*dx/(T(k)^(1/2)+v(k));
end
dt = min(dt_x);
% time loop
for j = 1:nt
% updating old values at the start of every time loop
rho_old = rho;
v_old = v;
T_old = T;
% predictor step
for i = 2:n-1
drho_dx = (rho(i+1)-rho(i))/dx;
dv_dx = (v(i+1)-v(i))/dx;
dT_dx = (T(i+1)-T(i))/dx;
dlnA_dx = (log(A(i+1))-log(A(i)))/dx;
% computing the derivatives using forward difference
% continuity equation
drho_dt_p(i) = -rho(i)*dv_dx-rho(i)*v(i)*dlnA_dx-v(i)*drho_dx;
% momentum equation
dv_dt_p(i) = -v(i)*dv_dx-(1/gamma)*(dT_dx+(T(i)/rho(i))*drho_dx);
% energy equation
dT_dt_p(i) = -v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlnA_dx);
% predictor step update
rho(i) = rho_old(i)+drho_dt_p(i)*dt;
v(i) = v_old(i)+dv_dt_p(i)*dt;
T(i) = T_old(i)+dT_dt_p(i)*dt;
end
% corrector step
for i = 2:n-1
drho_dx = (rho(i)-rho(i-1))/dx;
dv_dx = (v(i)-v(i-1))/dx;
dT_dx = (T(i)-T(i-1))/dx;
dlnA_dx = (log(A(i))-log(A(i-1)))/dx;
% computing the derivatives using rearward difference
% continuity equation
drho_dt_c(i) = -rho(i)*dv_dx-rho(i)*v(i)*dlnA_dx-v(i)*drho_dx;
% momentum equation
dv_dt_c(i) = -v(i)*dv_dx-(1/gamma)*(dT_dx+(T(i)/rho(i))*drho_dx);
% energy equation
dT_dt_c(i) = -v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlnA_dx);
end
% average value of time derivatives
drho_dt = 0.5*(drho_dt_p+drho_dt_c);
dv_dt = 0.5*(dv_dt_p+dv_dt_c);
dT_dt = 0.5*(dT_dt_p+dT_dt_c);
% corrector step update
for i = 2:n-1
rho(i) = rho_old(i)+drho_dt(i)*dt;
v(i) = v_old(i)+dv_dt(i)*dt;
T(i) = T_old(i)+dT_dt(i)*dt;
end
% Inlet boundary condition
v(1) = 2*v(2)-v(3);
% Outlet boundary conditions
v(n) = 2*v(n-1)-v(n-2);
rho(n) = 2*rho(n-1)-rho(n-2);
T(n) = 2*T(n-1)-T(n-2);
% Calculating pressure, mass flow rate and Mach number
p = rho.*T;
M = v./T.^0.5;
m = rho.*v.*A;
m_t(j,:) = rho.*v.*A; % capturing timewise variation of massflow at all nodes
% Capturing flow field variables at the throat for each time step
rho_throat(j) = rho(throat);
T_throat(j) = T(throat);
p_throat(j) = p(throat);
M_throat(j) = M(throat);
% checking convergence
error = max(abs(m-m_old));
if error>tol
n_iter = n_iter + 1;
m_old = m;
end
end
% plotting the steady state solutions
figure(1);
sgtitle(sprintf([\'Steady-state distribution of flow-field variables across nozzle\'...
\'\\n-Non-conservation form, No. of time steps = %g\'],nt));
subplot(4,1,1);
plot(x,M,\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'M\');
title(sprintf(\'Mach number\'));
subplot(4,1,2);
plot(x,rho,\'r\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'\\rho/\\rho_o\');
title(sprintf(\'Non-dimensional Density\'));
subplot(4,1,3);
plot(x,T,\'g\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'T/T_o\');
title(sprintf(\'Non-dimensional Temperature\'));
subplot(4,1,4);
plot(x,p,\'m\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'p/p_o\');
title(sprintf(\'Non-dimensional Pressure\'));
% Timewise variarion of flow field variables at nozzle throat
figure(2)
plot(nt_x,rho_throat,\'LineWidth\',2);
hold on
plot(nt_x,T_throat,\'LineWidth\',2);
hold on
plot(nt_x,p_throat,\'LineWidth\',2);
hold on
plot(nt_x,M_throat,\'LineWidth\',2);
title(sprintf([\'Timewise variation of non-dimensional flow-field variables\'...
\'\\nat nozzle throat- Non-conservation form\']));
xlabel(\'Number of time steps\');
legend(\'\\rho/\\rho_o\',\'T/T_o\',\'p/p_o\',\'M\');
% Mass flow variation across nozzle at different time steps
figure(3)
plot(x,m_initial,\'--\',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m_an,\':\',\'LineWidth\',2);
xlabel(\'Non-dimensional distance through nozzle(x)\');
ylabel(\'Non-dimensional mass flow (\\rhoVA)/(\\rho_oa_oA*)\');
legend(\'0\\Deltat\',\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\',\'Analytical solution\');
title(sprintf([\'Distribution of non-dimensional mass flow rate across the nozzle\'...
\'\\n-Non-conservation form\']));
end
Function to solve conservation form
% function to compute solution to quasi 1D subsonic-supersonic isentropic
% flow nozzle in conservation form
function [rho,T,p,M,m,n_iter]= conservation_form(L,n,nt,C,tol)
% Initialization
x = linspace(0,L,n);
dx = L/(n-1);
A = 1 + 2.2*(x-1.5).^2;
throat = (n+1)/2; % assuming equal no of grid points on either side of throat
nt_x = linspace(1,nt,nt);
gamma = 1.4;
error = 1;
n_iter = 1;
m_old = 0;
% Initial conditions
u2_i = 0.59; % initial mass flow rate, at t=0;
for i= 1:n
if ((x(i)>=0) && (x(i)<=0.5))
rho(i) = 1;
T(i) = 1;
elseif ((x(i)<1.5) && (x(i)>0.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
v(i) = u2_i/(rho(i)*A(i));
u1(i) = rho(i)*A(i);
u2(i) = rho(i)*A(i)*v(i);
u3(i) = rho(i)*A(i)*(T(i)/(gamma-1)+(gamma/2)*v(i)^2);
end
% analytical steady-state mass flow rate solution
m_an = 0.579*ones(n);
% copy of initial values
u1_initial = u1;
u2_initial = u2;
u3_initial = u3;
% loop to calculate minimum value of dt among all nodes
for k = 1:n
dt_x(k) = C*dx/(T(k)^(1/2)+v(k));
end
dt = min(dt_x);
%dt = 0.0267;
% time loop
for j = 1:nt
% updating old values at the start of every time loop
u1_old = u1;
u2_old = u2;
u3_old = u3;
% predictor step
for i = 2:n-1
% computing flux and source terms
f1_i = u2(i);
f1_ip1 = u2(i+1);
f2_i = u2(i)^2/u1(i)+((gamma-1)/gamma)*(u3(i)-(gamma/2)*u2(i)^2/u1(i));
f2_ip1 = u2(i+1)^2/u1(i+1)+((gamma-1)/gamma)*(u3(i+1)-(gamma/2)*u2(i+1)^2/u1(i+1));
j2 = (1/gamma)*rho(i)*T(i)*((A(i+1) - A(i))/dx);
f3_i = gamma*u2(i)*u3(i)/u1(i)-(gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_ip1 = gamma*u2(i+1)*u3(i+1)/u1(i+1)-(gamma*(gamma-1)/2)*u2(i+1)^3/u1(i+1)^2;
% computing the derivatives using forward difference
du1_dt_p(i) = -(f1_ip1-f1_i)/dx; % continuity eq
du2_dt_p(i) = j2-(f2_ip1-f2_i)/dx; % momentum eq
du3_dt_p(i) = -(f3_ip1-f3_i)/dx; % energy eq
% updating the predicted values
u1(i) = u1(i)+du1_dt_p(i)*dt;
u2(i) = u2(i)+du2_dt_p(i)*dt;
u3(i) = u3(i)+du3_dt_p(i)*dt;
end
% Updating primitive variables after predictor step
rho = u1./A;
v = u2./u1;
T = (gamma-1)*((u3./u1) - (gamma/2)*v.^2);
% corrector step
for i = 2:n-1
% computing new set of flux and source terms using predicted
% values for the dependent variables
f1_i = u2(i);
f1_im1 = u2(i-1);
f2_i = u2(i)^2/u1(i)+((gamma-1)/gamma)*(u3(i)-(gamma/2)*u2(i)^2/u1(i));
f2_im1 = u2(i-1)^2/u1(i-1)+((gamma-1)/gamma)*(u3(i-1)-(gamma/2)*u2(i-1)^2/u1(i-1));
j2 = (1/gamma)*rho(i)*T(i)*((A(i) - A(i-1))/dx);
f3_i = gamma*u2(i)*u3(i)/u1(i)-(gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_im1 = gamma*u2(i-1)*u3(i-1)/u1(i-1)-(gamma*(gamma-1)/2)*u2(i-1)^3/u1(i-1)^2;
% computing the derivatives using rearward difference
du1_dt_c(i) = -(f1_i-f1_im1)/dx;
du2_dt_c(i) = j2-(f2_i-f2_im1)/dx;
du3_dt_c(i) = -(f3_i-f3_im1)/dx;
end
% computing average value of time derivatives
du1_dt = 0.5*(du1_dt_p+du1_dt_c);
du2_dt = 0.5*(du2_dt_p+du2_dt_c);
du3_dt = 0.5*(du3_dt_p+du3_dt_c);
for i = 2:n-1
% computing corrected values
u1(i) = u1_old(i)+du1_dt(i)*dt;
u2(i) = u2_old(i)+du2_dt(i)*dt;
u3(i) = u3_old(i)+du3_dt(i)*dt;
end
% applying boundary conditions
% inlet
u1(1) = rho(1)*A(1);
u2(1) = 2*u2(2)-u2(3);
v(1) = u2(1)/(rho(1)*A(1));
u3(1) = u1(1)*(T(1)/(gamma-1)+(gamma/2)*v(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);
% computing the values of the primitive variables
for i= 1:n
rho(i) = u1(i)/A(i);
v(i) = u2(i)/u1(i);
T(i) = (gamma-1)*(u3(i)/u1(i)-gamma/2*v(i)^2);
end
% Calculating pressure, mass flow rate and Mach number
p = rho.*T;
M = v./T.^0.5;
m = rho.*v.*A;
m_t(j,:) = rho.*v.*A; % capturing timewise variation of massflow at all nodes
% Capturing flow field variables at the throat for each time step
rho_throat(j) = rho(throat);
T_throat(j) = T(throat);
p_throat(j) = p(throat);
M_throat(j) = M(throat);
% checking convergence
error = max(abs(m-m_old));
if error>tol
n_iter = n_iter + 1;
m_old = m;
end
end
% plotting steady state solutions of flow-field variables
figure(4);
sgtitle(sprintf([\'Steady-state distribution of flow-field variables across nozzle\'...
\'\\n-Conservation form, No. of time steps = %g\'],nt));
subplot(4,1,1);
plot(x,M,\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'M\');
title(sprintf(\'Mach number\'));
subplot(4,1,2);
plot(x,rho,\'r\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'\\rho/\\rho_o\');
title(sprintf(\'Non-dimensional Density\'));
subplot(4,1,3);
plot(x,T,\'g\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'T/T_o\');
title(sprintf(\'Non-dimensional Temperature\'));
subplot(4,1,4);
plot(x,p,\'m\',\'LineWidth\',2);
xlabel(\'x/L\');
ylabel(\'p/p_o\');
title(sprintf(\'Non-dimensional Pressure\'));
% Timewise variarion of flow field variables at nozzle throat
figure(5)
plot(nt_x,rho_throat,\'LineWidth\',2);
hold on
plot(nt_x,T_throat,\'LineWidth\',2);
hold on
plot(nt_x,p_throat,\'LineWidth\',2);
hold on
plot(nt_x,M_throat,\'LineWidth\',2);
title(sprintf([\'Timewise variation of non-dimensional flow-field variables\'...
\'\\nat nozzle throat- Conservation form\']));
xlabel(\'Number of time steps\');
legend(\'\\rho/\\rho_o\',\'T/T_o\',\'p/p_o\',\'M\');
% Mass flow variation across nozzle at different time steps
figure(6)
plot(x,u2_initial,\'--\',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m_an,\':\',\'LineWidth\',2);
xlabel(\'Non-dimensional distance through nozzle(x)\');
ylabel(\'Non-dimensional mass flow (\\rhoVA)/(\\rho_oa_oA*)\');
legend(\'0\\Deltat\',\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\',\'Analytical solution\');
title(sprintf([\'Distribution of non-dimensional mass flow rate across the nozzle\'...
\'\\n-Conservation form\']));
endResults
1. Qualitative aspects of quasi 1d nozzle flow
The steady-state solution profiles of all flow-field variables are the same for both conservation and non-conservation forms and concur with the reference profile in John D Anderson Jr\'s book.
3. Convergence rate
An idea of when the system reaches the steady state in both forms of the governing equations can be had by looking at the timewise variation of the flow-field variables at the throat of the nozzle.
It can be seen that, in both conservation and non-conservation forms, at around 400 iterations the profiles begin to straighten, meaning that they are approaching a steady-state solution. The oscillations seem to completely die out after 600 time-steps in the non-conservation form but in conservation form, small oscillations are observed. This means that the conservation form takes longer to converge.
The non-dimensional mass flow rate simulation across the nozzle can also be used to determine when the system reaches a steady-state.
The above graphs can be used to determine when the steady-state is reached and when the calculations can be stopped.
Alternatively, we can employ a convergence criterion inside the loop and stop the calculations as soon as it is met. To reach a tolerance of 1e-4, the number of iterations taken by both forms of governing equations are as follows;
Therefore, it can be inferred that the non-conservation is faster. The actual time taken by each form is determined by using the timeit function.
Comparing the normalized mass flow rate between conservation and non-conservation forms
Comparing the non-dimensional mass flow rate variation across the nozzle, the conservation form provides a much better mass flow rate distribution, that is closer to the analytical solution. Conservation form conserves the mass flow rate better as it is one of the dependent variables being directly solved by the governing equations, unlike the non-conservation form where the primitive variables are solved in the governing equations and mass flow is calculated.
2. Grid dependency test
The grid dependency test is done to determine the extent to which the solutions change, on changing the grid size. An ideal solution does not change with the number of grid points inside a domain, or gid size. To check our solutions for grid independence, let us compare the steady-state values of flow-field variables at nozzle throat for n=31 and n=61.
From the above throat value comparison of flow-field variables, it can be seen that for non-conservation form, the solutions change by the order of 1e-3. For conservation form. the change is of the order 1e-2. Since the change is very small, it can be concluded that the solutions are grid independent.
Conclusions
The quasi 1D subsonic-supersonic isentropic nozzle flow problem sourced from John D Anderson\'s \'Computer Fluid Dynamics- The basics with applications\' Chapter-7, Section- 7.3 & 7.5, was successfully simulated and the results concur with that of the author\'s mentioned in the same book. During this exercise the following conclusions were made;
1. Non-conservation form is easier to program.
2. Non-conservation form has faster convergence to steady-state.
3. Conservation form preserves mass better and has mass-flow rate distribution closer to the analytical value.
4. Both conservation and non-conservation forms are grid independent.