doc/user/beam.tex

% $Id$
\section{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	% $Id$
2	\section{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