弹性动力学平面应力问题

关键词： 瞬态 线性 弹性力学 平面应力 惯性力 Newmark法

本节以二维弹性动力学平面应力问题为例，讲解如何使用 FEtch 系统对线性二阶双曲型方程进行求解。

控制方程

微分形式

对于二维直角坐标系，弹性动力学问题的未知量服从如下偏微分方程

平衡方程 \(\rho\frac{\partial^2 u}{\partial t^2}+\mu\frac{\partial u}{\partial t}-\frac{\partial \sigma_{x x}}{\partial x}-\frac{\partial \sigma_{x y}}{\partial y}-f_{x} = 0 \ \left(\text{in } \Omega\right)\) \(\rho\frac{\partial^2 v}{\partial t^2}+\mu\frac{\partial v}{\partial t}-\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$ 为直角坐标系下的坐标分量，$t$ 为时间。

初始条件

初始位移 \(u=u_0,\ v=v_0 \ \left(\text{in } \Omega\right)\) 初始速度 \(\dot{u}=\dot{u_0},\ \dot{v}=\dot{v_0} \ \left(\text{in } \Omega\right)\)

初始加速度 \(\ddot{u}=\ddot{u_0},\ \ddot{v}=\ddot{v_0}\ \left(\text{in } \Omega\right)\)

这里 $\dot{}$ 和 $\ddot{}$ 表示对时间的一阶和二阶导数。

边界条件

考虑两类边界条件:

第一类边界条件，给定位移 \(u=\bar{u},\ v=\bar{v}\ \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)\)

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

弱形式

PDE 文件的编写依赖于待解偏微分方程的弱解方程。根据虚位移原理，对平衡方程两边分别乘以位移的变分并在求解区域内积分，得到其等效积分形式 \(\int_{\Omega}\left(\rho\frac{\partial^2 u}{\partial t^2}+\mu\frac{\partial u}{\partial t}-\frac{\partial \sigma_{x x}}{\partial x}-\frac{\partial \sigma_{x y}}{\partial y}+f_{x}\right) \delta ud\Omega\\+\int_{\Omega}\left(\rho\frac{\partial^2 v}{\partial t^2}+\mu\frac{\partial v}{\partial t}-\frac{\partial \sigma_{x y}}{\partial x}-\frac{\partial \sigma_{y y}}{\partial y}+f_{y}\right) \delta v d\Omega=0\) 对上式进行分部积分，得到弱形式 \(\int_{\Omega}\rho\left(\frac{\partial^2 u}{\partial t^2}\delta u+\frac{\partial^2 v}{\partial t^2}\delta v\right)d\Omega+\int_{\Omega}\mu\left(\frac{\partial u}{\partial t}\delta u+\frac{\partial v}{\partial t}\delta v\right)d\Omega\\+\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}\)

将本构方程代入上式，即可得到以位移为基本未知量的弱形式。于是，相关有限元模型的质量矩阵、阻尼矩阵和矩阵刚度矩阵分别由上式左端的体积分部分的 1、2、3 项形成，而载荷向量则包含上式右端的体积分和边界积分两部分的贡献。

在求得了位移之后，可以应用最小二乘法，继续求取应力 $\boldsymbol{\sigma}$ 。这里主要介绍位移场的求解，对应力场求解感有兴趣的用户可以参考弹性静力学章节的有关介绍。

脚本编写

求解策略

对于式 $\eqref{eq1}$，以位移作为基本未知量，位移场的有限元模型可以写成以下矩阵形式 \(M\ddot{U}+C\dot{U}+KU=F\) 其中，$M$ 为系统的质量矩阵，$C$ 为阻尼矩阵，$K$ 为刚度矩阵，$F$ 为载荷向量，$U$、$\dot{U}$ 和 $\ddot{U}$ 分别为为节点位移、速度和加速度向量。该式即为需要求解的半离散的线性方程组。对于时域积分，有多种方法可以选择，这里使用最常用的 Newmark 法。

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

PDE 和 FBC 文件

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

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

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

MDI 文件

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

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。

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

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

GCN 文件

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

NFE 文件

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

算例1

如图所示，边长为 $1.0\ \mathrm{m}\times0.1\ \mathrm{m}$ 的梁结构，左侧的位移完全约束，右侧和两侧边可自由移动，右侧施加均布的法向压力荷载 $F_N$，不计材料重力影响，求梁的位移的演化过程。具体参数值为 $E＝1.0×10^{3}\ \mathrm{Pa}$，$\nu = 0.3$ ，$f_{x} ＝ 0$， $f_{y} ＝ 0$，$F_N ＝ 10\sin(2\pi t)\ \mathrm{N/m}$。

前处理

填写材料参数文件

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

填写时间步长文件

由 dyn.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出时间步长和终止时间。文件内容如下：

0.02 2

可知，时间步长为 0.02 s，计算终止时间为 2 s。

赋单元材料号

材料以面单元的形式赋值

赋初始条件

初始位移

初始速度

初始加速度

施加固定边界

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

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

施加荷载

梁右侧的荷载以线单元的形式施加

网格剖分

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

计算结果

位移分布

以下为位移分布图，依次表示 $u$ 和 $v$。

速度分布

以下为速度分布图，依次表示 $u$ 和 $v$。

算例2

模型和具体参数值同算例 1。右侧施加分布荷载改为切向力 $F_T＝ 10\sin(2\pi t)\ \mathrm{N/m}$。求梁的位移的演化过程。

前处理

填写材料参数文件

因材料与算例 1 相同，只需修改与荷载相关的线形边界单元的参数值。修改后的 dyn.mat 文件内容如下。

其他处理

时间步长文件不变，赋单元材料号、赋初始条件、施加固定边界、施加荷载、网格剖分的操作过程均与算例1完全一样，不再赘述。

计算结果

位移分布

以下为位移分布图，依次表示 $u$ 和 $v$。

速度分布

以下为速度分布图，依次表示 $u$ 和 $v$。