ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

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