doc/user/beam.tex

% $Id$
\chapter{Bending Beam under Uniform Load}
\label{BEAM CHAP}
\sectionauthor{Jannine Eisenmann}{}

Given is a two-dimension bending beam which is fixed in all directions
at the left end and free at the other, see \fig{BEAM FIG 1}. Furthermore the beam is loaded
with a uniform load $p$.

\begin{figure}
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
\caption{Loaded Beam.}
\label{BEAM FIG 1}
\end{figure}


For emphasizing the displacement this picture is drawn (uebertrieben),
cause since we use the linear theory otherwise no displacements would
be visible. 
There are two ways of solving this problem: on the one hand one can
use the differential equation for the deformed neutral fibre of the
beam. This classical differential equation is based on several simplifying
assumptions concerning the statics and kinematics of the problem.
However the results are known to be highly accurate at least for slender
beams with length to hight ratios $> 10$. Alternatively, in connection
with finite element based differential equation toolkits one may simply
consider the beam as a special case of a 2D or 3D elastic continuum
problem and solve the stress equilibrium equations combined with Hooke's
law for the specific boundary conditions of a beam. Both cases assume
isotropic and linear elastic material properties.

The beam equation is easily solved analytically. The analytical solutions
can be used for benchmarkung finite element solutions. In section
1.2 we formluate a finite element code for general isotropic elasticity
problems including thin and deep beams under arbitrary loading conditiond
in 2D or 3D. 


\section{2-dimensional}
As the stress equilibrium equation is a partial differential equation
(PDE), we choose to use Finley to solve it, since Finley is a finite
element kernel library, that has been incorporated as a differential
equation solver into the python based numerical toolkit called escript.
We divided the beam into ten square elements of the size 1x1. Each
element consists of 8 nodes, which leads to a quadratic interpolation
of the model point displacements \\

The key ingredients is \textbf{Hooks Law}. We use Hooks Law because
we are dealing with \textbf{linear elastic material} \textbf{behaviour}.
We have \\


$\sigma_{ik}=2G\left(\varepsilon_{ik}+\frac{\nu}{1-2\nu}\cdot e\cdot\delta_{ik}\right)$\hfill{}(1)\\
where the engineering strain$\varepsilon_{ik}$is defined as:

$\varepsilon_{ik}=\frac{1}{2}\cdot\left(u_{k,i}+u_{i,k}\right)$\hfill{}(2)\\


with

\begin{enumerate}
\item e= Volume strain = $\varepsilon_{kk}$
\item $\delta_{ik}$= Kronecker symbol
\end{enumerate}
Inserting equation (2) in (1) (and with further mathematical conversions)
leads to the following partial differential equation :\\


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


For tension equilibrium we require:\\


$\sigma_{ij,j}=0$\\


which leads to the following expression:\\


\[
\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\]


with the natural boundary condition:

\[
n_{j}\sigma_{ij}=-p_{i}\]
\\
$p_{i}$ representing a uniform load on top of the beam.

A Dirichlet Boundary condition is assumed on the left end of the beam
for which no displacements can occure.\\
\\
\includegraphics[%
  width=0.60\linewidth,bb = 0 0 200 100, draft, type=eps]{/home/jeannine/sandbox/report/draws/dir_cond_beam.eps}\\
This is described in the code with the setting a mask for the left
end and setting values to that mask:

\begin{python}
q = xNodes{[}0{]}.whereZero(){*}{[}1.0,1.0{]}

r = Vector({[}0.0, 0.0{]}, where = nodes)
\end{python}
The Finley template PDE reads:\\


\[
-(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.

These equations are the base for the \textbf{Finley Script}:

\begin{python}
from ESyS import {*}
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

     ymask = xFaceElements{[}1{]}.whereEqualTo(h)

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

     y = y{*}ymask


     \#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()
\end{python}
The finer the mesh the more exact is the solution. E.g. a 20x2 mesh
is more exact than a 10x1 mesh. 

1	jgs	110	% $Id$
2			\chapter{Bending Beam under Uniform Load}
3			\label{BEAM CHAP}
4			\sectionauthor{Jannine Eisenmann}{}
5
6			Given is a two-dimension bending beam which is fixed in all directions
7			at the left end and free at the other, see \fig{BEAM FIG 1}. Furthermore the beam is loaded
8			with a uniform load $p$.
9
10			\begin{figure}
11			% \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
12			\caption{Loaded Beam.}
13			\label{BEAM FIG 1}
14			\end{figure}
15
16
17
18			For emphasizing the displacement this picture is drawn (uebertrieben),
19			cause since we use the linear theory otherwise no displacements would
20			be visible.
21			There are two ways of solving this problem: on the one hand one can
22			use the differential equation for the deformed neutral fibre of the
23			beam. This classical differential equation is based on several simplifying
24			assumptions concerning the statics and kinematics of the problem.
25			However the results are known to be highly accurate at least for slender
26			beams with length to hight ratios $> 10$. Alternatively, in connection
27			with finite element based differential equation toolkits one may simply
28			consider the beam as a special case of a 2D or 3D elastic continuum
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
31			isotropic and linear elastic material properties.
32
33			The beam equation is easily solved analytically. The analytical solutions
34			can be used for benchmarkung finite element solutions. In section
35			1.2 we formluate a finite element code for general isotropic elasticity
36			problems including thin and deep beams under arbitrary loading conditiond
37			in 2D or 3D.
38
39
40			\section{2-dimensional}
41			As the stress equilibrium equation is a partial differential equation
42			(PDE), we choose to use Finley to solve it, since Finley is a finite
43			element kernel library, that has been incorporated as a differential
44			equation solver into the python based numerical toolkit called escript.
45			We divided the beam into ten square elements of the size 1x1. Each
46			element consists of 8 nodes, which leads to a quadratic interpolation
47			of the model point displacements \\
48
49			The key ingredients is \textbf{Hooks Law}. We use Hooks Law because
50			we are dealing with \textbf{linear elastic material} \textbf{behaviour}.
51			We have \\
52
53
54			$\sigma_{ik}=2G\left(\varepsilon_{ik}+\frac{\nu}{1-2\nu}\cdot e\cdot\delta_{ik}\right)$\hfill{}(1)\\
55			where the engineering strain$\varepsilon_{ik}$is defined as:
56
57			$\varepsilon_{ik}=\frac{1}{2}\cdot\left(u_{k,i}+u_{i,k}\right)$\hfill{}(2)\\
58
59
60			with
61
62			\begin{enumerate}
63			\item e= Volume strain = $\varepsilon_{kk}$
64			\item $\delta_{ik}$= Kronecker symbol
65			\end{enumerate}
66			Inserting equation (2) in (1) (and with further mathematical conversions)
67			leads to the following partial differential equation :\\
68
69
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}$\\
71
72
73			For tension equilibrium we require:\\
74
75
76			$\sigma_{ij,j}=0$\\
77
78
79			which leads to the following expression:\\
80
81
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\]
84
85
86			with the natural boundary condition:
87
88			\[
89			n_{j}\sigma_{ij}=-p_{i}\]
90			\\
91			$p_{i}$ representing a uniform load on top of the beam.
92
93			A Dirichlet Boundary condition is assumed on the left end of the beam
94			for which no displacements can occure.\\
95			\\
96			\includegraphics[%
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
99			end and setting values to that mask:
100
101			\begin{python}
102			q = xNodes{[}0{]}.whereZero(){*}{[}1.0,1.0{]}
103
104			r = Vector({[}0.0, 0.0{]}, where = nodes)
105			\end{python}
106			The Finley template PDE reads:\\
107
108
109			\[
110			-(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}\]
111			\\
112			with the natural boundary condition:
113
114			\[
115			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}\]
116
117
118			Yields by comparing the coefficients :
119
120			\begin{enumerate}
121			\item $A_{ijkl}$= $\left[\mu\left(\delta_{ik}\delta_{ij}+\delta_{jl}\delta_{il}\right)+\lambda\left(\delta_{ij}\delta_{kl}\right)\right]$
122			\item $y_{i}$= $-p_{i}$
123			\item $u_{k}$= displacement $u$
124			\end{enumerate}
125			$B_{ijk,}=C_{ikl}=D_{ik}=X_{ij}=Y_{i}=d_{ik}=0$\\
126
127
128			Where 0 in the last line is taken as a scalar, vector or tensor, depending
129			on what the belonging coefficient is.
130
131			These equations are the base for the \textbf{Finley Script}:
132
133			\begin{python}
134			from ESyS import {*}
135			import Finley
136
137
138
139			\#mu lamda
140
141			def mu (E, nu): \#= shear modul G
142
143			return E/(2{*}(1+nu))
144
145			def lamda (E, nu):
146
147			return (nu{}E)/((1-2{}nu){*}(1+nu))
148
149
150
151			def main()
152
153			\#material, beam PROPERTIES
154
155			L = 10.0 \#length of beam {[}m{]}
156
157			h = 1 \#height of beam {[}m{]}
158
159			p = 1 \#outer uniform load {[}kN/m{]}
160
161			E0 = 210000 \#Young's modulus {[}kN/m\textasciicircum{}2{]}
162
163			nu = 0.4 \#Poisson ratio {[}-{]}
164
165
166
167			print \char`\"{}L=\char`\"{}, L
168
169			print \char`\"{}h=\char`\"{}, h
170
171			print \char`\"{}p=\char`\"{}, p
172
173			print \char`\"{}E=\char`\"{}, E0
174
175			print \char`\"{}nu=Poisson =\char`\"{}, nu
176
177			print \char`\"{}mu=\char`\"{}, mu (E0,nu)
178
179			print \char`\"{}lamda=\char`\"{}, lamda (E0,nu)
180
181
182
183			\#SET MESH for FE
184
185			mesh= Finley.Rectangle(n0=10 ,n1=1 ,order=2, l0=L, l1=h)
186
187			nodes = mesh.Nodes()
188
189			xNodes = nodes.getX()
190
191			elements = mesh.Elements()
192
193			faceElements = mesh.FaceElements()
194
195			xFaceElements = faceElements.getX()
196
197
198
199			\#DISPLACEMENT boundary
200
201			q = xNodes{[}0{]}.whereZero(){*}{[}1.,1.0{]} \#setting the mask for r
202
203			r = Vector({[}0.0, 0.0{]}, where = nodes) \#fixed end
204
205
206
207			\#STRESS boundary
208
209			ymask = xFaceElements{[}1{]}.whereEqualTo(h)
210
211			y = Vector({[}0, -p{]}, where = faceElements)
212
213			y = y{*}ymask
214
215
216
217			\#Finley coeff.
218
219			A = Tensor4(0, where = elements)
220
221			for i in range (2) :
222
223			for j in range (2) :
224
225			A{[}i,i,j,j{]}+= lamda (E0,nu)
226
227			A{[}j,i,j,i{]}+= mu (E0,nu)
228
229			A{[}j,i,i,j{]}+= mu (E0,nu)
230
231
232
233			M,b = mesh.assemble(A=A, B=B, q=q,
234
235			y=y,r=r,type=CSR,num\_equations=2)
236
237			u= M.iterative(b, tolerance=1e-8,iter\_max=50000,
238
239			iterative\_method=PCG)
240
241			print \char`\"{}u{[}0{]}=\char`\"{},u{[}0{]}
242
243			print \char`\"{}u{[}1{]}=\char`\"{},u{[}1{]}
244
245			main()
246			\end{python}
247			The finer the mesh the more exact is the solution. E.g. a 20x2 mesh
248			is more exact than a 10x1 mesh.
249
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision