// 27.10.17
// **********************************************************
// Sahib Tofig oglu Huseynov                                *
// National Aviation Academy (NAA) of Azerbaijan           *
// Department of aircraft and power systems                 *
// ************************************************ *********
//***********************************************************************
// SUMMARY OF THE COMPUTATIONAL PROGRAM, COMPOSED IN THE SCILAB LANGUAGE 
//***********************************************************************
// The algorithm of this program is composed according to the method 
// of characteristics. As a result of numerical simulation, we can 
// determine the optimal shape of the rocket nozzle.
// It is solving following optimization problem:
// 1) designing of an ideal nozzle shape with cosearning of pressure and
// some important parameters in the combustion chamber of solid propellant.
// 2) obtasining the minimum length of the nozzle and optimal dimensions
// of the cross-section of the nozzle (i.e., the optimal geometric 
// parameters of the nozzle profile and the diameters of the throat and
// outlet cross section of the nozzle).
// Main purpose of the algorithm is to achieve maximum value of the 
// objective function (thrust of the reactive engine).That is, in result
// of computation it will obtained the max value of thrust force of the 
// reactive engine with min values of length of nozzle for the given 
// values of the ambient pressure and  for main parameters of combustion
// chamber the jet engine.
// Also we can determine the Mach number along the nozzle
// ************************************************ ***********************
// *** Close all previous drowed windows with graphics *******************
close(); close(); close(); close()
//************************************************************************
// *** Main physical and geometrical input data of the problem ************
T_c = 3358;          // - temperature in the combustion chamber, [K];;
P_c = 9e6;           // - pressure in the combustion chamber, [Pa];
P_amb = 101e3;       // - ambient pressure, [Pa];
T_amb = 300;         // - ambient temperature, [K];
gamma = 1.6;         // - ratio of Specific Heats Cp/Cv (Gamma); 
W = 19.531;          // - molecular weight of gas (kg/kmol);
width = 0.3535;      // - nozzle width (meters)
h_th = 0.1945;       // - throat height (meters)
//************************************************************************
// ASSIGNED DATA FOR THE PARAMETERS OF THE METHOD OF CHARACTERISTICS
//************************************************************************
num = 100;             // - number of characteristic lines,
theta_i = .03;        //  - initial step in theta,           
max_iter = 10000;     //  - maximum value of iterations.                   
//************************************************************************
R = 8314/W;
dh = h_th/100;  // - it is the step size.
//*****    First Section   *********************************************
// Obtaing of the condition when P = u.
h(1) = h_th;
A_star = h_th*width;
M =1;
dM1 = .1;
for i=1: max_iter
    h(i) = h(1) + (i-1)*dh;
    Ae(i) = h(i)*width;
    A_Asq = (Ae(i)/A_star)^2;
    A_ratio(i)=sqrt(A_Asq);   
    //  Now is useing the algorithm of the Newton - Rhapson method.
    res = 1;
    if i > 1
        M = Ma(i-1);
 end   
   while res > .001
        M2 = M + dM1;
        funa1 = -A_Asq + (1/M^2)*((2/(gamma+1))*(1+(gamma-1)*M^2/2))^((gamma+1)/(gamma-1));
        funa2 = -A_Asq + (1/M2^2)*((2/(gamma+1))*(1+(gamma-1)*M2^2/2))^((gamma+1)/(gamma-1));
        dv_dm = (funa2-funa1)/dM1;        
        M = M - funa1/dv_dm;
        res = abs(funa1);    
    end
    Ma(i) = M;  
// Computation of preassure.
    P(i) = P_c*(1+(gamma-1)*Ma(i)^2/2)^(-gamma/(gamma-1));   
// Computation thrust force along nozzle. 
    Te(i) = T_c/(1+(gamma-1)*Ma(i)^2/2);
    Tt(i) = T_c/(1+(gamma-1)/2);
    Ve(i) = Ma(i)*sqrt(Te(i)*gamma*R);
    Vt(i) = sqrt(Tt(i)*gamma*R);
    rhot(i) = P(i)/(R*Te(i));
    mdot(i) = rhot(i)*Ve(i)*Ae(i);
    TT(i) = mdot(i)*Ve(i) + (P(i) - P_amb)*Ae(i);   
if P(i) < P_amb
        // Break of computation
        // Obtaining the thrust force in all points
        P_exit = P(i)*(1+(gamma*2/(gamma+1))*(Ma(i)^2-1));       
         if P_exit <= P_amb
             P(i) = P_exit;
             break
         else
         end  
      else
    end 
end
figure(1)
//ppp=plot(Ae,TT)
plot(Ae,TT,'Color','b','LineWidth',3)
title('Curve of thrust force of reactive engine')
xlabel('S_1, [m^2)')
ylabel('P_T, [N]')
//*********************************************************************** 
//   Maximum value of the engine thrust force 
//***********************************************************************
[a,b]=max(TT); 
b = b;
A_max = Ae(b);
Max_thrust = TT(b);
plot(A_max,Max_thrust,'ro','MarkerFaceColor','r','MarkerEdgeColor','b','MarkerSize',7)
legend('P_T, [N],','P_T_max , [N]')
//******************************** SECOND SECTION ***********************
//                            METHOD OF CHARACTERISTIC                     
//***********************************************************************
M_e = Ma(b);       // Mach number at the ideal nozzle exit  
// Computation of theta_max  
theta_max = (180/%pi)*(sqrt((gamma+1)/(gamma-1))*atan((sqrt((gamma-1)*(M_e^2-1)/(gamma+1))))-atan(sqrt(M_e^2-1)))/2;
//  Computetion of D_theta for each characteristic line.
del_theta = (theta_max - theta_i)/(num-1);
// Computation of theta_max
for i=1:num
 //   Prosedure of initialization of Mach numeber.    
    for j=1:num
        if i==1
 //Computation of theta for each char. line (both for first lines). 
            theta(i,j) = theta_i + del_theta*(j-1);
            nu(i,j) = theta(i,j);
            K_m(i,j) = theta(i,j) + nu(i,j);
            K_p(i,j) = theta(i,j) - nu(i,j);              
        elseif i > 1        
            K_p(i,j) = -K_m(1,i);
            
       // Define theta
            if j >= i
                theta(i,j) = del_theta*(j-i);
            else                 
                theta(i,j) = theta(j,i);                
            end
            nu(i,j) = theta(i,j) - K_p(i,j);
            K_m(i,j) = theta(i,j) + nu(i,j);
        end   
    // Computation of Prandtl-Meyer function by using of Newton - Rhapson
    // method.
    dM = .1; // When dM = 0.1 then end of iteration
    if j == 1
        M_ex(i,j) = 1.00;
    else
        M_ex(i,j) = M_ex(i,j-1);
    end
    M = M_ex(i,j);       
    res = 1;
    while res > .01
        M2 = M + dM;
        funv1 = (-nu(i,j)*(%pi/180)+(sqrt((gamma+1)/(gamma-1))*atan((sqrt((gamma-1)*(M^2-1)/(gamma+1))))-atan(sqrt(M^2-1))));
        funv2 = (-nu(i,j)*(%pi/180)+(sqrt((gamma+1)/(gamma-1))*atan((sqrt((gamma-1)*(M2^2-1)/(gamma+1))))-atan(sqrt(M2^2-1))));
        dv_dm = (funv2-funv1)/dM;        
        M = M - funv1/dv_dm;
        res = abs(funv1);       
    end
    M_ex(i,j) = M; 
 // Computation of the mu angle.
    mu(i,j) = (180/%pi)*asin(1/M_ex(i,j)); 
    end   
    // Procedure for adding the last point to the characteristic line
    theta(i,num+1) = theta(i,num);
    nu(i,num+1) = nu(i,num);
    K_m(i,num+1) = K_m(i,num);
    K_p(i,num+1) = K_p(i,num);
end
char = zeros(num,num+1,2);
for i=1:num    
    for j=1:num+1       
//   Computation of intersection pointsş    
//   It is a point with number 1 for all characteristic lines           
        if j == 1 
            char(i,j,1) = 0;
            char(i,j,2) = h_th/2;
        end        
// First line hits the symmetry line.
        if i == 1 & j==2            
            char(i,j,1) = (-h_th/2)/tan((%pi/180)*(theta(1,j-1)-mu(1,j-1)));
            char(i,j,2) = 0;
        end
//  All other lines hit the symmetry line.
        if j == i+1 & j>2            
              char(i,j,1) = -char(i-1,j,2)/tan((%pi/180)*(.5*theta(i,j-2)-.5*(mu(i,j-2)+mu(i,j-1)))) + char(i-1,j,1);
              char(i,j,2) = 0;
              test(i,j) = (theta(i,j-2)-.5*(mu(i,j-2)+mu(i,j-1)));
              testpty(i,j) = char(i-1,j,2);
              testptx(i,j) = char(i-1,j,1);                
        end            
//   Computation of all data points for charteristic line 1.
        if i ==1 & j>2 & j ~= i+1
            C_p = tan((%pi/180)*(.5*(theta(i,j-2)+theta(i,j-1))+.5*(mu(i,j-2)+mu(i,j-1))));
            C_m = tan((%pi/180)*(.5*(theta(j-1,1)+theta(i,j-1))-.5*(mu(j-1,1)+mu(i,j-1))));
            A = [1,-C_m;1,-C_p];
            B = [char(1,1,2) - char(1,1,1)*C_m;
                char(1,j-1,2) - char(1,j-1,1)*C_p];
                iterm(1,:)=inv(A)*B;
                char(i,j,1) = iterm(1,2);
                char(i,j,2) = iterm(1,1);
        end      
//   Computation of other points for all charteristic line.
        if i > 1 & j~=i+1 & j>2        
            C_p = tan((%pi/180)*(.5*(theta(i,j-2)+theta(i,j-1))+.5*(mu(i,j-2)+mu(i,j-1))));
            C_m = tan((%pi/180)*(.5*(theta(i-1,j-1)+theta(i,j-1))-.5*(mu(i-1,j-1)+mu(i,j-1))));
            A = [1,-C_m;1,-C_p];
            B = [char(i-1,j,2) - char(i-1,j,1)*C_m; char(i,j-1,2) - char(i,j-1,1)*C_p];   
            iterm(1,:) = inv(A)*B;
            char(i,j,1) = iterm(1,2);
            char(i,j,2) = iterm(1,1);  
        end
    end
end  
//  Fill in similar points to the points, where characteristic lines 
//  share the points.
for i = 2:num
    for j=2:num
        char(j,i,1) = char(i-1,j+1,1);
        char(j,i,2) = char(i-1,j+1,2);
    end
end        
// CREATION OF THE NOZZLE SHAPE AND OBTAINIG OF CHARACTERISTIC LINES
//                        ALONG NOZZLE SHAPE 
//  Assaining the initial point (start point). 
noz(1,1) = 0;
noz(1,2) = h_th/2;
//   Computation the all points of the nozzle.
for i = 2 : num
    // Computing of different slopes and points to intersect.
    m1 = tan((%pi/180)*(theta(i-1,num)+mu(i-1,num)));    
    if i ==2
        m2 = (%pi/180)*theta_max;
    else
        m2 = ((%pi/180)*(theta(i-1,num+1)));
    end
    m3 = ((%pi/180)*(theta(i-1,num)));
    m4 = tan((m2+m3)/2);  
    A = [1,-m4; 1,-m1];
    B = [noz(i-1,2) - noz(i-1,1)*m4; char(i-1,num+1,2) - char(i-1,num+1,1)*m1];                
    iterm(1,:) = inv(A)*B;
    noz(i,1) = iterm(1,2);
    noz(i,2) = iterm(1,1);    
// Extending of characteristic lines to the wall.
    char(i-1,num+2,1)= noz(i,1);
    char(i-1,num+2,2)= noz(i,2);
end
// Computetion of last line.
m1 = tan((%pi/180)*(theta(num,num)+ mu(num,num)));
m2 = ((%pi/180)*(theta(num-1,num)));
m3 = ((%pi/180)*(theta(num,num+1)));
m4 = tan((m2+m3)/2);
A = [1,-m4; 1,-m1];
B = [noz(num,2) - noz(num,1)*m4; char(num,num+1,2) - char(num,num+1,1)*m1];               
iterm(1,:) = inv(A)*B;
noz(num+1,1) = iterm(1,2);
noz(num+1,2) = iterm(1,1);  
//   Characteristic lines are drowing to the nozzle shape line.
char(num,num+2,1)= noz(num+1,1);
char(num,num+2,2)= noz(num+1,2);
// **** PLOTTING OF NOZZLE SHAPE BOTH WITH CHARACTERISTIC LINES ********
figure(2);clf;
//Ploting of characteristic lines itaration.
for i = 1 : num  
   plot(char(i,:,1),char(i,:,2), '.-g')      
   plot(char(i,:,1),-char(i,:,2), '.-g')
end
   //Drow of upper and lower lines of nozzle shape.
plot(noz(:,1),noz(:,2),'k','LineWidth',3)
plot(noz(:,1),-noz(:,2),'k','LineWidth',3)
title('Nozzle shape with characteristic lines  (by defining of maximum value of engine thrust with minimum nozzle length)')
xlabel('Nozzle length, [m]')
ylabel('Cross section diameter of the nozzle, [m]')
//  Computation of errors in A/A* and Mexit
error_Area = 100*(width*2*noz(num,2) - A_max)/(A_max)
error_Mach = 100*(M_e - M_ex(num,num))/M_e
// Plot Mach Number and pressure through nozzle using the quasi-1D
// area relations.  (Isentropic expansion through nozzle)
Mnoz(1) = 1.0;  //  Computing of choked flow.
M = Mnoz(1);
for i=1: size(noz,1)
    Ae(i) = 2*noz(i,2)*width;
    A_Asq = (Ae(i)/A_star)^2;
    A_ratio(i)=sqrt(A_Asq);    
    // Again application of Newton - Rhapson method. 
    res = 1;
    if i > 1
        M = Mnoz(i-1);    
     while res > .001
        M2 = M + dM1;
        funa1 = -A_Asq + (1/M^2)*((2/(gamma+1))*(1+(gamma-1)*M^2/2))^((gamma+1)/(gamma-1));
        funa2 = -A_Asq + (1/M2^2)*((2/(gamma+1))*(1+(gamma-1)*M2^2/2))^((gamma+1)/(gamma-1));
        dv_dm = (funa2-funa1)/dM1;      
        M = M - funa1/dv_dm;
        res = abs(funa1);       
    end
    Mnoz(i) = M;
    end
    // Computation of preassure.
    Pnoz(i) = P_c*(1+(gamma-1)*Mnoz(i)^2/2)^(-gamma/(gamma-1));
end
figure(3);
plot(noz(:,1),Mnoz,'-r','LineWidth',3,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',10)  
plot(noz(size(noz,1),1),M_e,'ro','MarkerFaceColor','y','MarkerEdgeColor','b','MarkerSize',10)  
plot(noz(size(noz,1),1),1) 
xlabel('Nozzle length, [m]')
ylabel('Mach number')
h1=legend(' -Mach number',' -M_outlet (predicted)',[2])
title('Mach number depend on nozzle length and predicted mach number on the nozzle outlet.')
figure(4);
plot(noz(:,1),Pnoz/P_amb,'.r','Color','black','LineWidth',3,'MarkerSize',7,'MarkerEdgeColor','r') 
plot(noz(size(noz,1),1),1,'ro','MarkerFaceColor','y','MarkerEdgeColor','b','MarkerSize',10)
xlabel('Nozzle length, [m]')
ylabel('P/P_c')
legend('Pnoz/P_amb','Pnoz/P_amb in outlet of the nozzle')
title('Pnoz/P_amb depend on nozzle length and its value on the nozzle outlet.')
