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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 4 months ago) by jfenwick
File MIME type: application/x-tex
File size: 6993 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ahallam 2401
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 6651 % Copyright (c) 2003-2018 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/LICENSE-2.0
9 ahallam 2401 %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 % Development 2012-2013 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 \verb|whereNegative| function we now
60 caltinay 2982 use \verb|whereZero| 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*(1-whereZero(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 \verb|whereZero| 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 non-positive, $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    

  ViewVC Help
Powered by ViewVC 1.1.26