vonMises_plasticity.py.rst 17.3 KB
 Jeremy BLEYER committed Apr 25, 2018 1   Jeremy BLEYER committed Oct 04, 2018 2 3 4 5 .. raw:: html
 Jeremy BLEYER committed Jan 22, 2019 6 7 .. _vonMisesPlasticity:  Jeremy BLEYER committed Apr 25, 2018 8 ==================================================  9 Elasto-plastic analysis of a 2D von Mises material  Jeremy BLEYER committed Apr 25, 2018 10 11 12 13 14 15 ================================================== ------------- Introduction -------------  16 This example is concerned with the incremental analysis of an elasto-plastic  Jeremy BLEYER committed Apr 25, 2018 17 18 19 20 von Mises material. The structure response is computed using an iterative predictor-corrector return mapping algorithm embedded in a Newton-Raphson global loop for restoring equilibrium. Due to the simple expression of the von Mises criterion, the return mapping procedure is completely analytical (with linear isotropic  21 hardening). The corresponding sources can be obtained from :download:vonMises_plasticity.zip.  Jeremy BLEYER committed Apr 25, 2018 22 Another implementation of von Mises plasticity can also be found at  Jeremy BLEYER committed Jan 22, 2019 23 24 https://bitbucket.org/fenics-apps/fenics-solid-mechanics as well as in the numerical tour :ref:PlasticityMFront.  Jeremy BLEYER committed Apr 25, 2018 25 26 27  We point out that the 2D nature of the problem will impose keeping track of the out-of-plane :math:\varepsilon_{zz}^p plastic strain and dealing with  28 representations of stress/strain states including the :math:zz component. Note  Jeremy BLEYER committed Apr 25, 2018 29 also that are not concerned with the issue of volumetric locking induced by  30 31 incompressible plastic deformations since quadratic triangles in 2D is enough to mitigate the locking phenomenon.  Jeremy BLEYER committed Apr 25, 2018 32 33 34 35 36 37 38 39 40 41 42  The plastic strain evolution during the cylinder expansion will look like this: .. image:: plastic_strain.gif :scale: 80% ----------------- Problem position -----------------  43 In FEniCS 2017.2, the FEniCS Form Compiler ffc now uses uflacs as a default  Jeremy BLEYER committed Apr 25, 2018 44 45 46 47 representation instead of the old quadrature representation. However, using Quadrature elements generates some bugs in this representation and we therefore need to revert to the old representation. Deprecation warning messages are also disabled. See this post _  48 and the corresponding issue _  Jeremy BLEYER committed Apr 25, 2018 49 50 51 for more information.:: from __future__ import print_function  Jeremy BLEYER committed Dec 16, 2018 52  from dolfin import *  Jeremy BLEYER committed Apr 25, 2018 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71  import numpy as np parameters["form_compiler"]["representation"] = 'quadrature' import warnings from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning) The material is represented by an isotropic elasto-plastic von Mises yield condition of uniaxial strength :math:\sigma_0 and with isotropic hardening of modulus :math:H. The yield condition is thus given by: .. math:: f(\boldsymbol{\sigma}) = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}} - \sigma_0 -Hp \leq 0 where :math:p is the cumulated equivalent plastic strain. The hardening modulus can also be related to a tangent elastic modulus :math:E_t = \dfrac{EH}{E+H}. The considered problem is that of a plane strain hollow cylinder of internal (resp. external)  72 radius :math:R_i (resp. :math:R_e) under internal uniform pressure :math:q.  Jeremy BLEYER committed Apr 25, 2018 73 Only a quarter of cylinder is generated using Gmsh and converted to .xml format. ::  74   Jeremy BLEYER committed Apr 25, 2018 75 76 77 78 79 80 81 82 83 84 85 86 87  # elastic parameters E = Constant(70e3) nu = Constant(0.3) lmbda = E*nu/(1+nu)/(1-2*nu) mu = E/2./(1+nu) sig0 = Constant(250.) # yield strength Et = E/100. # tangent modulus H = E*Et/(E-Et) # hardening modulus Re, Ri = 1.3, 1. # external/internal radius mesh = Mesh("thick_cylinder.xml") facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml") ds = Measure('ds')[facets]  88   Jeremy BLEYER committed Apr 25, 2018 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 Function spaces will involve a standard CG space for the displacement whereas internal state variables such as plastic strains will be represented using a **Quadrature** element. This choice will make it possible to express the complex non-linear material constitutive equation at the Gauss points only, without involving any interpolation of non-linear expressions throughout the element. It will ensure an optimal convergence rate for the Newton-Raphson method. See Chapter 26 of the FEniCS book _. We will need Quadrature elements for 4-dimensional vectors and scalars, the number of Gauss points will be determined by the required degree of the Quadrature element (e.g. degree=1 yields only 1 Gauss point, degree=2 yields 3 Gauss points and degree=3 yields 6 Gauss points (note that this is suboptimal)):: deg_u = 2 deg_stress = 2 V = VectorFunctionSpace(mesh, "CG", deg_u) We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default') W = FunctionSpace(mesh, We) W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default') W0 = FunctionSpace(mesh, W0e)  107   Jeremy BLEYER committed Apr 25, 2018 108 109 110 .. note:: In older versions, it was possible to define Quadrature function spaces directly using FunctionSpace(mesh, "Quadrature", 1). This is no longer the case since  111  FEniCS 2016.1 (see this issue _). Instead, Quadrature elements must first be defined  Jeremy BLEYER committed Apr 25, 2018 112 113  by specifying the associated degree and quadrature scheme before defining the associated FunctionSpace.  114 115 116   Jeremy BLEYER committed Apr 25, 2018 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 Various functions are defined to keep track of the current internal state and currently computed increments:: sig = Function(W) sig_old = Function(W) n_elas = Function(W) beta = Function(W0) p = Function(W0, name="Cumulative plastic strain") u = Function(V, name="Total displacement") du = Function(V, name="Iteration correction") Du = Function(V, name="Current increment") v = TrialFunction(V) u_ = TestFunction(V) Boundary conditions correspond to symmetry conditions on the bottom and left  132 parts (resp. numbered 1 and 3). Loading consists of a uniform pressure on the  Jeremy BLEYER committed Apr 25, 2018 133 134 135 136 137 138 internal boundary (numbered 4). It will be progressively increased from 0 to :math:q_{lim}=\dfrac{2}{\sqrt{3}}\sigma_0\log\left(\dfrac{R_e}{R_i}\right) which is the analytical collapse load for a perfectly-plastic material (no hardening):: bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]  139   Jeremy BLEYER committed Apr 25, 2018 140 141 142 143 144 145  n = FacetNormal(mesh) q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0) loading = Expression("-q*t", q=q_lim, t=0, degree=2) def F_ext(v): return loading*dot(n, v)*ds(4)  146   Jeremy BLEYER committed Apr 25, 2018 147 148 149 150 151 152 153 154 ----------------------------- Constitutive relation update ----------------------------- Before writing the variational form, we now define some useful functions which will enable performing the constitutive relation update using a return mapping procedure. This step is quite classical in FEM plasticity for a von Mises criterion with isotropic hardening and follow notations from [BON2014]_. First, the strain  155 tensor will be represented in a 3D fashion by appending zeros on the out-of-plane  Jeremy BLEYER committed Apr 25, 2018 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 components since, even if the problem is 2D, the plastic constitutive relation will involve out-of-plane plastic strains. The elastic consitutive relation is also defined and a function as_3D_tensor will enable to represent a 4 dimensional vector containing respectively :math:xx, yy, zz and :math:xy components as a 3D tensor:: def eps(v): e = sym(grad(v)) return as_tensor([[e[0, 0], e[0, 1], 0], [e[0, 1], e[1, 1], 0], [0, 0, 0]]) def sigma(eps_el): return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el def as_3D_tensor(X): return as_tensor([[X[0], X[3], 0], [X[3], X[1], 0], [0, 0, X[2]]]) The return mapping procedure consists in finding a new stress :math:\boldsymbol{\sigma}_{n+1}  174 175 and internal variable :math:p_{n+1} state verifying the current plasticity condition from a previous stress :math:\boldsymbol{\sigma}_{n} and internal variable :math:p_n state and  Jeremy BLEYER committed Apr 25, 2018 176 an increment of total deformation :math:\Delta \boldsymbol{\varepsilon}. An elastic  177 178 179 trial stress :math:\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon} is first computed. The plasticity criterion is then evaluated with the previous plastic strain :math:f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n where  Jeremy BLEYER committed Apr 25, 2018 180 181 182 :math:\sigma^{eq}_{\text{elas}} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}} with the deviatoric elastic stress :math:\boldsymbol{s} = \operatorname{dev}\boldsymbol{\sigma}_{\text{elas}}. If :math:f_{\text{elas}} < 0, no plasticity occurs during this time increment and  183 :math:\Delta p,\Delta \boldsymbol{\varepsilon}^p =0.  Jeremy BLEYER committed Apr 25, 2018 184 185 186 187  Otherwise, plasticity occurs and the increment of plastic strain is given by :math:\Delta p = \dfrac{f_{\text{elas}}}{3\mu+H}. Hence, both elastic and plastic evolution can be accounted for by defining the  188 plastic strain increment as follows:  Jeremy BLEYER committed Apr 25, 2018 189 190 191 192 193 194  .. math:: \Delta p = \dfrac{\langle f_{\text{elas}}\rangle_+}{3\mu+H} where :math:\langle \star \rangle_+ represents the positive part of :math:\star and is obtained by function ppos. Plastic evolution also requires the computation  195 of the normal vector to the final yield surface given by  Jeremy BLEYER committed Apr 25, 2018 196 :math:\boldsymbol{n}_{\text{elas}} = \boldsymbol{s}/\sigma_{\text{elas}}^{eq}. In the following,  197 this vector must be zero in case of elastic evolution. Hence, we multiply it by  Jeremy BLEYER committed Apr 25, 2018 198 :math:\dfrac{\langle f_{\text{elas}}\rangle_+}{ f_{\text{elas}}} to tackle  199 200 201 202 both cases in a single expression. The final stress state is corrected by the plastic strain as follows :math:\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} - \beta \boldsymbol{s} with :math:\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p. It can be observed that the last term vanishes in case of elastic evolution so  Jeremy BLEYER committed Apr 25, 2018 203 204 205 206 207 208 209 that the final stress is indeed the elastic predictor. :: ppos = lambda x: (x+abs(x))/2. def proj_sig(deps, old_sig, old_p): sig_n = as_3D_tensor(old_sig) sig_elas = sig_n + sigma(deps) s = dev(sig_elas)  210  sig_eq = sqrt(3/2.*inner(s, s))  Jeremy BLEYER committed Apr 25, 2018 211 212 213 214 215 216 217 218 219  f_elas = sig_eq - sig0 - H*old_p dp = ppos(f_elas)/(3*mu+H) n_elas = s/sig_eq*ppos(f_elas)/f_elas beta = 3*mu*dp/sig_eq new_sig = sig_elas-beta*s return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \ as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \ beta, dp  220 .. note::  Jeremy BLEYER committed Apr 25, 2018 221 222 223 224 225 226 227 228 229 230 231  We could have used conditionals to write more explicitly the difference between elastic and plastic evolution. In order to use a Newton-Raphson procedure to resolve global equilibrium, we also need to derive the algorithmic consistent tangent matrix given by: .. math:: \mathbf{C}_{\text{tang}}^{\text{alg}} = \mathbf{C} - 3\mu\left(\dfrac{3\mu}{3\mu+H}-\beta\right) \boldsymbol{n}_{\text{elas}} \otimes \boldsymbol{n}_{\text{elas}} - 2\mu\beta\mathbf{DEV} where :math:\mathbf{DEV} is the 4th-order tensor associated with the deviatoric  232 233 234 235 operator (note that :math:\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C} for elastic evolution). Contrary to what is done in the FEniCS book, we do not store it as the components of a 4th-order tensor but it will suffice keeping track of the normal vector and the :math:\beta parameter related to the plastic strains. We instead define a function  Jeremy BLEYER committed Apr 25, 2018 236 237 238 239 240 241 242 243 244 245 246 247 248 249 computing the tangent stress :math:\boldsymbol{\sigma}_\text{tang} = \mathbf{C}_{\text{tang}}^{\text{alg}} \boldsymbol{\varepsilon} as follows:: def sigma_tang(e): N_elas = as_3D_tensor(n_elas) return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e) -------------------------------------------- Global problem and Newton-Raphson procedure -------------------------------------------- We are now in position to derive the global problem with its associated Newton-Raphson procedure. Each iteration will require establishing equilibrium  250 by driving to zero the residual between the internal forces associated with the current  Jeremy BLEYER committed Apr 25, 2018 251 252 253 254 255 256 stress state sig and the external force vector. Because we use Quadrature elements a custom integration measure must be defined to match the quadrature degree and scheme used by the Quadrature elements:: metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"} dxm = dx(metadata=metadata)  257   Jeremy BLEYER committed Apr 25, 2018 258 259 260 261  a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)  262 The consitutive update defined earlier will perform nonlinear operations on  Jeremy BLEYER committed Apr 25, 2018 263 264 265 266 the stress and strain tensors. These nonlinear expressions must then be projected back onto the associated Quadrature spaces. Since these fields are defined locally in each cell (in fact only at their associated Gauss point), this projection can be performed locally. For this reason, we define a local_project function  Jeremy BLEYER committed Apr 25, 2018 267 that use the LocalSolver to gain in efficiency (see also :ref:TipsTricksProjection)  Jeremy BLEYER committed Apr 25, 2018 268 for more details::  269   Jeremy BLEYER committed Apr 25, 2018 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286  def local_project(v, V, u=None): dv = TrialFunction(V) v_ = TestFunction(V) a_proj = inner(dv, v_)*dxm b_proj = inner(v, v_)*dxm solver = LocalSolver(a_proj, b_proj) solver.factorize() if u is None: u = Function(V) solver.solve_local_rhs(u) return u else: solver.solve_local_rhs(u) return .. note:: We could have used the standard project if we are not interested in optimizing  287  the code. However, the use of Quadrature elements would have required telling  Jeremy BLEYER committed Apr 25, 2018 288  project to use an appropriate integration measure to solve the global :math:L^2  289  projection that occurs under the hood. This would have needed either redefining  Jeremy BLEYER committed Apr 25, 2018 290 291 292 293 294  explicitly the projection associated forms (as we just did) or specifiying the appropriate quadrature degree to the form compiler as follows :code:project(sig_, W, form_compiler_parameters={"quadrature_degree":deg_stress}) Before defining the Newton-Raphson loop, we set up the output file and appropriate  295 FunctionSpace (here piecewise constant) and Function for output of the equivalent  Jeremy BLEYER committed Apr 25, 2018 296 plastic strain since XDMF output does not handle Quadrature elements::  297   Jeremy BLEYER committed Apr 25, 2018 298 299 300 301 302 303 304 305 306 307 308 309 310  file_results = XDMFFile("plasticity_results.xdmf") file_results.parameters["flush_output"] = True file_results.parameters["functions_share_mesh"] = True P0 = FunctionSpace(mesh, "DG", 0) p_avg = Function(P0, name="Plastic strain") We now define the global Newton-Raphson loop. We will discretize the applied loading using Nincr increments from 0 up to 1.1 (we exclude zero from the list of load steps). A nonlinear discretization is adopted to refine the load steps during the plastic evolution phase. At each time increment, the system is assembled and the residual norm is computed. The incremental displacement Du is initialized to zero and the inner iteration loop performing the constitutive  311 update is initiated. Inside this loop, corrections du to the displacement  Jeremy BLEYER committed Apr 25, 2018 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 increment are computed by solving the Newton system and the return mapping update is performed using the current total strain increment deps. The resulting quantities are then projected onto their appropriate FunctionSpaces. The Newton system and residuals are reassembled and this procedure continues until the residual norm falls below a given tolerance. After convergence of the iteration loop, the total displacement, stress and plastic strain states are updated :: Nitermax, tol = 200, 1e-8 # parameters of the Newton-Raphson procedure Nincr = 20 load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5 results = np.zeros((Nincr+1, 2)) for (i, t) in enumerate(load_steps): loading.t = t A, Res = assemble_system(a_Newton, res, bc) nRes0 = Res.norm("l2") nRes = nRes0 Du.interpolate(Constant((0, 0))) print("Increment:", str(i+1)) niter = 0 while nRes/nRes0 > tol and niter < Nitermax: solve(A, du.vector(), Res, "mumps") Du.assign(Du+du) deps = eps(Du) sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p) local_project(sig_, W, sig) local_project(n_elas_, W, n_elas) local_project(beta_, W0, beta) A, Res = assemble_system(a_Newton, res, bc) nRes = Res.norm("l2") print(" Residual:", nRes) niter += 1 u.assign(u+Du) sig_old.assign(sig) p.assign(p+local_project(dp_, W0))  346   Jeremy BLEYER committed Apr 25, 2018 347 348 349 ---------------- Post-processing ----------------  350   Jeremy BLEYER committed Apr 25, 2018 351 352 Inside the incremental loop, the displacement and plastic strains are exported at each time increment, the plastic strain must first be projected onto the  353 previously defined DG FunctionSpace. We also monitor the value of the cylinder  Jeremy BLEYER committed Apr 25, 2018 354 355 356 357 358 359 displacement on the inner boundary. The load-displacement curve is then plotted:: file_results.write(u, t) p_avg.assign(project(p, P0)) file_results.write(p_avg, t) results[i+1, :] = (u(Ri, 0)[0], t)  360   Jeremy BLEYER committed Apr 25, 2018 361 362 363 364 365 366 367 368 369 370 371  import matplotlib.pyplot as plt plt.plot(results[:, 0], results[:, 1], "-o") plt.xlabel("Displacement of inner boundary") plt.ylabel(r"Applied pressure $q/q_{lim}$") plt.show() The load-displacement curve looks as follows: .. image:: cylinder_expansion_load_displ.png :scale: 20%  372 It can also be checked that the analytical limit load is also well reproduced  Jeremy BLEYER committed Apr 25, 2018 373 374 375 376 377 378 when considering a zero hardening modulus. ----------- References -----------  379 .. [BON2014] Marc Bonnet, Attilio Frangi, Christian Rey.  Jeremy BLEYER committed Apr 25, 2018 380 381  *The finite element method in solid mechanics.* McGraw Hill Education, pp.365, 2014