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

Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 9 months ago) by ksteube
File MIME type: application/x-tex
File size: 18888 byte(s)
Copyright updated in all files


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

## Properties

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