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

Revision 568 - (hide annotations)
Tue Feb 28 03:59:06 2006 UTC (13 years, 7 months ago) by gross
File MIME type: application/x-tex
File size: 18829 byte(s)
update
 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 107 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 27 jgs 102 in the interior of the domain is given by 28 \begin{equation} 29 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q 30 \label{DIFFUSION TEMP EQ 1} 31 \end{equation} 32 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 33 jgs 107 material the parameters depend on their location in the domain. $q$ is 34 jgs 102 a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 35 jgs 107 as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant heat production rate 36 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 37 jgs 102 \begin{equation} 38 q(x,t)= 39 \left\{ 40 \begin{array}{lcl} 41 q^c & & \|x-x^c\| \le r \\ 42 & \mbox{if} \\ 43 0 & & \mbox{else} \\ 44 \end{array} 45 \right. 46 \label{DIFFUSION TEMP EQ 1b} 47 \end{equation} 48 for all $x$ in the domain and all time $t>0$. 49 50 On the surface of the domain we are 51 jgs 107 specifying a radiation condition 52 which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 53 jgs 102 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 54 \begin{equation} 55 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T) 56 \label{DIFFUSION TEMP EQ 2} 57 \end{equation} 58 jgs 107 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 59 As usual $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} 60 jgs 102 at the surface of the domain. 61 62 jgs 107 To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 63 jgs 102 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 64 \begin{equation} 65 T(x,0)=T\hackscore{ref} 66 \label{DIFFUSION TEMP EQ 4} 67 \end{equation} 68 for all $x$ in the domain. It is pointed out that 69 jgs 107 the initial conditions satisfy the 70 jgs 102 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 71 72 jgs 107 The temperature is calculated at discrete time nodes $t^{(n)}$ where 73 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. 74 jgs 107 In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest 75 and most robust scheme to approximate the time derivative of the the temperature is 76 the backward Euler 77 \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler 78 scheme is based 79 on the Taylor expansion of $T$ at time $t^{(n)}$: 80 jgs 102 \begin{equation} 81 T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)}) 82 jgs 107 =T^{(n-1)} - h \cdot T\hackscore{,t}^{(n)} 83 jgs 102 \label{DIFFUSION TEMP EQ 6} 84 \end{equation} 85 jgs 107 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 86 jgs 102 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 87 \begin{equation} 88 \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)} 89 \label{DIFFUSION TEMP EQ 7} 90 \end{equation} 91 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 92 jgs 107 Together with the natural boundary condition 93 \begin{equation} 94 \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)}) 95 \label{DIFFUSION TEMP EQ 2222} 96 \end{equation} 97 taken from \eqn{DIFFUSION TEMP EQ 2} 98 jgs 102 this forms a boundary value problem that has to be solved for each time step. 99 As a first step to implement a solver for the temperature diffusion problem we will 100 first implement a solver for the boundary value problem that has to be solved at each time step. 101 102 jgs 121 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 103 jgs 102 The partial differential equation to be solved for $T^{(n)}$ has the form 104 \begin{equation} 105 \omega u - (\kappa u\hackscore{,i})\hackscore{,i} = f 106 \label{DIFFUSION HELM EQ 1} 107 \end{equation} 108 where $u$ plays the role of $T^{(n)}$ and we set 109 \begin{equation} 110 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 111 \label{DIFFUSION HELM EQ 1b} 112 \end{equation} 113 jgs 107 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 114 jgs 102 takes the form 115 \begin{equation} 116 \kappa u\hackscore{,i} n\hackscore{i} = g - \eta u\mbox{ on } \Gamma 117 \label{DIFFUSION HELM EQ 2} 118 \end{equation} 119 gross 568 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 120 jgs 107 is called the Helmholtz equation \index{Helmholtz equation}. 121 jgs 102 122 gross 568 We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the 123 Helmholtz equation. For a single PDE the \LinearPDE class supports the following form: 124 \begin{equation}\label{LINEARPDE.SINGLE.1} 125 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+(B\hackscore{j} u)\hackscore{,j}+C\hackscore{l} u\hackscore{,l}+D u =-X\hackscore{j,j}+Y \; . 126 jgs 102 \end{equation} 127 gross 568 The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the 128 \Function on the PDE or objects that can be converted into such \Data objects. 129 $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar. 130 The following natural 131 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 132 \begin{equation}\label{LINEARPDE.SINGLE.2} 133 n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y \;. 134 jgs 102 \end{equation} 135 gross 568 Notice that the coefficients $A$, $B$ and $X$ are the same like in the PDE~\eqn{LINEARPDE.SINGLE.1}. The coefficients $d$ and $y$ are 136 each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the 137 solution at certain locations in the domain. They have the form 138 \begin{equation}\label{LINEARPDE.SINGLE.3} 139 u=r \mbox{ where } q>0 140 \end{equation} 141 $r$ and $q$ are each \Scalar where $q$ is the characteristic function 142 \index{characteristic function} defining where the constraint is applied. 143 The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1} 144 or \eqn{LINEARPDE.SINGLE.2}. 145 The \Poisson class of the \linearPDEs module, 146 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general 147 \LinearPDE class. The \linearPDEs module provides a \Helmholtz class but 148 we will make direct use of the general \LinearPDE class. 149 jgs 102 150 jgs 107 By inspecting the Helmholtz equation \index{Helmholtz equation} 151 jgs 102 we can easily assign values to the coefficients in the 152 general PDE of the \LinearPDE class: 153 \begin{equation}\label{DIFFUSION HELM EQ 3} 154 \begin{array}{llllll} 155 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 156 d=\eta & y= g & \\ 157 \end{array} 158 \end{equation} 159 jgs 107 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 160 gross 568 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present\footnote{There is a difference 161 in \escript of being not present and set to zero. As not present coefficients are not processed, 162 it is more efficient to leave a coefficient undefined insted assigning zero to it.} 163 jgs 102 164 gross 568 Defining and solving the Helmholtz equation is very easy now: 165 \begin{python} 166 from esys.escript import * 167 from linearPDEs import LinearPDE 168 mypde=LinearPDE(mydomain) 169 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 170 u=mypde.getSolution() 171 \end{python} 172 where we assume that \code{mydomain} is a \Domain object and 173 \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values 174 typically \code{float} or \Data objects. The \method{setValue} method 175 assigns values to the coefficients of the general PDE. The \method{getSolution} method solves 176 the PDE and returns the solution \code{u} of the PDE. \code{kronecker} is \escript function 177 returning the Kronecker symbol. 178 179 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 180 If a value is assigned to a coefficint several times, the last assigned value is used when 181 the solution is calculated: 182 \begin{python} 183 mypde=LinearPDE(mydomain) 184 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 185 mypde.setValue(D=omega,Y=f,y=g) 186 mypde.setValue(d=2*eta) # overwrites d=eta 187 u=mypde.getSolution() 188 \end{python} 189 In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the 190 PDE is called symmetric if 191 \begin{equation}\label{LINEARPDE.SINGLE.4} 192 A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j} \;. 193 \end{equation} 194 Note that $D$ and $d$ may have any value and the right hand side $X$, $Y$, $y$ as well as the constraints 195 are not relevant. The Helmholtz problem is symmetric. 196 The \LinearPDE class provides the method \code{checkSymmetry} method to check if the given PDE is symmetric. 197 \begin{python} 198 mypde=LinearPDE(mydomain) 199 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 200 print mypde.checkSymmetry() # returns True 201 mypde.setValue(B=kronecker(mydomain)) 202 print mypde.checkSymmetry() # returns False 203 mypde.setValue(C=kronecker(mydomain)) 204 print mypde.checkSymmetry() # returns True 205 \end{python} 206 Unfortunately, a \code{checkSymmetry} is very expensive and is recommendable to use for 207 testing and debugging purposes only. The \code{setSymmetryOn} method is used to 208 declare a PDE symmetric: 209 \begin{python} 210 mypde = LinearPDE(mydomain) 211 mypde.setValue(A=kappa*kronecker(mydomain)) 212 mypde.setSymmetryOn() 213 \end{python} 214 215 216 217 218 219 220 221 ===== 222 223 Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation 224 and use the \method{getSolution} method of parent class to actually solve the problem. 225 226 jgs 102 We want to implement a 227 jgs 107 new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but 228 is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$, 229 jgs 102 $g$ rather than the general form given by \eqn{EQU.FEM.1}. 230 jgs 107 Python's mechanism of subclasses allows us to do this in a very simple way. 231 gross 568 That means that all methods (such as the \method{getSolution}) 232 jgs 107 from the parent class \LinearPDE are available for any \Poisson object. However, new 233 jgs 102 methods can be added and methods of the parent class can be redefined. In fact, 234 the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign 235 values to the coefficients of the PDE. This is exactly what we will do when we define 236 our new \class{Helmholtz} class: 237 238 jgs 107 To test our \class{Helmholtz} class on a rectangular domain 239 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we 240 we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that 241 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 242 jgs 102 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 243 jgs 107 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 244 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py} 245 gross 568 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 246 jgs 102 implements this test problem using the \finley PDE solver: 247 \begin{python} 248 gross 568 from esys.escript import * 249 from linearPDEs import LinearPDE 250 jgs 107 from esys.finley import Rectangle 251 jgs 102 #... set some parameters ... 252 jgs 107 kappa=1. 253 jgs 102 omega=0.1 254 eta=10. 255 #... generate domain ... 256 jgs 107 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 257 jgs 102 #... open PDE and set coefficients ... 258 gross 568 mypde=LinearPDE(mydomain) 259 mypde.setSymmetryOn() 260 jgs 102 n=mydomain.getNormal() 261 x=mydomain.getX() 262 gross 568 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x,d=eta,y=kappa*n+eta*x) 263 jgs 102 #... calculate error of the PDE solution ... 264 u=mypde.getSolution() 265 print "error is ",Lsup(u-x) 266 \end{python} 267 The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}. 268 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 269 jgs 107 is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument. 270 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 271 used to approximate the solution and our solution is a linear function of the spatial coordinates one might 272 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$). 273 However most PDE packages use an iterative solver which is terminated 274 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the 275 jgs 102 \method{setTolerance} of the \LinearPDE class. 276 277 jgs 121 \subsection{The Transition Problem} 278 jgs 102 \label{DIFFUSION TRANS SEC} 279 Now we are ready to solve the original time dependent problem. The main 280 jgs 107 part of the script is the loop over time $t$ which takes the following form: 281 jgs 102 \begin{python} 282 jgs 107 t=0 283 T=Tref 284 jgs 102 mypde=Helmholtz(mydomain) 285 while t

Properties

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