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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 625 - (show annotations)
Thu Mar 23 00:41:25 2006 UTC (13 years, 6 months ago) by gross
File MIME type: application/x-tex
File size: 18678 byte(s)
some updates and linearPDE.tex switched off
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 % http://www.opensource.org/licenses/osl-3.0.php
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, see~\cite{XXX} for alternatives. 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)[0])
218 print mypde.checkSymmetry() # returns False
219 mypde.setValue(C=kronecker(mydomain)[0])
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[0], \
257 d=eta,y=kappa*n[0]+eta*x[0])
258 #... calculate error of the PDE solution ...
259 u=mypde.getSolution()
260 print "error is ",Lsup(u-x[0])
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<t_end:
289 mypde.setValue(Y=q+rhocp/h*T)
290 T=mypde.getSolution()
291 t+=h
292 \end{python}
293 \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
294 in the domain and \var{h} is the time step size.
295 The variable \var{T}
296 holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
297 Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
298 temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
299 initial temperature distribution, which equal to $T\hackscore{ref}$.
300 The \LinearPDE class object \var{mypde}
301 and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient
302 $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information
303 from previous time steps as much as possible.
304
305 The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
306 in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
307 \var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired:
308 \begin{python}
309 from esys.escript import length,whereNegative
310 xc=[0.02,0.002]
311 r=0.001
312 x=mydomain.getX()
313 qH=q0*whereNegative(length(x-xc)-r)
314 \end{python}
315 \var{x} is a \Data class object of
316 the \escript module defining locations in the \Domain \var{mydomain}.
317 The \function{length()} imported from the \escript module returns the
318 Euclidean norm:
319 \begin{equation}\label{DIFFUSION HELM EQ 3aba}
320 \|y\|=\sqrt{
321 y\hackscore{i}
322 y\hackscore{i}
323 } = \function{esys.escript.length}(y)
324 \end{equation}
325 So \code{length(x-xc)} calculates the distances
326 of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
327 Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
328 converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative}
329 applied to
330 \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
331 zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
332
333 Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
334 \index{scripts!\file{diffusion.py}}:
335 \begin{python}
336 from esys.escript import *
337 from esys.escript.linearPDEs import LinearPDE
338 from esys.finley import Rectangle
339 #... set some parameters ...
340 xc=[0.02,0.002]
341 r=0.001
342 qc=50.e6
343 Tref=0.
344 rhocp=2.6e6
345 eta=75.
346 kappa=240.
347 tend=5.
348 # ... time, time step size and counter ...
349 t=0
350 h=0.1
351 i=0
352 #... generate domain ...
353 mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
354 #... open PDE ...
355 mypde=LinearPDE(mydomain)
356 mypde.setSymmetryOn()
357 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
358 # ... set heat source: ....
359 x=mydomain.getX()
360 qH=qc*whereNegative(length(x-xc)-r)
361 # ... set initial temperature ....
362 T=Tref
363 # ... start iteration:
364 while t<tend:
365 i+=1
366 t+=h
367 print "time step :",t
368 mypde.setValue(Y=qH+rhocp/h*T)
369 T=mypde.getSolution()
370 saveVTK("T.%d.xml"%i,temp=T)
371 \end{python}
372 The script will create the files \file{T.1.xml},
373 \file{T.2.xml}, $\ldots$, \file{T.50.xml} in the directory where the script has been started. The files give the
374 temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \VTK file format.
375
376 \begin{figure}
377 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes1.eps}}
378 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes16.eps}}
379 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes32.eps}}
380 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes48.eps}}
381 \caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
382 \label{DIFFUSION FIG 2}
383 \end{figure}
384 An easy way to visualize the results is the command
385 \begin{verbatim}
386 mayavi -d T.1.xml -m SurfaceMap &
387 \end{verbatim}
388 Use the \texttt{Configure Data}
389 to move forward and and backwards in time.
390 \fig{DIFFUSION FIG 2} shows the result for some selected time steps.

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26