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

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

Parent Directory Parent Directory | Revision Log Revision Log


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