From 855fb58f6373fd181bccc344108ceb9b9dafd470 Mon Sep 17 00:00:00 2001 From: frohne Date: Tue, 7 Aug 2012 14:59:24 +0000 Subject: [PATCH] a first version of the documentation git-svn-id: https://svn.dealii.org/trunk@25760 0785d39b-7218-0410-832d-ea1e28bc413d --- .../examples/step-42/doc/intro-step-42.tex | 182 ++++++++++++++++++ 1 file changed, 182 insertions(+) create mode 100644 deal.II/examples/step-42/doc/intro-step-42.tex diff --git a/deal.II/examples/step-42/doc/intro-step-42.tex b/deal.II/examples/step-42/doc/intro-step-42.tex new file mode 100644 index 0000000000..a7e1fce28d --- /dev/null +++ b/deal.II/examples/step-42/doc/intro-step-42.tex @@ -0,0 +1,182 @@ +\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 -- 2.39.5