CN107451371B - 提高三角单元计算精度的方法 - Google Patents

Abstract

本发明涉及一种提高三角单元计算精度的方法，本发明设计了一种全新的三角单元，作为传统有限单元法的补充。本发明方法包括：将几何体划分为若干个三角网格；与几何体边界有交集的单元定义为边界单元，反之则为体内单元，若该三角单元是体内单元W_e，则与该体内单元共边的相邻单元分别是

则将单元W_e的插值结点取为

这4个单元的所有节点的集合；若该三角单元是边界单元W_e，则将边界单元做镜像，形成了1～2个虚单元，产生了1～2个虚节点，边界单元W_e的插值节点选取边界单元W_e及其相邻单元的所有实节点和所有虚节点中的i₁～i₆即为边界单元W_e的插值节点；采用二阶点插值方法进行三角单元的插值。本发明在工程界有一定的应用前景。

Description

提高三角单元计算精度的方法

技术领域

本发明涉及一种提高三角单元计算精度的方法。

背景技术

当今世界最强大的数值计算工具非有限单元法莫属了，它广泛运用与土木工程、航空、航天工程、机械制造工程，它的灵活性、鲁棒性、高精度性等优良特点使之备受广大工程师以及科研工作者的青睐。经过几十年的发展，有限元已经衍生出了大量不同类型的方法，根据单元形状的不同，有三角单元法(简称T3单元)，四边形等参元(Q4单元)等。

但是，以上传统单元存在着各自的优缺点，其中包括：

对于三角形T3单元：

优点：拓扑适应性很高，对于实际工程中复杂的几何体都可以用三角单元进行离散；插值简单，由于每个三角单元仅有3个插值节点，可以使用最简单的一次多项式对场变量进行逼近模拟，计算效率较高；

缺点：传统的三角形单元插值精度非常低，三角形单元的应力与应变分量皆为常数，与计算点的位置无关，这与实际不符合，因此,极低的计算精度成了三角单元在工程中大量使用的最大“瓶颈”，目前突破这个“瓶颈”的最主要方法是增大网格划分密度，然而这样带来的后果就是，精度提高幅度有限，但计算负担确大幅增大。若不大幅增加三角单元的个数，提高三角单元计算精度的传统方法是增加单元节点总数，在这个背景下，T6单元应运而生，即对每个三角单元，增加三条边的中点作为插值节点，这样的做法对提高三角单元的计算精度确实有帮助，效果也相当明显，但不良后果就是“T6单元的节点数暴增，是原T3单元的3倍，由此形成的单元刚度矩阵也非常大，每个单元刚度矩阵的宽度是原T3单元的3倍，在求解方程的时候导致计算效率降低了60％左右。

对于四边形角形Q4单元：

优点：单元计算精度高于T3单元，每个Q4单元是由4个节点首尾相连形成的四边形单元，这种单元由于插值节点数增加到了4个点，根据经典的拉格朗日插值理论，Q4单元可以使用非完备二次多项式进行插值模拟，计算精度较T3单元的C0阶提高到了C1阶，因此，是目前工程界的“宠儿”；

Q4单元的刚度矩阵的带宽小于T6单元，但其计算精度接近T6单元；

缺点：几何拓扑适应性较差。对于“凸”形几何体，可以用四边形单元较好的进行离散，但对于“凹”形几何体而言，用四边形单元划分计算域就会遇到困难，四边形单元无法适应“凹”形几何体的边界，从而使数值模拟出现偏差甚至错误，若要很好的适应凹”形几何边界，就需要在边界区域细分网格，这为网格划分工作带来了巨大的困难，同时也增加了计算负担。

发明内容

为解决上述技术问题，本发明的目的是提供一种在不增加单元节点的基础上，提高三角单元计算精度的方法。

本发明提高三角单元计算精度的方法，包括：

将实际工程中的几何体划分为若干个三角单元；

判断三角单元是否为体内单元；

若为体内单元，可搜索出该体内单元共边的相邻单元

则这4个单元的所有节点的集合i₁,i₂,i₃,i₄,i₅,i₆取为该体内单元的插值节点；

若该三角单元非体内单元，则为边界单元，将该边界单元向几何体边界外做镜像，形成了1～2个虚单元，同时产生了1～2个虚节点，边界单元的插值节点选取如下6个节点：边界单元及其相邻单元的所有实节点和所有虚节点的集合i₁～i₆即为边界单元的插值节点；

当确定了每个三角单元的插值节点i₁～i₆后，采用二阶点插值方法进行三角单元的插值。

进一步地，所述的二阶点插值方法具体包括：

三角单元位移函数用选定的6个插值节点的节点位移u₁～u₆表示为

其中

N_i——由二次多项式点插值法构造的形函数，即

其中：

式中，

——x＝{x,y}为三角单元Ω_e内任一节点的坐标；

——x₁,…,x₆，为选取的三角单元Ω_e的6个插值节点横坐标；

——y₁,…,y₆，为选取的三角单元Ω_e的6个插值节点纵坐标。

借由上述方案，本发明至少具有以下优点：

本发明重新定义了三角单元的插值节点，由经典的线性单元的3个插值节点扩展到了相邻单元覆盖的所有节点，这样处理的优点就是不会额外增加求解域的节点总数，不会大幅增加控制方程组的大小，但计算精度要远高于经典的线性三角(T3单元，或简称LFEM),与T6单元一样具有2阶插值精度，能通过二阶小片检验，但计算效率又明显高于T6单元，因此，本发明是对传统有限单元法的补充和改进，在工程界有一定的应用前景。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。

附图说明

图1是传统的T3单元的插值节点的选取方案；

图2是传统的六节点三角形单元(T6)节点选择方案；

图3是新型三角三元插值节点选取方案；

图4是边界单元；

图5是二阶小片检验的计算模型；

图6是LFEM、T33单元的竖向应力解与精确解的对比图；

图7是受端部荷载作用的悬臂梁；

图8是经典的三角单元法与T33单元法的中轴线竖向位移对比图；

图9是中心有孔无限大平板受水平方向均布荷载作用；

图10是中心圆孔方板所使用的计算模型(取1/4板作为计算对象)；

图11是中心圆孔方板单元布置；

图12是经典三角单元法与本发明方法求得的1/4板底各节点(x＝0)竖向位移对比图；

图13是厚壁圆筒计算模型；

图14是厚壁圆筒的单元划分；

图15是T3、T6、T66等三角单元法与本发明方法(T33单元法)计算得到的1/4圆筒底端(y＝0)处沿水平方向的位移场的计算结果；

图16是均质边坡模型；

图17是均质边坡单元划分；

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

说明：为了与传统有限元做区隔，本发明提出的新型单元类型可简称为T33单元。

如图1所示，传统的T3单元的插值节点的选取方案，对于任一个三角T3单元Ω_e,其传统的插值节点选取方法是：

选取其三个角点i₁,i₂,i₃作为插值节点。

而对于传统的T6单元，其提高三角单元精度的方法是在原有角点i₁,i₂,i₃的基础上再额外增加每条边的中点c₁、c₂、c₃作为插值节点，如图2所示传统的六节点三角形单元(T6)节点选择方案。

这样，每个三角形单元的插值点增加到了6个，但节点总数也增加到6个。

可以看出，对于三角单元，增加插值节点总数，便可以提高其计算精度。

本发明提高三角单元计算精度的方法，对每个三角单元的插值节点选择方式进行变革，采用如图3所示的插值节点选择方式：本发明新型三角三元插值节点选取方案。具体方案介绍如下：

确定单元为体内单元或边界单元。

对于任一体内单元：

与其共边的相邻单元分别是

则将体内单元的插值节点取为

这4个单元的所有节点的集合，即i₁,i₂,i₃,i₄,i₅,i₆

对于边界单元：

该边界单元位于求解域的边界上，由图4边界单元：

可以看出，这种边界单元的相邻单元仅有1～2个，则该边界单元与相邻单元的节点总数仅有4～5个，不足6个，根据经典的数值逼近理论可知，若要使该边界单元具有二次插值精度，则插值节点总数需6个，这就意味着对于边界单元来说，目前的插值节点选取方案还不足以完成三角单元的二次插值模拟，这就需要对原方案进行修正，以便达到选取6个插值节点的目的。

本发明的修正方法是：

将边界单元沿求解域边界做镜像，形成了1～2个“虚”单元，同时产生了1～2个“虚”节点，如图4中用○标记的节点，则边界单元的插值节点选取如下6个节点：

边界单元及其相邻单元的所有实节点(图4●所示)+所有虚节点(图4○所示)，图4中的i₁～i₆即为边界单元的插值节点。

单元插值方案

当确定了每个三角单元的插值节点i₁～i₆后，可以选取二阶点插值作为新型单元的插值方案，该方案与经典的有限元方案有着较大的差别，具体插值方法如下：

三角单元位移函数可以用上述选定的6个插值节点的节点位移u₁～u₆近似表示为

其中

N_i——由二次多项式点插值法构造的形函数，即

其中：

p^T(x)＝(1,x,y,x²,xy,y²) (3)

式中，

——x＝{x,y}为三角单元Ω_e内任一节点的坐标

——x₁,…,x₆，按照上述方案选取的三角单元Ω_e的6个插值节点横坐标；

——y₁,…,y₆，按照上述方案选取的三角单元Ω_e的6个插值节点纵坐标；

由式(2)可以看出，本发明采用的点插值形函数为完备的二次多项式，其一阶偏导数为完备的一次多项式，因此该形函数具有C1阶插值精度，较传统的三角单元形函数仅具有C0阶精度提高了一大步。

另一方面，本发明三角单元在进度上与髙阶三角单元——T6单元在计算精度上是一致的，但节点总数却大幅低于T6单元，有限元系统方程的大小取决于求解域内节点总数求解域内节点总数越少，系统方程组越小，方程组的求解效率就越高，由此可见，本方案的计算效率高于传统的三角形T6单元法

本发明在计算精度与计算效率方面均优于传统的三角单元法，为与传统有限元法进行区隔，将本方案的三角单元命名为T33单元。本发明在少增加甚至不增加单元节点总数的前提下，增加每个三角单元的插值节点总数，从而达到既提高计算精度，又不增大计算负担目的。下面从经典算例与工程算例两方面验证T33单元的有效性与优越性。

具体实施例：

实施例1

本发明采用二阶小片检验具体包括如下：

如图5所示是一个矩形方板，边长L＝10m，方板的水平方向及底板均受到法向位移约束，这个方板仅受重力作用，容重γ＝9.8N/m³，弹性模量E＝9.8N/m³，泊松比ν＝0.3.

这个算例的竖向应力σ_y精确解为

σ_y＝γ(H-y) (5)

其中，γ——方板容重

H——板高

y——计算点纵坐标

将经典的线性三角单元(简写LFEM)及本发明单元(T33)与精确解(EXACT)做了个比较，图6展示的是三者的吻合情况。

图6展示的结果非常明显，本发明提出的T33单元与精确解较吻合，然而经典的LFEM与精确解偏差较多，表明本发明T33单元能通过著名的二次小片检验，具有二阶插值精度，经典的LFEM由于插值精度所限，是无法通过二次小片检验的。

实施例2

本发明采用受均布荷载作用的悬臂梁检验具体包括如下：

如图7所示，无重悬臂梁右端受均布荷载作用，其解析解为

其中梁的惯性矩为

并且

应力为

其中，

P——均布荷载；

E——弹性模量；

v——珀松比；

I——转动惯量；

L——悬臂梁长度；

c——悬臂梁半高；

x,y——计算点的横、纵坐标；

将本发明T33单元与经典的T3单元T6单元的结果进行比较，其中T6单元的计算结果取自ABAQUS，以上几种数值计算方法得到的悬臂梁中轴线上各点的竖向位移与精确解的比较结果如图8所示。

图8经典的三角单元法与本发明T33单元法的中轴线竖向位移对比图。

在图8中，虽然以上几种方法与精确解均能有较好的吻合，但相比之下，本发明T33单元的计算结果g更接近精确解，表明本发明T33单元的精度略高于其他解。

实施例3

本发明采用具有中心圆孔的无限大平板检验，具体包括如下：

考虑一个中心有孔的无限大平板，在无穷远处受水平方向均布荷载q作用，其中，q＝1，见图9所示。

将坐标原点取在圆孔中心，由于结构的对称性，仅取其中1/4进行研究，如图9所示，边长为D，中心圆孔半径为a，孔边缘为自由边界，左边界和下边界分别固定x方向位移和y方向位移，上表面和右表面为指定面力边界。

弹性力学给出的解析解为：

其中

E为材料的弹性模量，ν是泊松比。

(r,θ)是相对于原点在圆孔中心的坐标系的极坐标，r是极径，r²＝x²+y²，θ是从x轴正向沿逆时针方向算起，

中心孔径a＝1，材料参数取为E＝1000Pa，泊松比μ＝0.3，共布下361个三角单元，如图11所示。

分别用T6，T3及本发明T33单元求得的左边界(x＝0)各节点的竖向位移，如图12所示，图12各种三角单元法求得的1/4板底各节点(x＝0)竖向位移对比图。在图12中，T33单元较其他类型的三角单元具有更高的精度。

实施例4

本发明解决受均匀内压的厚壁圆筒问题具体包括如下：

考察一个厚壁圆筒，内径a＝1m，外径b＝5m，承受内压p＝1.5×10⁸N/m²，由于轴对称，则取1/4作为研究对象，平面应变问题，如图13所示，弹性模量E＝2.1×10⁹N/m²，泊松比v＝0.3，

应力和位移的精确解为

其中：

a——内径

b——外径

p——内压

在此1/4圆筒上布置了189个节点，划分了308个三角单元，如图14所示。

分别采用经典的三节点线性三角单元(T3单元)，6节点二次三角单元(T6单元)以及本方案的二阶点插值三角单元(T33单元)做结果比较。其中，在采用T6单元做计算分析时，在全域需要布置685个节点，而T33单元则在原有的三角单元的基础上仅需额外添加68个“虚”节点做场变量插值之用，节点数远少于T6单元。

图15展示的是上述T3、T6、T66三种三角单元计算得到的1/4圆筒底端(y＝0)处沿水平方向的位移场的计算结果。

图15表明，本方案提出的T33方法与精确解有着非常高的吻合性，而其他方法在两端的精度稍有逊色。

实施例5

本发明采用工程算例验证具体包括如下：

如图16所示一均质边坡模型，假定凝聚力c＝3kPa，摩擦角φ＝20ο，容重γ＝20kN/m³，杨氏模量E＝200MPa，泊松比v＝0.4。材料为各向同性，理想弹塑性，符合Mohr-Coulomb准则和关联流动法则。在边坡上布置了241个节点，同时被划分为206个单元，如图17所示均质边坡单元划分。

使用高精度的6节点T6单元与本发明的T33单元的塑性区的分布情况做了一个比较，其中塑性区均是采用强度折减法获取，安全系数取为1.39。

由结果可以知道，本发明T33单元的塑性区分布与商用软件的高精度单元得到是塑性区分布特征基本一致，验证了本方法的工程实用性。

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

Claims (1)

1.一种提高三角单元计算精度的方法，其特征在于，包括：

将实际工程中的二阶小片划分为若干个三角单元；

判断三角单元是否为体内单元；

若为体内单元，可搜索出该体内单元共边的相邻单元

则这4个单元的所有节点的集合i₁,i₂,i₃,i₄,i₅,i₆取为该体内单元的插值节点；

若该三角单元非体内单元，则为边界单元，将该边界单元向几何体边界外做镜像，形成了1～2个虚单元，同时产生了1～2个虚节点，边界单元的插值节点选取如下6个节点：边界单元及其相邻单元的所有实节点和所有虚节点的集合i₁～i₆即为边界单元的插值节点；

当确定了每个三角单元的插值节点i₁～i₆后，采用二阶点插值方法进行三角单元的插值；

所述的二阶点插值方法具体包括：

三角单元位移函数用选定的6个插值节点的节点位移u₁～u₆表示为

其中

N_i——由二次多项式点插值法构造的形函数，即

其中：

p^T(X)＝(1,x,y,x²,xy,y²) (3)

式中，

——X＝{x,y}为三角单元Ω_e内任一节点的坐标；

——x₁,…,x₆，为选取的三角单元Ω_e的6个插值节点横坐标；

——y₁,…,y₆，为选取的三角单元Ω_e的6个插值节点纵坐标；

所述二阶小片的竖向应力σ_y精确解为

σ_y＝γ(H-y) (5)

其中，γ——方板容重

H——板高

y——计算点纵坐标，即表示三角单元Ω_e内任一节点的纵坐标。

Images