# Diff of /trunk/doc/user/heatedblock.tex

trunk/esys2/doc/user/beam.tex revision 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/doc/user/heatedblock.tex revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC
# Line 1  Line 1
1  % $Id$  % $Id$
2  \section{Bending Beam under Uniform Load}  %
4  \sectionauthor{Jannine Eisenmann}{}  %               \url{http://www.access.edu.au
5    %         Primary Business: Queensland, Australia.
6  Given is a two-dimension bending beam which is fixed in all directions  %   Licensed under the Open Software License version 3.0
7  at the left end and free at the other, see \fig{BEAM FIG 1}. Furthermore the beam is loaded  %      http://www.opensource.org/licenses/osl-3.0.php
8  with a uniform load $p$.  %
9    \section{Elastic Deformation}
10  \begin{figure}  \label{ELASTIC CHAP}
11  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}  In this section we want to discuss the deformation of a linear elastic body caused by expansion through a heat distribution. We want
12  \caption{Loaded Beam.}  to displacement field $u\hackscore{i}$ which solves the momentum equation
13  \label{BEAM FIG 1}  \index{momentum equation}:
14  \end{figure}  \begin{eqnarray}\label{HEATEDBLOCK general problem}
15     - \sigma\hackscore{ij,j}=0
16    \end{eqnarray}
17    where the stress $\sigma$ is given by
18  For emphasizing the displacement this picture is drawn (uebertrieben),  \begin{eqnarray}\label{HEATEDBLOCK linear elastic}
19  cause since we use the linear theory otherwise no displacements would   \sigma\hackscore{ij}= \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
20  be visible.   - (\lambda+\frac{2}{3} \mu)  \; \alpha  \;  (T-T\hackscore{ref})\delta\hackscore{ij} \;.
21  There are two ways of solving this problem: on the one hand one can  \end{eqnarray}
22  use the differential equation for the deformed neutral fibre of the  In this formula $\lambda$ and $\mu$ are the Lame coefficients, $\alpha$ is the
23  beam. This classical differential equation is based on several simplifying  temperature expansion coefficient, $T$ is the temperature distribution and $T_{ref}$ a reference temperature. Note that
24  assumptions concerning the statics and kinematics of the problem.  \eqn{HEATEDBLOCK general problem} is similar to eqn{WAVE general problem} introduced in section~\Sec{WAVE CHAP} but the
25  However the results are known to be highly accurate at least for slender  inertia term $\rho u\hackscore{i,tt}$ has been dropped as we assume a static scenario here. Moreover, in
26  beams with length to hight ratios $> 10$. Alternatively, in connection  comparison to the \eqn{WAVE stress}
27  with finite element based differential equation toolkits one may simply  definition of stress $\sigma$ in \eqn{HEATEDBLOCK linear elastic} an extra term is introduced
28  consider the beam as a special case of a 2D or 3D elastic continuum  to bring in stress due to volume changes trough temperature dependent expansion.
29  problem and solve the stress equilibrium equations combined with Hooke's
30  law for the specific boundary conditions of a beam. Both cases assume  Our domain is the unit cube
31  isotropic and linear elastic material properties.  \begin{eqnarray} \label{HEATEDBLOCK natural location}
32    \Omega=\{(x\hackscore{i} | 0 \le x\hackscore{i} \le 1 \}
33  The beam equation is easily solved analytically. The analytical solutions  \end{eqnarray}
34  can be used for benchmarkung finite element solutions. In section  On the boundary the normal stress component is set to zero
35  1.2 we formluate a finite element code for general isotropic elasticity  \begin{eqnarray} \label{HEATEDBLOCK natural}
36  problems including thin and deep beams under arbitrary loading conditiond  \sigma\hackscore{ij}n\hackscore{j}=0
37  in 2D or 3D.  \end{eqnarray}
38    and on the face with $x\hackscore{i}=0$ we set the $i$-th component of the displacement to $0$
39    \begin{eqnarray} \label{HEATEDBLOCK constraint}
40  \section{2-dimensional}  u\hackscore{i}(x)=0 & \mbox{ where } & x\hackscore{i}=0 \;
41  As the stress equilibrium equation is a partial differential equation  \end{eqnarray}
42  (PDE), we choose to use Finley to solve it, since Finley is a finite  For the temperature distribution we use
43  element kernel library, that has been incorporated as a differential  \begin{eqnarray} \label{HEATEDBLOCK temperature}
44  equation solver into the python based numerical toolkit called escript.  T(x)= T\hackscore{0} e^{-\beta \|x-x^{c}\|};
45  We divided the beam into ten square elements of the size 1x1. Each  \end{eqnarray}
46  element consists of 8 nodes, which leads to a quadratic interpolation  with a given positive constant $\beta$ and location $x^{c}$ in the domain. Later in \Sec{MODELFRAME} we will use
47  of the model point displacements \\  $T$ from an time-dependent temperature diffusion problem as discussed in \Sec{DIFFUSION CHAP}.
48
49  The key ingredients is \textbf{Hooks Law}. We use Hooks Law because  When we insert~\eqn{HEATEDBLOCK linear elastic} we get a second oder system of linear PDEs for the displacements $u$ which is called
50  we are dealing with \textbf{linear elastic material} \textbf{behaviour}.  the Lame equation\index{Lame equation}. We want to solve
51  We have \\  this using the \LinearPDE class to this. For a system of PDEs and a solution with several components the \LinearPDE class
52    takes PDEs of the form
53    \label{LINEARPDE.SYSTEM.1 TUTORIAL}
54  $\sigma_{ik}=2G\left(\varepsilon_{ik}+\frac{\nu}{1-2\nu}\cdot e\cdot\delta_{ik}\right)$\hfill{}(1)\\  -A\hackscore{ijkl} u\hackscore{k,l}){,j}=-X\hackscore{ij,j} \; .
55  where the engineering strain$\varepsilon_{ik}$is defined as:
56    $A$ is of \RankFour and $X$ is of \RankTwo. We show here the coefficients relevant
57  $\varepsilon_{ik}=\frac{1}{2}\cdot\left(u_{k,i}+u_{i,k}\right)$\hfill{}(2)\\  for the we traying to solve. The full form is given in~\eqn{LINEARPDE.SYSTEM.1}.
58    The natural boundary conditions \index{boundary condition!natural} take the form:
59    \label{LINEARPDE.SYSTEM.2 TUTORIAL}
60  with  n\hackscore{j} A\hackscore{ijkl} u\hackscore{k,l}=n\hackscore{j}X\hackscore{ij} \;.
61
62  \begin{enumerate}  Constraints \index{constraint} take the form
63  \item e= Volume strain = $\varepsilon_{kk}$  \label{LINEARPDE.SYSTEM.3 TUTORIAL}
64  \item $\delta_{ik}$= Kronecker symbol  u\hackscore{i}=r\hackscore{i} \mbox{ where } q\hackscore{i}>0
65  \end{enumerate}
66  Inserting equation (2) in (1) (and with further mathematical conversions)  $r$ and $q$ are each \RankOne.
67  leads to the following partial differential equation :\\  We can easily identify the coefficients in~\eqn{LINEARPDE.SYSTEM.1 TUTORIAL}:
68    \begin{eqnarray}\label{LINEARPDE ELASTIC COEFFICIENTS}
69    A\hackscore{ijkl}=\lambda \delta\hackscore{ij} \delta\hackscore{kl} + \mu (
70  $\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}$\\  +\delta\hackscore{ik} \delta\hackscore{jl}
71    \delta\hackscore{il} \delta\hackscore{jk}) \\
72    X\hackscore{ij}=(\lambda+\frac{2}{3} \mu) \;  \alpha \; (T-T\hackscore{ref})\delta\hackscore{ij} \\
73  For tension equilibrium we require:\\  \end{eqnarray}
74    The characteristic function $q$ defining the locations and components where constraints are set is given by:
76  $\sigma_{ij,j}=0$\\  q\hackscore{i}(x)=\left\{
77    \begin{array}{cl}
78    1  & x\hackscore{i}=0  \\
79  which leads to the following expression:\\  0  & \mbox{otherwise}   \\
80    \end{array}
81    \right.
82  $83 \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$  Under the assumption that $\lambda$, $\mu$, $\beta$ and $T\hackscore{ref}$
84    are constant setting $Y\hackscore{i}=\lambda+\frac{2}{3} \mu) \; \alpha \; T\hackscore{i}$ seems to be also possible. However,
85    this choice would lead to a different natural boundary condition which does not set the normal stress component as defined
86  with the natural boundary condition:  in~\eqn{HEATEDBLOCK linear elastic} to zero.
87
88  $Analogously to concept of symmetry for a single PDE, we call the PDE defined by~\eqn{LINEARPDE.SYSTEM.1 TUTORIAL} symmetric if 89 n_{j}\sigma_{ij}=-p_{i}$  \index{symmetric PDE}
90  \\  \begin{eqnarray}\label{LINEARPDE.SYSTEM.SYMMETRY TUTORIAL}
91  $p_{i}$ representing a uniform load on top of the beam.  A\hackscore{ijkl} =A\hackscore{klij} \\
92    \end{eqnarray}
93  A Dirichlet Boundary condition is assumed on the left end of the beam  Note that different from the scalar case now the coefficients $D$ and $d$ have to be inspected. It is easy to see that
94  for which no displacements can occure.\\  the Lame equation in fact is symmetric. The \LinearPDE class is notified by this fact by calling its \method{setSymmetryOn} method.
95  \\
96  \includegraphics[%  After we have solved the Lame equation we want to analyse the actual stress distribution. Typically the von--Mises stress\index{von--Mises stress} defined by
97    width=0.60\linewidth,bb = 0 0 200 100, draft, type=eps]{/home/jeannine/sandbox/report/draws/dir_cond_beam.eps}\\
98  This is described in the code with the setting a mask for the left  \sigma\hackscore{mises} = \sqrt{
99  end and setting values to that mask:  \frac{1}{6} ((\sigma\hackscore{00}-\sigma\hackscore{11})^2+
100                   (\sigma\hackscore{11}-\sigma\hackscore{22})^2
101                   (\sigma\hackscore{22}-\sigma\hackscore{00})^2)
102    +  \sigma\hackscore{01}^2+\sigma\hackscore{12}^2+\sigma\hackscore{20}^2}
103
104    is used to detect material damages. Here we want to calculate the von--Mises and write the stress to a file for visualization.
105
106    \index{scripts!\file{diffusion.py}}:
107    The following script which is available in \file{heatedbox.py} in the \ExampleDirectory solves the Lame equation
108    and writes the displacements and the von--Mises stress\index{von--Mises stress} into a file \file{deform.xml} in the \VTK file format:
109  \begin{python}  \begin{python}
110  q = xNodes{[}0{]}.whereZero(){*}{[}1.0,1.0{]}  from esys.escript import *
111    from esys.escript.linearPDEs import LinearPDE
112  r = Vector({[}0.0, 0.0{]}, where = nodes)  from esys.finley import Brick
113    #... set some parameters ...
114    lam=1.
115    mu=0.1
116    alpha=1.e-6
117    xc=[0.3,0.3,1.]
118    beta=8.
119    T_ref=0.
120    T_0=1.
121    #... generate domain ...
122    mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
123    x=mydomain.getX()
124    #... set temperature ...
125    T=T_0*exp(-beta*length(x-xc))
126    #... open symmetric PDE ...
127    mypde=LinearPDE(mydomain)
128    mypde.setSymmetryOn()
129    #... set coefficients ...
130    C=Tensor4(0.,Function(mydomain))
131    for i in range(mydomain.getDim()):
132      for j in range(mydomain.getDim()):
133         C[i,i,j,j]+=lam
134         C[j,i,j,i]+=mu
135         C[j,i,i,j]+=mu
136    msk=whereZero(x[0])*[1.,0.,0.] \
137       +whereZero(x[1])*[0.,1.,0.] \
138       +whereZero(x[2])*[0.,0.,1.]
139    sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
140    mypde.setValue(A=C,X=sigma0,q=msk)
141    #... solve pde ...
142    u=mypde.getSolution()
143    #... calculate von-Misses stress
145    sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
146    sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \
147                      (sigma[2,2]-sigma[0,0])**2)/6. \
148                     +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
149    #... output ...
150    saveVTK("deform.xml",disp=u,stress=sigma_mises)
151  \end{python}  \end{python}

$-(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}$
\\
with the natural boundary condition:

$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}$

Yields by comparing the coefficients :

\begin{enumerate}
\item $A_{ijkl}$= $\left[\mu\left(\delta_{ik}\delta_{ij}+\delta_{jl}\delta_{il}\right)+\lambda\left(\delta_{ij}\delta_{kl}\right)\right]$
\item $y_{i}$= $-p_{i}$
\item $u_{k}$= displacement $u$
\end{enumerate}
$B_{ijk,}=C_{ikl}=D_{ik}=X_{ij}=Y_{i}=d_{ik}=0$\\

Where 0 in the last line is taken as a scalar, vector or tensor, depending
on what the belonging coefficient is.
152
153  These equations are the base for the \textbf{Finley Script}:  \begin{figure}
154    \centerline{\includegraphics[width=\figwidth]{figures/HeatedBlock.eps}}
155    \caption{von--Mises Stress and Displacement Vectors.}
156    \label{HEATEDBLOCK FIG 2}
157    \end{figure}
158
159    Finally the the results can be visualize by calling
160  \begin{python}  \begin{python}
161  from ESyS import {*}  mayavi -d deform.xml -f CellToPointData -m VelocityVector -m SurfaceMap &
import Finley

\#mu lamda

def mu (E, nu): \#= shear modul G

return E/(2{*}(1+nu))

def lamda (E, nu):

return (nu{*}E)/((1-2{*}nu){*}(1+nu))

def main()

\#material, beam PROPERTIES

L = 10.0         \#length of beam {[}m{]}

h = 1            \#height of beam {[}m{]}

p = 1            \#outer uniform load {[}kN/m{]}

E0 = 210000      \#Young's modulus {[}kN/m\textasciicircum{}2{]}

nu = 0.4         \#Poisson ratio {[}-{]}

print \char\"{}L=\char\"{}, L

print \char\"{}h=\char\"{}, h

print \char\"{}p=\char\"{}, p

print \char\"{}E=\char\"{}, E0

print \char\"{}nu=Poisson =\char\"{}, nu

print \char\"{}mu=\char\"{}, mu (E0,nu)

print \char\"{}lamda=\char\"{}, lamda (E0,nu)

\#SET MESH for FE

mesh= Finley.Rectangle(n0=10 ,n1=1 ,order=2, l0=L, l1=h)

nodes = mesh.Nodes()

xNodes = nodes.getX()

elements = mesh.Elements()

faceElements = mesh.FaceElements()

xFaceElements = faceElements.getX()

\#DISPLACEMENT boundary

q = xNodes{[}0{]}.whereZero(){*}{[}1.,1.0{]}    \#setting the                                           mask for r

r = Vector({[}0.0, 0.0{]}, where = nodes) \#fixed end

\#STRESS boundary

y = Vector({[}0, -p{]}, where = faceElements)

\#Finley coeff.

A = Tensor4(0, where = elements)

for i in range (2) :

for j in range (2) :

A{[}i,i,j,j{]}+= lamda (E0,nu)

A{[}j,i,j,i{]}+= mu (E0,nu)

A{[}j,i,i,j{]}+= mu (E0,nu)

M,b = mesh.assemble(A=A, B=B, q=q,

y=y,r=r,type=CSR,num\_equations=2)

u= M.iterative(b, tolerance=1e-8,iter\_max=50000,

iterative\_method=PCG)

print \char\"{}u{[}0{]}=\char\"{},u{[}0{]}

print \char\"{}u{[}1{]}=\char\"{},u{[}1{]}

main()
162  \end{python}  \end{python}
163  The finer the mesh the more exact is the solution. E.g. a 20x2 mesh  Note that the filter \text{CellToPointData} is applied to create smooth representation of the
164  is more exact than a 10x1 mesh.  von--Mises stress. \fig{HEATEDBLOCK FIG 2} shows the results where the vertical planes showing the
165    von--Mises stress and the horizontal plane shows the vector representing displacements.
166

Legend:
 Removed from v.121 changed lines Added in v.757