1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032012 by University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Open Software License version 3.0 
8 
% http://www.opensource.org/licenses/osl3.0.php 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development since 2012 by School of Earth Sciences 
12 
% 
13 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
14 

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

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

26 
\begin{figure}[ht] 
27 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}} 
28 
\caption{Example 2: One dimensional model of an Iron bar} 
29 
\label{fig:onedhdmodel} 
30 
\end{figure} 
31 

32 
This problem is very similar to the example of temperature diffusion in granite 
33 
blocks presented in the previous Section~\ref{Sec:1DHDv00}. Thus, it is possible 
34 
to modify the script we have already developed for the granite blocks to suit 
35 
the iron bar problem. 
36 
The obvious differences between the two problems are the dimensions of the 
37 
domain and different materials involved. This will change the time scale of the 
38 
model from years to hours. The new settings are 
39 
\begin{python} 
40 
#Domain related. 
41 
mx = 1*m #meters  model length 
42 
my = .1*m #meters  model width 
43 
ndx = 100 # mesh steps in x direction 
44 
ndy = 1 # mesh steps in y direction  one dimension means one element 
45 
#PDE related 
46 
rho = 7874. *kg/m**3 #kg/m^{3} density of iron 
47 
cp = 449.*J/(kg*K) # J/Kg.K thermal capacity 
48 
rhocp = rho*cp 
49 
kappa = 80.*W/m/K # watts/m.Kthermal conductivity 
50 
qH = 0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
51 
Tref = 20 * Celsius # base temperature of the rod 
52 
T0 = 100 * Celsius # temperature at heating element 
53 
tend= 0.5 * day #  time to end simulation 
54 
\end{python} 
55 
We also need to alter the initial value for the temperature. Now we need to set 
56 
the temperature to $T_{0}$ at the left end of the rod where we have 
57 
$x_{0}=0$ and 
58 
$T_{ref}$ elsewhere. Instead of \verbwhereNegative function we now 
59 
use \verbwhereZero which returns the value one for those sample points where 
60 
the argument (almost) equals zero and the value zero elsewhere. The initial 
61 
temperature is set to 
62 
\begin{python} 
63 
# ... set initial temperature .... 
64 
T= T0*whereZero(x[0])+Tref*(1whereZero(x[0])) 
65 
\end{python} 
66 

67 
\subsection{Dirichlet Boundary Conditions} 
68 
In the iron rod model we want to keep the initial temperature $T_0$ on 
69 
the left side of the domain constant with time. 
70 
This implies that when we solve the PDE~\refEq{eqn:hddisc}, the solution must 
71 
have the value $T_0$ on the left hand side of the domain. As mentioned 
72 
already in Section~\ref{SEC BOUNDARY COND} where we discussed boundary 
73 
conditions, this kind of scenario can be expressed using a 
74 
\textbf{Dirichlet boundary condition}. Some people also use the term 
75 
\textbf{constraint} for the PDE. 
76 

77 
To define a Dirichlet boundary condition we need to specify where to apply the 
78 
condition and determine what value the 
79 
solution should have at these locations. In \esc we use $q$ and $r$ to define 
80 
the Dirichlet boundary conditions for a PDE. The solution $u$ of the PDE is set 
81 
to $r$ for all sample points where $q$ has a positive value. 
82 
Mathematically this is expressed in the form; 
83 
\begin{equation} 
84 
u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0 
85 
\end{equation} 
86 
In the case of the iron rod we can set 
87 
\begin{python} 
88 
q=whereZero(x[0]) 
89 
r=T0 
90 
\end{python} 
91 
to prescribe the value $T_{0}$ for the temperature at the left end of 
92 
the rod where $x_{0}=0$. 
93 
Here we use the \verbwhereZero function again which we have already used to 
94 
set the initial value. 
95 
Notice that $r$ is set to the constant value $T_{0}$ for all sample 
96 
points. In fact, values of $r$ are used only where $q$ is positive. Where $q$ 
97 
is nonpositive, $r$ may have any value as these values are not used by the PDE 
98 
solver. 
99 

100 
To set the Dirichlet boundary conditions for the PDE to be solved in each time 
101 
step we need to add some statements; 
102 
\begin{python} 
103 
mypde=LinearPDE(rod) 
104 
A=zeros((2,2))) 
105 
A[0,0]=kappa 
106 
q=whereZero(x[0]) 
107 
mypde.setValue(A=A, D=rhocp/h, q=q, r=T0) 
108 
\end{python} 
109 
It is important to remark here that if a Dirichlet boundary condition is 
110 
prescribed on the same location as any Neumann boundary condition, the Neumann 
111 
boundary condition will be \textbf{overwritten}. This applies to Neumann 
112 
boundary conditions that \esc sets by default and those defined by the user. 
113 

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

128 
\begin{figure}[ht] 
129 
\begin{center} 
130 
\includegraphics[width=4in]{figures/ttrodpyplot150} 
131 
\caption{Example 2: Total Energy in the Iron Rod over Time (in seconds)} 
132 
\label{fig:onedheatout1 002} 
133 
\end{center} 
134 
\end{figure} 
135 

136 
\begin{figure}[ht] 
137 
\begin{center} 
138 
\includegraphics[width=4in]{figures/rodpyplot001} 
139 
\includegraphics[width=4in]{figures/rodpyplot050} 
140 
\includegraphics[width=4in]{figures/rodpyplot200} 
141 
\caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps 
142 
$1$, $50$ and $200$} 
143 
\label{fig:onedheatout 002} 
144 
\end{center} 
145 
\end{figure} 
146 

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