ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 107 - (show annotations)
Thu Jan 27 06:21:48 2005 UTC (18 years ago) by jgs
File MIME type: application/x-tex
File size: 18183 byte(s)
commit of branch dev-01 back to main trunk on 2005-01-27

1 % $Id$
2 \chapter{How to Solve A Diffusion Problem}
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}
11 \section{\label{DIFFUSION OUT SEC}Outline}
12 In this chapter we will discuss how to solve the time dependent-temperature diffusion\index{diffusion equation} for
13 a block of material. Within the block there is a heat source which drives the temperature diffusion.
14 On the surface, energy can radiate into the surrounding environment.
15 \fig{DIFFUSION FIG 1} shows the configuration.
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 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 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a
23 solver for the temperature diffusion problem.
25 \section{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
27 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
28 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 material the parameters depend on their location in the domain. $q$ is
35 a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}
36 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 \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$.
51 On the surface of the domain we are
52 specifying a radiation condition
53 which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional
54 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 $\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 at the surface of the domain.
63 To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time
64 $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 the initial conditions satisfy the
71 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.
73 The temperature is calculated at discrete time nodes $t^{(n)}$ where
74 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.
75 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 \begin{equation}
82 T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)})
83 =T^{(n-1)} - h \cdot T\hackscore{,t}^{(n)}
84 \label{DIFFUSION TEMP EQ 6}
85 \end{equation}
86 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at
87 $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 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 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.
103 \section{\label{DIFFUSION HELM SEC}Helmholtz Problem}
104 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 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}
115 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 is called the Helmholtz equation \index{Helmholtz equation}.
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 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.
130 The form of a single PDE that can be handled by the \LinearPDE class is
131 \begin{equation}\label{EQU.FEM.1}
132 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .
133 \end{equation}
134 We show here the terms which are relevant for the Helmholtz problem.
135 The general form and systems is discussed in \Sec{SEC LinearPDE}.
136 $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients}
137 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.
148 By inspecting the Helmholtz equation \index{Helmholtz equation}
149 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 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
158 $i=j$ and $0$ otherwise.
160 We want to implement a
161 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 $g$ rather than the general form given by \eqn{EQU.FEM.1}.
164 Python's mechanism of subclasses allows us to do this in a very simple way.
165 The \Poisson class of the \linearPDEsPack module,
166 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 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 class Helmholtz(LinearPDE):
177 def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
178 ndim=self.getDim()
179 kronecker=numarray.identity(ndim)
180 self._setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)
181 \end{python}
182 \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
183 of \LinearPDE which we have imported in the first line of the script.
184 We add the method \method{setValue} to the class which overwrites the
185 \method{setValue} method of the \LinearPDE class. The new method which has the
186 parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
187 maps the parameters of the coefficients of the general PDE defined
188 in \eqn{EQU.FEM.1}. We are actually using the \method{_setValue} of
189 the \LinearPDE class (notice the leeding underscoure).
190 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 You can use your favourite editor to create and edit the file.
198 An object of the \class{Helmholtz} class is created through the statments:
199 \begin{python}
200 from mytools import Helmholtz
201 mypde = Helmholtz(mydomain)
202 mypde.setValue(kappa=10.,omega=0.1,f=12)
203 u = mypde.getSolution()
204 \end{python}
205 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 assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
209 specified, the default values $0$ are used. \footnote{It would be better to use the default value
210 \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
211 would not be processed when the PDE is evaluated}. In the fourth statement the solution of the
212 PDE is returned.
214 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 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
219 $\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 \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 from mytools import Helmholtz
225 from esys.escript import Lsup
226 from esys.finley import Rectangle
227 #... set some parameters ...
228 kappa=1.
229 omega=0.1
230 eta=10.
231 #... generate domain ...
232 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
233 #... open PDE and set coefficients ...
234 mypde=Helmholtz(mydomain)
235 n=mydomain.getNormal()
236 x=mydomain.getX()
237 mypde.setValue(kappa,omega,omega*x[0],eta,kappa*n[0]+eta*x[0])
238 #... calculate error of the PDE solution ...
239 u=mypde.getSolution()
240 print "error is ",Lsup(u-x[0])
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 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 \method{setTolerance} of the \LinearPDE class.
252 \section{The Transition Problem}
254 Now we are ready to solve the original time dependent problem. The main
255 part of the script is the loop over time $t$ which takes the following form:
256 \begin{python}
257 t=0
258 T=Tref
259 mypde=Helmholtz(mydomain)
260 while t<t_end:
261 mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)
262 T=mypde.getSolution()
263 t+=h
264 \end{python}
265 \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
266 in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde}
267 is created before the loop over time is entered while in each time step only the coefficients
268 are reset in each time step. This way some information about the representation of the PDE can be reused
269 \footnote{The efficience can be improved further by setting the coefficients in the operator
270 \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only updating the coefficients
271 in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
272 method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator.}. The variable \var{T}
273 holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
274 Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
275 temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
276 initial temperature distribution, which equal to $T\hackscore{ref}$.
278 The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
279 in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
280 \var{q0} is a fixed constant. The following script defines \var{q} as desired:
281 \begin{python}
282 from esys.escript import length
283 xc=[0.02,0.002]
284 r=0.001
285 x=mydomain.getX()
286 q=q0*(length(x-xc)-r).whereNegative()
287 \end{python}
288 \var{x} is a \Data class object of
289 the \escript module defining locations in the \Domain \var{mydomain}.
290 The \function{length()} imported from the \escript module returns the
291 Euclidean norm:
292 \begin{equation}\label{DIFFUSION HELM EQ 3aba}
293 \|y\|=\sqrt{
294 y\hackscore{i}
295 y\hackscore{i}
296 } = \function{esys.escript.length}(y)
297 \end{equation}
298 So \code{length(x-xc)} calculates the distances
299 of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
300 Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
301 converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of
302 a \Data class object, in this case the result of the expression
303 \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
304 zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
306 Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
307 \index{scripts!\file{diffusion.py}}:
308 \begin{python}
309 from mytools import Helmholtz
310 from esys.escript import Lsup
311 from esys.finley import Rectangle
312 #... set some parameters ...
313 xc=[0.02,0.002]
314 r=0.001
315 qc=50.e6
316 Tref=0.
317 rhocp=2.6e6
318 eta=75.
319 kappa=240.
320 tend=5.
321 # ... time, time step size and counter ...
322 t=0
323 h=0.1
324 i=0
325 #... generate domain ...
326 mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
327 #... open PDE ...
328 mypde=Helmholtz(mydomain)
329 # ... set heat source: ....
330 x=mydomain.getX()
331 q=qc*(length(x-xc)-r).whereNegative()
332 # ... set initial temperature ....
333 T=Tref
334 # ... start iteration:
335 while t<tend:
336 i+=1
337 t+=h
338 print "time step :",t
339 mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
340 T=mypde.getSolution()
341 T.saveDX("T%d.dx"%i)
342 \end{python}
343 The script will create the files \file{T.1.dx},
344 \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the
345 temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.
347 \begin{figure}
348 \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
349 \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
350 \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
351 \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
352 \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
353 \label{DIFFUSION FIG 2}
354 \end{figure}
356 An easy way to visualize the results is the command
357 \begin{verbatim}
358 dx -edit diffusion.net &
359 \end{verbatim}
360 where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.
361 Use the \texttt{Sequencer} to move forward and and backwards in time.
362 \fig{DIFFUSION FIG 2} shows the result for some selected time steps.


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

  ViewVC Help
Powered by ViewVC 1.1.26