2D_elasticity.py.rst 6.55 KB
 Jeremy BLEYER committed Apr 17, 2018 1 2 3  .. # gedit: set fileencoding=utf8 :  Jeremy BLEYER committed Apr 23, 2018 4 .. _LinearElasticity2D:  Jeremy BLEYER committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 17, 2018 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 committed Apr 19, 2018 84 We can then have::  Jeremy BLEYER committed Apr 17, 2018 85 86 87 88  E = Constant(1e5) nu = Constant(0.3) model = "plane_stress"  89   Jeremy BLEYER committed Apr 17, 2018 90 91 92  mu = E/2/(1+nu) lmbda = E*nu/(1+nu)/(1-2*nu) if model == "plane_stress":  Jeremy BLEYER committed Apr 19, 2018 93  lmbda = 2*mu*lmbda/(lmbda+2*mu)  94   Jeremy BLEYER committed Apr 17, 2018 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 committed Apr 17, 2018 100  lambda functions of Python  101   Jeremy BLEYER committed Apr 17, 2018 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 committed Apr 17, 2018 104  will be provided in the :ref:OrthotropicElasticity example.  105 106   Jeremy BLEYER committed Apr 17, 2018 107 108 Variational formulation -----------------------  109   Jeremy BLEYER committed Apr 17, 2018 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 committed Apr 19, 2018 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 committed Apr 19, 2018 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 committed Apr 19, 2018 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.)