--- /dev/null
+\documentclass{article}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+
+\title{Documentation of step-42, An obstacle problem for elasto-plastic material behavior in three dimensions}
+\author{Joerg Frohne}
+\date{Juni, 2012}
+
+\begin{document}
+
+\section{Introduction}
+
+This example composes an advanced version of step-41 since it considers an
+elasto-plastic material behavior with isotropic hardening in three dimensions.
+That means that we have to take care of an additional nonlinearity by the
+material behavior. An other difference compared to step-41 is that
+the contact area is arranged at the boundary of the deformable body now.\\
+Since you can slightly reach a few million degrees of freedom in three dimensions,
+even with adaptive mesh refinement, we decided to use trilinos and p4est to run
+our code in parallel. On the other hand we have to deal with hanging nodes because
+of the adaptive mesh, which is an other advance of step-41.
+
+
+\section{Classical formulation}
+
+The classical formulation of the problem possesses the following form:
+\begin{align*}
+ \varepsilon(u) &= A\sigma + \lambda & &\quad\text{in } \Omega,\\
+ \lambda(\tau - \sigma) &\geq 0\quad\forall\tau\text{ mit }\mathcal{F}(\tau)\leq 0 & &\quad\text{in } \Omega,\\
+ -\textrm{div}\ \sigma &= f & &\quad\text{in } \Omega,\\
+ u(\mathbf x) &= 0 & &\quad\text{on }\Gamma_D,\\
+ \sigma_t(u) &= 0,\quad\sigma_n(u)\leq 0 & &\quad\text{on }\Gamma_C,\\
+\sigma_n(u)(u_n - g) &= 0,\quad u_n(\mathbf x) - g(\mathbf x) \leq 0 & &\quad\text{on } \Gamma_C
+\end{align*}
+with $u\in H^2(\Omega)$. The vector valued function $u$ denotes the
+displacement in the deformable body. The first two lines describe the elast-plastic
+material behaviour. Therein the equation shows the deformation $\varepsilon (u)$ as the additive
+decomposition of the elastic part $A\sigma$ and the plastic part $\lambda$. $A$ is defined as
+the compliance tensor of fourth order which contains some material constants and $\sigma$ as the
+symmetric stress tensor of second order. So we have to consider the inequality in the second
+row component-by-component and furthermore we have to distinguish two cases.\\
+The continuous and convex function $\mathcal{F}$ denotes the von mises flow function
+$$\mathcal{F}(\tau) = \vert\tau^D\vert - \sigma_0$$
+with $\sigma_0$ as yield stress. If there is no plastic deformation - that is $\lambda=0$ - this yields $\vert\sigma^D\vert < \sigma_0$
+and otherwise if $\lambda > 0$ it follows that $\vert\sigma^D\vert = \sigma_0$. That means if the stress is smaller as the yield stress
+there are only elastic deformations. Therein the Index $D$ denotes the deviator part of the stress $\sigma$ which
+is dedined as
+$$\sigma^D = \sigma - \dfrac{1}{3}tr(\sigma).$$
+It describes the hydrostatic part of the stress tensor in contrast to the volumetric part. For metal the hydrostatic
+stress composes the main indicator for the plastic deformation.\\
+The second equation is called equilibrium condition with a force of areal density $f$ which we will neglect in our example.
+The boundary of $\Omega$ separates as follows $\Gamma=\Gamma_D\bigcup\Gamma_C$ and $\Gamma_D\bigcap\Gamma_C=\emptyset$.
+At the boundary $\Gamma_D$ we have zero Dirichlet conditions. $\Gamma_C$ denotes the potential contact boundary.\\
+The last two lines decribe the so-called Signorini contact conditions. If there is no contact the normal stress
+$$ \sigma_n = \sigma n\cdot n$$
+is zero with the outward normal $n$. If there is contact ($u_n = g$) the tangential stress $\sigma_t = \sigma\cdot n - \sigma_n n$
+vanishes, because we consider a frictionless situation and the normal stress is negative.
+
+\section{Derivation of the variational inequality}
+
+As a starting point we want to minimise an energy functional:
+$$E(\tau) := \dfrac{1}{2}\int\limits_{\Omega}\tau A \tau d\tau,\quad \tau\in \Pi W^{div}$$
+with
+$$W^{div}:=\lbrace \tau\in L^2(\Omega,\mathbb{R}^{dim\times\dim}_{sym}),div(\tau)\in L^2(\Omega,\mathbb{R}^{dim})\rbrace$$
+and
+$$\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace$$
+as the set of admissible stresses which is defined
+by a continious, convex flow function $\mathcal{F}$.
+
+With the goal to derive the dual formulation of the minimisation problem, we define a lagrange function:
+$$L(\tau,\varphi) := E(\tau) + (\varphi, div(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{div}\times U$$
+with $U := \lbrace u\in H^1(\Omega), u = g \text{ on } \Gamma_D,u_n\leq 0 \text{ on } \Gamma_C \rbrace$.\\
+By building the fr\'echet derivatives of $L$ for both components we obtain the dual formulation for the stationary case
+which is known as \textbf{Hencky-Type-Model}:\\
+Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times U$ with
+$$\left(A\sigma,\tau - \sigma\right) + \left(u, div(\tau) - div(\sigma)\right) \geq 0,\quad \forall \tau\in \Pi W^{div}$$
+$$-\left(div(\sigma),\varphi - u\right) \geq 0,\quad \forall \varphi\in U.$$
+By integrating by parts and multiplying the first inequality by $C=A^{-1}$ we achieve the primal-mixed version of our problem:\\
+Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times U$ with
+$$\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W$$
+$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U.$$
+Therein $\varepsilon$ denotes the linearised deformation tensor with $\varepsilon(u) := \dfrac{1}{2}\left(\nabla u + \nabla u^T\right)$ for small deformations.\\
+Most materials - especially metals - have the property that they show some hardening effects during the forming process.
+There are different constitutive laws to describe those material behaviour. The most simple one is called linear isotropic hardening
+with the flow function $\mathcal{F}(\tau,\eta) = \vert\tau^D\vert - (\sigma_0 + \gamma\eta)$.
+It can be considered by establishing an additional term in our primal-mixed formulation:\\
+Find a pair $\lbrace(\sigma,\xi),u\rbrace\in \Pi (W\times L^2(\Omega,\mathbb{R}))\times U$ with
+$$\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) + \gamma\left( \xi, \eta - \xi\right) \geq 0,\quad \forall (\tau,\eta)\in \Pi (W,L^2(\Omega,\mathbb{R}))$$
+$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U,$$
+with the hardening parameter $\gamma > 0$.\\
+Now we want to derive a primal problem which only depends on the displacement $u$. For that purpose we
+set $\eta = \xi$ and eliminate the stress $\sigma$ by applying the projection theorem on\\
+$$\left(\sigma - C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W,$$
+which yields with the second inequality:\\
+Find the displacement $u\in U$ with
+$$\left(P_{\Pi}(C\varepsilon(u)),\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U,$$
+with the projection:
+$$P_{\Pi}(\tau):=\begin{cases}
+ \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0 + \gamma\xi,\\
+ \hat\alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0 + \gamma\xi,
+ \end{cases}$$
+with the radius
+$$\hat\alpha := \sigma_0 + \gamma\xi .$$
+With the relation $\xi = \vert\varepsilon(u) - A\sigma\vert$ it is possible to eliminate $\xi$ inside the projection $P_{\Pi}$:\\
+$$P_{\Pi}(\tau):=\begin{cases}
+ \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\
+ \alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0,
+ \end{cases}$$
+$$\alpha := \sigma_0 + \dfrac{\gamma}{2\mu+\gamma}\left(\vert\tau^D\vert - \sigma_0\right) ,$$
+with a further material parameter $\mu>0$ called shear modulus.\\
+So what we do is to calculate the stresses by using Hooke's law for linear elastic, isotropic materials
+$$\sigma = C \varepsilon(u) = 2\mu \varepsilon^D(u) + \kappa tr(\varepsilon(u))I = \left[2\mu\left(\mathbb{I} -\dfrac{1}{3} I\otimes I\right) + \kappa I\otimes I\right]\varepsilon(u)$$
+with the new material parameter $\kappa>0$ (bulk modulus). The variables $I$ and $\mathbb{I}$ denote the identity tensors of second and forth order.\\
+In the next step we test in a pointwise sense where the deviator part of the stress in a norm is bigger as the yield stress.
+If there are such points we project the deviator stress in those points back to the yield surface. Methods of this kind
+are called projections algorithm or radial-return-algorithm.\\
+Now we have a primal formulation of our elasto-plastic contact problem which only depends on the displacement $u$.
+It consists of a nonlinear variational inequality and has a unique solution as it shows the theorem of Lions and Stampaccia
+(A proof can be found in Rodrigues: Obstacle Problems in Mathematical Physics, North-Holland, Amsterdam, 1987).\\
+To handle the nonlinearity of the constitutive law we use a newton method and to deal with the contact we apply an
+active set method like in step-41. To be more concrete we combine both methods to an inexact semi smooth newton
+method - inexact since we use an iterative solver for the linearised problems in each newton step.
+
+\section{Linearisation of the constitutive law for the newton method}
+
+For the newton method we have to linearise the following semi-linearform
+$$a(\psi;\varphi) := \left(P_{\Pi}(C\varepsilon(\varphi)),\varepsilon(\varphi)\right).$$
+Becaus we have to find the solution $u$ in the convex set $U$, we have to apply an SQP-method (SQP: sequential quadratic
+programming). That means we have to solve a minimisation problem for a known $u^i$ in every SQP-step of the form
+\begin{eqnarray*}
+ & & a(u^{i};u^{i+1} - u^i) + \dfrac{1}{2}a'(u^i;u^{i+1} - u^i,u^{i+1} - u^i)\\
+ &=& a(u^i;u^{i+1}) - a(u^i;u^i) +\\
+ & & \dfrac{1}{2}\left( a'(u^i;u^{i+1},u^{i+1}) - 2a'(u^i;u^i,u^{i+1}) - a'(u^i;u^i,u^i)\right)\\
+ &\rightarrow& min,\quad u^{i+1}\in U.
+\end{eqnarray*}
+Neglecting the constant terms $ a(u^i;u^i)$ and $ a'(u^i;u^i,u^i)$ we obtain the following minimisation problem
+$$\dfrac{1}{2} a'(u^i;u^{i+1},u^{i+1}) - F(u^i)\rightarrow min,\quad u^{i+1}\in U$$
+with
+$$F(\varphi) := \left(a'(\varphi;\varphi,u^{i+1}) - a(\varphi;u^{i+1}) \right).$$
+In the case of our constitutive law the derivitive of the semi-linearform $a(.;.)$ at the point $u^i$ is
+
+$$a'(u^i;\psi,\varphi) =$$
+$$
+\begin{cases}
+\left(\left[2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes I\right) + \kappa I\otimes I\right]\varepsilon(\psi),\varepsilon(\varphi)\right), & \quad \vert\tau^D\vert \leq \sigma_0\\
+\left(\left[\dfrac{\alpha}{\vert\tau^D\vert}2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes I - \dfrac{\tau^D\otimes\tau^D}{\vert\tau^D\vert}\right) + \kappa I\otimes I\right]\varepsilon(\psi),\varepsilon(\varphi) \right), & \quad \vert\tau^D\vert > \sigma_0
+\end{cases}
+$$
+with
+$$\tau^D := C\varepsilon^D(u^i).$$
+Again the first case is for elastic and the second for plastic deformation.
+
+\section{Formulation as a saddle point problem}
+
+On the line of step-41 we compose a saddle point problem out of the minimisation problem. Again we do so to gain a formulation
+that allows us to solve a linear system of equations finally.
+
+\section{Active Set methods to solve the saddle point problem}
+
+\section{The primal-dual active set algorithm combined with the inexact semi smooth newton method}
+
+The inexact newton method works as follows:
+\begin{itemize}
+ \item[(0)] Initialize $\mathcal{A}_k$ and $\mathcal{F}_k$, such that $\mathcal{S} = \mathcal{A}_k \cup \mathcal{F}_k$ and $\mathcal{A}_k \cap \mathcal{F}_k = \emptyset$ and set $k = 1$.
+ \item[(1)] Assembel the newton matrix $a'(U^k;\varphi_i,\varphi_j)$ and the right-hand-side $F(U^k)$.
+ \item[(2)] Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies
+ \begin{align*}
+ AU^k + B\Lambda^k & = F, &\\
+ \left[BU^k\right]_i & = G & & \forall i\in\mathcal{A}_k\\
+ \Lambda^k_i & = 0 & & \forall i\in\mathcal{F}_k.
+ \end{align*}
+% Note that $\mathcal{S}$ contains only dofs related to the boundary $\Gamma_C$. So in contrast to step-41 there are much more than $\vert \mathcal{S}\vert$ equations necessary to determine $U$ and $\Lambda$.
+ \item[(3)] Define the new active and inactive sets by
+ $$\mathcal{A}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c\left(\left[BU^k\right]_i - G_i\left) < 0\rbrace,$$
+ $$\mathcal{F}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c\left(\left[BU^k\right]_i - G_i\left) \geq 0\rbrace.$$
+ \item[(4)] If $\mathcal{A}_{k+1} = \mathcal{A}_k$ and $\vert F(U^{k+1}\vert < \delta$ then stop, else set $k=k+1$ and go to step (1).
+\end{itemize}
+
+\section{Implementation}
+
+\end{document}
\ No newline at end of file
projects / dealii-svn.git / commitdiff