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

Annotation of /trunk/doc/user/heatedblock.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 568 - (hide annotations)
Tue Feb 28 03:59:06 2006 UTC (13 years, 8 months ago) by gross
Original Path: trunk/doc/user/beam.tex
File MIME type: application/x-tex
File size: 24299 byte(s)
update
1 jgs 110 % $Id$
2 gross 568
3     The \LinearPDE class is used to define a general linear, steady, second order PDE
4     for an unknown function $u$ on a given $\Omega$ defined through a \Domain object.
5     In the following $\Gamma$ denotes the boundary of the domain $\Omega$. $n$ denotes
6     the outer normal field on $\Gamma$.
7    
8     For a single PDE with a solution with a single component the linear PDE is defined in the
9     following form:
10     \begin{equation}\label{LINEARPDE.SINGLE.1}
11     -(A\hackscore{jl} u\hackscore{,l}){,j}+(B\hackscore{j} u)\hackscore{,j}+C\hackscore{l} u\hackscore{,l}+D u =-X\hackscore{j,j}+Y \; .
12     \end{equation}
13     $u_{,j}$ denotes the derivative of $u$ with respect to the $j$-th spatial direction. Einstein's summation convention, ie. summation over indexes appearing twice in a term of a sum is performed, is used.
14     The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the
15     \Function on the PDE or objects that can be converted into such \Data objects.
16     $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar.
17     The following natural
18     boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
19     \begin{equation}\label{LINEARPDE.SINGLE.2}
20     n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y \;.
21     \end{equation}
22     Notice that the coefficients $A$, $B$ and $X$ are defined in the PDE. The coefficients $d$ and $y$ are
23     each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the
24     solution at certain locations in the domain. They have the form
25     \begin{equation}\label{LINEARPDE.SINGLE.3}
26     u=r \mbox{ where } q>0
27     \end{equation}
28     $r$ and $q$ are each \Scalar where $q$ is the characteristic function
29     \index{characteristic function} defining where the constraint is applied.
30     The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1}
31     or \eqn{LINEARPDE.SINGLE.2}. The PDE is symmetrical \index{symmetrical} if
32     \begin{equation}\label{LINEARPDE.SINGLE.4}
33     A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j}
34     \end{equation}
35     For a system of PDEs and a solution with several components the PDE has the form
36     \begin{equation}\label{LINEARPDE.SYSTEM.1}
37     -(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u_k)\hackscore{,j}+C\hackscore{ikl} u\hackscore{k,l}+D\hackscore{ik} u_k =-X\hackscore{ij,j}+Y\hackscore{i} \; .
38     \end{equation}
39     $A$ is a \RankFour, $B$ and $C$ are each a \RankThree, $D$ and $X$ are each a \RankTwo and $Y$ is a \RankOne.
40     The natural boundary conditions \index{boundary condition!natural} take the form:
41     \begin{equation}\label{LINEARPDE.SYSTEM.2}
42     n\hackscore{j}(A\hackscore{ijkl} u\hackscore{k,l}){,j}+(B\hackscore{ijk} u_k)+d\hackscore{ik} u_k=n\hackscore{j}-X\hackscore{ij}+y\hackscore{i} \;.
43     \end{equation}
44     The coefficient $d$ is a \RankTwo and $y$ is a
45     \RankOne both in the \FunctionOnBoundary. Constraints \index{constraint} take the form
46     \begin{equation}\label{LINEARPDE.SYSTEM.3}
47     u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0
48     \end{equation}
49     $r$ and $q$ are each \RankOne. Notice that at some locations not necessarily all components must
50     have a constraint. The system of PDEs is symmetrical \index{symmetrical} if
51     \begin{eqnarray}\label{LINEARPDE.SYSTEM.4}
52     A\hackscore{ijkl}=A\hackscore{klij} \\
53     B\hackscore{ijk}=C\hackscore{kij} \\
54     D\hackscore{ik}=D\hackscore{ki} \\
55     d\hackscore{ik}=d\hackscore{ki} \
56     \end{eqnarray}
57     \LinearPDE also supports solution discontinuities \index{discontinuity} over contact region $\Gamma^{contact}$
58     in the domain $\Omega$. To specify the conditions across the discontinuity we are using the
59     generalised flux $J$ which is in the case of a systems of PDEs and several components of the solution
60     defined as
61     \begin{equation}\label{LINEARPDE.SYSTEM.5}
62     J\hackscore{ij}=A\hackscore{ijkl}u\hackscore{k,l}+B\hackscore{ijk}u\hackscore{k}-X\hackscore{ij}
63     \end{equation}
64     For the case of single solution component and single PDE $J$ is defined
65     \begin{equation}\label{LINEARPDE.SINGLE.5}
66     J\hackscore{j}=A\hackscore{jl}u\hackscore{,l}+B\hackscore{j}u\hackscore{k}-X\hackscore{j}
67     \end{equation}
68     In the context of discontinuities \index{discontinuity} $n$ denotes the normal on the
69     discontinuity pointing from side 0 towards side 1. For a system of PDEs
70     the contact condition takes the form
71     \begin{equation}\label{LINEARPDE.SYSTEM.6}
72     n\hackscore{j} J^{0}\hackscore{ij}=n\hackscore{j} J^{1}\hackscore{ij}=y^{contact}\hackscore{i} - d^{contact}\hackscore{ik} [u]\hackscore{k} \; .
73     \end{equation}
74     where $J^{0}$ and $J^{1}$ are the fluxes on side $0$ and side $1$ of the
75     discontinuity $\Gamma^{contact}$, respectively. $[u]$, which is the difference
76     of the solution at side 1 and at side 0, denotes the jump of $u$ across $\Gamma^{contact}$.
77     The coefficient $d^{contact}$ is a \RankTwo and $y^{contact}$ is a
78     \RankOne both in the \FunctionOnContactZero or \FunctionOnContactOne.
79     In case of a single PDE and a single component solution the contact condition takes the form
80     \begin{equation}\label{LINEARPDE.SINGLE.6}
81     n\hackscore{j} J^{0}\hackscore{j}=n\hackscore{j} J^{1}\hackscore{j}=y^{contact} - d^{contact}[u]
82     \end{equation}
83     In this case the the coefficient $d^{contact}$ and $y^{contact}$ are eaach \Scalar
84     both in the \FunctionOnContactZero or \FunctionOnContactOne.
85     ======================
86    
87    
88     We have used a special case of the \LinearPDE class, namely the
89     \Poisson class already in \Chap{FirstSteps}.
90     Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
91     and use the \method{getSolution} method of parent class to actually solve the problem.
92    
93     The form of a single PDE that can be handled by the \LinearPDE class is
94     \begin{equation}\label{EQU.FEM.1}
95     -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .
96     \end{equation}
97     We show here the terms which are relevant for the Helmholtz problem.
98     The general form and systems is discussed in \Sec{SEC LinearPDE}.
99     $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients}
100     Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar.
101     They may be constant or may depend on their
102     location in the domain but must not depend on the unknown solution $u$.
103     The following natural boundary conditions \index{boundary condition!natural} that
104     are used in the \LinearPDE class have the form
105     \begin{equation}\label{EQU.FEM.2}
106     n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y \;.
107     \end{equation}
108     where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that
109     the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients.
110    
111     By inspecting the Helmholtz equation \index{Helmholtz equation}
112     we can easily assign values to the coefficients in the
113     general PDE of the \LinearPDE class:
114     \begin{equation}\label{DIFFUSION HELM EQ 3}
115     \begin{array}{llllll}
116     A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\
117     d=\eta & y= g & \\
118     \end{array}
119     \end{equation}
120     $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
121     $i=j$ and $0$ otherwise.
122    
123     We want to implement a
124     new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
125     is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,
126     $g$ rather than the general form given by \eqn{EQU.FEM.1}.
127     Python's mechanism of subclasses allows us to do this in a very simple way.
128     The \Poisson class of the \linearPDEs module,
129     which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
130     \LinearPDE class. That means that all methods (such as the \method{getSolution})
131     from the parent class \LinearPDE are available for any \Poisson object. However, new
132     methods can be added and methods of the parent class can be redefined. In fact,
133     the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
134     values to the coefficients of the PDE. This is exactly what we will do when we define
135     our new \class{Helmholtz} class:
136     \begin{python}
137     from esys.linearPDEs import LinearPDE
138     import numarray
139     class Helmholtz(LinearPDE):
140     def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
141     ndim=self.getDim()
142     kronecker=numarray.identity(ndim)
143     self.setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)
144     \end{python}
145     \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
146     of \LinearPDE which we have imported in the first line of the script.
147     We add the method \method{setValue} to the class which overwrites the
148     \method{setValue} method of the \LinearPDE class. The new method which has the
149     parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
150     maps the parameters of the coefficients of the general PDE defined
151     in \eqn{EQU.FEM.1}. We are actually using the \method{_LinearPDE__setValue} of
152     the \LinearPDE class.
153     The coefficient \var{A} involves the Kronecker symbol. We use the
154     \numarray function \function{identity} which returns a square matrix with ones on the
155     main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix.
156     The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the
157     PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to
158     put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.
159     You can use your favourite editor to create and edit the file.
160    
161     An object of the \class{Helmholtz} class is created through the statments:
162     \begin{python}
163     from mytools import Helmholtz
164     mypde = Helmholtz(mydomain)
165     mypde.setValue(kappa=10.,omega=0.1,f=12)
166     u = mypde.getSolution()
167     \end{python}
168     In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure
169     that \file{mytools.py} is in the directory from where you started Python.
170     \var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are
171     assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
172     specified, the default values $0$ are used. \footnote{It would be better to use the default value
173     \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
174     would not be processed when the PDE is evaluated}. In the fourth statement the solution of the
175     PDE is returned.
176    
177     To test our \class{Helmholtz} class on a rectangular domain
178     of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
179     we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that
180     the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
181     an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
182     $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
183     So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}
184     \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory
185     implements this test problem using the \finley PDE solver:
186     \begin{python}
187     from mytools import Helmholtz
188     from esys.escript import Lsup
189     from esys.finley import Rectangle
190     #... set some parameters ...
191     kappa=1.
192     omega=0.1
193     eta=10.
194     #... generate domain ...
195     mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
196     #... open PDE and set coefficients ...
197     mypde=Helmholtz(mydomain)
198     n=mydomain.getNormal()
199     x=mydomain.getX()
200     mypde.setValue(kappa,omega,omega*x[0],eta,kappa*n[0]+eta*x[0])
201     #... calculate error of the PDE solution ...
202     u=mypde.getSolution()
203     print "error is ",Lsup(u-x[0])
204     \end{python}
205     The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.
206     \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
207     is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument.
208     The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
209     used to approximate the solution and our solution is a linear function of the spatial coordinates one might
210     expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
211     However most PDE packages use an iterative solver which is terminated
212     when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
213     \method{setTolerance} of the \LinearPDE class.
214    
215     \subsection{The Transition Problem}
216     \label{DIFFUSION TRANS SEC}
217     Now we are ready to solve the original time dependent problem. The main
218     part of the script is the loop over time $t$ which takes the following form:
219     \begin{python}
220     t=0
221     T=Tref
222     mypde=Helmholtz(mydomain)
223     while t<t_end:
224     mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)
225     T=mypde.getSolution()
226     t+=h
227     \end{python}
228     \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
229     in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde}
230     is created before the loop over time is entered while in each time step only the coefficients
231     are reset in each time step. This way some information about the representation of the PDE can be reused
232     \footnote{The efficience can be improved further by setting the coefficients in the operator
233     \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only updating the coefficients
234     in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
235     method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator.}. The variable \var{T}
236     holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
237     Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
238     temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
239     initial temperature distribution, which equal to $T\hackscore{ref}$.
240    
241     The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
242     in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
243     \var{q0} is a fixed constant. The following script defines \var{q} as desired:
244     \begin{python}
245     from esys.escript import length
246     xc=[0.02,0.002]
247     r=0.001
248     x=mydomain.getX()
249     q=q0*(length(x-xc)-r).whereNegative()
250     \end{python}
251     \var{x} is a \Data class object of
252     the \escript module defining locations in the \Domain \var{mydomain}.
253     The \function{length()} imported from the \escript module returns the
254     Euclidean norm:
255     \begin{equation}\label{DIFFUSION HELM EQ 3aba}
256     \|y\|=\sqrt{
257     y\hackscore{i}
258     y\hackscore{i}
259     } = \function{esys.escript.length}(y)
260     \end{equation}
261     So \code{length(x-xc)} calculates the distances
262     of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
263     Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
264     converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of
265     a \Data class object, in this case the result of the expression
266     \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
267     zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
268    
269     Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
270     \index{scripts!\file{diffusion.py}}:
271     \begin{python}
272     from mytools import Helmholtz
273     from esys.escript import Lsup
274     from esys.finley import Rectangle
275     #... set some parameters ...
276     xc=[0.02,0.002]
277     r=0.001
278     qc=50.e6
279     Tref=0.
280     rhocp=2.6e6
281     eta=75.
282     kappa=240.
283     tend=5.
284     # ... time, time step size and counter ...
285     t=0
286     h=0.1
287     i=0
288     #... generate domain ...
289     mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
290     #... open PDE ...
291     mypde=Helmholtz(mydomain)
292     # ... set heat source: ....
293     x=mydomain.getX()
294     q=qc*(length(x-xc)-r).whereNegative()
295     # ... set initial temperature ....
296     T=Tref
297     # ... start iteration:
298     while t<tend:
299     i+=1
300     t+=h
301     print "time step :",t
302     mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
303     T=mypde.getSolution()
304     T.saveDX("T%d.dx"%i)
305     \end{python}
306     The script will create the files \file{T.1.dx},
307     \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the
308     temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.
309    
310     \begin{figure}
311     \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
312     \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
313     \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
314     \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
315     \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
316     \label{DIFFUSION FIG 2}
317     \end{figure}
318    
319     An easy way to visualize the results is the command
320     \begin{verbatim}
321     dx -edit diffusion.net &
322     \end{verbatim}
323     where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.
324     Use the \texttt{Sequencer} to move forward and and backwards in time.
325     \fig{DIFFUSION FIG 2} shows the result for some selected time steps.
326     ====
327 jgs 121 \section{Bending Beam under Uniform Load}
328 jgs 110 \label{BEAM CHAP}
329     \sectionauthor{Jannine Eisenmann}{}
330    
331     Given is a two-dimension bending beam which is fixed in all directions
332     at the left end and free at the other, see \fig{BEAM FIG 1}. Furthermore the beam is loaded
333     with a uniform load $p$.
334    
335     \begin{figure}
336     % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
337     \caption{Loaded Beam.}
338     \label{BEAM FIG 1}
339     \end{figure}
340    
341    
342    
343     For emphasizing the displacement this picture is drawn (uebertrieben),
344     cause since we use the linear theory otherwise no displacements would
345     be visible.
346     There are two ways of solving this problem: on the one hand one can
347     use the differential equation for the deformed neutral fibre of the
348     beam. This classical differential equation is based on several simplifying
349     assumptions concerning the statics and kinematics of the problem.
350     However the results are known to be highly accurate at least for slender
351     beams with length to hight ratios $> 10$. Alternatively, in connection
352     with finite element based differential equation toolkits one may simply
353     consider the beam as a special case of a 2D or 3D elastic continuum
354     problem and solve the stress equilibrium equations combined with Hooke's
355     law for the specific boundary conditions of a beam. Both cases assume
356     isotropic and linear elastic material properties.
357    
358     The beam equation is easily solved analytically. The analytical solutions
359     can be used for benchmarkung finite element solutions. In section
360     1.2 we formluate a finite element code for general isotropic elasticity
361     problems including thin and deep beams under arbitrary loading conditiond
362     in 2D or 3D.
363    
364    
365     \section{2-dimensional}
366     As the stress equilibrium equation is a partial differential equation
367     (PDE), we choose to use Finley to solve it, since Finley is a finite
368     element kernel library, that has been incorporated as a differential
369     equation solver into the python based numerical toolkit called escript.
370     We divided the beam into ten square elements of the size 1x1. Each
371     element consists of 8 nodes, which leads to a quadratic interpolation
372     of the model point displacements \\
373    
374     The key ingredients is \textbf{Hooks Law}. We use Hooks Law because
375     we are dealing with \textbf{linear elastic material} \textbf{behaviour}.
376     We have \\
377    
378    
379     $\sigma_{ik}=2G\left(\varepsilon_{ik}+\frac{\nu}{1-2\nu}\cdot e\cdot\delta_{ik}\right)$\hfill{}(1)\\
380     where the engineering strain$\varepsilon_{ik}$is defined as:
381    
382     $\varepsilon_{ik}=\frac{1}{2}\cdot\left(u_{k,i}+u_{i,k}\right)$\hfill{}(2)\\
383    
384    
385     with
386    
387     \begin{enumerate}
388     \item e= Volume strain = $\varepsilon_{kk}$
389     \item $\delta_{ik}$= Kronecker symbol
390     \end{enumerate}
391     Inserting equation (2) in (1) (and with further mathematical conversions)
392     leads to the following partial differential equation :\\
393    
394    
395     $\sigma_{ij}=\left[\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)+\lambda\left(\delta_{ij}\delta_{kl}\right)\right]u_{k,l}$\\
396    
397    
398     For tension equilibrium we require:\\
399    
400    
401     $\sigma_{ij,j}=0$\\
402    
403    
404     which leads to the following expression:\\
405    
406    
407     \[
408     \left(\left[\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)+\lambda\left(\delta_{ij}\delta_{kl}\right)\right]u_{k,l}\right)_{,j}=0\]
409    
410    
411     with the natural boundary condition:
412    
413     \[
414     n_{j}\sigma_{ij}=-p_{i}\]
415     \\
416     $p_{i}$ representing a uniform load on top of the beam.
417    
418     A Dirichlet Boundary condition is assumed on the left end of the beam
419     for which no displacements can occure.\\
420     \\
421     \includegraphics[%
422     width=0.60\linewidth,bb = 0 0 200 100, draft, type=eps]{/home/jeannine/sandbox/report/draws/dir_cond_beam.eps}\\
423     This is described in the code with the setting a mask for the left
424     end and setting values to that mask:
425    
426     \begin{python}
427     q = xNodes{[}0{]}.whereZero(){*}{[}1.0,1.0{]}
428    
429     r = Vector({[}0.0, 0.0{]}, where = nodes)
430     \end{python}
431     The Finley template PDE reads:\\
432    
433    
434     \[
435     -(A_{ijkl}u_{k,l})_{,j}-(B_{ijk}u_{k})_{,j}+C_{ikl}u_{k,l}+D_{ik}u_{k}=-X_{ij,j}+Y_{i}\]
436     \\
437     with the natural boundary condition:
438    
439     \[
440     n_{j}(A_{ijkl}u_{k,l}+B_{ijkl}u_{k})+d_{ik}u_{k}=n_{j}X_{ij}+y_{i}on\Gamma_{i}^{D}\]
441    
442    
443     Yields by comparing the coefficients :
444    
445     \begin{enumerate}
446     \item $A_{ijkl}$= $\left[\mu\left(\delta_{ik}\delta_{ij}+\delta_{jl}\delta_{il}\right)+\lambda\left(\delta_{ij}\delta_{kl}\right)\right]$
447     \item $y_{i}$= $-p_{i}$
448     \item $u_{k}$= displacement $u$
449     \end{enumerate}
450     $B_{ijk,}=C_{ikl}=D_{ik}=X_{ij}=Y_{i}=d_{ik}=0$\\
451    
452    
453     Where 0 in the last line is taken as a scalar, vector or tensor, depending
454     on what the belonging coefficient is.
455    
456     These equations are the base for the \textbf{Finley Script}:
457    
458     \begin{python}
459     from ESyS import {*}
460     import Finley
461    
462    
463    
464     \#mu lamda
465    
466     def mu (E, nu): \#= shear modul G
467    
468     return E/(2{*}(1+nu))
469    
470     def lamda (E, nu):
471    
472     return (nu{*}E)/((1-2{*}nu){*}(1+nu))
473    
474    
475    
476     def main()
477    
478     \#material, beam PROPERTIES
479    
480     L = 10.0 \#length of beam {[}m{]}
481    
482     h = 1 \#height of beam {[}m{]}
483    
484     p = 1 \#outer uniform load {[}kN/m{]}
485    
486     E0 = 210000 \#Young's modulus {[}kN/m\textasciicircum{}2{]}
487    
488     nu = 0.4 \#Poisson ratio {[}-{]}
489    
490    
491    
492     print \char`\"{}L=\char`\"{}, L
493    
494     print \char`\"{}h=\char`\"{}, h
495    
496     print \char`\"{}p=\char`\"{}, p
497    
498     print \char`\"{}E=\char`\"{}, E0
499    
500     print \char`\"{}nu=Poisson =\char`\"{}, nu
501    
502     print \char`\"{}mu=\char`\"{}, mu (E0,nu)
503    
504     print \char`\"{}lamda=\char`\"{}, lamda (E0,nu)
505    
506    
507    
508     \#SET MESH for FE
509    
510     mesh= Finley.Rectangle(n0=10 ,n1=1 ,order=2, l0=L, l1=h)
511    
512     nodes = mesh.Nodes()
513    
514     xNodes = nodes.getX()
515    
516     elements = mesh.Elements()
517    
518     faceElements = mesh.FaceElements()
519    
520     xFaceElements = faceElements.getX()
521    
522    
523    
524     \#DISPLACEMENT boundary
525    
526     q = xNodes{[}0{]}.whereZero(){*}{[}1.,1.0{]} \#setting the mask for r
527    
528     r = Vector({[}0.0, 0.0{]}, where = nodes) \#fixed end
529    
530    
531    
532     \#STRESS boundary
533    
534     ymask = xFaceElements{[}1{]}.whereEqualTo(h)
535    
536     y = Vector({[}0, -p{]}, where = faceElements)
537    
538     y = y{*}ymask
539    
540    
541    
542     \#Finley coeff.
543    
544     A = Tensor4(0, where = elements)
545    
546     for i in range (2) :
547    
548     for j in range (2) :
549    
550     A{[}i,i,j,j{]}+= lamda (E0,nu)
551    
552     A{[}j,i,j,i{]}+= mu (E0,nu)
553    
554     A{[}j,i,i,j{]}+= mu (E0,nu)
555    
556    
557    
558     M,b = mesh.assemble(A=A, B=B, q=q,
559    
560     y=y,r=r,type=CSR,num\_equations=2)
561    
562     u= M.iterative(b, tolerance=1e-8,iter\_max=50000,
563    
564     iterative\_method=PCG)
565    
566     print \char`\"{}u{[}0{]}=\char`\"{},u{[}0{]}
567    
568     print \char`\"{}u{[}1{]}=\char`\"{},u{[}1{]}
569    
570     main()
571     \end{python}
572     The finer the mesh the more exact is the solution. E.g. a 20x2 mesh
573     is more exact than a 10x1 mesh.
574    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26