含非均匀体力机械结构弹性力学问题的双互易边界元法

曾 华， 周枫林， 余江鸿

（湖南工业大学 机械工程学院， 湖南 株洲412007）

0 引言

在工程实践中往往会遇到一些弹性力学方面的问题，而一般的弹性力学问题用解析方法求解就会非常麻烦，此时，采用数值方法来解决这类问题更为有效[1]。 到目前为止，人们所熟悉的工程数值方法主要有：有限单元法[2]、Ritz 法[3]、Galerkin 法[4]以及边界元法[5]。 近年来，随着计算机技术飞速发展，边界元法得到了广泛应用，边界元法最重要的特点是降低维度，从而提高计算效率。由于边界元法具有较高的计算精度， 在处理无限域问题或应力集中问题的精度比有限元分析要好， 所以这种方法受到了国内外许多研究者的重视。 边界元法作为一种常用的数值分析方法，在弹性力学、流体力学、热传导、声场，以及电磁场等研究领域得到了广泛应用[6]。 在边界元法中，由于边界积分方程对于含体力项的弹性问题求解起来十分困难，需要采用一种更有效的方法来解决体力项的问题，而双重互易法（DRM）的出现就很好地解决体力项的问题，近年来双重互易法被许多学者所研究，其中文献[7-10]研究较多， 双重互易法实现了将体力项产生的域内积分转化为等效的边界积分， 从而有效避免了域内积分的复杂计算， 本文通过编写边界元程序验证了双重互易法在求解弹性问题的有效性。

1 弹性力学问题的边界元法

描述线弹性问题的平衡方程、几何方程和本构方程为：

式中：i，j—方向；σij—应力；εkl—应变；bi—体积力；Eijkl—弹性模量；uij—位移，uij，j表示uij在j 方向上的导数；V—域内体积。 若要求出方程的解还需要包含一定的边界条件，可以写成：

式中：pi—面力；nj—边界外法线方向的法向量；S—边界面积。

通过Betti 功互等定理（简称Betti 定理）我们可以建立边界积分方程，形如：

式中：c—光滑系数，一般光滑边界点c=1/2， u*和p*分别表示位移基本解和面力基本解， 其值与节点位置、 泊松比、弹性模量有关。

2 双重互易法

边界积分方程中如果在无体力项的情况下， 只要边界条件中位移或者面力中的一项是已知的， 就可以通过边界积分方程求出另外一项， 但是由于在实际工程应用中，物体往往都有体力项存在，所以对于含体力项的弹性力学方程求解问题，原本的边界元方法就不太适用。

为了解决体力项的问题，采用了双重互易法。即当定义中考虑体力项的区域积分时， 可以将体力项积分近似的表示为：

式中：H 和G 矩阵分别由面力和位移基本解组成；F-1表示体积力的逆矩阵。

3 数值算例

为了验证双互易边界元法的有效性，在Visual studio的软件环境下设计两个边界元法分析程序。 第一个是在给定边界位移的条件下计算立方体模型域内的位移变化， 第二个是计算球体模型的域内位移变化。

3.1 立方体模型域内位移变化

建立一个边长为2 的立方体简易模型， 网格划分如图1 所示，该立方体的网格大小为0.1，采用的边界节点数目为792， 所构建的二次三节点单元为1332 个， 内部RBF 插值点个数为1895， 内部插值点位置分布如图2 所示，材料参数：弹性模量设置为1，泊松比0.25，材料密度1.14。 立方体的空间位置范围为：

图1 立方体网格划分

图2 插值点位置分布

其中网格大小与节点数量的关系如表2 所示，m 表示网格大小，N 表示节点数量，L 表示域内插值点数量。选取5 个域内节点作为参考点，见表1，x 方向上的计算结果如表2 所示。

表1 参考点坐标

表2 网格大小对应的节点数量

从表3 中可以看出，当固定y，z 方向的坐标值，随机选择域内五组参考点， 各点在x 方向上的位移计算结果与解析解的吻合度比较好，说明了采用此方法是有效的，为了进一步证明方法的准确性， 验证了在不同网格大小下计算结果与解析解的变化，计算A1 在x 方向上位移相对差值如图3 所示。

从图3 可以看出，随着插值点的不断增加，位移相对误差也会越来越小， 计算结果随着网格大小的增加呈单调递减的趋势，当网格大小为0.1 时，A1 的位移相对误差为3.1393%，由此可以看出，只要选择合适的网格大小，双互易边界元法具有良好的的精度。

表3 x 方向上的位移结果

图3 A1 在x 分向上的位移相对差值

3.2 球体模型域内位移变化

建立一个半径为2 的球体模型， 网格划分如图4 所示， 该球体的网格大小为0.3 采用的边界节点数目为581，所构建的二次三节点单元为1066 个，内部RBF 插值点个数为2003，内部插值点位置分布如图5 所示，材料参数与正方体相同，球体坐标原点为（0，0，0）。 空间坐标为

图4 球体网格划分

图5 插值点位置分布

边界条件：

同样选取5 个参考点， 见表4，y 方向的结果如表5所示，通过表5 中计算结果与解析解进行对比可知，此方法得到的结果与精确解吻合度很好， 两组实验都反映了采用双互易边界元法在求解弹性问题上是有效的。

表4 参考点坐标

表5 y 方向的位移结果

4 结论

运用双互易边界元法避免了复杂的域内积分计算；通过两组数值分析， 验证了双互易边界元法在求解弹性问题上的有效性，并随着网格划分得越来越细，计算结果具有更好的精度； 不足之处在于算法并未采用任何加速计算，因此提高计算效率是进一步需要完成的工作。