diff --git a/Linear_Elastoplasticity/CMakeLists.txt b/Linear_Elastoplasticity/CMakeLists.txt new file mode 100644 index 00000000..44168d3d --- /dev/null +++ b/Linear_Elastoplasticity/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the step-8 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "linear_elastoplasticity") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 9.1 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/Linear_Elastoplasticity/README.md b/Linear_Elastoplasticity/README.md new file mode 100644 index 00000000..0fb5b4f7 --- /dev/null +++ b/Linear_Elastoplasticity/README.md @@ -0,0 +1,68 @@ +## Overview +There are some well-established experimental procedures that are commonly employed to investigate the plastic reponse of metals and similar materials. They are typically setup such that the stress and strain field in a region of the domain is predictable, which means that (i) the experiments are reliable and reproducable, and (ii) that can be more easily correlated to microstructural studies and numerical simulations, or vice versa. This is particularly useful when the materials exhibit complex constitutive responses (such as anisotropic plastic behaviour, or viscoplastic behaviour), as a well defined set of experiments can be used to understand the nature of the material response on both the macroscopic and microscopic levels. + +In this code we have implemented basic linear elasticity with linear isotropic hardening. Although this is a very simple model of plasticity, it is useful as a stepping-stone to fortify one's understanding of the steps necessary to implement more elaborate models. The constiutive model(s) are used to examine the behaviour of two geometries, namely a notched cylindrical specimen and a tensile specimen that conforms to a certified standard. + +## Theory +The theory of linear elastoplasticity, starting from the absolute fundamentals moving all the way to the numerical model, involved several distinct topic that deserve thorough explanation. This includes: +* the governing equations of quasi-static elasticity, +* a generalised continuum framework for elastoplasticity, +* constiutive laws for elasticity and elastoplasticity with linear isotropic hardening (including the derivation of the elastoplastic tangent modulus in the continuous setting), +* time integration algorithms for elasto-plasticity (leading to the time-discrete setting and algorithmically correct elastoplastic tangent modulus), and finally +* the finite element discretisation itself. +Since its easier to do so, we have provided a [summary of the fundamental theory](./doc/theory/theory-linear_elastoplasticity.pdf) in the supplied PDF document. Unfortunately, though, the reader will have to satisfied with a marginally condensed (bullet point) form of the elastoplastic framework, constitutive laws and integration algorithms. + +## Recommended Literature +* J. C. Simo, and T. J. R. Hughes (2000), Computational Inelasticity. Springer Science & Business Media. ISBN: 978-0-387-97520-7. DOI: [10.1007/b98904](http://doi.org/10.1007/b98904) +* O. C. Zienkiewicz and R. L. Taylor (2005), The finite element method for solid and structural mechanics. Butterworth-Heinemann. ISBN: 978-1-85617-634-7. DOI: [10.1016/C2009-0-26332-X](https://doi.org/10.1016/C2009-0-26332-X) +* B. D. Reddy (2016) Theoretical and Numerical Elastoplasticity. In: J. Schroeder, P. Wriggers (eds) Advanced Finite Element Technologies. CISM International Centre for Mechanical Sciences (Courses and Lectures), vol 566, pp 177-194. Springer, Cham. ISBN: 978-3-319-31923-0. DOI: [10.1007/978-3-319-31925-4_7](https://doi.org/10.1007/978-3-319-31925-4_7) + +## Compiling and running +Similar to the example programs, run +``` +cmake -DDEAL_II_DIR=/path/to/deal.II . +``` +in this directory to configure the problem. +You can switch between debug and release mode by calling either +``` +make debug +``` +or +``` +make release +``` +The problem may then be run with +``` +make run +``` + +## Problem setup and results +In this code we have two toy test cases set up to demonstrate either purely elastic (fully recoverable) material behaviour or elastoplasticity (dissipative) with linear isotropic hardening. +The constitutive parameters are entirely fictitious, and were chosen such that plastic flow would be easily observed in both cases. + +### Tensile specimen +In the first example, we investigate plastic behaviour of a dumbbell-shaped tensile specimen, which is dimensioned according the DIN 53504 type S3 specification. The supplied geometry is shown in the figure below. +![Problem geometry](./doc/results/tensile_specimen-geometry.png) +Due to the symmetry of the geometry, only one-eigth of the geometry need be modelled. +Naturally, the cut planes around which the model is reflected must be treated as symmetry boundary conditions. Note that we also attach manifolds to the curved surface in order to maintain the accuracy of the geometry upon refinement. + +These specimens are typically clamped in a tensile test machine, and stretched according to some experimental protocol. Accordingly, a portion of the wider surface of the specimen has a prescribed displacement in the axial direction (aligned with the X-Cartesian basis vector), while the out of plane displacement over the same subregion of the boundary is not permitted. The second condition is a rudimentary approximation to the real physical conditions, as the thickness of the specimen would change during clamping and during the course of the experiment itself. + +The next figure shows the geometry after it is lengthened by 20%. The plotted contours depict the local values of the isotropic hardening variable, which is an indicator of how much (non-recoverable) plastic flow has occured. +![Displaced solution](./doc/results/tensile_specimen-isotropic_hardening.png) +There are three regions where distinct and interesting behaviour may be noted. The first is the thin part of the specimen, where necking occurs. The strain in this region is nearly constant (which is precisely the intention of the design of the specimen); this leads to near uniform plastic flow in this region, which is ideal for experimental observations. The second interesting region is where the specimen starts to thicken. The change in curvature of the geometry causes a slight stress concentration in this location, leading to increased plastic flow. Again, by design this occurs out of the region of experimental observation. The third place where plastic flow is observed is near the clamped surface. This is a numerical artefact caused by the simplified boundary conditions. In reality, the thickness of the sample would reduce here and the amount of out-of-plane deformation would be less than that predicted here. + +### Notched cylinder +In the second example we model cylindrical specimen that is also placed in tension. However, a defect in the form of a notch is introduced into the specimen, with the idea that this will induce a localised stress concentration that will, in turn, lead to plastic flow. +The geometry is shown below; we may again apply symmetry arguments to reduce the size of the problem to be modelled. In this case (for which the geometry is created inside `deal.II`), we model only a half of the specimen. With a little effort, this could also be reduced to a one-eigth model. +![Problem geometry](./doc/results/notched_cylinder-geometry.png) + +Like before, the specimen is stretched along its axis such that its overall length increases by 20%. +As anticipated, there is indeed a significant amount of plastic flow at the reentrant corner of the notch where the stress concentration is predicted to be. +![Displaced solution](./doc/results/notched_cylinder-isotropic_hardening.png) +We also observe that the majority of the specimen is predicted to undergo some plastic flow, although the larger-radius section adjacent to the notch undergoes no plastic deformation. This could be a real phenomenon, but could also be an experimental artefact. Due to the refinement strategy chosen for this case, the mesh around the notch is more refined than in the bulk of the specimen. This means that the boundary of the bulk geometry is captured with differing accuracy as one looks towards the notch, and may artifically be inducing some non-uniform stress response.Further investigation and mesh refinement studies could render the answer to this outstanding question. + +## Future work / possibilities for extension +* Implement and investigate other hardening rules, such as kinematic hardening and anisotropic hardening. +* Many materials also exhibit rate-dependent responses, so using an elasto-viscoplastic constitutive law might be more appropriate in those cases. +* Due to the large deformation, it may be necessary to extend this model from the linear response regimen to one that incorporates the effects of geometric stiffening, i.e. finite strain elasto-plasticity. \ No newline at end of file diff --git a/Linear_Elastoplasticity/doc/results/notched_cylinder-geometry.png b/Linear_Elastoplasticity/doc/results/notched_cylinder-geometry.png new file mode 100644 index 00000000..117583df Binary files /dev/null and b/Linear_Elastoplasticity/doc/results/notched_cylinder-geometry.png differ diff --git a/Linear_Elastoplasticity/doc/results/notched_cylinder-isotropic_hardening.png b/Linear_Elastoplasticity/doc/results/notched_cylinder-isotropic_hardening.png new file mode 100644 index 00000000..3221b4f1 Binary files /dev/null and b/Linear_Elastoplasticity/doc/results/notched_cylinder-isotropic_hardening.png differ diff --git a/Linear_Elastoplasticity/doc/results/tensile_specimen-geometry.png b/Linear_Elastoplasticity/doc/results/tensile_specimen-geometry.png new file mode 100644 index 00000000..78b4748c Binary files /dev/null and b/Linear_Elastoplasticity/doc/results/tensile_specimen-geometry.png differ diff --git a/Linear_Elastoplasticity/doc/results/tensile_specimen-isotropic_hardening.png b/Linear_Elastoplasticity/doc/results/tensile_specimen-isotropic_hardening.png new file mode 100644 index 00000000..84cc7ce9 Binary files /dev/null and b/Linear_Elastoplasticity/doc/results/tensile_specimen-isotropic_hardening.png differ diff --git a/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_1-Radial_projection.png b/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_1-Radial_projection.png new file mode 100644 index 00000000..f6e94e99 Binary files /dev/null and b/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_1-Radial_projection.png differ diff --git a/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_2-General_closest_point_projection.png b/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_2-General_closest_point_projection.png new file mode 100644 index 00000000..5baae29a Binary files /dev/null and b/Linear_Elastoplasticity/doc/theory/Figures/Sec_3_5_2-General_closest_point_projection.png differ diff --git a/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.bib b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.bib new file mode 100644 index 00000000..024e7d4b --- /dev/null +++ b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.bib @@ -0,0 +1,21 @@ +%% LaTeX2e file `theory-linear_elastoplasticity.bib' +%% generated by the `filecontents' environment +%% from source `theory-linear_elastoplasticity' on 2020/05/25. +%% +@Misc{Mergheim2018a, + author = {Julia Mergheim and Jean-Paul Pelteret and Ester Comellas}, + title = {Lecture notes: Material Modelling and Simulation}, + year = {2018}, + publisher = {Friedrich--Alexander Uniersit\"{a}t Erlangen--N\"{u}rnberg} +} +@Book{Simo2006a, + title = {Computational Inelasticity}, + publisher = {Springer Science \& Business Media}, + year = {2000}, + author = {Simo, J. C. and Hughes, Thomas J. R.}, + editor = {Marsden, J. E. and Sirovich, L. and Wiggins, S.}, + volume = {7}, + series = {Interdisciplinary Applied Mathematics}, + isbn = {978-0-387-97520-7}, + doi = {10.1007/b98904}, +} diff --git a/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.pdf b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.pdf new file mode 100644 index 00000000..2947b9d2 Binary files /dev/null and b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.pdf differ diff --git a/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.tex b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.tex new file mode 100644 index 00000000..bea7a5d1 --- /dev/null +++ b/Linear_Elastoplasticity/doc/theory/theory-linear_elastoplasticity.tex @@ -0,0 +1,803 @@ +\documentclass[]{scrartcl} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{graphicx} +\usepackage{hyperref} +\usepackage[sharp]{easylist} +\usepackage{cleveref} +\usepackage{subcaption} +\usepackage{filecontents} +\usepackage[sort]{natbib} + +\title{Theory: Linear elastoplasicity} +\author{Jean-Paul Pelteret} + +\begin{filecontents}{\jobname.bib} +@Misc{Mergheim2018a, + author = {Julia Mergheim and Jean-Paul Pelteret and Ester Comellas}, + title = {Lecture notes: Material Modelling and Simulation}, + year = {2018}, + publisher = {Friedrich--Alexander Uniersit\"{a}t Erlangen--N\"{u}rnberg} +} +@Book{Simo2006a, + title = {Computational Inelasticity}, + publisher = {Springer Science \& Business Media}, + year = {2000}, + author = {Simo, J. C. and Hughes, Thomas J. R.}, + editor = {Marsden, J. E. and Sirovich, L. and Wiggins, S.}, + volume = {7}, + series = {Interdisciplinary Applied Mathematics}, + isbn = {978-0-387-97520-7}, + doi = {10.1007/b98904}, +} +\end{filecontents} + +\begin{document} + +\maketitle + +\begin{abstract} +An introduction to the theory applied for elastoplasticity. +\end{abstract} + +\section{Governing equations for quasi-static linear elasticity} +The strong statement of the balance of linear momentum reads +\begin{gather} +\nabla \cdot \boldsymbol{\sigma} + \mathbf{b} + = \mathbf{0} +\quad \text{on} \quad \Omega \quad , +\end{gather} +where $\nabla = \frac{\partial}{\partial x}$ is a differential operator, +$\boldsymbol{\sigma}$ is the Cauchy stress tensor and +$\mathbf{b} = \rho \mathbf{g}$ is the body force density vector. +This is expressed in index notation as +\begin{gather} +\frac{\partial \sigma_{ij}}{\partial x_{j}} + b_{i} + = 0 +\quad \text{on} \quad \Omega \quad . +\end{gather} +Pre-multiplying the above by test function $\delta \mathbf{v}$ and integrating over the domain $\Omega$ renders +\begin{gather} +- \int\limits_{\Omega} \delta v_{i} \, \frac{\partial \sigma_{ij}}{\partial x_{j}} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv +\end{gather} +that, by using the product rule for derivatives (i.e. integration by parts), becomes +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \sigma_{ij} \, dv +- \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \delta v_{i} \, \sigma_{ij} \right] \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv +\quad . +\end{gather} +Finally, by applying divergence theorem to the second term in the above, we attain the weak form of the balance of linear momentum +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \sigma_{ij} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + + \int\limits_{\partial\Omega} \delta v_{i} \, \underbrace{\sigma_{ij} \, n_{j}}_{\bar{t}_{i}} \, da +\quad , +\label{equ: Weak form of the balance of linear momentum} +\end{gather} +wherein $\mathbf{n}$ represents the outward facing normal on $\partial\Omega$, the boundary of the domain, and $\bar{\mathbf{t}}$ the prescribed traction on the Neumann boundary. + +\clearpage +\section{General framework for elastoplasticity} +\begin{easylist}[itemize] +# Kinematics +## Small strain tensor +\begin{gather} +\boldsymbol{\varepsilon} + = \frac{1}{2} \left[ \nabla \mathbf{u} + \left[ \nabla \mathbf{u} \right]^{T} \right] + = \nabla\mathbf{u}^{sym} +\end{gather} +## Additive volumetric-isochoric split of the strain +\begin{gather} +\boldsymbol{\varepsilon} + = \boldsymbol{\varepsilon}^{vol} + \boldsymbol{\varepsilon}^{dev} +\quad ; \quad +\boldsymbol{\varepsilon}^{vol} = \dfrac{1}{3} \textrm{ tr}(\boldsymbol{\varepsilon}) \mathbf{I} +\end{gather} +## Additive split of strain into elastic and plastic contributions +\begin{gather*} +\boldsymbol{\varepsilon} + = \boldsymbol{\varepsilon}^{e} + \boldsymbol{\varepsilon}^{p} +\end{gather*} +# General free energy function +\begin{gather*} +\psi = \psi (\boldsymbol{\varepsilon}^{e}, \alpha, \boldsymbol{\beta}) +\end{gather*} +## Two new internal variables: +### $\alpha$ describes the relative increase of the elastic region (isotropic hardening variable) +### $\boldsymbol{\beta}$ describes the kinematic hardening (a rank-2 tensor; related to the centre point of the elastic region and introduces anisotropy due to plastic flow) +# Principle of irreversibility: +## Dissipation inequality +\begin{align*} +\mathcal{D} + &= \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} - \dot\psi (\boldsymbol{\varepsilon}^{e}, \alpha, \boldsymbol{\beta}) \\ + &= \boldsymbol{\sigma} : \left[ \dot{\boldsymbol{\varepsilon}}^{e} + \dot{\boldsymbol{\varepsilon}}^{p} \right] + - \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^{e}} : \dot{\boldsymbol{\varepsilon}^{e}} + \underbrace{- \frac{\partial \psi}{\partial \alpha}}_{R} \dot{\alpha} + \underbrace{- \frac{\partial \psi}{\partial \boldsymbol{\beta}}}_{\boldsymbol{B}} : \dot{\boldsymbol{\beta}} \\ + &= \left[ \boldsymbol{\sigma} - \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^{e}} \right] : \dot{\boldsymbol{\varepsilon}^{e}} + + \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}}^{p} + + R \dot{\alpha} + + \boldsymbol{B} : \dot{\boldsymbol{\beta}} + \\ + &\geq 0 +\end{align*} +### Hardening stress $R$ +### Back stress $\boldsymbol{B}$ +## Definition of Cauchy stress +\begin{gather} +\boldsymbol{\sigma} + := \frac{\partial \psi}{\partial \boldsymbol{\varepsilon}^{e}} +\end{gather} +## Reduced dissipation inequality (using definition of Cauchy stress) +\begin{gather*} +\mathcal{D}^{red} + = \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}}^{p} + + R \dot{\alpha} + + \boldsymbol{B} : \dot{\boldsymbol{\beta}} + \geq 0 +\end{gather*} +# Evolution equations derived by the postulate of maximum dissipation [note: this ensures that hardening occurs when accommodated by the model, thus dissipating more energy than perfect plastic flow] +## Elastic region is restricted by yield function $\Phi \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right)$: +\begin{gather} +\mathbb{E} := \left\lbrace \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \vert \Phi \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \leq 0 \right\rbrace +\end{gather} +## Conditions: +### Elastic state: $\Phi < 0$ +### Plastic state: $\Phi = 0$ +# Define Lagrange multiplier problem to maximise the dissipation, i.e. +\begin{gather*} +\textnormal{maximise } \mathcal{D}^{red} \textnormal{ subject to } \Phi \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \leq 0 \\ +\mathcal{L} \left( \boldsymbol{\sigma},R,\boldsymbol{B}, \dot{\lambda} \right) + = -\mathcal{D}^{red} + \dot{\lambda} \Phi + \quad \rightarrow \quad \textnormal{stationary} +\end{gather*} +### $\dot{\lambda}$ is known as the consistency parameter +## Taking directional derivatives, we can identify the following evolution equations for the internal variables: +\begin{align*} +\frac{\partial \mathcal{L}}{\partial \boldsymbol{\sigma}} = \boldsymbol{0} +\quad &\Rightarrow \quad +\dot{\boldsymbol{\varepsilon}}^{p} = \dot{\lambda} \frac{\partial \Psi}{\partial \boldsymbol{\sigma}} +\quad \textnormal{(flow rule)} +\\ +\frac{\partial \mathcal{L}}{\partial R} = 0 +\quad &\Rightarrow \quad +\dot{\alpha} = \dot{\lambda} \frac{\partial \Psi}{\partial R} +\quad \textnormal{(evolution equation for isotropic hardening)} +\\ +\frac{\partial \mathcal{L}}{\partial \boldsymbol{B}} = \boldsymbol{0} +\quad &\Rightarrow \quad +\dot{\boldsymbol{\beta}} = \dot{\lambda} \frac{\partial \Psi}{\partial \boldsymbol{B}} +\quad \textnormal{(evolution equation for kinematic hardening)} +\\ +%\frac{\partial \mathcal{L}}{\partial \dot{\lambda}} +% = \\ +\end{align*} +% http://www.engr.colostate.edu/~thompson/hPage/CourseMat/Courses/CE696DV/plasticity.pdf +%https://www.ethz.ch/content/dam/ethz/special-interest/baug/ifb/ifb-dam/homepage-IfB/Education/msc_courses/msc-mechanics-building-materials/documents/plastizitaet.pdf +% p 36 +## Loading / unloading condition: +\begin{align*} +\frac{\partial \Psi}{\partial \boldsymbol{\sigma}} \dot{\boldsymbol{\sigma}} +\begin{cases} +< 0 \quad & \textnormal{Unloading (elastic)} \\ += 0 \quad & \textnormal{Loading (plastic flow)} \\ +> 0 \quad & \textnormal{Loading (plastic hardening)} +\end{cases} +%\textnormal{Loading (plastic hardening)} &\frac{\partial \Psi}{\partial \boldsymbol{\sigma}} \dot{\boldsymbol{\sigma}} > 0 \\ +%\textnormal{Loading (plastic flow): } &\frac{\partial \Psi}{\partial \boldsymbol{\sigma}} \dot{\boldsymbol{\sigma}} = 0 \\ +%\textnormal{Unloading: } &\frac{\partial \Psi}{\partial \boldsymbol{\sigma}} \dot{\boldsymbol{\sigma}} < 0 +\end{align*} +### The loading-unloading conditions can be shown to be equivalent to the Karush-Kuhn-Tucker (KKT) conditions (an important statement in optimisation theory): +% https://en.wikipedia.org/wiki/Flow_plasticity_theory +% Simo, Computational elasticity, p77 +\begin{gather} +\dot{\lambda} \geq 0 +\quad;\quad +\Phi \leq 0 +\quad;\quad +\dot{\lambda} \Phi = 0 +\end{gather} +#### Situations and definitions for any $\left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \in \mathbb{E}$ +\begin{align*} +\Phi < 0 + \quad &\Rightarrow \quad + \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \in \textnormal{int}(\mathbb{E}) + \quad \Rightarrow \quad + \dot{\lambda} = 0 \quad \textnormal{(Elastic response)} \\ +\Phi = 0 + \quad &\Rightarrow \quad + \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) \in \partial\mathbb{E} + \quad \Rightarrow \quad + \begin{cases} + \dot{\Phi} < 0 \quad \Rightarrow \quad \dot{\lambda} = 0 \quad \textnormal{(Elastic unloading)} \\ + \dot{\Phi} = 0 \quad \textnormal{and} \quad \dot{\lambda} = 0 \quad \textnormal{(Neutral loading)} \\ + \dot{\Phi} = 0 \quad \textnormal{and} \quad \dot{\lambda} > 0 \quad \textnormal{(Plastic loading)} + \end{cases} +\end{align*} +## Hardening / softening condition: +\begin{align*} +\frac{\partial \Psi}{\partial \alpha} \dot{\alpha} +\begin{cases} +< 0 \quad & \textnormal{Softening} \\ +> 0 \quad & \textnormal{Hardening} +\end{cases} +%\textnormal{Softening: } &\frac{\partial \Psi}{\partial \alpha} \dot{\alpha} > 0 \\ +%\textnormal{Hardening: } &\frac{\partial \Psi}{\partial \alpha} \dot{\alpha} < 0 +\end{align*} +# Example of Von-Mises yield criterion +## Convex yield surface defined by +\begin{gather*} +\Phi \left( \boldsymbol{\sigma} \right) + = \sqrt{\frac{3}{2}} \Vert \boldsymbol{\sigma}^{dev} \Vert - \sigma_{y} + \leq 0 +\end{gather*} +### Cylindrical yield-surface as observed in the principal stress space +### Norm of deviatoric stress contribution: +\begin{gather*} +\Vert \boldsymbol{\sigma}^{dev} \Vert + = \sqrt{\boldsymbol{\sigma}^{dev} : \boldsymbol{\sigma}^{dev}} +\end{gather*} +\end{easylist} + +\clearpage +\section{Constitutive law: Linear elasticity} + +\begin{easylist}[itemize] +# Linear elasticity (a special case of hyperelasticity) +## Dissipation (in-)equality +\begin{gather} +\mathcal{D} + = \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} - \dot\psi \left( \boldsymbol{\varepsilon} \right) + = 0 +\quad \Rightarrow \quad +\boldsymbol{\sigma} + = \frac{\partial \psi \left( \boldsymbol{\varepsilon} \right)}{\partial \boldsymbol{\varepsilon}} +\end{gather} +## Free/Strain energy function (Hooke's law) +\begin{gather} +\psi(\boldsymbol{\varepsilon}) + = \dfrac{\lambda}{2} \left[\textrm{tr}(\boldsymbol{\varepsilon})\right]^2 + + \mu \textrm{ tr}(\boldsymbol{\varepsilon}^2) +\end{gather} +## Cauchy stress +\begin{gather} +\boldsymbol{\sigma}(\boldsymbol{\varepsilon}) = \dfrac{\partial \psi}{\partial \boldsymbol{\varepsilon}} = \lambda \textrm{ tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2\mu\boldsymbol{\varepsilon} +\end{gather} +## Elastic tangent +\begin{gather} +\mathbb{C} = \dfrac{{\partial}^2 \psi}{\partial \boldsymbol{\varepsilon}^2}: \dot{\boldsymbol{\varepsilon}} = \lambda \mathbf{I} \otimes \mathbf{I} + 2\mu\;\mathbb{I}^{sym} +\end{gather} +### with the linear relationship +\begin{gather} +\boldsymbol{\sigma} + = \boldsymbol{\mathbb{C}} : \boldsymbol{\varepsilon} +\end{gather} +\end{easylist} + +\clearpage +\section{Constitutive law: Elastoplasticity} + +The linear elastoplastic constitutive laws and framework described here are described in detail by \citep{Simo2006a,Mergheim2018a}. + +\subsection{Associated Von Mises elastoplasticity with linear isotropic and kinematic hardening} +\begin{easylist} +# Kinematics +\begin{gather*} +\varepsilon + = \varepsilon^{e} + \varepsilon^{p} +\end{gather*} +# Free energy (additive decomposition into elastic and plastic parts) +\begin{align*} +\psi (\boldsymbol{\varepsilon}^{e}, \alpha, \boldsymbol{\beta}) + &= \psi^{e} (\boldsymbol{\varepsilon}^{e}) + + \psi^{p} (\alpha, \boldsymbol{\beta}) \\ +\psi^{e} (\boldsymbol{\varepsilon}^{e}) + &= \frac{\lambda}{2} \left[ \textnormal{tr} (\boldsymbol{\varepsilon}^{e}) \right]^{2} + + \mu \textnormal{tr} \left( \left[\boldsymbol{\varepsilon}^{e} \right]^{2} \right) \\ +\psi^{p} (\alpha, \boldsymbol{\beta}) + &= \frac{1}{2} k r \alpha^{2} + \frac{1}{2} k \left[ 1 - r \right] \Vert \boldsymbol{\beta} \Vert^{2} + \quad , \quad r \in \left[ 0,1 \right] +\end{align*} +# Stresses +\begin{align*} +\boldsymbol{\sigma} + &= \dfrac{\partial \psi}{\partial \boldsymbol{\varepsilon^{e}}} + = \lambda \textrm{ tr}(\boldsymbol{\varepsilon}^{e}) \mathbf{1} + 2\mu\boldsymbol{\varepsilon}^{e} \\ +R + &= - \dfrac{\partial \psi}{\partial \alpha} + = - k r \alpha \\ +\boldsymbol{B} + &= - \dfrac{\partial \psi}{\partial \boldsymbol{\beta}} + = -k \left[ 1 - r \right] \boldsymbol{\beta} +\end{align*} +# Yield function (Von Mises) +\begin{gather*} +\Phi \left( \boldsymbol{\sigma},R,\boldsymbol{B} \right) + = \Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert - \sqrt{\frac{2}{3}} \left[ \sigma_{y} - R \right] + \leq 0 +\end{gather*} +# Evolution equations +\begin{align*} +\dot{\boldsymbol{\varepsilon}}^{p} + &= \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} + = \dot{\lambda} \dfrac{\boldsymbol{\sigma}^{dev} + \boldsymbol{B}}{\Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert} + = \dot{\lambda} \boldsymbol{n} +\\ +\dot{\alpha} + &= \dot{\lambda} \frac{\partial \Phi}{\partial R} + = \sqrt{\frac{2}{3}} \dot{\lambda} +\\ +\dot{\boldsymbol{\beta}} + &= \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{B}} + = \dot{\lambda} \dfrac{\boldsymbol{\sigma}^{dev} + \boldsymbol{B}}{\Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert} + = \dot{\lambda} \boldsymbol{n} +\end{align*} +## Must satisfy KTT conditions +\begin{gather*} +\dot{\lambda} \geq 0 +\quad;\quad +\Phi \leq 0 +\quad;\quad +\dot{\lambda} \Phi = 0 +\end{gather*} +# Principle of irreversibility +\begin{align*} +\mathcal{D}^{red} + &= \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}}^{p} + + R \dot{\alpha} + + \boldsymbol{B} : \dot{\boldsymbol{\beta}} \\ + &= \boldsymbol{\sigma} : \dot{\lambda} \boldsymbol{n} + + R \sqrt{\frac{2}{3}} \dot{\lambda} + + \boldsymbol{B} : \dot{\lambda} \boldsymbol{n} \\ + &= \dot{\lambda} \left[ \left[ + \boldsymbol{\sigma} + \boldsymbol{B} \right] : \boldsymbol{n} + + R \sqrt{\frac{2}{3}} + \right] \\ + &= \dot{\lambda} \left[ \left[ + \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \right] : \dfrac{\boldsymbol{\sigma}^{dev} + \boldsymbol{B}}{\Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert} + + R \sqrt{\frac{2}{3}} + \right] \\ + &= \dot{\lambda} \left[ \ + \Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert + + R \sqrt{\frac{2}{3}} + \right] \\ + &\geq 0 +\end{align*} +## If $\Phi < 0 \quad \Rightarrow \quad \dot{\lambda} = 0$ +## If $\Phi = 0 \quad \Rightarrow \quad \dot{\lambda} > 0$ + and $\Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert + R \sqrt{\dfrac{2}{3}} = \sqrt{\dfrac{2}{3}} \sigma_{y} \geq 0$ +# Consistency condition (evolution of yield surface) +\begin{align*} +\dot{\Phi} + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + + \frac{\partial \Phi}{\partial R} \dot{R} + + \frac{\partial \Phi}{\partial \boldsymbol{B}} : \dot{\boldsymbol{B}} \\ + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \underbrace{ \left[ \dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon^{p}}} \right]}_{\dot{\boldsymbol{\varepsilon^{e}}}} + - \frac{\partial \Phi}{\partial R} \dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} \dot{\alpha} + - \frac{\partial \Phi}{\partial \boldsymbol{B}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} : \dot{\boldsymbol{\beta}} \\ +% &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dot{\boldsymbol{\varepsilon}} +% - \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dot{\boldsymbol{\varepsilon^{p}}} +% - \frac{\partial \Phi}{\partial R} \dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} \dot{\alpha} +% - \frac{\partial \Phi}{\partial \boldsymbol{B}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} : \dot{\boldsymbol{\beta}} \\ + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dot{\boldsymbol{\varepsilon}} + - \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \dot{\lambda} + - \frac{\partial \Phi}{\partial R} \dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} \frac{\partial \Phi}{\partial R} \dot{\lambda} + - \frac{\partial \Phi}{\partial \boldsymbol{B}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} : \frac{\partial \Phi}{\partial \boldsymbol{B}} \dot{\lambda} \\ + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dot{\boldsymbol{\varepsilon}} + - \dot{\lambda} \left[ + \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} + + \frac{\partial \Phi}{\partial R} \dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} \frac{\partial \Phi}{\partial R} + + \frac{\partial \Phi}{\partial \boldsymbol{B}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} : \frac{\partial \Phi}{\partial \boldsymbol{B}} + \right] \\ + &= 0 \quad \textnormal{(stay on yield surface during plastic deformation)} +\end{align*} +\begin{align*} +\Rightarrow \quad \dot{\lambda} + = \dfrac { \dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dot{\boldsymbol{\varepsilon}} }{ \dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} : \dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} + + \dfrac{\partial \Phi}{\partial R} \dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} \dfrac{\partial \Phi}{\partial R} + + \dfrac{\partial \Phi}{\partial \boldsymbol{B}} : \dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} : \dfrac{\partial \Phi}{\partial \boldsymbol{B}} } +\end{align*} +with +\begin{align*} +\dfrac{\partial^{2} \psi}{\partial \boldsymbol{\varepsilon^{e}} \otimes \partial \boldsymbol{\varepsilon^{e}}} + &= \mathbb{C} \\ +\dfrac{\partial^{2} \psi}{\partial \alpha \partial \alpha} + &= -k r \\ +\dfrac{\partial^{2} \psi}{\partial \boldsymbol{\beta} \otimes \partial \boldsymbol{\beta}} + &= -k \left[ 1-r \right] \mathbb{I}^{sym} \\ +\dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} + &= \dfrac{\partial \Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert}{\partial \left( \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \right)} : \dfrac{\partial \left( \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \right)}{\partial \boldsymbol{\sigma}} + = \dfrac{\boldsymbol{\sigma}^{dev} + \boldsymbol{B}}{\Vert \boldsymbol{\sigma}^{dev} + \boldsymbol{B} \Vert} : \mathbb{I}^{sym} + = \boldsymbol{n} \\ +\dfrac{\partial \Phi}{\partial R} + &= \sqrt{\frac{2}{3}} \\ +\dfrac{\partial \Phi}{\partial \boldsymbol{B}} + &= \boldsymbol{n} +\end{align*} +Then +\begin{align*} +\Rightarrow \quad \dot{\lambda} + &= \dfrac { \boldsymbol{n} : \mathbb{C} : \dot{\boldsymbol{\varepsilon}} }{ \boldsymbol{n} : \mathbb{C} : \boldsymbol{n} + + \sqrt{\frac{2}{3}} \left[ -k r \right] \sqrt{\frac{2}{3}} + + \boldsymbol{n} : \left[ -k \left[ 1-r \right] \mathbb{I}^{sym} \right] : \boldsymbol{n} } \\ + &= \dfrac { \boldsymbol{n} : \mathbb{C} : \dot{\boldsymbol{\varepsilon}} }{ \boldsymbol{n} : \mathbb{C} : \boldsymbol{n} + - \frac{2}{3} k r + -k \left[ 1-r \right] } \\ + &= \dfrac { \boldsymbol{n} : \mathbb{C} : \dot{\boldsymbol{\varepsilon}} }{ \boldsymbol{n} : \mathbb{C} : \boldsymbol{n} + + k \left[\frac{1}{3} r - 1 \right] } \\ +\end{align*} +\end{easylist} + +\subsection{Elasto-plastic tangent modulus (continuous setting)} +\begin{easylist} +# Goal +\begin{gather*} +\dot{\boldsymbol{\sigma}} + = \boldsymbol{\mathcal{C}}^{ep} : \dot{\boldsymbol{\varepsilon}} +\end{gather*} +# Stress rate +\begin{gather*} +\boldsymbol{\sigma} + = \frac{\partial \Psi}{\partial \boldsymbol{\varepsilon}^{e}} \\ +\Rightarrow \quad +\dot{\boldsymbol{\sigma}} + = \frac{\partial^{2} \Psi}{\partial \boldsymbol{\varepsilon}^{e} \otimes \partial \boldsymbol{\varepsilon}^{e}} : \dot{\boldsymbol{\varepsilon}}^{e} + = \boldsymbol{\mathcal{C}} : \left[ \dot{\boldsymbol{\varepsilon}} - \dot{\boldsymbol{\varepsilon}}^{p} \right] + = \boldsymbol{\mathcal{C}} : \left[ \dot{\boldsymbol{\varepsilon}} - \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \right] +\end{gather*} +# Remember: Plastic potential $\Phi \left( \boldsymbol{\sigma},\boldsymbol{S} \right)$ with evolution laws +$ +\dot{\boldsymbol{\varepsilon}}^{p} + = \dot{\lambda} \dfrac{\partial \Phi}{\partial \boldsymbol{\sigma}} +$, $ +\dot{\boldsymbol{d}} + = \dot{\lambda} \dfrac{\partial \Phi}{\partial \boldsymbol{S}} +$ +## Definition and derivatives of stress-like terms $\boldsymbol{S} = \left[ R,\boldsymbol{B} \right]$ with internal variables $\boldsymbol{d} = \left[ \alpha,\boldsymbol{\beta} \right]$ +\begin{gather*} +\boldsymbol{S} + = - \frac{\partial \Psi}{\partial \boldsymbol{d}} \\ +\dot{\boldsymbol{S}} + = - \frac{\partial^{2} \Psi}{\partial \boldsymbol{d} \otimes \partial \boldsymbol{d}} \circ \dot{\boldsymbol{d}} + = - \frac{\partial^{2} \Psi}{\partial \boldsymbol{d} \otimes \partial \boldsymbol{d}} \circ \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{S}} \\ +\end{gather*} +# Use consistency condition to compute $\dot{\lambda}$: +\begin{gather*} +\dot{\lambda} \dot{\Phi} = 0 +\quad \Rightarrow \quad +\dot{\Phi} = 0 \quad \textnormal{for plastic loading} +\end{gather*} +\begin{align*} +\dot{\Phi} + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \dot{\boldsymbol{\sigma}} + + \frac{\partial \Phi}{\partial \boldsymbol{S}} \circ \dot{\boldsymbol{S}} \\ + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \boldsymbol{\mathcal{C}} : \left[ \dot{\boldsymbol{\varepsilon}} - \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \right] + - \frac{\partial \Phi}{\partial \boldsymbol{S}} \circ \frac{\partial^{2} \Psi}{\partial \boldsymbol{d} \otimes \partial \boldsymbol{d}} \circ \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{S}} \\ + &= 0 \\ +\Rightarrow \quad +\dot{\lambda} + &= \frac{1}{D} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \boldsymbol{\mathcal{C}} : \dot{\boldsymbol{\varepsilon}} +\quad \textnormal{with} \\ +D + &= \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \boldsymbol{\mathcal{C}} : \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} + + \frac{\partial \Phi}{\partial \boldsymbol{S}} \circ \frac{\partial^{2} \Psi}{\partial \boldsymbol{d} \otimes \partial \boldsymbol{d}} \circ \frac{\partial \Phi}{\partial \boldsymbol{S}} +\end{align*} +# Elasto-plastic tangent +\begin{gather*} +\dot{\boldsymbol{\sigma}} + = \boldsymbol{\mathcal{C}}^{ep} : \dot{\boldsymbol{\varepsilon}} + = \underbrace{\left[ \boldsymbol{\mathcal{C}} - \frac{1}{D} \left[ \boldsymbol{\mathcal{C}} : \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} \right] \otimes \left[ \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} : \boldsymbol{\mathcal{C}} \right] \right]}_{\boldsymbol{\mathcal{C}}^{ep} \textnormal{ (symmetric)}} : \dot{\boldsymbol{\varepsilon}} +\end{gather*} +# Result for isotropic hardening +\begin{subequations} +\begin{align} +\boldsymbol{\mathcal{C}}^{ep, dev} + &= 2 \mu \boldsymbol{\mathcal{I}}^{dev} - \frac{\left[ 2\mu \right]^{2}}{2\mu + \frac{2}{3}k} \boldsymbol{n} \otimes \boldsymbol{n} \\ +\boldsymbol{\mathcal{C}}^{ep, vol} + &= \kappa \boldsymbol{1} \otimes \boldsymbol{1} + \label{equ: Constutitve law: Algorithmically correct tangent: Volumetric} +\end{align} +\end{subequations} +with +\begin{gather*} +\boldsymbol{n} + = \frac{\boldsymbol{\sigma}^{dev}}{\Vert \boldsymbol{\sigma}^{dev} \Vert} +\end{gather*} +\end{easylist} + +\section{Integration algorithms for elasto-plasticity} +\subsection{Radial return algorithm for von Mises plasticity with linear isotropic hardening} +\begin{easylist} +# Constitutive setting +\begin{gather*} +\boldsymbol{\varepsilon} + = \boldsymbol{\varepsilon}^{e} + \boldsymbol{\varepsilon}^{p} + = \boldsymbol{\varepsilon}^{vol} + \boldsymbol{\varepsilon}^{dev} +\\ +\psi + = \frac{1}{2} \kappa \left[ \textnormal{tr} \left( \boldsymbol{\varepsilon}^{vol} \right) \right]^{2} + + \mu \textnormal{tr} \left( \left[\boldsymbol{\varepsilon}^{e, dev} \right]^{2} \right) + + \frac{1}{2} k \alpha^{2} +\\ +\Phi \left( \boldsymbol{\sigma}, R \right) + = \Vert \boldsymbol{\sigma}^{dev} \Vert - \sqrt{\frac{2}{3}} \left[ \sigma_{y} - R \right] + \leq 0 +\\ +\dot{\boldsymbol{\varepsilon}}^{p} + = \dot{\lambda} \frac{\partial \Phi}{\partial \boldsymbol{\sigma}} +% = \dot{\lambda} \dfrac{\boldsymbol{\sigma}^{dev}}{\Vert \boldsymbol{\sigma}^{dev} \Vert} + = \dot{\lambda} \boldsymbol{n} +\quad , \quad +\boldsymbol{n} + = \dfrac{\boldsymbol{\sigma}^{dev}}{\Vert \boldsymbol{\sigma}^{dev} \Vert} +\\ +\dot{\alpha} + = \dot{\lambda} \frac{\partial \Phi}{\partial R} + = \sqrt{\frac{2}{3}} \dot{\lambda} +\end{gather*} +# Apply backward Euler method to evolution equations +\begin{gather*} +\boldsymbol{\varepsilon}^{p}_{n+1} + = \boldsymbol{\varepsilon}^{p}_{n} + + \Delta t \, \dot{\boldsymbol{\varepsilon}}^{p}_{n+1} + = \boldsymbol{\varepsilon}^{p}_{n} + + \Delta t \, \dot{\lambda}_{n+1} \boldsymbol{n}_{n+1} +\\ +\alpha_{n+1} + = \alpha_{n} + + \Delta t \, \dot{\alpha}_{n+1} + = \alpha_{n} + + \Delta t \, \dot{\lambda}_{n+1} \sqrt{\frac{2}{3}} +\end{gather*} +# Deviatoric stress +\begin{align*} +\boldsymbol{\sigma}^{dev}_{n+1} + & = 2 \mu \boldsymbol{\varepsilon}^{e, dev}_{n+1} + = 2 \mu \left[ \boldsymbol{\varepsilon}^{dev}_{n+1} - \boldsymbol{\varepsilon}^{p, dev}_{n+1} \right] \\ + &= 2 \mu \boldsymbol{\varepsilon}^{dev}_{n+1} + - 2 \mu \left[ \boldsymbol{\varepsilon}^{p}_{n} + + \Delta t \, \dot{\lambda}_{n+1} \boldsymbol{n}_{n+1} \right] \\ +% &= \underbrace{2 \mu \left[ \boldsymbol{\varepsilon}^{dev} - \boldsymbol{\varepsilon}^{p}_{n} \right]}_{\boldsymbol{\sigma}^{dev, trial}_{n+1}} +% - 2 \mu \Delta t \, \dot{\lambda}_{n+1} \boldsymbol{n}_{n+1} \\ + &= \underbrace{2 \mu \left[ \boldsymbol{\varepsilon}^{dev}_{n+1} - \boldsymbol{\varepsilon}^{p}_{n} \right]}_{\boldsymbol{\sigma}^{dev, trial}_{n+1}} + - 2 \mu \Delta \lambda_{n+1} \dfrac{\boldsymbol{\sigma}^{dev}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} +\quad (\textnormal{Note: } \Delta t \, \dot{\lambda}_{n+1} = \Delta \lambda_{n+1}) +%\\ +%\Rightarrow \quad +%\boldsymbol{\sigma}^{dev, trial}_{n+1} +% &= \left[ 1 + \dfrac{2 \mu \Delta \lambda_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} \right] \boldsymbol{\sigma}^{dev}_{n+1} +\\ +\Rightarrow \quad +\boldsymbol{\sigma}^{dev}_{n+1} + &= \left[ 1 + \dfrac{2 \mu \Delta \lambda_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} \right]^{-1} \boldsymbol{\sigma}^{dev, trial}_{n+1} +\end{align*} +\clearpage +### The final stress state is obtained by projecting the trial stress state onto the current yield surface. +\begin{figure}[!htb] +\centering +\includegraphics[width=0.7\linewidth]{Figures/Sec_3_5_1-Radial_projection} +\caption{Radial projection \citep{Mergheim2018a}} +\label{fig:sec351-radialprojection} +\end{figure} +## Still need to determine $\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert$. Use coaxiality of the stresses: +\begin{gather*} +\dfrac{\boldsymbol{\sigma}^{dev, trial}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} + = \dfrac{\boldsymbol{\sigma}^{dev}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} + = \boldsymbol{n}_{n+1} +\end{gather*} +### Contract both sides of previous equation by $\boldsymbol{n}_{n+1}$ and use definition of tensor norm $\Vert \boldsymbol{\sigma} \Vert$: +\begin{gather*} +\boldsymbol{\sigma}^{dev}_{n+1} : \dfrac{\boldsymbol{\sigma}^{dev}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} + = \boldsymbol{\sigma}^{dev, trial}_{n+1} : \dfrac{\boldsymbol{\sigma}^{dev, trial}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} + - 2 \mu \Delta \lambda_{n+1} \dfrac{\boldsymbol{\sigma}^{dev}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} : \dfrac{\boldsymbol{\sigma}^{dev}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert} \\ +\Rightarrow \quad +\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert + = \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert + - 2 \mu \Delta \lambda_{n+1} +\end{gather*} +## Yield function (plastic flow: $\Phi = 0$) +\begin{gather*} +\Phi_{n+1} + = \Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert - \sqrt{\frac{2}{3}} \left[ \sigma_{y} + \underbrace{k \alpha_{n+1}}_{- R_{n+1}} \right] + = 0 \\ +\Rightarrow \quad +\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert + = \sqrt{\frac{2}{3}} \left[ \sigma_{y} + k \left[ \alpha_{n} + \Delta \lambda_{n+1} \sqrt{\frac{2}{3}} \right] \right] +\end{gather*} +## Plastic multiplier update computed by equating the two above equations +\begin{gather*} +\Vert \boldsymbol{\sigma}^{dev}_{n+1} \Vert + = \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert + - 2 \mu \Delta \lambda_{n+1} + = \sqrt{\frac{2}{3}} \left[ \sigma_{y} + k \left[ \alpha_{n} + \Delta \lambda_{n+1} \sqrt{\frac{2}{3}} \right] \right] \\ +\Rightarrow \quad +\underbrace{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert - \sqrt{\frac{2}{3}} \left[ \sigma_{y} + k \alpha_{n} \right]}_{\Phi^{trial}_{n+1}} + = \Delta \lambda_{n+1} \left[ 2 \mu + \frac{2}{3} k \right] \\ +\Rightarrow \quad +\Delta \lambda_{n+1} + = \frac{\Phi^{trial}_{n+1}}{2 \mu + \frac{2}{3} k} +\end{gather*} +# Remarks +## Within ideal von Mises plasticity or with linear hardening, $\dot{\lambda}_{n+1}$ can be computed directly. + Otherwise, $\Phi_{n+1}$ has to be solver iteratively for $\dot{\lambda}_{n+1}$ (e.g. using Newton's method) +## For general yield functions, the closest point projection is the extension of the radial return algorithm. +%# Algorithm +%## Show slides +\end{easylist} + +\subsection{General closest-point projection (for general plasticity and hardening laws)} +\begin{easylist} +# Difference between radial and closest point projection +\begin{figure}[!htb] +\centering +\begin{subfigure}[t]{0.45\textwidth} +\includegraphics[width=0.8\linewidth]{Figures/Sec_3_5_1-Radial_projection} +\caption{Radial projection \citep{Mergheim2018a}} +\label{fig:sec351-radialprojection} +\end{subfigure} +\quad +\begin{subfigure}[t]{0.45\textwidth} +\includegraphics[width=0.8\linewidth]{Figures/Sec_3_5_2-General_closest_point_projection} +\caption{Closest-point projection \citep{Simo2006a}} +\label{fig:sec352-generalclosestpointprojection} +\end{subfigure} +\end{figure} +\end{easylist} + +\subsection{Consistent elastoplastic tangent modulus for von Mises plasticity with only linear isotropic hardening (time discrete setting)} +\begin{easylist} +# Note: The definition of the consistent elastoplastic tangent depends on the algorithms used to update stresses, internal variables +# Goal: Tangent to compute +\begin{gather} + \boldsymbol{\mathcal{C}}^{ep} + = \frac{\partial \boldsymbol{\sigma}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \boldsymbol{\mathcal{C}}^{ep, vol} + \boldsymbol{\mathcal{C}}^{ep, dev} +\label{equ: Constutitve law: Algorithmically correct tangent} +\end{gather} +# Stress tensors +\begin{gather*} +\boldsymbol{\sigma}_{n+1} + = \boldsymbol{\sigma}^{vol}_{n+1} + \boldsymbol{\sigma}^{dev}_{n+1} \\ +\boldsymbol{\sigma}^{vol}_{n+1} + = \kappa \boldsymbol{\varepsilon}^{vol}_{n+1} \\ +\boldsymbol{\sigma}^{dev}_{n+1} + = \boldsymbol{\sigma}^{dev, trial}_{n+1} + - 2 \mu \Delta \lambda_{n+1} \boldsymbol{n}_{n+1} \\ +\boldsymbol{\sigma}^{dev, trial}_{n+1} + = 2 \mu \left[ \boldsymbol{\varepsilon}^{dev}_{n+1} - \boldsymbol{\varepsilon}^{p}_{n} \right] +\end{gather*} +# Derivatives (for tangent stiffness contributions for trial solution) +\begin{gather*} +\boldsymbol{\mathcal{C}}^{ep, vol} + = \frac{\partial \boldsymbol{\sigma}^{vol}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \kappa \boldsymbol{1} \otimes \boldsymbol{1} \\ +\boldsymbol{\mathcal{C}}^{ep, dev, trial} + = \frac{\partial \boldsymbol{\sigma}^{dev, trial}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = 2 \mu \boldsymbol{\mathcal{I}}^{dev} +\end{gather*} +# Derivatives (for tangent stiffness contributions when plastic flow) +## Remember: +\begin{gather*} +\Delta \lambda_{n+1} + = \frac{\Phi^{trial}_{n+1}}{2 \mu + \frac{2}{3} k} + = \frac{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert - \sqrt{\frac{2}{3}} \left[ \sigma_{y} - R_{n} \right]}{2 \mu + \frac{2}{3} k} +\quad \textnormal{(Note: $R_{n}$ is fixed)} \\ +\boldsymbol{n}_{n+1} + = \dfrac{\boldsymbol{\sigma}^{dev, trial}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} +\end{gather*} +## Therefore: +\begin{gather*} +\frac{\partial \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \frac{\partial \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert}{\partial \boldsymbol{\sigma}^{dev, trial}_{n+1}} : \frac{\partial \boldsymbol{\sigma}^{dev, trial}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \underbrace{\frac{\boldsymbol{\sigma}^{dev, trial}_{n+1}}{\Vert\boldsymbol{\sigma}^{dev, trial}_{n+1}\Vert}}_{\boldsymbol{n}_{n+1}} : 2 \mu \boldsymbol{\mathcal{I}}^{dev} + = 2 \mu \boldsymbol{n}_{n+1} +\\ +\frac{\partial \Delta \lambda_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \frac{1}{{2 \mu + \frac{2}{3} k}} \frac{\partial \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \frac{2 \mu}{{2 \mu + \frac{2}{3} k}} \boldsymbol{n}_{n+1} \\ +\end{gather*} +\begin{align*} +\frac{\partial \boldsymbol{n}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} +% &= 2 \mu \boldsymbol{\mathcal{I}}^{dev}_{n+1} - \frac{2\mu}{2\mu + \frac{2}{3}k} \boldsymbol{n}_{n+1} \otimes \boldsymbol{n}_{n+1} \quad \textnormal{(from previous section)} \\ + &\equiv \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert^{-1} \frac{\partial \boldsymbol{\sigma}^{dev, trial}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + + \boldsymbol{\sigma}^{dev, trial}_{n+1} \otimes \frac{\partial \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert^{-1}}{\partial \boldsymbol{\varepsilon}_{n+1}} \\ + &= \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert^{-1} 2 \mu \boldsymbol{\mathcal{I}}^{dev} + - \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert^{-2} \boldsymbol{\sigma}^{dev, trial}_{n+1} \otimes \frac{\partial \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert}{\partial \boldsymbol{\varepsilon}_{n+1}} \\ + &= \Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert^{-1} \left[ 2 \mu \boldsymbol{\mathcal{I}}^{dev} + - \frac{\boldsymbol{\sigma}^{dev, trial}_{n+1}}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} \otimes 2 \mu \boldsymbol{n}_{n+1} \right] \\ + &= \frac{2\mu}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} \left[ \boldsymbol{\mathcal{I}}^{dev} + - \boldsymbol{n}_{n+1} \otimes \boldsymbol{n}_{n+1} \right] +\end{align*} +## Resulting deviatoric part of the elasto-plastic tangent +\begin{align} +\boldsymbol{\mathcal{C}}^{ep, dev} + &= \frac{\partial \boldsymbol{\sigma}^{dev}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + = \frac{\partial \boldsymbol{\sigma}^{dev, trial}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + - 2 \mu \left[ \boldsymbol{n}_{n+1} \otimes \frac{\partial \Delta \lambda_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} + + \Delta \lambda_{n+1} \frac{\partial \boldsymbol{n}_{n+1}}{\partial \boldsymbol{\varepsilon}_{n+1}} \right] +\\ + &= 2 \mu \boldsymbol{\mathcal{I}}^{dev} + - 2 \mu \left[ \boldsymbol{n}_{n+1} \otimes \left[ \frac{2 \mu}{{2 \mu + \frac{2}{3} k}} \boldsymbol{n}_{n+1} \right] + \right.\\ &\left. + + \Delta \lambda_{n+1} \left[ \frac{2\mu}{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert} \left[ \boldsymbol{\mathcal{I}}^{dev} - \boldsymbol{n}_{n+1} \otimes \boldsymbol{n}_{n+1} \right] \right] \right] +\\ + &= 2 \mu \boldsymbol{\mathcal{I}}^{dev} + - \frac{ \left[2 \mu \right]^{2}}{2 \mu + \frac{2}{3} k} \boldsymbol{n}_{n+1} \otimes \boldsymbol{n}_{n+1} \\ + & + \left[\Delta \lambda_{n+1} \frac{\left[ 2 \mu \right]^{2}}{{\Vert \boldsymbol{\sigma}^{dev, trial}_{n+1} \Vert}} \left[ \boldsymbol{\mathcal{I}}^{dev}_{n+1} - \boldsymbol{n}_{n+1} \otimes \boldsymbol{n}_{n+1} \right] \right] + \label{equ: Constutitve law: Algorithmically correct tangent: Deviatoric} +\end{align} +## Compared to continuous case, the tangents are not equal but they converge to the same expression when $\Delta \lambda_{n+1} \rightarrow 0$ +\end{easylist} + +\clearpage +\section{Finite element discretisation +\label{sec: FE discretisation} +} +For clarity, it is worth mentioning up front that in this \namecref{sec: FE discretisation} we will denote the general elasto-plastic tangent by $\boldsymbol{\mathbb{C}}^{ep}$; +so, in the elastic region $\boldsymbol{\mathbb{C}}^{ep} \equiv \boldsymbol{\mathbb{C}}$ while in the plastic region $\boldsymbol{\mathbb{C}}^{ep}$ represents the consistent tangent. + +Combining \cref{equ: Weak form of the balance of linear momentum,equ: Constutitve law: Algorithmically correct tangent} (specifically, \cref{equ: Constutitve law: Algorithmically correct tangent: Volumetric,equ: Constutitve law: Algorithmically correct tangent: Deviatoric}) renders the complete expression of the balance of linear momentum, namely +\begin{gather} +\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \left[ \boldsymbol{\mathbb{C}}^{ep} \right]_{ijkl} \varepsilon_{kl} \, dv + = \int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + + \int\limits_{\partial\Omega} \delta v_{i} \, \underbrace{\sigma_{ij} \, n_{j}}_{\bar{t}_{i}} \, da +\quad , +\label{equ: Weak form of the balance of linear momentum: Muscle model} +\end{gather} + +We discretise the trial solution and test function using vector-valued finite element shape functions (ansatz) +\begin{gather} +\mathbf{u} \left( \mathbf{x} \right) + \approx \sum\limits_{I} \mathbf{N}^{I} \left( \mathbf{x} \right) u^{I} +\quad , \quad +\mathbf{v} \left( \mathbf{x} \right) + \approx \sum\limits_{I} \mathbf{N}^{I} \left( \mathbf{x} \right) v^{I} +\end{gather} +where $\mathbf{N}^{I} \left( \mathbf{x} \right)$ is the (position-dependent) vector-valued finite element shape function corresponding to the $I^{\textnormal{th}}$ degree-of-freedom, and $u^{I}, v^{I}$ are coefficients of the solution and trial function. + +We now use these shape functions to discretise the weak expression for the balance of linear momentum. +Starting on the right-hand side of \cref{equ: Weak form of the balance of linear momentum: Muscle model}, the body force and traction contributions are computed by +\begin{align} +\int\limits_{\Omega} \delta v_{i} \, b_{i} \, dv + &= \int\limits_{\Omega} \left[ \sum\limits_{I} \mathbf{N}^{I} \left( \mathbf{x} \right) \delta v^{I} \right]_{i} \, b_{i} \, dv + = \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \mathbf{N}^{I}_{i} \, b_{i} \, dv +\label{equ: Discretisation: Body force} +\\ +\int\limits_{\Omega} \delta v_{i} \, t_{i} \, dv + &= \sum\limits_{I} \delta v^{I} \int\limits_{\Omega} \mathbf{N}^{I}_{i} \, t_{i} \, dv +\label{equ: Discretisation: Traction} +\quad . +\end{align} +The last component of \cref{equ: Weak form of the balance of linear momentum: Muscle model} that we wish to express in discrete form is the left-hand side of the equation. +Before we do, we observe that using the minor symmetry of the material stiffness tensor we can re-express the contraction of it and the small strain tensor as +\begin{gather} +\boldsymbol{\mathbb{C}} : \boldsymbol{\varepsilon} + = \boldsymbol{\mathbb{C}} : \frac{1}{2} \left[ \nabla \mathbf{u} + \left[ \nabla \mathbf{u} \right]^{T} \right] + \equiv \boldsymbol{\mathbb{C}} : \nabla \mathbf{u} +\end{gather} +from which we can similarly deduce that +\begin{gather} +\nabla \mathbf{u} : \boldsymbol{\mathbb{C}} + \equiv \boldsymbol{\varepsilon} : \boldsymbol{\mathbb{C}} +\end{gather} +Therefore, this contribution written in discrete form is +\begin{align} +&\int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \left[ \boldsymbol{\mathbb{C}}^{ep} \right]_{ijkl} \varepsilon_{kl} \, dv + \equiv \int\limits_{\Omega} \frac{\partial \delta v_{i}}{\partial x_{j}} \, \mathbb{C}^{ep}_{ijkl} \, \frac{\partial \delta u_{k}}{\partial x_{l}} \, dv \notag\\ + &= \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \sum\limits_{I} \mathbf{N}^{I} \left( \mathbf{x} \right) \delta v^{I} \right]_{i} \, \mathbb{C}^{ep}_{ijkl} \, \frac{\partial}{\partial x_{l}} \left[ \sum\limits_{J} \mathbf{N}^{J} \left( \mathbf{x} \right) \delta u^{J} \right]_{k} \, dv \notag\\ + &= \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \frac{\partial}{\partial x_{j}} \left[ \mathbf{N}^{I} \left( \mathbf{x} \right) \right]_{i} \, \mathbb{C}^{ep}_{ijkl} \, \frac{\partial}{\partial x_{l}} \left[ \mathbf{N}^{J} \left( \mathbf{x} \right) \right]_{k} \, dv \right] \delta u^{J} \notag\\ + &= \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \left[ \mathbf{N}^{I} \left( \mathbf{x} \right) \right]_{i,j} \, \mathbb{C}^{ep}_{ijkl} \, \left[ \mathbf{N}^{J} \left( \mathbf{x} \right) \right]_{k,l} \, dv \right] \delta u^{J} \notag\\ + &\equiv \sum\limits_{I,J} \delta v^{I} \left[ \int\limits_{\Omega} \left[\left[ \mathbf{N}^{I} \left( \mathbf{x} \right) \right]_{i,j} \right]^{S} \, \mathbb{C}^{ep}_{ijkl} \, \left[\left[ \mathbf{N}^{J} \left( \mathbf{x} \right) \right]_{k,l}\right]^{S} \, dv \right] \delta u^{J} +\label{equ: Discretisation: Material tangent} +\quad . +\end{align} + +\Cref{equ: Discretisation: Body force,equ: Discretisation: Traction,equ: Discretisation: Material tangent} are collectively used to develop the system of linear equations that are solved at each time step. + +\bibliographystyle{plain} +\bibliography{\jobname} + +\end{document} diff --git a/Linear_Elastoplasticity/linear_elastoplasticity.cc b/Linear_Elastoplasticity/linear_elastoplasticity.cc new file mode 100644 index 00000000..38efc322 --- /dev/null +++ b/Linear_Elastoplasticity/linear_elastoplasticity.cc @@ -0,0 +1,1613 @@ +/* --------------------------------------------------------------------- + * Copyright (C) 2018 by the deal.II authors and + * Jean-Paul Pelteret + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +/* + * Authors: Jean-Paul Pelteret + * University of Erlangen-Nuremberg, 2018 + */ + + +// @sect3{Include files} + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace LinearElastoplasticity +{ + using namespace dealii; + + +// @sect3{Run-time parameters} +// +// There are several parameters that can be set in the code so we set up a +// ParameterHandler object to read in the choices at run-time. + namespace Parameters + { + +// @sect4{Finite Element system} + +// Here we specify the polynomial order used to approximate the solution. +// The quadrature order should be adjusted accordingly. + struct FESystem + { + unsigned int poly_degree; + unsigned int quad_order; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void FESystem::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + prm.declare_entry("Polynomial degree", "1", + Patterns::Integer(0), + "Displacement system polynomial order"); + + prm.declare_entry("Quadrature order", "2", + Patterns::Integer(0), + "Gauss quadrature order"); + } + prm.leave_subsection(); + } + + void FESystem::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + poly_degree = prm.get_integer("Polynomial degree"); + quad_order = prm.get_integer("Quadrature order"); + } + prm.leave_subsection(); + } + + +// @sect4{Boundary conditions} + +// The boundary conditions to be applied to the problem. + + struct BoundaryConditions + { + double final_stretch; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void BoundaryConditions::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Boundary conditions"); + { + prm.declare_entry("Final stretch", "2.0", + Patterns::Double(0), + "The amount by which the specimen is to be stretched"); + } + prm.leave_subsection(); + } + + void BoundaryConditions::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Boundary conditions"); + { + final_stretch = prm.get_double("Final stretch"); + } + prm.leave_subsection(); + } + + +// @sect4{Geometry} + +// Make adjustments to the geometry of the notched specimen. + + struct Geometry + { + std::string geometry_type; + double scale; + unsigned int n_global_refinement_steps; + unsigned int n_local_refinement_steps; + + struct NotchedCylinder + { + double length; + double radius; + double notch_length; + double notch_radius; + } notched_cylinder; + + struct TensileSpecimen + { + const double length = 35.0/2; + + const Tensor<1,3> manifold_direction = Tensor<1,3>({0,0,1}); + const Point<3> centre_manifold_id_10 = {6.0,4.0,0.0}; + const Point<3> centre_manifold_id_11 = {10.472135955,0.0,0.0}; + } tensile_specimen; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Geometry::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + prm.declare_entry("Geometry type", "Notched cylinder", + Patterns::Selection("Notched cylinder|Tensile specimen"), + "The geometry to be modelled"); + + prm.declare_entry("Grid scale", "1e-3", + Patterns::Double(0), + "Global grid scaling factor"); + + prm.declare_entry("Global refinement steps", "1", + Patterns::Integer(0), + "Number of global refinement steps"); + + prm.declare_entry("Local refinement steps", "1", + Patterns::Integer(0), + "Number of initial local refinement cycles"); + + prm.enter_subsection("Notched cylinder"); + { + prm.declare_entry("Length", "100", + Patterns::Double(0), + "Overall length of the specimen"); + + prm.declare_entry("Radius", "10", + Patterns::Double(0), + "Overall radius of the specimen"); + + prm.declare_entry("Notch length", "10", + Patterns::Double(0), + "Overall length of the notch in the specimen"); + + prm.declare_entry("Notch radius", "1", + Patterns::Double(0), + "Overall radius of the notch in the specimen"); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + void Geometry::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + geometry_type = prm.get("Geometry type"); + scale = prm.get_double("Grid scale"); + n_global_refinement_steps = prm.get_integer("Global refinement steps"); + n_local_refinement_steps = prm.get_integer("Local refinement steps"); + + prm.enter_subsection("Notched cylinder"); + { + notched_cylinder.length = prm.get_double("Length"); + notched_cylinder.radius = prm.get_double("Radius"); + notched_cylinder.notch_length = prm.get_double("Notch length"); + notched_cylinder.notch_radius = prm.get_double("Notch radius"); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + AssertThrow(notched_cylinder.length > notched_cylinder.notch_length, + ExcMessage("Cannot create geometry with given dimensions")); + AssertThrow(notched_cylinder.radius > notched_cylinder.notch_radius, + ExcMessage("Cannot create geometry with given dimensions")); + } + + +// @sect4{Material properties} + +// Make adjustments to the material of the notched specimen. + + struct MaterialProperties + { + std::string mat_type; + double mu_e; // Shear modulus + double nu_e; // Poisson ratio + double k_p; // Isotropic hardening constant + double sigma_y_p; // Yield stress + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void MaterialProperties::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + prm.declare_entry("Type", "Elastic", + Patterns::Selection("Elastic|Elastoplastic (isotropic hardening)"), + "Type of material that composes the specimen"); + + prm.declare_entry("Shear modulus", "100", + Patterns::Double(0), + "Shear modulus of the specimen"); + + prm.declare_entry("Poisson ratio", "0.3", + Patterns::Double(-1,0.5), + "Poisson ratio of the specimen"); + + prm.declare_entry("Isotropic hardening constant", "10", + Patterns::Double(0), + "Isotropic hardening constant"); + + prm.declare_entry("Yield stress", "20", + Patterns::Double(0), + "Yield stress"); + } + prm.leave_subsection(); + } + + void MaterialProperties::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + mat_type = prm.get("Type"); + mu_e = prm.get_double("Shear modulus"); + nu_e = prm.get_double("Poisson ratio"); + k_p = prm.get_double("Isotropic hardening constant"); + sigma_y_p = prm.get_double("Yield stress"); + } + prm.leave_subsection(); + } + + +// @sect4{Time} + +// Set the timestep size $ \varDelta t $ and the simulation end-time. + struct Time + { + double delta_t; + double end_time; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "3", + Patterns::Double(0), + "End time"); + + prm.declare_entry("Time step size", "0.1", + Patterns::Double(0), + "Time step size"); + } + prm.leave_subsection(); + } + + void Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + } + prm.leave_subsection(); + } + + +// @sect4{All parameters} + +// Finally we consolidate all of the above structures into a single container +// that holds all of our run-time selections. + struct AllParameters : public FESystem, + public BoundaryConditions, + public Geometry, + public MaterialProperties, + public Time + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + } + + void AllParameters::declare_parameters(ParameterHandler &prm) + { + FESystem::declare_parameters(prm); + BoundaryConditions::declare_parameters(prm); + Geometry::declare_parameters(prm); + MaterialProperties::declare_parameters(prm); + Time::declare_parameters(prm); + } + + void AllParameters::parse_parameters(ParameterHandler &prm) + { + FESystem::parse_parameters(prm); + BoundaryConditions::parse_parameters(prm); + Geometry::parse_parameters(prm); + MaterialProperties::parse_parameters(prm); + Time::parse_parameters(prm); + } + } + + +// @sect3{Time incrementation} + + class Time + { + public: + Time (const double time_end, + const double delta_t) + : + timestep(0), + time_current(0.0), + time_end(time_end), + delta_t(delta_t) + {} + virtual ~Time() + {} + double current() const + { + return time_current; + } + double end() const + { + return time_end; + } + double get_delta_t() const + { + return delta_t; + } + unsigned int get_timestep() const + { + return timestep; + } + void increment() + { + time_current += delta_t; + ++timestep; + } + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; +}; + + +// @sect3{Constitutive laws} + + template + class Material_Base + { + public: + Material_Base() + {} + + virtual ~Material_Base() + {} + + // Cauchy stress + virtual SymmetricTensor<2,dim> + get_sigma(const SymmetricTensor<2, dim> &epsilon) const = 0; + + // Elastic / Elastoplastic tangent + virtual SymmetricTensor<4,dim> + get_K(const SymmetricTensor<2, dim> &epsilon) const = 0; + + virtual void + update_internal_equilibrium(const SymmetricTensor<2, dim> &epsilon) = 0; + + virtual void + update_end_timestep() = 0; + }; + + + /** + * Linear elastic material + */ + template + class Material_Linear_Elastic : public Material_Base + { + public: + Material_Linear_Elastic(const double mu_e, + const double nu_e) + : + kappa_e((2.0 * mu_e * (1.0 + nu_e)) / (3.0 * (1.0 - 2.0 * nu_e))), + mu_e(mu_e) + { + Assert(kappa_e > 0, ExcInternalError()); + } + + virtual ~Material_Linear_Elastic() + {} + + // Cauchy stress + SymmetricTensor<2,dim> + get_sigma(const SymmetricTensor<2, dim> &epsilon) const override + { + // Assume isochoric/deviatoric split + const SymmetricTensor<2, dim> epsilon_vol = (trace(epsilon)/dim)*Physics::Elasticity::StandardTensors::I; + const SymmetricTensor<2, dim> epsilon_dev = epsilon - epsilon_vol; + + return dim*kappa_e*epsilon_vol + + 2.0*mu_e*epsilon_dev; + } + + // Elastic / Elastoplastic tangent + SymmetricTensor<4,dim> get_K(const SymmetricTensor<2, dim> &epsilon) const override + { + (void)epsilon; + + // Elastic part of tangent, assuming an isochoric/deviatoric split + const SymmetricTensor<4,dim> K_e + = (kappa_e - 2.0/dim*mu_e)*Physics::Elasticity::StandardTensors::IxI + + (2.0*mu_e)*Physics::Elasticity::StandardTensors::S; + + // Check that stress-strain relationship really is linear and + // correctly implemented + Assert((get_sigma(epsilon) - K_e*epsilon).norm() < 1e-6, ExcInternalError()); + + return K_e; + } + + void + update_internal_equilibrium(const SymmetricTensor<2, dim> &epsilon) override + { + (void)epsilon; + } + + void + update_end_timestep() override + {} + + protected: + const double kappa_e; // bulk modulus + const double mu_e; // shear modulus + }; + + + /** + * Linear elastoplastic material: von Mises plasticity with linear isotropic hardening + */ + template + class Material_Linear_Elastoplastic_Isotropic_Hardening : public Material_Base + { + public: + Material_Linear_Elastoplastic_Isotropic_Hardening(const double mu_e, + const double nu_e, + const double k_p, + const double sigma_y_p, + const Time &time) + : + kappa_e((2.0 * mu_e * (1.0 + nu_e)) / (3.0 * (1.0 - 2.0 * nu_e))), + mu_e(mu_e), + k_p (k_p), + sigma_y_p (sigma_y_p), + time(time), + phi_t(std::numeric_limits::min()), + n_t(), + delta_lambda_t(0.0), + sigma_dev_trial_norm(1.0), + epsilon_p_t(), + epsilon_p_t1(), + alpha_p_t (0.0), + alpha_p_t1 (0.0) + { + Assert(kappa_e > 0, ExcInternalError()); + } + + virtual ~Material_Linear_Elastoplastic_Isotropic_Hardening() + {} + + // Cauchy stress + SymmetricTensor<2,dim> + get_sigma(const SymmetricTensor<2, dim> &epsilon) const override + { + return get_sigma(epsilon, epsilon_p_t); + } + + // Isotropic hardening stress + double + get_R() const + { + return get_R(alpha_p_t); + } + + // Elastic / Elastoplastic tangent + SymmetricTensor<4,dim> + get_K(const SymmetricTensor<2, dim> &epsilon) const override + { + (void)epsilon; + + // Elastic part of tangent, assuming an isochoric/deviatoric split + const SymmetricTensor<4,dim> K_e + = (kappa_e - 2.0/dim*mu_e)*Physics::Elasticity::StandardTensors::IxI + + (2.0*mu_e)*Physics::Elasticity::StandardTensors::S; + + if (phi_t <= 0) + { + // Check that stress-strain relationship really is linear and + // correctly implemented + Assert((get_sigma(epsilon) - K_e*epsilon).norm() < 1e-6, ExcInternalError()); + + return K_e; + } + else + { + // Compute correction for radial projection during plastic flow + static const SymmetricTensor<4,dim> I_dev + = Physics::Elasticity::StandardTensors::S + - (1.0/dim)*Physics::Elasticity::StandardTensors::IxI; + const SymmetricTensor<4,dim> n_t_x_n_t = outer_product(n_t,n_t); + + Assert(sigma_dev_trial_norm != 0.0, ExcInternalError()); + const SymmetricTensor<4,dim> K_p_correction + = - (Utilities::fixed_power<2>(2.0*mu_e))/(2.0*mu_e + 2.0/dim*kappa_e)*n_t_x_n_t + - (Utilities::fixed_power<2>(2.0*mu_e)*delta_lambda_t)/(sigma_dev_trial_norm)*(I_dev - n_t_x_n_t); + + return K_e + K_p_correction; + } + } + + void + update_internal_equilibrium(const SymmetricTensor<2, dim> &epsilon) override + { + // Trial stress + const SymmetricTensor<2, dim> sigma_trial = get_sigma(epsilon,epsilon_p_t1); + // Deviatoric part of trial stress + const SymmetricTensor<2, dim> sigma_dev_trial = deviator(sigma_trial); + // Trial hardening stress + const double R_trial = get_R(alpha_p_t1); + + // Update yield function value (von Mises plasticity for linear + // isotropic hardening) + { + phi_t = sigma_dev_trial.norm() - std::sqrt(2.0/dim)*(sigma_y_p - R_trial); + + // Update internal variables based condition of + // elasticity or plastic flow + if (phi_t <= 0.0) + { + // Elastic update + epsilon_p_t = epsilon_p_t1; + alpha_p_t = alpha_p_t1; + } + else + { + // Plastic update: + // Linear isotropic hardening -> use radial return algorithm + + // Internal variables for radial return algorithm + sigma_dev_trial_norm = sigma_dev_trial.norm(); + n_t = sigma_dev_trial / sigma_dev_trial_norm; + delta_lambda_t = phi_t / (2.0*mu_e + 2.0/dim*kappa_e); + + // Update internal plastic variables + epsilon_p_t = epsilon_p_t1 + delta_lambda_t*n_t; + alpha_p_t = alpha_p_t1 + delta_lambda_t*std::sqrt(2.0/dim); + } + } + } + + void + update_end_timestep() override + { + // Record history of plastic variables + epsilon_p_t1 = epsilon_p_t; + alpha_p_t1 = alpha_p_t; + } + + // --- Post-processing --- + + // Plastic variables: + // Isotropic hardening variable + double + get_alpha_p() const + { + return alpha_p_t; + } + + protected: + const double kappa_e; // buld modulus + const double mu_e; // shear modulus + const double k_p; // isotropic hardening modulus + const double sigma_y_p; // yield stress + + const Time &time; + + // Variables related to yield function at current timestep + double phi_t; // yield function + SymmetricTensor<2,dim> n_t; // radial return direction + double delta_lambda_t;// increment in plastic multiplier + double sigma_dev_trial_norm; // norm of deviatoric part of trial stress + + // Internal variables: + // Plastic strain + SymmetricTensor<2,dim> epsilon_p_t; // current timestep + SymmetricTensor<2,dim> epsilon_p_t1; // previous timestep + // Isotropic hardening variable + double alpha_p_t; // current timestep + double alpha_p_t1; // previous timestep + + SymmetricTensor<2,dim> + get_sigma(const SymmetricTensor<2, dim> &epsilon, + const SymmetricTensor<2, dim> &epsilon_p) const + { + // Assume isochoric/deviatoric split + const SymmetricTensor<2, dim> epsilon_vol = (trace(epsilon)/dim)*Physics::Elasticity::StandardTensors::I; + const SymmetricTensor<2, dim> epsilon_dev = epsilon - epsilon_vol; + + const SymmetricTensor<2,dim> sigma + = dim*kappa_e*epsilon_vol + + 2.0*mu_e*(epsilon_dev - epsilon_p); + return sigma; + } + + double + get_R(const double alpha_p) const + { + return -k_p*alpha_p; + } + }; + + +// @sect3{Quadrature point data} + + template + class PointHistory + { + public: + PointHistory() + {} + virtual ~PointHistory() + {} + void + setup_lqp (const types::material_id mat_id, + const Parameters::AllParameters ¶meters, + const Time &time) + { + if (mat_id == 0) + { + material.reset(new Material_Linear_Elastic( + parameters.mu_e, parameters.nu_e)); + } + else if (mat_id == 1) + { + material.reset(new Material_Linear_Elastoplastic_Isotropic_Hardening( + parameters.mu_e, parameters.nu_e, + parameters.k_p, parameters.sigma_y_p, + time)); + } + else + AssertThrow(mat_id <= 1, ExcMessage("Unknown material id")); + } + + SymmetricTensor<2, dim> + get_sigma(const SymmetricTensor<2, dim> &epsilon) const + { + return material->get_sigma(epsilon); + } + SymmetricTensor<4, dim> + get_K(const SymmetricTensor<2, dim> &epsilon) const + { + return material->get_K(epsilon); + } + void + update_internal_equilibrium(const SymmetricTensor<2, dim> &epsilon) + { + material->update_internal_equilibrium(epsilon); + } + void + update_end_timestep() + { + material->update_end_timestep(); + } + const Material_Base * const + get_material() const + { + Assert(material, ExcInternalError()); + return &(*material); + } + private: + std::shared_ptr< Material_Base > material; +}; + + +// @sect3{The LinearElastoplasticProblem class template} + + template + class LinearElastoplasticProblem + { + public: + LinearElastoplasticProblem (const std::string &input_file); + ~LinearElastoplasticProblem (); + void run (); + + private: + void make_grid (); + void refine_grid (); + void setup_system (); + void setup_qph(); + void make_constraints (); + void assemble_system (); + void solve (); + void output_results () const; + void update_end_timestep (); + + Parameters::AllParameters parameters; + + Triangulation triangulation; + DoFHandler dof_handler; + + Time time; + + FESystem fe; + QGauss qf_cell; + QGauss qf_face; + + AffineConstraints hanging_node_constraints; + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + // Local data + CellDataStorage::active_cell_iterator, + PointHistory > quadrature_point_history; + + const FEValuesExtractors::Vector u_fe; + static const unsigned int n_components = dim; + static const unsigned int first_u_component = 0; + }; + + +// @sect4{LinearElastoplasticProblem::LinearElastoplasticProblem} + + template + LinearElastoplasticProblem::LinearElastoplasticProblem (const std::string &input_file) + : + parameters(input_file), + dof_handler (triangulation), + time (parameters.end_time, parameters.delta_t), + fe (FE_Q(parameters.poly_degree), dim), + qf_cell (parameters.quad_order), + qf_face (parameters.quad_order), + u_fe(first_u_component) + { + Assert(dim==3, ExcNotImplemented()); + } + + +// @sect4{LinearElastoplasticProblem::~LinearElastoplasticProblem} + + template + LinearElastoplasticProblem::~LinearElastoplasticProblem () + { + dof_handler.clear (); + } + + +// @sect4{LinearElastoplasticProblem::make_grid} + + template + void LinearElastoplasticProblem::make_grid () + { + Assert (dim == 3, ExcNotImplemented()); + + if (parameters.geometry_type == "Notched cylinder") + { + const double notch_radius_inner = parameters.notched_cylinder.notch_radius - (parameters.notched_cylinder.radius-parameters.notched_cylinder.notch_radius); + + Triangulation<2> tria_cross_section; + { + const Point<2> centre; + Triangulation<2> tria_inner_1; + GridGenerator::hyper_ball (tria_inner_1,centre, + notch_radius_inner,false); + + Triangulation<2> tria_inner_2; + GridGenerator::hyper_shell (tria_inner_2,centre, + notch_radius_inner, + parameters.notched_cylinder.notch_radius, + 4); + GridTools::rotate(M_PI/4,tria_inner_2); + + Triangulation<2> tria_outer; + GridGenerator::hyper_shell (tria_outer,centre, + parameters.notched_cylinder.notch_radius, + parameters.notched_cylinder.radius, + 4); + GridTools::rotate(M_PI/4,tria_outer); + + Triangulation<2> tria_inner; + GridGenerator::merge_triangulations(tria_inner_1,tria_inner_2,tria_inner); + GridGenerator::merge_triangulations(tria_inner,tria_outer,tria_cross_section); + } + + Triangulation tria_unnotched; + { + std::vector slice_coordinates; + slice_coordinates.push_back(0.0); + slice_coordinates.push_back(parameters.notched_cylinder.notch_length); + slice_coordinates.push_back(2.0*parameters.notched_cylinder.notch_length); + // Biased slices with function taken from GridGenerator::concentric_hyper_shells + const double bias_first_slice = 2.0*parameters.notched_cylinder.notch_length; + const double bias_last_slice = parameters.notched_cylinder.length; + const double skewness = 2.0; + const unsigned int n_slices = 5; + for (unsigned int s=0; s(n_slices)))) / std::tanh(skewness); + AssertThrow(f >= 0.0 && f <= 1.0, ExcMessage("Bias function not in correct range.")); + const double r = bias_first_slice + (bias_last_slice-bias_first_slice)*f; + AssertThrow(r >= bias_first_slice && r <= bias_last_slice, ExcMessage("New slice location not in correct range.")); + slice_coordinates.push_back(r); + } + slice_coordinates.push_back(parameters.notched_cylinder.length); + GridGenerator::extrude_triangulation (tria_cross_section, + slice_coordinates, + tria_unnotched); + GridTools::rotate(M_PI/2,1,tria_unnotched); + } + + std::set::active_cell_iterator> cells_to_remove; + for (typename Triangulation::active_cell_iterator + cell = tria_unnotched.begin_active(); + cell!=tria_unnotched.end(); ++cell) + { + // Remove all cells within the notch length and greater than + // the notch radius. + if (cell->center()[0] < parameters.notched_cylinder.notch_length) + { + // Outer layer of cells correspond to the notch to be removed + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() == true && + std::abs(cell->face(face)->center()[0] - parameters.notched_cylinder.notch_length/2.0) < 1e-12) + cells_to_remove.insert(cell); + + // Additional verification + Point rad = cell->center(); + rad[0] = 0.0; + const double radius = rad.norm(); + if (radius > parameters.notched_cylinder.notch_radius) + cells_to_remove.insert(cell); + } + } + + AssertThrow(cells_to_remove.empty() == false, ExcMessage("Found no cells to remove.")); + GridGenerator::create_triangulation_with_removed_cells(tria_unnotched,cells_to_remove,triangulation); + + // Set boundary and manifold IDs + triangulation.set_all_manifold_ids(2); // TFI manifold + + const double radius_for_cyl_cells = notch_radius_inner*std::cos(M_PI/4); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + { + Point rad = cell->center(); + rad[0] = 0.0; + const double radius = rad.norm(); + if (radius > radius_for_cyl_cells) + { + cell->set_all_manifold_ids(1); // Cylindrical manifold + } + } + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + { + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true) + { + if (std::abs(cell->face(face)->center()[0] - 0.0) < 1e-6) + cell->face(face)->set_boundary_id(1); + else if (std::abs(cell->face(face)->center()[0] - parameters.notched_cylinder.length) < 1e-6) + cell->face(face)->set_boundary_id(2); + else + cell->face(face)->set_all_manifold_ids(1); // Cylindrical manifold + } + } + } + + CylindricalManifold cylindrical_manifold (0 /*axis*/); + CylindricalManifold cylindrical_manifold_tfi (0 /*axis*/); + TransfiniteInterpolationManifold tfi_manifold; + + triangulation.set_manifold (1, cylindrical_manifold); + tfi_manifold.initialize(triangulation); + triangulation.set_manifold (2, tfi_manifold); + + triangulation.refine_global(parameters.n_global_refinement_steps); + GridTools::scale (parameters.scale, triangulation); + } + else if (parameters.geometry_type == "Tensile specimen") + { + const std::string filename = "tensile_specimen.inp"; + std::ifstream input (filename.c_str()); + + GridIn gridin; + gridin.attach_triangulation(triangulation); + gridin.read_abaqus(input, false); + +// GridTools::copy_boundary_to_manifold_id(triangulation,false); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + { + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() == true) + { + if (cell->face(face)->boundary_id() == 10 || cell->face(face)->boundary_id() == 11) + cell->face(face)->set_all_manifold_ids(cell->face(face)->boundary_id()); + } + } + } + + // The triangulation scaling needs to happen before the manifolds + // with a non-trivial centre point are defined / applied. + GridTools::scale (parameters.scale, triangulation); + + CylindricalManifold cylindrical_manifold_id_10 (parameters.tensile_specimen.manifold_direction,parameters.tensile_specimen.centre_manifold_id_10*parameters.scale); + CylindricalManifold cylindrical_manifold_id_11 (parameters.tensile_specimen.manifold_direction,parameters.tensile_specimen.centre_manifold_id_11*parameters.scale); + triangulation.set_manifold (10, cylindrical_manifold_id_10); + triangulation.set_manifold (11, cylindrical_manifold_id_11); + + triangulation.refine_global(parameters.n_global_refinement_steps); + } + else + AssertThrow(false, ExcMessage("Unknown geometry")); + + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell!=triangulation.end(); ++cell) + { + if (parameters.mat_type == "Elastic") + cell->set_material_id(0); + else if (parameters.mat_type == "Elastoplastic (isotropic hardening)") + cell->set_material_id(1); + else + AssertThrow(false, ExcMessage("Unknown material id")); + } + } + + +// @sect4{LinearElastoplasticProblem::refine_grid} + + template + void LinearElastoplasticProblem::refine_grid () + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.03); + triangulation.execute_coarsening_and_refinement(); + } + + +// @sect4{LinearElastoplasticProblem::setup_system} + + template + void LinearElastoplasticProblem::setup_system () + { + dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + hanging_node_constraints, + /*keep_constrained_dofs = */ true); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit (sparsity_pattern); + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + // Setup quadrature point history + setup_qph(); + } + + +// @sect4{LinearElastoplasticProblem::setup_qph} + + template + void LinearElastoplasticProblem::setup_qph() + { + std::cout << "Setting up quadrature point data..." << std::endl; + const unsigned int n_q_points_cell = qf_cell.size(); + + quadrature_point_history.clear(); + quadrature_point_history.initialize(triangulation.begin_active(), + triangulation.end(), + n_q_points_cell); + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + { + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points_cell, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points_cell; ++q_point) + lqph[q_point]->setup_lqp(cell->material_id(), parameters, time); + } +} + + // @sect4{LinearElastoplasticProblem::assemble_system} + + template + void LinearElastoplasticProblem::assemble_system () + { + system_matrix = 0; + system_rhs = 0; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<1, 1, 1>; + using CellIteratorType = decltype(dof_handler.begin_active()); + + const ScratchData sample_scratch_data (fe, qf_cell, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + CopyData sample_copy_data (fe.dofs_per_cell); + + auto cell_worker = [this] (const CellIteratorType &cell, + ScratchData &scratch_data, + CopyData ©_data) + { + const FEValues &fe_values = scratch_data.reinit(cell); + FullMatrix &cell_matrix = copy_data.matrices[0]; + + std::vector &local_dof_indices = copy_data.local_dof_indices[0]; + cell->get_dof_indices(local_dof_indices); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell; + const unsigned int n_q_points_cell = fe_values.n_quadrature_points; + + std::vector > solution_grads_u_total (n_q_points_cell); + fe_values[u_fe].get_function_gradients( + solution, + solution_grads_u_total); + + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + + for (unsigned int q_point_cell=0; q_point_cell &Grad_u = solution_grads_u_total[q_point_cell]; + const SymmetricTensor<2,dim> epsilon = Physics::Elasticity::Kinematics::epsilon(Grad_u); + lqph[q_point_cell]->update_internal_equilibrium(epsilon); + const SymmetricTensor<4,dim> K = lqph[q_point_cell]->get_K(epsilon); + + const double JxW = fe_values.JxW(q_point_cell); + + // Precompute the shape function symmetric gradients + // (related to the variation of the small strain) + std::vector< SymmetricTensor<2,dim> > d_epsilon (dofs_per_cell); + for (unsigned int K=0; K &d_epsilon_I = d_epsilon[I]; + + for (unsigned int J=I; J &d_epsilon_J = d_epsilon[J]; + + cell_matrix(I,J) + += contract3(d_epsilon_I,K,d_epsilon_J) * JxW; + } + } + } + + for (unsigned int I=0; I &cell_matrix = copy_data.matrices[0]; + const Vector &cell_rhs = copy_data.vectors[0]; + const std::vector &local_dof_indices = copy_data.local_dof_indices[0]; + + constraints.distribute_local_to_global(cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + }; + + MeshWorker::mesh_loop(dof_handler.active_cell_iterators(), + cell_worker, copier, + sample_scratch_data, + sample_copy_data, + MeshWorker::assemble_own_cells); + } + + +// @sect4{LinearElastoplasticProblem::make_constraints} + + template + void LinearElastoplasticProblem::make_constraints () + { + constraints.clear(); + + if (parameters.geometry_type == "Notched cylinder") + { + // Fully constraints (pinned) vertex at centre of notched (-X) face + { + Point pinned_point; + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (cell->vertex(v).distance(pinned_point) < 1e-9) + { + for (unsigned int d=0; dvertex_dof_index(v,d)); + } + } + } + // No X displacement on constraint on -X faces + { + ComponentMask component_mask_x (n_components, false); + component_mask_x.set(0, true); + VectorTools::interpolate_boundary_values (dof_handler, + 1, + ZeroFunction(dim), + constraints, + component_mask_x); + } + // Prescribed horizontal displacement on +X faces + { + ComponentMask component_mask_x (n_components, false); + component_mask_x.set(0, true); + const double total_displacement + = (parameters.final_stretch-1.0)*parameters.notched_cylinder.length + * parameters.scale + * (time.current()/time.end()); + VectorTools::interpolate_boundary_values (dof_handler, + 2, + ConstantFunction(total_displacement,dim), + constraints, + component_mask_x); + } + } + else if (parameters.geometry_type == "Tensile specimen") + { + // No X displacement on constraint on -X faces + { + ComponentMask component_mask_x (n_components, false); + component_mask_x.set(0, true); + VectorTools::interpolate_boundary_values (dof_handler, + 1, + ZeroFunction(dim), + constraints, + component_mask_x); + } + // No Y displacement on constraint on -Y faces + { + ComponentMask component_mask_y (n_components, false); + component_mask_y.set(1, true); + VectorTools::interpolate_boundary_values (dof_handler, + 3, + ZeroFunction(dim), + constraints, + component_mask_y); + } + // No Z displacement on constraint on -Z faces + { + ComponentMask component_mask_z (n_components, false); + component_mask_z.set(2, true); + VectorTools::interpolate_boundary_values (dof_handler, + 5, + ZeroFunction(dim), + constraints, + component_mask_z); + } + // Prescribed horizontal displacement on clamped surface + { + ComponentMask component_mask_x (n_components, false); + component_mask_x.set(0, true); + const double total_displacement + = (parameters.final_stretch-1.0)*parameters.tensile_specimen.length + * parameters.scale + * (time.current()/time.end()); + VectorTools::interpolate_boundary_values (dof_handler, + 6, + ConstantFunction(total_displacement,dim), + constraints, + component_mask_x); + } + // No out-of-plane displacement on clamped surface + { + ComponentMask component_mask_yz (n_components, true); + component_mask_yz.set(0, false); + VectorTools::interpolate_boundary_values (dof_handler, + 6, + ZeroFunction(dim), + constraints, + component_mask_yz); + } + } + else + AssertThrow(false, ExcMessage("Unknown geometry")); + + constraints.merge(hanging_node_constraints, AffineConstraints::right_object_wins); + constraints.close(); + } + + +// @sect4{LinearElastoplasticProblem::solve} + + template + void LinearElastoplasticProblem::solve () + { + SolverControl solver_control (system_matrix.m(), 1e-12); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + constraints.distribute(solution); + } + + +// @sect4{LinearElastoplasticProblem::update_end_timestep} + + template + void LinearElastoplasticProblem::update_end_timestep () + { + const unsigned int n_q_points_cell = qf_cell.size(); + + for (typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell!=dof_handler.end(); ++cell) + { + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + + for (unsigned int q_point_cell=0; q_point_cellupdate_end_timestep(); + } + } + + + // @sect4{LinearElastoplasticProblem::output_results} + + template + class PostProcessIsotropicHardening : public DataPostprocessorScalar + { + public: + PostProcessIsotropicHardening(const CellDataStorage::active_cell_iterator, + PointHistory > &quadrature_point_history, + const QGauss &qf_cell); + virtual ~PostProcessIsotropicHardening() = default; + virtual void evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> &computed_quantities) const override; + + private: + const CellDataStorage::active_cell_iterator, + PointHistory > &quadrature_point_history; + + FullMatrix projection_matrix_qpoint_to_support_point; + }; + + + template + PostProcessIsotropicHardening::PostProcessIsotropicHardening( + const CellDataStorage::active_cell_iterator, + PointHistory > &quadrature_point_history, + const QGauss &quadrature_cell) + : DataPostprocessorScalar("isotropic_hardening", update_values), + quadrature_point_history (quadrature_point_history) + { + const FE_Q fe_projection (1); + const QTrapez quadrature_projection; + projection_matrix_qpoint_to_support_point = FullMatrix ( + quadrature_projection.size(), + quadrature_cell.size()); + + FETools::compute_projection_from_quadrature_points_matrix + (fe_projection, + quadrature_projection, + quadrature_cell, + projection_matrix_qpoint_to_support_point); + } + + + template + void PostProcessIsotropicHardening::evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> & computed_quantities) const + { + Assert(computed_quantities.size() == inputs.solution_values.size(), + ExcDimensionMismatch(computed_quantities.size(), + inputs.solution_values.size())); + + const typename DoFHandler::active_cell_iterator cell + = inputs.template get_cell>(); + + // number of support points (nodes) to project to + const unsigned int n_support_points = projection_matrix_qpoint_to_support_point.n_rows(); + // number of quadrature points to project from + const unsigned int n_quad_points = projection_matrix_qpoint_to_support_point.n_cols(); + + // component projected to the nodes + Vector component_at_node(n_support_points); + // component at the quadrature point + Vector component_at_qp(n_quad_points); + + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_quad_points, ExcDimensionMismatch(lqph.size(), n_quad_points)); + + // populate the vector of components at the qps + for (unsigned int q_point = 0; q_point < n_quad_points; ++q_point) + { + const Material_Base * const material = lqph[q_point]->get_material(); + if (const Material_Linear_Elastoplastic_Isotropic_Hardening* const material_ep + = dynamic_cast* const>(material)) + component_at_qp(q_point) = material_ep->get_alpha_p(); + else + component_at_qp(q_point) = 0.0; + } + + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix_qpoint_to_support_point.vmult(component_at_node, component_at_qp); + + Assert(computed_quantities.size() == n_support_points, + ExcDimensionMismatch(computed_quantities.size(), + n_support_points)); + for (unsigned int i = 0; i < n_support_points; i++) + { + Assert(inputs.solution_values[i].size() == dim, + ExcDimensionMismatch(inputs.solution_values[i].size(), dim)); + (void)inputs; + + computed_quantities[i](0) = component_at_node(i); + } + } + + + template + void LinearElastoplasticProblem::output_results () const + { + std::string filename = "solution-"; + filename += Utilities::int_to_string(time.get_timestep(),4); + filename += ".vtu"; + std::ofstream output (filename.c_str()); + + PostProcessIsotropicHardening pp_isotropic_hardening (quadrature_point_history, qf_cell); + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + std::vector + data_component_interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + std::vector solution_name(dim, "displacement"); + + data_out.add_data_vector (solution, solution_name, + DataOut::type_dof_data, + data_component_interpretation); + data_out.add_data_vector(solution, pp_isotropic_hardening); + + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags (flags); + data_out.build_patches (StaticMappingQ1::mapping, fe.degree, DataOut::curved_inner_cells); + data_out.write_vtu (output); + + static std::vector< std::pair< double, std::string >> times_and_names; + times_and_names.emplace_back (time.current(), filename); + std::ofstream pvd_output ("solution.pvd"); + DataOutBase::write_pvd_record (pvd_output, times_and_names); + } + + + // @sect4{LinearElastoplasticProblem::run} + + template + void LinearElastoplasticProblem::run () + { + make_grid(); + setup_system (); + output_results (); + + update_end_timestep(); + time.increment(); + + auto solve_timestep = [&]() + { + // Next we assemble the system and enforce boundary + // conditions. + make_constraints(); + assemble_system (); + + // Then we solve the linear system + solve (); + }; + + while (time.current() < time.end()+0.01*time.get_delta_t()) + { + std::cout + << "Timestep " << time.get_timestep() + << " @ time " << time.current() + << std::endl; + + // Compute the solution at the current timestep + solve_timestep(); + + // Perform local refinement + if (time.get_timestep() == 1) + for (unsigned int cycle = 0; cycle < parameters.n_local_refinement_steps; ++cycle) + { + std::cout + << "Executing refinement cycle " << cycle + << " of " << parameters.n_local_refinement_steps + << "..." << std::endl; + refine_grid(); + setup_system(); + solve_timestep(); + } + + // Output some values to file + output_results (); + + update_end_timestep(); + time.increment(); + } + } +} // namespace LinearElastoplasticity + + +// @sect3{The main function} + +int main (int argc, char *argv[]) +{ + try + { + dealii::deallog.depth_console (0); + const unsigned int dim = 3; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + + LinearElastoplasticity::LinearElastoplasticProblem elastoplastic_problem ("parameters.prm"); + elastoplastic_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/Linear_Elastoplasticity/mesh/tensile_specimen.inp b/Linear_Elastoplasticity/mesh/tensile_specimen.inp new file mode 100755 index 00000000..8f037c58 --- /dev/null +++ b/Linear_Elastoplasticity/mesh/tensile_specimen.inp @@ -0,0 +1,422 @@ +*HEADING +cubit(ems/elasticity/tensile_specimen/mesh/tensile_specimen.inp): 10/23/2012: 16 +version: 13.2 +** +********************************** P A R T S ********************************** +*PART, NAME=Part-Default +** +********************************** N O D E S ********************************** +*NODE, NSET=ALLNODES + 1, 0.000000e+00, 1.000000e+00, 0.000000e+00 + 2, 0.000000e+00, 5.000000e-01, 0.000000e+00 + 3, 5.000000e-01, 5.000000e-01, 0.000000e+00 + 4, 5.000000e-01, 1.000000e+00, 0.000000e+00 + 5, 0.000000e+00, 1.000000e+00, 2.500000e-01 + 6, 0.000000e+00, 5.000000e-01, 2.500000e-01 + 7, 5.000000e-01, 5.000000e-01, 2.500000e-01 + 8, 5.000000e-01, 1.000000e+00, 2.500000e-01 + 9, 0.000000e+00, 0.000000e+00, 0.000000e+00 + 10, 5.000000e-01, 0.000000e+00, 0.000000e+00 + 11, 0.000000e+00, 0.000000e+00, 2.500000e-01 + 12, 5.000000e-01, 0.000000e+00, 2.500000e-01 + 13, 1.000000e+00, 5.000000e-01, 0.000000e+00 + 14, 1.000000e+00, 1.000000e+00, 0.000000e+00 + 15, 1.000000e+00, 5.000000e-01, 2.500000e-01 + 16, 1.000000e+00, 1.000000e+00, 2.500000e-01 + 17, 1.000000e+00, 0.000000e+00, 0.000000e+00 + 18, 1.000000e+00, 0.000000e+00, 2.500000e-01 + 19, 1.500000e+00, 5.000000e-01, 0.000000e+00 + 20, 1.500000e+00, 1.000000e+00, 0.000000e+00 + 21, 1.500000e+00, 5.000000e-01, 2.500000e-01 + 22, 1.500000e+00, 1.000000e+00, 2.500000e-01 + 23, 1.500000e+00, 0.000000e+00, 0.000000e+00 + 24, 1.500000e+00, 0.000000e+00, 2.500000e-01 + 25, 2.000000e+00, 5.000000e-01, 0.000000e+00 + 26, 2.000000e+00, 1.000000e+00, 0.000000e+00 + 27, 2.000000e+00, 5.000000e-01, 2.500000e-01 + 28, 2.000000e+00, 1.000000e+00, 2.500000e-01 + 29, 2.000000e+00, 0.000000e+00, 0.000000e+00 + 30, 2.000000e+00, 0.000000e+00, 2.500000e-01 + 31, 2.500000e+00, 5.000000e-01, 0.000000e+00 + 32, 2.500000e+00, 1.000000e+00, 0.000000e+00 + 33, 2.500000e+00, 5.000000e-01, 2.500000e-01 + 34, 2.500000e+00, 1.000000e+00, 2.500000e-01 + 35, 2.500000e+00, 0.000000e+00, 0.000000e+00 + 36, 2.500000e+00, 0.000000e+00, 2.500000e-01 + 37, 3.000000e+00, 5.000000e-01, 0.000000e+00 + 38, 3.000000e+00, 1.000000e+00, 0.000000e+00 + 39, 3.000000e+00, 5.000000e-01, 2.500000e-01 + 40, 3.000000e+00, 1.000000e+00, 2.500000e-01 + 41, 3.000000e+00, 0.000000e+00, 0.000000e+00 + 42, 3.000000e+00, 0.000000e+00, 2.500000e-01 + 43, 3.500000e+00, 5.000000e-01, 0.000000e+00 + 44, 3.500000e+00, 1.000000e+00, 0.000000e+00 + 45, 3.500000e+00, 5.000000e-01, 2.500000e-01 + 46, 3.500000e+00, 1.000000e+00, 2.500000e-01 + 47, 3.500000e+00, 0.000000e+00, 0.000000e+00 + 48, 3.500000e+00, 0.000000e+00, 2.500000e-01 + 49, 4.000000e+00, 5.000000e-01, 0.000000e+00 + 50, 4.000000e+00, 1.000000e+00, 0.000000e+00 + 51, 4.000000e+00, 5.000000e-01, 2.500000e-01 + 52, 4.000000e+00, 1.000000e+00, 2.500000e-01 + 53, 4.000000e+00, 0.000000e+00, 0.000000e+00 + 54, 4.000000e+00, 0.000000e+00, 2.500000e-01 + 55, 4.500000e+00, 5.000000e-01, 0.000000e+00 + 56, 4.500000e+00, 1.000000e+00, 0.000000e+00 + 57, 4.500000e+00, 5.000000e-01, 2.500000e-01 + 58, 4.500000e+00, 1.000000e+00, 2.500000e-01 + 59, 4.500000e+00, 0.000000e+00, 0.000000e+00 + 60, 4.500000e+00, 0.000000e+00, 2.500000e-01 + 61, 5.000000e+00, 5.000000e-01, 0.000000e+00 + 62, 5.000000e+00, 1.000000e+00, 0.000000e+00 + 63, 5.000000e+00, 5.000000e-01, 2.500000e-01 + 64, 5.000000e+00, 1.000000e+00, 2.500000e-01 + 65, 5.000000e+00, 0.000000e+00, 0.000000e+00 + 66, 5.000000e+00, 0.000000e+00, 2.500000e-01 + 67, 0.000000e+00, 1.000000e+00, 5.000000e-01 + 68, 0.000000e+00, 5.000000e-01, 5.000000e-01 + 69, 5.000000e-01, 5.000000e-01, 5.000000e-01 + 70, 5.000000e-01, 1.000000e+00, 5.000000e-01 + 71, 0.000000e+00, 0.000000e+00, 5.000000e-01 + 72, 5.000000e-01, 0.000000e+00, 5.000000e-01 + 73, 1.000000e+00, 5.000000e-01, 5.000000e-01 + 74, 1.000000e+00, 1.000000e+00, 5.000000e-01 + 75, 1.000000e+00, 0.000000e+00, 5.000000e-01 + 76, 1.500000e+00, 5.000000e-01, 5.000000e-01 + 77, 1.500000e+00, 1.000000e+00, 5.000000e-01 + 78, 1.500000e+00, 0.000000e+00, 5.000000e-01 + 79, 2.000000e+00, 5.000000e-01, 5.000000e-01 + 80, 2.000000e+00, 1.000000e+00, 5.000000e-01 + 81, 2.000000e+00, 0.000000e+00, 5.000000e-01 + 82, 2.500000e+00, 5.000000e-01, 5.000000e-01 + 83, 2.500000e+00, 1.000000e+00, 5.000000e-01 + 84, 2.500000e+00, 0.000000e+00, 5.000000e-01 + 85, 3.000000e+00, 5.000000e-01, 5.000000e-01 + 86, 3.000000e+00, 1.000000e+00, 5.000000e-01 + 87, 3.000000e+00, 0.000000e+00, 5.000000e-01 + 88, 3.500000e+00, 5.000000e-01, 5.000000e-01 + 89, 3.500000e+00, 1.000000e+00, 5.000000e-01 + 90, 3.500000e+00, 0.000000e+00, 5.000000e-01 + 91, 4.000000e+00, 5.000000e-01, 5.000000e-01 + 92, 4.000000e+00, 1.000000e+00, 5.000000e-01 + 93, 4.000000e+00, 0.000000e+00, 5.000000e-01 + 94, 4.500000e+00, 5.000000e-01, 5.000000e-01 + 95, 4.500000e+00, 1.000000e+00, 5.000000e-01 + 96, 4.500000e+00, 0.000000e+00, 5.000000e-01 + 97, 5.000000e+00, 5.000000e-01, 5.000000e-01 + 98, 5.000000e+00, 1.000000e+00, 5.000000e-01 + 99, 5.000000e+00, 0.000000e+00, 5.000000e-01 + 100, 5.500000e+00, 0.000000e+00, 5.000000e-01 + 101, 5.500000e+00, 0.000000e+00, 2.500000e-01 + 102, 5.500000e+00, 5.000000e-01, 5.000000e-01 + 103, 5.500000e+00, 5.000000e-01, 2.500000e-01 + 104, 6.000000e+00, 0.000000e+00, 5.000000e-01 + 105, 6.000000e+00, 0.000000e+00, 2.500000e-01 + 106, 6.000000e+00, 5.000000e-01, 5.000000e-01 + 107, 6.000000e+00, 5.000000e-01, 2.500000e-01 + 108, 5.500000e+00, 0.000000e+00, 0.000000e+00 + 109, 5.500000e+00, 5.000000e-01, 0.000000e+00 + 110, 6.000000e+00, 0.000000e+00, 0.000000e+00 + 111, 6.000000e+00, 5.000000e-01, 0.000000e+00 + 112, 5.500000e+00, 1.000000e+00, 5.000000e-01 + 113, 5.500000e+00, 1.000000e+00, 2.500000e-01 + 114, 6.000000e+00, 1.000000e+00, 5.000000e-01 + 115, 6.000000e+00, 1.000000e+00, 2.500000e-01 + 116, 5.500000e+00, 1.000000e+00, 0.000000e+00 + 117, 6.000000e+00, 1.000000e+00, 0.000000e+00 + 118, 6.847692e+00, 0.000000e+00, 5.000000e-01 + 119, 6.847692e+00, 0.000000e+00, 2.500000e-01 + 120, 6.736928e+00, 5.330372e-01, 5.000000e-01 + 121, 6.736928e+00, 5.330372e-01, 2.500000e-01 + 122, 7.695384e+00, 0.000000e+00, 5.000000e-01 + 123, 7.695384e+00, 0.000000e+00, 2.500000e-01 + 124, 7.460065e+00, 6.306936e-01, 5.000000e-01 + 125, 7.460065e+00, 6.306936e-01, 2.500000e-01 + 126, 8.543076e+00, 0.000000e+00, 5.000000e-01 + 127, 8.543076e+00, 0.000000e+00, 2.500000e-01 + 128, 8.156227e+00, 7.886674e-01, 5.000000e-01 + 129, 8.156227e+00, 7.886674e-01, 2.500000e-01 + 130, 9.390769e+00, 0.000000e+00, 5.000000e-01 + 131, 9.390769e+00, 0.000000e+00, 2.500000e-01 + 132, 8.813418e+00, 1.000000e+00, 5.000000e-01 + 133, 8.813418e+00, 1.000000e+00, 2.500000e-01 + 134, 6.847692e+00, 0.000000e+00, 0.000000e+00 + 135, 6.736928e+00, 5.330372e-01, 0.000000e+00 + 136, 7.695384e+00, 0.000000e+00, 0.000000e+00 + 137, 7.460065e+00, 6.306936e-01, 0.000000e+00 + 138, 8.543076e+00, 0.000000e+00, 0.000000e+00 + 139, 8.156227e+00, 7.886674e-01, 0.000000e+00 + 140, 9.390769e+00, 0.000000e+00, 0.000000e+00 + 141, 8.813418e+00, 1.000000e+00, 0.000000e+00 + 142, 6.626164e+00, 1.066074e+00, 5.000000e-01 + 143, 6.626164e+00, 1.066074e+00, 2.500000e-01 + 144, 7.224745e+00, 1.261387e+00, 5.000000e-01 + 145, 7.224745e+00, 1.261387e+00, 2.500000e-01 + 146, 7.769377e+00, 1.577335e+00, 5.000000e-01 + 147, 7.769377e+00, 1.577335e+00, 2.500000e-01 + 148, 8.236068e+00, 2.000000e+00, 5.000000e-01 + 149, 8.236068e+00, 2.000000e+00, 2.500000e-01 + 150, 6.626164e+00, 1.066074e+00, 0.000000e+00 + 151, 7.224745e+00, 1.261387e+00, 0.000000e+00 + 152, 7.769377e+00, 1.577335e+00, 0.000000e+00 + 153, 8.236068e+00, 2.000000e+00, 0.000000e+00 + 154, 1.047214e+01, 0.000000e+00, 0.000000e+00 + 155, 1.047214e+01, 1.500000e+00, 0.000000e+00 + 156, 9.589422e+00, 1.369306e+00, 0.000000e+00 + 157, 9.931452e+00, 0.000000e+00, 0.000000e+00 + 158, 1.047214e+01, 0.000000e+00, 2.500000e-01 + 159, 1.047214e+01, 1.500000e+00, 2.500000e-01 + 160, 9.589422e+00, 1.369306e+00, 2.500000e-01 + 161, 9.931452e+00, 0.000000e+00, 2.500000e-01 + 162, 1.047214e+01, 3.000000e+00, 0.000000e+00 + 163, 9.247391e+00, 2.738613e+00, 0.000000e+00 + 164, 1.047214e+01, 3.000000e+00, 2.500000e-01 + 165, 9.247391e+00, 2.738613e+00, 2.500000e-01 + 166, 1.047214e+01, 0.000000e+00, 5.000000e-01 + 167, 1.047214e+01, 1.500000e+00, 5.000000e-01 + 168, 9.589422e+00, 1.369306e+00, 5.000000e-01 + 169, 9.931452e+00, 0.000000e+00, 5.000000e-01 + 170, 1.047214e+01, 3.000000e+00, 5.000000e-01 + 171, 9.247391e+00, 2.738613e+00, 5.000000e-01 + 172, 1.148607e+01, 0.000000e+00, 5.000000e-01 + 173, 1.148607e+01, 0.000000e+00, 2.500000e-01 + 174, 1.148607e+01, 1.500000e+00, 5.000000e-01 + 175, 1.148607e+01, 1.500000e+00, 2.500000e-01 + 176, 1.250000e+01, 0.000000e+00, 5.000000e-01 + 177, 1.250000e+01, 0.000000e+00, 2.500000e-01 + 178, 1.250000e+01, 1.500000e+00, 5.000000e-01 + 179, 1.250000e+01, 1.500000e+00, 2.500000e-01 + 180, 1.148607e+01, 0.000000e+00, 0.000000e+00 + 181, 1.148607e+01, 1.500000e+00, 0.000000e+00 + 182, 1.250000e+01, 0.000000e+00, 0.000000e+00 + 183, 1.250000e+01, 1.500000e+00, 0.000000e+00 + 184, 1.148607e+01, 3.000000e+00, 5.000000e-01 + 185, 1.148607e+01, 3.000000e+00, 2.500000e-01 + 186, 1.250000e+01, 3.000000e+00, 5.000000e-01 + 187, 1.250000e+01, 3.000000e+00, 2.500000e-01 + 188, 1.148607e+01, 3.000000e+00, 0.000000e+00 + 189, 1.250000e+01, 3.000000e+00, 0.000000e+00 + 190, 1.375000e+01, 0.000000e+00, 5.000000e-01 + 191, 1.375000e+01, 0.000000e+00, 2.500000e-01 + 192, 1.375000e+01, 1.500000e+00, 5.000000e-01 + 193, 1.375000e+01, 1.500000e+00, 2.500000e-01 + 194, 1.500000e+01, 0.000000e+00, 5.000000e-01 + 195, 1.500000e+01, 0.000000e+00, 2.500000e-01 + 196, 1.500000e+01, 1.500000e+00, 5.000000e-01 + 197, 1.500000e+01, 1.500000e+00, 2.500000e-01 + 198, 1.625000e+01, 0.000000e+00, 5.000000e-01 + 199, 1.625000e+01, 0.000000e+00, 2.500000e-01 + 200, 1.625000e+01, 1.500000e+00, 5.000000e-01 + 201, 1.625000e+01, 1.500000e+00, 2.500000e-01 + 202, 1.750000e+01, 0.000000e+00, 5.000000e-01 + 203, 1.750000e+01, 0.000000e+00, 2.500000e-01 + 204, 1.750000e+01, 1.500000e+00, 5.000000e-01 + 205, 1.750000e+01, 1.500000e+00, 2.500000e-01 + 206, 1.375000e+01, 0.000000e+00, 0.000000e+00 + 207, 1.375000e+01, 1.500000e+00, 0.000000e+00 + 208, 1.500000e+01, 0.000000e+00, 0.000000e+00 + 209, 1.500000e+01, 1.500000e+00, 0.000000e+00 + 210, 1.625000e+01, 0.000000e+00, 0.000000e+00 + 211, 1.625000e+01, 1.500000e+00, 0.000000e+00 + 212, 1.750000e+01, 0.000000e+00, 0.000000e+00 + 213, 1.750000e+01, 1.500000e+00, 0.000000e+00 + 214, 1.375000e+01, 3.000000e+00, 5.000000e-01 + 215, 1.375000e+01, 3.000000e+00, 2.500000e-01 + 216, 1.500000e+01, 3.000000e+00, 5.000000e-01 + 217, 1.500000e+01, 3.000000e+00, 2.500000e-01 + 218, 1.625000e+01, 3.000000e+00, 5.000000e-01 + 219, 1.625000e+01, 3.000000e+00, 2.500000e-01 + 220, 1.750000e+01, 3.000000e+00, 5.000000e-01 + 221, 1.750000e+01, 3.000000e+00, 2.500000e-01 + 222, 1.375000e+01, 3.000000e+00, 0.000000e+00 + 223, 1.500000e+01, 3.000000e+00, 0.000000e+00 + 224, 1.625000e+01, 3.000000e+00, 0.000000e+00 + 225, 1.750000e+01, 3.000000e+00, 0.000000e+00 +** +********************************** E L E M E N T S **************************** +*ELEMENT, TYPE=C3D8R, ELSET=EB1 + 1, 1, 2, 3, 4, 5, 6, 7, 8 + 2, 2, 9, 10, 3, 6, 11, 12, 7 + 3, 4, 3, 13, 14, 8, 7, 15, 16 + 4, 3, 10, 17, 13, 7, 12, 18, 15 + 5, 14, 13, 19, 20, 16, 15, 21, 22 + 6, 13, 17, 23, 19, 15, 18, 24, 21 + 7, 20, 19, 25, 26, 22, 21, 27, 28 + 8, 19, 23, 29, 25, 21, 24, 30, 27 + 9, 26, 25, 31, 32, 28, 27, 33, 34 + 10, 25, 29, 35, 31, 27, 30, 36, 33 + 11, 32, 31, 37, 38, 34, 33, 39, 40 + 12, 31, 35, 41, 37, 33, 36, 42, 39 + 13, 38, 37, 43, 44, 40, 39, 45, 46 + 14, 37, 41, 47, 43, 39, 42, 48, 45 + 15, 44, 43, 49, 50, 46, 45, 51, 52 + 16, 43, 47, 53, 49, 45, 48, 54, 51 + 17, 50, 49, 55, 56, 52, 51, 57, 58 + 18, 49, 53, 59, 55, 51, 54, 60, 57 + 19, 56, 55, 61, 62, 58, 57, 63, 64 + 20, 55, 59, 65, 61, 57, 60, 66, 63 + 21, 5, 6, 7, 8, 67, 68, 69, 70 + 22, 6, 11, 12, 7, 68, 71, 72, 69 + 23, 8, 7, 15, 16, 70, 69, 73, 74 + 24, 7, 12, 18, 15, 69, 72, 75, 73 + 25, 16, 15, 21, 22, 74, 73, 76, 77 + 26, 15, 18, 24, 21, 73, 75, 78, 76 + 27, 22, 21, 27, 28, 77, 76, 79, 80 + 28, 21, 24, 30, 27, 76, 78, 81, 79 + 29, 28, 27, 33, 34, 80, 79, 82, 83 + 30, 27, 30, 36, 33, 79, 81, 84, 82 + 31, 34, 33, 39, 40, 83, 82, 85, 86 + 32, 33, 36, 42, 39, 82, 84, 87, 85 + 33, 40, 39, 45, 46, 86, 85, 88, 89 + 34, 39, 42, 48, 45, 85, 87, 90, 88 + 35, 46, 45, 51, 52, 89, 88, 91, 92 + 36, 45, 48, 54, 51, 88, 90, 93, 91 + 37, 52, 51, 57, 58, 92, 91, 94, 95 + 38, 51, 54, 60, 57, 91, 93, 96, 94 + 39, 58, 57, 63, 64, 95, 94, 97, 98 + 40, 57, 60, 66, 63, 94, 96, 99, 97 + 41, 99, 100, 101, 66, 97, 102, 103, 63 + 42, 100, 104, 105, 101, 102, 106, 107, 103 + 43, 66, 101, 108, 65, 63, 103, 109, 61 + 44, 101, 105, 110, 108, 103, 107, 111, 109 + 45, 97, 102, 103, 63, 98, 112, 113, 64 + 46, 102, 106, 107, 103, 112, 114, 115, 113 + 47, 63, 103, 109, 61, 64, 113, 116, 62 + 48, 103, 107, 111, 109, 113, 115, 117, 116 + 49, 104, 118, 119, 105, 106, 120, 121, 107 + 50, 118, 122, 123, 119, 120, 124, 125, 121 + 51, 122, 126, 127, 123, 124, 128, 129, 125 + 52, 126, 130, 131, 127, 128, 132, 133, 129 + 53, 105, 119, 134, 110, 107, 121, 135, 111 + 54, 119, 123, 136, 134, 121, 125, 137, 135 + 55, 123, 127, 138, 136, 125, 129, 139, 137 + 56, 127, 131, 140, 138, 129, 133, 141, 139 + 57, 106, 120, 121, 107, 114, 142, 143, 115 + 58, 120, 124, 125, 121, 142, 144, 145, 143 + 59, 124, 128, 129, 125, 144, 146, 147, 145 + 60, 128, 132, 133, 129, 146, 148, 149, 147 + 61, 107, 121, 135, 111, 115, 143, 150, 117 + 62, 121, 125, 137, 135, 143, 145, 151, 150 + 63, 125, 129, 139, 137, 145, 147, 152, 151 + 64, 129, 133, 141, 139, 147, 149, 153, 152 + 65, 154, 155, 156, 157, 158, 159, 160, 161 + 66, 155, 162, 163, 156, 159, 164, 165, 160 + 67, 157, 156, 141, 140, 161, 160, 133, 131 + 68, 156, 163, 153, 141, 160, 165, 149, 133 + 69, 158, 159, 160, 161, 166, 167, 168, 169 + 70, 159, 164, 165, 160, 167, 170, 171, 168 + 71, 161, 160, 133, 131, 169, 168, 132, 130 + 72, 160, 165, 149, 133, 168, 171, 148, 132 + 73, 166, 172, 173, 158, 167, 174, 175, 159 + 74, 172, 176, 177, 173, 174, 178, 179, 175 + 75, 158, 173, 180, 154, 159, 175, 181, 155 + 76, 173, 177, 182, 180, 175, 179, 183, 181 + 77, 167, 174, 175, 159, 170, 184, 185, 164 + 78, 174, 178, 179, 175, 184, 186, 187, 185 + 79, 159, 175, 181, 155, 164, 185, 188, 162 + 80, 175, 179, 183, 181, 185, 187, 189, 188 + 81, 176, 190, 191, 177, 178, 192, 193, 179 + 82, 190, 194, 195, 191, 192, 196, 197, 193 + 83, 194, 198, 199, 195, 196, 200, 201, 197 + 84, 198, 202, 203, 199, 200, 204, 205, 201 + 85, 177, 191, 206, 182, 179, 193, 207, 183 + 86, 191, 195, 208, 206, 193, 197, 209, 207 + 87, 195, 199, 210, 208, 197, 201, 211, 209 + 88, 199, 203, 212, 210, 201, 205, 213, 211 + 89, 178, 192, 193, 179, 186, 214, 215, 187 + 90, 192, 196, 197, 193, 214, 216, 217, 215 + 91, 196, 200, 201, 197, 216, 218, 219, 217 + 92, 200, 204, 205, 201, 218, 220, 221, 219 + 93, 179, 193, 207, 183, 187, 215, 222, 189 + 94, 193, 197, 209, 207, 215, 217, 223, 222 + 95, 197, 201, 211, 209, 217, 219, 224, 223 + 96, 201, 205, 213, 211, 219, 221, 225, 224 +** +********************************** S I D E S E T S ********************************** +*ELSET, ELSET=SS1_S3 + 1, 2, 21, 22 +*SURFACE, NAME=SS1 +SS1_S3, S3 +*ELSET, ELSET=SS3_S1 + 41, 42, 43, 44, 49, 50, 51, 52, 53, 54, + 55, 56, 73, 74, 75, 76, 81, 82, 83, 84, + 85, 86, 87, 88 +*ELSET, ELSET=SS3_S4 + 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, + 22, 24, 26, 28, 30, 32, 34, 36, 38, 40 +*ELSET, ELSET=SS3_S6 + 65, 67, 69, 71 +*SURFACE, NAME=SS3 +SS3_S1, S1 +SS3_S4, S4 +SS3_S6, S6 +*ELSET, ELSET=SS5_S1 + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, + 65, 66, 67, 68 +*ELSET, ELSET=SS5_S5 + 43, 44, 47, 48, 53, 54, 55, 56, 61, 62, + 63, 64, 75, 76, 79, 80, 85, 86, 87, 88, + 93, 94, 95, 96 +*SURFACE, NAME=SS5 +SS5_S1, S1 +SS5_S5, S5 +*ELSET, ELSET=SS6_S3 + 81, 82, 83, 84, 89, 90, 91, 92 +*SURFACE, NAME=SS6 +SS6_S3, S3 +*ELSET, ELSET=SS10_S2 + 57, 58, 59, 60, 61, 62, 63, 64 +*SURFACE, NAME=SS10 +SS10_S2, S2 +*ELSET, ELSET=SS11_S4 + 66, 68, 70, 72 +*SURFACE, NAME=SS11 +SS11_S4, S4 +** +********************************** P R O P E R T I E S ************************ +*SOLID SECTION, ELSET=EB1, MATERIAL=Default-Steel +** +*END PART +** +** +** +********************************** E N D P A R T S ********************************** +** +** +********************************** A S S E M B L Y ************************************ +** +*ASSEMBLY, NAME=ASSEMBLY1 +** +*INSTANCE, NAME=Part-Default_1, PART=Part-Default +*END INSTANCE +** +*END ASSEMBLY +** +** +** +*MATERIAL, NAME = Default-Steel +*ELASTIC, TYPE=ISOTROPIC +2.068000e+05, 2.900000e-01 +*DENSITY +7.000000e-06 +*CONDUCTIVITY,TYPE=ISO +4.500000e-02 +*SPECIFIC HEAT +5.000000e+02 +** +** +************************************** H I S T O R Y ************************************* +** +*PREPRINT +** +**************************************** S T E P 1 *************************************** +*STEP,INC=100,NAME=Default Set +** +*STATIC +1, 1, 1e-05, 1 +** +** +** +** +*END STEP diff --git a/Linear_Elastoplasticity/parameters.prm b/Linear_Elastoplasticity/parameters.prm new file mode 100644 index 00000000..540a3cc1 --- /dev/null +++ b/Linear_Elastoplasticity/parameters.prm @@ -0,0 +1,73 @@ +# Listing of Parameters +# --------------------- +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Boundary conditions + # The amount by which the specimen is to be stretched + set Final stretch = 1.2 +end + + +subsection Geometry + # Geometry to simulate + # Options: Notched cylinder ; Tensile specimen + set Geometry type = Tensile specimen + + # Global grid scaling factor + set Grid scale = 1e-3 + + # Number of global refinement steps + set Global refinement steps = 2 + + # Number of initial local refinement cycles + set Local refinement steps = 1 + + subsection Notched cylinder + # Overall length of the specimen + set Length = 100 + + # Overall radius of the specimen + set Radius = 10 + + # Overall length of the notch in the specimen + set Notch length = 1 + + # Overall radius of the notch in the specimen + set Notch radius = 9 + end +end + + +subsection Material properties + # Type of material that composes the specimen + # Options: Elastic ; Elastoplastic (isotropic hardening) + set Type = Elastoplastic (isotropic hardening) + + # Shear modulus + set Shear modulus = 100 + + # Poisson ratio + set Poisson ratio = 0.3 + + # Isotropic hardening constant + set Isotropic hardening constant = 10 + + # Yield stress + set Yield stress = 20 +end + + +subsection Time + # End time + set End time = 1.0 + + # Time step size + set Time step size = 0.1 +end