% Calculate the temperature of a cup of water being heated in a microwave
% oven assuming uniform heat generation in the water using:
%
% 1. One-dimensional model for steay state temperature distribution
% 2. Lumped capcitance method for bulk water temperture vs. time
% 3. PDE toolbox for two-dimensional temperature distribution vs. time for:
% a. Cartesian coordinate system with insulation on left side and
% convection on top and right side
% b. Cylindrical coordinate system with convection on top and side
%
% NOTE: For comparison to 1-D steady state solution, PDE toolbox solution
% can be obtained for no cup bottom and insulation on top and bottom.
%% Define geometry, properties, initial, and boundary condition values
clear, clc, close all
coord = input('Enter 1 for Cartesian or 2 for cylindrical coordinates: ');
bottom = input('Enter 1 for 1-D or 2 for 2-D solution with bottom: ');
q_dot = 1.5e4; % W/m^3, uniform heat generation rate in water
L = 0.080; % m, length of outside of cup
r_1 = 0.035; % m, inner radius of cup
r_2 = 0.040; % m, outer radius of cup
L_c = r_2-r_1; % m, thickness of ceramic cup
k_w = 0.600; % W/m-K, thermal conductivity of water
c_w = 4180; % J/kg-K, heat capcity of water
rho_w = 1000; % kg/m^3, density of water
k_c = 2.0; % W/m-K, thermal conductivity of ceramic cup
c_c = 1080; % J/kg-K, heat capacity of ceramic cup
rho_c = 2300; % kg/m^3, density of ceramic cup
T_i = 20; % deg C, initial temperature
T_inf = T_i; % deg C, temperature of surrounding air
h = 10; % W/m^2-K, convection coefficient for surrounding air
if bottom == 1
t = 0; % m, no bottom
y = L/2; % m, y location at middle of water
edges = 1; % Convection on sides, insulated on top and bottom
else
t = L_c; % m, thickness of bottom
y = (L + t)/2; % m, y location at middle of water
edges = [2,5,6]; % Convection on top and sides, insulated on bottom
end
%% 1. One-dimensional model for steay state temperature distribution
r_w = 0:r_1/20:r_1; % m, x or r locations for water
r_c = r_1:L_c/20:r_2; % m, x or r locations for cup
r = [r_w r_c];
if coord == 1 % Cartesian coordinates
As = 2*pi*r_1*L; % m^2, surface area for convection
Rt_cond = L_c/(k_c*As); % K/W, conduction thermal resistance
Rt_conv = 1/(h*As); % K/W, convection thermal resistance
V_c = L_c*As + 2*pi*r_1^2*t; % m^3, volume of cup
V_w = r_2*As - V_c; % m^3, volume of water
q = q_dot*V_w; % W, total heat transfer
T_2 = T_inf + q*Rt_conv; % deg. C, temperature at r_2
T_1 = T_2 + q*Rt_cond; % deg. C, temperature at r_1
T_0 = T_1 + (q_dot*r_1^2)/...
(2*k_w); % deg. C, temperature at r_0
T_vs_r_c = T_1 + (T_2 - T_1)*(r_c - r_1)/L_c;
T_vs_r_w = T_1 + (T_0 - T_1)*(1 - r_w.^2/r_1^2);
x_str = 'x (m)';
else % Cylindrical coordinates
As = 2*pi*r_2*L; % m^2, surface area for convection
Rt_conv = 1/(h*As); % K/W, conduction thermal resistance
Rt_cond = log(r_2/r_1)/...
(2*pi*k_c*L); % K/W, convection thermal resistance
V_c = pi*(r_2^2-r_1^2)*(L-t)...
+ pi*r_2^2*t; % m^3, volume of cup
V_w = pi*r_1^2*(L-t); % m^3, volume of water
q = q_dot*V_w; % W, total heat transfer
T_2 = T_inf + q*Rt_conv; % deg. C, temperature at r_2
T_1 = T_2 + q*Rt_cond; % deg. C, temperature at r_1
T_0 = T_1 + (q_dot*r_1^2)/...
(4*k_w); % deg. C, temperature at r_0
T_vs_r_c = T_2 + (T_1 - T_2)*log(r_c/r_2)/log(r_1/r_2);
T_vs_r_w = T_1 + (T_0 - T_1)*(1 - r_w.^2/r_1^2);
x_str = 'r (m)';
end
T_vs_r_1D = [T_vs_r_w T_vs_r_c];
%% 2. Lumped capcitance method for bulk water temperture vs. time
tau_w = (rho_w*V_w*c_w)/(h*As); % sec, time constant for water
tau_c = (rho_c*V_c*c_c)/(h*As); % sec, time constant for cup
tau = tau_w + tau_c;
tfinal = 5*tau; % sec, final time
Nt = 60; % number of solution times
tlist = 0:tfinal/Nt:tfinal; % sec, solution times
b_tau = (q_dot*V_w)/(h*As);
T_vs_t_0D = T_i + b_tau*(1 - exp(-tlist/tau));
%% 3. PDE toolbox for two-dimensional temperature distribution vs. time
% Create thermal models for steady and transient caculations
thermalmodelS = createpde('thermal','steadystate');
thermalmodelT = createpde('thermal','transient');
% Create geometry for water and cup and append to model
R1 = [3;4;0;r_2;r_2;0;0;0;L;L];
R2 = [3;4;0;r_1;r_1;0;t;t;L;L];
gd = [R1,R2]; % geometry description matrix
sf = 'R1+R2'; % set formula
ns = char('R1','R2')'; % name space matrix
dl = decsg(gd,sf,ns); % decomposed geometry matrix
geometryFromEdges(thermalmodelS,dl);
geometryFromEdges(thermalmodelT,dl);
% Specify thermal properties of the material and heat source
if coord == 1 % Cartesian coordinates
thermalProperties(thermalmodelS,'Face',1,'ThermalConductivity',k_c);
thermalProperties(thermalmodelS,'Face',2,'ThermalConductivity',k_w);
thermalProperties(thermalmodelT,'Face',1,'ThermalConductivity',k_c, ...
'MassDensity',rho_c,...
'SpecificHeat',c_c);
thermalProperties(thermalmodelT,'Face',2,'ThermalConductivity',k_w,...
'MassDensity',rho_w,...
'SpecificHeat',c_w);
internalHeatSource(thermalmodelS,q_dot,'Face',2);
internalHeatSource(thermalmodelT,q_dot,'Face',2);
else % Cylindrical coordinates with variable transformation
kFunc_w = @(region,state) k_w*region.x;
cFunc_w = @(region,state) c_w*region.x;
kFunc_c = @(region,state) k_c*region.x;
cFunc_c = @(region,state) c_c*region.x;
qFunc = @(region,state) q_dot*region.x;
thermalProperties(thermalmodelS,'Face',1,'ThermalConductivity',kFunc_c);
thermalProperties(thermalmodelS,'Face',2,'ThermalConductivity',kFunc_w);
thermalProperties(thermalmodelT,'Face',1,'ThermalConductivity',kFunc_c, ...
'MassDensity',rho_c,...
'SpecificHeat',cFunc_c);
thermalProperties(thermalmodelT,'Face',2,'ThermalConductivity',kFunc_w,...
'MassDensity',rho_w,...
'SpecificHeat',cFunc_w);
internalHeatSource(thermalmodelS,qFunc,'Face',2);
internalHeatSource(thermalmodelT,qFunc,'Face',2);
end
% Set initial and boundary conditions (default is zero heat flux)
IC = thermalIC(thermalmodelT,T_i);
if coord == 1 % Cartesian coordinates
BCS = thermalBC(thermalmodelS,'Edge',edges,...
'ConvectionCoefficient',h,...
'AmbientTemperature',T_inf);
BCT = thermalBC(thermalmodelT,'Edge',edges,...
'ConvectionCoefficient',h,...
'AmbientTemperature',T_inf);
else % Cylindrical coordinates with varaiable transformation
hFunc = @(region,~) h*region.x;
BCS = thermalBC(thermalmodelS,'Edge',edges,...
'ConvectionCoefficient',hFunc,...
'AmbientTemperature',T_inf);
BCT = thermalBC(thermalmodelT,'Edge',edges,...
'ConvectionCoefficient',hFunc,...
'AmbientTemperature',T_inf);
end
% Generate mesh for maxiumum spacing Hmax
mesh = generateMesh(thermalmodelS,'Hmax',L_c/4);
thermalmodelT
