# Contents of /trunk/doc/user/beam.tex

Revision 568 - (show annotations)
Tue Feb 28 03:59:06 2006 UTC (15 years, 6 months ago) by gross
File MIME type: application/x-tex
File size: 24299 byte(s)
update

 1 % $Id$ 2 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 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