#
# .. # gedit: set fileencoding=utf8 :
# .. raw:: html
#
#

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

#
# .. _OrthotropicElasticity:
#
# ===============================
# Orthotropic linear elasticity
# ===============================
#
#
# Introduction
# ------------
#
# 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
# :download:`orthotropic_elasticity.py`.
#
# We consider here the case of a square plate perforated by a circular hole of
# 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
domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R)
mesh = generate_mesh(domain, N)
# Constitutive relation
# ---------------------
#
# Constitutive relations will be defined using an engineering (or Voigt) notation (i.e.
# 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\\
# -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix}
# \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
# relation symmetry :math:`\nu_{yx}=\nu_{xy}E_y/E_x`) and :math:`G_{xy}` being the
# 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
# --------------------------------
#
# Different parts of the quarter plate boundaries are now defined as well as the
# 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)
# Symmetric boundary conditions are applied on the ``Top`` and ``Left`` boundaries
# 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
#