doc/cookbook/example02.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2010 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht]
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}}
\caption{Example 2: One dimensional model of an Iron bar.}
\label{fig:onedhdmodel}
\end{figure}

\section{Example 2: One Dimensional Heat Diffusion in an Iron Rod}
\sslist{example02.py}
\label{Sec:1DHDv0}

Our second example is a cold iron bar at a constant temperature of $T\hackscore{ref}=20^{\circ} C$, see \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end keeping the the temperature at a constant level $T\hackscore0=100^{\circ} C$.  As heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source.

This problem is very similar to the example of temperature diffusion in granite blocks presented in the previous section~\ref{Sec:1DHDv00}. So we modify the script we have already developed for the granite blocks to adjust 
it to the iron bar problem.  
The obvious difference between the two problems are the dimensions of the domain and different materials involved. This will change the time scale of the model from years to hours. 
The new settings are;
\begin{python}
#Domain related.
mx = 1*m #meters - model length
my = .1*m #meters - model width
ndx = 100 # mesh steps in x direction 
ndy = 1 # mesh steps in y direction - one dimension means one element
#PDE related
rho = 7874. *kg/m**3 #kg/m^{3} density of iron
cp = 449.*J/(kg*K) # J/Kg.K thermal capacity
rhocp = rho*cp 
kappa = 80.*W/m/K   # watts/m.Kthermal conductivity
qH = 0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source
Tref = 20 * Celsius  # base temperature of the rod
T0 = 100 * Celsius # temperature at heating element
tend= 0.5 * day # - time to end simulation
\end{python}
We also need to alter the initial value for the temperature. Now we need to set the 
temperature to $T\hackscore{0}$ at the left end of the rod where we have $x\hackscore{0}=0$ and 
$T\hackscore{ref}$ elsewhere. Instead of \verb|whereNegative| function we use now the 
\verb|whereZero| which returns the value one for those sample points where 
the argument (almost) equals zero and the value zero elsewhere. The initial
temperature is set to;
\begin{python}
# ... set initial temperature ....
T= T0*whereZero(x[0])+Tref*(1-whereZero(x[0]))
\end{python}

\subsection{Dirichlet Boundary Conditions}
In iron rod model  we want to keep the initial temperature $T\hackscore0$ on the left side of the domain over time. 
So when we solve the PDE~\refEq{eqn:hddisc} the solution must have the value $T\hackscore0$ on the left hand
side of the domain. As mentioned already in Section~\ref{SEC BOUNDARY COND} where we discussed
boundary condition this kind of condition are called a \textbf{Dirichlet boundary condition}. Some people also
use the term \textbf{constraint} for the PDE. 

To define a Dirichlet boundary condition we need to define where to apply the condition and what value the 
solution should have at these locations. In \esc we use $q$ and $r$ to define the Dirichlet boundary conditions
for a PDE. The solution $u$ of the PDE is set to $r$ for all sample points where $q$ has a positive value.
Mathematically this is expressed in the form;
\begin{equation}
  u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0
\end{equation} 
In the case of the iron rod 
we can set;
\begin{python}
q=whereZero(x[0])
r=T0
\end{python}
to prescribe the value $T0$ for the temperature at the left end of the rod where $x\hackscore{0}=0$. 
Here we use the \verb|whereZero| function again which we have already used to set the initial value.
Notice that $r$ is set to the constant value $r$ for all sample points. In fact, 
values of $r$ are used only where $q$ is positive. Where $q$ is non-positive,
$r$ may have any value as these values are not used by the PDE solver. 

To set the Dirichlet boundary conditions for the PDE to be solved in each time step we need
to add some statements;
\begin{python}
mypde=LinearPDE(rod)
A=zeros((2,2)))
A[0,0]=kappa
q=whereZero(x[0])
mypde.setValue(A=A, D=rhocp/h, q=q, r=T0)
\end{python}
It is important to remark here that the Dirichlet condition \textbf{overwrites} any Neuman boundary 
condition \esc sets by default (or you may set).  

\begin{figure}
\begin{center}
\includegraphics[width=4in]{figures/ttrodpyplot150}
\caption{Example 2: Total Energy in the Iron Rod over Time (in seconds).}
\label{fig:onedheatout1 002} 
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=4in]{figures/rodpyplot001}
\includegraphics[width=4in]{figures/rodpyplot050}
\includegraphics[width=4in]{figures/rodpyplot200}
\caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps $1$, $50$ and $200$.}
\label{fig:onedheatout 002} 
\end{center}
\end{figure}

Besides some cosmetic modification this all we need to change. The total energy over time is shown in \reffig{fig:onedheatout1 002}. As heat
is transfered into the rod by the heater the total energy is growing over time but reaches a plateau 
when the temperature is constant is the rod, see \reffig{fig:onedheatout 002}. 
You will notice that the time scale of this model is several order of magnitudes faster than
for the granite rock problem due to the different length scale and material parameters. 
In practice it can take a few models run before the right time scale has been chosen\footnote{An estimate of the
time scale for a diffusion problem is given by the formula $\frac{\rho c\hackscore{p} L\hackscore{0}^2}{4 \kappa}$, see
\url{http://en.wikipedia.org/wiki/Fick\%27s_laws_of_diffusion}}.


\section{For the Reader}
\begin{enumerate}
 \item Move the boundary line between the two granite blocks to another part of the domain.
 \item Split the domain into multiple granite blocks with varying temperatures.
 \item Vary the mesh step size. Do you see a difference in the answers? What does happen with the compute time?
 \item Insert an internal heat source (Hint: The internal heat source is given by $q\hackscore{H}$.)
 \item Change the boundary condition for iron rod example such that the temperature 
 at the right end is kept at a constant level $T\hackscore{ref}$, which corresponds to the installation of a cooling element (Hint: Modify $q$ and $r$). 
\end{enumerate}

1
2	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3	%
4	% Copyright (c) 2003-2010 by University of Queensland
5	% Earth Systems Science Computational Center (ESSCC)
6	% http://www.uq.edu.au/esscc
7	%
8	% Primary Business: Queensland, Australia
9	% Licensed under the Open Software License version 3.0
10	% http://www.opensource.org/licenses/osl-3.0.php
11	%
12	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13	\begin{figure}[ht]
14	\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}}
15	\caption{Example 2: One dimensional model of an Iron bar.}
16	\label{fig:onedhdmodel}
17	\end{figure}
18
19	\section{Example 2: One Dimensional Heat Diffusion in an Iron Rod}
20	\sslist{example02.py}
21	\label{Sec:1DHDv0}
22
23	Our second example is a cold iron bar at a constant temperature of $T\hackscore{ref}=20^{\circ} C$, see \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end keeping the the temperature at a constant level $T\hackscore0=100^{\circ} C$. As heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source.
24
25	This problem is very similar to the example of temperature diffusion in granite blocks presented in the previous section~\ref{Sec:1DHDv00}. So we modify the script we have already developed for the granite blocks to adjust
26	it to the iron bar problem.
27	The obvious difference between the two problems are the dimensions of the domain and different materials involved. This will change the time scale of the model from years to hours.
28	The new settings are;
29	\begin{python}
30	#Domain related.
31	mx = 1*m #meters - model length
32	my = .1*m #meters - model width
33	ndx = 100 # mesh steps in x direction
34	ndy = 1 # mesh steps in y direction - one dimension means one element
35	#PDE related
36	rho = 7874. kg/m*3 #kg/m^{3} density of iron
37	cp = 449.J/(kgK) # J/Kg.K thermal capacity
38	rhocp = rho*cp
39	kappa = 80.*W/m/K # watts/m.Kthermal conductivity
40	qH = 0 * J/(secm*3) # J/(sec.m^{3}) no heat source
41	Tref = 20 * Celsius # base temperature of the rod
42	T0 = 100 * Celsius # temperature at heating element
43	tend= 0.5 * day # - time to end simulation
44	\end{python}
45	We also need to alter the initial value for the temperature. Now we need to set the
46	temperature to $T\hackscore{0}$ at the left end of the rod where we have $x\hackscore{0}=0$ and
47	$T\hackscore{ref}$ elsewhere. Instead of \verb\|whereNegative\| function we use now the
48	\verb\|whereZero\| which returns the value one for those sample points where
49	the argument (almost) equals zero and the value zero elsewhere. The initial
50	temperature is set to;
51	\begin{python}
52	# ... set initial temperature ....
53	T= T0whereZero(x[0])+Tref(1-whereZero(x[0]))
54	\end{python}
55
56	\subsection{Dirichlet Boundary Conditions}
57	In iron rod model we want to keep the initial temperature $T\hackscore0$ on the left side of the domain over time.
58	So when we solve the PDE~\refEq{eqn:hddisc} the solution must have the value $T\hackscore0$ on the left hand
59	side of the domain. As mentioned already in Section~\ref{SEC BOUNDARY COND} where we discussed
60	boundary condition this kind of condition are called a \textbf{Dirichlet boundary condition}. Some people also
61	use the term \textbf{constraint} for the PDE.
62
63	To define a Dirichlet boundary condition we need to define where to apply the condition and what value the
64	solution should have at these locations. In \esc we use $q$ and $r$ to define the Dirichlet boundary conditions
65	for a PDE. The solution $u$ of the PDE is set to $r$ for all sample points where $q$ has a positive value.
66	Mathematically this is expressed in the form;
67	\begin{equation}
68	u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0
69	\end{equation}
70	In the case of the iron rod
71	we can set;
72	\begin{python}
73	q=whereZero(x[0])
74	r=T0
75	\end{python}
76	to prescribe the value $T0$ for the temperature at the left end of the rod where $x\hackscore{0}=0$.
77	Here we use the \verb\|whereZero\| function again which we have already used to set the initial value.
78	Notice that $r$ is set to the constant value $r$ for all sample points. In fact,
79	values of $r$ are used only where $q$ is positive. Where $q$ is non-positive,
80	$r$ may have any value as these values are not used by the PDE solver.
81
82	To set the Dirichlet boundary conditions for the PDE to be solved in each time step we need
83	to add some statements;
84	\begin{python}
85	mypde=LinearPDE(rod)
86	A=zeros((2,2)))
87	A[0,0]=kappa
88	q=whereZero(x[0])
89	mypde.setValue(A=A, D=rhocp/h, q=q, r=T0)
90	\end{python}
91	It is important to remark here that the Dirichlet condition \textbf{overwrites} any Neuman boundary
92	condition \esc sets by default (or you may set).
93
94	\begin{figure}
95	\begin{center}
96	\includegraphics[width=4in]{figures/ttrodpyplot150}
97	\caption{Example 2: Total Energy in the Iron Rod over Time (in seconds).}
98	\label{fig:onedheatout1 002}
99	\end{center}
100	\end{figure}
101
102	\begin{figure}
103	\begin{center}
104	\includegraphics[width=4in]{figures/rodpyplot001}
105	\includegraphics[width=4in]{figures/rodpyplot050}
106	\includegraphics[width=4in]{figures/rodpyplot200}
107	\caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps $1$, $50$ and $200$.}
108	\label{fig:onedheatout 002}
109	\end{center}
110	\end{figure}
111
112	Besides some cosmetic modification this all we need to change. The total energy over time is shown in \reffig{fig:onedheatout1 002}. As heat
113	is transfered into the rod by the heater the total energy is growing over time but reaches a plateau
114	when the temperature is constant is the rod, see \reffig{fig:onedheatout 002}.
115	You will notice that the time scale of this model is several order of magnitudes faster than
116	for the granite rock problem due to the different length scale and material parameters.
117	In practice it can take a few models run before the right time scale has been chosen\footnote{An estimate of the
118	time scale for a diffusion problem is given by the formula $\frac{\rho c\hackscore{p} L\hackscore{0}^2}{4 \kappa}$, see
119	\url{http://en.wikipedia.org/wiki/Fick\%27s_laws_of_diffusion}}.
120
121
122
123
124
125
126	\section{For the Reader}
127	\begin{enumerate}
128	\item Move the boundary line between the two granite blocks to another part of the domain.
129	\item Split the domain into multiple granite blocks with varying temperatures.
130	\item Vary the mesh step size. Do you see a difference in the answers? What does happen with the compute time?
131	\item Insert an internal heat source (Hint: The internal heat source is given by $q\hackscore{H}$.)
132	\item Change the boundary condition for iron rod example such that the temperature
133	at the right end is kept at a constant level $T\hackscore{ref}$, which corresponds to the installation of a cooling element (Hint: Modify $q$ and $r$).
134	\end{enumerate}
135