clear all;
clc;
 
T_amb=20;                  %Ambient Temperature [C]
 
%Device Properties
P=5;                       %Power Dissipated [W]
rd=1.95e-3;                %Device Radius [m]
Ad=pi*rd^2;                %Device Area [m^2]
 
%Fluid Thermo-Physical Properties
Fluid='Air';               %Fluid Variable (Air or Water)
 
if strcmp(Fluid,'Air')==1
% Air Fluid Properties 
rho=1.1614;                %Density [kg/m^3]
nu=15.84e-6;               %Kinematic Viscosity [m^2/s]
alpha_f=28.86e-6;          %Thermal diffusivity [m^2/s]
k_f=26.1e-3;               %Thermal conductivity [W/(m.K)]
Pr=nu/alpha_f;             %Prandtl number [-]
beta_f=1/(273+20);         %CTE [1/K]
g=9.81;                    %Gravitational Acceleration [m/s^2]
SD_fact=1;                 %Factor for Search Function
 
elseif strcmp(Fluid,'Water')==1
rho=1000;                  %Density [kg/m^3]
nu=1e-6;                   %Kinematic Viscosity [m^2/s]
alpha_f=1.435e-7;          %Thermal diffusivity [m^2/s]
k_f=0.6;                   %Thermal conductivity [W/(m.K)]
Pr=nu/alpha_f;             %Prandtl number [-]
beta_f=0.000276;           %CTE [1/K]
g=9.81;                    %Gravitational Acceleration [m/s^2]
SD_fact=10;                %Factor for Search Function
 
end
 
%Coating Spreader Properties
t=[20 50 (100:100:1000)];  %Coating Thickness [microns]
k=387.6;                   %Coating Thermal Conduc. [W/(m.K)]
kt=t*387.6*1e-6;           %k multiplied by t [W/K]
 
%PCB Properties
tf=1.6e-3;                 %FR-4 Thickness [m]
kf=0.343;                  %FR-4 Thermal Conduc. [W/(m.K)]
 
%Parylene C Properties
R_Par = 120*1e-6;         %Parylene Thermal Insul. [m^2.K/W]
 
%Heat Transfer Initialization
h=100;                     %HTC Guess [W/(m^2.K)]
sqrt_alpha=sqrt(1+kf/(h*kf*R_Par+tf*h+kf)); %sqrt(alpha)
L_char=1e-2;               %Characteristic Length Guess [m]

%Solution Initialization
n=length(kt);              %Number of Solutions
SD=zeros(1,n);             %Spreading Depth Vector
CSD=zeros(1,n);            %Characteristic SD Vector
R=zeros(1,n);              %Thermal Resistance Vector
h_tot=zeros(1,n);          %HTC Vector
Ra_tot=zeros(1,n);         %Rayleigh Number Vector
sqrt_alpha_tot=zeros(1,n); %sqrt(alpha) Vector 
 
%Thickness Loop
for i=1:n
    
    kt_i=kt(i);
    t_i=t(i);
    CSD_i=(kt_i/h)^0.5;    
    error_h=2;             %Error on HTC Initialization
 
%Iterative Solution Loop
while(error_h>0.01)
    
    R_i=besselk(0,sqrt_alpha*rd/CSD_i)/...
        besselk(1,sqrt_alpha*rd/CSD_i)...
        /(2*sqrt_alpha*pi*rd*h*CSD_i);  %Spreader Thermal Res.
    Rd=1/h/(pi*rd^2);      %Device Area Thermal Resistance
    R_i=Rd*R_i/(Rd+R_i);   %Parallel Equivalent Resistance
 
    T_max=T_amb+R_i*P;     %Maximum Temperature [C] (Device)
 
    %HTC Calculation
    Ra=g*beta_f*(T_max-T_amb)*(L_char)^3....
        /nu/alpha_f;       %Rayleigh Number
    
    h_new=k_f/(L_char)*(0.825+0.387*Ra^(1/6)/...
        (1+(0.492/Pr)^(9/16))^(8/27))^2; %New HTC 
   
    error_h=abs(h-h_new)/h*100;          %HTC Error [%]
  
    %Update HTC and Other Variables
    h=h_new;
    sqrt_alpha=sqrt(1+kf/(h*kf*R_Par+tf*h+kf));
    CSD_i=(kt_i/h)^0.5;
    
    %Calculating Spreading Depth
    f=@(x)(besselk(0,sqrt_alpha*x/CSD_i)-...
        0.1*besselk(0,sqrt_alpha*rd/CSD_i));
    roots=fzero(f,CSD_i*SD_fact);
    SD_i=roots(1);
    
    %Updating the Characteristic Length
    L_char=2*SD_i;
        
end
 
%Appending Solutions
h_tot(1,i)=h;
SD(1,i)=SD_i;
CSD(1,i)=CSD_i;
R(1,i)=R_i;
Ra_tot(1,i)=Ra;
sqrt_alpha_tot(1,i)=sqrt_alpha;
 
end
