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

Revision 3370 - (hide annotations)
Sun Nov 21 23:22:25 2010 UTC (9 years ago) by ahallam
File MIME type: application/x-tex
File size: 6815 byte(s)
Rearranged figures for release.

 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 2975 14 15 \section{Example 2: One Dimensional Heat Diffusion in an Iron Rod} 16 \sslist{example02.py} 17 \label{Sec:1DHDv0} 18 ahallam 2606 19 ahallam 2979 Our second example is of a cold iron bar at a constant temperature of 20 jfenwick 3308 $T_{ref}=20^{\circ} C$, see \reffig{fig:onedhdmodel}. The bar is 21 ahallam 2979 perfectly insulated on all sides with a heating element at one end keeping the 22 jfenwick 3308 temperature at a constant level $T_0=100^{\circ} C$. As heat is 23 caltinay 2982 applied energy will disperse along the bar via conduction. With time the bar 24 ahallam 2979 will reach a constant temperature equivalent to that of the heat source. 25 gross 2905 26 ahallam 3370 \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 ahallam 2979 This problem is very similar to the example of temperature diffusion in granite 33 caltinay 2982 blocks presented in the previous Section~\ref{Sec:1DHDv00}. Thus, it is possible 34 ahallam 2979 to modify the script we have already developed for the granite blocks to suit 35 the iron bar problem. 36 caltinay 2982 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 ahallam 2801 \begin{python} 40 gross 2905 #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 ahallam 2606 #PDE related 46 gross 2905 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 ahallam 2801 \end{python} 55 ahallam 2979 We also need to alter the initial value for the temperature. Now we need to set 56 jfenwick 3308 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 caltinay 2982 use \verb|whereZero| which returns the value one for those sample points where 60 gross 2905 the argument (almost) equals zero and the value zero elsewhere. The initial 61 caltinay 2982 temperature is set to 62 gross 2905 \begin{python} 63 # ... set initial temperature .... 64 T= T0*whereZero(x[0])+Tref*(1-whereZero(x[0])) 65 \end{python} 66 ahallam 2401 67 gross 2953 \subsection{Dirichlet Boundary Conditions} 68 jfenwick 3308 In the iron rod model we want to keep the initial temperature $T_0$ on 69 ahallam 2979 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 jfenwick 3308 have the value $T_0$ on the left hand side of the domain. As mentioned 72 caltinay 2982 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 gross 2905 77 ahallam 2979 To define a Dirichlet boundary condition we need to specify where to apply the 78 caltinay 2982 condition and determine what value the 79 ahallam 2979 solution should have at these locations. In \esc we use $q$ and $r$ to define 80 caltinay 2982 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 gross 2905 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 caltinay 2982 In the case of the iron rod we can set 87 ahallam 2801 \begin{python} 88 gross 2905 q=whereZero(x[0]) 89 r=T0 90 ahallam 2801 \end{python} 91 jfenwick 3308 to prescribe the value $T_{0}$ for the temperature at the left end of 92 the rod where $x_{0}=0$. 93 ahallam 2979 Here we use the \verb|whereZero| function again which we have already used to 94 set the initial value. 95 jfenwick 3308 Notice that $r$ is set to the constant value $T_{0}$ for all sample 96 caltinay 2982 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 ahallam 2401 100 ahallam 2979 To set the Dirichlet boundary conditions for the PDE to be solved in each time 101 caltinay 2982 step we need to add some statements; 102 ahallam 2801 \begin{python} 103 gross 2905 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 ahallam 2801 \end{python} 109 ahallam 2979 It is important to remark here that if a Dirichlet boundary condition is 110 caltinay 2982 prescribed on the same location as any Neumann boundary condition, the Neumann 111 boundary condition will be \textbf{overwritten}. This applies to Neumann 112 ahallam 2979 boundary conditions that \esc sets by default and those defined by the user. 113 ahallam 2401 114 ahallam 3370 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 gross 2905 \begin{center} 130 \includegraphics[width=4in]{figures/ttrodpyplot150} 131 caltinay 2982 \caption{Example 2: Total Energy in the Iron Rod over Time (in seconds)} 132 gross 2905 \label{fig:onedheatout1 002} 133 \end{center} 134 \end{figure} 135 136 ahallam 3370 \begin{figure}[ht] 137 gross 2905 \begin{center} 138 \includegraphics[width=4in]{figures/rodpyplot001} 139 \includegraphics[width=4in]{figures/rodpyplot050} 140 \includegraphics[width=4in]{figures/rodpyplot200} 141 ahallam 2979 \caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps 142 caltinay 2982 $1$, $50$ and $200$} 143 gross 2905 \label{fig:onedheatout 002} 144 \end{center} 145 \end{figure} 146 147 \section{For the Reader} 148 ahallam 2401 \begin{enumerate} 149 ahallam 2979 \item Move the boundary line between the two granite blocks to another part of 150 the domain. 151 gross 2905 \item Split the domain into multiple granite blocks with varying temperatures. 152 ahallam 2979 \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 jfenwick 3308 by $q_{H}$.) 156 caltinay 2982 \item Change the boundary condition for the iron rod example such that the 157 ahallam 2979 temperature 158 jfenwick 3308 at the right end is kept at a constant level $T_{ref}$, which 159 ahallam 2979 corresponds to the installation of a cooling element (Hint: Modify $q$ and 160 $r$). 161 ahallam 2401 \end{enumerate} 162