弹性力学平面应力问题

关键词： 稳态 线性 弹性力学 平面应力 应力场 耦合

本节以二维弹性力学平面应力问题为例，讲解如何使用 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}\) 本构方程（平面应力） \(\left(\begin{array}{c}\sigma_{x x} \\ \sigma_{y y} \\ \sigma_{x y} \end{array}\right)=\frac{E}{(1+\nu)(1-\nu)}\left(\begin{array}{ccc} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & (1-\nu) / 2 \end{array}\right)\left(\begin{array}{l} \varepsilon_{x x} \\ \varepsilon_{y y} \\ \varepsilon_{x y} \end{array}\right)\) 以上式子可以简记为 \(\nabla\boldsymbol{\cdot \sigma}+\boldsymbol{f}=\boldsymbol{0}\) \(\boldsymbol{\varepsilon}=\left(\nabla \boldsymbol{u}+\nabla^{T} \boldsymbol{u}\right) / 2\) \(\boldsymbol{\sigma}=\mathbf{D} \boldsymbol{\varepsilon}\) 其中， $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$ 为梯度算子，$\mathbf{D}$ 为刚度矩阵。参数 $E$ 为弹性模量， $\nu$ 为泊松比。 $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}$ 分别为边界力在 $x$ 和 $y$ 方向的分量。$\mathbf{n}$ 为边界上的单位外法向量。

弱形式

PDE 文件的编写依赖于待解偏微分方程的弱解方程。根据虚位移原理，对平衡方程两边分别乘以位移的变分并在求解区域内积分，得到其等效积分形式 \(\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] 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) d\Omega \\=\int_{\Omega}\left(f_{x} \delta u+f_{y} \delta v\right) d\Omega+\int_{\Gamma_2}\left(T_{x} \delta u+T_{y} \delta v\right) d\Gamma \tag{1}\label{eq1}\)

将本构方程代入上式，即可得到以位移为基本未知量的弱形式。于是，此有限元问题的单元刚度矩阵由上式左边体积分部分形成，而单元载荷则包含上式右端体积分和边界积分两部分的贡献。

在求得了位移之后，利用本构方程和几何方程可以得到应力与位移的关系 \(\boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{D} \boldsymbol{\varepsilon}(\boldsymbol{u})\) 应用最小二乘法，以待求应力 $\boldsymbol{\sigma}$ 与由位移求得的应力 $\boldsymbol{\sigma}_0(\boldsymbol{u})$ 的偏差为目标泛函，即 \(F(\boldsymbol{\sigma})=\int_{\Omega}[\boldsymbol{\sigma-\sigma_0}(\boldsymbol{u})]^{2} d \Omega\) 对 $F(\boldsymbol{\sigma})$ 取极值，由变分原理知需满足 $\delta F(\sigma) = 0$。于是有 \(\delta F(\boldsymbol{\sigma})=\int_{\Omega} 2 \delta \sigma^{\mathrm{T}}[\boldsymbol{\sigma}-\boldsymbol{\sigma}_0(\boldsymbol{u})] \mathrm{d} \Omega=0\) 整理得 \(\int_{\Omega} \boldsymbol{\sigma} \delta \boldsymbol{\sigma} d\Omega=\int_{\Omega} \boldsymbol{\sigma_{0}}(\boldsymbol{u}) \delta \boldsymbol{\sigma} d\Omega \tag{2}\label{eq2}\)

脚本编写

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

PDE 和 FBC 文件

位移场

对于式 $\eqref{eq1}$，以位移作为基本未知量，位移场的控制方程可以写成以下矩阵形式 \(KU=F \tag{3}\label{eq3}\) 其中，$K$ 为刚度矩阵，$U$ 为位移向量，$F$ 为载荷向量。该式即为最终需要求解的线性方程组。

基于有限元语言语法，将式 $\eqref{eq1}$ 的体积分部分，编写成了disp.pde 文件，具体内容如下。

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

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

应力场

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

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

MDI 文件

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

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

对于弹性力学问题，理论上四边形单元的计算精度要明显好于三角形单元。因此，这里配合 2 点高斯积分，采用了四边形单元。

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

GCN 文件

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

注意：stress 求解器是一种特殊化了的求解器，专门用于求解已知场的梯度场。它与 str.nfe 算法模板文件互相匹配。stress 求解器不必初始化（即不使用 #start*），而是直接使用与之耦合的 a 场的网格和材料数据作为输入数据。因此，后面我们会看到，GiD 前处理模板中没有为 b 场添加菜单选项。

NFE 文件

位移场

位移场使用了 ell 算法，系统会基于模板文件自动生成以下的 elasa.nfe 文件。如果用户有其他需要，可以在此基础上进一步地修改。

应力场

通过已知位移求应力，应力场使用了 str 算法，系统会基于模板文件自动生成以下的 elasb.nfe 文件。

注意，NFE 文件的 coef 行，要保证与 PDE文件中 coef 行耦合场变量的数目和顺序保持一致，否则会出现信息传递错误。

算例1

边长为 $0.5\ \mathrm{m}\times1.0\ \mathrm{m}$ 的长方形薄板，底边的位移竖向固定，横向自由，两侧边自由，顶边施加分布荷载 $P$，不计材料重力影响，求薄板的变形和应力分布。具体参数值为 $E＝1.0×10^{10}\ \mathrm{Pa}$，$\nu = 0.3$ ，$f_{x} ＝ 0$， $f_{y} ＝ 0$，$P ＝ 100\ \mathrm{N/m}$（压力）。

前处理

填写材料参数文件

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

赋单元材料号

材料以面单元的形式赋值

施加固定边界

底部左侧端点 $u$ 和 $v$ 全部约束

底部边界 $u$ 自由，$v$ 约束

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

施加荷载

物体上部的荷载以线单元的形式施加

网格剖分

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

计算结果

位移分布

以下为位移分布图，依次表示 $u$ 和 $v$。均为线性变化，与理论值完全相符。

应力分布

以下为应力分布图，依次表示 $\sigma_{x x}$ 、$\sigma_{y y}$ 和 $\sigma_{x y}$ 。$\sigma_{y y}=-100\ \mathrm{Pa}$，呈均匀分布，和理论解完全相同。$\sigma_{x x}$ 和 $\sigma_{x y}$ 理论值为 0，而计算结果为 $10^{-11}$ 数量级，逼近于 0，误差极小。

算例2

正方形薄板，尺寸为 $2\ \mathrm{m}\times2\ \mathrm{m}$，中间有一小圆孔，半径为 $0.2\ \mathrm{m}$。在薄板四周施加分布荷载 $P$，求薄板的变形和应力分布。具体参数值同算例 1。

前处理

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

填写材料参数文件

因材料和荷载与算例 1 相同，所以仍使用算例 1 的 elas.mat 文件。

赋单元材料号

材料以面单元的形式赋值

施加固定边界

底部边界 $u$ 自由，$v$ 约束；左侧边界 $v$ 自由，$u$ 约束。

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

施加荷载

物体上部的荷载以线单元的形式施加

网格剖分

根据所开发程序的计算要求，采用四边形线性单元进行网格剖分。共使用了 441 个单元，400 个节点。为提高计算精度，按 1.1:1 的比例，沿径向对小孔周围的网格进行了局部加密。

计算结果

位移分布

以下为位移分布图，依次表示 $u$ 和 $v$，以及径向位移 $r=\sqrt{u^2+v^2}$。结果与理论值完全相符。

应力分布

以下为应力分布图，依次表示 $\sigma_{x x}$ 、$\sigma_{y y}$ 和 $\sigma_{x y}$ 。其中法向应力的最大值为 $201\ \mathrm{Pa}$，和理论解 $2P=200\ \mathrm{Pa}$ 非常接近。