1 
ahallam 
2401 

2 
jfenwick 
3989 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
jfenwick 
6651 
% Copyright (c) 20032018 by The University of Queensland 
4 
jfenwick 
3989 
% http://www.uq.edu.au 
5 
ahallam 
2401 
% 
6 


% Primary Business: Queensland, Australia 
7 
jfenwick 
6112 
% Licensed under the Apache License, version 2.0 
8 


% http://www.apache.org/licenses/LICENSE2.0 
9 
ahallam 
2401 
% 
10 
jfenwick 
3989 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
jfenwick 
4657 
% Development 20122013 by School of Earth Sciences 
12 


% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
jfenwick 
3989 
% 
14 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 
ahallam 
2975 

16 


\section{Example 2: One Dimensional Heat Diffusion in an Iron Rod} 
17 


\sslist{example02.py} 
18 


\label{Sec:1DHDv0} 
19 
ahallam 
2606 

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

27 
ahallam 
3370 
\begin{figure}[ht] 
28 


\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}} 
29 


\caption{Example 2: One dimensional model of an Iron bar} 
30 


\label{fig:onedhdmodel} 
31 


\end{figure} 
32 



33 
ahallam 
2979 
This problem is very similar to the example of temperature diffusion in granite 
34 
caltinay 
2982 
blocks presented in the previous Section~\ref{Sec:1DHDv00}. Thus, it is possible 
35 
ahallam 
2979 
to modify the script we have already developed for the granite blocks to suit 
36 


the iron bar problem. 
37 
caltinay 
2982 
The obvious differences between the two problems are the dimensions of the 
38 


domain and different materials involved. This will change the time scale of the 
39 


model from years to hours. The new settings are 
40 
ahallam 
2801 
\begin{python} 
41 
gross 
2905 
#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 
ahallam 
2606 
#PDE related 
47 
gross 
2905 
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 
ahallam 
2801 
\end{python} 
56 
ahallam 
2979 
We also need to alter the initial value for the temperature. Now we need to set 
57 
jfenwick 
3308 
the temperature to $T_{0}$ at the left end of the rod where we have 
58 


$x_{0}=0$ and 
59 


$T_{ref}$ elsewhere. Instead of \verbwhereNegative function we now 
60 
caltinay 
2982 
use \verbwhereZero which returns the value one for those sample points where 
61 
gross 
2905 
the argument (almost) equals zero and the value zero elsewhere. The initial 
62 
caltinay 
2982 
temperature is set to 
63 
gross 
2905 
\begin{python} 
64 


# ... set initial temperature .... 
65 


T= T0*whereZero(x[0])+Tref*(1whereZero(x[0])) 
66 


\end{python} 
67 
ahallam 
2401 

68 
gross 
2953 
\subsection{Dirichlet Boundary Conditions} 
69 
jfenwick 
3308 
In the iron rod model we want to keep the initial temperature $T_0$ on 
70 
ahallam 
2979 
the left side of the domain constant with time. 
71 


This implies that when we solve the PDE~\refEq{eqn:hddisc}, the solution must 
72 
jfenwick 
3308 
have the value $T_0$ on the left hand side of the domain. As mentioned 
73 
caltinay 
2982 
already in Section~\ref{SEC BOUNDARY COND} where we discussed boundary 
74 


conditions, this kind of scenario can be expressed using a 
75 


\textbf{Dirichlet boundary condition}. Some people also use the term 
76 


\textbf{constraint} for the PDE. 
77 
gross 
2905 

78 
ahallam 
2979 
To define a Dirichlet boundary condition we need to specify where to apply the 
79 
caltinay 
2982 
condition and determine what value the 
80 
ahallam 
2979 
solution should have at these locations. In \esc we use $q$ and $r$ to define 
81 
caltinay 
2982 
the Dirichlet boundary conditions for a PDE. The solution $u$ of the PDE is set 
82 


to $r$ for all sample points where $q$ has a positive value. 
83 
gross 
2905 
Mathematically this is expressed in the form; 
84 


\begin{equation} 
85 


u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0 
86 


\end{equation} 
87 
caltinay 
2982 
In the case of the iron rod we can set 
88 
ahallam 
2801 
\begin{python} 
89 
gross 
2905 
q=whereZero(x[0]) 
90 


r=T0 
91 
ahallam 
2801 
\end{python} 
92 
jfenwick 
3308 
to prescribe the value $T_{0}$ for the temperature at the left end of 
93 


the rod where $x_{0}=0$. 
94 
ahallam 
2979 
Here we use the \verbwhereZero function again which we have already used to 
95 


set the initial value. 
96 
jfenwick 
3308 
Notice that $r$ is set to the constant value $T_{0}$ for all sample 
97 
caltinay 
2982 
points. In fact, values of $r$ are used only where $q$ is positive. Where $q$ 
98 


is nonpositive, $r$ may have any value as these values are not used by the PDE 
99 


solver. 
100 
ahallam 
2401 

101 
ahallam 
2979 
To set the Dirichlet boundary conditions for the PDE to be solved in each time 
102 
caltinay 
2982 
step we need to add some statements; 
103 
ahallam 
2801 
\begin{python} 
104 
gross 
2905 
mypde=LinearPDE(rod) 
105 


A=zeros((2,2))) 
106 


A[0,0]=kappa 
107 


q=whereZero(x[0]) 
108 


mypde.setValue(A=A, D=rhocp/h, q=q, r=T0) 
109 
ahallam 
2801 
\end{python} 
110 
ahallam 
2979 
It is important to remark here that if a Dirichlet boundary condition is 
111 
caltinay 
2982 
prescribed on the same location as any Neumann boundary condition, the Neumann 
112 


boundary condition will be \textbf{overwritten}. This applies to Neumann 
113 
ahallam 
2979 
boundary conditions that \esc sets by default and those defined by the user. 
114 
ahallam 
2401 

115 
ahallam 
3370 
Besides some cosmetic modification this is all we need to change. The total 
116 


energy over time is shown in \reffig{fig:onedheatout1 002}. As heat 
117 


is transferred into the rod by the heater the total energy is growing over time 
118 


but reaches a plateau when the temperature is constant in the rod, see 
119 


\reffig{fig:onedheatout 002}. 
120 


You will notice that the time scale of this model is several order of 
121 


magnitudes faster than for the granite rock problem due to the different length 
122 


scale and material parameters. 
123 


In practice it can take a few model runs before the right time scale has been 
124 


chosen\footnote{An estimate of the 
125 


time scale for a diffusion problem is given by the formula $\frac{\rho 
126 


c_{p} L_{0}^2}{4 \kappa}$, see 
127 


\url{http://en.wikipedia.org/wiki/Fick\%27s_laws_of_diffusion}}. 
128 



129 


\begin{figure}[ht] 
130 
gross 
2905 
\begin{center} 
131 


\includegraphics[width=4in]{figures/ttrodpyplot150} 
132 
caltinay 
2982 
\caption{Example 2: Total Energy in the Iron Rod over Time (in seconds)} 
133 
gross 
2905 
\label{fig:onedheatout1 002} 
134 


\end{center} 
135 


\end{figure} 
136 



137 
ahallam 
3370 
\begin{figure}[ht] 
138 
gross 
2905 
\begin{center} 
139 


\includegraphics[width=4in]{figures/rodpyplot001} 
140 


\includegraphics[width=4in]{figures/rodpyplot050} 
141 


\includegraphics[width=4in]{figures/rodpyplot200} 
142 
ahallam 
2979 
\caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps 
143 
caltinay 
2982 
$1$, $50$ and $200$} 
144 
gross 
2905 
\label{fig:onedheatout 002} 
145 


\end{center} 
146 


\end{figure} 
147 



148 


\section{For the Reader} 
149 
ahallam 
2401 
\begin{enumerate} 
150 
ahallam 
2979 
\item Move the boundary line between the two granite blocks to another part of 
151 


the domain. 
152 
gross 
2905 
\item Split the domain into multiple granite blocks with varying temperatures. 
153 
ahallam 
2979 
\item Vary the mesh step size. Do you see a difference in the answers? What 
154 


does happen with the compute time? 
155 


\item Insert an internal heat source (Hint: The internal heat source is given 
156 
jfenwick 
3308 
by $q_{H}$.) 
157 
caltinay 
2982 
\item Change the boundary condition for the iron rod example such that the 
158 
ahallam 
2979 
temperature 
159 
jfenwick 
3308 
at the right end is kept at a constant level $T_{ref}$, which 
160 
ahallam 
2979 
corresponds to the installation of a cooling element (Hint: Modify $q$ and 
161 


$r$). 
162 
ahallam 
2401 
\end{enumerate} 
163 


