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

Revision 565 - (hide annotations)
Mon Feb 27 00:04:17 2006 UTC (15 years, 7 months ago) by gross
File MIME type: application/x-tex
File size: 18169 byte(s)
some updates

 1 jgs 102 % $Id$ 2 jgs 121 \section{The Diffusion Problem} 3 jgs 102 \label{DIFFUSION CHAP} 4 5 \begin{figure} 6 \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}} 7 \caption{Temperature Diffusion Problem with Circular Heat Source} 8 \label{DIFFUSION FIG 1} 9 \end{figure} 10 11 jgs 121 \subsection{\label{DIFFUSION OUT SEC}Outline} 12 jgs 107 In this chapter we will discuss how to solve the time dependent-temperature diffusion\index{diffusion equation} for 13 jgs 102 a block of material. Within the block there is a heat source which drives the temperature diffusion. 14 jgs 107 On the surface, energy can radiate into the surrounding environment. 15 jgs 102 \fig{DIFFUSION FIG 1} shows the configuration. 16 17 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A 18 time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 19 jgs 107 We will see that at each time step a Helmholtz equation \index{Helmholtz equation} 20 must be solved. 21 The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 22 jgs 102 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a 23 solver for the temperature diffusion problem. 24 25 jgs 121 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 26 jgs 102 27 jgs 107 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 28 jgs 102 in the interior of the domain is given by 29 \begin{equation} 30 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q 31 \label{DIFFUSION TEMP EQ 1} 32 \end{equation} 33 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 34 jgs 107 material the parameters depend on their location in the domain. $q$ is 35 jgs 102 a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 36 jgs 107 as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant heat production rate 37 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 38 jgs 102 \begin{equation} 39 q(x,t)= 40 \left\{ 41 \begin{array}{lcl} 42 q^c & & \|x-x^c\| \le r \\ 43 & \mbox{if} \\ 44 0 & & \mbox{else} \\ 45 \end{array} 46 \right. 47 \label{DIFFUSION TEMP EQ 1b} 48 \end{equation} 49 for all $x$ in the domain and all time $t>0$. 50 51 On the surface of the domain we are 52 jgs 107 specifying a radiation condition 53 which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 54 jgs 102 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 55 \begin{equation} 56 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T) 57 \label{DIFFUSION TEMP EQ 2} 58 \end{equation} 59 jgs 107 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 60 As usual $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} 61 jgs 102 at the surface of the domain. 62 63 jgs 107 To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 64 jgs 102 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 65 \begin{equation} 66 T(x,0)=T\hackscore{ref} 67 \label{DIFFUSION TEMP EQ 4} 68 \end{equation} 69 for all $x$ in the domain. It is pointed out that 70 jgs 107 the initial conditions satisfy the 71 jgs 102 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 72 73 jgs 107 The temperature is calculated at discrete time nodes $t^{(n)}$ where 74 jgs 102 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant. 75 jgs 107 In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest 76 and most robust scheme to approximate the time derivative of the the temperature is 77 the backward Euler 78 \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler 79 scheme is based 80 on the Taylor expansion of $T$ at time $t^{(n)}$: 81 jgs 102 \begin{equation} 82 T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)}) 83 jgs 107 =T^{(n-1)} - h \cdot T\hackscore{,t}^{(n)} 84 jgs 102 \label{DIFFUSION TEMP EQ 6} 85 \end{equation} 86 jgs 107 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 87 jgs 102 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 88 \begin{equation} 89 \frac{\rho c\hackscore p}{h} T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q + \frac{\rho c\hackscore p}{h} T^{(n-1)} 90 \label{DIFFUSION TEMP EQ 7} 91 \end{equation} 92 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 93 jgs 107 Together with the natural boundary condition 94 \begin{equation} 95 \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)}) 96 \label{DIFFUSION TEMP EQ 2222} 97 \end{equation} 98 taken from \eqn{DIFFUSION TEMP EQ 2} 99 jgs 102 this forms a boundary value problem that has to be solved for each time step. 100 As a first step to implement a solver for the temperature diffusion problem we will 101 first implement a solver for the boundary value problem that has to be solved at each time step. 102 103 jgs 121 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 104 jgs 102 The partial differential equation to be solved for $T^{(n)}$ has the form 105 \begin{equation} 106 \omega u - (\kappa u\hackscore{,i})\hackscore{,i} = f 107 \label{DIFFUSION HELM EQ 1} 108 \end{equation} 109 where $u$ plays the role of $T^{(n)}$ and we set 110 \begin{equation} 111 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 112 \label{DIFFUSION HELM EQ 1b} 113 \end{equation} 114 jgs 107 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 115 jgs 102 takes the form 116 \begin{equation} 117 \kappa u\hackscore{,i} n\hackscore{i} = g - \eta u\mbox{ on } \Gamma 118 \label{DIFFUSION HELM EQ 2} 119 \end{equation} 120 The partial differential 121 \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 122 jgs 107 is called the Helmholtz equation \index{Helmholtz equation}. 123 jgs 102 124 We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the 125 Helmholtz equation. We have used a special case of the \LinearPDE class, namely the 126 \Poisson class already in \Chap{FirstSteps}. 127 jgs 107 Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation 128 and use the \method{getSolution} method of parent class to actually solve the problem. 129 jgs 102 130 jgs 107 The form of a single PDE that can be handled by the \LinearPDE class is 131 jgs 102 \begin{equation}\label{EQU.FEM.1} 132 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; . 133 \end{equation} 134 jgs 107 We show here the terms which are relevant for the Helmholtz problem. 135 jgs 102 The general form and systems is discussed in \Sec{SEC LinearPDE}. 136 jgs 107 $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients} 137 jgs 102 Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar. 138 They may be constant or may depend on their 139 location in the domain but must not depend on the unknown solution $u$. 140 The following natural boundary conditions \index{boundary condition!natural} that 141 are used in the \LinearPDE class have the form 142 \begin{equation}\label{EQU.FEM.2} 143 n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y \;. 144 \end{equation} 145 where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that 146 the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients. 147 148 jgs 107 By inspecting the Helmholtz equation \index{Helmholtz equation} 149 jgs 102 we can easily assign values to the coefficients in the 150 general PDE of the \LinearPDE class: 151 \begin{equation}\label{DIFFUSION HELM EQ 3} 152 \begin{array}{llllll} 153 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 154 d=\eta & y= g & \\ 155 \end{array} 156 \end{equation} 157 jgs 107 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 158 jgs 102 $i=j$ and $0$ otherwise. 159 160 We want to implement a 161 jgs 107 new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but 162 is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$, 163 jgs 102 $g$ rather than the general form given by \eqn{EQU.FEM.1}. 164 jgs 107 Python's mechanism of subclasses allows us to do this in a very simple way. 165 gross 565 The \Poisson class of the \linearPDEs module, 166 jgs 107 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general 167 \LinearPDE class. That means that all methods (such as the \method{getSolution}) 168 from the parent class \LinearPDE are available for any \Poisson object. However, new 169 jgs 102 methods can be added and methods of the parent class can be redefined. In fact, 170 the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign 171 values to the coefficients of the PDE. This is exactly what we will do when we define 172 our new \class{Helmholtz} class: 173 \begin{python} 174 from esys.linearPDEs import LinearPDE 175 import numarray 176 jgs 107 class Helmholtz(LinearPDE): 177 jgs 102 def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0) 178 jgs 107 ndim=self.getDim() 179 kronecker=numarray.identity(ndim) 180 jgs 121 self._LinearPDE_setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g) 181 jgs 102 \end{python} 182 \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass 183 jgs 107 of \LinearPDE which we have imported in the first line of the script. 184 jgs 102 We add the method \method{setValue} to the class which overwrites the 185 jgs 107 \method{setValue} method of the \LinearPDE class. The new method which has the 186 jgs 102 parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments 187 maps the parameters of the coefficients of the general PDE defined 188 jgs 121 in \eqn{EQU.FEM.1}. We are actually using the \method{_LinearPDE__setValue} of 189 the \LinearPDE class. 190 jgs 107 The coefficient \var{A} involves the Kronecker symbol. We use the 191 \numarray function \function{identity} which returns a square matrix with ones on the 192 main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix. 193 The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the 194 PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to 195 put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory. 196 jgs 102 You can use your favourite editor to create and edit the file. 197 198 An object of the \class{Helmholtz} class is created through the statments: 199 \begin{python} 200 jgs 107 from mytools import Helmholtz 201 mypde = Helmholtz(mydomain) 202 jgs 102 mypde.setValue(kappa=10.,omega=0.1,f=12) 203 jgs 107 u = mypde.getSolution() 204 jgs 102 \end{python} 205 jgs 107 In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure 206 that \file{mytools.py} is in the directory from where you started Python. 207 \var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are 208 jgs 102 assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are 209 jgs 107 specified, the default values $0$ are used. \footnote{It would be better to use the default value 210 jgs 102 \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and 211 jgs 107 would not be processed when the PDE is evaluated}. In the fourth statement the solution of the 212 jgs 102 PDE is returned. 213 214 jgs 107 To test our \class{Helmholtz} class on a rectangular domain 215 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we 216 we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that 217 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 218 jgs 102 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 219 jgs 107 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 220 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py} 221 jgs 102 \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory 222 implements this test problem using the \finley PDE solver: 223 \begin{python} 224 jgs 107 from mytools import Helmholtz 225 from esys.escript import Lsup 226 from esys.finley import Rectangle 227 jgs 102 #... set some parameters ... 228 jgs 107 kappa=1. 229 jgs 102 omega=0.1 230 eta=10. 231 #... generate domain ... 232 jgs 107 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 233 jgs 102 #... open PDE and set coefficients ... 234 mypde=Helmholtz(mydomain) 235 n=mydomain.getNormal() 236 x=mydomain.getX() 237 jgs 107 mypde.setValue(kappa,omega,omega*x,eta,kappa*n+eta*x) 238 jgs 102 #... calculate error of the PDE solution ... 239 u=mypde.getSolution() 240 print "error is ",Lsup(u-x) 241 \end{python} 242 The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}. 243 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 244 jgs 107 is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument. 245 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 246 used to approximate the solution and our solution is a linear function of the spatial coordinates one might 247 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$). 248 However most PDE packages use an iterative solver which is terminated 249 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the 250 jgs 102 \method{setTolerance} of the \LinearPDE class. 251 252 jgs 121 \subsection{The Transition Problem} 253 jgs 102 \label{DIFFUSION TRANS SEC} 254 Now we are ready to solve the original time dependent problem. The main 255 jgs 107 part of the script is the loop over time $t$ which takes the following form: 256 jgs 102 \begin{python} 257 jgs 107 t=0 258 T=Tref 259 jgs 102 mypde=Helmholtz(mydomain) 260 while t

## Properties

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