Film_Cooling_SIm_V0_1

%----
%Film Cooling Sim
%V0.1
%----

close all
clear

%Variables --------
%Constants ----
R_univ = 8314;  %J/kg,K  Universal gas constant
SBC = 5.670373*10^-8;   %W/m2,K4  Stefan–Boltzmann constant
%----

%Simulation ----
n_Sim = 10^3;   %Number of steps in the sim
%----

%Engine Properties ----
%Fluid Properties (Taken from RPA)
P_c = 2*10^6;   %Pa Chamber pressure (stagnation)
T_c = 2948; %K  Chamber temperature (stagnation)
gamma_c = 1.156; %Heat capacity ratio
R_c = 321.1;    %J/kg,K  Gas constant
Lchar = 1.0;  %m  Characteristic length

%Physical Properties
D_t = 30/1000;  %m  Diameter of throat
D_c = D_t*sqrt(9);  %m  Diameter of chamber
D_e = D_t*sqrt(3.845);    %m  Diameter of nozzle exit
theta_Frus = 45;    %Degrees  Contraction angle of chamber frustrum
theta_Nzzl = 15;    %Degrees  Contraction anfle of nozzle (frustrum)
T_Wall_Max = 1000;  %K Maximum engine wall temperature
%----

%Coolant Properties ----
%Methanol

%Base Properties
MW_FC = 32.04;  %kg/kmol  Molar weight of coolant
rho_l_FC = 792; %kg/m3  Density of coolant as liquid
CP_l_FC = 2490; %J/kg,K  Specific heat of coolant as liquid
CP_g_FC = 1920; %J/kg,K  Specific heat of coolant as gas
k_l_FC = 0.202; %W/m,K  Thermal conductivity of coolant as liquid
k_g_FC = 0.0262; %W/m,K  Thermal conductivity of coolant as gas
mu_l_FC = 6.05*10^-4; %N,s/m2  Dynamic viscosity of coolant as liquid
mu_g_FC = 1.38*10^-5; %N,s/m2  Dynamic viscosity of coolant as gas
T_boil_FC = 337.85;  %K  Boiling temperature of coolant
DeltaHv_FC = 1.17*10^6;  %J/kg,K Specific heat of vaporization of coolant
epsilon_FC = 0.8;   %Emissivity of coolant

%Injection Properties
f_FC = 7/100;   %Fraction of total engine mass flux as coolant
T_inj_FC = 288;  %K  Temperature of coolant injection
u_inj_FC = 16;  %m/s  Coolant injection velocity
%----
%--------

%Preliminary Math ----
mDot_Engine = 0.25*pi*(D_t^2)*((P_c/sqrt(T_c))*sqrt(gamma_c/R_c)* ...
    (0.5*gamma_c+0.5)^(-1*(gamma_c+1)/(2*gamma_c-2))); %kg/s
mDot_FC = mDot_Engine*f_FC;   %kg/s  Mass flux of coolant

L_Frus = abs(D_c - D_t)/tand(theta_Frus)/2;    %m  Length of chamber frustrum
L_Nzzl = abs(D_e - D_t)/tand(theta_Nzzl)/2;    %m  Length of nozzle

V_c = 0.25*pi*(D_t^2)*Lchar;    %m3  Chamber volume
V_Frus = (pi*L_Frus/12)*(D_c*D_c+D_c*D_t+D_t*D_t);  %m3  Frustrum volume
V_Cyl = V_c - V_Frus;   %m3  Chamber cylinder volume
L_Cyl = V_Cyl/(0.25*pi*D_c^2);  %m  Chamber cylinder length
L_c = L_Cyl + L_Frus;   %m  Length of chamber
L_Engine = L_Cyl + L_Frus + L_Nzzl;  %m  Total length of engine

dD_Frus = (D_t - D_c)/L_Frus;
dD_Nzzl = (D_e - D_t)/L_Nzzl;

PR_l_FC = CP_l_FC*mu_l_FC/k_l_FC;   %Prandtl number of coolant as liquid
PR_g_FC = CP_g_FC*mu_g_FC/k_g_FC;   %Prandtl number of coolant as gas
%----

%Modeling the Engine ----
dL = L_Engine/(n_Sim - 1);   %m
L_Array = 0:dL:L_Engine;    %m
D_Array = zeros(1,n_Sim); %m
M_s_Array = zeros(1,n_Sim); %Mach
T_s_Array = zeros(1,n_Sim); %K
P_s_Array = zeros(1,n_Sim); %Pa
u_s_Array = zeros(1,n_Sim); %m/s

for i = 1:1:n_Sim
    L_i = L_Array(i);   %m
    D_i = 0;    %m
    SuperSonic_i = 0;   %'Is the flow supersonic?' bool

    if (L_i <= L_Cyl)
        D_i = D_c;  %m
    elseif (L_i > L_Cyl)&&(L_i <= L_c)
        D_i = dD_Frus*(L_i - L_Cyl) + D_c;  %m
    else
        D_i = dD_Nzzl*(L_i - L_c) + D_t;    %m
        SuperSonic_i = 1;
    end

    G_i = (D_i/D_t)^2;
    M_i = AreaMachFunc(G_i,gamma_c,SuperSonic_i,4);
    T_i = T_c*(1+0.5*(gamma_c-1)*M_i^2)^-1; %K
    P_i = P_c*(1+0.5*(gamma_c-1)*M_i^2)^(-1*gamma_c/(gamma_c-1));   %Pa
    u_i = M_i*sqrt(gamma_c*R_c*T_i);    %m/s

    if (L_i <= L_Cyl)
        u_i = (L_i/L_Cyl)*u_i;  %m/s
    end

    D_Array(i) = D_i;   %m
    M_s_Array(i) = M_i;   %Mach
    T_s_Array(i) = T_i; %K
    P_s_Array(i) = P_i; %Pa
    u_s_Array(i) = u_i; %m/s
end
%----

%Film Cooling Modeling --------
T_FC_Array = zeros(1,n_Sim);    %K
qa_Array = zeros(1,n_Sim);  %W/m2

%Liquid Phase ----
%It is assumed that the coolant will quickly vaporize.
%Therefore the cooling contribution during liquid phase will be...
%calculated in one step.
RE_D_l_FC = rho_l_FC*D_c*u_inj_FC/mu_l_FC;   %Reynold's number of coolant as liquid

%Gnielinski correlation for Nusselt number
DFF_l_FC = (0.79*log(RE_D_l_FC)-1.64)^-2;  %Darcy friction factor
NU_D_l_FC = ((DFF_l_FC/8)*(RE_D_l_FC-1000)*PR_l_FC)/...
    (1+12.7*sqrt(DFF_l_FC/8)*(-1 + PR_l_FC^(2/3))); %Nusselt number
h_l_FC = NU_D_l_FC*k_l_FC/D_c;  %W/m2,K  Heat transfer coeff.

%Heat Transfer
QDot_l_FC = mDot_FC*(CP_l_FC*(T_boil_FC - T_inj_FC) + DeltaHv_FC);  %W
T_l_FC_AVG = (T_boil_FC + T_inj_FC)/2;  %K

%Length of liquid film cooling
L_l_FC = QDot_l_FC/(h_l_FC*pi*D_c*(T_c - T_l_FC_AVG));  %m

%Heat Transfer per Area
A_l_FC = pi*D_c*L_l_FC; %m2  Area covered by liquid film cooling
qa_l_FC = QDot_l_FC/A_l_FC; %W/m2

%Where did the liquid film cooling end?
Delta_L_Array = abs(L_Array - L_l_FC);
[Delta_min,n_l_End] = min(Delta_L_Array);
T_FC_Array(n_l_End) = T_boil_FC;    %K

T_FC_Array(1,1:n_l_End) = ...
    T_inj_FC:(T_boil_FC-T_inj_FC)/(n_l_End-1):T_boil_FC;  %K
qa_Array(1,1:n_l_End) = qa_l_FC/(n_l_End); %W/m2
%----

%Gas Phase ----
for j = (n_l_End+1):1:n_Sim
    D_j = D_Array(j);   %m
    D_jp = D_Array(j-1);    %m
    T_s_jp = T_s_Array(j);   %K
    P_s_jp = P_s_Array(j);   %Pa
    u_s_jp = u_s_Array(j);    %m/s
    T_FC_jp = T_FC_Array(j-1); %K

    A_j = 0.25*pi*(D_j+D_jp)*sqrt(((abs(D_j-D_jp))^2)+4*dL^2);  %m2
    %A_j = pi*D_jp*dL;   %m2

    %Convective Heat Transfer with Exhaust
    u_FC_j = u_inj_FC; %m/s
    if (u_s_jp > u_FC_j)
        u_FC_j = u_s_jp;    %m/s
    end
    rho_jp = P_s_jp/((R_univ/MW_FC)*T_FC_jp);    %kg/m3
    RE_D_jp = rho_jp*u_FC_j*D_jp/mu_g_FC;
    DFF_jp = (0.79*log(RE_D_jp)-1.64)^-2;
    NU_jp = ((DFF_jp/8)*(RE_D_jp-1000)*PR_g_FC)/...
    (1+12.7*sqrt(DFF_jp/8)*(-1 + PR_g_FC^(2/3)));
    h_jp = k_g_FC*NU_jp/D_jp;   %W/m2,K

    QDot_Conv = h_jp*A_j*(T_s_jp - T_FC_jp);    %W

    %Convective Heat Transfer with Wall(add when wall temp can be modeled)
    QDot_Wall = 0;  %W
%     if (T_FC_jp > T_Wall_Max)
%         QDot_Wall = h_jp*A_j*(T_Wall_Max - T_FC_jp);   %W
%     end

    %Radiative Heat Transfer(add when wall temp can be modeled)
    %Also it's contribution is very small
    %QDot_Radi = -1*A_j*epsilon_FC*SBC*T_FC_jp^4;  %W
    QDot_Radi = 0;  %W

    %Summing
    QDot_j = QDot_Conv + QDot_Wall + QDot_Radi; %W
    DeltaT = QDot_j/(mDot_FC*CP_g_FC);  %K
    T_FC_j = DeltaT + T_FC_jp;  %K
    qa_j = QDot_j/A_j; %W/m2

    T_FC_Array(j) = T_FC_j; %K
    qa_Array(j) = qa_j; %W/m2
end
%----

T_FC_Max = max(T_FC_Array); %K  Maximum film coolant temperaure

[~,n_t] = min(abs(D_Array - D_t));
T_FC_t = T_FC_Array(n_t);   %K  Film coolant temperature at throat
%--------

%Plotting ----
% %Engine Shape
% figure
% plot(L_Array,D_Array/2,'k-')
% grid on
% axis equal
% title('Engine Shape')
% xlabel('Length (m)')
% ylabel('Radius (m)')
% xlim([0,L_Engine])
% ylim([0,1*D_c])

%Heat Flux per Area 3D
figure
[X,Y,Z] = cylinder(0.5*D_Array);
[m,n] = size(Z);
C = zeros(m,n);
for i = 1:1:n
    C(:,i) = qa_Array;
end
eng = surf(X,Y,L_Engine*(-1*Z+1),C);
eng.EdgeColor = 'none';
axis equal
cbar = colorbar;
cbar.Label.String = 'Watts-per-meter-squared';
title('Heat Flux per Area in Engine')
xlabel('Meters')
ylabel('Meters')
zlabel('Meters')

%Film Cooling Temperature 3D
figure
[X,Y,Z] = cylinder(0.5*D_Array);
[m,n] = size(Z);
C = zeros(m,n);
for i = 1:1:n
    C(:,i) = T_FC_Array;
end
eng = surf(X,Y,L_Engine*(-1*Z+1),C);
eng.EdgeColor = 'none';
axis equal
cbar = colorbar;
cbar.Label.String = 'Kelvin';
title('Film Coolant Temperature in Engine')
xlabel('Meters')
ylabel('Meters')
zlabel('Meters')

%Heat Flux per Area
figure
yyaxis left
plot(L_Array,qa_Array,'b-');
ylabel('Heat Flux per Area (W/m2)')
grid on
yyaxis right
plot(L_Array,D_Array/2,'r-');
ylabel('Engine Radius (m)')
title('Heat Flux per Area VS. Length')
xlabel('Length (m)')
xlim([0,L_Engine])

%Film Cooling Temperature
figure
yyaxis left
plot(L_Array,T_FC_Array,'b-');
ylabel('Film Coolant Temperature (K)')
grid on
yyaxis right
plot(L_Array,D_Array/2,'r-');
ylabel('Engine Radius (m)')
title('Film Coolant Temperature VS. Length')
xlabel('Length (m)')
xlim([0,L_Engine])
%----

%Display ----
Display = ['Maximum Film Coolant Temperature: ',num2str(T_FC_Max),' K'];
disp(Display)
Display = ['Film Coolant Temperature at Throat: ',num2str(T_FC_t),' K'];
disp(Display)
%----

disp('Done!')

%----------------
%Functions
%----------------
function [M_Guess] = AreaMachFunc(G,gamma,SuperSonic,Accuracy)
M_Step = 10^(-1*Accuracy);
M_Guess = 1;
Delta1 = 10^99;
Loop = 1;

if SuperSonic == 0
    M_Guess = M_Step;
end

G1 = ((gamma+1)/2)^(-1*(gamma+1)/(2*gamma-2));
while Loop > 0
    G2 = (1+0.5*(gamma-1)*M_Guess^2)^((gamma+1)/(2*gamma-2));
    G_guess = G1*G2/M_Guess;
    Delta2 = abs(G_guess - G);
    if (Delta2 > Delta1)
        M_Guess = M_Guess - M_Step;
        Loop = 0;
    else
        M_Guess = M_Guess + M_Step;
        Delta1 = Delta2;
        if (SuperSonic == 0)&&(M_Guess >= 1)
            %disp('Error! Left subsonic bounds!')
            M_Guess = 1;
            Loop = 0;
        end
    end
end

end