# Diff of /trunk/doc/cookbook/example02.tex

trunk/doc/cookbook/TEXT/onedheatdiff002.tex revision 2620 by gross, Thu Aug 20 06:24:00 2009 UTC trunk/doc/cookbook/example02.tex revision 2953 by gross, Thu Feb 25 10:00:23 2010 UTC
# Line 1  Line 1
1
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2009 by University of Queensland  % Copyright (c) 2003-2010 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  %  %
# Line 10  Line 10
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 \verb|whereNegative| function we use now the
48  \end{verbatim}  \verb|whereZero| 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*(1-whereZero(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=(tend-t)/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{x-axis} defined by \verb boundloc  .
68  The PDE can then be solved as before.    u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0
69
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 \verb|whereZero| 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 non-positive,
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

Legend:
 Removed from v.2620 changed lines Added in v.2953