AME Fall '17 Senior Design
\begin
Discover why 18 million people worldwide trust Overleaf with their work.
%\documentclass[openany]{book}
\documentclass[12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage{amssymb} %for fancy L
\usepackage{calrsfs} %for fancy L
\usepackage{cancel}
\usepackage{tabularx}
\usepackage[hyphens]{url}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage[titletoc,title]{appendix}
\usepackage{subfig}
\DeclareMathAlphabet{\pazocal}{OMS}{zplm}{m}{n} %for fancy L
\usepackage{epsfig, float,array,tabu,longtable,}
\usepackage{hyperref,wrapfig}
\usepackage{enumerate}
\usepackage{graphicx,psfrag}
\usepackage{cite}
\usepackage{sectsty}
\usepackage{epstopdf}
\usepackage{amsmath,esint, setspace, fancyhdr, amsfonts, bookmark, blindtext}
\usepackage[normalem]{ulem}
\usepackage{tikz}
\usepackage{rotating}
\usepackage[americanvoltages,fulldiodes,siunitx]{circuitikz}
\usepackage{stackengine}
\usetikzlibrary{matrix}
\usepackage{multirow}
\usepackage{multicol}
\usetikzlibrary{shapes,backgrounds,patterns}
\usetikzlibrary{mindmap,trees,decorations.markings}
\usetikzlibrary{quotes,angles}
\usepackage{verbatim}
\renewcommand{\baselinestretch}{1}
\setlength{\textheight}{8in}
\setlength{\textwidth}{6.5in}
\setlength{\headheight}{0in}
\setlength{\headsep}{0.25in}
\usepackage{graphicx}
\setlength{\topmargin}{0in}
\setlength{\oddsidemargin}{0in}
\setlength{\evensidemargin}{0in}
\setlength{\parindent}{.3in}
\usepackage{listings}
\usepackage{color} %red, green, blue, yellow, cyan, magenta, black, white
\definecolor{mygreen}{RGB}{28,172,0} % color values Red, Green, Blue
\definecolor{mylilas}{RGB}{170,55,241}
\doublespacing
\begin{document}
\begin{titlepage}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}} % Defines a new command for the horizontal lines, change thickness here
\center % Center everything on the page
%---------------------------------------------------------
% HEADING SECTIONS
%---------------------------------------------------------
\textsc{\LARGE University of Notre Dame}\\[1.5cm] % Name of your university/college
\textsc{\Large AME 40463}\\[0.5cm] % Major heading such as course name
\textsc{\large ME Senior Design}\\[0.5cm] % Minor heading such as course title
%---------------------------------------------------------
% TITLE SECTION
%---------------------------------------------------------
\HRule \\[0.6cm]
{ \huge \bfseries Design Proposal}\\[0.4cm] % Title of your document
\HRule \\[1.0cm]
%---------------------------------------------------------
% AUTHOR SECTION
%---------------------------------------------------------
%\begin{minipage}{0.4\textwidth}
\begin{center} \large
% \emph{Authors:}
\medskip
{\textsc{\textbf{Cole Mitchell} }} \quad {\textsc{\textbf{Matt Policelli}}} \quad
{\textsc{\textbf{Maurcio Segovia} }} \quad
{\textsc{\textbf{Kyra Twohy} }} % Your name, bold and small caps
\end{center}
%\end{minipage}
\begin{center} \Large
{\textsc{\textbf{Group 7} }}
\end{center}
~
%\begin{minipage}{0.4\textwidth}
%\begin{flushright} \large
%\emph{Supervisor:} \\
%Dr. James \textsc{Smith} % Supervisor's Name
%\end{flushright}
%\end{minipage}\\[4cm]
% If you don't want a supervisor, uncomment the two lines below and remove the section above
%\Large \emph{Author:}\\
%John \textsc{Smith}\\[3cm] % Your name
%---------------------------------------------------------
% DATE SECTION
%---------------------------------------------------------
\begin{center}
{\large \today}
\end{center}
% Date, change the \today to a set date if you want to be precise
%---------------------------------------------------------
% LOGO SECTION
%---------------------------------------------------------
%\vfill
\newcommand*{\plogo}{\includegraphics[width=0.25\textwidth]{nd-seal2.png}}
\plogo\\[1cm] % Include a department/university logo - this will require the graphicx package
%---------------------------------------------------------
\vfill % Fill the rest of the page with whitespace
\end{titlepage}
\newpage
\tableofcontents
\newpage
% EXECUTIVE SUMMARY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Executive Summary}
The ultimate goal of this project was to design a Thermal Energy Control System (TECS) that will heat, maintain, and cool the core temperature of a workpiece, specifically an aluminum cylinder, for a user-defined amount of time following a user-defined temperature profile. The final design, detailed in this report, is centered around a Thermoelectric Module (TEM) that is used to both heat and cool the workpiece. Additional features include a fan, used to bring the workpiece to room temperature, and a housing built to protect the user from high temperature surfaces. The total system cost $\$$565 including all materials for the physical build, as well as materials used to test during the design process. This process was informed by the constraints given by the customer or industry standards. These along with the resulting critical design decisions made to create the TECS are discussed in this report.
% DESIGN CONSTRAINTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Design Constraints}
The ultimate design was originally guided by the constraints set by the customer. The overall core temperature profile, shown in Figure 1, which is dependent on user inputted values is the fundamental constraint.
\begin{figure}[H]
\begin{center}
\includegraphics[width=5in]{tempprofile.png}
\caption{The ideal transient core temperature profile of the workpiece.}
\label{1.fig}
\end{center}
\end{figure}
The user-defined values include the maximum and final temperatures, as well as the time needed to reach and hold those given temperatures. Another guiding constraint was the size of the build volume, which would allow for three dimension printing in the future, and this was set as a 8 inch cube. The total system could not cost more than $\$$600 and could not jeopardize the user's safety in anyway. All other constraints are included in Appendix A.
% HEATING SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Heating Regimen}
The guiding principles for selecting a heating method were controllability and effectiveness. Two methods of heating were considered, a heating coil and a thermoelectric module (TEM).
The initial design proposed to heat the workpiece using a TEM, as well as a heating coil in front of a fan blowing on the workpiece. Two separate trade studies investigated capabilities of each method in meeting the customer requirements for heating the workpiece.
The method of heat transfer from the heating coil and fan combination was modeled as forced convection on a cylinder in cross-flow, where the cylinder was modeled using lumped capacitance. Correlations and equations defining these relationships are contained in Appendix B. Two extreme heating cases were considered, a rise of 10$^o$F within 3 minutes and a rise of 1$^o$F within 5 minutes both in a linear fashion, in the parametric studies to find the heat transfer coefficient to achieve the desired temperature profile. Once the heat transfer coefficient was found, the dynamics of the fluid needed from the fan could be calculated and these are shown in Figures \ref{fastCFM} and \ref{slowCFM}.
\begin{figure}[H]
\centering
\subfloat[Rise: $10^{\circ}$F in 3 min. $C$=$0.26$, $m$=$0.60$ for Eq. \ref{Nusselt} ]{\includegraphics[width=0.49\textwidth]{fast_CFM.eps}\label{fastCFM}}
\hfill
\subfloat[Rise:$1^{\circ}$F in 5 min. $C$=$0.51$, $m$=$0.50$ for Eq. \ref{Nusselt} ]{\includegraphics[width=0.49\textwidth]{slow_CFM.eps}\label{slowCFM}}
\caption{Fluid dynamics the fan needed to generate to heat the workpiece.}
\end{figure}
For the largest change in temperature over the shortest time period, the minimum flow rate for the fan was 450 cubic feet per minute, CFM. Multiple fans were found that achieved this constraint; however, the smallest change in temperature over the longest time period presented some complications. The heat transfer coefficient needed for this case resulted in flow rates under 10 CFM. Although it is theoretically possible to use pulse width modulation, PWM, to achieve these fluid dynamics, the changes between necessary values would be too fine to reliably achieve.
The heating coil was also investigated as part of this heating system. It can be modeled as a resistor, which generates heat via Joule heating. An energy balance can be done with the fluid as it enters and exits the area in which the fluid comes in contact with the coil. This would couple with the fan to create a heat transfer
model of both the air and the coil making it a complex system to model, instead the TEM was selected as a simpler system to model.
As detailed in the following section, the TEM was selected to be the cooling method as well. %spoilers
Therefore, as $\frac{|\Delta T_{max}|}{\Delta t_{min}}$ as stated in the customer requirements is the same for both heating and cooling, the criterion used to select the model of TEM was the module's ability to meet the customer cooling requirements. Equations \ref{TEMcooling} through \ref{TEMdifference} in Appendix D show that this was a valid method of selecting a TEM to be used for both heating and cooling an object, as by definition a TEM will always move more heat into the hot side than it will move out of the cold side for a power input $>0$.
% COOLING SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Cooling Regimen}
\subsection{Fan Cooling}
The method for cooling the workpiece also was selected based on the guiding principles of controllability and effectiveness. Two methods of cooling were considered, a fan and a TEM.
It was found that a standalone fan would not be able to cool the workpiece by the necessary 10$^o$F in 3 minutes, the largest decrease in temperature over the shortest amount of time, colloquially referred as the "worst case scenario" for cooling.
With \$60 allotted to the budget for cooling fans, two different fans from supplier Allied Electronics were considered, the circular Sunon A1175-HBT.TC.R.GN
and the square Ebm Papst 4414/2H \cite{supplier}.
The fans' performance in cooling the cylindrical aluminum workpiece was modeled as a cylinder in cross flow as depicted in Figure \ref{cyl}.
\begin{figure}[H]
\begin{center}
\includegraphics[width=5in]{cyl_cross_flow.png}
\caption{A diagram of a cylinder in cross flow taken from a heat transfer text book.\cite{heattransfer}}
\label{cyl}
\end{center}
\end{figure}
The Nusselt number and, therefore, the heat transfer coefficient of the cooling fans on the workpiece was taken from the average of the Hilpert, Zukauskas, and Churchill \& Bernstein correlations which correspond to Equations \ref{hilp}, \ref{zuka} and \ref{cb} respectively in Appendix C.
With forced convection correlation accounted for and a lumped capacitance model (Equation \ref{lumped}) in place, the two fans were evaluated on their ability to cool the workpiece in accordance with the worst case scenario. Particularly, the fan specifications were used to back out a heat transfer coefficient by converting a volumetric flow rate (often given in units of ft$^3$/min) into a free stream velocity in order to yield a Reynolds number.
The heat transfer coefficient of each fan as a function of volumetric flow rate (reported in CFM) can be seen as Figures \ref{h47} and \ref{circh}.
\begin{figure}[H]
\centering
\subfloat[4.7 in by 4.7 in by 1.5 in Ebm-Papst fan. ]{\includegraphics[width=0.49\textwidth]{47_h.eps}\label{h47}}
\hfill
\subfloat[6.8 in diameter, 2 in thick Sunon fan.]{\includegraphics[width=0.49\textwidth]{circ_h.eps}\label{circh}}
\caption{Heat transfer coefficient as a function cubic feet per minute for the two fans.}
\end{figure}
Next, a 10-minute snapshot of the transient temperature profile of the workpiece was parametrically determined in MATLAB for different flowrates for each fan. These curves can be seen in Figures \ref{inf47} and \ref{inf_circ}.
\begin{figure}[H]
\centering
\subfloat[4.7 in by 4.7 in by 1.5 in Ebm-Papst fan. ]{\includegraphics[width=0.49\textwidth]{to_inf_47in.eps}\label{inf47}}
\hfill
\subfloat[6.8 in diameter, 2 in thick Sunon fan.]{\includegraphics[width=0.49\textwidth]{time_to_inf_circ.eps}\label{inf_circ}}
\caption{A 10 minute snapshot of workpiece temperature demonstrating the cooling capabilities each of the fans.}
\end{figure}
% \begin{figure}[H]
% \centering
% \subfloat[4.7 in by 4.7 in by 1.5 in Ebm-Papst fan. ]{\includegraphics[width=0.49\textwidth]{to_inf_47in.eps}\label{inf47}}
% \hfill
% \subfloat[6.8 in diameter, 2 in thick Sunon fan.]{\includegraphics[width=0.49\textwidth]{{time_to_inf_circ.eps}\label{inf_circ}}
% \caption{A 10 minute snapshot of the cooling capabilities each of the fans.}
% \end{figure}
Close inspection of Figures \ref{inf47} and \ref{inf_circ} revealed that both of the fans considered were unable to decrease the temperature of the workpiece by 10 $^o$F within 3 minutes per the "worst cast scenario". The implication of this finding was that cooling the workpiece with solely a fan blowing room temperature air was not robust enough to satisfy the entire range of the design constraint.
Table \ref{tab1} below represents the figure of merit matrix for the two fans considered for the cooling regime summarizing key specifications.
\begin{table}[H]
\caption{Figure of Merit Matrix for Fans}
\centering
\label{tab1}
\begin{tabular}{ c c c } \toprule
Manufacturer & Ebm-Papst & Sunon \\ \midrule
Model & 4414/2H & A1175-HBT.TC.R.GN \\ \midrule
Geometry & Square [LxWxT] & Circular[DxT] \\ \midrule
Size (in) & 4.7x4.7x1.5 & 6.7x2 \\ \midrule
Cost ($\$$) & 40.83 & 57.18 \\ \midrule
Max air flow (CFM) & 295 & 240 \\ \midrule
Upper limit h (W/m$^2$K) & 100 & 80 \\ \bottomrule
\end{tabular}
\end{table}
The Ebm-Papst square fan was purchased to help bring the workpiece down to room temperature as support to the primary cooling system, the TEM, which is discussed in detail in the following section.
\subsection{Thermoelectric Cooling}
With a budget of $\$100$ to spend on TEMs, four TEM models were studied to determine their suitability for both heating and cooling the module in the worst-case scenario: the TE-127-1.4-1.15, TE-127-1.4-1.5, TE-127-2.0-2.5, and TE-127-1.4-2.5. However, as stated in Section 3, the heating requirements will be met by definition if the module is capable of meeting the cooling requirements.
The heat transfer into the workpiece through conduction was investigated to determine which model of TEM could best handle the worst case scenario. It was assumed that if the module could meet the heating requirements for the worst case scenario, it could also meet it for the opposite extreme case of 1$^o$F in 5 minutes by controlling the module with PWM or, if necessary, rapidly switching between heating and cooling. Because the performance of the TEM is inversely proportional to the temperature difference between the two sides, a variable duty cycle would be needed to properly control the operation of the module such that a linear temperature profile could be achieved. A model was created in MATLAB to determine the duty cycle needed to satisfied the worst case scenario based on the plots of power absorbed by the module ($Q_{module}$) against the temperature gradient of the module ($DT$) provided in manufacturer specification sheets as a function of time. A detailed discussion of how this model was constructed can be found in Appendix D. Figures \ref{fig:f1} and \ref{fig:f2} below demonstrate the worst case temperature profile being met by each of the TEMs, and the duty cycle as a function of time required to meet it, respectively.
\begin{figure}[H]
\centering
\subfloat[Temperature vs. Time]{\includegraphics[width=.5\textwidth]{Tvst.jpg}\label{fig:f1}}
\hfill
\subfloat[Duty Cycle vs. Time]{\includegraphics[width=.5\textwidth]{Duty_Cycle_vs__Time.jpg}\label{fig:f2}}
\caption{Figure 6b. shows the duty cycle required as a function of time required to match the temperature profile shown in (a). Figure 6b. shows that all four TEM models considered will be able to match the temperature profile, as the duty cycle never exceeds 1 or goes below 0.}
\end{figure}
In deciding a model, several quantities were important to consider. First, the module to must be able to meet the required cooling performance, the duty cycle $D \leq 1$ at all times. Furthermore, the lower the maximum duty cycle, $D_{max}$, was the more precise the lower the input duty cycle could be; for example, if $D_{max} = .5$, the input could be scaled such that an input of 1 resulted in a duty cycle of .5, effectively doubling the precision. Second, the temperature difference between the two sides could never be more than the burnout temperature listed in the specification sheets, which was around $70^o$C for all models that were examined. The most important characteristics considered when choosing TEMs to compare were cost and size. A figure of merit matrix, as seen on Table \ref{TEM_table}, shows different characteristics and specification for each of the TEMs studied:
\begin{table}[H]
\caption{Figure of Merit Matrix for TEM Module}
\centering
\begin{tabular}{ c c c c c } \toprule
TEM Model & TE-127-1.4-1.15 & TE-127-1.4-1.5 & TE-127-2.0-2.5 & TE-127-1.4-2.5 \\ \midrule
$D_{min}-D_{max}$ & $0.012-0.304 $ & $0.016-0.354 $ & $0.012-0.268 $ & $0.023-0.580 $\\ \midrule
Cost ($\$$) & 27.00 & 27.50 & 59.00 & 32.30 \\ \midrule
Size [L x W x T] (in) & 1.57x1.57x0.13 & 1.57x1.57x0.15 & 1.89x1.89x0.19 & 1.57x1.57x0.19 \\ \midrule
$(T_h- T_c)_{max}$ ($^o$F) & 52.75 & 53.37 & 53.29 & 52.66\\ \bottomrule
\end{tabular}
\label{TEM_table}
\end{table}
% Tables don't have captions do we still need this info: \caption{Figure of Merit matrix for the TEM models studied. The table compares the minimum and maximum fractional values of the duty cycle $D$, the cost per module (shipping and tax not included), and the dimensions of the module.}
From Table \ref{TEM_table}, it is clear that all four models studied would be capable of meeting the cooling requirements; however, the TE-127-1.4-2.5 was eliminated due to the maximum duty cycle being over 0.50. From the three remaining modules, the only noticeable difference was cost, and so the most inexpensive model was selected, the TE-127-1.4-1.15.
It is important to indicate that the figure of merit matrix gave cost per module before any tax, shipping, or miscellaneous fees. Thus, an unforeseen problem arose when procuring the selected module, it was found that the supplier would significantly up-charge orders under $\$85.00$. While the cost of purchasing two modules with the fee was still under the target budget of $\$100$, the fee was a waste of money, especially on a limited budget. As a result, the decision was made to purchase a module with similar specifications to the chosen module from a different company to avoid any upcharges. The module purchased was the Laird PC6,12,F1.4040,TA,W6 and key module specifications used to verify it as a satisfactory replacement for the TE-127-1.4-1.15 can be found in Appendix D.
% HOUSING SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Housing Design}
The housing for TECS is the third fundamental part of the heat transfer model. The housing is used to help control the environment inside the TECS and to protect the user from coming in contact with any safety hazards, including heating elements that could burn the user, which occurs at 120$^o$F \cite{MeasureMan}. The material selected for the TECS could not discolor or deform and at least one section of the the housing must allow the user to see the workpiece.
The heat transfer through the housing material was investigated to ensure that the material picked would prevent the user from being burned by the TECS. The system was modeled with radiation from the heating coil, using Equation \ref{rad}, and forced convection of warm air perpendicular to the wall using Equation \ref{convec}. Losses due to radiation and convection were neglected to ensure the material performed even in a "worst case scenario". The correlations and calculations used to obtain the outer surface temperature of the housing material are located in Appendix E. A figure of merit matrix, given in Table \ref{table3} was used to compare the different materials and inform the decision.
\begin{table}[H]
\caption{Figure of Merit Matrix for Housing Material \cite{heattransfer}}
\centering
\begin{tabular}{ c c c c } \toprule
State Variable & \shortstack{Clear UV Resistant\\ Polycarbonate\cite{poly}} & \shortstack{Yellow Pine \\ Plywood \cite{wood}} & \shortstack{Enhanced
Temperature \\UHMW \cite{umhw}}\\ \midrule
Size [LxWxT] (in) & 12x12x$\frac{1}{4}$& 96x48x$\frac{11}{32}$ & 12x12x$\frac{1}{4}$ \\ \midrule
Cost ($\$$) & 7.90 & 13.73 & 10.85 \\ \midrule
\shortstack{Thermal \\ Conductivity (W/mK)} & 0.19 & 0.13 & 0.42 \\ \midrule
\shortstack{Melting Point or \\Combustion Temperature \\($^o$F)} & 297 & 230 & 277\\ \midrule
Weight per piece (lbs) & 0.375 & 33.91 & 0.612 \\ \midrule
\shortstack{Calculated Outer \\ Surface Temperature \\($^o$F)} & 94.17 & 70.16 & 107.20 \\ \bottomrule
\end{tabular}
\label{table3}
\end{table}
Because the temperature of the surface is less than the maximum of 120$^o$F, all materials satisfy the temperature requirement. Two of the three materials considered also were below 96$^o$F, where humans feel discomfort \cite{MeasureMan}. However some performed better than others, the minimum requirement was met making the choice of the housing material based on the other material properties and the cost. The most economical material was the plywood and the size of the piece of wood was enough to create the entire housing. However, the material was not translucent, so another material was needed to create a viewing window. This was accomplished by buying one sheet of the polycarbonate to create the window. The total weight of the material was also checked to ensure the user could lift the TECS. The plywood and polycarbonate weighed just under 10 lbs, making it well under the maximum 51 lbs and the total cost was substantially under budget at a total of $\$$21.63.
% CONTROL SYSTEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Control System}
A control system was used to synchronize the performance of heating and cooling elements with the desired core temperature profile, as well as plotting the core temperature profile for the user to see. To be able to accomplish both these functions two micro controllers were used, one to control the temperature of the workpiece and one to report the temperature of the workpiece to the user. The heating and cooling of the workpiece utilizes proportional and integral (PI) control to match the core temperature. The implementation of the PI control for the governing equations of this system is detailed in the Appendix F. This type of controller was selected to minimize overshoot, without decreasing rise time. The control system takes the temperature of the object and compares it to the desired temperature at that given time on the core temperature profile. Then it outputs the power the TEM needs to put into the system to achieve the desired temperature. This can be seen in Figure \ref{logic}, shown below.
\begin{figure}[H]
\begin{center}
\includegraphics[width=6.5in]{LogicFlow.jpg}
\caption{A diagram of the logic used in the control system.}
\label{logic}
\end{center}
\end{figure}
The temperature of the workpiece was taken from an infrared temperature (IR) sensor and then fed into the control system. This was to avoid violating the constraint that contact sensors for the control system are not allowed in the build volume. The actual core temperature was found using a resistance temperature detector (RTD) and was not part of the temperature control system. This type of sensor was chosen because it functioned over the desired range of temperatures with a high level of accuracy. Discussion on the calibration of the temperature sensors is including in Appendix F. Arduinos were selected as the micro controllers because they had enough analog and digital pins for the TEM, large fan, heat sink fan, RTD, and IR sensor.
% Conclusion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
The overall product design satisfied the customer specifications and performed the task of heating, cooling, and maintaining the temperature of the workpiece in accordance with the time constraints. The total budget, found in Appendix H, is under the alloted $\$$600, but could be optimized to be less if only the material used in the final product design was purchased. The TEM and large fan integrate into a control system that follows the user profile with minimal overshoot. The controller dynamically adapts to changing conditions in room temperature or rates of desired change in temperature. The completed design can be visualized below, in Figure \ref{section} and \ref{external}.
\begin{figure}[H]
\centering
\subfloat[Internal section of the TECS.]{\includegraphics[width=0.49\textwidth]{SECTION_TECS.PNG}\label{section}}
\hfill
\subfloat[External view of the TECS.]{\includegraphics[width=0.49\textwidth]{TECS.png}\label{external}}
\caption{Computer generated models of the TECS.}
\end{figure}
All drawings of individual parts and the completed model are contained in the Appendix G. Ideally, this system will be integrated into a 3D printing system to ensure acceptable printing onto preexisting metal pieces.
% BIBLIOGRAPHY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{References}
\bibliographystyle{unsrt}
\bibliography{bibliography}
\fancyhead[R,OL]{bibliography}
% APPENDIX A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix A: Design Specifications}
A list of Design Requirements is given below.
\begin{enumerate}
\item $T_{initial}$ will be room temperature, assumed to be 68$^o$F. Actual room temperature may vary based on the given room.
\item $T_{max}$ will be targeted to be 78$^o$F, which is 10$^o$F above $T_{initial}$. $T_{max}$ will function for the entire range given in the customer requirement.
\item $T_{final}$ will be targeted to be 69$^o$ or 9$^o$ below $T_{max}$. $T_{final}$ will function for the entire range between $T_{initial}$ and $T_{max}$.
\item 5 minutes will be the targeted $t_{rise}$, but it will also function over the entire range given in the customer requirement
\item 5 minutes will be the targeted $t_{dwell_{1}}$, but it will also function over the entire range given in the customer requirement
\item 5 minutes will be the targeted $t_{fall}$, but it will also function over the entire range given in the customer requirement
\item 5 minutes will be the targeted $t_{dwell_{2}}$, but it will also function over the entire range given in the customer requirement
\item ($t_{rise}$ + $t_{dwell_{1}}$ + $t_{fall}$ + $t_{dwell_{2}}$) will be 20 minutes or less given user input
\item $T_{max}$, $T_{final}$, $t_{rise}$, $t_{dwell_1}$, $t_{fall}$, and $t_{dwell_2}$ will be entered in touchscreen monitor external to the TECS.
\item The operating voltage of the TECS cannot exceed 120V AC for safety precautions \cite{Powerstream}
\item All power sources will include fuses which limit current flow to the level(s) which corresponds to the maximum required power to prevent overheating
\item Power consumption must be minimized. Target set at 1.5 kW for the entire system\cite{kettle}.
\item The TECS should have a “build volume” that is an 8 in. x 8 in. x 8 in. open space (comparable to the build volume of the MakerSpace 3D printer).
\item The workpiece should be centrally located in the build volume and will be supported by a platform.
\item Only the workpiece and sensors can be situated in the build volume.
\item The workpiece will be aluminum based on more suitable material properties\cite{heattransfer}.
\item The workpiece must be a volume between 3 in.$^3$ and 4 in.$^3$ and can be no longer than 3 in. in any representative dimension (diameter, length, etc.)
\item The workpiece must be representative of a potentially “useful” part. Target is to use a bolt as the workpiece.
\item The TECS will display a time history of the temperature of the core on the external display to validate the temperature profile.
\item No internal sensors or monitors may be used to control the internal temperature of the workpiece. The temperature of the core will be taken using a thermocouple in the core that is not part of a control loop.
\item There is no limit (other than budgetary) for sensors on or within the TECS or non-contact sensors for the workpiece.
\item A user must be able to transport the TECS without assistance. It should be 12 lbs. to allow for repetitive lifting for a short duration \cite{USHealth}.
\item The TECS should never suffer discoloration or deformation.
\item A user of the TECS should never suffer a burn. Tissue damage is caused at 120$^o$F and discomfort is felt at 96$^o$F. Surface temperatures will not exceed 96$^o$F and insulation will be used when necessary to prevent surfaces from being exceedingly hot \cite{MeasureMan}.
\item The TECS will have a window with an unobstructed line of site to the workpiece for viewing purposes.
\item The TECS must withstand a drop to the floor from a table top.
\item The cost of the TECS cannot exceed \$600. This does not include the cost of any laptop
that may be used as a part of the TECS.
\item The TECS design will be completed by November 29th.
\item Overshoot of the exact core temperature profile will be within 5$\%$.
\end{enumerate}
% APPENDIX B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix B: Heating}
% Mauricio's equations
% Maybe put reason for TEM being able to heat workpiece if it can cool it here instead of in main text?
The aluminum workpiece was modeled using the {\it{Lumped Capacitance}} model with material properties given in \cite{heattransfer}. The Biot number, $Bi$, was found to be less than 0.10 if the heat transfer coefficient, $h_{f}$ is less than 2500$\frac{W}{m^{2}K}$. Therefore, the temperature of the workpiece was modeled by the following equation:
\begin{equation}\label{lump}
\rho{V}C_{p}\frac{dT}{dt}=q_{cond}+q_{conv}+q_{rad,net} \;
\begin{cases}
q_{cond}=-\frac{k_{s}}{L_{c}}A_{cond}(T-T_{s})\\
q_{conv}=-h_{f}A_{conv}(T-T_{f})\\
q_{rad,net}=\sigma{A_{rad}}(\alpha{G}-\epsilon{E})
\end{cases} \;
\begin{cases}
A_{cond}=\frac{\pi{D^{2}}}{4}\\
A_{conv}=\pi{D}L\\
A_{rad}=\pi{D}L+\frac{\pi{D^{2}}}{4}
\end{cases}
\end{equation}
where $\rho$ is the density of the workpiece, $V$ is the volume of the workpiece, $C_{p}$ is the specific heat capacity at constant pressure, $k_{s}$ is the thermal conductivity, $L_{c}$ is the characteristic length, $L$ is the true length, $D$ is the diameter, $T_{s}$ is the temperature of the platform, $T_{f}$ is the temperature of the fluid, $\sigma$ is Boltzman's constant, $\alpha$ is the absorptivity, $G$ is the radiation flux from the surroundings, $\epsilon$ is the emissivity, and $E$ is the radiation flux from the workpiece. The temperature of the surface of the workpiece was assumed to be at a constant room temperature. The radiation was treated such that the workpiece behaved like a grey body and the temperature of the surroundings was kept at a constant room temperature. The workpiece was modeled as a solid aluminum cylinder with with a diameter and length of 1.6 inches in cross-flow with air (properties given in \cite{heattransfer}).
Heat transfer via conduction was modeled using contact resistance between the workpiece and the platform. A value was the contact resistance, $R_{cond}$, was chosen to be 0.90$\frac{Km^{2}}{W}$ to model interaction between aluminum-aluminum surfaces \cite{heattransfer}. Therefore, conduction took the form of the following equation:
\begin{equation}\label{qcond}
q_{cond}=\frac{A_{cond}(T-T_{s})}{R_{cond}}
\end{equation}
Heat transfer via radiation was modeled using the following equation:
\begin{equation}\label{qrad}
q_{rad,net}=\epsilon\sigma{A_{rad}}(T_{surr}^{4}-T^{4})
\end{equation}
where $T_{surr}$ is the temperature of the surroundings. The radiation environment was taken to be as a large enclosure and the workpiece was assumed to be a grey body. Therefore, the emissivity and absorptivity were taken to be the same value of 0.90.
Convection was modeled using a cylinder in cross-flow, seen in Figure \ref{cyl}.
\begin{equation}\label{Nusselt}
Nu_{D}=\frac{h_{f}D}{k_{f}}=C\;Re_{D}^{m}\;Pr^{1/3}
\end{equation}
\bigskip
\noindent
Because the surrounding air completely enveloped the workpiece, free convection was not taken into account. The Nusselt number, $Nu_{D}$, given in Equation \ref{Nusselt} has certain values for the coefficient and exponent of the Reynolds number, $Re_{D}$, and is multiplied by the Prandlt number, $Pr$ \cite{heattransfer}. The values used were $C=0.26$, $m=0.60$ for a fast temperature rise and $C=0.75$, $m=0.40$ for slow temperature rise.
If a heating coil were to be used, then the dimensions would be such that it can be assumed to be lumped. Since the heating coil is modeled as a resistor, it will have heat generation via Joule heating. Thus, the equation that describes the thermal energy of the heating coil is as follows:
\begin{equation}
\rho_{c}{V_{c}}C_{c_{p}}\frac{dT_{c}}{dt}=i_{c}V_{c}-h_{f}A_{conv}(T_{c}-T_{f})
\end{equation}
where the subscript c represents the coil, $i_{c}$ is the current through the heating coil, and $V_{c}$ is the voltage drop across the heating coil. Note that the temperature of the fluid, $T_{f}$, is to be maintained constant. However, the temperature of the coil and current of the coil must be controlled as well. An energy balance can be done with the fluid as it enters and exits the area in which the fluid comes in contact with the coil. This would couple the heat transfer model of both the air and the coil making it a complex system to model. All the calculations for this were done using the code given below.
\lstset{language=matlab,frame=single}
\begin{lstlisting}
%%Trade Study: Mauricio Segovia
%Will convection with hot air be enough to increase
% the temperature of the workpiece within the
% alloted time using small fans and a heating coil
%Fan
close all; clear all; clc;
sigma=5.67*10^-8;
%%%%%%%%%% Material Properties of Aluminum %%%%%%%%%%%%
Cp=901; %J/kg-K
rho=2712; %kg/m3
ks=205; %W/m-K
e=0.90; %Emissivity of Aluminum Commercial Sheet
%Absopritivty can be 0.10-0.65 painted black 0.85-1.00
a=0.90; %Since the surface will be painted black
%%%%%%%%% Material Properties of Air%%%%%%%%%%
Pr=0.710;
kf=2.5*(10^-2); %W/m-K (thermal conductivity of FLUID)
vf=1.4*(10^-5); %m2/s
rho_f=1.2; %kg/m3
%%%%%%%%% Geometry %%%%%%%%%
%Volume 3in^3-4in^3
%Heat transfer coefficient for convection (free and forced)
%Time span of tempt change 3min<t<6min
%Max change in temperature is 10F(5.1C)
%Cylinder:
r=0.80; L=1.60; %inches
V=pi*r^2*L; %volume check in inches
r=r*0.0254; D=2*r;
L=D;
A_cond=pi*r^2;
A_conv=2*pi*r*L;
A_rad=pi*r*(r+2*L);
%%%%%%%% Ensure Lumped Capacitance%%%%%%%%%%%%
V=pi*r^2*L;
SA=2*pi*r*(r+L);
Bi=0.10; %in order to assume LCM
hf=Bi*ks*SA/V; %As long has hf<3000, we can assume LCM
%%%%%%%%Set Desired Temperature Profile%%%%%%%%%%
Tsurr=293; %(20C) the wall and platform will be at room temp
Rcond=0.90; %Contact resistance
% Choice of temperature rise and time
% deltaT=10; DELTAt=3; %change of (10F) in 3 min
% deltaT=(deltaT)*5/9;
%DELTAt=DELTAt*60; %convert to celcius and seconds
% dTdt=deltaT/DELTAt;
deltaT=1; DELTAt=5; %change of (1F) in 5 min
deltaT=(deltaT)*5/9; DELTAt=DELTAt*60;
dTdt=deltaT/DELTAt;
time=linspace(0,DELTAt,1000);
T=dTdt*time+Tsurr;
q_total=rho*V*Cp*dTdt;
q_rad_net=e*sigma*A_rad*(Tsurr^4-T.^4);
q_cond=A_cond*(T-Tsurr)/Rcond;
%%%%%%%%% Calculate the needed convection coefficient %%%%%%%%%%%%
Tf=70+[2 4 6 8]; %Temperature of the ambient air: needs
% to be higher than that of the final temperature of the
% object or else the object will never be
% heated to that temperature
Tf=convtemp(Tf,'F','K');
for i=1:length(Tf)
disp(['T_{f}= ',num2str(convtemp(Tf(i),'K','F')),'F'])
end
h=zeros(length(T),length(Tf));
figure(1)
for i=1:length(Tf)
h(:,i) = (q_total+q_cond-q_rad_net)./(A_conv*(Tf(i)-T));
plot(time,h(:,i))
hold on
end
xlabel('Time (sec)')
ylabel('Heat Transfer ^{W}/_{m^{2}K}')
legend(['T_{f}= ',num2str(convtemp(Tf(1),'K','F')),'F'],...
['T_{f}= ',num2str(convtemp(Tf(2),'K','F')),'F'],...
['T_{f}= ',num2str(convtemp(Tf(3),'K','F')),'F'],...
['T_{f}= ',num2str(convtemp(Tf(4),'K','F')),'F'])
xlim([0 DELTAt])
set(gca,'FontSize',14)
%Find the Nusselt number associated with these h-coeff.
%We will have to look for the correlation for air over a cylinder
%It is the hilbert correlation but the constants depend on the Re
%number so we are going to see how much Re changes with the
%constants
%REMEMBER FOR EXTERNAL FLOW Re_Cr=2(10^5)
ReLim=[40,1000,2*10^5,10^6];
C=[0.75 0.51 0.26 0.076]; %Re of 40,1000,10^5,10^6
m=[0.40 0.50 0.60 0.70];
%Choose an h
disp(['Choose T_{f}= ',num2str(convtemp(Tf(4),'K','F')),'F'])
hf=h(:,4);
Nu=hf*D/kf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%Get Reynolds%%%%%%%%%%%%%%%%%%%%%%%%%
Re=zeros(length(T),length(C));
figure(2)
for j=1:length(ReLim)
Re(:,j)=(Nu/(C(j)*Pr^(1/3))).^(1/m(j));
plot(time,Re(:,j))
hold on;
end
xlabel('Time (sec)')
ylabel('Reynolds Number')
legend(['Re_{lim}= ',num2str(ReLim(1))],...
['Re_{lim}= ',num2str(ReLim(2))],...
['Re_{lim}= ',num2str(ReLim(3))],...
['Re_{lim}= ',num2str(ReLim(4))])
xlim([0 DELTAt])
set(gca,'FontSize',14)
%Choose appropriate Re whose solution falls in the acceptable
%range of the coefficients for the correlations (Hilpert)
choice=2;
ReD=Re(:,choice);
disp(['Hilpert Correlation with C= ',num2str(C(choice)),...
' and m= ',num2str(m(choice))])
disp(['Limits of Correlation are Re(min)= ',num2str(0),...
' to Re(max)= ',num2str(ReLim(choice)),])
%%%%%%%%%%%%%%%Get free stream velocity%%%%%%%%%%%%%%%%
figure(3)
uf=ReD*vf/D;
plot(time,uf)
xlabel('Time (sec)')
ylabel('Free Stream Velocity (^{m}/_{s})')
legend(['Re_{lim}= ',num2str(ReLim(choice))])
xlim([0 DELTAt])
% ylim([0 30]) %cap velocity at speed of sound (300 m/s)
set(gca,'FontSize',14)
disp(['Max free stream velocity of ',num2str(max(uf)),...
' meters per sec'])
%%%%%%%%%%Get volumetric flowrate%%%%%%%%%%%%
%choose a cross-section read of flow
Aflow1=16; %inches squared
Aflow=Aflow1*(0.0254^2); %meters squared
Vflow=Aflow*uf;
disp(['Max volume flow rate of ',num2str(max(Vflow)),...
' meters cubed per sec of and area of ',num2str(Aflow1),'in2'])
%Now find the volumetric flowrate in cubic feet per minut3
%because this is what fans will use
CFM=Vflow*2118.88; %cubic feet per minute
figure(4)
plot(time,CFM)
xlabel('Time (sec)')
ylabel('Volume Flow Rate (CFM)')
legend(['Re_{lim}= ',num2str(ReLim(choice))])
xlim([0 DELTAt])
% ylim([0 300]) %cap velocity at speed of sound
set(gca,'FontSize',14)
% figure(4)
% mFlow=rho_f*Vflow;
% plot(time,mFlow)
% xlabel('Time (sec)')
% ylabel('Mass Flow (^{kg}/_{s})')
% legend(['Re_{lim}= ',num2str(ReLim(4))])
% xlim([0 DELTAt])
% % ylim([0 300]) %cap velocity at speed of sound
% set(gca,'FontSize',14)
\end{lstlisting}
% APPENDIX C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix C: Cooling}
The correlations and formulas used to model the fan cooling system are given below. The workpiece temperature was modeled by Equation \ref{lump} which makes the lumped capacitance assumption for a sufficiently small Biot number ($Bi \leq$ 0.1) \cite{heattransfer}.
\begin{equation}
\rho C_p V\frac{dT}{dt} = h(T_s - T_{\infty}),
\label{lumped}
\end{equation}
The Hilpert correlation (Equation \ref{hilp}), the Zukauskas (Equation \ref{zuka}), and Churchill \& Bernstein correlation (Equation \ref{cb}) are all used to estimate the surface measured Nusselt number of a cylinder in crossflow \cite{heattransfer}.
\begin{align}
\overline{Nu} &= C Re^m Pr^{1/3}, \label{hilp} \\
\overline{Nu} &= C Re^m Pr^{0.37} \left(\frac{Pr}{Pr_s}\right)^{1/4}, \label{zuka} \\
\overline{Nu} &= 0.3 + \frac{0.62 Re^{1/2}Pr^{1/3}}{[1 + (0.4 \,Pr)^{2/3}]^{1/4}}\left[1 + \frac{Re}{282,000}^{5/8}\right]^{4/5}, \label{cb}
\end{align}
The following MATLAB code used to parametrically study the effect of different fan speed and their cooling capabilities on the workpiece based on example problem from the heat transfer text.
\newpage
\begin{lstlisting}
%% Cooling Fan (Forced Convection) Trade Study
% Author: Cole Mitchell (901847474)
clc; clear;
%% Initial set up
Afan = 0.12^2; % area of fan [m^2]
CFM = 295/2; % enter this from spec sheet [ft^3/m]
vdot = 0.0004719474.*CFM; % conversion from CFM to m^3/s
u = vdot/Afan; % free stream velocity [m/s]
rho = 1.2; % kg/m^3
mdot = rho.*vdot; % mass flow rate
Prs = 0.71435; % Prandtl #
Pr = 0.71538; % Prandtl # of surface
R = 0.041/2; % radius of the WP
D = 0.041; % diameter of wP
L = 0.041; % length of WP
Lc = pi*R^2*L/(2*pi*R*(R+L)); % critical length calc for cyl
%% Zukauskas Correlation:
Re = Lc*u/(1.5111E-5) % Reynolds #
% Correlation constants
c = 0.193;
n = 0.37;
m = 0.618;
Nu_z = c.*Re.^m*Pr^n*(Pr/Prs)^.25; % Nusselt #
k = 0.025679; % thermal conductivity of air
h_z = Nu_z.*k./Lc;
%% Churchill & Bernstein Correlation
Nu_cb = 0.3 + ((0.62.*Re.^.5.*Pr^(1/3))./(1+
(0.4*Pr)^(2/3))^.25).*(1 + (Re./282000).^(5/8)).^.8;
h_cb = Nu_cb*k/Lc;
%% Hilpert Correlation
Nu_h = 0.683.*Re.^.466*Pr^(1/3);
h_h = Nu_h*k/Lc;
%% Average and plot h
h = (h_cb + h_z + h_h)./3; % average all three HT coeffs
figure(4)
f1 = plot(CFM, h_z, CFM, h_cb, CFM, h_h, CFM, h, 'b--')
set(f1, 'linewidth', 1.5)
xlabel('Volumetric Flow Rate [ft^3/min]' )
ylabel('Heat Transfer Coefficient [W/m^2-k]')
title('120 mm x 120 mm Fan')
legend('Zukauskas', 'Churchill & Bernstein', 'Hilpert',...
'Average', 'Location', 'Southeast')
set(gca, 'FontSize', 14)
%% Force Analysis & Plot
P = 1.2/2*u.^2; % Bernoulli's Principle
Area = pi*R*L; % m^2
F_bern = P.*Area;
mu = 1.05; % coef. of friction Al on Al
w = 1.39; % N
figure(2)
f2 = plot([100 300], [w*mu w*mu], 'k--', CFM, F_bern);
set(f2(1), 'linewidth', 1.5);
set(f2(2), 'linewidth', 1.5);
hold on
loglog(CFM, F_bern)
xlabel('Volumetric Flow Rate [ft^3/min]' )
ylabel('Force [N]')
legend('Friction force','Drag force on WP')
set(gca, 'FontSize', 14)
clear;
R = 0.041/2; % radius of the WP
D = 0.041; % diammeter of the WP
L = 0.041; % length of the WP
rho = 2.7e3; % density of Al
Cp = 900; % specific heat of Al
V = pi*R^2*L; % volume of the WP
h = 90; % pick off a heat transfer coeff from
A = pi*D*L; % cross-sectional area
Ts(1) = 298.71; % starting temp
Tinf(1) = 293.15; % ambient temp
i = 1; % iterator
dt = 10; % time step
time(1) = 0
while (time < 60*10)
dT(i) = (h.*A*(Ts(i)-Tinf)/(rho*Cp*V))*dt; % change in T
i = i +1;
time(i) = time(i-1) + dt;
Ts(i)= Ts(i-1) - dT(i-1);
if (dT(i-1) < 0) % stops if dT goes below 0
break
end
end
Tsc1 = Ts -273.15; % convert to degrees C
Tsf1 = Tsc1*(9/5) + 32 % convert to degrees F
%% Plot
plot(time/60,Tsf1, '-')
hold on
set(gca, 'FontSize', 14)
xlabel('Time (min)')
ylabel('Temperature of workpiece (^oF)')
set('linewidth', 1.5)
\end{lstlisting}
\newpage
% APPENDIX D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Appendix D: TEM Cooling}
\subsection{Operation of TEMs}
Thermoelectric modules utilize the Peltier effect to create a temperature gradient across semiconductors of two different materials when a current is passed through the semiconductors. The direction of heat flow depends on the polarity of the current, and is fully reversible. The heat moved through the cold side of the TEM, $Q_c$ is given by the equation
\begin{equation}\label{TEMcooling}
Q_{c} = \frac{\beta V_{in}}{R_{TEM}} = \beta I_{in},
\end{equation}
where $\beta$ is the Peltier coefficient of the module, and depends on the temperature difference between the two sides as well as the types of semiconductors used. The temperature dependence of the Peltier coefficient is seen in manufacturer specification sheets, from which the following relationship can be extrapolated for a given voltage (or current depending on the manufacturer) applied to the module:
\begin{equation}\label{TEMresponse}
Q_{C}(T_{hot}-T_{cold}) = \alpha (T_{hot} - T_{cold}) + Q_{C|T_{hot} = T_{cold}}
\end{equation}
This gives the heat which is passed through as a function of the temperature difference between the sides, which is important for controlling the heat being removed or added to the workpiece.
Further, the heat transferred through the hot side of the TEM is given by the equation
\begin{equation}\label{TEMheating}
Q_{h} = Q_{c} + P_{in} = Q_c + V^2_{in}/R_{TEM},
\end{equation}
where $P_{in}$ represents the joule heating of the module, as would be seen in a heating coil. Because for $Q_c = 0 if V_{in} = 0$, for any $Q_{c} > 0$,
\begin{equation}\label{TEMdifference}
Q_h > Q_c
\end{equation}
\subsection{Construction of Mathematical Model And Selected TEM Specifications}
The TEM module was selected based on a mathematical model that computed $T_{upper}$, $T_{lower}$, radiative and convective heat loss out of the workpiece, $Q_{top}$, $Q_{bottom}$, and duty cycle all as a function of time. In doing this, it was key to assume lumped capacitance, which was a valid assumption because the thermal mass of the TEM plates are extremely low and because the workpiece itself was designed such that lumped capacitance could be assumed. Because of this assumption, it could be assumed that $T_{upper}(t) = T_{required}(t)$.
The temperature of the bottom side of the EM can be calculated using the equation
\begin{equation}\label{tlower}
T_{lower} = T_{\infty} + \Theta Q_{h},
\end{equation}
where $\Theta$ is the thermal resistance of the heat sink\footnote{The heat sink chosen (and purchased) was CTS Electronic Components Heat Sink Model APF40-40-13CB, which has a thermal resistance of $1.9\frac{^oC}{W}$}. A fan is used to maintain the thermal resistance of the heat sink, and due to this $T_{\infty}$ is assumed to remain constant, and is not incremented as the upper temperature is. Once the lower temperature is calculated, the temperature difference between the upper and lower sides can be computed:
\begin{equation}
DT(t) = |T_{upper}-T_{lower}|
\end{equation}
This allows the performance of the TEM to be computed:
\begin{equation}
Q_{h}(t) = \alpha_h DT + Q_{h|T_{hot}-T{cold} = 0}
Q_{c}(t) = \alpha_c DT + Q_{c|T_{hot}-T{cold} = 0}
\end{equation}
This allows the duty cycle to be computed for both heating and cooling regimes:
\begin{equation}
D = \frac{Q_{needed}}{Q_{TEM}},
\end{equation}
where $Q_{TEM} = Q_{h}$ for heating and $Q_{c}$ for cooling, with losses accounted for in $Q_{needed}$.
From this, the temperature change of the workpiece can be calculated over time step $\delta t$:
\begin{equation}
\delta T = \frac{D*Q_{TEM}*\delta t}{\rho C_p V}
\end{equation}
The temperature at time $t+\delta t$ can then be computed:
\begin{equation}
T(t+\delta t) = T(t) + \delta T
\end{equation}
Below is the MATLAB code constructed from these equations which was used to test the suitability of different TEM models for use as a heat source and cooling device:
\begin{lstlisting}
% Trade Study: TEM Module Choice
% Matthew Policelli
% 14 September 2017
% Startup Stuff
clear
close all
clc
global k_al Cp_al rho_al D_wp H_wp A_wp As_wp V_wp
global W_TEM L_TEM A_TEM Aexp_TEM
global k_air Pr v
test = input('Type "1" to test, Type "2" to input values\n');
% Workpiece Properties
%material properties
k_al = 205; %Thermal conductivity, W/mK
Cp_al = 0.91e3; %J/kgK
rho_al = 2800; %kg/m^3
%dimensions
D_wp = 38.1e-3; %m (Diameter of the workpiece = 1.5 in)
H_wp = 44.45e-3; %m (Height of the workpiece = 1.75 in)
%Areas and Volumes
A_wp = pi*power(D_wp,2)/4;
%Area of top or bottom surface of the workpiece, m^2
As_wp = A_wp + pi*D_wp*H_wp;
%Surface area of workpiece exposed to air, m^2
V_wp = A_wp*H_wp; %Volume of the workpiece, m^3
% Air Properties
Pr = 0.707;
k_air = 2.57e-2; %W/mK
v = 15.11e-6; %kinematic viscosity, m^2/s
% TEM Properties
% Need equations for DT response of current and
% Qc for lowest voltage of each TEM
% input dimensions of TEM
if test == 1
W_TEM = 40e-3;
L_TEM = 40e-3;
end
if test == 2
W_TEM = 1e-3*input('Width of TEM (mm)\n'); %Width of TEM, m
L_TEM = 1e-3*input('Length of TEM (mm)\n'); %Length of TEM, m
end
% Calculate Surface areas for TEM
A_TEM = W_TEM*L_TEM; %surface area of TEM, m^2
Aexp_TEM = A_TEM - A_wp; %Area of TEM exposed to air, m^2
if test == 1
a1 = 1;
b1 = 23;
a2 = 1;
b2 = 27;
end
if test == 2
% Input DT response equations
% -Qc = -a1*DT+b1
a1 = input('Slope of TEM Qc line\n');
b1 = input('Intercept of TEM Qc line\n');
% Qh = -a2*DT+b2
a2 = input('Slope of TEM Qh line\n');
b2 = input('Intercept of TEM Qh line\n');
end
% Temperature and Time inputs
dT_const = 0;
dt_const = 600; %time held constant, seconds (10 mins)
if test == 1
dT_h = 5.5556;
dt_h = 180;
dT_c = dT_h;
dt_c = 180;
end
if test == 2
dT_h = input('Heating Temperature Change, deg C\n');
dt_h = 60*input('Time to heat, minutes\n'); %seconds
dT_c = input('Cooling Temperature Change, deg C\n');
dt_c = 60*input('Time to cool, minutes\n'); %seconds
end
% From Temp and Time inputs, calculate Q_needed
Qh_needed = rho_al*Cp_al*V_wp*dT_h/dt_h;
Qconst_needed = 0;
% Calculate Duty Cycle as a function of time
dt = .005; % 10 Hz frequency
T = 20; %room temperature, Celcius
i = 1;
n = 0;
% Heating
Tl = 20;
Tlstore(1) = Tl;
DT = 0;
for t = 0:dt:dt_h
% Store Values
Tstore(i) = T;
tstore(i) = t;
% Compute Tl
if t == 0
Q_exh = 0;
end
if t>0
Q_exh = a1*DT-b1;
SinkResist = 1.9;
Tl = 20 + SinkResist*(D*Q_exh);
Tlstore(i) = Tl;
end
% Compute DT
DT = T-Tl;
DTstore(i) = DT;
% Compute Q_TEM
Q_TEM = -a2*DT+b2;
% Compute Q_losses
Tfilm = (T+20)/2;
beta = 1/Tfilm;
Gr = 9.81*beta*(Tfilm-20)*power(D_wp,3)/power(v,2);
Ra = Gr*Pr;
% Convention: Le Fevre and Ede
NuL = (4/3)*power((7*Ra*Pr)/(5*(20+21*Pr)),1/4)+
(4*(272+315*Pr)*H_wp)/(35*(64+63*Pr)*D_wp);
h = NuL*k_air/H_wp;
Q_conv = h*As_wp*(T-20);
% Radiative Losses
Q_rad = 5.67e-8*(power(T+273,4)-power(293,4))*As_wp;
% Total losses = Sum of Radiative and Convective Losses
Qlosses = Q_conv+Q_rad;
%Qlosses_store(i) = Qlosses;
% Compute D
if n == 10
n = 0;
end
if n == 0
D = (Qh_needed+Qlosses)/Q_TEM;
end
Dstore(i) = D;
Qstore(i) = D*Q_TEM;
Pstore(i) = D*(Q_TEM-Q_exh);
% compute temperature change over sensor period
dT = 1/(rho_al*Cp_al*V_wp)*(D*Q_TEM-Qlosses)*dt;
% increment
T = T+dT;
i = i+1;
n = n+1;
end
% Constant Part
t1 = t;
for t = t1:dt:t1+dt_const
% Store Values
Tstore(i) = T;
tstore(i) = t;
% Compute Tl
Q_exh = a1*DT-b1;
SinkResist = 1.9;
Tl = 20+ SinkResist*(D*Q_exh);
Tlstore(i) = Tl;
% Compute DT
DT = T - Tl;
DTstore(i) = DT;
% Compute Q_TEM
Q_TEM = -a2*DT+b2;
% Compute Q_losses
Tfilm = (T+20)/2;
beta = 1/Tfilm;
Gr = 9.81*beta*(Tfilm-20)*power(D_wp,3)/power(v,2);
Ra = Gr*Pr;
% Convention: Le Fevre and Ede
NuL = (4/3)*power((7*Ra*Pr)/(5*(20+21*Pr)),1/4)
+(4*(272+315*Pr)*H_wp)/(35*(64+63*Pr)*D_wp);
h = NuL*k_air/H_wp;
Q_conv = h*As_wp*(T-20);
% Radiative Losses
Q_rad = 5.67e-8*(power(T+273,4)-power(293,4))*As_wp;
% Total losses = Sum of Radiative and Convective Losses
Qlosses = Q_conv+Q_rad;
% Compute D
if n == 10
n = 0;
end
if n == 0
D = Qlosses/Q_TEM;
end
Dstore(i) = D;
Qstore(i) = D*Q_TEM;
Pstore(i) = D*(Q_TEM - Q_exh);
% compute temperature change over sensor period
dT = 1/(rho_al*Cp_al*V_wp)*(D*Q_TEM-Qlosses)*dt;
%increment
T = T+dT;
i = i+1;
n = n+1;
end
t2 = t;
% Cooling Part
Qc_needed = rho_al*Cp_al*V_wp*-dT_c/dt_c;
for t = t2:dt:960
% Store Values
Tstore(i) = (T);
tstore(i) = t;
% Compute Tl
Q_exh = -a2*DT+b2;
SinkResist = 1.9;
Tl = 20 + SinkResist*(D*Q_exh);
Tlstore(i) = (Tl);
% Compute DT
DT = Tl-T;
if DT < 0
DT = 0;
end
DTstore(i) = DT;
% Compute Q_TEM
Q_TEM = a1*DT-b1;
% Compute Q_losses
Tfilm = (T+20)/2;
beta = 1/Tfilm;
Gr = 9.81*beta*(Tfilm-20)*power(D_wp,3)/power(v,2);
Ra = Gr*Pr;
% Convention: Le Fevre and Ede
NuL = (4/3)*power((7*Ra*Pr)/(5*(20+21*Pr)),1/4)
+(4*(272+315*Pr)*H_wp)/(35*(64+63*Pr)*D_wp);
h = NuL*k_air/H_wp;
Q_conv = h*As_wp*(T-20);
% Radiative Losses
Q_rad = 5.67e-8*(power(T+273,4)-power(293,4))*As_wp;
% Total losses = Sum of Radiative and Convective Losses
Qlosses = Q_conv+Q_rad;
% Compute D
if n == 10
n = 0;
end
if n == 0
D = (Qc_needed+Qlosses)/Q_TEM;
end
Dstore(i) = D;
Qstore(i) = (D*Q_TEM);
Pstore(i) = (D*(Q_exh-Q_TEM));
% compute temperature change over sensor
dT = 1/(rho_al*Cp_al*V_wp)*(D*Q_TEM-Qlosses)*dt;
%increment
T = T+dT;
i = i+1;
n = n+1;
end
\end{lstlisting}
\subsection{TEM Specifications for Chosen Module}
%Link: https://assets.lairdtech.com/home/brandworld/files/THR-DS-PC612F14040TAW6_092215.pdf
Important specifications considered when choosing the module after deciding to go with an alternative were $DT_{max}, I_{max}, V_{max}$, and the Q vs. DT plot. These are listed below for a hot side temperature of $25^o C$, with Figure \ref{TEMspec} showing the Q vs. DT plot \cite{Laird}.
\begin{center}
$DT_{max} = 67^o C$
\newline $I_{max} = 6A$
\newline $V_{max} = 14.5A$ \newline
\end{center}
\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{Capture.PNG}
\caption{Q vs. DT plot for different current imputs}
\label{TEMspec}
\end{center}
\end{figure}
% APPENDIX E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix E: Housing}
Correlations and equations used to model the heat transfer through the TECS housing material are given below. The radiation of the heating coil, shown in Equation \ref{rad} \cite{heattransfer}, and increased air temperature inside the housing were considered to find the inner surface temperature, $T_{s1}$.
\begin{eqnarray} \label{rad}
{q_{rad}}&=&{A\sigma\epsilon({T_{coil}}^4-{T_{s1}}^4)},
\end{eqnarray}
The convection was modeled using the correlation for a thin vertical plate perpendicular to the flow given in Equation \ref{convec} \cite{heattransfer}. It should be noted that for the fans being considered, the Reynold's number would not be high enough to achieve forced convection for the entire length of the plate. Instead, the Reynolds number was assumed to be 10,000, which is the lowest value needed to have forced convection. Although this is an inaccurate model of the actual system, it approximates the temperature to be higher, allowing for the actual system to meet all constraints by wider margins in practice.
\begin{eqnarray} \label{convec}
{\frac{\bar{h}D}{k}}&=&{C{Re_{D}}^{m}Pr^{1/3}}
\end{eqnarray}
The outer surface temperature was then found using the one dimensional conduction equation, shown in Equation \ref{conduc}, below \cite{heattransfer}.
\begin{eqnarray} \label{conduc}
{q_{cond}}&=&{\frac{kA}{L}(T_{s1}-T_{s2})}
\end{eqnarray}
Below is the MATLAB code used to approximate the outer temperature of the housing and the example problem from the heat transfer book.
\begin{lstlisting}
% Trade Study: TECS Housing Material
% Kyra Twohy
% 14 Spetember 2017
%% Set Material Properties
% get meterial properties for plastics, wood, metal, and air
% Clear UV Resistant Polycarbonate
% (https://www.mcmaster.com/#standard-plastic-sheets/=19dukul)
cp_p=1100; % J/kg-C from Engineering Toolbox
k_p=0.19; % W/m-K from Engineering Toolbox
t_p=0.003175*2; % m, +/- 0.006"
tensile_p=8000; % psi
epsilon_p=0.9;
% Enhanced Temperature UHMW
k_u=0.42;
t_u=0.003175*2;
% Marine-Grade Plywood Sheets
% (https://www.mcmaster.com/#wood/=19dv51h)
t_w=0.008731255; % m
k_w=0.13; % W/m-K from Engineering Toolbox
cp_w=2.5; % kJ/kg-K
epsilon_w=0.885;
% Properties of Air
Pr=0.710;
kf=2.5*(10^-2); % W/m-K (thermal conductivity of FLUID)
v=1.4*(10^-5); % m^2/s
rho_f=1.2; % kg/m^3
A=0.3556^2; % m
C=0.667; % from table 7.3 in HT book vertical thin wall
m=0.5; % from table 7.3 in HT book
D=0.3556; % m
V=0.05522; % m^3/s
nu=15.89e-6; % m^2/s
Re=10000; %V*D/nu % Reynold's number for air,
this isn't high enough for this correlation, assume 10000
h=(C*Re^(m)*Pr^(1/3)*kf)/D;
% set for type of convection, vertical thin wall
Tair=40+273.15; % K
% Properties of Resistance Heater
Pwr=100; % W
alpha=0.82;
Tcoil=102+273.15; % K
Ac=0.03; % m^2
% Properties of Radiation
sigma=5.67e-8; % W/m^2-K
%% Heat Transfer Analysis
% radiation equation
% outer temperature can't be more than 96F
% Find radiation from the Heating Coil
syms Ts1
eq1=-sigma*alpha*A*(Tcoil^4-Ts1^4)-h*A*(Tair-Ts1)+Pwr*0.5;
% find the inside temperature of the wall
s=solve(eq1,Ts1);
Ts11=s(2,1);
T_s1=vpa(Ts11,4);
T_s2_p=double((((-0.5*Pwr)*t_p)/(k_p*A))+T_s1);
% estimate outer surface of the wall, polycarbonate
T_s2_u=double((((-0.5*Pwr)*t_u)/(k_u*A))+T_s1);
% estimate outer surface of the wall, uhmw
T_s2_w=double((((-0.5*Pwr)*t_w)/(k_w*A))+T_s1);
% estimate outer surface of the wall, plywood
T_s2_pF=convtemp(T_s2_p,'K','F');
T_s2_uF=convtemp(T_s2_u,'K','F');
T_s2_wF=convtemp(T_s2_w,'K','F');
%% Example Problem 3.1, using conduction equation
Rtot=(35-10)/(100);
Lins_a=0.014*(1.8*0.25-((3e-3)/(0.3))-(1/(2+5.9)));
Lins_w=0.014*(1.8*0.25-((3e-3)/(0.3))-(1/(200)));
q=35-100*((3e-3)/(0.3*1.8)); % correct answer is 34.4C -> verified
%% Example Problem 7.25, using convection correlation
% Solving for the one h using the same correlation
h_test=(0.0263/1)*0.667*((5*(1/15.89e-6))^.5)*(0.707)^(1/3);
% correct value is 8.7 W/(m^2 K)-> verified
%% Example Problem 12.12, using radiation equation
qflux=(0.95*750)-0.22*(120-30)^(4/3)-(0.1*5.67e-8)*(393^4-263^4);
% correct answer is 516 W/m^2,
% *note no area in this equation it is just flux
\end{lstlisting}
% APPENDIX F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix F: Control System}
The following governing equation was used to describe the transient temperature profile of the workpiece:
\begin{equation}\label{Control1}
\rho{V}C_{p}\frac{dT}{dt}=Q_{TEM}-h_{forced}A_{conv}(T-T_{\infty})+q_{rad,net}+q_{free}
\end{equation}
From the preliminary heat transfer calculations, it was found that the losses due to free convection and net radiation range from $0.10-0.30$W. For the purpose of deriving a control law, those modes of heat transfer are assumed constant and are ignored. Therefore, Equation \ref{Control1} reduces to the following expression:
\begin{equation}\label{Control2}
\rho{V}C_{p}\frac{dT}{dt}=Q_{TEM}-h_{forced}A_{conv}(T-T_{\infty})
\end{equation}
To further simplify the problem, only one input will be considered: the power of the TEM. The heat transfer coefficient for forced convection will be made constant. Additionally, a non-dimension variables is introduced in order to simplify the equation when taking the Laplace transform:
\begin{equation}\label{nondim}
\Theta=\frac{T-T_{initial}}{T_{max}-T_{initial}} \qquad\hbox{therefore, $\Theta(t=0)=0$.}
\end{equation}
Introduction Equation \ref{nondim} into Equation \ref{Control2}, defining $\Delta{T_{max}}=T_{max}-T_{initial}$, and noting that $T_{initial}=T_{\infty}=T_{room}$ gives the following expression:
\begin{equation}\label{Control3}
\rho{V}C_{p}\Delta{T_{max}}\frac{d\Theta}{dt}=Q_{TEM}-h_{forced}A_{conv}(\Delta{T_{max}}\Theta)
\end{equation}
The power of the TEM is made proportional to the error in the nondimensionalized temperature and made proportional to the integral of the error. Therefore,
\begin{equation}\label{Control4}
Q_{TEM}=k_{p}(\Theta_{des}-\Theta)+k_{i}\int_{0}^{\hat{t}}(\Theta_{des}-\Theta)d\hat{t}
\end{equation}
Combing Equations \ref{Control3} and \ref{Control4}, applying the Laplace transform and solving for the ratio of $\Theta$ to $\Theta_{des}$ yields
\begin{equation}\label{transferfunction}
\frac{\Theta}{\Theta_{des}}=\frac{sk_{p}+k_{i}}{ \rho{V}C_{p}\Delta{T_{max}}s^{2}+(k_{p}+h_{forced}A_{conv}\Delta{T_{max}})s+k_{i}}
\end{equation}
Classical control theory states that the zeros and poles of the transfer function (as in Eq. \ref{transferfunction}), dictate the response given an input. To have a stable response, the poles of the transfer function \ref{transferfunction} need to be in the left region of the complex plane. In other words, by letting $a=\rho{V}C_{p}\Delta{T_{max}}$, $b=h_{forced}A_{conv}\Delta{T_{max}}$
\begin{equation}\label{conditions}
s=\frac{-(b+k_{p})\pm\sqrt{(b+k_{p})^{2}-4ak_{i}}}{2a}\quad\hbox{such that}\quad\begin{cases}
(b+k_{p})^{2}-4ak_{i}<0 \qquad\hbox{or}\\
-(b+k_{p})\pm\sqrt{(b+k_{p})^{2}-4ak_{i}}<0
\end{cases}
\end{equation}
In order to see for what values of $k_{p}$ and $k_{i}$ the conditions in \ref{conditions} are satisfied, a Mathematica notebook was written. The following conditions were found:
\begin{equation}\label{kilimit}
k_{i_{limit}}=\frac{a^{2}+2ak_{p}+k^{2}_{p}}{4b} \qquad\begin{cases}
k_{i}\leq{k_{i_{limit}}}\implies\hbox{poles are complex} \\
k_{i}>k_{i_{limit}}\implies\hbox{poles are negative and real}
\end{cases}
\end{equation}
Note that poles closest to the imaginary line will more or less dominate the system response. Thus, the condition that made the poles complex was chosen. A parametric study was done in MATLAB to see the possible pairs of $k_{p}$ and $k_{i}$ that make the poles complex and satisfy the necessary condition give different values for the forced heat transfer coefficient.
The calibration of the temperature sensors was done to ensure accurate integration with the control system. For the RTD, the manufacturer data sheet was used to to get the resistance of the sensor as a function of temperature. For the range of temperatures the workpiece would reach during the testing phase, the relationship between resistance and temperature was extrapolated using the curve fitting tool in Matlab. The third order function was used in the code to accurately report the core temperature.
%% APPENDIX G %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Appendix G: Design Drawings}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{workpiece1.pdf}
\label{workpiece}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{heatsink.pdf}
\label{heatsink}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{SMALL_FAN.png}
\label{smallfan}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{TEM_MOUNT.jpg}
\label{temmount}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{LARGE_FAN.png}
\label{largefan}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{windowplate.pdf}
\label{window}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{frontdoor.pdf}
\label{frontdoor}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{LEFT_PANEL.png}
\label{left}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{RIGHT_PANEL.png}
\label{right}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{TOP_PANEL.png}
\label{top}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{BOTTOM_PANEL.png}
\label{bottom}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{BACK_PANEL.png}
\label{bottom}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{IRsensor.pdf}
\label{irsensor}
\end{center}
\end{figure}
\begin{figure}[H]
\begin{center}
\includegraphics[width=7.5in,angle=90]{TECS.pdf}
\label{TECS}
\end{center}
\end{figure}
% APPENDIX H: BUDGET %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Appendix H: Complete Budget}
\begin{table}[H]
\centering
\caption{Cost of Individual Components for TECS}
\label{my-label}
\begin{tabular}{ccc}
\textbf{Item Description} & \textbf{Total Price [USD]} \\
Hermetically Sealed Flexible Thermistor Sensors, PFA Jacketed (K21) & 64.00 \\
1.75in x 6in Aluminum 6061 Rod & 18.57 \\
Amazon: Electrical Tape & 8.99 \\
Amazon: Resistor Set & 10.99 \\
Amazon: RedBoard & 19.95 \\
Amazon: 5PCS Thermistor Anycubic NTC 3950 100K & 8.99 \\
Mouser: Thermopile Sensor & 0.00 \\
Amazon: 5pc N Mosfet & 10.98 \\
Amazon: 3pcs Solderless Breadboard & 7.99 \\
Amazon: 60V 5A Rectifier Diodes & 9.99 \\
McMasterCarr: Clear UV-Resistant Polycarbonate 12x12x0.125 & 7.90 \\
DigiKey: Thermoelectic Module & 85.73 \\
DigiKey: Heat Sink & 8.09 \\
DigiKey: Small Fan & 18.35 \\
Adafruit: Thermopile 5V Sensor & 41.85 \\
Home Depot: BC Sanded Pine Plywood & 21.75 \\
Amazon: AC-DC 12V 5A & 9.98 \\
Amazon: Prime & 11.00 \\
McMasterCarr: Surface Mount Hinge & 4.50 \\
McMasterCarr: Plastic Pull Handle & 5.82 \\
McMasterCarr: 10-24 0.50 Length Head Screw & 6.31 \\
McMasterCarr: Steel Corner Bracket & 20.00 \\
AlliedDelec: 24V Fan 4414/2H & 54.80 \\
McMasterCarr: 10-24 1.00 Length Head Screw for Wood & 10.52 \\
Brownell (mini Cord and INA122PA Amp) & 35.00 \\
McMasterCarr: Steel Hex Nut 10-24 & 1.84 \\
McMaster Screw 10-24 1.0 Length Phillips & 7.28 \\
Amazon: 180kOhm Resistors & 9.77 \\
Kabelin Ace Hardware Store & 3.00 \\
Amazon: Protoboard Prototyping PCB Diy Soldering Universal & 10.00 \\
Stores & 15.00 \\
Amazon: ACDC 12V 5A & 16.00 \\
\textbf{Total} &
\textbf{564.94}
\end{tabular}
\end{table}
\end{document}