1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2003-2018 by The University of Queensland
% http://www.uq.edu.au
%
% Primary Business: Queensland, Australia
% Licensed under the Apache License, version 2.0
% http://www.apache.org/licenses/LICENSE-2.0
%
% Development until 2012 by Earth Systems Science Computational Center (ESSCC)
% Development 2012-2013 by School of Earth Sciences
% Development from 2014 by Centre for Geoscience Computing (GeoComp)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Example 2: One Dimensional Heat Diffusion in an Iron Rod}
\sslist{example02.py}
\label{Sec:1DHDv0}
Our second example is of a cold iron bar at a constant temperature of
$T_{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
temperature at a constant level $T_0=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.
\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}
This problem is very similar to the example of temperature diffusion in granite
blocks presented in the previous Section~\ref{Sec:1DHDv00}. Thus, it is possible
to modify the script we have already developed for the granite blocks to suit
the iron bar problem.
The obvious differences 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_{0}$ at the left end of the rod where we have
$x_{0}=0$ and
$T_{ref}$ elsewhere. Instead of \verb|whereNegative| function we now
use \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 the iron rod model we want to keep the initial temperature $T_0$ on
the left side of the domain constant with time.
This implies that when we solve the PDE~\refEq{eqn:hddisc}, the solution must
have the value $T_0$ on the left hand side of the domain. As mentioned
already in Section~\ref{SEC BOUNDARY COND} where we discussed boundary
conditions, this kind of scenario can be expressed using 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 specify where to apply the
condition and determine 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 $T_{0}$ for the temperature at the left end of
the rod where $x_{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 $T_{0}$ 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 if a Dirichlet boundary condition is
prescribed on the same location as any Neumann boundary condition, the Neumann
boundary condition will be \textbf{overwritten}. This applies to Neumann
boundary conditions that \esc sets by default and those defined by the user.
Besides some cosmetic modification this is all we need to change. The total
energy over time is shown in \reffig{fig:onedheatout1 002}. As heat
is transferred into the rod by the heater the total energy is growing over time
but reaches a plateau when the temperature is constant in 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 model runs 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_{p} L_{0}^2}{4 \kappa}$, see
\url{http://en.wikipedia.org/wiki/Fick\%27s_laws_of_diffusion}}.
\begin{figure}[ht]
\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}[ht]
\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}
\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_{H}$.)
\item Change the boundary condition for the iron rod example such that the
temperature
at the right end is kept at a constant level $T_{ref}$, which
corresponds to the installation of a cooling element (Hint: Modify $q$ and
$r$).
\end{enumerate}
|
