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

Revision 660 - (show annotations)
Fri Mar 24 08:43:21 2006 UTC (14 years, 8 months ago) by gross
File MIME type: application/x-tex
File size: 18645 byte(s)
that does not look to bad now although more stuff could be added.

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

## Properties

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