2D_elasticity.py.rst 6.55 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3

..    # gedit: set fileencoding=utf8 :

Jeremy BLEYER's avatar
Jeremy BLEYER committed
4
.. _LinearElasticity2D:
Jeremy BLEYER's avatar
Jeremy BLEYER committed
5 6 7 8 9 10 11 12 13 14


=========================
 2D linear elasticity
=========================


Introduction
------------

15
In this first numerical tour, we will show how to compute a small strain solution for
Jeremy BLEYER's avatar
Jeremy BLEYER committed
16
a 2D isotropic linear elastic medium, either in plane stress or in plane strain,
17 18
in a tradtional displacement-based finite element formulation. The corresponding
file can be obtained from :download:`2D_elasticity.py`. Extension to 3D
Jeremy BLEYER's avatar
Jeremy BLEYER committed
19 20 21 22 23 24 25 26 27
is straightforward and an example can be found in the :ref:`ModalAnalysis` example.

We consider here the case of a cantilever beam modeled as a 2D medium of dimensions
:math:`L\times  H`. Geometrical parameters and mesh density are first defined
and the rectangular domain is  generated using the ``RectangleMesh`` function.
We also choose a criss-crossed structured mesh::

 from __future__ import print_function
 from fenics import *
28

Jeremy BLEYER's avatar
Jeremy BLEYER committed
29 30 31 32 33 34 35 36 37 38 39
 L = 25.
 H = 1.
 Nx = 250
 Ny = 10
 mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")


Constitutive relation
---------------------

We now define the material parameters which are here given in terms of a Young's
40 41 42 43
modulus :math:`E` and a Poisson coefficient :math:`\nu`. In the following, we will
need to define the constitutive relation between the stress tensor :math:`\boldsymbol{\sigma}`
and the strain tensor :math:`\boldsymbol{\varepsilon}`. Let us recall
that the general expression of the linear elastic isotropic constitutive relation
Jeremy BLEYER's avatar
Jeremy BLEYER committed
44 45 46 47 48 49 50 51 52 53 54
for a 3D medium is given by:

.. math::
 \boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon}
 :label: constitutive_3D

for a natural (no prestress) initial state where the Lamé coefficients are given by:

.. math::
 \lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \dfrac{E}{2(1+\nu)}
 :label: Lame_coeff
55 56

In this demo, we consider a 2D model either in plane strain or in plane stress conditions.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
57 58
Irrespective of this choice, we will work only with a 2D displacement vector :math:`\boldsymbol{u}=(u_x,u_y)`
and will subsequently define the strain operator ``eps`` as follows::
59

Jeremy BLEYER's avatar
Jeremy BLEYER committed
60 61 62
 def eps(v):
     return sym(grad(v))

63
which computes the 2x2 plane components of the symmetrized gradient tensor of
Jeremy BLEYER's avatar
Jeremy BLEYER committed
64 65 66 67 68
any 2D vectorial field. In the plane strain case, the full 3D strain tensor is defined as follows:

.. math::
 \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & 0\\
 \varepsilon_{xy} & \varepsilon_{yy} & 0 \\ 0 & 0 & 0\end{bmatrix}
69 70

so that the 2x2 plane part of the stress tensor is defined in the same way as for the 3D case
Jeremy BLEYER's avatar
Jeremy BLEYER committed
71 72
(the out-of-plane stress component being given by :math:`\sigma_{zz}=\lambda(\varepsilon_{xx}+\varepsilon_{yy})`.

73
In the plane stress case, an out-of-plane strain component :math:`\varepsilon_{zz}`
Jeremy BLEYER's avatar
Jeremy BLEYER committed
74
must be considered so that :math:`\sigma_{zz}=0`. Using this condition in the
75
3D constitutive relation, one has :math:`\varepsilon_{zz}=-\dfrac{\lambda}{\lambda+2\mu}(\varepsilon_{xx}+\varepsilon_{yy})`.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
76 77 78 79 80 81
Injecting into :eq:`constitutive_3D`, we have for the 2D plane stress relation:

.. math::
 \boldsymbol{\sigma} = \lambda^* \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon}

where :math:`\boldsymbol{\sigma}, \boldsymbol{\varepsilon}, \mathbf{1}` are 2D tensors and with
82 83
:math:`\lambda^* = \dfrac{2\lambda\mu}{\lambda+2\mu}`. Hence, the 2D constitutive relation
is identical to the plane strain case by changing only the value of the Lamé coefficient :math:`\lambda`.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
84
We can then have::
Jeremy BLEYER's avatar
Jeremy BLEYER committed
85 86 87 88

 E = Constant(1e5)
 nu = Constant(0.3)
 model = "plane_stress"
89

Jeremy BLEYER's avatar
Jeremy BLEYER committed
90 91 92
 mu = E/2/(1+nu)
 lmbda = E*nu/(1+nu)/(1-2*nu)
 if model == "plane_stress":
Jeremy BLEYER's avatar
Jeremy BLEYER committed
93
     lmbda = 2*mu*lmbda/(lmbda+2*mu)
94

Jeremy BLEYER's avatar
Jeremy BLEYER committed
95 96 97 98
 def sigma(v):
     return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

.. note::
99
 Note that we used the variable name ``lmbda`` to avoid any confusion with the
Jeremy BLEYER's avatar
Jeremy BLEYER committed
100
 lambda functions of Python
101

Jeremy BLEYER's avatar
Jeremy BLEYER committed
102
 We also used an intrinsic formulation of the constitutive relation. Example of
103
 constitutive relation implemented with a matrix/vector engineering notation
Jeremy BLEYER's avatar
Jeremy BLEYER committed
104
 will be provided in the :ref:`OrthotropicElasticity` example.
105 106


Jeremy BLEYER's avatar
Jeremy BLEYER committed
107 108
Variational formulation
-----------------------
109

Jeremy BLEYER's avatar
Jeremy BLEYER committed
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
For this example, we consider a continuous polynomial interpolation of degree 2
and a uniformly distributed loading :math:`\boldsymbol{f}=(0,-f)` corresponding
to the beam self-weight. The continuum mechanics variational formulation (obtained
from the virtual work principle) is given by:

.. math::
 \text{Find } \boldsymbol{u}\in V \text{ s.t. } \int_{\Omega}
 \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v}) d\Omega
 = \int_{\Omega} \boldsymbol{f}\cdot\boldsymbol{v}  d\Omega \quad \forall\boldsymbol{v} \in V

which translates into the following FEniCS code::

 rho_g = 1e-3
 f = Constant((0,-rho_g))

 V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
 du = TrialFunction(V)
 u_ = TestFunction(V)
 a = inner(sigma(du), eps(u_))*dx
 l = inner(f, u_)*dx


Resolution
----------

Fixed displacements are imposed on the left part of the beam, the ``solve``
function is then called and solution is plotted by deforming the mesh::

 def left(x, on_boundary):
     return near(x[0],0.)

 bc = DirichletBC(V, Constant((0.,0.)), left)

 u = Function(V, name="Displacement")
 solve(a == l, u, bc)

 plot(1e3*u, mode="displacement")

The (amplified) solution should look like this:

.. image:: cantilever_deformed.png
   :scale: 15%


Validation and post-processing
------------------------------

The maximal deflection is compared against the analytical solution from
Euler-Bernoulli beam theory which is here :math:`w_{beam} = \dfrac{qL^4}{8EI}`::

 print("Maximal deflection:", -u(L,H/2.)[1])
 print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))

163
One finds :math:`w_{FE} = 5.8638\text{e-3}` against :math:`w_{beam} = 5.8594\text{e-3}`
Jeremy BLEYER's avatar
Jeremy BLEYER committed
164 165 166 167 168 169 170 171 172 173 174 175
that is a 0.07% difference.


The stress tensor must be projected on an appropriate function space in order to
evaluate pointwise values or export it for Paraview vizualisation. Here we choose
to describe it as a (2D) tensor and project it onto a piecewise constant function
space::

 Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
 sig = Function(Vsig, name="Stress")
 sig.assign(project(sigma(u), Vsig))
 print("Stress at (0,H):", sig(0, H))
176

Jeremy BLEYER's avatar
Jeremy BLEYER committed
177 178 179
Fields can be exported in a suitable format for vizualisation using Paraview.
VTK-based extensions (.pvd,.vtu) are not suited for multiple fields and parallel
writing/reading. Prefered output format is now .xdmf::
180

Jeremy BLEYER's avatar
Jeremy BLEYER committed
181 182 183 184 185
 file_results = XDMFFile("elasticity_results.xdmf")
 file_results.parameters["flush_output"] = True
 file_results.parameters["functions_share_mesh"] = True
 file_results.write(u, 0.)
 file_results.write(sig, 0.)