弹塑性力学问题的有限元求解

学术科技 2024-09-24 09:28 山东

点击上方「蓝字」关注我们

有限元方法是研究材料力学行为的一个有效工具。本文以二维理想弹塑性力学平面应变问题为例，介绍有限元方法和 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)

几何方程

\begin{aligned} \varepsilon_{x x}&=\frac{\partial u}{\partial x}, \\ \varepsilon_{y y}&=\frac{\partial v}{\partial y}, \\ \varepsilon_{x y}&=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \end{aligned}

以上式子可以简记为

\nabla\boldsymbol{\cdot \sigma}+\boldsymbol{f}=\boldsymbol{0}

\boldsymbol{\varepsilon}=\left(\nabla \boldsymbol{u}+\nabla^{T} \boldsymbol{u}\right) / 2

其中，、为位移的分量，、、为应变的分量，、、为应力的分量，、为体力的分量。为梯度算子，、为直角坐标系下的坐标分量。

边界条件

考虑两类边界条件：

第一类边界条件，给定位移

u=u_0,\ v=v_0 \quad \left(\text{ on } \Gamma_1\right)

第二类边界条件，给定外力

\begin{aligned} 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 \quad \left(\text{ on } \Gamma_2\right) \end{aligned}

其中，、为边界上的位移，、分别为边界力在和方向的分量。为边界上的单位外法向量。

本构关系

在弹塑性问题中，应力与应变之间一般不再存在一一对应的关系，本构方程只能用增量的形式表出，这里采用应变空间表述的弹塑性本构理论。考虑各向同性的理想弹塑性材料，塑性本构关系采用相关联的塑性流动法则。此时的塑性势函数等于屈服函数，塑性本构矩阵是对称矩阵。在小应变的情况下，应变增量可以分解为弹性应变增量和塑性应变增量两部分：

\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)

其中，是弹性矩阵。

由塑性势理论

\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

得到应力应变关系为

\begin{aligned} \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} \end{aligned}

K_p=\frac{\partial F}{\partial \boldsymbol{\sigma}}^T \mathbf{D}_e \frac{\partial F}{\partial \boldsymbol{\sigma}}

其中，为非负比例系数，是弹塑性矩阵。

这里考虑理想弹塑性的 Mises 模型，其屈服面函数为

F(\boldsymbol{\sigma})=\sigma_e-\sigma_y

其中为等效应力，，；为屈服应力。弹性参数包括弹性模量和泊松比。

弱形式

根据虚位移原理，对平衡方程两边分别乘以位移的变分并在求解区域内积分，得到其等效积分形式

\begin{aligned} &\int_{\Omega}\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}+f_{x}\right) \delta u\mathrm{d}\Omega\\+&\int_{\Omega}\left(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}+f_{y}\right) \delta v \mathrm{d}\Omega=0 \end{aligned}

对上式进行分部积分，得到弱形式

\begin{aligned} &\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 \end{aligned}

或

\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

计算步骤

在弹塑性问题中，材料的性质与应力和变形的历史有关，本构方程用增量形式表达。这就需要按载荷作用的实际情况，在小的载荷增量下逐步地计算求解。这里将加载过程划分为若干个增量步。对于每一个增量步，继续引入迭代过程，

\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

其中，上标代表迭代步数。将和本构方程代入上式，得到如下增量格式

\begin{aligned} &\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 \end{aligned}

选用标准的有限单元对上式进行空间离散化，得到以位移增量为基本未知量的代数方程组

K_{ep}\Delta U^{n+1}=\Delta F^{n+1}

求解该方程组，即可得到当前迭代步的位移增量。重复上述过程，直至位移增量小于某一阈值时，确认收敛，停止迭代。

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

算例

一个内、外半径分别为和的厚壁圆筒，受内压的作用，不计体力。设圆筒的长度远大于筒的直径，材料是理想弹塑性的，满足 Mises 屈服条件，具体参数值为，，屈服强度。圆筒的理论屈服极限为。为了模拟屈服过程，采用线性加载至，求该圆筒内部的应力场和位移场。

前处理

前处理采用 GiD 软件。根据圆截面的对称性，取四分之一的区域进行模拟。在两个侧边施加法向约束，内压荷载以线单元的形式施加。根据所开发程序的计算要求，采用四边形二次单元进行网格剖分。共使用了 1600 个单元，6561 个节点。时间步长取 0.2 s，计算终止时间为 19.2 s。

计算结果

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

位移分布

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

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

应力分布

以下为应力分布的演化过程，依次表示、和。

下图给出了当时，环向和轴向应力的分布图。

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

等效塑性应变

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

参考文献

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

往期推荐

一图读懂 FEtch 的工作流程

谈谈有限元程序开发，兼谈 FEtch 入门

通用前后处理软件 GiD 简介