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

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


 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2018 by The University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Apache License, version 2.0 8 9 % 10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development 2012-2013 by School of Earth Sciences 12 % Development from 2014 by Centre for Geoscience Computing (GeoComp) 13 % 14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 16 \section{Example 2: One Dimensional Heat Diffusion in an Iron Rod} 17 \sslist{example02.py} 18 \label{Sec:1DHDv0} 19 20 Our second example is of a cold iron bar at a constant temperature of 21 $T_{ref}=20^{\circ} C$, see \reffig{fig:onedhdmodel}. The bar is 22 perfectly insulated on all sides with a heating element at one end keeping the 23 temperature at a constant level $T_0=100^{\circ} C$. As heat is 24 applied energy will disperse along the bar via conduction. With time the bar 25 will reach a constant temperature equivalent to that of the heat source. 26 27 \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 This problem is very similar to the example of temperature diffusion in granite 34 blocks presented in the previous Section~\ref{Sec:1DHDv00}. Thus, it is possible 35 to modify the script we have already developed for the granite blocks to suit 36 the iron bar problem. 37 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 \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 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 use \verb|whereZero| which returns the value one for those sample points where 61 the argument (almost) equals zero and the value zero elsewhere. The initial 62 temperature is set to 63 \begin{python} 64 # ... set initial temperature .... 65 T= T0*whereZero(x)+Tref*(1-whereZero(x)) 66 \end{python} 67 68 \subsection{Dirichlet Boundary Conditions} 69 In the iron rod model we want to keep the initial temperature $T_0$ on 70 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 have the value $T_0$ on the left hand side of the domain. As mentioned 73 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 78 To define a Dirichlet boundary condition we need to specify where to apply the 79 condition and determine what value the 80 solution should have at these locations. In \esc we use $q$ and $r$ to define 81 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 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 In the case of the iron rod we can set 88 \begin{python} 89 q=whereZero(x) 90 r=T0 91 \end{python} 92 to prescribe the value $T_{0}$ for the temperature at the left end of 93 the rod where $x_{0}=0$. 94 Here we use the \verb|whereZero| function again which we have already used to 95 set the initial value. 96 Notice that $r$ is set to the constant value $T_{0}$ for all sample 97 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 101 To set the Dirichlet boundary conditions for the PDE to be solved in each time 102 step we need to add some statements; 103 \begin{python} 104 mypde=LinearPDE(rod) 105 A=zeros((2,2))) 106 A[0,0]=kappa 107 q=whereZero(x) 108 mypde.setValue(A=A, D=rhocp/h, q=q, r=T0) 109 \end{python} 110 It is important to remark here that if a Dirichlet boundary condition is 111 prescribed on the same location as any Neumann boundary condition, the Neumann 112 boundary condition will be \textbf{overwritten}. This applies to Neumann 113 boundary conditions that \esc sets by default and those defined by the user. 114 115 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 \begin{center} 131 \includegraphics[width=4in]{figures/ttrodpyplot150} 132 \caption{Example 2: Total Energy in the Iron Rod over Time (in seconds)} 133 \label{fig:onedheatout1 002} 134 \end{center} 135 \end{figure} 136 137 \begin{figure}[ht] 138 \begin{center} 139 \includegraphics[width=4in]{figures/rodpyplot001} 140 \includegraphics[width=4in]{figures/rodpyplot050} 141 \includegraphics[width=4in]{figures/rodpyplot200} 142 \caption{Example 2: Temperature ($T$) distribution in the iron rod at time steps 143 $1$, $50$ and $200$} 144 \label{fig:onedheatout 002} 145 \end{center} 146 \end{figure} 147 148 \section{For the Reader} 149 \begin{enumerate} 150 \item Move the boundary line between the two granite blocks to another part of 151 the domain. 152 \item Split the domain into multiple granite blocks with varying temperatures. 153 \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 by $q_{H}$.) 157 \item Change the boundary condition for the iron rod example such that the 158 temperature 159 at the right end is kept at a constant level $T_{ref}$, which 160 corresponds to the installation of a cooling element (Hint: Modify $q$ and 161 $r$). 162 \end{enumerate} 163