CN105574281A - 一种自适应二维有限元疏密网格界面过渡方法 - Google Patents

Abstract

本发明公开了一种自适应二维有限元疏密网格界面过渡方法，应用于工程结构分析领域。该方法通过惩罚因子以最小二乘形式强迫不同疏密网格子区域间的有限元节点位移连续，根据最小势能原理，获得界面过渡元和系统的刚度矩阵。将刚度矩阵通过用户单元子程序接口，引入到商用有限元软件中。用户只需设定单元的材料属性、几何特性和载荷等性质来选取合适的惩罚因子，即可实现二维有限元网格界面过渡。本发明方法能有效地解决有限元网格划分时局部细化问题和不同地域、不同分析者采用不同有限元软件建立的有限元模型之间远程对接问题，具有自适应性，简单易行，计算精度高，占用内存小，计算速度快。

Description

一种自适应二维有限元疏密网格界面过渡方法

技术领域

本发明涉及工程结构分析领域，具体涉及一种自适应二维有限元疏密网格界面过渡方法。

背景技术

作为一种对连续弹性体进行离散化的近似数值解法，有限元法因其能够模拟几何形状复杂的结构，计算复杂结构的场问题，得出近似解；同时因其解题步骤已经系统化、标准化，有相当数量灵活通用的计算机程序，可以解决各种不同的工程问题，成为目前工程计算中应用最广泛最有效的数值计算方法之一。

经典有限元网格离散单元的基本思想是通过相互协调的节点传递位移和力，要求划分网格时单元之间的节点必须一致，否则会引起较大的计算误差。然而，在工程结构的有限元分析中，常常存在两种情况：一是构件内部存在稀疏和致密两种网格。如航空发动机机匣、火焰筒等为多孔结构，有限元分析时为了保证一定的计算精度，同时提高计算效率，建模时孔周围要有一套比较致密的网格，其他部分则采用较稀疏的网格。二是随着模型共享和大规模分析日益增加，现代复杂问题多由不同的分析者或组织在不同的地理位置应用不用软件建立不同构件的模型，如航空发动机的涡轮盘和涡轮叶片。这些模型在彼此界面处很可能节点不协调，需要花费大量的精力和时间来实现两种网格的过渡。因此，研究有限元网格过渡方法对有限元法的进一步应用有着重要的意义。

国内外对有限元网格过渡方法开展了一定的研究，提出了主-从边界法、特殊单元法、位移函数法等。这些方法的基本思想均是调整常规有限元方程中的刚度项，在刚度里体现位移不协调而产生的附加能量，然后进行整体求解来解决网格间的不协调过渡，但由于这些方法需要花费大量的时间和精力来选取界面的虚拟节点数，且不同的方法选取的虚拟节点数是不同的，当虚拟节点数选择不合适时，造成刚度矩阵奇异。另外，计算结果的精度要受到虚拟节点处插值函数的限制，给计算带来了很大的不便。目前，国内外还没有较为高效的有限元网格过渡方法，也尚未见公开的发明专利。因此，提出一种降低有限元计算时间，同时又能保证足够的计算精度的自适应二维有限元疏密网格界面过渡方法非常有必要。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种自适应二维有限元疏密网格界面过渡方法，通过划分区域的计算，实现了不同模型之间的拼接，解决了现有技术的不足。

技术方案：为实现上述目的，本发明采用的技术方案为：一种自适应二维有限元疏密网格界面过渡方法，其特征在于，该方法包括以下步骤：

1)建立带垂直界面的二维平板的几何模型，将该模型划分为两个区域，并对两个区域进行划分网格，网格相交处为节点；该二维平板的一端固定，另一端受均布拉伸载荷f；固定端的节点为固定节点，其余节点为自由节点；

2)以惩罚因子γ表示系统的总势能；

3)采用传统有限元方法确定不含虚拟节点的罚函数的有限元方程；

4)计算界面单元刚度矩阵；

5)计算系统刚度矩阵；

6)确定惩罚因子γ；

7)采用常规方法计算出位移和应力应变。

进一步的，所述步骤1)中，所述两个区域分别为区域一、区域二；采用四边形网格进行划分，每个网格为一个单元；区域一和区域二之间为界面，区域一的另一侧为固定端。

进一步的，步骤2)中，系统的总势能为：

$\mathrm{\π}={\mathrm{\π}}_{{\mathrm{\Ω}}_1}+{\mathrm{\π}}_{{\mathrm{\Ω}}_2}+\frac{1}{2}\mathrm{\γ}\underset{S}{\∫}{(u_1-u_2)}^{2}ds$

π对各自由度求一阶变分，并置为零，可得：

$\mathrm{\δ}\mathrm{\π}{|}_{q_1^i,q_1^o,q_2^i,q_2^o}=0$

$\mathrm{\δ}\mathrm{\π}={\mathrm{\δ\π}}_{{\mathrm{\Ω}}_1}+{\mathrm{\δ\π}}_{{\mathrm{\Ω}}_2}+\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_2\right)ds+\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_1\right)ds-\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_1\right)ds-\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_2\right)ds;$

其中，区域一的积分为：

$\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_1\right)ds=\mathrm{\γ}{\∫}_0^{2/3}\left({\mathrm{\δu}}_1{u}_1\right)ds+\mathrm{\γ}{\∫}_{2/3}^{4/3}\left({\mathrm{\δu}}_1{u}_1\right)ds+\mathrm{\γ}{\∫}_{4/3}^2\left({\mathrm{\δu}}_1{u}_1\right)ds=\mathrm{\γ}{\∫}_0^{2/3}({\mathrm{\δq}}_{11}{N}_1+{\mathrm{\δq}}_{12}{N}_2)(q_{11}{N}_1+q_{12}{N}_2)ds+$

$\mathrm{\γ}{\∫}_{2/3}^{4/3}({\mathrm{\δq}}_{12}{N}_1+{\mathrm{\δq}}_{13}{N}_2)(q_{12}{N}_1+q_{13}{N}_2)ds+\mathrm{\γ}{\∫}_{4/3}^2({\mathrm{\δq}}_{13}{N}_1+{\mathrm{\δq}}_{14}{N}_2)(q_{13}{N}_1+q_{14}{N}_2)ds;$

区域二的积分为：

$\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_2\right)ds=\mathrm{\γ}{\∫}_0^1\left({\mathrm{\δu}}_2{u}_2\right)ds+\mathrm{\γ}{\∫}_1^2\left({\mathrm{\δu}}_2{u}_2\right)ds=\mathrm{\γ}{\∫}_0^1({\mathrm{\δq}}_{21}{N}_1+{\mathrm{\δq}}_{22}{N}_2)(q_{21}{N}_1+q_{22}{N}_2)ds+$ $\mathrm{\γ}{\∫}_1^2({\mathrm{\δq}}_{22}{N}_1+{\mathrm{\δq}}_{23}{N}_2)(q_{22}{N}_1+q_{23}{N}_2)ds;$

其中：

δ为一阶变分标志；

和π_Ω2分别为区域一和区域二的势能；

μ₁和μ₂分别为区域一和区域二相邻界面处的位移；

N₁和N₂分别为区域一和区域二的位移插值函数；

区域一积分中的自由q₁₁、q₁₂、q₁₃、q₁₄分别是区域一与区域二相邻界面处区域一上的节点的轴向位移；

区域二积分中的自由q₂₁、q₂₂、q₂₃分别是区域二与区域一相邻界面处区域二上的节点的轴向位移；

进一步的，步骤3)中，不含虚拟节点的罚函数的有限元方程为：

$\left[\begin{array}{cccc}{K}_1^{oo}& K_1^{oi}& 0& 0\\ K_1^{io}& K_1^{ii}+G_1^{ii}& -G& 0\\ 0& -G^T& G_2^{ii}+K_2^{ii}& K_2^{io}\\ 0& 0& K_2^{oi}& K_2^{oo}\end{array}\right]\left\{\begin{array}{c}{q}_1^0\\ q_1^i\\ q_2^i\\ q_2^0\end{array}\right\}=\left[\begin{array}{c}{f}_1^i\\ f_1^o\\ f_2^i\\ f_2^o\end{array}\right]$

其中：

K₁为区域一所有单元的刚度矩阵；K₂为区域二所有单元的刚度矩阵；

上标i表示界面上的单元的自由度；上标o表示不在界面上的单元的自由度；

G由下列积分来计算：

$G_1^{ii}=\mathrm{\γ}{\∫}_s{N}_1^T{N}_{1}ds$

$G_2^{ii}=\mathrm{\γ}{\∫}_s{N}_2^T{N}_{2}ds$

$G=\mathrm{\γ}{\∫}_s{N}_1^T{N}_{2}ds.$

进一步的，步骤4)中，区域一中的固定节点自由度等于0；界面单元的刚度矩阵为：

$\left[\begin{array}{ccccccc}\frac{2\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& 0& -\frac{7\mathrm{\γ}}{27}& -\frac{2\mathrm{\γ}}{27}& 0\\ \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& -\frac{25\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{\mathrm{\γ}}{108}\\ 0& \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& -\frac{\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{25\mathrm{\γ}}{108}\\ 0& 0& \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& 0& -\frac{2\mathrm{\γ}}{27}& -\frac{7\mathrm{\γ}}{27}\\ -\frac{7\mathrm{\γ}}{27}& -\frac{25\mathrm{\γ}}{108}& -\frac{\mathrm{\γ}}{108}& 0& \frac{\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& 0\\ -\frac{2\mathrm{\γ}}{27}& -\frac{23\mathrm{\γ}}{54}& -\frac{23\mathrm{\γ}}{54}& -\frac{2\mathrm{\γ}}{27}& \frac{\mathrm{\γ}}{6}& \frac{2\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}\\ 0& -\frac{\mathrm{\γ}}{108}& -\frac{25\mathrm{\γ}}{108}& -\frac{7\mathrm{\γ}}{27}& 0& \frac{\mathrm{\γ}}{6}& \frac{\mathrm{\γ}}{3}\end{array}\right].$

进一步的，步骤5)中，系统整体刚度矩阵为：

$K=\left[\begin{array}{cccccccccc}\frac{1}{2}+\frac{2\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& 0& -\frac{7\mathrm{\γ}}{27}& -\frac{2\mathrm{\γ}}{27}& 0& 0& 0& 0\\ \frac{\mathrm{\γ}}{9}& 1+\frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& -\frac{25\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{\mathrm{\γ}}{108}& 0& 0& 0\\ 0& \frac{\mathrm{\γ}}{9}& 1+\frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& -\frac{\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{25\mathrm{\γ}}{108}& 0& 0& 0\\ 0& 0& \frac{\mathrm{\γ}}{9}& \frac{1}{2}+\frac{2\mathrm{\γ}}{9}& 0& -\frac{2\mathrm{\γ}}{27}& -\frac{7\mathrm{\γ}}{27}& 0& 0& 0\\ -\frac{7\mathrm{\γ}}{27}& -\frac{25\mathrm{\γ}}{108}& -\frac{\mathrm{\γ}}{108}& 0& \frac{1}{2}+\frac{\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& 0& -\frac{1}{4}& -\frac{1}{4}& 0\\ -\frac{2\mathrm{\γ}}{27}& -\frac{23\mathrm{\γ}}{54}& -\frac{23\mathrm{\γ}}{54}& -\frac{2\mathrm{\γ}}{57}& \frac{\mathrm{\γ}}{6}& 1+\frac{2\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& -\frac{1}{4}& -\frac{1}{2}& -\frac{1}{4}\\ 0& -\frac{\mathrm{\γ}}{108}& -\frac{25\mathrm{\γ}}{108}& -\frac{7\mathrm{\γ}}{27}& 0& \frac{\mathrm{\γ}}{6}& \frac{1}{2}+\frac{\mathrm{\γ}}{3}& 0& -\frac{1}{4}& -\frac{1}{4}\\ 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{4}& 0& \frac{1}{2}& 0& 0\\ 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{2}& -\frac{1}{4}& 0& 1& 0\\ 0& 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{4}& 0& 0& \frac{1}{2}\end{array}\right].$

进一步的，步骤6)中，惩罚因子与单元材料、几何特性间的关系为：

其中：β为常数；a为单元长度；E为材料的弹性模量；h为高度。

有益效果：本发明提供了一种自适应二维有限元网格界面过渡方法，该方法不需要选取虚拟节点，采用该法对复杂系统进行分析时能采用弹性建模方式，即采用该法时不需要进行繁琐的网格过渡划分，只需要在网格间的界面处运用该法即可，与现有方法相比，本发明的二维有限元网格界面过渡方法主要有以下优势：

1)适用范围广：二维有限元疏密网格界面过渡方法适用于所有二维平面模型中网格疏密差异的协调，以及带孔板模型中网格的致密化；同时，适用于不同地域，不同分析者采用不同软件建立的有限元模型的远程对接。此外，对于任意载荷(集中载荷和分布载荷，拉伸载荷和弯曲载荷)都能适用，适用范围广；

2)占用内存小：二维有限元疏密网格界面过渡方法界面单元刚度矩阵规模较小，且具有对称性、带状性，在施加了必要的边界条件后为正定的。其有限元程序简单，占用的计算内存小；

3)计算精度高：二维有限元疏密网格界面过渡方法可以很好地消除由于疏密网格界面位移不协调所产生的误差，其计算结果与基准解吻合，计算精度高；

4)计算效率高：采用二维有限元疏密网格界面过渡方法可以消除了繁琐复杂的模型转换过程，使得有限元网格的局部加密变得易于实现。同时不需要选取虚拟节点数以及虚拟节点的插值函数，计算更加方便，计算速度加快，计算效率提高；

5)操作简便：将二维有限元疏密网格界面过渡方法通过有限元程序来实现，只需确定每种网格的材料属性、几何特性和载荷等参数即可进行计算，具有自适应性，步骤简单，操作方便；

附图说明

图1是带垂直界面的二维平板模型。

图2是拉伸中心带孔平板模型。

图3是包含过渡界面的有限元模型。

图4是常规有限元模型。

图5是包含过渡界面的有限元模型水平位移分布云图。

图6是常规有限元模型水平位移分布云图。

图7是包含过渡界面的有限元模型垂直位移分布云图。

图8是常规有限元模型垂直位移分布云图。

图9是包含过渡界面的有限元模型方向应力σ₁分布云图。

图10是常规有限元模型方向应力σ₁分布云图。

图11是包含过渡界面的有限元模型方向应力σ₂分布云图。

图12是常规有限元模型方向应力σ₂分布云图。

图13是沿y轴应力集中系数曲线图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示，一种自适应二维有限元疏密网格界面过渡方法，其特征在于，该方法包括以下步骤：

1)建立带垂直界面的二维平板的几何模型，将该模型划分为两个区域，并对两个区域进行划分网格，网格相交处为节点；该二维平板的一端固定，另一端受均布拉伸载荷f；固定端的节点为固定节点，其余节点为自由节点；所述两个区域分别为区域一、区域二；采用四边形网格进行划分，每个网格为一个单元；区域一和区域二之间为界面，区域一的另一侧为固定端。本例中采用两组四节点四边形作为所述单元来进行描述，

2)以惩罚因子γ表示系统的总势能；系统的总势能为：

$\mathrm{\π}={\mathrm{\π}}_{{\mathrm{\Ω}}_1}+{\mathrm{\π}}_{{\mathrm{\Ω}}_2}+\frac{1}{2}\mathrm{\γ}\underset{S}{\∫}{(u_1-u_2)}^{2}ds$

π对各自由度求一阶变分，并置为零，可得：

$\mathrm{\δ}\mathrm{\π}={\mathrm{\δ\π}}_{{\mathrm{\Ω}}_1}+{\mathrm{\δ\π}}_{{\mathrm{\Ω}}_2}+\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_2\right)ds+\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_1\right)ds-\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_1\right)ds-\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_2\right)ds;$

其中，区域一的积分为：

$\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_1{u}_1\right)ds=\mathrm{\γ}{\∫}_0^{2/3}\left({\mathrm{\δu}}_1{u}_1\right)ds+\mathrm{\γ}{\∫}_{2/3}^{4/3}\left({\mathrm{\δu}}_1{u}_1\right)ds+\mathrm{\γ}{\∫}_{4/3}^2\left({\mathrm{\δu}}_1{u}_1\right)ds=\mathrm{\γ}{\∫}_0^{2/3}({\mathrm{\δq}}_{11}{N}_1+{\mathrm{\δq}}_{12}{N}_2)(q_{11}{N}_1+q_{12}{N}_2)ds+$

$\mathrm{\γ}{\∫}_{2/3}^{4/3}({\mathrm{\δq}}_{12}{N}_1+{\mathrm{\δq}}_{13}{N}_2)(q_{12}{N}_1+q_{13}{N}_2)ds+\mathrm{\γ}{\∫}_{4/3}^2({\mathrm{\δq}}_{13}{N}_1+{\mathrm{\δq}}_{14}{N}_2)(q_{13}{N}_1+q_{14}{N}_2)ds;$

区域二的积分为：

$\mathrm{\γ}\underset{S}{\∫}\left({\mathrm{\δu}}_2{u}_2\right)ds=\mathrm{\γ}{\∫}_0^1\left({\mathrm{\δu}}_2{u}_2\right)ds+\mathrm{\γ}{\∫}_1^2\left({\mathrm{\δu}}_2{u}_2\right)ds=\mathrm{\γ}{\∫}_0^1({\mathrm{\δq}}_{21}{N}_1+{\mathrm{\δq}}_{22}{N}_2)(q_{21}{N}_1+q_{22}{N}_2)ds+$ $\mathrm{\γ}{\∫}_1^2({\mathrm{\δq}}_{22}{N}_1+{\mathrm{\δq}}_{23}{N}_2)(q_{22}{N}_1+q_{23}{N}_2)ds;$

其中：

δ为一阶变分标志；

和π_Ω2分别为区域一和区域二的势能；

μ₁和μ₂分别为区域一和区域二相邻界面处的位移；

N₁和N₂分别为区域一和区域二的位移插值函数；

区域一积分中的自由q₁₁、q₁₂、q₁₃、q₁₄分别是区域一与区域二相邻界面处区域一上的节点的轴向位移；本实例中具体与图1中节点的对应关系为：节点2、4、6、8分别对应q₁₁、q₁₂、q₁₃、q₁₄

区域二积分中的自由q₂₁、q₂₂、q₂₃分别是区域二与区域一相邻界面处区域二上的节点的轴向位移；本实例中具体与图1中节点的对应关系为：节点9、11、13分别对应q₂₁、q₂₂、q₂₃。

3)采用传统有限元方法确定不含虚拟节点的罚函数的有限元方程；进一步的，不含虚拟节点的罚函数的有限元方程为：

$\left[\begin{array}{cccc}{K}_1^{oo}& K_1^{oi}& 0& 0\\ K_1^{io}& K_1^{ii}+G_1^{ii}& -G& 0\\ 0& -G^T& G_2^{ii}+K_2^{ii}& K_2^{io}\\ 0& 0& K_2^{oi}& K_2^{oo}\end{array}\right]\left\{\begin{array}{c}{q}_1^0\\ q_1^i\\ q_2^i\\ q_2^0\end{array}\right\}=\left[\begin{array}{c}{f}_1^i\\ f_1^o\\ f_2^i\\ f_2^o\end{array}\right]$

其中：

K₁为区域一所有单元的刚度矩阵；K₂为区域二所有单元的刚度矩阵；

上标i表示界面上的单元的自由度；上标o表示不在界面上的单元的自由度；

G由下列积分来计算：

$G_1^{ii}=\mathrm{\γ}{\∫}_s{N}_1^T{N}_{1}ds$

$G_2^{ii}=\mathrm{\γ}{\∫}_s{N}_2^T{N}_{2}ds$

$G=\mathrm{\γ}{\∫}_s{N}_1^T{N}_{2}ds.$

4)计算界面单元刚度矩阵；进一步的，区域一中的固定节点自由度等于0；界面单元的刚度矩阵为：

$\left[\begin{array}{ccccccc}\frac{2\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& 0& -\frac{7\mathrm{\γ}}{27}& -\frac{2\mathrm{\γ}}{27}& 0\\ \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& -\frac{25\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{\mathrm{\γ}}{108}\\ 0& \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& -\frac{\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{25\mathrm{\γ}}{108}\\ 0& 0& \frac{\mathrm{\γ}}{9}& \frac{4\mathrm{\γ}}{9}& 0& -\frac{2\mathrm{\γ}}{27}& -\frac{7\mathrm{\γ}}{27}\\ -\frac{7\mathrm{\γ}}{27}& -\frac{25\mathrm{\γ}}{108}& -\frac{\mathrm{\γ}}{108}& 0& \frac{\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& 0\\ -\frac{2\mathrm{\γ}}{27}& -\frac{23\mathrm{\γ}}{54}& -\frac{23\mathrm{\γ}}{54}& -\frac{2\mathrm{\γ}}{27}& \frac{\mathrm{\γ}}{6}& \frac{2\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}\\ 0& -\frac{\mathrm{\γ}}{108}& -\frac{25\mathrm{\γ}}{108}& -\frac{7\mathrm{\γ}}{27}& 0& \frac{\mathrm{\γ}}{6}& \frac{\mathrm{\γ}}{3}\end{array}\right].$

5)计算系统刚度矩阵；系统整体刚度矩阵为：

$K=\left[\begin{array}{cccccccccc}\frac{1}{2}+\frac{2\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& 0& -\frac{7\mathrm{\γ}}{27}& -\frac{2\mathrm{\γ}}{27}& 0& 0& 0& 0\\ \frac{\mathrm{\γ}}{9}& 1+\frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& 0& -\frac{25\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{\mathrm{\γ}}{108}& 0& 0& 0\\ 0& \frac{\mathrm{\γ}}{9}& 1+\frac{4\mathrm{\γ}}{9}& \frac{\mathrm{\γ}}{9}& -\frac{\mathrm{\γ}}{108}& -\frac{23\mathrm{\γ}}{54}& -\frac{25\mathrm{\γ}}{108}& 0& 0& 0\\ 0& 0& \frac{\mathrm{\γ}}{9}& \frac{1}{2}+\frac{2\mathrm{\γ}}{9}& 0& -\frac{2\mathrm{\γ}}{27}& -\frac{7\mathrm{\γ}}{27}& 0& 0& 0\\ -\frac{7\mathrm{\γ}}{27}& -\frac{25\mathrm{\γ}}{108}& -\frac{\mathrm{\γ}}{108}& 0& \frac{1}{2}+\frac{\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& 0& -\frac{1}{4}& -\frac{1}{4}& 0\\ -\frac{2\mathrm{\γ}}{27}& -\frac{23\mathrm{\γ}}{54}& -\frac{23\mathrm{\γ}}{54}& -\frac{2\mathrm{\γ}}{57}& \frac{\mathrm{\γ}}{6}& 1+\frac{2\mathrm{\γ}}{3}& \frac{\mathrm{\γ}}{6}& -\frac{1}{4}& -\frac{1}{2}& -\frac{1}{4}\\ 0& -\frac{\mathrm{\γ}}{108}& -\frac{25\mathrm{\γ}}{108}& -\frac{7\mathrm{\γ}}{27}& 0& \frac{\mathrm{\γ}}{6}& \frac{1}{2}+\frac{\mathrm{\γ}}{3}& 0& -\frac{1}{4}& -\frac{1}{4}\\ 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{4}& 0& \frac{1}{2}& 0& 0\\ 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{2}& -\frac{1}{4}& 0& 1& 0\\ 0& 0& 0& 0& 0& -\frac{1}{4}& -\frac{1}{4}& 0& 0& \frac{1}{2}\end{array}\right].$

6)确定惩罚因子γ；惩罚因子与单元材料、几何特性间的关系为：

其中：β为常数；a为单元长度；E为材料的弹性模量；h为高度。这个表达式使解的精确性只依赖于参数β的值。界面两侧区域可得到两个不同的惩罚因子，但是无虚拟节点罚函数位移协调法中只需要一个惩罚因子，故取两侧单元惩罚因子的平均值作为惩罚因子。

7)采用常规方法计算出位移和应力应变：具体的，举例来说可以通过商用有限元软件ABAQUS中的用户单元子程序UEL接口来实现。

子程序从输入文件中读取连接的两个网格的材料、几何特性等必需信息，包括：界面上的节点、两个区域的厚度、杨氏模量。根据这些参数，程序可计算出G矩阵，并组合成界面单元刚度矩阵。从而，可以进行各种计算，验证界面元的有效性。

根据ABAQUS用户手册的要求，UEL子程序用FORTRAN77编写。为尽量减小数值截断误差，所有的变量定义为双精度型。

用户单元子程序格式如下：

实施例：

如附图2所示几何尺寸平板两端承受均布拉伸载荷f，板的长度和宽度分别为：

b＝200，h＝100，孔的半径r＝2.5；界面到圆心的距离为5。

利用对称性，只对板的1/4建模；在孔与界面间采用较细致的网格，而其它部分采用较稀疏的网格。使用有限元软件分别建立了包含节点的常规有限元模型(见附图3)和采用了不包含节点的界面过渡方法的模型(见附图4)进行对比分析。网格划分采用平面应力四节点四边形单元CSP4。界面的每一部分应用5个界面单元。

两种方法的计算结果基本一致。越过界面水平位移U₁(见附图5和附图6)和垂直位移U₂(见附图7和附图8)不显示任何变化。因此，界面单元不影响位移的分布。比较方向应力σ₁的分布云图(见附图9和附图10)和方向应力σ₂的分布云图(见附图11和附图12)可知，不同的模型间只有很小的差别。

中心带圆孔的无限大平板在均布拉伸载荷下的弹性解表明孔上边缘处的应力集中系数K_t＝3.0。附图13中，横坐标为y轴坐标，纵坐标为应力集中系数，将弹性精确解与不包含节点的界面元模型和常规有限元模型进行对比。结果非常接近，确认了该有限元界面过渡方法的精确性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (7)

1.一种自适应二维有限元疏密网格界面过渡方法，其特征在于，该方法包括以下步骤：

1)建立带垂直界面的二维平板的几何模型，将该模型划分为两个区域，并对两个区域进行划分网格，网格相交处为节点；该二维平板的一端固定，另一端受均布拉伸载荷f；固定端的节点为固定节点，其余节点为自由节点；

2)以惩罚因子γ表示系统的总势能；

3)采用传统有限元方法确定不含虚拟节点的罚函数的有限元方程；

4)计算界面单元刚度矩阵；

5)计算系统刚度矩阵；

6)确定惩罚因子γ；

7)采用常规方法计算出位移和应力应变。

2.如权利要求1所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，所述步骤1)中，所述两个区域分别为区域一、区域二；采用四边形网格进行划分，每个网格为一个单元；区域一和区域二之间为界面，区域一的另一侧为固定端。

3.如权利要求2所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，步骤2)中，系统的总势能为：

π对各自由度求一阶变分，并置为零，可得：

其中，区域一的积分为：

区域二的积分为：

其中：

δ为一阶变分标志；

和π_Ω2分别为区域一和区域二的势能；

μ₁和μ₂分别为区域一和区域二相邻界面处的位移；

N₁和N₂分别为区域一和区域二的位移插值函数；

区域一积分中的自由q₁₁、q₁₂、q₁₃、q₁₄分别是区域一与区域二相邻界面处区域一上的节点的轴向位移；

区域二积分中的自由q₂₁、q₂₂、q₂₃分别是区域二与区域一相邻界面处区域二上的节点的轴向位移。

4.如权利要求2所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，步骤3)中，不含虚拟节点的罚函数的有限元方程为：

其中：

K₁为区域一所有单元的刚度矩阵；

K₂为区域二所有单元的刚度矩阵；

上标i表示界面上的单元的自由度；上标o表示不在界面上的单元的自由度；

G由下列积分来计算：

5.如权利要求2所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，步骤4)中，区域一中的固定节点自由度等于0；界面单元的刚度矩阵为：

6.如权利要求2所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，步骤5)中，系统整体刚度矩阵为：

7.如权利要求2所述的一种自适应二维有限元疏密网格界面过渡方法，其特征在于，步骤6)中，惩罚因子与单元材料、几何特性间的关系为：

其中：

β为常数；

a为单元长度；

E为材料的弹性模量；

h为高度。