
% 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}
\renewcommand{\baselinestretch}{1.3}
\newcommand{\pwrt}[2]{$\frac{\partial #1}{\partial #2}$}
\newtheorem{theorem}{Theorem}
\newcommand{\dt}{\Delta t}
%
% 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{hyperref}
\usepackage{url}
\usepackage{epsfig}
\usepackage{float}
%
\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\\Skin Burn Model - Numerical Methods}




\markboth{Skin Burn Model - Numerical Methods}{Skin Burn Model - Numerical Methods} 

\author{Suruchi Singh\\
Department of Mathematics\\ 
Aditi Mahavidayalaya\\
University of Delhi\\
Delhi 110039 INDIA\\
ssuruchi2005@yahoo.co.in}

\makeStitlePDFLaTex

\begin{abstract}
The heat equation is an important partial differential equation (PDE) which describes the distribution of heat in a given region  over time. Here we learn to solve a heat equation numerically. It is difficult to study the behavior of temperature in problems with interfaces analytically so numerical methods play an important role in solving these. The idea presented here can be used to solve a wide variety of models for non-homogenous inner structure models with partial differential equations. 
\end{abstract}



\subsection*{SCENARIO DESCRIPTION - BACKGROUND READING}  %%%%%%%%%%%%%%%%%%%
The heat equation arises in various fields such as image processing, hyperthermia treatment, chemical reactions, etc.and it can be applied to the study of skin burns.  The  skin is the largest organ of the human body and  plays an important role such as defense, sensory effects, and thermoregulation. Skin burns are one of the most devastating injuries encountered in a human's daily life. These injuries usually result from heat, electricity, sunburn, radiation or chemicals. Skin tissue has a complicated multilayer structure namely epidermis, dermis and subcutaneous (see Figure 1) and the complex boundary and interfacial conditions raise the difficulty of solving the problem analytically. Each layer possesses different thermo-physical properties. Using Numerical Methods we develop algorithms which enable a computer to find the solution of complicated problems.
\begin{figure}[H]
\centering
\includegraphics[width=1\textwidth]{skinModel.pdf}
\caption{Triple layered skin structure}
\label{Figure3}
\end{figure}

\subsection*{HEAT ENERGY MODEL}%%%%%%%%%%%%%%%%%%%
In any material the thermal energy transfers from the region of higher temperature to the region of lower temperature. Consider a uniform rod of length $l$ with density $\rho$ and  specific heat $c$ with cross-sectional area of the rod to be $A$.  $\kappa$ is called the thermal conductivity.

Let $u(x, t)$ be the temperature of the thin slice of the rod between $x$ and $x +\Delta x$. We offer a derivation of  the heat equation which is based on Fourier's law and the Conservation of Energy Principle. Basically, we will compute the change in heat in a small section of the rod in two different manners and equate these to form a partial differential equation.

Now, the heat energy of the segment between $x$ and $x +\Delta x$ is computed by
\begin{equation}
c \rho A \Delta  x u(x, t) \\
\end{equation}

\noindent while Fourier's law of heat transfer gives the following:


\begin{equation}
\frac{\text {Rate of heat transfer}}{\text{area}}=-\kappa\Delta u(x,t)\, .\\
\end{equation}

Using the Conservation of Energy Principle,
\begin{equation}
       \text {change of heat energy in the segment in time $\Delta t$ }= \text {energy in} - \text{energy out}\, ,
 \end{equation}
 
\noindent from  (1)-(3), we obtain a difference equation which is then simplified with division of both sides by the factor $ A \Delta t \Delta x$.
\begin{equation}
c \rho A \Delta  x (u(x, t+\Delta t)-u(x, t))=\Delta t A(-\kappa\Delta u(x,t)+\kappa\Delta u(x+\Delta x,t))
\end{equation}
\begin{equation}
c \rho  \frac{(u(x, t+\Delta t)-u(x, t))}{\Delta t }=\frac{(-\kappa\Delta u(x,t)+\kappa\Delta u(x+\Delta x,t))}{\Delta  x}
\end{equation}
Taking the limit as  $\Delta t$ and  $\Delta x$ both tend to $0$ gives the Heat Equation in our case, a partial differential equation:
 \begin{equation}
(\kappa u_{x})_{x}= \rho c\frac{\partial u(x,t)}{\partial t}
\end{equation}

\subsection*{SKIN BURN MODEL}
To predict the behavior of elevated temperature in skin burn treatment, a model is as follows (see \cite{klinger1974heat,liuc2008finite,pennes1948analysis,shih2007analytical,weinbaum1985new}):
\begin{equation}
(\kappa u_{x})_{x}-w_bc_bu(x,t)+f= \rho c\frac{\partial u(x,t)}{\partial t}
\end{equation}
where $\rho$ and $c$  denote density and specific heat of the tissue; $c_b$ is the specific heat of the blood; $w_b$ blood perfusion rate; $f$ is internal heat source. ($f$ is the heat generated due to metabolism)

\subsection*{NUMERICAL APPROACH TO SOLVE THE MODEL}%%%%%%%%%%%%%%%%%%%

First we discretize the region $ \{(x, t) : a \le x\le b, t\ge 0\}$ into $N$ subintervals in space direction with spacing of $h =\frac{(b-a)}{N}> 0$, so that $x_{l+1}-x_{l}=h$. Let $\Delta t$ be the time step. For $l =0,1,2,...,N$  and $j > 0$, we let the solution $u$ at the grid point $(x_l, t_j)$ be denoted by $u^j_l$.
We can now approximate the  time derivative using  Forward Euler's method with,
\begin{equation}
u_{{t}_{l}}^j=\frac{u^{j+1}_l-u^{j}_l}{\Delta t}
\end{equation} 
and  spatial derivative $u_{xx}$ by central differencing 
 \begin{equation}
u_{{xx}_{l}}^j=\frac{u^j_{l+1}-2u^j_l+u^j_{l-1}}{h^2}\, .
\end{equation} 
Substituting (8) and (9) in PDE (7),\\
\begin{equation}
\rho c\frac{u^{j+1}_l-u^{j}_l}{\Delta t}=\kappa \frac{u^j_{l+1}-2u^j_l+u^j_{l-1}}{h^2}-w_bc_bu^j_l+f^j_l
\end{equation} 
\noindent we obtain a system of algebraic equations:
\begin{equation}
u^{j+1}_l=u^{j}_l+\frac{\kappa\Delta t}{\rho c} \frac{u^j_{l+1}-2u^j_l+u^j_{l-1}}{h^2}-\frac{w_bc_b\Delta t}{\rho c}u^j_l+\frac{\Delta t}{\rho c}f^j_l\, .
\end{equation} 
Suppose we have initial temperature as\\
 \begin{equation}
u(x,0)=T_0
\end{equation} 
\begin{equation}
u^{0}_l=T_0  \text{ for all }  l
\end{equation} 
and temperature at the boundaries as
\begin{equation}
u(a,t)=T_1(t),
u(b,t)=T_2(t)
\end{equation} 
\begin{equation}
u^{j}_o=T_1(t_j) ,u^{j}_N=T_2(t_j)  \text{ for all }  j
\end{equation} 

\subsection*{Initial Actions}

For  numerical experiments we take thickness of the tissue as 12.08 mm, density $\rho =0.001g/mm^3$,  specific heat $c= 4.2J/g^oC$ , $c_b = 4.2J/g^oC$, perfusion rate  $w_b =5*10^{-7} g/mm^3/sec$,  thermal conductivity $\kappa=5.2*10^{-4}$ and metabolic heat $f$=0

\begin{enumerate}
\item Actions 1: When skin comes in contact with hot steel plate, skin surface temperature is constant. Taking elevated temperature $12^o$C on the skin boundary, solve the system (11) using  MATLAB or any other software.   
\item Actions 2: Plot the graph of temperature versus depth at time 30 sec as shown in Figure 2.
\item  Actions 3: Observe how temperature behaves at different depths. Figure 3 shows the behavior of temperature  at depths 0.72mm and 2.08mm as time progresses.
\end{enumerate}
\begin{figure}[H]
\centering
\includegraphics[width=1\textwidth]{newfig2.jpg}
\caption{Temperature versus depth graph }
\label{Figure2}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=1\textwidth]{Newfig3.jpg}
\caption{Time versus temperature graph at two different depths}
\label{Figure3}
\end{figure}

\subsection*{FURTHER MODIFICATIONS IN THE MODEL}%%%%%%%%%%%%%%%%%%%
There are different layers of skin with different physical and thermal properties. Suppose two layers meet at grid point $x=x_{l_1}$. The interface separates the region into two parts $\Omega^-=(a, x_{l_1})$ and $\Omega^+=(x_{l_1}, b)$. 
Let $\kappa$, $\varepsilon $, and $\sigma$ be piecewise constant functions, i.e. 
\begin{equation*}
\kappa=
\begin{cases}
\kappa^- &\text{in}\ \Omega ^-
\\
\kappa^+ & \text{in}\ \Omega ^+
\end{cases}
, \quad
\varepsilon = 
\begin{cases}
\varepsilon ^- & \text{in}\ \Omega ^-
\\
\varepsilon ^+ &\text{in}\ \Omega ^+
\end{cases}
, \quad
\sigma = 
\begin{cases}
\sigma ^-  &\text{in}\ \Omega ^-
\\
\sigma ^+  &\text{in}\ \Omega ^+ \, .
\end{cases}
\end{equation*}
Physically the temperature should be continuous at the interface also which gives  
\begin{equation}
u^+-u^-=\lim_{x\to \ {l_1}^+} u(x)-\lim_{x\to \ {l_1}^-} u(x)=0\, .
\end{equation}
In addition, continuity of heat flux is applied to the boundary of adjacent layers,
\begin{equation}
\kappa^+ u_{x}^+=\kappa^- u_{x}^-\, .
\end{equation}
Now using  Taylor's expansion about $x_l$, we have 
\begin{eqnarray}
u_{l_1-1}&=&u^--u_x^-h+\frac{1}{2}u_{xx}^-h^2+O(h^3)\, \quad \text{and} \\
u_{l_1+1}&=&u^++u_x^+h+\frac{1}{2}u_{xx}^+h^2+O(h^3)\, .
\end{eqnarray}
Solving equations (7), (16)-(19) for $u_x^-$ and  $u_x^+$, we obtain
\begin{eqnarray}
({\rho^+ c^+}+{\rho^- c^-})u^{j+1}_{l_1} & = &    ({\rho^+ c^+}+{\rho^- c^-})u^{j}_{l_1}+2\Delta t \frac{\kappa^+u^j_{{l_1}+1}-(\kappa^++\kappa^-)u^j_{l_1}+\kappa^-u^j_{{l_1}-1}}{h^2}\\ \nonumber
&& -\Delta t (w_b^-c_b^-+w_b^+c_b^+)u^j_{l_1}+\Delta t (f^++f^-)
\end{eqnarray}


\subsection*{Activities}
\begin{enumerate}
\item Solve the model numerically using MATLAB or any other software and predict the behavior of temperature at different layers when constant temperature  $12^o$C  is given to the skin boundary. The values of the parameters is given in Table 1.
\begin{table}[!htbp]
\begin{center}

\begin{tabular}{|l|c|c|c|c|}
 \hline
 & Blood & Epidermis &Dermis &Subcutaneous \\
 \hline
density [g/mm$^3$] & &0.0012&0.0012&0.001\\
\hline
specific heat [J/g$^o$C] &4.2&3.6&3.4&3.06\\
\hline
thermal conductivity [w/mm$^o$C] & &$2.6*10^{-4}$ &$5.2*10^{-4}$ &$2.1*10^{-4}$\\
\hline
perfusion rate [g/mm$^3$/sec] & &0 &$5*10^{-7}$ &$5*10^{-7}$\\
\hline
thickness [mm] & &0.08&2&10\\
\hline
\end{tabular}

\caption{Values of the parameters used for the numerical experiments}
\end{center}
\end{table}
\item Plot the graph of temperature versus depth.
\item Observe how temperature behaves at different depths.
\item When skin gets burned one can get relief by putting the skin under the cold water tap or under any cold beverage available at the moment. What conditions  need to be imposed so that the model predicts the effect of cool pack on the each layer.
\item Examine the case when thermal conductivity is variable. 
\item  Predict the behavior of temperature at different depths when linear temperature is imposed on the skin boundary.
\item Apply this technique on other non-homogeneous models.

\end{enumerate}

\begin{thebibliography}{99}

\bibitem{klinger1974heat} Klinger, H.G. 1974. 
 Heat transfer in perfused biological tissue-I: General theory.
{\it Bulletin of Mathematical Biology\/}. 36: 403--415.

\bibitem{liuc2008finite} Liu, K.C.  and P.J. Cheng. 2008.
Finite Propagation of Heat Transfer in a Multilayer Tissue.
{\it Journal of Thermophysics and Heat Transfer\/}.  22(4):  775--782.

\bibitem{pennes1948analysis} Pennes, H.H. 1948.
Analysis of tissue and arterial blood temperatures in the resting human forearm.
{\it Journal of Applied Physiology\/}.  1(2): 93--122.

\bibitem{shih2007analytical} Shih, T.C., P. Yuan, W.L. LIn and H.S. Kou. 2007.
 Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface.
{\it Medical Engineering \& Physics\/}. 29(9): 946--953.


\bibitem{weinbaum1985new} Weinbaum, S.  and L.M. Jiji. 1985. 
A new simplified bioheat equation for the effect of blood flow on local average tissue temperature.
{\it ASME J. Biomech. Eng\/}. 107(2): 131--139.





\end{thebibliography}


\teacher{
\section*{NOTES FOR TEACHERS}

\subsection*{Prerequisites}
Before starting this activity, students should be able to
\begin{itemize}
\item have some knowledge of partial differential equations,
\item write Taylor's expansion, and
\item write code for numerical methods.
\end{itemize}

\subsection*{Objectives}
Upon successful completion of this activity, students will be able to
\begin{itemize}
\item solve the heat equation,
\item solve the problems with interface numerically, 
\item predict the behavior of temperature at the different depths of the skin, and  
\item read and explain behavior of the temperature using the graphs.
\end{itemize}

A more realistic model of skin is 3-dimensional. Students can extend the same technique to solve the 3-dimensional case. Other  similar  situations can be handled. This numerical technique is of lower order. Higher order accurate methods can also be obtained to predict the temperature distribution more accurately. 
}


\end{document}