orthotropic_elasticity.py 5.78 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2
#
# ..    # gedit: set fileencoding=utf8 :
3 4 5
# .. raw:: html
#
#  <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><p align="center"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png"/></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>
Jeremy BLEYER's avatar
Jeremy BLEYER committed
6 7 8 9 10 11 12 13 14 15 16
#
# .. _OrthotropicElasticity:
#
# ===============================
#  Orthotropic linear elasticity
# ===============================
#
#
# Introduction
# ------------
#
17
# In this numerical tour, we will show how to tackle the case of orthotropic elasticity (in a 2D setting). The corresponding file can be obtained from
Jeremy BLEYER's avatar
Jeremy BLEYER committed
18 19
# :download:`orthotropic_elasticity.py`.
#
20
# We consider here the case of a square plate perforated by a circular hole of
Jeremy BLEYER's avatar
Jeremy BLEYER committed
21 22 23 24 25 26 27 28 29 30 31 32
# radius :math:`R`, the plate dimension is :math:`2L\times 2L` with :math:`L \gg R`
# Only the top-right quarter of the plate will be considered. Loading will consist
# of a uniform traction on the top/bottom boundaries, symmetry conditions will also
# be applied on the correponding symmetry planes. To generate the perforated domain
# we use here the ``mshr`` module and define the boolean "*minus*" operation
# between a rectangle and a circle::

from fenics import *
from mshr import *

L, R = 1., 0.1
N = 50 # mesh density
33

Jeremy BLEYER's avatar
Jeremy BLEYER committed
34 35 36 37 38 39 40
domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R)
mesh = generate_mesh(domain, N)


# Constitutive relation
# ---------------------
#
41
# Constitutive relations will be defined using an engineering (or Voigt) notation (i.e.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
42 43 44 45 46 47 48 49 50 51
# second order tensors will be written as a vector of their components) contrary
# to the :ref:`LinearElasticity2D` example which used an intrinsic notation. In
# the material frame, which is assumed to coincide here with the global :math:`(Oxy)`
# frame, the orthotropic constitutive law writes :math:`\boldsymbol{\varepsilon}=\mathbf{S}
# \boldsymbol{\sigma}` using the compliance matrix
# :math:`\mathbf{S}` with:
#
# .. math::
#   \begin{Bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ 2\varepsilon_{xy}
#   \end{Bmatrix} = \begin{bmatrix} 1/E_x & -\nu_{xy}/E_x & 0\\
52
#   -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix}
Jeremy BLEYER's avatar
Jeremy BLEYER committed
53 54 55 56 57
#   \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy}
#   \end{Bmatrix}
#
# with :math:`E_x, E_y` the two Young's moduli in the orthotropy directions, :math:`\nu_{xy}`
# the in-plane Poisson ration (with the following relation ensuring the constitutive
58
# relation symmetry :math:`\nu_{yx}=\nu_{xy}E_y/E_x`) and :math:`G_{xy}` being the
Jeremy BLEYER's avatar
Jeremy BLEYER committed
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
# shear modulus. This relation needs to be inverted to obtain the stress components as a function
# of the strain components :math:`\boldsymbol{\sigma}=\mathbf{C}\boldsymbol{\varepsilon}` with
# :math:`\mathbf{C}=\mathbf{S}^{-1}`::

Ex, Ey, nuxy, Gxy = 100., 10., 0.3, 5.
S = as_matrix([[1./Ex,nuxy/Ex,0.],[nuxy/Ex,1./Ey,0.],[0.,0.,1./Gxy]])
C = inv(S)

# .. note::
#  Here we used the ``inv`` opertor to compute the elasticity matrix :math:`\mathbf{C}`.
#  We could also have computed analytically the inverse relation. Note that the ``inv``
#  operator is implemented only up to 3x3 matrices. Extension to the 3D case yields 6x6
#  matrices and therefore requires either analytical inversion or numerical inversion
#  using Numpy for instance (assuming that the material parameters are constants).
#
# We define different functions for representing the stress and strain either as
# second-order tensor or using the Voigt engineering notation::

def eps(v):
    return sym(grad(v))
def strain2voigt(e):
    """e is a 2nd-order tensor, returns its Voigt vectorial representation"""
    return as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
    """
    s is a stress-like vector (no 2 factor on last component)
    returns its tensorial representation
    """
    return as_tensor([[s[0],s[2]],[s[2],s[1]]])
def sigma(v):
    return voigt2stress(dot(C,strain2voigt(eps(v))))


# Problem position and resolution
# --------------------------------
#
95
# Different parts of the quarter plate boundaries are now defined as well as the
Jeremy BLEYER's avatar
Jeremy BLEYER committed
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
# exterior integration measure ``ds``::

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],L) and on_boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0) and on_boundary
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1],0) and on_boundary

# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
Left().mark(facets, 2)
Bottom().mark(facets, 3)
ds = Measure('ds')[facets]

# We are now in position to define the variational form which is given as in :ref:`LinearElasticity2D`,
# the linear form now contains a Neumann term corresponding to a uniform vertical traction :math:`\sigma_{\infty}`
# on the top boundary::

# Define function space
V = VectorFunctionSpace(mesh, 'Lagrange', 2)

# Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name='Displacement')
a = inner(sigma(du), eps(u_))*dx

# uniform traction on top boundary
T = Constant((0, 1e-3))
l = dot(T, u_)*ds(1)

133
# Symmetric boundary conditions are applied on the ``Top`` and ``Left`` boundaries
Jeremy BLEYER's avatar
Jeremy BLEYER committed
134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
# and the problem is solved::

# symmetry boundary conditions
bc = [DirichletBC(V.sub(0), Constant(0.), facets, 2),
      DirichletBC(V.sub(1), Constant(0.), facets, 3)]

solve(a == l, u, bc)

import matplotlib.pyplot as plt
p = plot(sigma(u)[1,1]/T[1], mode='color')
plt.colorbar(p)
plt.title(r"$\sigma_{yy}$",fontsize=26)

# The :math:`\sigma_{xx}` and :math:`\sigma_{yy}` components should look like
# that:
#
# .. image:: circular_hole_sigxx.png
#    :scale: 11 %
#    :align: left
# .. image:: circular_hole_sigyy.png
#    :scale: 11 %
#    :align: right
#