ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 3306 - (show annotations)
Mon Oct 25 05:09:13 2010 UTC (10 years, 10 months ago) by caltinay
File MIME type: application/x-tex
File size: 18295 byte(s)
Commented declaremodule and modulesynopsis.

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


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

  ViewVC Help
Powered by ViewVC 1.1.26