单晶塑性abaqus-umat自学

本文以翻译黄永刚老师的umat文章为主， 这个颜色的为补充材料和个人的理解。

摘要

包含单晶塑性的abaqus用户材料子程序。本文回顾有限元形式的单晶弹塑性和黏弹塑性，包含小变形理论，以及严格的有限变形和有限旋转理论。单晶的非弹性变形由于晶体滑移，假定遵从Schmid准则。提出多种多样的分切应力与剪切变形的自硬化和潜硬化关系，并且包含在子程序内。

1. 简介

  有限元程序ABAQUS已经广泛应用在固体的变形与应力分析。除了自身大量的本构模型，ABAQUS还提供用户子程序接口。应力、应变以及求解依赖的状态变量以增量形式求出。当umat子程序被访问时，提供增量步开始的状态（应力以及求解依赖的状态变量）和应变和时间增量。umat子程序执行二段程序：更新应力和状态变量使其等于增量步结束时对应的值，提供材料的Jacobian矩阵，$\partial \Delta \sigma /\partial \Delta \epsilon$，为Newton-Rhapson迭代提供相应的本构模型。

  本文的主要目标是提供单晶塑性的本构的框架。运动学部分严格满足有限变形。塑性变形仅仅考虑位错滑移，不考虑扩散变形，孪晶和晶界滑移。Schmid应力，或滑移系上分切应力假定是滑移的驱动力。完整的子程序用于单晶和双晶的变形和断裂分析。

2. 单晶弹塑性本构模型回顾

2.1运动学

  晶体力学的运动学理论的先驱是Taylor（1938），精准的力学理论是Hill（1966），Rice（1971）和Hill and Rice（1972）搭建的。下面是关于Asaro and Rice（1977）和Asaro（1983）简单的理论总结。
  晶体材料的弹性变形和旋转是由晶格承担的。单晶的非弹性变形源于晶体滑移。材料流动是通过位错在晶格上的运动。总的变形梯度$\mathbf{F}$:$$\mathbf{F}={\mathbf{F}}^{\ast }\cdot {\mathbf{F}}^{\mathbf{P}}$$
  其中 $F^P$ 是中间构型的塑性剪切，$F^{\ast }$ 为晶格拉伸和旋转。假设弹性特征不受滑移的影响，应力由$F^{\ast }$ 单独决定。$F^P$ 的变化率与 $\alpha$ 滑移系的滑移率 ${\dot{\gamma }}^{(\alpha )}$ 有关。
$$\dot{{\mathbf{F}}^P}\cdot \dot{{\mathbf{F}}^{P-1}}=\sum _{\alpha }{\dot{\gamma }}^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}$$
  对所有激活的滑移系求和，$s^{(\alpha )}$ 和 $m^{(\alpha )}$ 分别是参考构型下滑移方向和滑移面的法向。

---

补充：单滑移变形梯度描述

速度梯度张量 $\mathbf{L}$
$$\begin{array}{rl}\mathbf{L}& =\dot{\mathbf{F}}\cdot {\mathbf{F}}^{-1}={\dot{\mathbf{F}}}^{\ast }\cdot {\mathbf{F}}^{\mathbf{P}}\cdot {\mathbf{F}}^{\mathbf{P}-1}\cdot {\mathbf{F}}^{\ast -1}+{\mathbf{F}}^{\ast }\cdot {\dot{\mathbf{F}}}^{\mathbf{P}}\cdot {\mathbf{F}}^{\mathbf{P}-1}\cdot {\mathbf{F}}^{\ast -1}\\ & ={\dot{\mathbf{F}}}^{\ast }\cdot {\mathbf{F}}^{\ast -1}+{\mathbf{F}}^{\ast }\cdot {\dot{\mathbf{F}}}^{\mathbf{P}}\cdot {\mathbf{F}}^{\mathbf{P}-1}\cdot {\mathbf{F}}^{\ast -1}\end{array}$$
求 ${\mathbf{F}}^{\mathbf{P}-1}$ 由于：
$$\begin{array}{rl}& [\mathbf{I}-{\gamma }^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}]\cdot [\mathbf{I}+{\gamma }^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}]\\ & \\ =& \mathbf{I}-{\gamma }^{(\alpha {)}^2}{s}^{(\alpha )}{m}^{(\alpha )}\cdot s^{(\alpha )}{m}^{(\alpha )}\end{array}$$
滑移方向和滑移面法向垂直：$m^{(\alpha )}\cdot s^{(\alpha )}=0$

因此： $${\mathbf{F}}^{\mathbf{P}-1}=\mathbf{I}-{\gamma }^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}$$
$$\dot{{\mathbf{F}}^{\mathbf{P}}}\cdot {\mathbf{F}}^{\mathbf{P}-1}=[{\dot{\gamma }}^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}]\cdot [\mathbf{I}-{\gamma }^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}]={\dot{\gamma }}^{(\alpha )}{s}^{(\alpha )}{m}^{(\alpha )}$$

---

  为了方便定义$s^{\ast (\alpha )}$ 为 $\alpha$ 滑移系在变形后的构型中的滑移方向：
$$s^{\ast (\alpha )}=F^{\ast }\cdot s^{(\alpha )}$$
  定义$m^{\ast (\alpha )}$ 为 $\alpha$ 滑移系在变形后的构型中的滑移面的法向：
$$m^{\ast (\alpha )}=m^{(\alpha )}\cdot F^{\ast -1}$$

---

因此：
$$\begin{array}{rl}\mathbf{L}& ={\dot{\mathbf{F}}}^{\ast }\cdot {\mathbf{F}}^{\ast -1}+{\dot{\gamma }}^{(\alpha )}{\mathbf{F}}^{\ast }\cdot s^{(\alpha )}{m}^{(\alpha )}\cdot {\mathbf{F}}^{\ast -1}\\ & \\ & ={\dot{\mathbf{F}}}^{\ast }\cdot {\mathbf{F}}^{\ast -1}+{\dot{\gamma }}^{(\alpha )}{\mathbf{s}}^{\ast (\alpha )}{\mathbf{m}}^{\ast (\alpha )}\\ & \\ & ={\mathbf{L}}^{\ast }+{\mathbf{L}}^{\mathbf{P}}\end{array}$$

---

  当前状态的速度梯度矩阵：
$$\mathbf{L}=\dot{\mathbf{F}}\cdot {\mathbf{F}}^{-1}=\mathbf{D}+\mathbf{\Omega }$$
  其中，$D$ 为对称的拉伸率，$\Omega$ 为反对称的自旋张量。可分解为晶格部分（上标*）和塑性部分（上标P）：
$$\mathbf{D}={\mathbf{D}}^{\ast }+{\mathbf{D}}^{\mathbf{P}}\mathbf{\Omega }={\mathbf{\Omega }}^{\ast }+{\mathbf{\Omega }}^{\mathbf{P}}$$
  满足：$${\mathbf{D}}^{\ast }+{\mathbf{\Omega }}^{\ast }={\dot{\mathbf{F}}}^{\ast }\cdot {\mathbf{F}}^{\ast -1}{\mathbf{D}}^{\mathbf{P}}+{\mathbf{\Omega }}^{\mathbf{P}}=\sum _{\alpha }{\dot{\gamma }}^{(\alpha )}{\mathbf{s}}^{\ast (\alpha )}{\mathbf{m}}^{\ast (\alpha )}$$

2.2 本构

  依据 Hill and Rice（1972），存在弹性势能，$\Phi =\Phi (F^{\ast })$ ，保证对称的晶格变形率张量 ${\mathbf{D}}^{\ast }$, 和Cauchy应力的Jaumann客观应变率满足如下关系：
$${\stackrel{\nabla }{\sigma }}^{\ast }+\mathbf{\sigma }\left(\mathbf{I}:{\mathbf{D}}^{\ast }\right)=\mathbf{L}:{\mathbf{D}}^{\ast }\text{\hspace{1em}(2.2.1)}$$
其中$I$是二阶单位张量，$\mathbf{L}$ 是弹性张量具有如下对称性$L_{ijkl}=L_{jikl}=L_{ijlk}=L_{klij}$，Jaumann率 ${\stackrel{\nabla }{\sigma }}^{\ast }$ 是晶格坐标轴上的共转应力率，Jaumann率 $\stackrel{\nabla }{\sigma }$ 是材料坐标轴上的共转应力率，二者之间的关系：
$${\stackrel{\nabla }{\sigma }}^{\ast }=\stackrel{\nabla }{\sigma }+(\Omega -{\Omega }^{\ast })\cdot \sigma -\sigma \cdot (\Omega -{\Omega }^{\ast })\text{\hspace{1em}(2.2.2)}$$
其中，$\stackrel{\nabla }{\sigma }=\dot{\mathbf{\sigma }}-\Omega \cdot \sigma +\sigma \cdot \Omega$

---

补充：

在Lagrange坐标下，变形率 $D^{\ast }$ 和 Kirchhoff应力 ${\tau }^{\ast }=J{\sigma }^{\ast }$ 共轭
$${\stackrel{\nabla }{\tau }}^{\ast }=\mathbf{L}:{\mathbf{D}}^{\ast }$$
其中：
$$\begin{array}{rl}{\stackrel{\nabla }{\tau }}^{\ast }& =\dot{{\tau }^{\ast }}-\Omega \cdot {\tau }^{\ast }+{\tau }^{\ast }\cdot \Omega \\ & =\dot{J}{\sigma }^{\ast }+J\dot{{\sigma }^{\ast }}-J\Omega \cdot {\sigma }^{\ast }+J{\sigma }^{\ast }\cdot \Omega \\ & =J{\sigma }^{\ast }(I:D)+J\dot{{\sigma }^{\ast }}-J\Omega \cdot {\sigma }^{\ast }+J{\sigma }^{\ast }\cdot \Omega \\ & =J\left[{\sigma }^{\ast }(I:D)+\stackrel{\nabla }{{\sigma }^{\ast }}\right]\end{array}$$
考虑到从初始构型到中间构型只发生滑移，忽略体积变化$J=\det F\approx 1$
因此：
$${\stackrel{\nabla }{\sigma }}^{\ast }+\mathbf{\sigma }\left(\mathbf{I}:{\mathbf{D}}^{\ast }\right)=\mathbf{L}:{\mathbf{D}}^{\ast }$$

---

  晶体滑移遵循Schmid准则，任意滑移系 $\alpha$的滑移率 ${\dot{\gamma }}^{(\alpha )}$假设只与当前应力 $\sigma$ 的Schmid应力 ${\tau }^{(\alpha )}$。Schmid应力是忽略晶格畸变后的分切应力。Asaro and Rice（1977）讨论了现有有限弹性畸变的一些可能的一般化。本文基于Rice（1971），热力学应力与滑移共轭。定义：
$${\tau }^{(\alpha )}=m^{\ast (\alpha )}\cdot \frac{{\rho }_0}{\rho }\sigma \cdot s^{\ast (\alpha )}\text{\hspace{1em}(2.2.3)}$$
其中 ${\rho }_0$ 和 $\rho$ 参考构型和当前构型的密度；Hill and Rice(1972)认为${\tau }^{(\alpha )}$等于随晶格旋转的Kirchhoff应力的最大剪切分量。Schmid应力的变化率如下：
$${\dot{\tau }}^{(\alpha )}=m^{\ast (\alpha )}\cdot \left[{\stackrel{\nabla }{\sigma }}^{\ast }+\mathbf{\sigma }\left(\mathbf{I}:{\mathbf{D}}^{\ast }\right)-{\mathbf{D}}^{\ast }\cdot \mathbf{\sigma }+\mathbf{\sigma }\cdot {\mathbf{D}}^{\ast }\right]\cdot s^{\ast (\alpha )}\text{\hspace{1em}(2.2.4)}$$

---

补充：

式（2.2.3）中：$$\frac{{\rho }_0}{\rho }=\frac{dv}{dV}=J$$
因此：
$${\dot{\tau }}^{(\alpha )}={\dot{m}}^{\ast (\alpha )}\cdot J\sigma \cdot s^{\ast (\alpha )}+m^{\ast (\alpha )}\cdot \dot{(J\sigma )}\cdot s^{\ast (\alpha )}+m^{\ast (\alpha )}\cdot J\sigma \cdot {\dot{s}}^{\ast (\alpha )}$$
其中：
$$\begin{array}{rl}{m}^{\ast (\alpha )}\cdot \dot{(J\sigma )}\cdot s^{\ast (\alpha )}& =m^{\ast (\alpha )}\cdot [\stackrel{\nabla }{(J\sigma )}+{\Omega }^{\ast }\cdot \sigma -\sigma \cdot {\Omega }^{\ast }]\cdot s^{\ast (\alpha )}\\ & =m^{\ast (\alpha )}\cdot [{\stackrel{\nabla }{\sigma }}^{\ast }+\mathbf{\sigma }\left(\mathbf{I}:{\mathbf{D}}^{\ast }\right)-{\Omega }^{\ast }\cdot \sigma +\sigma \cdot {\Omega }^{\ast }]\cdot s^{\ast (\alpha )}\end{array}$$
忽略体积变化$J=\det F\approx 1$
$$\begin{array}{rl}{\dot{m}}^{\ast (\alpha )}\cdot \sigma \cdot s^{\ast (\alpha )}& =m^{(\alpha )}\cdot {\dot{F}}^{\ast -1}\cdot \sigma \cdot s^{\ast (\alpha )}\\ & =m^{\ast (\alpha )}\cdot F^{\ast }\cdot {\dot{F}}^{\ast -1}\cdot \sigma \cdot s^{\ast (\alpha )}\\ & =m^{\ast (\alpha )}\cdot -L^{\ast }\cdot \sigma \cdot s^{\ast (\alpha )}\end{array}$$
注：$$\dot{\mathbf{I}}=\dot{({\mathbf{F}}^{\ast }\cdot {\mathbf{F}}^{\ast -1})}=\dot{{\mathbf{F}}^{\ast }}\cdot {\mathbf{F}}^{\ast -1}+{\mathbf{F}}^{\ast }\cdot {\dot{\mathbf{F}}}^{\ast -1}=0$$
其中：$\dot{{\mathbf{F}}^{\ast }}\cdot {\mathbf{F}}^{\ast -1}={\mathbf{L}}^{\ast }$ ,因此${\mathbf{F}}^{\ast }\cdot {\dot{\mathbf{F}}}^{\ast -1}=-{\mathbf{L}}^{\ast }$
同理：
$$\begin{array}{rl}{m}^{\ast (\alpha )}\cdot \sigma \cdot {\dot{s}}^{\ast (\alpha )}& =m^{(\ast \alpha )}\cdot \sigma \cdot \dot{F^{\ast }}\cdot s^{(\alpha )}\\ & =m^{\ast (\alpha )}\cdot \sigma \cdot \dot{F^{\ast }}\cdot F^{\ast -1}{s}^{\ast (\alpha )}\\ & =m^{\ast (\alpha )}\cdot \sigma \cdot L^{\ast }\cdot s^{\ast (\alpha )}\end{array}$$
因此：
$$\begin{array}{rl}{\dot{\tau }}^{(\alpha )}& =\\ & \end{array}$$