ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (17 years, 6 months ago) by jgs
File MIME type: application/x-tex
File size: 17701 byte(s)
*** empty log message ***

1 jgs 102 % $Id$
2     \chapter{How to Solve The Diffusion Equation}
3     \label{DIFFUSION CHAP}
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     \begin{figure}
12     \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
13     \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
14     \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
15     \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
16     \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
17     \label{DIFFUSION FIG 2}
18     \end{figure}
21     \section{\label{DIFFUSION OUT SEC}Outline}
22     In this chapter we will discuss how to solve the time depeneded-temperature diffusion\index{diffusion equation} within
23     a block of material. Within the block there is a heat source which drives the temperature diffusion.
24     On the surface energy can radiate into the surrounding environment.
25     \fig{DIFFUSION FIG 1} shows the configuration.
27     In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A
28     time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$.
29     We will see that at time step a so-called Helmholtz equation \index{Helmholtz equation} has to be solved.
30     The implementation iof a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}.
31     In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a
32     solver for the temperature diffusion problem.
34     \section{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
36     The temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
37     in the interior of the domain is given by
38     \begin{equation}
39     \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q
40     \label{DIFFUSION TEMP EQ 1}
41     \end{equation}
42     where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite
43     material the parameters are depending on the their location in the domain. $q$ is
44     a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}
45     as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant $q^{c}$
46     on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:
47     \begin{equation}
48     q(x,t)=
49     \left\{
50     \begin{array}{lcl}
51     q^c & & \|x-x^c\| \le r \\
52     & \mbox{if} \\
53     0 & & \mbox{else} \\
54     \end{array}
55     \right.
56     \label{DIFFUSION TEMP EQ 1b}
57     \end{equation}
58     for all $x$ in the domain and all time $t>0$.
60     On the surface of the domain we are
61     are specifying a radiation condition
62     which precribes the normal component of the flux $J\hackscore i= \kappa T\hackscore{,i}$ to be proportional
63     to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$:
64     \begin{equation}
65     \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T)
66     \label{DIFFUSION TEMP EQ 2}
67     \end{equation}
68     $\eta$ is a given material coefficient depending on the material and the surrounding medium.
69     As usual $n_i$ is the $i$-th component of the outer normal field \index{outer normal field}
70     at the surface of the domain.
72     To solve the the time depended \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time
73     $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature:
74     \begin{equation}
75     T(x,0)=T\hackscore{ref}
76     \label{DIFFUSION TEMP EQ 4}
77     \end{equation}
78     for all $x$ in the domain. It is pointed out that
79     the initial conditions is fullfilling the
80     boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.
82     The temperature is calculated discrete time nodes $t^{(n)}$ where
83     $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.
84     In the following the upper index ${(n)}$ is refering to a value at time $t^{(n)}$. The simplest
85     and most robust scheme to approximate the time derivative of the the temperature is the
86     \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler scheme bases
87     on the Taylor expansion
88     \begin{equation}
89     T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)})
90     =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)}
91     \label{DIFFUSION TEMP EQ 6}
92     \end{equation}
93     which 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 + \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 from \eqn{DIFFUSION TEMP EQ 2}.
101     this forms a boundary value problem that has to be solved for each time step.
102     As a first step to implement a solver for the temperature diffusion problem we will
103     first implement a solver for the boundary value problem that has to be solved at each time step.
105     \section{\label{DIFFUSION HELM SEC}Helmholtz Problem}
106     The partial differential equation to be solved for $T^{(n)}$ has the form
107     \begin{equation}
108     \omega u - (\kappa u\hackscore{,i})\hackscore{,i} = f
109     \label{DIFFUSION HELM EQ 1}
110     \end{equation}
111     where $u$ plays the role of $T^{(n)}$ and we set
112     \begin{equation}
113     \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.
114     \label{DIFFUSION HELM EQ 1b}
115     \end{equation}
116     With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2}
117     takes the form
118     \begin{equation}
119     \kappa u\hackscore{,i} n\hackscore{i} = g - \eta u\mbox{ on } \Gamma
120     \label{DIFFUSION HELM EQ 2}
121     \end{equation}
122     The partial differential
123     \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
124     is called a Helmholtz equation \index{Helmholtz equation}.
126     We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the
127     Helmholtz equation. We have used a special case of the \LinearPDE class, namely the
128     \Poisson class already in \Chap{FirstSteps}.
129     Here we will write our own specialized class of the \LinearPDE to solve the Helmholtz equation.
131     The general form of a single PDE that can be handeled by the \LinearPDE class is
132     \begin{equation}\label{EQU.FEM.1}
133     -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .
134     \end{equation}
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 inpecting 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 like the \LinearPDE class but
162     is defined through the coefficeints $\kappa$, $\omega$, $f$, $\eta$,
163     $g$ rather than the general form given by \eqn{EQU.FEM.1}.
164     Python's
165     mechanism of inhertence allows doing this in a very easy way.
166     The advantage is that our new \class{Helmholtz} can be used in any context
167     that works with a \LinearPDE class but with an easier interface to define the PDE.
168     This improves reuasablity as well as maintainability of program codes.
170     We want to implement a
171     new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but
172     is defined through the coefficeints $\kappa$, $\omega$, $f$, $\alpha$,
173     $g$ rather than the general form given by \eqn{EQU.FEM.1}.
174     Python's mechanism of subclasses allows doing this in a very easy way.
175     The \Poisson class of the \linearPDEsPack module,
176     which we have already used in \Chap{FirstSteps}, is in fact a subclass of the
177     \LinearPDE class. That means that it all methods (such as the \method{getSolution})
178     from the parent class \LinearPDE are defined for any \Poisson object. However, new
179     methods can be added and methods of the parent class can be redefined. In fact,
180     the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
181     values to the coefficients of the PDE. This is exactly what we will do when we define
182     our new \class{Helmholtz} class:
183     \begin{python}
184     from esys.linearPDEs import LinearPDE
185     import numarray
186     class Helmholtz(LinearPDE)
187     def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
188     self._setValue(A=kappa*numarray.identity(self.getDim()),D=omega,Y=f,d=eta,y=g)
189     \end{python}
190     \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
191     of the \LinearPDE which have imported in the first line of the script.
192     We add the method \method{setValue} to the class which overwrites the
193     \method{setValue} method of the \LinearPDE class. The new methods which has the
194     parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
195     maps the parameters of the coefficients of the general PDE defined
196     in \eqn{EQU.FEM.1}. The coefficient \var{A} is defined by Kroneckers symbol. we use the
197     \numarray function \function{identity} which return a square matrix which has ones on the
198     main diagonal and zeros out side the main diagonal. The argument of \function{identity} gives the order of the matrix.
199     Here we use
200     the \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimension of the domain of the
201     PDE. As we will make use of the \class{Helmholtz} class several times, it is convient to
202     put is definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.
203     You can use your favourite editor to create and edit the file.
205     An object of the \class{Helmholtz} class is created through the statments:
206     \begin{python}
207     from mytools import *
208     mypde=Helmholtz(mydomain)
209     mypde.setValue(kappa=10.,omega=0.1,f=12)
210     u=mypde.getSolution()
211     \end{python}
212     In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure
213     that \file{mytools.py} is in directory from where you have started Python.
214     \var{mydomain} is the \Domain of the PDE. In the third statment values are
215     assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
216     specified the default values $0$ are used \footnote{It would be better to use the default value
217     \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
218     would not be processed when the PDE is evaluated.}. In the forth statement the solution of the
219     PDE is returned.
221     We want to test our \class{Helmholtz} class on a rectangular domain
222     of length $l$ and height $h$. We do this by choosing a simple test solution,
223     here we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that
224     the test solution becomes the solution of the problem. If we take $\kappa=1$
225     an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
226     $\kappa n\hackscore{i} u\hackscore{,i}=n\hackscore{0}$.
227     So we have to set $g=n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}
228     \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory
229     implements this test problem using the \finley PDE solver:
230     \begin{python}
231     from mytools import *
232     from esys.escript import *
233     import esys.finley
234     #... set some parameters ...
235     omega=0.1
236     eta=10.
237     #... generate domain ...
238     mydomain = esys.finley.Rectangle(l0=5.,l1=1.,n0=50, n1=10)
239     #... open PDE and set coefficients ...
240     mypde=Helmholtz(mydomain)
241     n=mydomain.getNormal()
242     x=mydomain.getX()
243     mypde.setValue(1,omega,omega*x[0],eta,n[0]+eta*x[0])
244     #... calculate error of the PDE solution ...
245     u=mypde.getSolution()
246     print "error is ",Lsup(u-x[0])
247     \end{python}
248     The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.
249     \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
250     is imported by the \code{from escript import *} statement and returns the maximum absulute value of it argument. To run
251     The error shown by the print statement should be in the order of $10^{-7}$. As piecewise linear interpolation is
252     used to approximate the solution and our solution is a linear function of the spatial coordinates one may
253     expect that the error is zero. However, as most PDE packages uses an iterative solver which is terminated
254     when a given tolerance has been reached. The default tolerance is $10^{-8}$. Thsi value can be altered by using the
255     \method{setTolerance} of the \LinearPDE class.
257     \section{The Transition Problem}
258     \label{DIFFUSION TRANS SEC}
259     Now we are ready to solve the original time dependent problem. The main
260     part of the script is the loop over the time $t$ which takes the following form:
261     \begin{python}
262     mypde=Helmholtz(mydomain)
263     while t<t_end:
264     mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)
265     T=mypde.getSolution()
266     t+=h
267     \end{python}
268     \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
269     in the domain and \var{h} is the time step size which has to be chosen. Notice that the \class{Hemholtz}
270     is created before the loop over time is entered while in each time step only the coefficients
271     are reset in each time step. This way some information about the reperesentation of the PDE can be reused
272     \footnote{The efficience can be improved further by setting the coefficients in the operator
273     \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only update the coefficients
274     in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
275     method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator}. The variable \var{T}
276     holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
277     Helmholtz \eqn{DIFFUSION HELM EQ 1}. Statement \code{T=mypde.getSolution()} overwrites \var{T} with the
278     temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to sepcify the
279     initial temperature distribution, which equal to $T\hackscore{ref}$.
281     The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} shall be \var{q0}
282     at an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
283     \var{q0} is a fixed constant. The following script defines \var{q} as desired:
284     \begin{python}
285     xc=[0.02,0.002]
286     r=0.001
287     x=mydomain.getX()
288     q=q0*(length(x-xc)-r).whereNegative()
289     \end{python}
290     \var{x} is a \Data class object of
291     the \escript module defining the locations of points in the \Domain \var{mydomain}.
292     \code{length(x-xc)} calculates the distances in the Euclidean norm
293     of the locations \var{x} to the center of the circle \var{xc} where the heat source is acting.
294     Notice that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
295     converted into a \Data class object before subtracted from \var{x}. The method \method{whereNegative} of
296     a \Data class object, in this case the result of the expression
297     \code{length(x-xc)-r}, returns a \Data class which is equal one where the object negative and
298     zero elsewhere. After multiplication with \var{q0} we get a function with the deired property.
300     Now we can put the components together to the script \file{diffusion.py} which is available in the \ExampleDirectory:
301     \index{scripts!\file{diffusion.py}}:
302     \begin{python}
303     from mytools import *
304     from esys.escript import *
305     import esys.finley
306     #... set some parameters ...
307     x_c=[0.02,0.002]
308     r=0.001
309     q0=50.e6
310     Tref=0.
311     rhocp=2.6e6
312     eta=75.
313     kappa=240.
314     t_end=5.
315     # ...time step size and counter ...
316     h=0.1
317     i=0
318     t=0
319     #... generate domain ...
320     mydomain = esys.finley.Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
321     #... open PDE ...
322     mypde=Helmholtz(mydomain)
323     # ... set heat source: ....
324     x=mydomain.getX()
325     q=q0*(length(x-x_c)-r).whereNegative()
326     # ... set initial temperature ....
327     T=Tref
328     # ... start iteration:
329     while t<t_end:
330     i+=1
331     t+=h
332     print "time step :",t
333     mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
334     T=mypde.getSolution()
335     T.saveDX("T%d.dx"%i)
336     \end{python}
337     The script will create the files \file{T.1.dx},
338     \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the
339     temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.
340     An easy way to visualize the results is the command
341     \begin{verbatim}
342     dx -edit diffusion.net
343     \end{verbatim}
344     where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.
345     \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