1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 

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

18 
\label{Sec:1DHDv0} 
19 
\begin{figure}[ht] 
20 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}} 
21 
\caption{Example 2: One dimensional model of an Iron bar.} 
22 
\label{fig:onedhdmodel} 
23 
\end{figure} 
24 

25 
Our second example is of a cold iron bar at a constant temperature of 
26 
$T\hackscore{ref}=20^{\circ} C$, see \reffig{fig:onedhdmodel}. The bar is 
27 
perfectly insulated on all sides with a heating element at one end keeping the 
28 
the temperature at a constant level $T\hackscore0=100^{\circ} C$. As heat is 
29 
applied; energy will disperse along the bar via conduction. With time the bar 
30 
will reach a constant temperature equivalent to that of the heat source. 
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 difference between the two problems are the dimensions of the domain 
37 
and different materials involved. This will change the time scale of the model 
38 
from years to hours. 
39 
The new settings are; 
40 
\begin{python} 
41 
#Domain related. 
42 
mx = 1*m #meters  model length 
43 
my = .1*m #meters  model width 
44 
ndx = 100 # mesh steps in x direction 
45 
ndy = 1 # mesh steps in y direction  one dimension means one element 
46 
#PDE related 
47 
rho = 7874. *kg/m**3 #kg/m^{3} density of iron 
48 
cp = 449.*J/(kg*K) # J/Kg.K thermal capacity 
49 
rhocp = rho*cp 
50 
kappa = 80.*W/m/K # watts/m.Kthermal conductivity 
51 
qH = 0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
52 
Tref = 20 * Celsius # base temperature of the rod 
53 
T0 = 100 * Celsius # temperature at heating element 
54 
tend= 0.5 * day #  time to end simulation 
55 
\end{python} 
56 
We also need to alter the initial value for the temperature. Now we need to set 
57 
the 
58 
temperature to $T\hackscore{0}$ at the left end of the rod where we have 
59 
$x\hackscore{0}=0$ and 
60 
$T\hackscore{ref}$ elsewhere. Instead of \verbwhereNegative function we use 
61 
now the 
62 
\verbwhereZero which returns the value one for those sample points where 
63 
the argument (almost) equals zero and the value zero elsewhere. The initial 
64 
temperature is set to; 
65 
\begin{python} 
66 
# ... set initial temperature .... 
67 
T= T0*whereZero(x[0])+Tref*(1whereZero(x[0])) 
68 
\end{python} 
69 

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

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

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

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

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

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

155 

156 

157 

158 

159 

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