/[escript]/trunk/doc/cookbook/example02.tex
ViewVC logotype

Contents of /trunk/doc/cookbook/example02.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3989 - (show annotations)
Tue Sep 25 02:21:54 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: application/x-tex
File size: 6931 byte(s)
More copyright fixes.
pyvisi traces removed.
Some install doco
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2012 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/osl-3.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 \verb|whereNegative| function we now
59 use \verb|whereZero| 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*(1-whereZero(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 \verb|whereZero| 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 non-positive, $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

  ViewVC Help
Powered by ViewVC 1.1.26