\documentclass[11pt,twoside]{article}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{array,longtable}
\usepackage{theorem}

\parindent =0pt
\parskip = 12pt



\newcommand{\pr}[2]{\ensuremath{\frac{\partial #1}{\partial #2}}}
\newcommand{\bo}[1]{ \mbox{ \boldmath $#1$ } }
\newcommand{\pdt}[1]{{#1}^{\prime}(t)}
\newcommand{\dt}[1]{\frac{\partial #1}{\partial t}}
\newcommand{\Dt}{\frac{D}{D t}}
\newcommand{\intO}[1]{\int_{\Omega_t}{#1}\, d\Omega_t}
\newcommand{\intdO}[1]{\int_{\partial \Omega_t}{#1}\, dS}
\newcommand{\divd}[1]{\nabla \cdot \! #1 \! }
\newcommand{\ir}{\mathbf{i}_r}
\newcommand{\itheta}{\mathbf{i}_{\theta}}
\newcommand{\ds}{\displaystyle}
\renewcommand{\vec}[1]{{\bf {#1}} }
\renewcommand{\descriptionlabel}[1]{{\em{#1}}}

\begin{document}
\bibliographystyle{plain}
\bibliographystyle{unsrt}


\pagestyle{plain}

\date{}



{\huge \bf Mandatory Assignment No.2 \\for INF-MAT2351\\ \large
  Deadline Noon May 7, 2015.}
\\ \ \\

{\it Every student should do the programming by herself/himself and
  must submit her/his own report as a PDF file. Use LaTeX, or a
  similar tool, to typeset your report. The program codes should also
  be submitted as an appendix to the report. Submit the assignment in Devilry.}



\section{Introduction}%\begin{abstract}

In this project we study the cooling of a can of beer in a refrigerator. We assume that the can of beer has a temperature of $24$ degree Celsius ($^oC$) originally, and the fridge is $4 ^o C$. We will try to simulate how the heat diffuses out of the can. We are interested in answering the following question: How long time does the beer have to stay in the fridge to be cold enough to be enjoyed.



To do this simulation accurate, will require knowledge to physical parameters in the beer, air and aluminum (in the can). We also need to know the flow of air in the fridge and the flow of beer in the can.  This would be a very complicated simulation and we will therefore in this project only simulate the heat flow in the beer, where we assume that the beer is not moving.

%\end{abstract}



\section{Model}



A can of beer is usually almost like a cylinder. The Norwegian $\frac 12$-liter cans have radius $r_c=32 mm$ and the height is about $h_c=160 mm$, see Figure \ref{fig:beer:1}.



\begin{figure}[ptb]
  \begin{center}
    %\includegraphics[height=2.367in, width=2.2618in]{can.eps}
\includegraphics[height=2.367in, width=2.2618in]{can.png}
\caption{\label{fig:beer:1} The can of beer is approximately a cylinder with height $h_c=160mm$ and radius $r_c=32mm$}
  \end{center}
\end{figure}


The three dimensional heat conduction equation reads

\begin{align}
  \varrho c_v {\partial T\over\partial t} = \nabla\cdot\left( k\nabla
  T\right) + f(x,y,z,t),
\label{beer1}
\end{align}

where $T(x,y,z,t)$ is the temperature, $\varrho$ is mass density, $c_v$ is the specific heat capacity, $k$ is heat conductivity and $f(x,y,z,t)$ is the generated heat per unit time and per unit volume, see Equation (7.25) in \cite{A}.



Since most of the beer is water, we assume that the physical properties of beer are equal to water. The physical properties of water are as follows:

\begin{itemize}

\item $\varrho = 1 000 \frac {kg}{m^3}$ ~~(The mass density of water)

\item $c_v = 4200 \frac J{kg\, K}$ ~~(The specific heat capacity)

\item $k = 0.58 \frac W {m\, K}$ ~~(The approximate heat conductivity) 

      %(k = 0.556-0.602 nr t = 5-25 ^oC)

\end{itemize}

All these numbers will vary with pressure and temperature, but these are the approximate numbers at room temperature and with earth pressure.



\section{Simplifications}

We assume that the physical properties, $\varrho$, $c_v$ and $k$, are constant and that there is no generated heat in the beer, $f(x,y,z,t)=0$.

Equation \eqref{beer1} can therefore be simplified to  


\begin{align}
{\partial T\over\partial t} =  \frac k{\varrho c_v}  \nabla\cdot\left( \nabla  T\right),
\label{beer2}
\end{align}


where 

\[
\nabla\cdot\left( \nabla T\right) = \frac {\partial^2 T}{\partial x^2}
+ \frac {\partial^2 T}{\partial y^2} + \frac {\partial^2 T}{\partial z^2}.
\]




\subsection{Invariance in $z$-direction}

We assume there is no variation in $z$. This means that $\frac {\partial
T}{\partial z}=0$ or

$$T=T(x,y,t).$$

Equation \eqref{beer2} can then be simplified to

\begin{align}
   {\partial T(x,y,t) \over\partial t} =  \frac k{\varrho c_v}  \nabla\cdot\left( \nabla  T(x,y,t) \right),
\label{beer3}
\end{align}

where 

$$\nabla\cdot\left( \nabla T(x,y,t)\right) = \frac {\partial^2
  T(x,y,t)}{\partial x^2} + \frac {\partial^2 T(x,y,t)}{\partial y^2}.$$



\begin{figure}
  [ptb]
  \begin{center}
    \includegraphics[height=2.367in, width=2.2618in]{circ.png}
    \caption{We simplify our simulations to a disc with radius $r_c=32mm$} 
    \label{fig:beer:2}
  \end{center}
\end{figure}





In practice the assumption $\frac {\partial T}{\partial z}=0$, will not be true. However, if the bottom and the top of the can of beer are insulated, it will be a reasonable assumption. In that case, no heat can flow out of the bottom and the top and physically it seems reasonable that all changes in temperature are in $(x,y)$-directions.



Another situation where the above assumption is reasonable is if the height is much larger than the radius. Then the temperature variations in z-direction will be small if you are far away from the bottom and the top.



%\begin{itemize}

%\item{\bf a)} Give an example of a physical situation where invariance in

%  $z$-direction is reasonable.

%\end{itemize}



\subsection{Polar coordinates}



We now have a two dimensional problem with circular geometry. The solution is therefore assumed to be radially symmetric, i.e.

$$T=T(r,t),$$

where the coordinate system is placed with origo at the center of the circle. 







Physically this is reasonable since: initially the temperature is constant, which is a radially symmetric function, and the boundary ($r=r_c$) is symmetric and therefore heat flow will be symmetric.





It can be shown that

\begin{align}
  \nabla \cdot (\nabla T(r,t))=\frac 1 r \frac {\partial }{\partial r}(r
  \frac {\partial T(r,t)}{\partial r}), 
\end{align}

when $T$ is a radially symmetric function $T=T(r,t)$. Thus \eqref{beer3} can be written

\begin{align}
   {\partial T(r,t) \over\partial t} =  \frac k{\varrho c_v}  \frac 1r \frac {\partial }{\partial r}(r  \frac {\partial T(r,t)}{\partial r}).
\label{beer4}
\end{align}









\subsection{Initial conditions}

We assume that the temperature is $24^oC$ initially, i.e.

\begin{align}
  T(r,0)=T_0 \quad (=24^oC).
\end{align}


\subsection{Boundary conditions}



If the fridge is large compared to the can of beer, we can assume that the air in the fridge is constant at $4^oC$, $T_f=4$. 

Newton's law of cooling says that the heat flow out of the boundary is proportional to the difference between the temperature in the beer and the temperature in the fridge (the aluminium in the can is very thin), i.e.

\begin{align}
  -k \frac {\partial T(r_c,t)}{\partial r} = h_T(T(r_c,t)-T_f).
  \label{Bcond1}
\end{align}



Here the heat transfer coefficient will usually be approximately

\begin{align}
  h_T=5\frac W{m^2\, K}, %(3-25 ved stille luft)
\end{align}

when the air is not moving. However, if the air in the fridge starts to move the coefficient will increase. If the air moves very much, we get the boundary condition: 

\begin{align}
  T(r_c,t)=T_f.
  \label{Bcond2}
\end{align}



\begin{description}

\item{\bf a)} Show that \eqref{Bcond2} follows from \eqref{Bcond1} when

  $h_T \rightarrow \infty$. (Hint: assume that $|\frac {\partial T(r_c,t)}{\partial r}| < \infty$)

\end{description}







It can be shown that the boundary condition for $r=0$ is

\begin{align}
  \frac {\partial T}{\partial r}(0,t)=0.
\end{align}





We assume that \eqref{Bcond2} is valid and write the full system

\begin{align}
  {\partial T(r,t) \over\partial t} &=\frac k{\varrho c_v} \frac 1r   \frac {\partial }{\partial r}(r  \frac {\partial T(r,t)}{\partial r}),\quad (r,t)\in (0,r_c)\times (0,t_{end}) \label{model1}\\
  \frac {\partial T}{\partial r}(0,t)&=0\label{model2}\\
  T(r_c,t)&=T_f \label{model3}\\
  T(r,0)&=T_0\label{model4}
\end{align}

If Newton's law of cooling is used for the boundary, then \eqref{model3} is replaced with 

\begin{align}
  -k \frac {\partial T(r_c,t)}{\partial r} = h_T(T(r_c,t)-T_f).
\end{align}

Above $r_c=0.032 \, m$, $T_f = 4 ^oC$, $T_0 = 24 ^oC$ and $t_{end}$ is
the end time.



\section{Scaling}



The problem is scaled in the following way

\begin{align}
  \bar {r} = \frac r r_c, \quad \bar {T} = \frac {T-T_f}{T_0-T_f} 
  \quad \text{and} \quad \bar {t} = \frac t t_c ,
\end{align}

where $r_c=0.032 \, m$, $T_f = 4 ^oC$, $T_0 = 24 ^oC$ and $t_c$ is the characteristic time (see question {\bf e)}).



\begin{description}

\item{\bf b)} Show that 
  \begin{align}
    \frac {\partial T}{\partial t} = \frac {T_0-T_f}{t_c} \frac{\partial \bar {T}}{\partial \bar {t}}.
  \end{align}

\item{\bf c)} Show that 

  \begin{align}
    \frac 1r \frac {\partial }{\partial r} \left( r \frac {\partial T}{\partial r} \right) = \frac {T_0-T_f}{r_c^2} \frac 1{\bar{r}} \frac {\partial} {\partial \bar {r}} \left(\bar {r}  \frac {\partial \bar {T}} {\partial \bar {r}} \right).
\end{align}

\item{\bf d)} Show that the Newton's law of cooling on scaled form can be written
  \begin{align}
    -\frac {k}{h_T r_c} \frac {\partial \bar {T}(1,\bar {t})} {\partial \bar {r}}= \bar {T} (1,\bar {t}).
  \end{align}



\item{\bf e)} Choose $t_c=\frac {\varrho c_v r_c^2}k$ and show that \eqref{model1}-\eqref{model4} can be written

  \begin{align}
    \frac {\partial \bar{T}}{\partial \bar{t}}&=\frac 1{\bar{r}} \frac {\partial }{\partial \bar{r}} \left( \bar{r} \frac {\partial \bar{T}}{\partial \bar{r}} \right),\quad (\bar{r},\bar{t})\in (0,1)\times (0,\bar{t}_{end}) \label{scaled1} \\
  \frac {\partial \bar{T}}{\partial \bar{r}}(0,\bar{t})&=0 \label{scaled2}\\
  \bar{T}(1,\bar{t})&=0 \label{scaled3}\\
  \bar{T}(\bar{r},0)&=1.  \label{scaled4}
\end{align}

%where $u=\bar{T}$, $x=\bar {r}$ and $t=\bar{t}$

\end{description}

%\item{\bf f)} Calculate $\bar{t}_{end}$ expressed with 







\subsection{Discretization}

In the following we define $u=\bar{T}$ and  $x=\bar {r}$, and write the system on scaled form:

  \begin{align}
    \frac {\partial u}{\partial \bar{t}}&=\frac 1x \frac {\partial }{\partial x} \left( x \frac {\partial u}{\partial x} \right),\quad (x,\bar{t})\in (0,1)\times (0,\bar{t}_{end}) \label{model:eq1} \\
  \frac {\partial u}{\partial x}(0,\bar{t})&=0 \label{model:eq2}\\
  u(1,\bar{t})&=0 \label{model:eq3}\\
  u(x,0)&=1\label{model:eq4},
\end{align}

when $h_T=\infty$. When $h_T < \infty$, \eqref{model:eq3} is replaced with

\begin{align}
  -\frac {k}{h_T r_c} \frac {\partial u}{\partial x}(1,\bar{t}) = u (1,\bar{t}).  \label{model:eq5}
\end{align}

Note that these equations are on dimensionless form (only numbers)



To discretize this problem we choose a number of spatial discretization points, $n$, and a number of temporal discretization points, $m$. We define the discretization parameters

\begin{align*}
  \Delta x = \frac 1 {n-1} \quad \text{and} \quad \Delta t = \frac 1 m.
\end{align*}



The discretization points are

\begin{align*}
  (x_i,t_\ell)=((i-1) \Delta x, \ell \Delta  t), \quad i=1,\dots ,n, \quad \ell = 0,\dots, m.
\end{align*}





The numerical solution $u_i^\ell$ approximates the exact solution $u(x_i,t_\ell)$ in the discretization points, i.e.

\begin{align}
  u_i^\ell \approx u(x_i,t_\ell), \quad i=1,\dots ,n, \quad \ell = 0,\dots, m.  
\end{align}





\begin{description}

\item{\bf f)} Write a numerical scheme for solving Equation \eqref{model:eq1} in the grid points. (The scheme should be a finite difference scheme on explicit Euler form, as described in Chapter 7.4 in \cite{A}.)







The boundary condition \eqref{model:eq2} shall in this project be discretized by

\begin{align}
  \frac {u^\ell_1-u^\ell_2}{\Delta x}=0
\end{align}

which gives

\begin{align}
  u^\ell_1=u^\ell_2, \quad \text{for} \quad \ell=0,\dots ,m.
\end{align}



\item{\bf g)} Show how the boundary conditions \eqref{model:eq3} and \eqref{model:eq5} and the initial condition \eqref{model:eq4} can be incorporated into the finite difference scheme.



\item{\bf h)} Write an algorithm for solving \eqref{model:eq1}-\eqref{model:eq4} with the finite difference method.



\item{\bf i)} Write an algorithm for solving \eqref{model:eq1}-\eqref{model:eq2} and \eqref{model:eq4}-\eqref{model:eq5} with the finite difference method.

\end{description}





\subsection{Implementation}



In the simulations we do not need to store all the values for $u_i^\ell$. We are just interested in the temperature in the last time step.



\begin{description}

\item{\bf j)}  Write the algorithms in {\bf h)} and {\bf i)}, in such a way that only the storage of $2n$ numbers are needed for $u$.



\item{\bf k)} Implement the algorithms from {\bf j)} in one of the following languages: Matlab, Python, Java, C++, C or Fortran.



\item{\bf l)} Based on numerical experiments, discuss stability. Discuss

  briefly accuracy (we define the solution to be accurate enough when

  decreasing $\Delta x$ and $\Delta t$ gives small changes in the solution).



\end{description}



For the next questions, use at least $40$ spatial discretization points

($n>40$), in order to get sufficient accurate solutions. Answer the

following questions for both two codes. 



\begin{description}

\item{\bf m)} Make the code such that it stops the simulations when

  $u_1^\ell<0.1$. What is the smallest $\bar{t}=\Delta t \ell$ such that

  $u_1^\ell < 0.1$.



\item{\bf n)} Make the code such that it stops the simulations when

  $\frac 1n \sum_{k=1}^n 2 \, u_k^\ell \, x_k < 0.1$. What is the

  smallest $\bar{t}=\Delta t \ell$ such that $\frac 1n \sum_{k=1}^n  2\, u_k^\ell \, x_k < 0.1$.  

\end{description}





Now we go back to the original problem.



\begin{description}

\item{\bf o)} At what time (nonscaled time) is the temperature at the centre of the can of beer $6^oC$?



\item{\bf p)} At what time is the average temperature in the beer $8^oC$?



\item{\bf q)} At what time is the average temperature in the beer $6^oC$?


\item{\bf r)} In practice the value of the heat transfer coefficient,

  $h_T$, is hard to determine. Solve {\bf o)}, {\bf p)} and {\bf q)} when $h_T=10 \, \frac {W}{m^2 K}$ and $h_T=18 \, \frac {W}{m^2 K}$.



\end{description}











\begin{thebibliography}{99}



\bibitem{A} A. Tveito, H. P. Langtangen, B. F. Nielsen and X. Cai. Elements of Scientific Computing.



\end{thebibliography}







\end{document}
























