 # Contents of /trunk/doc/user/diffusion.tex

Revision 1316 - (show annotations)
Tue Sep 25 03:18:30 2007 UTC (13 years, 2 months ago) by ksteube
File MIME type: application/x-tex
File size: 18910 byte(s)
Quickly edited chapters 1 and 2 of the User Guide, but it needs more work.
Ran entire document through spell checker.


 1 % 2 % $Id$ 3 % 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5 % 6 % Copyright 2003-2007 by ACceSS MNRF 7 % Copyright 2007 by University of Queensland 8 % 9 10 % Primary Business: Queensland, Australia 11 % Licensed under the Open Software License version 3.0 12 13 % 14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 % 16 17 \section{The Diffusion Problem} 18 \label{DIFFUSION CHAP} 19 20 \begin{figure} 21 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain.eps}} 22 \caption{Temperature Diffusion Problem with Circular Heat Source} 23 \label{DIFFUSION FIG 1} 24 \end{figure} 25 26 \subsection{\label{DIFFUSION OUT SEC}Outline} 27 In this chapter we will discuss how to solve a time-dependent temperature diffusion\index{diffusion equation} PDE for 28 a given block of material. Within the block there is a heat source which drives the temperature diffusion. 29 On the surface, energy can radiate into the surrounding environment. 30 \fig{DIFFUSION FIG 1} shows the configuration. 31 32 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A 33 time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 34 We will see that at each time step a Helmholtz equation \index{Helmholtz equation} 35 must be solved. 36 The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 37 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a 38 solver for the temperature diffusion problem. 39 40 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 41 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 42 in the interior of the domain is given by 43 \begin{equation} 44 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H 45 \label{DIFFUSION TEMP EQ 1} 46 \end{equation} 47 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 48 material the parameters depend on their location in the domain. $q\hackscore H$ is 49 a heat source (or sink) within the domain. We are using the Einstein summation convention \index{summation convention} 50 as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate 51 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 52 \begin{equation} 53 q\hackscore H(x,t)= 54 \left\{ 55 \begin{array}{lcl} 56 q^c & & \|x-x^c\| \le r \\ 57 & \mbox{if} \\ 58 0 & & \mbox{else} \\ 59 \end{array} 60 \right. 61 \label{DIFFUSION TEMP EQ 1b} 62 \end{equation} 63 for all $x$ in the domain and all time $t>0$. 64 65 On the surface of the domain we are 66 specifying a radiation condition 67 which prescribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 68 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 69 \begin{equation} 70 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T) 71 \label{DIFFUSION TEMP EQ 2} 72 \end{equation} 73 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 74 $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} 75 at the surface of the domain. 76 77 To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 78 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 79 \begin{equation} 80 T(x,0)=T\hackscore{ref} 81 \label{DIFFUSION TEMP EQ 4} 82 \end{equation} 83 for all $x$ in the domain. It is pointed out that 84 the initial conditions satisfy the 85 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 86 87 The temperature is calculated at discrete time nodes $t^{(n)}$ where 88 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant. 89 In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest 90 and most robust scheme to approximate the time derivative of the the temperature is 91 the backward Euler 92 \index{backward Euler} scheme. The backward Euler 93 scheme is based 94 on the Taylor expansion of $T$ at time $t^{(n)}$: 95 \begin{equation} 96 T^{(n)}\approx T^{(n-1)}+T\hackscore{,t}^{(n)}(t^{(n)}-t^{(n-1)}) 97 =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)} 98 \label{DIFFUSION TEMP EQ 6} 99 \end{equation} 100 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 101 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 102 \begin{equation} 103 \frac{\rho c\hackscore p}{h} T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n-1)} 104 \label{DIFFUSION TEMP EQ 7} 105 \end{equation} 106 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 107 Together with the natural boundary condition 108 \begin{equation} 109 \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)}) 110 \label{DIFFUSION TEMP EQ 2222} 111 \end{equation} 112 taken from \eqn{DIFFUSION TEMP EQ 2} 113 this forms a boundary value problem that has to be solved for each time step. 114 As a first step to implement a solver for the temperature diffusion problem we will 115 first implement a solver for the boundary value problem that has to be solved at each time step. 116 117 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 118 The partial differential equation to be solved for $T^{(n)}$ has the form 119 \begin{equation} 120 \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 121 \label{DIFFUSION HELM EQ 1} 122 \end{equation} 123 and we set 124 \begin{equation} 125 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 126 \label{DIFFUSION HELM EQ 1b} 127 \end{equation} 128 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 129 takes the form 130 \begin{equation} 131 \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}\mbox{ on } \Gamma 132 \label{DIFFUSION HELM EQ 2} 133 \end{equation} 134 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 135 is called the Helmholtz equation \index{Helmholtz equation}. 136 137 We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the 138 Helmholtz equation. For a single PDE the \LinearPDE class supports the following form: 139 \begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} 140 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; . 141 \end{equation} 142 where we show only the coefficients relevant for the problem discussed here. For the general form of 143 single PDE see \eqn{LINEARPDE.SINGLE.1}. 144 The coefficients $A$, and $Y$ have to be specified through \Data objects in the 145 \Function on the PDE or objects that can be converted into such \Data objects. 146 $A$ is a \RankTwo and $D$ and $Y$ are scalar. 147 The following natural 148 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 149 \begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} 150 n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;. 151 \end{equation} 152 Notice that the coefficient $A$ is the same like in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}. 153 The coefficients $d$ and $y$ are 154 each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the 155 solution at certain locations in the domain. They have the form 156 \begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} 157 u=r \mbox{ where } q>0 158 \end{equation} 159 $r$ and $q$ are each \Scalar where $q$ is the characteristic function 160 \index{characteristic function} defining where the constraint is applied. 161 The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any other condition set by 162 \eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}. 163 The \Poisson class of the \linearPDEs module, 164 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general 165 \LinearPDE class. The \linearPDEs module provides a \Helmholtz class but 166 we will make direct use of the general \LinearPDE class. 167 168 By inspecting the Helmholtz equation \index{Helmholtz equation} 169 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and 170 substituting $u$ for $T^{(n)}$ 171 we can easily assign values to the coefficients in the 172 general PDE of the \LinearPDE class: 173 \begin{equation}\label{DIFFUSION HELM EQ 3} 174 \begin{array}{llllll} 175 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 176 d=\eta & y= g & \\ 177 \end{array} 178 \end{equation} 179 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 180 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference 181 in \escript of being not present and set to zero. As not present coefficients are not processed, 182 it is more efficient to leave a coefficient undefined instead of assigning zero to it.} 183 In this diffusion example we do not need to define a characteristic function $q$ because the 184 boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary 185 conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). 186 187 Defining and solving the Helmholtz equation is very easy now: 188 \begin{python} 189 from esys.escript import * 190 from linearPDEs import LinearPDE 191 mypde=LinearPDE(mydomain) 192 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 193 u=mypde.getSolution() 194 \end{python} 195 where we assume that \code{mydomain} is a \Domain object and 196 \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values 197 typically \code{float} or \Data objects. The \method{setValue} method 198 assigns values to the coefficients of the general PDE. The \method{getSolution} method solves 199 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 200 returning the Kronecker symbol. 201 202 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 203 If a value is assigned to a coefficient several times, the last assigned value is used when 204 the solution is calculated: 205 \begin{python} 206 mypde=LinearPDE(mydomain) 207 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 208 mypde.setValue(D=omega,Y=f,y=g) 209 mypde.setValue(d=2*eta) # overwrites d=eta 210 u=mypde.getSolution() 211 \end{python} 212 In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the 213 PDE is called symmetric if 214 \begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} 215 A\hackscore{jl}=A\hackscore{lj}\;. 216 \end{equation} 217 Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints 218 are not relevant. The Helmholtz problem is symmetric. 219 The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 220 \begin{python} 221 mypde=LinearPDE(mydomain) 222 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 223 print mypde.checkSymmetry() # returns True 224 mypde.setValue(B=kronecker(mydomain)) 225 print mypde.checkSymmetry() # returns False 226 mypde.setValue(C=kronecker(mydomain)) 227 print mypde.checkSymmetry() # returns True 228 \end{python} 229 Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for 230 testing and debugging purposes only. The \method{setSymmetryOn} method is used to 231 declare a PDE symmetric: 232 \begin{python} 233 mypde = LinearPDE(mydomain) 234 mypde.setValue(A=kappa*kronecker(mydomain)) 235 mypde.setSymmetryOn() 236 \end{python} 237 Now we want to see how we actually solve the Helmholtz equation. 238 on a rectangular domain 239 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we 240 can verify the returned solution against the exact answer. Actually, we 241 take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that 242 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 243 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 244 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 245 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 246 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 247 implements this test problem using the \finley PDE solver: 248 \begin{python} 249 from esys.escript import * 250 from esys.escript.linearPDEs import LinearPDE 251 from esys.finley import Rectangle 252 #... set some parameters ... 253 kappa=1. 254 omega=0.1 255 eta=10. 256 #... generate domain ... 257 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 258 #... open PDE and set coefficients ... 259 mypde=LinearPDE(mydomain) 260 mypde.setSymmetryOn() 261 n=mydomain.getNormal() 262 x=mydomain.getX() 263 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x, \ 264 d=eta,y=kappa*n+eta*x) 265 #... calculate error of the PDE solution ... 266 u=mypde.getSolution() 267 print "error is ",Lsup(u-x) 268 saveVTK("x0.xml",sol=u) 269 \end{python} 270 To visualise the solution `x0.~xml' just use the command 271 \begin{python} 272 mayavi -d u.xml -m SurfaceMap & 273 \end{python} 274 and it is easy to see that the solution $T=x\hackscore{0}$ is calculated. 275 276 The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}. 277 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 278 imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 279 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 280 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 281 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$). 282 However most PDE packages use an iterative solver which is terminated 283 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the 284 \method{setTolerance} of the \LinearPDE class. 285 286 \subsection{The Transition Problem} 287 \label{DIFFUSION TRANS SEC} 288 Now we are ready to solve the original time-dependent problem. The main 289 part of the script is the loop over time $t$ which takes the following form: 290 \begin{python} 291 t=0 292 T=Tref 293 mypde=LinearPDE(mydomain) 294 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 295 while t

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision