## Thermal analysis of box-type solar cookers (2/2) The temperature plot of solar cooker to validate the mathematical model. In a practical situation, once the cooking fluid reaches 100 ̊ C, cooking fluid temperature remains constant because of latent heat of vaporization.

## MATLAB script

```% Author: Jonah Kadoko
% Project: Solar Fundamentals, Fall 2014
% Instructor:

% Description: This is the main script for modelling a solar box-type
% cooker. In this mathematical model, several coatings of commercially
% avialable coatings will be i nvestigated.

% T : is a 5x1 array of temperatures in the following order:
% T(1): Temperature of the lid of vessel, Tc
% T(2): Temp of the sides and base of the vessel, Tsb
% T(3): Temp of the cooking fluid, Tf
% T(4): Temp of the lower glass cover, Tgl
% T(5): Temp of upper glass cover, Tgu

% T(1) -
% Input Variables
clc
t_S=6*3600; % [s] starting time
t_F=12*3600; % [s] time to stop the iterations, 4.770 seems to be max time
% t_Span=linspace(0,t_F,floor(t_F/t_Step));
% Calculations
options=odeset('Stats','off'); % set the end conditions
[t, T]=ode45(@temp_Profile,[0 t_F],[75 80 80 25 25],options,t_S);
```

## Plot Results

```clf; % clear all figures

axes('FontSize',14);
T_max=max(max(T));
axis([t_S/3600 (t_F+t_S)/3600 0 T_max])
plot((t+t_S)/3600,T(:,1),':r','LineWidth',2);
hold on;
plot((t+t_S)/3600,T(:,2),'--g','LineWidth',2);
plot((t+t_S)/3600,T(:,3),'-.k','LineWidth',2);
plot((t+t_S)/3600,T(:,4),'-','LineWidth',2);
plot((t+t_S)/3600,T(:,5),'-k','LineWidth',2);
grid on;
title('Temperature Analaysis of Solar Cooker: Model Validation','FontWeight','bold','FontSize',20)
xlabel('Time (hrs)','FontWeight','bold','FontSize',16)
ylabel('Temperature (deg C)','FontWeight','bold','FontSize',16)

legend('Temp of the lid of vessel, Tc','Temp of the sides and base of the vessel, Tsb',...
'Temp of the cooking fluid, Tf','Temp of the lower glass cover, Tgl',...
'Temp of upper glass cover, Tgu');
hold off;
```

## Calculate the glazing properties

```a_0=0.05;
theta_1=60*2*pi/360;
theta_2=asin(sin(theta_1)/1.52); % the angle of the refraction
a=a_0/cos(theta_2);
rho=0.5*((sin(theta_2-theta_1)/sin(theta_2+theta_1))^2)*(1+cos(theta_1+theta_2)^2/cos(theta_1-theta_2)^2)
trans=(1-rho)^2*(1-a)^2/(1-rho^2*(1-a)^2)
```

## ODE Solver

```function dT_dt=temp_Profile(t,T,t_S)
% this is the ode for temperature distribution in the box-type solar cooker
% Variables:
% T : is a 5x1 array of temperatures in the following order:
% T(1): Temperature of the lid of vessel, Tc
% T(2): Temp of the sides and base of the vessel, Tsb
% T(3): Temp of the cooking fluid, Tf
% T(4): Temp of the lower glass cover, Tgl
% T(5): Temp of upper glass cover, Tgu
% t_S : Starting time in seconds, variable used to offset time from the ode
% time to solar time

% Initialization
dT_dt=zeros(5,1);

% Temperature
load 'Lawrence_TMY3_June_21'; % Site data, TMY(1)=time(hrs); TMY(2)=Global horizontal Irradiation
% TMY(3)=diffuse horizontal irradiation, TMY(4)=dry bulb temperature
if ((t+t_S)/3600)<1
t_i=1; % [hrs] time of day
else
t_i=(t_S+t)/3600;
end
T_a=interp1(Lawrence_TMY3(:,1),Lawrence_TMY3(:,4),t_i); %[C] temperature of the ambient air

% Absorber plate: gray body base, 'base' is the same as 'absorber'
% Box material ????
k=0.03; %[W/mK] thermal conductivity of insulating material
l_s=30*.0254; %[m] length of one of the sides of box
w_s=30*0.0254; %[m] width of one of the sides of box
h_s=12*0.0254; %[m] height of box
x_b=1*0.0254; %[m] thickness of bottom insulating material
x_s=3.25*0.0254; %[m] average thickness of insulating material on the sides
n_opt=0.6; %[] efficience of inner reflector
n_u=0.3; %[] overal utilization of cooker

% Cooking Vessel: assum a cylindrical vessel
% User entered data
a_p=0.09; %[] absorptivity of cooking vessel

E_c=1; %[] emissivity of the cooking vessel cover/lid
r_v=3.75*0.0254; %[m] radius of the cooking vessel
h_v=5*0.0254; %[m] height of the cooking vessel
x_t=0.125*0.0254; %[m] thickness of the cooking vessel
A_b=pi*r_v^2; %[m2] area of the base
A_sb=2*pi*r_v*h_v+A_b; %[m2] area of the sides + base of cooking vessel
A_s=2*pi*r_v*h_v; %[m2] area of the sides cooking vessel
A_t=A_b; %[m2] area of lid of vessel
h_ccl=convection_Air(T(1),T(4),(r_v/2),'Horizontal-Upper'); %(estimated=4)[W/m2K] convective heat transfer: wind
%(estimated=3.777; todo list)[W/m2K] convective heat transfer: vessel to lower glass
L_c_sbf=38*0.0254/(pi*3.75*2)+4.875*0.0254; % this is the height of vessel + charecteristic length of the base
h_c_sbf=500;
% h_c_sbf=convection_Water(T(2),T(3),L_c_sbf,'Vertical'); %(estimated=5)[W/m2K] convective heat transfer: walls to cooking fluid
% ToDo - this convective heat coefficient does not take into account the
% fact that some heat is transmitted through the base and some through the
% side walls of the vessel. The characteristic length is calculated by
% deviding the total surface area of the cooking vessel by the total
% perimeter

k_t=80; %(cast iron)[W/mK] thermal conductivity of cooking vessel material
U_b=k/x_b; %[W/m2K] bottom loss coefficient
% U_s=1/((x_s/k)+1/h_w); %[W/m2K] side loss coefficient
U_t=5; %(estimated)[W/m2K] top loss coefficient

m_c=0.531; %[kg] mass of the lid (added by JK)
c_c=506; % [J/kgK] heat capacity of lid (added by JK)
m_sb=5.675*0.454; %[kg] mass of sides and base of vessel(added by JK)
c_sb=c_c; %[J/kgK] heat capacity of base+surfaces (added by JK)

% Cooking fluid properties: Assume the food cooked is just water
h_f=3*0.0254; %[m] height of cooking fluid in the cooking vessel
rho_f=1000; %[kg/m^3] 1000 density of the cooking fluid
c_f=4170; %[J/(kg.K)], 4170 specific heat capacity of cooking fluid

m_f=(pi*(r_v-x_t)^2)*h_f*rho_f; % [kg] mass of cooking fluid
A_sf=2*pi*(r_v-x_t)*h_f; %[m2] surface area of vessel covered by cooking fluid
h_gap=3.17; %[W/m2K] convective heat transfer: fluid to ambient
R_t=(1/h_gap)+(x_t/k_t)+(1/U_t); % thermal resistance: fluid surface and ambient

% Cover/Aperture: assume double glaze
% User entered variablesf
x_g=0.1875*0.0254; %[m] 0.1875*0.0254 thickness of the cover
a_g=0.05; %[] absorptivity of glass cover
t_g=0.95; %(estimated)[] 0.95 transmissivity of glass
density_g=2180; %[kg/m3] 2180 density of glass

r=30*0.0254; %[m] width of the cover
l_g=30*0.0254; %[m] length of the cover
x_lu=0.5*0.0254; %[m] the air gap between upper and lower glass
E_gl=1; % [] emissivity of lower glass
E_gu=1; % [] emissivity of lower glass
A_g=l_g*r; % [m2] area of glass
h_clu=convection_Air(T(4),T(5),(r^2/(4*r)),'Horizontal-Upper'); %(estimated=3)[W/m2K] convective heat transfer: lower glass to upper glass
h_w=convection_Air(T(5),T_a,(r^2/(4*r)),'Horizontal-Upper'); %(estimated=4)[W/m2K] convective heat transfer: wind
h_rlu=4; %(estimated)[W/m2K] radiative heat transfer: lower glass to upper glass
h_rua=4; %(estimated)[W/m2K] radiative heat transfer: upper glass to ambient
k_g=1.38; %[W/mK] thermal conductivity of glazing material (added by JK)
m_g=(x_g*r*l_g)*density_g; %[kg] mass of glass
c_g=750; %[J/(kg.K)], specific heat capacity of glass

% Reflector: the reflector outside the box
% user entered values
rho_o=.85; %[] 0.85 outer mirror reflectivity
tilt=80*2*pi/360; %[degrees] 50 angle of outer reflector, convert to radians for

% compatibility  with MATLAB 2013+
c=30*0.0254; %[m] length of reflector
w_r=c; %[m] width of reflector
t_r=0.05;%[m] (estimated)thickness of reflector (added by JK)
A_r=c*w_r; %[m] area of reflector
s=sqrt(c^2+r^2-2*c*r*cos(tilt)); %[m] distance from the outer mirror to outer edge of cooker
F_rg=(c+r-s)/2*r; %[] this is the view factor of the reflector-glass
A_c=A_r+A_g; %[m2] aperture size is the total of A_reflector + A_window

% Irradiance and environment: Assume that intensity is constant throughout the day
h_angle=(t_i-12)*(2*pi/24); %[degrees] hour angle, starting at 12 noon

Steph_Boltz=5.67*10^-8; %[W/m^2*K^4] Stephano-boltzmann constant
I_h=interp1(Lawrence_TMY3(:,1),Lawrence_TMY3(:,2),t_i); %[W/m2] global irradiance on a horizontal surface, Lawrence, MA, June 21, TMY3
I_d=interp1(Lawrence_TMY3(:,1),Lawrence_TMY3(:,3),t_i); %(estimated)[W/m2] diffuse irradiance on a horizontal surface, Lawrence, MA, June 21, TMY3
I_b=I_h-I_d; % (estimated)[W/m2] beam irradiance on a horizontal surface, Lawrence, MA, June 21, TMY3

rho=0.2; %[] ground reflectivity
rho_i=0.85; %[] inner mirror reflectivity
lat=42.717*2*pi/360; %[degrees] latitude of site
decl=23.45*2*pi/360; %[degrees] solar declination angle, June 21, Lawrence, MA

R_b=(cos(lat-tilt)*cos(decl)*cos(h_angle)+sin(lat-h_angle)*sin(decl))/(cos(lat)*cos(decl)*cos(h_angle)+sin(lat)*sin(decl));
I_t=I_b*R_b+0.5*I_d*(1+cos(tilt))+0.5*rho*(I_b+I_d)*(1-cos(tilt)); %[W/m2] global irradiance on outer reflector
I_ro=I_t*rho_o*F_rg*A_r/A_g; %[W/m2] irradiance reflected by outer reflector
I_ri=(I_ro+I_h)*(t_g^2)*rho_i*n_opt; %[W/m2] irradiance reflected by inner reflector

% Calculations of variables - these are the variables that depend on other
% variables and temperature
U_s=1/((x_s/k)+1/h_w); %[W/m2K] side loss coefficient

% Calculating the view factor from cooking vessel cover to fluid
% Assume circular disks, page 817 of Fundamentals of Heat....
L_cf=h_v-h_f; % [m] the distance between the cooking fluid and the cooking vessel cover
r_c=r_v/L_cf;r_f=(r_v-x_t)/L_cf; % [] exactly from fundamentals of heat transfer book (check references)
S=1+(1+r_f^2)/r_c^2;
F_cf=0.5*(S-sqrt(S^2-4*(r_f/r_c)^2)); % [] view factor F

% Calculating the view factor from cooking vessel cover to lower glass
% Assume circular disks, page 817 of Fundamentals of Heat....
L_cl=7*0.0254; % [m] the distance between the top of the vessel to the lower glass
r_c2=r_v/L_cl;r_l=r/L_cl;
S=1+(1+r_l^2)/r_c2^2;
F_cl=0.5*(S-sqrt(S^2-4*(r_l/r_c2)^2)); % [] view factor F

% Calculating the view factor from the lower to upper glass
% assume two rectangles parallel, page 817 of Fundamentals of Heat ..
x_gl=r/x_lu;y_gu=r/x_lu;
F_lu=(2/(pi*x_gl*y_gu))*(log(sqrt((1+x_gl^2)*(1+y_gu^2)))+...
x_gl*sqrt(1+y_gu^2)*atan(x_gl/(sqrt(1+y_gu^2)))+...
y_gu*sqrt(1+x_gl^2)*atan(y_gu/(sqrt(1+x_gl^2)))+....
-x_gl*atan(x_gl)-y_gu*atan(y_gu));

% Calculations of the defferential equations
dT_dt(1)=((I_ro+I_h)*(t_g^2)*a_p*A_t-U_t*A_t*(T(1)-T_a)-Steph_Boltz*E_c*F_cf*A_t*(T(1)-T(3)))/(m_c*c_c);
dT_dt(2)=(I_ri*a_p*(A_c-A_t)-h_c_sbf*A_sf*(T(2)-T(3))-(U_b*A_b+U_s*A_s)*(T(2)-T_a))/(m_sb*c_sb);
dT_dt(3)=(Steph_Boltz*E_c*F_cf*A_t*(T(1)-T(3))+h_c_sbf*A_sf*(T(2)-T(3))-(A_t/R_t)*(T(3)-T_a))/(m_f*c_f);
dT_dt(4)=((I_ro+I_h)*t_g*a_p*A_g+Steph_Boltz*E_c*F_cl*A_t*(T(1)^4-T(4)^4)+h_ccl*A_t*(T(1)-T(4))-Steph_Boltz*E_gl*F_lu*A_g*(T(4)^4-T(5)^4)-h_clu*A_g*(T(4)-T(5)))/(m_g*c_g);
dT_dt(5)=((I_ro+I_h)*a_g*A_g+Steph_Boltz*E_gl*F_lu*A_g*(T(4)^4-T(5)^4)+h_clu*A_g*(T(4)-T(5))-Steph_Boltz*E_gu*A_g*(T(5)^4-T_a^4)-h_w*A_g*(T(5)-T_a))/(m_g*c_g);
```