/[escript]/trunk/doc/user/diffusion.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1387 - (show annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 9 months ago) by trankine
Original Path: temp/doc/user/diffusion.tex
File MIME type: application/x-tex
File size: 18910 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
1 %
2 % $Id$
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 %
6 % Copyright 2003-2007 by ACceSS MNRF
7 % Copyright 2007 by University of Queensland
8 %
9 % http://esscc.uq.edu.au
10 % Primary Business: Queensland, Australia
11 % Licensed under the Open Software License version 3.0
12 % http://www.opensource.org/licenses/osl-3.0.php
13 %
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 %
16
17 \section{The Diffusion Problem}
18 \label{DIFFUSION CHAP}
19
20 \begin{figure}
21 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain.eps}}
22 \caption{Temperature Diffusion Problem with Circular Heat Source}
23 \label{DIFFUSION FIG 1}
24 \end{figure}
25
26 \subsection{\label{DIFFUSION OUT SEC}Outline}
27 In this chapter we will discuss how to solve a time-dependent temperature diffusion\index{diffusion equation} PDE for
28 a given block of material. Within the block there is a heat source which drives the temperature diffusion.
29 On the surface, energy can radiate into the surrounding environment.
30 \fig{DIFFUSION FIG 1} shows the configuration.
31
32 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A
33 time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$.
34 We will see that at each time step a Helmholtz equation \index{Helmholtz equation}
35 must be solved.
36 The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}.
37 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a
38 solver for the temperature diffusion problem.
39
40 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
41 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
42 in the interior of the domain is given by
43 \begin{equation}
44 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H
45 \label{DIFFUSION TEMP EQ 1}
46 \end{equation}
47 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite
48 material the parameters depend on their location in the domain. $q\hackscore H$ is
49 a heat source (or sink) within the domain. We are using the Einstein summation convention \index{summation convention}
50 as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate
51 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:
52 \begin{equation}
53 q\hackscore H(x,t)=
54 \left\{
55 \begin{array}{lcl}
56 q^c & & \|x-x^c\| \le r \\
57 & \mbox{if} \\
58 0 & & \mbox{else} \\
59 \end{array}
60 \right.
61 \label{DIFFUSION TEMP EQ 1b}
62 \end{equation}
63 for all $x$ in the domain and all time $t>0$.
64
65 On the surface of the domain we are
66 specifying a radiation condition
67 which prescribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional
68 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$:
69 \begin{equation}
70 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T)
71 \label{DIFFUSION TEMP EQ 2}
72 \end{equation}
73 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium.
74 $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field}
75 at the surface of the domain.
76
77 To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time
78 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature:
79 \begin{equation}
80 T(x,0)=T\hackscore{ref}
81 \label{DIFFUSION TEMP EQ 4}
82 \end{equation}
83 for all $x$ in the domain. It is pointed out that
84 the initial conditions satisfy the
85 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.
86
87 The temperature is calculated at discrete time nodes $t^{(n)}$ where
88 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.
89 In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest
90 and most robust scheme to approximate the time derivative of the the temperature is
91 the backward Euler
92 \index{backward Euler} scheme. The backward Euler
93 scheme is based
94 on the Taylor expansion of $T$ at time $t^{(n)}$:
95 \begin{equation}
96 T^{(n)}\approx T^{(n-1)}+T\hackscore{,t}^{(n)}(t^{(n)}-t^{(n-1)})
97 =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)}
98 \label{DIFFUSION TEMP EQ 6}
99 \end{equation}
100 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at
101 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$
102 \begin{equation}
103 \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)}
104 \label{DIFFUSION TEMP EQ 7}
105 \end{equation}
106 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.
107 Together with the natural boundary condition
108 \begin{equation}
109 \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)})
110 \label{DIFFUSION TEMP EQ 2222}
111 \end{equation}
112 taken from \eqn{DIFFUSION TEMP EQ 2}
113 this forms a boundary value problem that has to be solved for each time step.
114 As a first step to implement a solver for the temperature diffusion problem we will
115 first implement a solver for the boundary value problem that has to be solved at each time step.
116
117 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem}
118 The partial differential equation to be solved for $T^{(n)}$ has the form
119 \begin{equation}
120 \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f
121 \label{DIFFUSION HELM EQ 1}
122 \end{equation}
123 and we set
124 \begin{equation}
125 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.
126 \label{DIFFUSION HELM EQ 1b}
127 \end{equation}
128 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}
129 takes the form
130 \begin{equation}
131 \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}\mbox{ on } \Gamma
132 \label{DIFFUSION HELM EQ 2}
133 \end{equation}
134 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
135 is called the Helmholtz equation \index{Helmholtz equation}.
136
137 We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the
138 Helmholtz equation. For a single PDE the \LinearPDE class supports the following form:
139 \begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL}
140 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; .
141 \end{equation}
142 where we show only the coefficients relevant for the problem discussed here. For the general form of
143 single PDE see \eqn{LINEARPDE.SINGLE.1}.
144 The coefficients $A$, and $Y$ have to be specified through \Data objects in the
145 \Function on the PDE or objects that can be converted into such \Data objects.
146 $A$ is a \RankTwo and $D$ and $Y$ are scalar.
147 The following natural
148 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
149 \begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL}
150 n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;.
151 \end{equation}
152 Notice that the coefficient $A$ is the same like in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}.
153 The coefficients $d$ and $y$ are
154 each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the
155 solution at certain locations in the domain. They have the form
156 \begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL}
157 u=r \mbox{ where } q>0
158 \end{equation}
159 $r$ and $q$ are each \Scalar where $q$ is the characteristic function
160 \index{characteristic function} defining where the constraint is applied.
161 The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any other condition set by
162 \eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}.
163 The \Poisson class of the \linearPDEs module,
164 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
165 \LinearPDE class. The \linearPDEs module provides a \Helmholtz class but
166 we will make direct use of the general \LinearPDE class.
167
168 By inspecting the Helmholtz equation \index{Helmholtz equation}
169 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and
170 substituting $u$ for $T^{(n)}$
171 we can easily assign values to the coefficients in the
172 general PDE of the \LinearPDE class:
173 \begin{equation}\label{DIFFUSION HELM EQ 3}
174 \begin{array}{llllll}
175 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\
176 d=\eta & y= g & \\
177 \end{array}
178 \end{equation}
179 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
180 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference
181 in \escript of being not present and set to zero. As not present coefficients are not processed,
182 it is more efficient to leave a coefficient undefined instead of assigning zero to it.}
183 In this diffusion example we do not need to define a characteristic function $q$ because the
184 boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary
185 conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}).
186
187 Defining and solving the Helmholtz equation is very easy now:
188 \begin{python}
189 from esys.escript import *
190 from linearPDEs import LinearPDE
191 mypde=LinearPDE(mydomain)
192 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g)
193 u=mypde.getSolution()
194 \end{python}
195 where we assume that \code{mydomain} is a \Domain object and
196 \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values
197 typically \code{float} or \Data objects. The \method{setValue} method
198 assigns values to the coefficients of the general PDE. The \method{getSolution} method solves
199 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function
200 returning the Kronecker symbol.
201
202 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily.
203 If a value is assigned to a coefficient several times, the last assigned value is used when
204 the solution is calculated:
205 \begin{python}
206 mypde=LinearPDE(mydomain)
207 mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
208 mypde.setValue(D=omega,Y=f,y=g)
209 mypde.setValue(d=2*eta) # overwrites d=eta
210 u=mypde.getSolution()
211 \end{python}
212 In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the
213 PDE is called symmetric if
214 \begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL}
215 A\hackscore{jl}=A\hackscore{lj}\;.
216 \end{equation}
217 Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints
218 are not relevant. The Helmholtz problem is symmetric.
219 The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric.
220 \begin{python}
221 mypde=LinearPDE(mydomain)
222 mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
223 print mypde.checkSymmetry() # returns True
224 mypde.setValue(B=kronecker(mydomain)[0])
225 print mypde.checkSymmetry() # returns False
226 mypde.setValue(C=kronecker(mydomain)[0])
227 print mypde.checkSymmetry() # returns True
228 \end{python}
229 Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for
230 testing and debugging purposes only. The \method{setSymmetryOn} method is used to
231 declare a PDE symmetric:
232 \begin{python}
233 mypde = LinearPDE(mydomain)
234 mypde.setValue(A=kappa*kronecker(mydomain))
235 mypde.setSymmetryOn()
236 \end{python}
237 Now we want to see how we actually solve the Helmholtz equation.
238 on a rectangular domain
239 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we
240 can verify the returned solution against the exact answer. Actually, we
241 take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that
242 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
243 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
244 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
245 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py}
246 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
247 implements this test problem using the \finley PDE solver:
248 \begin{python}
249 from esys.escript import *
250 from esys.escript.linearPDEs import LinearPDE
251 from esys.finley import Rectangle
252 #... set some parameters ...
253 kappa=1.
254 omega=0.1
255 eta=10.
256 #... generate domain ...
257 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
258 #... open PDE and set coefficients ...
259 mypde=LinearPDE(mydomain)
260 mypde.setSymmetryOn()
261 n=mydomain.getNormal()
262 x=mydomain.getX()
263 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \
264 d=eta,y=kappa*n[0]+eta*x[0])
265 #... calculate error of the PDE solution ...
266 u=mypde.getSolution()
267 print "error is ",Lsup(u-x[0])
268 saveVTK("x0.xml",sol=u)
269 \end{python}
270 To visualise the solution `x0.~xml' just use the command
271 \begin{python}
272 mayavi -d u.xml -m SurfaceMap &
273 \end{python}
274 and it is easy to see that the solution $T=x\hackscore{0}$ is calculated.
275
276 The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}.
277 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
278 imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument.
279 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
280 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might
281 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
282 However most PDE packages use an iterative solver which is terminated
283 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
284 \method{setTolerance} of the \LinearPDE class.
285
286 \subsection{The Transition Problem}
287 \label{DIFFUSION TRANS SEC}
288 Now we are ready to solve the original time-dependent problem. The main
289 part of the script is the loop over time $t$ which takes the following form:
290 \begin{python}
291 t=0
292 T=Tref
293 mypde=LinearPDE(mydomain)
294 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
295 while t<t_end:
296 mypde.setValue(Y=q+rhocp/h*T)
297 T=mypde.getSolution()
298 t+=h
299 \end{python}
300 \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
301 in the domain and \var{h} is the time step size.
302 The variable \var{T}
303 holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
304 Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
305 temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
306 initial temperature distribution, which equal to $T\hackscore{ref}$.
307 The \LinearPDE class object \var{mypde}
308 and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient
309 $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information
310 from previous time steps as much as possible.
311
312 The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
313 in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
314 \var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired:
315 \begin{python}
316 from esys.escript import length,whereNegative
317 xc=[0.02,0.002]
318 r=0.001
319 x=mydomain.getX()
320 qH=q0*whereNegative(length(x-xc)-r)
321 \end{python}
322 \var{x} is a \Data class object of
323 the \escript module defining locations in the \Domain \var{mydomain}.
324 The \function{length()} imported from the \escript module returns the
325 Euclidean norm:
326 \begin{equation}\label{DIFFUSION HELM EQ 3aba}
327 \|y\|=\sqrt{
328 y\hackscore{i}
329 y\hackscore{i}
330 } = \function{esys.escript.length}(y)
331 \end{equation}
332 So \code{length(x-xc)} calculates the distances
333 of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
334 Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is automatically
335 converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative}
336 applied to
337 \code{length(x-xc)-r}, returns a \Data object which is equal to one where the object is negative (inside the circle) and
338 zero elsewhere. After multiplication with \var{qc} we get a function with the desired property of having value \var{qc} inside
339 the circle and zero elsewhere.
340
341 Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
342 \index{scripts!\file{diffusion.py}}:
343 \begin{python}
344 from esys.escript import *
345 from esys.escript.linearPDEs import LinearPDE
346 from esys.finley import Rectangle
347 #... set some parameters ...
348 xc=[0.02,0.002]
349 r=0.001
350 qc=50.e6
351 Tref=0.
352 rhocp=2.6e6
353 eta=75.
354 kappa=240.
355 tend=5.
356 # ... time, time step size and counter ...
357 t=0
358 h=0.1
359 i=0
360 #... generate domain ...
361 mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
362 #... open PDE ...
363 mypde=LinearPDE(mydomain)
364 mypde.setSymmetryOn()
365 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
366 # ... set heat source: ....
367 x=mydomain.getX()
368 qH=qc*whereNegative(length(x-xc)-r)
369 # ... set initial temperature ....
370 T=Tref
371 # ... start iteration:
372 while t<tend:
373 i+=1
374 t+=h
375 print "time step :",t
376 mypde.setValue(Y=qH+rhocp/h*T)
377 T=mypde.getSolution()
378 saveVTK("T.%d.xml"%i,temp=T)
379 \end{python}
380 The script will create the files \file{T.1.xml},
381 \file{T.2.xml}, $\ldots$, \file{T.50.xml} in the directory where the script has been started. The files give the
382 temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \VTK file format.
383
384 \begin{figure}
385 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes1.eps}}
386 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes16.eps}}
387 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes32.eps}}
388 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes48.eps}}
389 \caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
390 \label{DIFFUSION FIG 2}
391 \end{figure}
392 \fig{DIFFUSION FIG 2} shows the result for some selected time steps.
393 An easy way to visualize the results is the command
394 \begin{verbatim}
395 mayavi -d T.1.xml -m SurfaceMap &
396 \end{verbatim}
397 Use the \texttt{Configure Data} window in mayavi
398 to move forward and and backwards in time.

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26