弹塑性力学问题

关键词：弹塑性力学 Mises 模型 平面应变 非线性

有限元方法是研究材料力学行为的一个有效工具。本节以二维理想弹塑性力学平面应变问题为例，讲解如何使用 FEtch 系统对弹塑性问题进行求解。

控制方程

微分形式

二维直角坐标系下，不考虑惯性力的作用，弹塑性平面问题的未知量服从如下偏微分方程。

平衡方程 \(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}+f_{x} = 0 \left(\text{ in } \Omega\right)\) \(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}+f_{y} = 0 \left(\text{ in } \Omega\right)\) 几何方程 \(\varepsilon_{x x}=\frac{\partial u}{\partial x}, \quad \varepsilon_{y y}=\frac{\partial v}{\partial y}, \quad \varepsilon_{x y}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\) 以上式子可以简记为 \(\nabla\boldsymbol{\cdot \sigma}+\boldsymbol{f}=\boldsymbol{0}\) \(\boldsymbol{\varepsilon}=\left(\nabla \boldsymbol{u}+\nabla^{T} \boldsymbol{u}\right) / 2\) 其中， $u$、$v$ 为位移 $\boldsymbol{u}$ 的分量，$\varepsilon_{x x}$、$\varepsilon_{y y}$ 、 $\varepsilon_{x y}$为应变 $\boldsymbol{\varepsilon}$ 的分量，$\sigma_{x x}$ 、$\sigma_{y y}$ 、$\sigma_{x y}$ 为应力 $\boldsymbol{\sigma}$ 的分量，$f_{x}$ 、$f_{y}$ 为体力 $\boldsymbol{f}$ 的分量。$\nabla$ 为梯度算子， $x$、$y$ 为直角坐标系下的坐标分量。

边界条件

考虑两类边界条件:

第一类边界条件，给定位移 \(u=u_0,\ v=v_0 \left(\text{ on } \Gamma_1\right)\)

第二类边界条件，给定外力 \(T_{x}=\left(\sigma_{x x},\sigma_{x y}\right)\cdot \mathbf{n}=t_x,\ T_{y}=\left(\sigma_{x y},\sigma_{y y}\right)\cdot \mathbf{n}=t_y \left(\text{ on } \Gamma_2\right)\)

其中，$u_0$ 、$v_0$为边界上的位移，$T_{x}$、$T_{y}$ 分别为边界力 $\boldsymbol{T}$ 在 $x$ 和 $y$ 方向的分量。$\mathbf{n}$ 为边界上的单位外法向量。

本构关系

在弹塑性问题中，应力与应变之间一般不再存在一一对应的关系，本构方程只能用增量的形式表出，这里采用应变空间表述的弹塑性本构理论。考虑各向同性的理想弹塑性材料，塑性本构关系采用相关联的塑性流动法则。此时的塑性势函数等于屈服函数，塑性本构矩阵是对称矩阵。 在小应变的情况下，应变增量 $\mathrm{d}\boldsymbol{\varepsilon}$ 可以分解为弹性应变增量 $\mathrm{d}\boldsymbol{\varepsilon}^e$ 和塑性应变增量 $\mathrm{d}\boldsymbol{\varepsilon}^p$ 两部分：
\(\mathrm{d}\boldsymbol{\varepsilon}=\mathrm{d}\boldsymbol{\varepsilon}^e+\mathrm{d}\boldsymbol{\varepsilon}^p\) 利用弹性应力应变关系，可得 \(\mathrm{d}\boldsymbol{\sigma}=\mathbf{D}_e \mathrm{d}\boldsymbol{\varepsilon}^e=\mathbf{D}_e\left(\mathrm{d}\boldsymbol{\varepsilon}-\mathrm{d}\boldsymbol{\varepsilon}^p\right)\)

其中， $\mathbf{D}_e$ 是弹性矩阵。

由塑性势理论 \(\mathrm{d}\boldsymbol{\varepsilon}^p=\mathrm{d}\lambda\frac{\partial F}{\partial \sigma}\)

和一致性条件 \(\mathrm{d}F=\frac{\partial F}{\partial \boldsymbol{\sigma}}\mathrm{d}\boldsymbol{\sigma}=0\)

得到应力应变关系为 \(\mathrm{d}\boldsymbol{\sigma}=\mathbf{D}_{ep}\mathrm{d}\boldsymbol{\varepsilon}=\left(\mathbf{D}_e-\frac{1}{K_p}\frac{\partial F}{\partial \boldsymbol{\sigma}}^T\mathbf{D}_e^T\mathbf{D}_e\frac{\partial F}{\partial \boldsymbol{\sigma}}\right)\mathrm{d}\boldsymbol{\varepsilon}\) \(K_p=\frac{\partial F}{\partial \boldsymbol{\sigma}}^T \mathbf{D}_e \frac{\partial F}{\partial \boldsymbol{\sigma}}\)

其中，$\mathrm{d}\lambda$ 为非负比例系数，$\mathbf{D}_{ep}$ 是弹塑性矩阵。

这里考虑理想弹塑性的 Mises 模型，其屈服面函数为 \(F(\boldsymbol{\sigma})=\sigma_e-\sigma_y\)

其中 $\sigma_e$ 为等效应力，$\sigma_e=\left(\frac{3}{2} \boldsymbol{\sigma}^{\prime}: \boldsymbol{\sigma}^{\prime}\right)^{1 / 2}$，$\boldsymbol{\sigma}^{\prime}=\boldsymbol{\sigma}-\frac{1}{3} \operatorname{tr}(\boldsymbol{\sigma}) \mathbf{1}$；$\sigma_y$ 为屈服应力。弹性参数包括弹性模量 $E$ 和泊松比 $\nu$ 。

弱形式

根据虚位移原理，对平衡方程两边分别乘以位移的变分并在求解区域内积分，得到其等效积分形式 \(\int_{\Omega}\left[\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}+f_{x}\right) \delta u+\left(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}+f_{y}\right) \delta v\right] \mathrm{d}\Omega=0\) 对上式进行分部积分，得到弱形式 \(\int_{\Omega}\left(\sigma_{x x} \delta \varepsilon_{x x}+\sigma_{y y} \delta \varepsilon_{y y}+\sigma_{x y} \delta \varepsilon_{x y}\right) \mathrm{d}\Omega \\=\int_{\Omega}\left(f_{x} \delta u+f_{y} \delta v\right) \mathrm{d}\Omega+\int_{\Gamma_2}\left(T_{x} \delta u+T_{y} \delta v\right) \mathrm{d}\Gamma\)

或 \(\int_{\Omega}\boldsymbol{\sigma}\delta \boldsymbol{\varepsilon}\mathrm{d}\Omega =\int_{\Omega}\boldsymbol{f} \delta \boldsymbol{u} \mathrm{d}\Omega+\int_{\Gamma_2}\boldsymbol{T}\delta \boldsymbol{u} \mathrm{d}\Gamma \tag{1}\label{eq1}\)

计算步骤

在弹塑性问题中，材料的性质与应力和变形的历史有关，本构方程用增量形式表达。这就需要按载荷作用的实际情况，在小的载荷增量下逐步地计算求解。这里将加载过程划分为若干个增量步。对于每一个增量步，继续引入迭代过程， \(\int_{\Omega}\boldsymbol{\sigma}^{n+1}\delta \boldsymbol{\varepsilon}\mathrm{d}\Omega =\int_{\Omega}\boldsymbol{f}^{n+1} \delta \boldsymbol{u} \mathrm{d}\Omega+\int_{\Gamma_2}\boldsymbol{T}^{n+1}\delta \boldsymbol{u} \mathrm{d}\Gamma\) 其中上标 $n$ 代表迭代步数。将 $\boldsymbol{\sigma}^{n+1}=\boldsymbol{\sigma}^{n}+\mathrm{\Delta}\boldsymbol{\sigma}^{n+1}$ 和本构方程 $\mathrm{\Delta}\boldsymbol{\sigma}^{n+1}=\mathbf{D}_{ep}\mathrm{\Delta}\boldsymbol{\varepsilon}^{n+1}$ 代入上式，得到如下增量格式 \(\int_{\Omega}\mathbf{D}_{ep}\mathrm{\Delta}\boldsymbol{\varepsilon}^{n+1}\delta \boldsymbol{\varepsilon}\mathrm{d}\Omega =\int_{\Omega}\boldsymbol{f}^{n+1} \delta \boldsymbol{u} \mathrm{d}\Omega+\int_{\Gamma_2}\boldsymbol{T}^{n+1}\delta \boldsymbol{u} \mathrm{d}\Gamma-\int_{\Omega}\boldsymbol{\sigma}^{n}\delta \boldsymbol{\varepsilon}\mathrm{d}\Omega \tag{2}\label{eq2}\) 选用标准的有限单元对控制方程 $\eqref{eq2}$ 进行空间离散化，得到以位移增量 $\Delta U$ 为基本未知量的代数方程组 \(K_{ep}\Delta U^{n+1}=\Delta F^{n+1}\) 求解该方程组，即可得到当前迭代步的位移增量。重复上述过程，直至位移增量小于某一阈值时，确认收敛，停止迭代。

可以看出，整个求解过程的关键是开发用于本构积分的子程序，实现 $\mathbf{D}_{ep}$ 和 $\boldsymbol{\sigma}^{n}$ 的计算。这会涉及到在循环迭代时对应力和内变量等状态变量的计算、存储和调用。这个过程在每个高斯点上都要执行，是一个相对较为复杂的问题。

脚本编写

一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。下面依次给出详细介绍。

PDE 和 FBC 文件

位移场

对于式 $\eqref{eq2}$，以位移作为基本未知量，基于有限元语言语法，将式 $\eqref{eq2}$ 的体积分部分，编写成了disp.pde 文件，具体内容如下。

由于弹塑性本构方程和积分算法是千变万化的，这里其实仅通过 PDE 文件给出了单元计算程序的整体结构。脚本中使用的 plasmod 模块和 initializer、updater 两个子程序均来自于独立的 plastic.f90 源文件。它专门用于更新高斯点上的应力状态和与弹塑性模型相关的状态变量，需要依靠用户自己来开发。

这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。

对于单元载荷向量中的边界积分部分，用 FBC 文件来描述（FBC 的命名习惯上常取对应 PDE 的文件名，disp.fbc）。将方程 $\eqref{eq1}$ 的边界积分项写入 FBC 文件如下。

注意到这里的单元载荷采用了线性加载的方式，随时间的增长而线性增加。

应力场

应力场只有体积分项而没有边界积分项，所以只需填写 PDE 文件。对应式 $\eqref{eq2}$，利用有限元语言语法编写的 PDE 文件如下（命名为 sdisp.pde）

这里用到的 getstate 子程序和 PDE 文件中的 initializer、updater 一样，同样来自于独立的 plastic.f90 源文件，用于获取应力和塑性应变的值。

注意，后面会看到，由于应力场计算直接使用了位移场的数据作为输入数据，所以，这里应力场 PDE 文件中 mate 行材料参数的定义必须和位移场 PDE 文件 mate 行的定义完全相同。

MDI 文件

MDI 文件负责组织不同场的 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 plas.mdi，具体内容如下。

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。由于应力场计算使用了位移场的计算结果，所以单元类型和积分方法需要和位移场保持一致。

对于弹塑性力学问题，为了保证应力场的计算精度，这里配合 3 点高斯积分，采用了 9 节点四边形单元。

MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 和 FBC 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

在这个例子中，有位移场和应力场两个场。除了对应的 PDE 和 FBC 文件，还要给出每个场使用的算法文件和计算的命令流程，这些信息由 GCN 文件给出。系统要求 GCN 与 MDI 文件同名，具体内容如下。

注意：与 str 算法相对应的求解器为 stress，这是为其专门设置的。stress 求解器不必初始化（即不使用 #start*），直接使用已知耦合场的输入数据即可。因此， GiD 前处理模板中也没有为其添加菜单选项。

NFE 文件

位移场

位移场使用了 nell 算法，系统会基于模板文件自动生成以下的 plasa.nfe 文件。根据弹塑性问题的需要，我们在原有文件的基础上作了一定的修改。

应力场

通过已知位移求应力，应力场使用了 str 算法，系统会基于模板文件自动生成以下的 plasb.nfe 文件。根据弹塑性问题的需要，我们在原有文件的基础上作了一定的修改，去除了关于耦合变量读取的内容。

算例

一个内、外半径分别为 $100\ \mathrm{mm}$ 和 $200\ \mathrm{mm}$ 的厚壁圆筒，受内压 $P$ 的作用，不计体力。设圆筒的长度远大于筒的直径，材料是理想弹塑性的，满足 Mises 屈服条件，具体参数值为 $E＝2.1×10^{11}\ \mathrm{Pa}$，$\nu = 0.3$ ，屈服强度 $\sigma_y = 2.4 ×10^{8}\ \mathrm{Pa}$ 。圆筒的理论屈服极限为 $P ＝ 1.9209×10^{8}\ \mathrm{Pa}$。为了模拟屈服过程，采用线性加载至 $P ＝ 1.92×10^{8}\ \mathrm{Pa}$，求该圆筒内部的应力场和位移场。

前处理

根据圆截面的对称性，取四分之一的区域进行模拟。

填写材料参数文件

系统自动生成了 plas.mat 文件，用于设置材料参数。参数初始值全部为 PDE 和 FBC 文件给出的默认值。可根据具体问题进行修改。对于本算例，修改完成后的文件内容如下。由于应力场计算直接使用了位移场的数据作为输入数据，所以没有相关的单元和参数选项。

填写时间步长文件

由 plas.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出虚拟的时间步长和终止时间，用于实现线性加载。文件内容如下：

0.2 19.2

可知，时间步长为 0.2 s，计算终止时间为 19.2 s。

赋单元材料号

材料以面单元的形式赋值

施加固定边界

在两个侧边施加法向约束

其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点力

施加荷载

内压荷载以线单元的形式施加

网格剖分

根据所开发程序的计算要求，采用四边形二次单元进行网格剖分。共使用了 1600 个单元，6561 个节点。

计算结果

对计算结果稍加整理，效果如下。

位移分布

以下为轴向位移分布的演化过程。

下图给出了随着荷载的线性施加，圆筒外壁轴向位移的变化情况，可以看出与理论值完全相符。

应力分布

以下为应力分布的演化过程，依次表示 $\sigma_{x x}$ 、$\sigma_{y y}$ 和 $\sigma_{x y}$ 。

下图给出了当 $P ＝ 1.8×10^{8}\ \mathrm{Pa}$ 时，环向和轴向应力的分布图。

可以发现数值解与解析解符合得很好，进而证明了算法和程序的有效性。

等效塑性应变

以下为等效塑性应变和塑性区的演化情况。可以清楚地看到随着荷载的线性施加，圆筒逐渐屈服的全过程。

参考文献

de Souza Neto E A, Peric D, Owen D R J. Computational methods for plasticity: theory and applications[M]. John Wiley & Sons, 2011.