% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{url}
\usepackage{float}
\usepackage{romannum}
\usepackage{textcomp} 
\usepackage{hyperref}

%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION 
\else STUDENT VERSION \fi\\Temperature Distribution in a Uniform Slender Bar}





\markboth{Temperature Distribution in a Uniform Slender Bar}{Temperature Distribution in a Uniform Slender Bar}

\author{Yuxin Zhang\\
Washington State University\\
Richland WA 99354 USA\\
yx.zhang@wsu.edu
}

\makeStitlePDFLaTex


\begin{abstract}
The temperature distribution along a uniform slender bar due to conduction and convection is investigated through experimental, analytical, and numerical approaches. A series of experiments are designed to study the effects of materials, ambient fluid flows, geometric characteristics, and boundary conditions on the temperatures at prescribed locations of the bar. A steady-state temperature distribution model is developed based on the principle of conservation of energy. The resulting second-order ordinary differential equation is solved analytically with specified boundary conditions. A numerical scheme using the finite volume method on a uniform grid is developed. Numerical simulations are performed in accordance with various physical scenarios in the experiments. The comparison of analytical, numerical solutions, and the experimental data provides insights on the heat transfer process and show correlation between modeling and experiments.

\end{abstract}


\section*{SCENARIO DESCRIPTION}

Conduction and convection are two important modes of heat transfer. To understand the mechanisms of these two types of heat transfer processes, we study the temperature distribution in uniform slender bars though experiments, modeling, and numerical simulations. Various boundary conditions are investigated to demonstrate the effects of conduction and convection on the temperature distribution. Some of the related engineering applications of this study include the design of a heat exchanger with fins or a thermal radiator.



\section{Experiments}

We investigate the steady-state temperature distribution within a uniform slender bar, such as a thin rod or a slender rectangular bar. 

\subsection{Preparation}

We use the following experimental equipment and materials: $1$m-long, thin metal bars (circular and rectangular), thermocouple, infrared thermometer, electric heater, water bath, insulation layer, and tape.
To account for the effect of thermal conductivities \cite{Bergman2011} on conductive heat transfer, we use bars made of three different materials, i.e. aluminum, steel and copper. In addition, two ambient conditions - air and water - will be considered in order to investigate the effect of convective heat transfer on the temperature distribution. 

In the following experiments, we will set and maintain the left end of the metal bar at a constant high temperature $T_b$, for example, $T_b=80^{\circ}\mathrm{C}$, while the lateral surface and the right tip of the bar are subject to varying conditions. We will measure the temperature distribution along the bar at different locations $x=0.2, 0.4, 0.6, 0.8 ~\mbox{and}~ 1$m.

\subsection{Conduction Only}\label{sec12}

To understand heat transfer due to conduction, we will start with a simple scenario. 

(\romannum{1}) We take a copper circular rod, attach wired thermocouples to the lateral surface of the rod at prescribed locations, wrap the lateral surface with an insulation layer, and maintain the right tip of the rod at a  fixed temperature, for example $T_{tip}=40^{\circ}\mathrm{C}$ such that $T_{tip}<T_b$. Then we measure the steady-state temperature distribution along the bar at prescribed locations. 

Record and plot the steady-state temperature distribution along the rod. Find out whether the temperature has a linear distribution. Explain the mode of heat transfer in this case.

(\romannum{2}) For the same setting, replace the rod with other materials - aluminum and steel, and observe temperature changes due to various thermal conductivities. 

\subsection{Conduction and Convection}\label{sec13}

Now we remove the insulation layer around the lateral surface of the bar (circular rod) and expose it to the air. In this case, in addition to conductive heat transfer within the bar, there is convective heat transfer between the metal bar and the surrounding air flow. Again, we will measure and study the temperature distribution. There are several experimental conditions we can set up including replacing the ambient fluids from air to water, changing the tip conditions (prescribed temperature, prescribed flux or convective condition) of the bar, etc. Please conduct the following experiments:

(\romannum{1}) Keep the same tip condition as in the previous experiment, for example at a constant temperature $T_{tip}=40^{\circ}\mathrm{C}$, and investigate the temperature distribution due to different ambient fluid flows - air and water.

(\romannum{2}) When the lateral surface is exposed to air, investigate the effects of the following tip conditions on the temperature distribution within the bar. 

\begin{itemize}
\item A prescribed temperature: keep the tip at a constant temperature.
\item An insulated or adiabatic tip (zero flux): wrap the tip with an insulation layer.
\item A convective condition: expose the tip to air.
\end{itemize}


(\romannum{3}) Replace the circular rod by a slender rectangular bar made of the same material. See how the distribution of temperature changes due to varying geometries.

Record and tabulate the above experimental results from \S \ref{sec12} - \S\ref{sec13} in tables. The data will be used for further analysis and comparison.



\section{Modeling}

The above heat transfer process is governed by two physical laws: Fourier\textquotesingle s law and Newton\textquotesingle s law of cooling \cite{Bergman2011}. Assume the uniform slender bar has a length of $L$ and a constant cross-sectional area $A$. The temperature at the left end $(x=0)$ of the bar is fixed at $T(0)=T_b$. The right end $(x=L)$ of the bar is subject to varying boundary conditions. The lateral surface of the bar is exposed to a moving fluid at temperature $T_{fluid}$ with a convective heat transfer coefficient $h$. If we take a finite element of length $dx$ from the bar, this small slice of material undergoes conduction along the axial direction and convective cooling/heating on the external surface. Assume the conductive heat transfer rates on the two sides of the element are $q_x$ and $q_{x+dx}$ and the differential convection rate through the side surface is $dq_s$ (see Figure~\ref{fig1}).

\begin{figure}[H]
\centering
\includegraphics[width=5in,height=6in,keepaspectratio]{fig1}
\caption{\label{fig1}Schematic of energy balance for an element of the bar}
\end{figure} 

Applying the principle of conservation of energy to the finite element gives:

\begin{equation}\label{eq1}
q_x=q_{x+dx}+dq_s.
\end{equation}

From Fourier\textquotesingle s law, we know that $q_x=-kA\frac{dT}{dx}$ where $k$ is the thermal conductivity of the material. The heat rate at the right side of the element can be approximated using a 1st-order Taylor series expansion $q_{x+dx}=q_x+\frac{dq_x}{dx}dx$. The differential heat transfer rate $dq_s$ due to convection is determined by the Newton\textquotesingle s law of cooling $dq_s=hdA_s(T-T_{fluid})$ where $dA_s$ is the lateral surface area of the finite element.

Substituting the above quantities into \eqref{eq1}, we can derive the governing equation for temperature distribution $T(x)$ along a uniform slender bar as follows:

\begin{equation}\label{eq2}
\frac{d^2T}{dx^2}-\frac{hp}{kA}(T-T_{fluid})=0,
\end{equation}
where $p$ is the perimeter of the cross section and $p=\frac{dA_s}{dx}$. The temperature distribution $T(x)$ is affected by the thermal conductivity $k$ of the material, the perimeter $p$ and area $A$ of the cross section, the ambient fluid temperature $T_{fluid}$, and the convective heat transfer coefficient $h$.

If we define a new temperature variable $\Phi(x)\equiv T(x)-T_{fluid}$, \eqref{eq2} reduces to 

\begin{equation}\label{eq3}
\frac{d^2\Phi}{dx^2}-C^2\Phi=0,
\end{equation}
where $C^2=\frac{hp}{kA}$. It can be verified that the solution of \eqref{eq3} takes the form of $\Phi(x)=Ae^{Cx}+Be^{-Cx}$. Here the coefficients $A$ and $B$ can be evaluated using the boundary conditions at the two ends of the bar.

Now we need to specify the boundary conditions. According to the experiments, the temperature on the left $(x=0)$ is fixed by $T(0)=T_b$. Thus, $\Phi(0)\equiv T(0)-T_{fluid}=T_b-T_{fluid}\equiv \Phi_b$. In accordance with the experimental setups, the boundary condition on the right $(x=L)$ can be specified by one of the following:

\begin{itemize}
\item A prescribed temperature $\Phi(L)=\Phi_L$.
\item A prescribed flux $\frac{d\Phi}{dx}_{|{x=L}}=q_f$. If $q_f=0$, it is an insulated (adiabatic) tip.
\item A convective condition $-k\frac{d\Phi}{dx}_{|{x=L}}=h\Phi(L)$.
\end{itemize}

Then, the temperature distribution $\Phi(x)$ can be obtained from solving \eqref{eq3} with boundary conditions.



\subsection{Example}
Consider a slender, finite-length rod with a convective tip. Derive the expression of the temperature distribution along the rod.

From the above analysis, we know that the solution can be represented by $\Phi(x)=Ae^{Cx}+Be^{-Cx}$ with unknown coefficients $A$ and $B$. The two boundary conditions are $\Phi(0)= \Phi_b$ and $-k\frac{d\Phi}{dx}_{|{x=L}}=h\Phi(L)$. Substituting $\Phi(x)$ into boundary conditions gives:

\begin{eqnarray} \label{eq4}
\Phi_b &=&  A+B  \nonumber \\
-kC(Ae^{CL}-Be^{-CL}) &=& h(Ae^{CL}+Be^{-CL}).
\end{eqnarray}

The coefficients $A$ and $B$ can be determined by solving the equation system \eqref{eq4}. Then, the temperature distribution $\Phi(x)$ along the rod is 

\begin{equation}\label{eq5}
\Phi(x)=\frac{\cosh{C(L-x)}+\frac{h}{Ck}\sinh{C(L-x)}}{\cosh{CL}+\frac{h}{Ck}\sinh{CL}}\Phi_b.
\end{equation}
A plot of the temperature $T(x)=T_{fluid}+\Phi(x)$ along a rod with specified geometric and thermophysical conditions can be found by the curve of the analytical solution in \S \ref{sec3}.


\subsection{Exercises}\label{sec22}

~ ~ ~(\romannum{1}) Solve the temperature distribution along the rod when the tip temperature is prescribed by $T(L)=T_{tip}$.

(\romannum{2}) Solve the temperature distribution along the rod when the tip is insulated, i.e., $\frac{d\Phi}{dx}_{|{x=L}}=0$.

(\romannum{3}) Replace the above circular rod by a rectangular bar with the same cross-sectional area. Discuss how this shape change affects the temperature distribution.



\section{Numerical Simulations}\label{sec3}

The steady-state equation \eqref{eq2} with boundary conditions is solved numerically using the finite volume method \cite{Ferziger2002}. The computational domain is divided into $N$ node-centered control volumes or cells, as shown in Figure~\ref{fig2}. The grid has a uniform mesh size of $\Delta x=\frac{L}{N}$. Integration of \eqref{eq2} over a control volume gives 

\begin{figure}[H]
\centering
\includegraphics[width=5in,height=6in,keepaspectratio]{fig2}
\caption{\label{fig2}The grid}
\end{figure} 

\begin{equation}\label{eq6}
\int_S \frac{dT}{dx}\overrightarrow{n}dS-\int_{\mathcal{V}}C^2(T-T_{fluid})d\mathcal{V}=0,
\end{equation}
where the integrals in \eqref{eq6} are evaluated locally within each control volume. The first surface integral represents the heat flux through the surfaces of a control volume. The second volume integral can be regarded as a source term within the control volume. For example, evaluating \eqref{eq6} over control volume $3$ gives

\begin{equation}\label{eq7}
\big(\frac{dT}{dx}\big)_e-\big(\frac{dT}{dx}\big)_w-C^2(T_P-T_{fluid})\Delta x=0.
\end{equation}
Here the subscripts $e$ and $w$ refer to the east and west surfaces of the control volume. $P$ represents the center of the control volume. The two temperature gradients in \eqref{eq7} can be evaluated using linear approximations with temperatures at neighboring cells $2$ and $4$. Therefore, for internal cells from $2$ to $N-1$, replacing the temperature	 gradients in  \eqref{eq7} gives 

\begin{equation}\label{eq8}
\frac{T_E-T_P}{\Delta x}-\frac{T_P-T_W}{\Delta x}-C^2(T_P-T_{fluid})\Delta x=0,
\end{equation}
where the subscripts $E$ and $W$ refer to the center nodes of the east and west cells around the $P$ cell. For boundary cells $1$ and $N$, the temperature gradients in \eqref{eq7} need to be evaluated using the given boundary conditions. Then, we can establish a system of $N$ discretized equations. The temperature distribution $T(x)$ at each cell node can be calculated by solving the system of linear equations.

Figure~\ref{fig3} shows an interactive user interface for simulating the temperature distribution along a bar with varying tip conditions. On the left panel, the dimensions (length $L$, cross-sectional area $A$ and perimeter $p$) of the bar and the thermophysical properties (thermal conductivity $k$ and convective heat transfer coefficient $h$) can be entered. The ambient temperature condition is determined by the fluid temperature (Fluid T). For boundary conditions, the temperature on the left is prescribed at a base temperature (Base T). The condition on the right can be set up from the dropdown list and the corresponding tip temperature (Tip T) or flux (Tip Flux) is to be entered. On the right panel, users can input data from experiments. The middle panel shows the result of temperature distribution along the bar. The top plot shows the comparison of results from simulation (blue line with square markers), analytical solution (red line), and experiment (yellow circles). The error indicates the difference between the numerical solution and the exact solution. The bottom scales present the temperature distribution from the simulation.

\begin{figure}[H]
\centering
\includegraphics[width=5.5in,height=6in,keepaspectratio]{fig3}
\caption{\label{fig3}Temperature distribution MATLAB App (downloadable on SIMIODE website)}
\end{figure} 

\subsection{Numerical Exercises}\label{sec31}
~ ~ ~(\romannum{1}) Investigate the interaction between conduction and convection by changing the values of thermal conductivity $k$ and convective heat transfer coefficient $h$. If $k\gg h$, what is the temperature distribution when the tip is maintained at a constant room temperature? Explain the corresponding physical condition for this case. What happens if $k\ll h$? How is heat transferred in the bar?

(\romannum{2}) Assume that the ambient fluid temperature remains constant. If a circular rod is under the convective tip condition, how could one increase the heat transfer rate by changing other parameters? Then, replace the above circular rod by a square bar with the same cross-sectional area. Discuss how this shape change affects the temperature distribution, especially the tip temperature.


(\romannum{3}) What is the difference between the adiabatic and the convective tip conditions? To reach the same tip temperature, how are these two conditions related to the length of the bar?

(\romannum{4}) How could one improve the accuracy of a numerical solution, especially when sharp gradients exist in the solution? 



\section{Comparison of Results}

We investigate the steady-state temperature distribution along slender metal bars through experimental, analytical, and numerical approaches. Now we want to compare the analytical and numerical solutions with the temperature data collected from selected experiments. 
Use the MATLAB App discussed in \S \ref{sec3} to compare results in the following scenarios:
\begin{itemize}
\item \S \ref{sec12} (\romannum{1}) The temperature distribution in a copper rod due to conduction.
\item \S \ref{sec13} (\romannum{2}) The temperature distribution in a rod with varying tip conditions.
\item \S \ref{sec13} (\romannum{3}) The change of temperature distribution due to varying geometries.
\end{itemize}
The comparison of results will provide more insights on the heat transfer process and show correlation between the modeling and the experiments. 




\begin{thebibliography}{98}

\bibitem{Bergman2011} Bergman, T. L., F.P.~Incropera, D.P.~DeWitt, and A.S.~Lavine. 2011. {\it Fundamentals of heat and mass transfer}. New York: John Wiley \& Sons.

\bibitem{Ferziger2002}Ferziger, J. H., M.~Peri{\'c}, and R.L.~Street. 2002. {\it Computational methods for fluid dynamics}. New York: Springer.


\end{thebibliography}

\teacher{
\section*{NOTES FOR TEACHERS}

This project can be used in engineering courses such as heat transfer, thermal fluids lab, and computational fluid dynamics. The full version requires students to investigate the temperature distribution along a uniform slender bar through experimental, analytical and numerical approaches. Several parameters and boundary conditions are to be investigated, which makes the study comprehensive and challenging. Students are expected to complete the entire project within one week including experimental observation, data collection, modeling and simulation. A written report is required to address all questions and exercises in the project. Instructors also have the freedom to choose certain modules in this project such as only to teach the analytical and numerical studies and ignore the experiments. 

The prerequisites of this project are some familiarity with undergraduate-level heat transfer, some exposure to ordinary differential equations, and numerical analysis. Specifically, students need to understand the basics of conduction and convection, conservation of energy, and Fourier\textquotesingle s law and Newton\textquotesingle s law of cooling. Students are expected to know some common techniques in solving ordinary differential equations analytically. Also, students need to understand the basic procedures in solving diffusion equations with source, numerically using the finite volume method. In this project, students are not required to develop a code to conduct simulations. Instead, simulations can be done using the MATLAB App directly even for users without MATLAB coding experience. The MATLAB App used in this project can be downloaded from the SIMIODE website for this Modeling Scenario. Instructors may also ask the students to develop their own versions of codes to have a deeper understanding of the numerical schemes.

\setcounter{equation}{0}
\renewcommand{\theequation}{C.\arabic{equation}}

\subsection*{Hints for Selected Exercises}

\S \ref{sec22}(\romannum{1})
We know that the solution can be represented by $\Phi(x)=Ae^{Cx}+Be^{-Cx}$ with unknown coefficients $A$ and $B$. The two boundary conditions are $\Phi(0)= \Phi_b$ and $\Phi(L)\equiv T(L)-T_{fluid}=T_{tip}-T_{fluid}\equiv \Phi_L$. Substituting $\Phi(x)$ into boundary conditions gives:

\begin{eqnarray} \label{eq9}
\Phi_b &=&  A+B  \nonumber \\
\Phi_L &=& Ae^{CL}+Be^{-CL}.
\end{eqnarray}

The coefficients $A$ and $B$ can be determined by solving the equation system \eqref{eq9}. Then, the temperature distribution along the rod is 

\begin{equation}\label{eq10}
\Phi(x)=\frac{\sinh{C(L-x)}+\frac{\Phi_L}{\Phi_b}\sinh{Cx}}{\sinh{CL}}\Phi_b.
\end{equation}

\S \ref{sec22}(\romannum{2})
This problem can be solved in a similar fashion. The two boundary conditions are $\Phi(0)= \Phi_b$ and $\frac{d\Phi}{dx}_{|{x=L}}=0$. Substitute the general expression of $\Phi(x)$ into boundary conditions:

\begin{eqnarray} \label{eq11}
\Phi_b &=&  A+B  \nonumber \\
0 &=& C(Ae^{CL}-Be^{-CL}).
\end{eqnarray}

The coefficients $A$ and $B$ can be determined by solving the equation system \eqref{eq11}. Then, the temperature distribution along the rod is 

\begin{equation}\label{eq12}
\Phi(x)=\frac{\cosh{C(L-x)}}{\cosh{CL}}\Phi_b.
\end{equation}

\S \ref{sec31}(\romannum{1})
For a bar with fixed temperatures on both ends, the effect of conduction dominates if $k\gg h$. The convective heat transfer between the lateral surface of the bar and the ambient fluid flow can be ignored. It is similar to the insulated condition described in \S \ref{sec12}. The temperature distribution has a linear trend. If $k\ll h$, it indicates the material has a very low thermal conductivity and there is a strong convection between the lateral surface of the bar and the ambient fluid flow. 

\S \ref{sec31}(\romannum{2})
The heat transfer rate can be increased by replacing the ambient fluid by fluids with higher convection heat transfer coefficients. Changing the circular rod to a square bar with the same cross-sectional area will slightly lower the temperature at the tip. This can be verified from the MATLAB App.

\S \ref{sec31}(\romannum{3})
The convective tip condition assumes an energy balance between conduction and convection at the tip, while the adiabatic tip condition assumes convective heat transfer at the tip is negligible. Therefore, a bar with the adiabatic tip needs to be longer than a bar with the convective tip to reach the same tip temperature.

\S \ref{sec31}(\romannum{4})
The accuracy of the numerical solution for this one-dimensional problem can be improved by using a finer uniform mesh or a non-uniform mesh where the grid is finer when sharp gradients exist.


\subsection*{Discretized Equations for Boundary Cells in the Finite Volume Scheme}

For cell $1$, the temperature is fixed at $T_b$ on the left:

\begin{equation}\label{eq13}
\frac{T_E-T_P}{\Delta x}-\frac{T_P-T_b}{\Delta x/2}-C^2(T_P-T_{fluid})\Delta x=0.
\end{equation}

For cell $N$, the discretized equation depends on the boundary condition on the tip:
(\romannum{1}) a prescribed temperature $T_{tip}$,
\begin{equation}\label{eq14}
\frac{T_{tip}-T_P}{\Delta x/2}-\frac{T_P-T_W}{\Delta x}-C^2(T_P-T_{fluid})\Delta x=0;
\end{equation}
(\romannum{2}) a prescribed flux $q_f$,
\begin{equation}\label{eq15}
\frac{q_f}{k}+\frac{T_P-T_W}{\Delta x}+C^2(T_P-T_{fluid})\Delta x=0;
\end{equation}
(\romannum{3}) a convective condition $-k\frac{T_{tip}-T_P}{\Delta x/2}=h(T_{tip}-T_{fluid})$,
\begin{equation}\label{eq16}
\frac{2h}{h\Delta x+2k}(T_{fluid}-T_P)-\frac{T_P-T_W}{\Delta x}-C^2(T_P-T_{fluid})\Delta x=0.
\end{equation}

}





\end{document}
