1 


2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
% 
4 
% Copyright (c) 20032009 by University of Queensland 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
% http://www.uq.edu.au/esscc 
7 
% 
% 
10 
% http://www.opensource.org/licenses/osl3.0.php 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

\begin{figure}[ht] 
14 
\section{One Dimensional Heat Diffusion accross an Interface} 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}} 
15 
%\label{Sec:1DHDv1} 
\caption{Example 2: One dimensional model of an Iron bar.} 
16 
It is quite simple to now expand upon the 1D heat diffusion problem we just tackled. Suppose we have two blocks of isotropic material which are very large in all directions to the point that the interface between the two blocks appears infinite in length compared to the distance we are modelling perpendicular to the interface and accross the two blocks. If \textit{Block 1} is of a temperature \verb 0 and \textit{Block 2} is at a temperature \verb T what would happen to the temperature distribution in each block if we placed them next to each other. This problem is very similar to our Iron Rod but instead of a constant heat source we instead have a heat disparity with a fixed amount of energy. In such a situation it is common knowledge that the heat energy in the warmer block will gradually conduct into the cooler block until the temperature between the blocks is balanced. 
\label{fig:onedhdmodel} 




\begin{figure}[h!] 


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


\caption{Temperature differential along a single interface between two granite blocks.} 


\label{fig:onedgbmodel} 

17 
\end{figure} 
\end{figure} 
18 


19 
By modifying our previous code it is possible to solve this new problem. In doing so we will also try to tackle a real world example and as a result, introduce and discuss some new variables. The linear model of the two blocks is very similar to the effect a large magmatic intrusion would have on a cold country rock. It is however, simpler at this stage to have both materials the same and for this example we will use granite \reffig{fig:onedgbmodel}. The intrusion will have an initial temperature defined by \verb Tref and the granite properties required are: 
\section{Example 2: One Dimensional Heat Diffusion in an Iron Rod} 
20 
\begin{verbatim} 
\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 will 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 
#PDE related 
36 
mx = 500*m #meters  model length 
rho = 7874. *kg/m**3 #kg/m^{3} density of iron 
37 
my = 100*m #meters  model width 
cp = 449.*J/(kg*K) # J/Kg.K thermal capacity 
38 
ndx = 500 # mesh steps in x direction 
rhocp = rho*cp 
39 
ndy = 1 # mesh steps in y direction 
kappa = 80.*W/m/K # watts/m.Kthermal conductivity 
40 
boundloc = mx/2 # location of boundary between two blocks 
qH = 0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
41 
q=0.*Celsius #our heat source temperature is now zero 
Tref = 20 * Celsius # base temperature of the rod 
42 
Tref=2273.*Celsius # Kelvin the starting temperature of our RHS Block 
T0 = 100 * Celsius # temperature at heating element 
43 
rho = 2750*kg/m**3 #kg/m^{3} density of granite 
tend= 0.5 * day #  time to end simulation 
44 
cp = 790.*J/(kg*K) #j/Kg.K thermal capacity 
\end{python} 
45 
rhocp = rho*cp #DENSITY * SPECIFIC HEAT 
We also need to alter the initial value for the temperature. Now we need to set the 
46 
eta=0. # RADIATION CONDITION 
temperature to $T\hackscore{0}$ at the left end of the rod where we have $x\hackscore{0}=0$ and 
47 
kappa=2.2*W/m/K #watts/m.K thermal conductivity 
$T\hackscore{ref}$ elsewhere. Instead of \verbwhereNegative function we use now the 
48 
\end{verbatim} 
\verbwhereZero which returns the value one for those sample points where 
49 

the argument (almost) equals zero and the value zero elsewhere. The initial 
50 
Since the scale and values involved in our problem have changed, the length and step size of the iteration must be considered. Instead of seconds which our units are in, it may be more prudent to decide the number of days or years we would like to run the simulation over. These can then be converted accordingly to SI units \editor{lutz new schema in here}. 
temperature is set to; 
51 
\begin{verbatim} 
\begin{python} 
52 
#Script/Iteration Related 
# ... set initial temperature .... 
53 
t=0. #our start time, usually zero 
T= T0*whereZero(x[0])+Tref*(1whereZero(x[0])) 
54 
tday=10*365. #the time we want to end the simulation in days 
\end{python} 
55 
tend=tday*24*60*60 

56 
outputs = 400 # number of time steps required. 
\subsection{Dirichlet Boundary Conditions} 
57 
h=(tendt)/outputs #size of time step 
In iron rod model we want to keep the initial temperature $T\hackscore0$ on the left side of the domain over time. 
58 
\end{verbatim} 
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 
If we assume that the dimensions of the blocks are continuous and extremely large compared with the model size then we need only model a small proportion of the boundary. It is practical to locate the boundary between the two blocks at the center of our model. As there is no heat source our \verb q variable can be set to zero. The new initial conditions are defined using the following: 
boundary condition this kind of condition are called a \textbf{Dirichlet boundary condition}. Some people also 
61 
\begin{verbatim} 
use the term \textbf{constraint} for the PDE. 
62 
#establish location of boundary between two blocks 

63 
bound = x[0]boundloc 
To define a Dirichlet boundary condition we need to define where to apply the condition and what value the 
64 
#set initial temperature 
solution should have at these locations. In \esc we use $q$ and $r$ to define the Dirichlet boundary conditions 
65 
T= 0*Tref*whereNegative(bound)+Tref*wherePositive(bound) 
for a PDE. The solution $u$ of the PDE is set to $r$ for all sample points where $q$ has a positive value. 
66 
\end{verbatim} 
Mathematically this is expressed in the form; 
67 
The \verb bound statement sets the boundary to the location along the \textit{xaxis} defined by \verb boundloc . 
\begin{equation} 
68 
The PDE can then be solved as before. 
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 \verbwhereZero 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 nonpositive, 
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 
FOR THE READER: 
\section{For the Reader} 
127 
\begin{enumerate} 
\begin{enumerate} 
128 
\item Move the boundary line between the two blocks to another part of the domain. 
\item Move the boundary line between the two granite blocks to another part of the domain. 
129 
\item Try splitting the domain in to multiple blocks with varying temperatures. 
\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} 
\end{enumerate} 
135 

