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

Revision 2905 - (hide annotations)
Tue Feb 2 04:26:42 2010 UTC (9 years, 9 months ago) by gross
Original Path: trunk/doc/cookbook/onedheatdiff002.tex
File MIME type: application/x-tex
File size: 6675 byte(s)
some missing files.

 1 ahallam 2401 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 jfenwick 2881 % Copyright (c) 2003-2010 by University of Queensland 5 ahallam 2401 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 ahallam 2606 \begin{figure}[h!] 14 ahallam 2645 \centerline{\includegraphics[width=4.in]{figures/onedheatdiff002}} 15 gross 2905 \caption{One dimensional model of an Iron bar.} 16 \label{fig:onedhdmodel} 17 ahallam 2606 \end{figure} 18 19 gross 2905 \section{One Dimensional Heat Diffusion in an Iron Rod} 20 \sslist{onedheatdiff002.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 granit blocks presented in the previous section~\ref{Sec:1DHDv00}. So we will modify the script we have already developed for the granit 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 ahallam 2801 \begin{python} 30 gross 2905 #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 ahallam 2606 #PDE related 36 gross 2905 rho = 7874. *kg/m**3 #kg/m^{3} density of iron 37 cp = 449.*J/(kg*K) # J/Kg.K thermal capacity 38 rhocp = rho*cp 39 kappa = 80.*W/m/K # watts/m.Kthermal conductivity 40 qH = 0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 41 Tref = 20 * Celsius # base temperature of the rod 42 T0 = 100 * Celsius # temperature at heating element 43 tend= 0.5 * day # - time to end simulation 44 ahallam 2801 \end{python} 45 gross 2905 We also need to alter the initial value for the temperature. Now we need to set the 46 temperature to $T\hackscore{0}$ at the left end of the rod where we have $x\hackscore{0}=0$ and 47 $T\hackscore{ref}$ elsewhere. Instead of \verb|whereNegative| function we use now the 48 \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 temperature is set to; 51 \begin{python} 52 # ... set initial temperature .... 53 T= T0*whereZero(x[0])+Tref*(1-whereZero(x[0])) 54 \end{python} 55 ahallam 2401 56 gross 2905 \subsection{Dirchlet Boundary Conditions} 57 In iron rod model we want to keep the initial temperature $T\hackscore0$ on the left side of the domain over time. 58 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 boundary condition this kind of condition are called a \textbf{Dirichlet boundary condition}. Some people also 61 use the term \textbf{constraint} for the PDE. 62 63 To define a Dirichlet boundary condition we need to define where to apply the condition and what value the 64 solution should have at these locations. In \esc we use $q$ and $r$ to define the Dirichlet boundary conditions 65 for a PDE. The solution $u$ of the PDE is set to $r$ for all sample points where $q$ has a positive value. 66 Mathematically this is expressed in the form; 67 \begin{equation} 68 u(x) = r(x) \mbox{ for any } x \mbox{ with } q(x) > 0 69 \end{equation} 70 In the case of the iron rod 71 we can set; 72 ahallam 2801 \begin{python} 73 gross 2905 q=whereZero(x[0]) 74 r=T0 75 ahallam 2801 \end{python} 76 gross 2905 to prescibe 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 alread 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 ahallam 2401 82 gross 2905 To set the Dirichlet boundary conditions for the PDE to be solved in each time step we need 83 to add some statements; 84 ahallam 2801 \begin{python} 85 gross 2905 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 ahallam 2801 \end{python} 91 gross 2905 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 ahallam 2401 94 gross 2905 \begin{figure} 95 \begin{center} 96 \includegraphics[width=4in]{figures/ttrodpyplot150} 97 \caption{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{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 \section{For the Reader} 127 ahallam 2401 \begin{enumerate} 128 gross 2905 \item Move the boundary line between the two granite blocks to another part of the domain. 129 \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 ahallam 2401 \end{enumerate} 135