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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2979 - (show annotations)
Tue Mar 9 02:54:32 2010 UTC (9 years, 8 months ago) by ahallam
File MIME type: application/x-tex
File size: 6953 byte(s)
cookbook review final final 3.1 - artak and tony corrections
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 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/osl-3.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 \verb|whereNegative| function we use
61 now the
62 \verb|whereZero| 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*(1-whereZero(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 \verb|whereZero| 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 non-positive,
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

  ViewVC Help
Powered by ViewVC 1.1.26