%% Project Info
%Mason Averill
%ME-545: Advanced FEA
%Spring 2022
%Main Project: Radiator Heat Transfer Analysis 
%% Clear the Workspace and Command Window Before Running
clear;
clc;
%% User Input: Volume Flow Rate and Radiator Configuration
total_volume_flow_rate=35*0.2*6.309*10^-5; %units: m^3/s, estimating based on empirical evidence for automotive engines. #hp*0.2=GPM
total_number_of_channels_per_radiator=10;
total_number_of_radiators=2;
%% User Input: Geometry of Radiator
%Channel Parameters
channel_width=1.6/1000; %units: m (as viewed from top, x-direction) X
channel_length=26/1000; %units: m (as viewed from top, y-direction) X
channel_depth=198.6/1000; %units: m (as viewed from top, z-direction) X
channel_wall_thickness=0.2/1000; %units: m X
distance_between_channels_on_center=8.5/1000; %units: m 

%Fin Parameters
fin_thickness=0.1/1000; %units: m 
distance_between_fins_on_center=0.65/1000; %units: m 

%% User Input: Fin and Channel Thermal Conductivity
k_fin=170; %units: W/m*k
k_channel=170; %units: W/m*k
%% User Input: Ambient Air Properites
air_speed=13; %units: m/s
T_inf=25; %units: C

%All Properites currently evaluated at average of 80C inlet temperature and T_inf of 25C (film temperature)
air_prandtl_number=0.705; %
k_air=28.5/1000; %units: W/m*k
air_density=1.093; %units: kg/m^3
air_dynamic_viscosity=19.7*10^-6; %units: Pa*s
%% User Input: Water Inlet Properties
%All Properites currently evaluated at average of 80C inlet temperature and T_inf of 25C (film temperature)
water_inlet_temperature=80; %units: C
water_dynamic_viscosity=0.355*10^-3; %units: Pa*s
water_density=970; %units: kg/m^3
k_water=667.02/1000; %units: W/m*k
c_p_water=4186; %units: J/kg*k

%% Calculated Fin Geometry
fin_width=(distance_between_channels_on_center-2*channel_wall_thickness-channel_width)/2; %units: m
fin_length=channel_length+2*channel_wall_thickness; %units: m
distance_between_fins_face_to_face=distance_between_fins_on_center-fin_thickness;
fin_cross_sectional_area=fin_thickness*fin_length;
fin_perimeter=2*(fin_thickness+fin_length);
%% Channel Interior Calculated Parameters
%Channel Areas and Perimeters
channel_cross_sectional_area=channel_length*channel_width; %units: m^2
channel_perimeter=2*(channel_width+channel_length); %units: m

%Flow Properties
%water_volume_flow_rate_per_channel=total_volume_flow_rate/(total_number_of_channels_per_radiator*total_number_of_radiators); %units: m^3/s
water_volume_flow_rate_per_channel=5*10^-6;

water_mass_flow_rate_per_channel=water_volume_flow_rate_per_channel*water_density; %units: kg/s
water_average_velocity=water_volume_flow_rate_per_channel/channel_cross_sectional_area; %units m/s
water_time_in_channel=channel_depth/water_average_velocity; %units: seconds
%% Interior of Channel Convection Coefficient Calculations
%Interior of Channel-Treating as Rectangular Duct
channel_hydraulic_diameter=4*channel_cross_sectional_area/channel_perimeter; 
Re_d_channel_interior=(water_density*water_average_velocity*channel_hydraulic_diameter)/water_dynamic_viscosity;
Pr_channel_interior=water_dynamic_viscosity*c_p_water/k_water;
hydrodynamic_entry_length_channel_interior=0.05*Re_d_channel_interior*channel_hydraulic_diameter;

%Laminar
if(Re_d_channel_interior<10000)
     
    if(channel_depth>hydrodynamic_entry_length_channel_interior)
     Nu_d_channel_interior=3.66; %Equation 8.55, Intro to Ht Transfer, Sixth Edition, Bergman constant temp
     h_channel_interior=Nu_d_channel_interior*k_water/channel_hydraulic_diameter;
    end
    
    if(channel_depth<=hydrodynamic_entry_length_channel_interior)
     Nu_d_channel_interior=2.25; %Slug flow, approximation of Nud
     h_channel_interior=Nu_d_channel_interior*k_water/channel_hydraulic_diameter;
    end

end 

%Turbulent
if(Re_d_channel_interior>=10000)
Nu_d_channel_interior=0.023*Re_d_channel_interior^(4/5)*Pr_channel_interior^0.3; %Equation 8.60, Intro to Ht Transfer, Sixth Edition, Bergman
h_channel_interior=Nu_d_channel_interior*k_water/channel_hydraulic_diameter;
end 

%% Exterior of Channel Convection Coefficient Calculations

%Treating Front/Rear of Fins and Channel as Modified Cylinder in Cross Flow: Thin Plate Perpendicular to Flow (Equation 7.44, Table 7.3 Intro to Ht Transfer, Sixth Edition, Bergman) 
Re_d_channel_exterior_cross_flow=(air_density*air_speed*channel_depth/air_dynamic_viscosity);
Re_d_fin_cross_flow=(air_density*air_speed*(fin_width)/air_dynamic_viscosity);

Nu_d_channel_exterior_cross_flow_front=0.667*Re_d_channel_exterior_cross_flow^0.5*air_prandtl_number^(1/3);
Nu_d_channel_exterior_cross_flow_rear=0.191*Re_d_channel_exterior_cross_flow^0.667*air_prandtl_number^(1/3);

Nu_d_fin_cross_flow_front=0.667*Re_d_fin_cross_flow^0.5*air_prandtl_number^(1/3);
Nu_d_fin_cross_flow_rear=0.191*Re_d_fin_cross_flow^0.667*air_prandtl_number^(1/3);

h_channel_exterior_cross_flow_front=k_air*Nu_d_channel_exterior_cross_flow_front/channel_depth;
h_channel_exterior_cross_flow_rear=k_air*Nu_d_channel_exterior_cross_flow_rear/channel_depth;

h_fin_cross_flow_front=k_air*Nu_d_fin_cross_flow_front/fin_width;
h_fin_cross_flow_rear=k_air*Nu_d_fin_cross_flow_rear/fin_width;

%Treating Exterior Edges and Fin as Flat Plate
x_vector_flat_plate=fin_length*0.1:fin_length*10^-3:fin_length*0.99;

Re_x_flat_plate=(air_density*air_speed*x_vector_flat_plate)/air_dynamic_viscosity; 

%Laminar Flow
if(max(Re_x_flat_plate)<5*10^5)
Nu_x_flat_plate=0.332.*Re_x_flat_plate.^0.5.*air_prandtl_number^(1/3); %Equation 7.21, Intro to Heat Transfer 6th ed. Bergman
h_x_flat_plate=k_air.*Nu_x_flat_plate./x_vector_flat_plate;
h_flat_plate_bar=trapz(x_vector_flat_plate,h_x_flat_plate)/(max(x_vector_flat_plate)-min(x_vector_flat_plate));
end

%Turbulent Flow
if(max(Re_x_flat_plate)>=5*10^5)
Nu_x_flat_plate=0.0296.*(Re_x_flat_plate).^(4/5).*air_prandtl_number^(1/3); %Equation 7.30, Intro to Heat Transfer 6th ed. Bergman
h_x_flat_plate=k_air.*Nu_x_flat_plate./x_vector_flat_plate;
h_flat_plate_bar=trapz(x_vector_flat_plate,h_x_flat_plate)/(max(x_vector_flat_plate)-min(x_vector_flat_plate));
end


h_channel_sides=h_flat_plate_bar;
h_fin=h_flat_plate_bar;
%% Discretization of the Domain in Z (along water flow direction): Can adjust resolution here
%User Input
number_of_elements_between_fins=20;

%Calculated Parameters
delta_z_no_fin=(distance_between_fins_on_center-fin_thickness)/number_of_elements_between_fins;
delta_z_fin=fin_thickness;
%% Initialization of Global Information Matrix
global_information_matrix=[];
global_information_matrix(1,1)=0;
global_information_matrix(1,2)=0;
global_information_matrix(1,3)=0;
fin_case=-1;
i=1;
number_of_fins=0;
while(max(global_information_matrix(:,2))<(channel_depth-distance_between_fins_on_center))
    
    if(fin_case==-1 && max(global_information_matrix(:,2))<(channel_depth-distance_between_fins_on_center))
    [row,~]=size(global_information_matrix);
        while(i<number_of_elements_between_fins)
        global_information_matrix(row+i,2)=global_information_matrix(row,2)+i*delta_z_no_fin;
        global_information_matrix(row+i,3)=0;
        i=i+1;
        end
        fin_case=1;
        i=1;
    end

    if(fin_case==0 && max(global_information_matrix(:,2))<(channel_depth-distance_between_fins_on_center))
    [row,~]=size(global_information_matrix);
        while(i<=number_of_elements_between_fins)
        global_information_matrix(row+i,2)=global_information_matrix(row,2)+i*delta_z_no_fin;
        global_information_matrix(row+i,3)=0;
        i=i+1;
        end
        fin_case=1;
        i=1;
    end
   
    if(fin_case==1 && max(global_information_matrix(:,2))<(channel_depth-distance_between_fins_on_center))
    [row,~]=size(global_information_matrix);
        while(i<=1)
        global_information_matrix(row+i,2)=global_information_matrix(row,2)+i*delta_z_fin;
        global_information_matrix(row+i,3)=1;
        i=i+1;
        end
        fin_case=0;
        i=1;
        number_of_fins=number_of_fins+1;
    end

    if(max(global_information_matrix(:,2))>=(channel_depth-distance_between_fins_on_center))
    [row,~]=size(global_information_matrix);
    i=1;
        while(max(global_information_matrix(:,2))<=channel_depth-delta_z_no_fin)
        global_information_matrix(row+i,2)=global_information_matrix(row,2)+i*delta_z_no_fin;
        global_information_matrix(row+i,3)=0;
        i=i+1;
        end
    end

end 
[row,~]=size(global_information_matrix);
for i=1:row
    global_information_matrix(i,1)=i;
end 

%% Thermal Resistance Calculations
%Rear
R_th_rear_no_fin_combined=(1/(h_channel_interior*channel_width*delta_z_no_fin))+(1/(h_channel_exterior_cross_flow_rear*channel_width*delta_z_no_fin));
R_th_rear_fin_channel_only=(1/(h_channel_interior*channel_width*delta_z_fin));
R_th_rear_fin_combined=(1/(h_channel_interior*channel_width*delta_z_fin))+(1/(h_fin_cross_flow_rear*channel_width*delta_z_fin));

%Front
R_th_front_no_fin_combined=(1/(h_channel_interior*channel_width*delta_z_no_fin))+(1/(h_channel_exterior_cross_flow_front*channel_width*delta_z_no_fin));
R_th_front_fin_channel_only=(1/(h_channel_interior*channel_width*delta_z_fin));
R_th_front_fin_combined=(1/(h_channel_interior*channel_width*delta_z_fin))+(1/(h_fin_cross_flow_front*channel_width*delta_z_fin));

%Side
R_th_side_no_fin=(1/(h_channel_interior*channel_length*delta_z_no_fin))+(1/(h_channel_sides*channel_length*delta_z_no_fin));
%% Fin Heat Transfer Preparation Calculations
m=sqrt((h_fin*fin_perimeter)/(k_fin*fin_cross_sectional_area));
%Assuming adiabatic at end of fin (should be used for symmetry exploited)
fin_ht_without_M=tanh(m*fin_width);
%% Global Information Matrix Heat Transfer and Temperature Calculations
Q_front=0;
Q_rear=0;
Q_sides=0;
Q_rejected=0;
T_channel_front=0;
T_channel_rear=0;
M=0;
global_information_matrix(1,4)=water_inlet_temperature;
[row,~]=size(global_information_matrix);
i=1;
while(i<row)

    if(global_information_matrix(i,3)==0) %no fin ht case
    %Compute Q's:
    Q_front=(global_information_matrix(i,4)-T_inf)/R_th_front_no_fin_combined;
    Q_rear=(global_information_matrix(i,4)-T_inf)/R_th_rear_no_fin_combined;
    Q_sides=2*(global_information_matrix(i,4)-T_inf)/R_th_side_no_fin;
    Q_rejected=Q_front+Q_rear+Q_sides;

    %Update Global Information Matrix
    global_information_matrix(i+1,4)=global_information_matrix(i,4)-Q_rejected/(water_mass_flow_rate_per_channel*c_p_water);
    global_information_matrix(i,5)=Q_front;
    global_information_matrix(i,6)=Q_rear;
    global_information_matrix(i,7)=Q_sides;
    global_information_matrix(i,8)=Q_rejected;
    end

    if(global_information_matrix(i,3)==1) %fin ht case
    Q_front=(global_information_matrix(i,4)-T_inf)/R_th_front_fin_combined;
    Q_rear=(global_information_matrix(i,4)-T_inf)/R_th_rear_fin_combined;
    T_channel_rear=global_information_matrix(i,4)-Q_rear*R_th_rear_fin_channel_only;
    T_channel_front=global_information_matrix(i,4)-Q_front*R_th_front_fin_channel_only;


    M=sqrt(h_fin*fin_perimeter*k_fin*fin_cross_sectional_area)*((T_channel_rear+T_channel_front)/2-T_inf);
    if(~isreal(M))
    M=0;
    end
    Q_sides=2*M*fin_ht_without_M;


    if(global_information_matrix(i,4)<T_inf)
    Q_sides=-2*M*fin_ht_without_M;
    end
    Q_rejected=Q_front+Q_rear+Q_sides;

    %Update Global Information Matrix
    global_information_matrix(i+1,4)=global_information_matrix(i,4)-Q_rejected/(water_mass_flow_rate_per_channel*c_p_water);
    global_information_matrix(i,5)=Q_front;
    global_information_matrix(i,6)=Q_rear;
    global_information_matrix(i,7)=Q_sides;
    global_information_matrix(i,8)=Q_rejected;

    end
i=i+1;
end

%% Global Information Matrix Labeling (String Version for Printing)
global_information_matrix_label_first_row=["Element Number" "Distance Along Channel (m)" "Fin Status" "Bulk Temperature of Water (C)" "Q Front (W)" "Q Rear (W)" "Q Sides (W)" "Q Total (W)"];
global_information_matrix_labeled=[global_information_matrix_label_first_row;cellstr(string(global_information_matrix))];
%% Format User Output

[row,~]=size(global_information_matrix);
water_outlet_temperature=global_information_matrix(row, 4);
water_delta_T=water_inlet_temperature-water_outlet_temperature;

total_heat_loss_per_channel=water_mass_flow_rate_per_channel*c_p_water*(water_inlet_temperature-water_outlet_temperature);
total_heat_loss_per_channel_check=sum(global_information_matrix(:,8));


%Gather Relevant Model Parameters and Categorize
system_calculation_parameters=[total_number_of_radiators,total_number_of_channels_per_radiator, total_volume_flow_rate, water_inlet_temperature, air_speed, T_inf];
geometry_calculation_parameters=[channel_width*1000, channel_length*1000, channel_wall_thickness*1000, channel_depth, fin_width*1000, distance_between_fins_on_center*1000, fin_thickness*1000, number_of_fins];
flow_calculation_parameters=[water_mass_flow_rate_per_channel, water_volume_flow_rate_per_channel, water_average_velocity, water_time_in_channel];
convection_coefficient_calculation_parameters=[h_channel_interior, h_channel_exterior_cross_flow_front, h_channel_exterior_cross_flow_rear, h_fin_cross_flow_front,h_fin_cross_flow_rear, h_channel_sides, h_fin];
other_thermo_calculation_parameters=[k_channel, k_fin];
resolution_calculation_parameters=[number_of_elements_between_fins,row];
final_result_calculation_parameters=[water_inlet_temperature, water_outlet_temperature, water_delta_T, water_inlet_temperature*(9/5)+32, water_outlet_temperature*(9/5)+32, water_delta_T*1.8;];

%Construct Formatted Output for each Categorized Model Parameter
system_calculation_parameters_formatted=sprintf('Total Number of Radiators=%i \nTotal Number of Channels=%i \nTotal Volume Flow Rate=%.3e m^3/s \nWater Inlet Temperature=%.1f C\nVehicle Speed=%.1f m/s \nAmbient Temperature=%.1f C\n',system_calculation_parameters);
geometry_calculation_parameters_formatted=sprintf('Coolant Passage Channel Width=%.1f mm \nCoolant Passage Channel Length=%.1f mm \nCoolant Passage Wall Thickness=%.2f mm\nCoolant Passage Channel Total Depth=%.1f m\nFin Width=%.2f mm\nDistance Between Fins on Center=%.2f mm\nFin Thickness=%.2f mm\nNumber of Fins on Each Side of the Channel=%i\n',geometry_calculation_parameters);
flow_calculation_parameters_formatted=sprintf('Water Mass Flow Rate Per Channel=%.3e kg/s \nWater Volume Flow Rate Per Channel=%.3e m^3/s \nWater Average Velocity=%.3f m/s \nTime for Water to Pass Through Channel=%.2f s\n',flow_calculation_parameters);
convection_coefficient_calculation_parameters_formatted=sprintf('Interior Convection Coefficient of Channel=%.1f W/m^2*k \nFront Exterior Convection Coefficient of Channel=%.1f W/m^2*k \nRear Exterior Convection Coefficient of Channel=%.1f W/m^2*k \nFront Convection Coefficient of Fins=%.1f W/m^2*k \nRear Convection Coefficient of Fins=%.1f W/m^2*k \nSides Exterior Convection Coefficient of Channel=%.1f W/m^2*k \nFin Convection Coefficient=%.1f W/m^2*k\n',convection_coefficient_calculation_parameters);
other_thermo_calculation_parameters_formatted=sprintf('Thermal Conductivity of the Channel Walls=%.1f W/m*k \nThermal Conductivity of the Fins=%.1f W/m*k\n',other_thermo_calculation_parameters);
resolution_calculation_parameters_formatted=sprintf('Number of Elements Considered Between Each Fin=%i \nNumber of Nodes Considered Along Water Flow Direction=%i\n',resolution_calculation_parameters);
final_result_calculation_parameters_formatted=sprintf('OVERALL RESULTS:\nWater Inlet Temperature=%.2f C \nWater Outlet Temperature=%.2f C \nWater Change in Temperature=%.2f C \nWater Inlet Temperature=%.2f F \nWater Outlet Temperature=%.2f F \nWater Change in Temperature=%.2f F\n',final_result_calculation_parameters);
output_label=sprintf('%s\n','Model Inputs and Results (can also be viewed in global_information_matrix_labeled variable):');

assembled_formatted_model_parameters_and_overall_results_1=[output_label system_calculation_parameters_formatted geometry_calculation_parameters_formatted flow_calculation_parameters_formatted convection_coefficient_calculation_parameters_formatted other_thermo_calculation_parameters_formatted resolution_calculation_parameters_formatted final_result_calculation_parameters_formatted];
assembled_formatted_model_parameters_and_overall_results_2=splitlines(convertCharsToStrings(assembled_formatted_model_parameters_and_overall_results_1));
assembled_formatted_model_parameters_and_overall_results_3=[system_calculation_parameters_formatted geometry_calculation_parameters_formatted flow_calculation_parameters_formatted convection_coefficient_calculation_parameters_formatted other_thermo_calculation_parameters_formatted resolution_calculation_parameters_formatted final_result_calculation_parameters_formatted];

[row2,~]=size(assembled_formatted_model_parameters_and_overall_results_2);

for i=1:row2
    global_information_matrix_labeled(i,9)=assembled_formatted_model_parameters_and_overall_results_2(i,1);
end

for i=row2:row+1
    global_information_matrix_labeled(i,9)="";
end

bulk_temperature_poly_coeffs=polyfit(global_information_matrix(:,2),global_information_matrix(:,4),2);

bulk_temperature_poly=polyval(bulk_temperature_poly_coeffs,global_information_matrix(:,2));

%% Display User Output 
%Print Model Parameters and Results to the Command Window
disp(assembled_formatted_model_parameters_and_overall_results_1);

% %{

figure(1)
plot(global_information_matrix(:,2),bulk_temperature_poly)
%plot(global_information_matrix(:,2),global_information_matrix(:,4))
xlabel('Distance Along Channel In Direction of Water Flow (m)')
ylabel('Bulk Temperature of Water (C)')
title('Bulk Temperature of Water vs Distance Along Channel')
xlim([0 max(global_information_matrix(:,2))])
ylim([(floor(min(global_information_matrix(:,4)))) water_inlet_temperature])
text(0.05*max(global_information_matrix(:,2)),water_inlet_temperature-0.65*(water_inlet_temperature-min(global_information_matrix(:,4))),assembled_formatted_model_parameters_and_overall_results_3)

% %}

%{
hold on
plot(global_information_matrix(:,2),global_information_matrix(:,4))
hold off
%}


 % %{
figure(2)
area(global_information_matrix(:,2),global_information_matrix(:,3)*fin_width)
xlabel('Distance Along Channel In Direction of Water Flow (m)')
ylabel('Fin Width Extending Towards Other Channels (m)')
title('Geometric Representation of Fins')
xlim([0 channel_depth])

 % %}