% 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,empheq,mathrsfs}
\usepackage{xcolor,verbatim}
\usepackage{graphicx}
\usepackage{framed}
\usepackage{tasks}
\usepackage{exsheets}
\SetupExSheets[question]{type=exam}
\DeclareMathOperator{\csch}{csch}
%\usepackage{hyperref}
\usepackage{url}
%
\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 Chemical and biomedical applications of differential equations %\footnote{DELETE THIS FOOTNOTE BEFORE SUBMISSION.  In the \TeX \, source code there is a toggle/option to produce a Teacher Version or a Student Version by ``carving'' out with  $\backslash$teacher\{\} that material which should only appear in the Teacher version. See the source code.} 
\else STUDENT VERSION \fi\\Chemical and biomedical applications of differential equations}



\teacher{
    \keywords{Catalysis, second order equation, chemical engineering, cyclohexane, hydrogel
    }%Invasive species; fundamental solution; traveling waves; dispersion; migration}%Enter all words that could be considered key, to include words as general as differential equations or as specific as  kinetics. This list is to have no more than  five keywords. These will serve as tags for colleagues to use in searching SIMIODE for location of material. The Editor may suggest/change/add/delete some for more refinement as necessary.}

\tagwords{dehydrogenation, hyperbolic trig functions; porous medium, chemistry
}%Partial differential equation model; dispersive migration model; logistic growth; traveling wave solutions; fundamental solutions} %Here you can enter a wider set of words which serve as tags (as do keywords) in the SIMIODE web site. Keep this list to ten words. The Editor may suggest/change/add/delete some for more refinement as necessary.}
}

\markboth{Chemical and biomedical applications of ODEs}{Chemical and biomedical applications of ODEs}
%\markboth{Dehydrogenation of Cyclohexane}{From Dehydrogenation of Cyclohexane to tissue growth: applying ODEs in real life}

\author{Eric Stachura\\
Department of Mathematics\\ 
Kennesaw State University\\
Marietta, GA \teacher{30060} USA\\
\teacher{eric.stachura@kennesaw.edu}
\and
Tamara Lozano\\
Yokogawa Corporation\\
Newnan, GA \teacher{30265} USA\\
\teacher{tamara.lozano@us.yokogawa.com}
}

%}

\makeStitlePDFLaTex

\teacher{
\begin{abstract}
Students will be lead through a classical chemical engineering problem: to calculate the concentration profile of cyclohexane within a catalyst pellet by solving a second order linear differential equation. Then students will analyze the concentration as the radius of the catalyst shrinks to zero. Finally, they will explore a Mathematica Module where they can explore the relationships between the various parameters within the problem. In particular, they will see that decreasing the radius of the catalyst pellet results in a lower chemical concentration at the center of the pellet than at the surface. An optional set of questions related to a hydrogel model for knee joint replacements, using similar mathematics, is also included.
\end{abstract}
}

\teacher{
\vspace*{-.1 true cm}
\listkeywords

\vspace*{-.9 true cm}
\listtags
}

\teacher{
\section*{NOTES}
Solutions to the questions are in blue text beneath each question. Additionally, the last few exercises on the hydrogel are optional and left at the discretion of the instructor.
}

\section*{STATEMENT}

A \emph{catalyst} is a material which is capable of accelerating a chemical reaction without being consumed during the reaction process. Catalysts often consist of a porous material with high surface area, in which particles of a precious metal have been dispersed. Catalytic particles are often referred to as pellets. Surrounding each pellet, there is a thin gas film that contains a mix of reactants and products; see Figure \ref{pellet} and Figure \ref{catsurf1} for a 2D version.

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.2]{pellet7.pdf}
\caption{Catalyst pellet with a thin gas film surrounding it.}
\label{pellet}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.3]{recatantgraphic1neww}
\caption{2D slice of catalyst surface with reactants and products.}
\label{catsurf1}
\end{center}
\end{figure}


First, what does it mean to be porous?  Intuitively, such a medium has ``pores", or voids, which are typically filled with a fluid. Examples include sponges, rocks, soil, biological tissues, etc. You can imagine something as in Figure \ref{porousmedium}. 

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.1]{pore_scale_flow_structure}
\caption{A porous medium.}
\label{porousmedium}
\end{center}
\end{figure}


\teacher{For a more in-depth analysis of porous media and the differential equations involved in them, see \cite{ES19porosity}.} Then, the reactant must diffuse through the pellet pores to reach the metal atoms that are dispersed through the pellet. See Figure \ref{catsurf2}.

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.3]{reactantgraphic2neww}
\caption{Reactants must penetrate the pores of the surface in order for the reaction to occur.}
\label{catsurf2}
\end{center}
\end{figure}

In order for a reaction to occur, the reactant must reach the pellet surface after being transferred through the gas film. Similarly, after the reaction takes place and the product is formed, that product must transfer from the metal back to the pellet surface through the pore, and then out into the reaction medium through the gas film. See Figure \ref{catsurf3}.

 
 
 \begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.3]{reactantgraphic3neww}
\caption{Products must be transferred back to the surface through the pore.}
\label{catsurf3}
\end{center}
\end{figure}

 
 


An important problem in chemical engineering is to predict the diffusion and reaction in a porous catalyst pellet. The goal is to predict the overall reaction rate on the catalyst pellet surface \cite{Davis2013}.  \\









%\textcolor{red}{INSERT PICTURE}
%Thus, for a given chemical component, the sum of the masses entering the system and being generated by the chemical reaction must be equal to the sum of masses leaving the system and being consumed in the reaction. T
Conservation of mass is a scientific principle stating that in a chemical reaction, matter can neither be created nor destroyed. This can be expressed mathematically for $c$, the concentration of a given chemical, in a spherical pellet with radius $r_p$ by the following second order differential equation:
\begin{align} \label{conservationofmass1}
D \left( \dfrac{1}{r^2} \dfrac{d}{dr} \left( r^2 \dfrac{dc}{dr} \right) \right) = k R(c), \quad 0 < r< r_p 
\end{align}
where:
\begin{itemize}
\item $D$ is the diffusivity constant, in units of $cm^2/s$;
%\item $c$ is the concentration of a given chemical;
\item $k$ is the rate constant; 
\item $R(c)$ is the reaction rate function (could be nonlinear), in units of moles per volume per time (mol/L/s)
\end{itemize}
The boundary conditions are:
\[ \dfrac{dc}{dr} = 0 \quad \text{ at } \quad r=0\]
and
\[c=c_0 \quad \text{ at } \quad r=r_p \quad (\text{concentration is fixed at the surface})\]
Note that the units of $k$ actually depend on the order of the reaction! For example, in a first order reaction, the units of $k$ are $1/s$. %The units of $R(c)$ further depend on which function one chooses; for instance, if $R(c)=c$, the units are number of moles per $dm^3$ of open pore volume.

\begin{framed}
\begin{question}
Define the \emph{effectiveness factor} $E$ as the average reaction rate in the pellet divided by the average reaction rate on the surface:
\[ E \coloneqq \dfrac{\int_0^{r_p} R(c(r)) r^2 dr}{\int_0^{r_p} R(c_0) r^2 dr }.\]
Integrate the conservation of mass equation (\ref{conservationofmass1}) to obtain an expression for the effectiveness factor $E$ in terms of $r_p, D, \dfrac{dc}{dr}, k$, and $R(c_0)$. 
\end{question}
\end{framed}

%\newpage

\teacher{
\textcolor{blue}{Integrating the conservation of mass equation yields
\[ D \int_0^{r_p} \left( \dfrac{1}{r^2} \dfrac{d}{dr} \left( r^2 \dfrac{dc}{dr} \right)\right) r^2 dr = k \int_0^{r_p} R(c(r)) r^2 dr = D r_p^2 \left. \dfrac{dc}{dr}\right|_{r=r_p},\]
which in terms of $E$ yields
\[ E= \dfrac{3 r_p^2 D \left. \dfrac{dc}{dr} \right|_{r=r_p} }{k R(c_0)}. \]
}
}

\begin{framed} 
\begin{question}
Interpret what $E=1$ and $E<1$ would mean in terms of the reaction rate in the pellet. 
\end{question}
\end{framed}

\teacher{
\textcolor{blue}{If $E=1$ then the overall reaction rate in the pellet equals the surface value of the reaction rate, so mass transfer has no limiting effects on the reaction rate. If $E<1$ then mass transfer has limited the overall rate in the pellet.  
}
}


%\textcolor{red}{Explain below}
A very useful catalytic chemical reaction is the dehydrogenation of cyclohexane. See Figure \ref{Lewis} for the Lewis structure of cyclohexane (remember your chemistry class?). This industrial process requires the use of $\gamma-$alumina, a porous catalyst pellet with spherical shape (Google it--it's a real thing!). On this sphere, having diameter $5 \text{ mm}$, particles of platinum (a precious and very expensive metal) have been dispersed to catalyze chemical reaction.



\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.7]{Lewisstructure.png}
\caption{Lewis structure of cyclohexane $C_6H_{12}$.}
\label{Lewis}
\end{center}
\end{figure}

As the name suggests, dehydrogenation involves the removal of hydrogen atoms, a process that usually requires high temperatures. The catalyst $\gamma$-alumina is a popular choice of catalyst for this reaction due to its chemical properties, which make it resistant to the extreme reaction conditions. Figure \ref{process} provides a depiction of this chemical reaction process. It is a 3 stage process, and in each stage a hydrogen molecule is released. 

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.7]{reactionpic.png}
\caption{Dehydrogenation process.}
\label{process}
\end{center}
\end{figure}


Suppose that at $700\; \text{K}$, the rate constant for this reaction is $k=4 \; \text{s}^{-1}$ and the diffusivity $D=5\times 10^{-2} \text{cm}^2/s$. \textbf{Your goal is to calculate the concentration profile of cyclohexane within the pellet, as well as the effectiveness factor.}% $E$ for various reaction rate functions $R(c)$. }

The quantity of interest to be calculated is the concentration profile of cyclohexane, which we define as 
\[ C \coloneqq \dfrac{\text{concentration of cyclohexane inside the pellet}}{\text{concentration of cyclohexane at the surface of the pellet}}\]
Mass conservation for cyclohexane gives the following model (see \cite[Chapter 12]{Fogler2010} for a full derivation):
\begin{align} \label{cyclomodel}
 \dfrac{d^2 C}{d R^2} + \dfrac{2}{R} \dfrac{dC}{dR} = \Phi^2 \dfrac{R(c)}{c_0}, \quad 0<R<1
 \end{align}
where $R$ is the re-scaled radial coordinate, i.e $R=r/r_p$. The constant $\Phi$ is called the \emph{Thiele modulus} and is given by
\[ \Phi = r_p \sqrt{\dfrac{k}{D}}.\]
This quantity was introduced by E. W. Thiele in \cite{Thiele1939}, and came to describe the relationship between diffusion and reaction rates in porous catalyst pellets.\\
The boundary conditions are similar to the ones stated before, except scaled appropriately: 
\[\dfrac{dC}{dR} = 0 \quad \text{ at } \quad R=0 \]
and
\[ C=1 \quad \text{ at } \quad R=1 \quad (\text{this is by definition--we take $c_0=1$})\]


\begin{framed} 
\begin{question}
Suppose $R(c)=kc$, so you have a linear reaction rate function. Use the boundary conditions to solve and provide an explicit solution of (\ref{cyclomodel}) . HINT: first make the substitution $z= C R$ and solve the equation satisfied by $z$. 
\end{question}
\end{framed}
\teacher{
\textcolor{blue}{Making the substitution yields the equation
\[ \dfrac{d^2 z}{dR^2} = \Phi^2 z\]
whose solution is
\[ z(R)= c_1 e^{\Phi R} +c_2 e^{-\Phi R}\]
Now we apply the boundary conditions: the condition on the derivative of $C$ yields $c_1=-c_2$, and so 
\[z=c_1\left( e^{\Phi R} -e^{-\Phi R}\right) = 2c_1 \sinh(\phi R)\]
by definition of $\sinh(x)$. The second boundary condition yields
\[1=2c_1 \sinh(\Phi)\]
and hence
\[C=\dfrac{z}{R} = \dfrac{\sinh(\Phi R)}{R\sinh(\Phi)}.\]
}
}

\begin{framed} 
\begin{question}
Using your answer from the previous question, what is the effectiveness factor $E$ just in terms of the Thiele modulus $\Phi$? See Question 1.
\end{question}
\end{framed}

\teacher{
\textcolor{blue}{We integrate
\[C(R)= \dfrac{\sinh(\Phi R)}{R\sinh(\Phi)}\]
with respect to $R^2 dR$ from $R=0$ to $R=1$ (since we have rescaled), and use the definition of $E$ to obtain
\[E = \dfrac{3}{\Phi}\left( \dfrac{1}{\tanh(\Phi)} -\dfrac{1}{\Phi}\right)\]
See in particular \cite{Carberry2001}.
}
}

\begin{framed} 
\begin{question}
Analyze the limit of $C(R)$ as $R\to 0$. Interpret this in terms of the reaction in the center of the pellet. Further, with the value of $E$ that you computed in the last question, what can you say about the average reaction rate in the pellet?
\end{question}
\end{framed}

\teacher{
\textcolor{blue}{With the help of a computer, we find that
\[\lim\limits_{R\to 0} C(R)= \Phi \csch(\Phi)\]
A plot of this function can be found in Figure \ref{phiplot}. In particular, since $\Phi >0$, we have that $\lim\limits_{R\to 0}C(R) \neq 0$ and is finite, which means the concentration is finite as $R\to 0$.
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{phi_plot}
\caption{Plot of $\Phi \csch(\Phi)$. Note that only the right half of the graph is relevant here since $\Phi>0$.}
\label{phiplot}
\end{center}
\end{figure}
Moreover, by plotting $E$ as a function of $\Phi$, we see that $E<1$. Thus by Question 2, we see that mass transfer has limited the overall reaction rate in the pellet. 
}
}

\begin{framed} 
\begin{question}
Suppose you have a nonlinear reaction rate now, so that $R(c) = kc^2$. How might you go about computing the concentration profile of cyclohexane and calculating the effectiveness factor in this case? 
\end{question}
\end{framed}
\teacher{
\textcolor{blue}{The material balance equation is now
\[ \dfrac{d^2C}{dR^2} + \dfrac{2}{R} \dfrac{dC}{dR} = \Phi^2 C^2, \quad 0 < R< 1\]
This is a nonlinear second order equation! Our only hope is to solve this numerically. 
}
}
\begin{framed} 
\begin{question}
Explore the catalyst analysis module, adapted from \cite{Baumann2017}. You will see the relationship between the concentration of the reactant concentration inside the pellet and the pellet radius. What happens as you decrease the pellet radius? What happens if you decrease the rate constant? Why do you think that is?
\end{question}
\end{framed}
\teacher{
\textcolor{blue}{A decrease in the pellet radius permits the reactant to diffuse towards the center more quickly. Since the reactant is being consumed in the reaction, the concentration in the center of the particle is lower than in the surface ($C_{a,s}$). Thus, as the pellet radius decreases, a decrease in the reactant concentration is observed. 
On the other hand, the rate constant is directly proportional to the rate of the reaction, which describes how fast is the reactant being consumed. A high rate constant would imply fast consumption of reactant, therefore its concentration decreases faster. When the rate constant is low, the reaction takes place at a slower pace and the decrease in concentration occurs more smoothly. \\
}
}
\teacher{
%These last exercises are optional. \\

The equations describing diffusion and reaction in porous catalysts can also be used to derive rates of tissue growth. Over 200,000 patients each year receive knee joint replacements. One approach to growing cartilage to repair damaged knees is to deliver cartilage forming cells in a hydrogel at the damaged area. The design of the gel must be so that the gel can maintain the necessary rates of diffusion of nutrients in the hydrogel. In particular, the gel thickness needs to be designed in such a way so as to allow for rapid transport of oxygen \cite[Chapter 12]{Fogler2010}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{cylinderpic}
\caption{Hydrogel model.}
\label{default}
\end{center}
\end{figure}


\begin{framed} 
\begin{question}
Your goal is to find the gel thickness at which the minimum oxygen consumption rate is $10^{-13}$ mol/cell/h, where $h$ stands for hours. The cell density is $10^{10}$ cells/$\text{dm}^3$, the bulk concentration of oxygen is $2\times 10^{-4}$ mol/$\text{dm}^3$, and the diffusivity is $10^{-5} \text{cm}^2/s$. Physical considerations (again, see \cite[Chapter 12]{Fogler2010} for a full derivation of this result) imply that the equation satisfied by the dimensionless concentration $C$ satisfies
\[ \dfrac{d^2 C}{dz^2} =2\Phi_0, \quad \Phi_0 = \dfrac{k}{2 DC_0} L^2\]
with boundary conditions
\[ C(z)=1 \quad \text{ at } \quad z=0\]
and
\[\dfrac{dC}{dz} =0 \quad \text{ at } \quad z=1.\]
Above, $L$ denotes the thickness of the gel. 
Solve for $C(z)$. It turns out, your solution is only valid for $\Phi_0 \leq 1$. To see why, try setting $\Phi_0 = 10$ and then calculate $C(0.1)$. What happens?
\end{question}
\end{framed}
\textcolor{blue}{Integrating the second order equation once yields
\[\dfrac{dC}{dz} = 2\Phi_0 z + K_1\]
for some constant $K_1$. The second boundary condition implies that $K_1 = -2\Phi_0$. Integrate the resulting first order equation once more to get
\[C(z) = \Phi_0 z^2 -2\Phi_0 z + K_2\]
for some constant $K_2$. Using the other boundary condition, we get $K_2 = 1$, so the concentration profile is 
\begin{align} \label{cartilageprofile}
C(z) = \Phi_0 z(z-2) + 1
\end{align}
Note that the dimensionless concentration profile (\ref{cartilageprofile}) is only valid for values of $\Phi_0 \leq 1$. Indeed, you can see this by setting $\Phi_0=10$ and $z=0.1$ we find that $C(z) = -0.9$, a negative concentration profile! Thus, mathematically we should actually restrict ourselves to \emph{positive} solutions. Physically, negative concentration values are not useful.
}
\begin{framed} 
\begin{question}
With the given data from Question 8, calculate the rate constant $k$ (your units should be in mole/$\text{dm}^3 \cdot \text{h}$). Then, calculate the Thiele modulus $\Phi_0$. 
\end{question}
\end{framed}
\textcolor{blue}{We calculate
\[ k = \dfrac{10^{10} \text{cells}}{\text{dm}^3} \cdot \dfrac{10^{-13} \text{mole}\;  O_2}{\text{cell} \cdot \text{h}} = 10^{-3} \text{mole}/\text{dm}^3 \; \text{h}\]
and 
\[\Phi_0 = \dfrac{10^{-3} \text{mol}/\text{dm}^3 \; h }{2 \times 0.2 \times 10^{-3} \text{mol}/\text{dm}^3 \cdot 10^{-5} \text{cm}^2/s \times (3600 s)/h} \cdot L^2 = 70 \; L^2 \; \text{cm}^{-2}\]
Notice the conversion from seconds to hours in the last step. 
}
\begin{framed} 
\begin{question}
Suppose the gel is completely effective such that the concentration of oxygen is reduced to zero by the time the oxygen reaches the center of the gel, i.e. $C=0$ at $z=1$. What is the gel thickness in this case? 
\end{question}
\end{framed}
\textcolor{blue}{From (\ref{cartilageprofile}) we find that $\Phi_0=1$, so that now using the definition of $\Phi_0$:
\[ \Phi_0 =1 =70 \; \text{cm}^{-2}\; L^2\]
and hence $L=0.12$ cm. 
}
\begin{framed} 
\begin{question}
Next suppose the minimum oxygen concentration for the cells to survive is $0.1$ mm/$\text{dm}^3$, which is one half that at the surface, i.e. $C(1)=0.5$. What is the gel thickness in this case?
\end{question}
\end{framed}
\textcolor{blue}{ Again from (\ref{cartilageprofile}) we have that
\[\Phi_0 = 0.5 = 70 \; \text{cm}^{-2}\; L^2\]
which gives $L=0.085$ cm or $850$ microns. 
}
}


\begin{thebibliography}{9}

\bibitem{Davis2013}
M. E. Davis, \textit{Numerical Methods and modeling for chemical engineers}. Mineola, New York: Courier Corporation, 2013. 

\bibitem{Carberry2001}
J. J. Carberry, \textit{Chemical and catalytic reaction engineering}. Mineola, New York: Courier Corporation, 2001. 


\bibitem{Baumann2017}
R. L. Baumann, \textit{Diffusion and Reaction in a Catalyst Pellet}. Wolfram Demonstrations Project, \url{http://demonstrations.wolfram.com/DiffusionAndReactionInACatalystPellet/}. Published October 31, 2017.

\bibitem{Fogler2010}
H. Scott Fogler, \textit{Elements of Chemical Reaction Engineering}. Westford, Massachusetts: Pearson Education, 2010.  

\bibitem{Thiele1939}
E. W. Thiele, \textit{Relation between catalytic activity and size of particle}. Industrial and Engineering Chemistry, 31, 916--920, 1939.

\teacher{
\bibitem{ES19porosity}
E. Stachura (2019), ``9-010-S-TravelingWaves," \url{https://www.simiode.org/resources/6296}.
}
\end{thebibliography}


\teacher{
\section*{COMMENTS}

We recommend having students work on this problem after they have been exposed to second order equations. We expect this problem to be of interest for chemistry or chemical engineering students. Numerical analysis can be added on if desired, in particular in Question 6. Additionally, the last exercises are somewhat involved, so we have left it as optional, and won't include it in the student version. Students may need help in integrating hyperbolic trig functions, but we suggest letting them do this on the computer. There is a lot of information for students to parse through in this project. In particular, we have designed the questions in such a way so that students need to \emph{read with meaning} to find the useful data. We have also decided not to include a derivation of the conservation of mass differential equation. The reasoning for this is that the problem is already technical and rather involved. However, it would be an interesting addition to this project to have students derive the differential equation (see \cite{Fogler2010}). 

%Finally, we have included another problem on calculating gel thickness for knee joint replacements. This can be included at the discretion of the instructor. We have left it out though due to the already complicated nature of the problem students are solving in this activity. 


%\begin{thebibliography}{99}
%
%
%\bibitem{ES19porosity}
%E. Stachura (2019), ``9-010-S-TravelingWaves," \url{https://www.simiode.org/resources/6296}.
%
%
%\end{thebibliography}
%\bibitem{2012Kerckhove} Kerckhove, M. 2012. {\it From Population Dynamics to Partial Differential Equations}. The Mathematica Journal, 14, 18pp.
%
%\bibitem{1994Holmes} Holmes, E. E., Lewis, M. A., Banks, J. E., and Veit, R. R. 1994. {\it Partial Differential Equations in Ecology: Spatial  Interactions and Population Dynamics}. Ecology, 75 (1), pp. 17--29. 
%
%\bibitem{2014Roos} Roos, A. 2012. {\it Modeling Population Dynamics}. Course notes. 222pp. 

%\bibitem{Gause1934} Gause, G.~F. 1971. {\it The Struggle for Existence}. New York:  Dover Publications, Inc. First published in 1934 by The Williams \& Wilkins Company and available completely on the world wide web at \url{http://www.ggause.com/Contgau.htm}. Accessed  31 March 2008
%
%
%\bibitem{Hutchinson1978a} Hutchinson, G. Evelyn. 1978. {\it Introduction to Population Ecology\/}. New Haven CT:  Yale University Press.


}






\end{document}
