orthotropic_elasticity.py 5.78 KB
 Jeremy BLEYER committed May 28, 2018 1 2 # # .. # gedit: set fileencoding=utf8 :  Jeremy BLEYER committed Oct 04, 2018 3 4 5 # .. raw:: html # #
 Jeremy BLEYER committed May 28, 2018 6 7 8 9 10 11 12 13 14 15 16 # # .. _OrthotropicElasticity: # # =============================== # Orthotropic linear elasticity # =============================== # # # Introduction # ------------ #  Jeremy BLEYER committed Oct 04, 2018 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 committed May 28, 2018 18 19 # :download:orthotropic_elasticity.py. #  Jeremy BLEYER committed Oct 04, 2018 20 # We consider here the case of a square plate perforated by a circular hole of  Jeremy BLEYER committed May 28, 2018 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  Jeremy BLEYER committed Oct 04, 2018 33   Jeremy BLEYER committed May 28, 2018 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 # --------------------- #  Jeremy BLEYER committed Oct 04, 2018 41 # Constitutive relations will be defined using an engineering (or Voigt) notation (i.e.  Jeremy BLEYER committed May 28, 2018 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\\  Jeremy BLEYER committed Oct 04, 2018 52 # -\nu_{yx}/E_y & 1/E_y & 0 \\ 0 & 0 & 1/G_{xy} \end{bmatrix}\begin{Bmatrix}  Jeremy BLEYER committed May 28, 2018 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  Jeremy BLEYER committed Oct 04, 2018 58 # relation symmetry :math:\nu_{yx}=\nu_{xy}E_y/E_x) and :math:G_{xy} being the  Jeremy BLEYER committed May 28, 2018 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 # -------------------------------- #  Jeremy BLEYER committed Oct 04, 2018 95 # Different parts of the quarter plate boundaries are now defined as well as the  Jeremy BLEYER committed May 28, 2018 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)  Jeremy BLEYER committed Oct 04, 2018 133 # Symmetric boundary conditions are applied on the Top and Left boundaries  Jeremy BLEYER committed May 28, 2018 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 #