Supersonic nozzle

Calculations were carried out using conservation and non-conservation form of Navier-Stokes Equation.

-Non-conservation form:

Number of cycles: most of the variations in values occurs in the initial number of time steps. But at most till 400-500 iterations, values dont change that much. And after about 1000 iterations, values almost remains constant.

Mass flow rate - for different time steps:

1. 100 steps

2. after 200 steps:

3. after 500 steps

4. after 700 steps

So basically after 700 time steps, value for mass flow rate seems to be stable.

-Grid dependence test: number of points were increased to 61,91 and it didnt effect the solution to much extent. Hence this code is not that much grid dependent.

Conservation form:

Code and Decode functions were used to calculate F using solution matrix, and primitive variables.

Number of cycles: In this case, convergence is reached faster in about 200-300 time steps as compared to that of non-conservation form, in which convergence was around 500 time steps.

Mass flow rate:

1. at 100 time steps

2. at 200 time steps:

3. at 500 time steps

4. at 700 time steps:

so basically convergence is reached faster in around 200-300 time steps, as one can judge from these graphs.

As through the course of reaching convergence , variations in solution for conservation form of equations is less as compared to that of non-conservation form. As judging from the mass flow rate graph, non-conservation form has large variations in mass-flow rate as compared to conservation form, in which mass flow rate remain closer to solution till convergence is reached.

Hence based upon results, conservation form appears to be faster and less varying through the course of convergence as compared to non-conservation form.

clear all
close all
clc

%input parameters-->using anderson data
n=31; % number of grid points
x=linspace(0,3,n); % grid space
dx=x(2)-x(1); % grid spacing
gamma=1.4;

%given functions of x
dens=1-0.3146*x; %density
T=1-0.2314*x;%temperature
v=(0.1+1.09*x).*T.^0.5; %velocity

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

nt=700; % number of time steps
dt=1e-2;%time step

method=2; %method =1 non-conservative form of solution  % anything else- conservation form

if method ==1
%%%%%Time Loop%%%%%
for k=1:nt
dens_old=dens;
T_old=T;
v_old=v;
      %%%Predictor step
      
      for j=2:n-1
          dv_dx=(v(j+1)-v(j))/dx; %forward diff on v
          ddens_dx=(dens(j+1)-dens(j))/dx; %forward diff on density
          dT_dx=(T(j+1)-T(j))/dx; %forward diff on temperature
          dlogA_dx=(log(A(j+1))-log((A(j))))/dx; %forward diff on area
          
          %continuity equation
          ddens_dt_p(j)=-dens(j)*dv_dx-dens(j)*v(j)*dlogA_dx-v(j)*ddens_dx; 
          
          %momentum equation
          dv_dt_p(j)=-v(j)*dv_dx-(1/gamma)* (dT_dx + (T(j)/dens(j))*ddens_dx);
          
          %energy equation
          dT_dt_p(j)=-v(j)*dT_dx-(gamma-1)*T(j)*(dv_dx+v(j)*dlogA_dx);
          
          %solution update
          v(j)=v(j)+dv_dt_p(j)*dt;
          dens(j)=dens(j)+ddens_dt_p(j)*dt;
          T(j)=T(j)+dT_dt_p(j)*dt;
      end
    
      %%%corrector step
      
      for j=2:n-1
          dv_dx=(v(j)-v(j-1))/dx; %forward diff on v
          ddens_dx=(dens(j)-dens(j-1))/dx; %forward diff on density
          dT_dx=(T(j)-T(j-1))/dx; %forward diff on temperature
          dlogA_dx=(log(A(j))-log((A(j-1))))/dx; %forward diff on area
          
          %continuity equation
          ddens_dt_c(j)=-dens(j)*dv_dx-dens(j)*v(j)*dlogA_dx-v(j)*ddens_dx; 
          
          %momentum equation
          dv_dt_c(j)=-v(j)*dv_dx-(1/gamma)* (dT_dx + (T(j)/dens(j))*ddens_dx);
          
          %energy equation
          dT_dt_c(j)=-v(j)*dT_dx-(gamma-1)*T(j)*(dv_dx+v(j)*dlogA_dx);
          
      end
    %%%average time derivative
    dv_dt=0.5*(dv_dt_p+ dv_dt_c);
    dens_dt=0.5*(ddens_dt_p+ ddens_dt_c);
    dT_dt=0.5*(dT_dt_p+ dT_dt_c);
    
    %final solution update
    for i=2:n-1
    dens(i)=dens_old(i)+dens_dt(i)*dt;
    v(i)=v_old(i)+dv_dt(i)*dt;
    T(i)=T_old(i)+dT_dt(i)*dt;
    end

    %%%applying boundary conditions
    v(1)=2*v(2)-v(3);
    T(1)=2*T(2)-T(3);
    dens(1)=2*dens(2)-dens(3);
    
    v(n)=2*v(n-1)-v(n-2);
    T(n)=2*T(n-1)-T(n-2);
    dens(n)=2*dens(n-1)-dens(n-2);
    plot(x,dens.*v.*A)
    title(sprintf(\'Time step=%d\',k))
    pause(0.01)
    axis([0 3 0 4])
end

else
    %%%%%%%%%conservative form of solution
    U=zeros(3,n);
    F=zeros(3,n);
    J=zeros(3,n);
    v=zeros(1,n);
    T=zeros(1,n);
    
    %%%%%%%%%%%%%%%%%%%%Boundary conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%inflow
    U(1,1)=A(1); % at inflow dens=1;T=1
    U(2,1)=2*U(2,2)-U(2,3);% then using this v is calculated
    v(1)=U(2,1)/U(1,1);
    U(3,1)=U(1,1)*((T(1,1)/(gamma-1))+gamma*0.5*v(1,1)^2);
    
    F(1,1)=U(2,1);
    F(2,1)=((U(2,1).^2)./U(1,1))+((gamma-1)/gamma)*(U(3,1)-(gamma*0.5*((U(2,1).^2)./U(1,1))));
    F(3,1)=gamma*(U(2,1).*U(3,1)./U(1,1))-gamma*0.5*(gamma-1)*(U(2,1).^3)./U(1,1).^2;
    
    %%outflow
    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);
     
    F(1,n)=U(2,n);
    F(2,n)=((U(2,n).^2)./U(1,n))+((gamma-1)/gamma)*(U(3,n)-(gamma*0.5*((U(2,n).^2)./U(1,n))));
    F(3,n)=gamma*(U(2,n).*U(3,n)./U(1,n))-gamma*0.5*(gamma-1)*(U(2,n).^3)./U(1,n).^2;
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %%%%%%%%%%%%%%%%%%%%%%initial conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 n1=(0.5/dx)+1; % for position at where x=0.5
 n2=(1.5/dx)+1; %for position at where x=1.5
 for i =1:n1
     dens(i)=1;
     T(i)=1;
 end
 
 for i =n1+1 : n2
     dens(i)=1-0.366*(x(i)-0.5);
     T(i)=1-0.167*(x(i)-0.5);
 end
 
 for i = n2+1 : n
     dens(i)=0.634-0.3879*(x(i)-1.5);
     T(i)=0.833-0.3507*(x(i)-1.5);
 end
 v=0.59./(dens.*A);
 
 U(1,:)=dens.*A;
 U(2,:)=dens.*A.*v;
 e=T;
 U(3,:)=dens.*((e/(gamma-1))+gamma*0.5*v.^2).*A;
[F]=Code(U,n);
 for i=1:n-1
     J(2,i)=(1/gamma)*dens(i)*T(i)*(A(i+1)-A(i))/dx;
 end
 J(2,n)=(1/gamma)*dens(n)*T(n)*(A(n)-A(n-1))/dx;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   
   dt=0.0267; % time step size
   nt=700;
   Up=zeros(3,n);
   dUp_dt=zeros(3,n);
   dUc_dt=zeros(3,n);
   vp=zeros(1,n);
   Tp=zeros(1,n);
   Jp=zeros(3,n);
  for m= 1:nt
%   %%%%predictor step%%%%

    for j=1:n-1
        dUp_dt(:,j)=J(:,j)-((F(:,j+1)-F(:,j))/dx);
    end
    dUp_dt(:,n)=J(:,n)-((F(:,n)-F(:,n-1))/dx);
    Up=U+dUp_dt*dt;
   %%%%%%%%%%%Setting boundary conditions for predictor step%%%%%%%%
   %%inflow
      Up(1,1)=A(1); % at inflow dens=1;T=1
      Up(2,1)=2*Up(2,2)-Up(2,3);% then using this v is calculated
      vp(1)=Up(2,1)/Up(1,1);
      Tp(1)=1;
      Up(3,1)=Up(1,1)*((Tp(1)/(gamma-1))+gamma*0.5*vp(1)^2);
%    
%    %%Outflow
     Up(1,n)=2*Up(1,n-1)-Up(1,n-2);
     Up(2,n)=2*Up(2,n-1)-Up(2,n-2);
     Up(3,n)=2*Up(3,n-1)-Up(3,n-2);
     vp(n)=Up(2,n)/Up(1,n);
% % %   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     [densp, Tp]=Decode(Up,x); %decoding primitive variables
     [Fp]=Code(Up,n);
      densp(1)=1;
      Tp(1)=1;

     for i=1:n-1
         Jp(2,i)=(1/gamma)*densp(i)*Tp(i)*(A(i+1)-A(i))/dx;
     end
    Jp(2,n)=(1/gamma)*densp(n)*Tp(n)*(A(n)-A(n-1))/dx;

% % % % % % %%%%%%%%%%%%%%%%%%%%%%%corrector step%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % 
      dUc_dt(:,1)=-(((Fp(:,2)-Fp(:,1))/dx)-Jp(:,1));
      for j=2:n
       dUc_dt(:,j)=-(((Fp(:,j)-Fp(:,j-1))/dx)-Jp(:,j));
      end
      dU_av_dt=0.5*(dUc_dt+dUp_dt);
      U_new=U+dU_av_dt*dt;
      %U=U_new;  
      %%inflow
      U_new(1,1)=A(1); % at inflow dens=1;T=1
     U_new(2,1)=2*U_new(2,2)-U_new(2,3);% then using this v is calculated
      v(1)=U(2,1)/U(1,1);
      T(1)=(gamma-1)*((U(3,1)/U(1,1))-gamma*0.5*((U(2,1)/U(1,1))^2));
      U_new(3,1)=U_new(1,1)*((T(1)/(gamma-1))+gamma*0.5*v(1)^2);
%   
%    %%Outflow
      U_new(1,n)=2*U_new(1,n-1)-U_new(1,n-2);
      U_new(2,n)=2*U_new(2,n-1)-U_new(2,n-2);
      U_new(3,n)=2*U_new(3,n-1)-U_new(3,n-2);
% % %   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       F=Code(U_new,n);
       [dens,T]=Decode(U_new,x); 
       v=U(2,:)./U_new(1,:);
       dens(1)=1;
        T(1)=1;
         for t=1:n-1
           J(2,t)=(1/gamma)*dens(t)*T(t)*(A(t+1)-A(t))/dx;
         end
         J(2,n)=(1/gamma)*dens(n)*T(n)*(A(n)-A(n-1))/dx;
         U=U_new;
         plot(x,U(2,:))
         title(sprintf(\'Time step=%d\',m))
         axis([0 3 0.5 1])
         pause(0.01)
  end
end


%%%%functions%%%%%
function [dens,T]=Decode(U,x)
A=1+2.2*(x-1.5).^2; %area
gamma=1.4;
dens=U(1,:)./A;
T=(gamma-1)*((U(3,:)./U(1,:))-gamma*0.5*((U(2,:)./U(1,:)).^2));
end
function [F]=Code(U,n)
gamma=1.4;
F=zeros(3,n);
 F(1,:)=U(2,:);
 F(2,:)=((U(2,:).^2)./U(1,:))+((gamma-1)/gamma)*(U(3,:)-(gamma*0.5*((U(2,:).^2)./U(1,:))));
 F(3,:)=gamma*(U(2,:).*U(3,:)./U(1,:))-gamma*0.5*(gamma-1)*(U(2,:).^3)./U(1,:).^2;
 end