一种抗畸变的增强型八节点等参元

乔海青; 鲍四元; 沈峰

doi:10.13471/j.cnki.acta.snus.ZR20230003

欢迎访问《中山大学学报(自然科学版)(中英文)》！ English Version 维普资讯中国知网万方数据

您当前的位置：

首页 >

文章列表页 >

一种抗畸变的增强型八节点等参元

研究论文 | 更新时间：2024-05-30

- 一种抗畸变的增强型八节点等参元
- An enhanced 8-node isoparametric element with anti-distortion
- 角色： First author , 第一作者 ,
  
  机构：
  苏州科技大学土木工程学院，江苏苏州 215011
  
  邮箱： 2640403550@qq.com
  
  简介：乔海青（1999年生），男；研究方向：有限元方法及理论；E-mail：2640403550@qq.com
  
  暂无本作者相关信息
  乔海青
  ，
  角色： Corresponding author , 通讯作者 ,
  
  机构：
  苏州科技大学土木工程学院，江苏苏州 215011
  
  邮箱： bsiyuan@126.com
  
  简介：鲍四元（1980年生），男；研究方向：辛方法，动力学中的数值方法；E-mail：bsiyuan@126.com
  
  暂无本作者相关信息
  鲍四元
  ，
  角色：
  
  机构：
  苏州科技大学土木工程学院，江苏苏州 215011
  
  邮箱：
  
  暂无本作者相关信息
  沈峰
  ，
- 中山大学学报(自然科学版)(中英文) 2024年63卷第3期页码：163-171
- 作者机构：
  
  苏州科技大学土木工程学院，江苏苏州 215011
- 作者简介：
  
  乔海青（1999年生），男；研究方向：有限元方法及理论；E-mail：2640403550@qq.com
  鲍四元（1980年生），男；研究方向：辛方法，动力学中的数值方法；E-mail：bsiyuan@126.com
- 基金信息：
  
  江苏省高校青蓝工程(2022GX05);江苏省研究生科研创新与实践计划项目(KYCX23_3332)
- DOI：10.13471/j.cnki.acta.snus.ZR20230003
  中图分类号： O34
- 纸质出版日期：2024-05-25，
  
  网络出版日期：2024-01-24，
扫描看全文
引用本文

阅读全文PDF
乔海青,鲍四元,沈峰.一种抗畸变的增强型八节点等参元[J].中山大学学报(自然科学版)(中英文),2024,63(03):163-171.

QIAO Haiqing,BAO Siyuan,SHEN Feng.An enhanced 8-node isoparametric element with anti-distortion[J].Acta Scientiarum Naturalium Universitatis Sunyatseni,2024,63(03):163-171.
乔海青,鲍四元,沈峰.一种抗畸变的增强型八节点等参元[J].中山大学学报(自然科学版)(中英文),2024,63(03):163-171. DOI： 10.13471/j.cnki.acta.snus.ZR20230003.

QIAO Haiqing,BAO Siyuan,SHEN Feng.An enhanced 8-node isoparametric element with anti-distortion[J].Acta Scientiarum Naturalium Universitatis Sunyatseni,2024,63(03):163-171. DOI： 10.13471/j.cnki.acta.snus.ZR20230003.

摘要

为了提高平面问题的求解精度，基于勒让德正交函数族引入泡状函数来增强映射关系的思想，构造Q8_Legendre单元。针对单元的位移场中含有泡状函数导致单元刚度矩阵规模变大的问题，在得到单元刚度矩阵后附加自由度凝聚，以减小整体刚度矩阵规模、节省计算机内存空间。算例表明，随着泡状函数项数的不断增加，Q8_Legendre单元能够通过分片试验，使得单元的位移精度得到改善，而应力精度保持相当。当出现畸变网格时，Q8_Legendre单元表现出抗畸变的性能，且不会出现闭锁现象。

译

Abstract

In order to improve the solution precision of the plane problem， the Q8_Legendre element is constructed based on the idea of introducing Bubble function to enhance the mapping relationship based on Legendre orthogonal function. Because the displacement field of the element contains Bubble function， the size of the element stiffness matrix becomes larger. In order to reduce the size of the overall stiffness matrix and save the computer memory space， the additional degrees of freedom are condensed after the stiffness matrix of the element is obtained. Some numerical examples show that the Legendre element can pass the patch test， and the displacement accuracy of the element is improved with the increasing of Bubble function number， while the stress accuracy remains almost the same. When the distortion mesh is encountered， the Q8_Legendre element exhibits the performance of anti-distortion ability， and locking phenomenon does not appear.

译

关键词

平面问题; 泡状函数; 位移修正模型; 高斯积分; 抗畸变

译

Keywords

plane problem; Bubble function; displacement correction model; Gauss quadrature; anti-distortion

译

有限单元法（

龙志飞等，2001；Zienkiewicz et al， 2005）是20世纪模拟工程问题时提出的强有力的工具。这种方法扩大了可求解问题的范围，但是也引入了一些新的问题，如网格质量和求解精度等。学者们积极开展高性能单元的研究工作，并发展了非协调元（Long et al.， 2009）、谱元法（Pozrikidis， 2014；Han et al.， 2020）、光滑有限元法（Liu et al.， 2007）等理论。龙驭球等（2001）对有限元的发展问题做了总结；岑松等（2017）概述了形状自由高性能有限元法研究进展。近年来，在有限元插值方法、位移场和新形状的单元方面有众多的成果，如Tatyana et al.（2022）提出一种类似于在标准有限元法中插入边界条件的思想研究泊松方程问题；Hau et al.（2021）提出连续插值多面体有限元法并应用于三维固体力学问题的分析中；Essongue et al.（2021）研究了增广有限元法的精度和收敛性；Ding et al.（2021）用超收敛恢复方法改进多边形线性有限元理论，建立多边形单元用于泊松方程的研究；Su et al.（2018）提出具有一般边界条件的功能梯度阶梯梁位移的有效公式，其在处理单元位移场时利用泡状函数进行修正；另外Przemieniecki（2009）在处理平面问题时提到用泡状函数修正单元的位移场。

译

以平面问题为例，采用8结点位移元结果可达到二次精度，但其抗畸变能力不强，所以本论文主要对常规的8节点位移元提出改进方法。由文献知，在矩形网格下8节点单元位移场从自然坐标系映射到物理坐标系下的多项式最高阶次是三次，其中所含完全多项式阶次是二次（

Lee et al.， 1993）。但当划分网格出现角度畸变时，其可以映射到物理坐标系下的最高多项式次数和完全多项式次数均是一次；当发生二次曲边畸变时，同样是一次多项式。实际问题中，在物理坐标系下平面问题的理论解往往不止是一阶多项式，因此畸变网格情况下平面问题的求解精度会降低。为了增加单元从自然坐标系映射到物理坐标下的位移场多项式最高次数和完全多项式阶次，本文对传统有限元理论中的8节点位移元的位移场提出改进方法，在单元的自然坐标系下增加泡状函数，以提高单元的性能。

译

正交函数族（

Houmat， 1997， 2000；He et al.， 2021）可以通过线性组合逼近任意连续函数，所以在工程计算中有着较广泛的应用；结合有限元理论发现以正交函数族为基础构造泡状函数是可行的。因此，本文选Legendre函数族构造泡状函数，以提升单元模拟精度。

译

本文的Q8_Legendre单元（以下简记Q8_L）研究思路为：基于有限元理论增加泡状函数建立8节点等参元的增强型位移场，按照传统有限元求解过程形成刚度矩阵和等效力矩阵；为了节省计算机内存空间，附加泡状函数系数采用凝聚法消去；并通过测试分片试验，以确认单元是否收敛；最后，通过算例测试所提出新单元的精度和抗畸变性能。

译

1 力学模型

1.1　单元的位移场

如图1所示的悬臂梁端部受集中力作用时，在物理坐标系下的解含有三次多项式；畸变网格出现时，如果采用8节点等参元，则物理坐标系中的可表达多项式阶次是一次，从而造成单元模拟精度下降（

Lee et al.， 1993）。通过增加一些泡状函数项使单元在自然坐标下完全多项式项数和阶次增加；当再映射到物理坐标系时，就会使单元的映射关系有所改善。

译

图1 受集中力的悬臂梁

Fig.1 Cantilever beam subjected to concentrated force

下载: 原图 | 高精图 | 低精图

对于平面问题，等参元Q8_Legendre在自然坐标系上父单元到物理坐标系下子单元的映射关系，如图2所示。

译

图2 父单元与子单元的映射

Fig.2 Mapping from the parent element to the sub-element

下载: 原图 | 高精图 | 低精图

图2中，自然坐标为 $ξ 和 η$ 的定义区间均为 $[- 1,1]$ 。Q8_L单元的坐标映射关系为

译

$\{\begin{array}{l} x = \sum_{i = 1}^{8} N_{i} x_{i}, \\ y = \sum_{i = 1}^{8} N_{i} y_{i} . \end{array}$

（1）

泡状函数基于Legendre正交函数族构造。以水平位移 $u$ 为例，Q8_L单元的位移场为

译

$u = \sum_{i = 1}^{8} N_{i} (ξ, η) u_{i} + \sum_{k = 1}^{n} \sum_{j = 1}^{m} c_{j k} Φ_{j} (ξ) Φ_{k} (η)$ ，

（2）

其中 $N_{i}$ （i=1，…，8）为常规8结点等参元的8个形状函数，m和n表示在自然坐标系每个维度的附加泡状函数的项数。自然坐标系中Legendre函数 $Φ_{j} (η)$ 和 $Φ_{k} (ξ)$ 的统一表达式（

Houmat， 1997）为

译

$Φ_{l} (x) = \sum_{p = 0}^{1 + [l / 2]} \frac{{(- 1)}^{p} (2 l - 5 p - 1)!}{2^{p} p! (l - 2 p + 1)} x^{l - 2 p + 1},$

（3）

其中 $l = 1,2, 3, \dots$ ， $[l / 2]$ 表示取 $l / 2$ 的整数部分， $0! = (- 1)! = 1$ .泡状函数 $Φ_{j} (η)$ 和 $Φ_{k} (ξ)$ 在单元域的边界值为0，单元域内值不为0.

译

1.2　静态凝聚及求解

本节讨论Q8_L单元的应变场和静态凝聚过程。所提单元的位移场可以统一写为

译

$U = [\begin{matrix} u \\ v \end{matrix}] = [\begin{matrix} N_{F E M} & N_{c} \end{matrix}] {[\begin{matrix} q & C_{c} \end{matrix}]}^{T}$ ，

（4）

其中N_FEM是普通有限元的形状函数矩阵，由N₁，N₂，…，N₈组成，N_c是式（1）中附加泡状函数组成的矩阵，C_c是对应的系数向量。

译

基于式（4），Q8_L单元对应的应变列式（

曾攀，2009）为

译

$[\begin{matrix} ε_{x} \\ ε_{y} \\ ε_{x y} \end{matrix}] = Q {[\begin{matrix} N_{F E M} \\ N_{c} \end{matrix}]}^{T} [\begin{matrix} q \\ C_{c} \end{matrix}]$ ，

（5）

其中

译

$Q = \frac{1}{|J|} [\begin{matrix} J_{22} & - J_{12} & 0 & 0 \\ 0 & 0 & - J_{21} & J_{11} \\ - J_{21} & J_{11} & J_{22} & - J_{12} \end{matrix}] [\begin{matrix} \frac{\partial}{\partial ξ} & 0 \\ \frac{\partial}{\partial η} & 0 \\ 0 & \frac{\partial}{\partial ξ} \\ 0 & \frac{\partial}{\partial η} \end{matrix}]$ ，

这里J是雅可比矩阵， $J_{i j} (i, j = 1,2)$ 是J的相关元素。故单元的应变关系为

译

$ε = [\begin{matrix} B_{F E M} & B_{c} \end{matrix}] {[\begin{matrix} q & C_{c} \end{matrix}]}^{T} = B {[\begin{matrix} q & C_{c} \end{matrix}]}^{T}$ ，

（6）

其中 $B_{F E M}, B_{c}$ 分别是传统有限元理论中的应变矩阵和附加泡状函数对应的应变矩阵。

译

通过虚功原理得到的刚度矩阵为

译

$\begin{array}{l} K^{e} = \iint t B^{T} D B |J| d ξ d η, \\ = \iint t [\begin{matrix} B_{F E M}^{T} D B_{F E M} & B_{F E M}^{T} D B_{c} \\ B_{c}^{T} D B_{F E M} & B_{c}^{T} D B_{c} \end{matrix}] |J| d ξ d η \\ = [\begin{matrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{matrix}], \end{array}$

（7）

其中 $K_{11}, K_{12}, K_{21}, K_{22}$ 分别为

译

$K_{11} = \iint t B_{F E M}^{T} D B_{F E M} |J| d ξ d η$ ，

$K_{12} = \iint t B_{F E M}^{T} D B_{c} |J| d ξ d η$ ，

$K_{21} = \iint t B_{c}^{T} D B_{F E M} |J| d ξ d η$ ，

$K_{22} = \iint t B_{c}^{T} D B_{c} |J| d ξ d η$ .

等效力向量的计算公式为

译

$F^{e} = \iint N^{T} [\begin{matrix} f_{x} \\ f_{y} \end{matrix}] |J| d ξ d η + \int_{S_{σ}} N^{T} [\begin{matrix} T_{x} \\ T_{y} \end{matrix}] d S_{σ} = [\begin{matrix} F_{F E M} \\ F_{c} \end{matrix}]$ ，

（8）

其中 $f_{x}$ ， $f_{y}$ 是体积分布力， $T_{x}$ 和 $T_{y}$ 是作用在单元边界的面力。每个单元的刚度矩阵方程为

译

$[\begin{matrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{matrix}] [\begin{matrix} q \\ C_{c} \end{matrix}] = [\begin{matrix} F_{F E M} \\ F_{c} \end{matrix}] + [\begin{matrix} F_{1} \\ 0 \end{matrix}]$ .

（9）

处理过程中，该单元的节点集中力和该单元的分布力分开处理，所以公式（9）的右边由两个向量组成，其中 $F_{1}$ 是作用在单元结点上的集中力， $F_{F E M}$ 是传统有限元的等效力， $F_{c}$ 是与泡状函数项相对应的分布力等效力。

译

因附加泡状项的加入，使得该单元的刚度矩阵不再是16阶方阵。为缩小整体线性方程组的规模，本文将附加系数用静态凝聚法凝聚。该单元节点位移列阵求解为

译

$(K_{11} - K_{12} K_{22}^{- 1} K_{21}) q = F_{_{F E M}} + F_{1} - K_{12} K_{22}^{- 1} F_{c}$ .

（10）

同时，新增的系数求解为

译

$C_{c} = K_{22}^{- 1} (F_{c} - K_{21} q)$ .

（11）

经过以上的凝聚过程后，可继续求解出节点位移和应力等。为了满足积分精度要求，求刚度矩阵时，子矩阵K₁₁使用5×5点的Gauss积分点，而子矩阵K₁₂，K₂₁，K₂₂使用20×20的Gauss积分点。

译

2 算例分析

2.1　收敛性测试

分片试验是检验一个单元是否收敛和健壮的数值方法。为检验本文提出的Q8_L单元是否能收敛于精确解，在测试分片试验时，给定位移为

译

$\{\begin{array}{l} u = 5 \times 10^{- 4} (x + y), \\ v = 5 \times 10^{- 4} (x + y) . \end{array}$

（12）

模型材料参数为 $E = 10^{6}$ ， $μ = 0.25$ ，厚度 $t = 0.001$ .用图3（a）和图3（b）两种不同的网格划分模式测试，相应的常应力解为

译

$\{\begin{array}{l} σ_{x} = σ_{y} = 666.666 6, \\ τ_{x y} = 400 . \end{array}$

（13）

图3 四边形网格

Fig.3 The quadrilateral mesh

下载: 原图 | 高精图 | 低精图

以平面应力问题为例，检测后发现所提单元能通过分片试验。

译

2.2　抗畸变性能和精度研究

本节用两个模型测试所提单元的性能。第一个模型注重Q8_L单元的收敛性能测试，结构及网格划分如图4所示。图4为一简支梁模型，参数 $E = 100$ ， $v = 0$ ，厚度 $t = 1$ ，集中力 $P = 1$ 。分别对三种网格应用Q8_L单元，结果表明：当单元数量固定后，随着 $m$ 、 $n$ 的增加，单元收敛于精确解。其中，Q8_L单元在图4（c）的网格划分情况下，当 $m$ 、 $n$ 分别取1、2和6时，计算精度相对于Q8单元提升了46.67%、54.40%和56.00%。图 $4 (c)$ 的网格划分情况下，所得简支梁底边中点的竖向位移 $v (5,0)$ ，如表1所示。

译

图4 简支梁不同种类的网格

Fig.4 A simply supported beam with different kinds of grids

下载: 原图 | 高精图 | 低精图

表1 均匀弯矩M作用下的结果^1）

Table 1 Numerical results under uniform bending moment M

Mesh

Element

Results

0.630

Q8_L

1.334

（m=n=1）

1.446

（m=n=2）

1.470

（m=n=6）

1） m=n=1表示Q8_L的附加项参数值取1，其他类似。

下载: CSV

下载: 表格图片

第二个模型是端部受集中力或集中力偶的悬臂梁（如图1所示）。该模型注重测试模型的抗畸变性能。模型参数 $E = 1 \times 10^{7}$ ， $v = 0.3$ ，厚度 $t = 1$ ，梁长 $L$ 取100和20两组值，梁高 $c$ 取10，按照平面应力问题考虑。结构受力分为两种情况：（1） $P = \frac{20 c^{2}}{L}$ ， $y$ 方向上有分布力 $f_{y} = \frac{120 y}{L} - \frac{120 y^{2}}{c L}$ 考虑；（2） $M = 20 c^{2}$ ， $y$ 方向上有分布力 $f_{y} = \frac{240 y}{c} - 120$ 考虑。对于 $L$ 取100和20的情况，采用的网格不同，具体如图5所示。

译

图5 悬臂梁的网格种类

Fig.5 Grid kinds of cantilever beam

下载: 原图 | 高精图 | 低精图

第二个模型主要研究当单元出现斜线和曲边时，Q8_L单元是否会有闭锁现象，以及计算精度随 $m$ 、 $n$ 的变化情况。计算所得的位移结果如表2-3所示。

译

表2 力P作用下的数值结果

Table 2 Numerical results under the load P

Displacement	Mesh	Element	Results
$v (100,0) \times 10^{3}$	1	Q8	6.107
	1	Q8_L	7.633 （m=n=1）	7.682 （m=n=2）	7.711 （m=n=3）	7.723 （m=n=4）	7.732 （m=n=6）
	2	Q12^1）/Q8	3.444/1.289
	2	Q8_L	6.335 （m=n=1）	6.430 （m=n=2）	6.474 （m=n=3）	6.499 （m=n=4）	6.518 （m=n=6）
$v (20,0) \times 10^{4}$	3	Q8	3.339
	3	Q8_L	3.502 （m=n=1）	3.599 （m=n=2）	3.585 （m=n=3）	3.605 （m=n=4）	3.612 （m=n=6）
	4	Q8	3.344
	4	Q8_L	3.335 （m=n=1）	3.619 （m=n=2）	3.646 （m=n=3）	3.656 （m=n=4）	3.658（m=n=6）

1）参见

Lee et al.（1993）。

下载: CSV

下载: 表格图片

表3 力偶M作用下的结果

Table 3 Numerical results under the coupling M

Displacement	Mesh	Element	Results
$v (100,0) \times 10^{3}$	1	Q8	-11.714
	1	Q8_L	-11.898 （m=n=1）	-11.950 （m=n=2）	-11.976 （m=n=3）	-11.986 （m=n=4）	-11.993 （m=n=6）
	2	Q12^1）/Q8	（-2.33）/（-1.896）
	2	Q8_L	-11.934 （m=n=1）	-11.974 （m=n=2）	-11.997 （m=n=3）	-12.009 （m=n=4）	-12.018 （m=n=6）
$v (20,0) \times 10^{4}$	3	Q8	-4.412
	3	Q8_L	-4.749 （m=n=1）	-4.782 （m=n=2）	-4.779 （m=n=3）	-4.785 （m=n=4）	-4.785 （m=n=5）
	4	Q12^1）/Q8	-4.412
	4	Q8_L	-4.393 （m=n=1）	-4.652 （m=n=2）	-4.681 （m=n=3）	-4.691 （m=n=4）	-4.694 （m=n=6）

1）参见

Lee et al.（1993）。

下载: CSV

下载: 表格图片

2.3　Macneal问题

Macneal问题（

Macneal，1989）是著名的有限元难题，如图6所示。参数

$E = 10^{7}$ ，

$t = 0.1$ ，

$v = 0.3$ 时，A点的竖向位移测试结果如表4-5所示。结果表明：（1）图6（c）网格情况下，使用Q8单元且计算使用3×3的Gauss积分点时出现自锁现象。对于本文提出的Q8_L单元，即使

$m$ 和

$n$ 很小也不会出现自锁现象；（2）随着

$m$ 和

$n$ 的增加，计算精度越来越高且逐步收敛于精确解。与所引用的

$P S (α)$ 等单元（Wu et al.， 1995； Huang et al.， 2004）相比，Q8_L单元计算结果的精度基本相当，甚至有些情况下超过了

$P S (α)$ 等单元。

译

图6 Macneal梁

Fig.6 Macneal beam

下载: 原图 | 高精图 | 低精图

Q8单元用2×2的Gauss积分能得到较为理想结果的原因是：缩减积分方法使得单元的刚度矩阵中元素值减小，缓解了自锁现象，但这仅仅是从积分方案上缓解。而本文提出的Q8_L单元能够从根本上改善单元的自锁现象。本算例还考虑了平面应变问题接近不可压缩的情况，篇幅原因只给出Q8_L单元当 $m = n = 6$ 时的结果，如表6-7所示。结果表明，畸变网格（图6（b）和图6（c））网格情况下，Q8_L受影响较小，同时计算精度上大大优于PS单元（

Pian et al.， 1984），且不存在闭锁现象。

译

表4 力偶M作用下Macneal问题结果^1）

Table 4 The results of Macneal problem under the couple M

Mesh	Element
	Q8（3×3）	Q8（2×2）	PS/P2-Q6-Ⅲ^2）	$P S (α)$ ^3）	Q8_L
	Q8（3×3）	Q8（2×2）	PS/P2-Q6-Ⅲ^2）	$P S (α)$ ^3）	m=n=1	m=n=3	m=n=5	m=n=6
a	0.616 8	0.989 6	1.000	1.000	0.995 9	0.998 6	0.999 2	0.999 3
b	0.564 3	0.998 1	0.163	0.999	0.996 8	0.998 8	0.999 3	0.999 3
c	0.089 1	0.984 8	0.846	1.000	0.996 9	0.998 8	0.999 3	0.999 3

1）Q8（2×2）表示Q8单元中采用2×2Gauss积分；2）参见

Huang et al.（2004）；3）参见Wu et al.（1995）。

下载: CSV

下载: 表格图片

表5 集中力P作用下Macneal问题结果^1）

Table 5 The results of Macneal problem under concentrated force P

Mesh	Element
	Q8 （3×3）	Q8（2×2）	PS/P2-Q6-Ⅲ1^1）	$P S (α$ ）^2）	Q8_L
	Q8 （3×3）	Q8（2×2）	PS/P2-Q6-Ⅲ1^1）	$P S (α$ ）^2）	m=n=1	m=n=3	m=n=5	m=n=6
a	0.550 722	0.984 218	0.993	0.993	0.981 3	0.985 8	0.986 7	0.986 9
b	0.488 326	0.997 077	0.221	0.996	0.981 1	0.987 1	0.988 2	0.988 3
c	0.0863 37	0.969 473	0.796	0.989	0.986 6	0.990 4	0.991 3	0.991 4

1）参见

Huang et al.（2004）；2）参见Wu et al.（1995）。

下载: CSV

下载: 表格图片

表6 力偶M作用下Macneal问题结果

Table 6 The results of Macneal problem under the action of couple M

Case	Mesh	Element
Case	Mesh	Q8（3×3）	PS^1）	Q8_L（m=n=6）
$μ$ =0.49	a	0.003 455	0.004 10	0.004 009
	b	0.003 598	0.000 69	0.004 024
	c	0.001 896	0.000 348	0.004 024
$μ$ =0.499	a	0.003 331 0.003 047 0.000 481	0.004 05	0.003 787
	b		0.000 68	0.003 830
	c		0.003 44	0.003 831

1）参见Pian et al.（1987）。

下载: CSV

下载: 表格图片

表7 集中力

$P$ 作用下的结果

Table 7 The results of the incompressibility under concentrated force

$P$

Case	Mesh	Element
Case	Mesh	Q8（3×3）	PS^1）	Q8_L（m=n=6）
$μ = 0.49$	a	0.062 66	0.08159	0.078 85
	b	0.065 58	0.018 65	0.078 71
	c	0.035 83	0.065 76	0.078 82
$μ$ =0.499	a	0.059 53	0.080 53	0.072 85
	b	0.052 79	0.018 44	0.073 31
	c	0.009 33	0.064 93	0.073 30

1）参见文献（Pian et al.， 1987）。

下载: CSV

下载: 表格图片

2.4　Cook梁问题

Cook梁如图7所示，材料参数 $E = 1$ ， $μ = 1 / 3$ ，厚度为 $t = 1$ 。通过计算A点的竖向位移，进一步考察Q8_L单元的性能。

译

图7 Cook斜梁

Fig.7 Cook skew beam

下载: 原图 | 高精图 | 低精图

原结构分为2×2、4×4和8×8共三种网格，其中斜梁的四条边分别均分为2、4和8段。表8再次显示，随着网格细化和m与n值的增大，A点位移解精度变的较为理想。由表8知，在网格划分为2×2情况下，Q8_L单元（m=n=1）A点竖向位移和对应的Q8单元在网格划分为4×4的情况下精度相差-1.752 9%；Q8_L单元（m=n=6）在网格划分为4×4的情况下，A点竖向位移结果和对应Q8单元在网格划分为8×8的情况相比提升了0.208 7%。图8展示了Q8_L（m=n=6）单元在两种不同网格下的竖向位移云图。

译

表8 受集中力时Cook梁的位移解

Table 8 The displacement solution of Cook beam under a concentrated force

Element	$v_{A}$
Element	2×2	4×4	8×8
Q8	23.17	23.71	23.88
AQ8-Ⅰ/Ⅱ^1）	22.98	23.74	23.89
Q8_L（m=n=1）	23.29	23.84	23.93
Q8_L（m=n=6）	23.70	23.93	23.95

1）参见文献（

Soh et al.， 2000）。

下载: CSV

下载: 表格图片

图8 竖向位移云图

Fig.8 Cloud image of vertical displacement

下载: 原图 | 高精图 | 低精图

2.5　应力集中问题

考虑平面应力集中问题（

Zienkiewicz et al.， 1992），参数E=1 000，

$μ = 0.3$ 。取该无限空间的部分区域计算对应点P（直角坐标

$x = y = 1.060 660 20$ 或极坐标

$r = 1.5, θ = 45 °$ ）的应力情况，如图9所示。

译

图9 应力集中

Fig.9 Stress concentration

下载: 原图 | 高精图 | 低精图

本例圆孔半径a取1，点P应力的理论解为 $σ_{x} = 1.148 1$ ， $σ_{y} = - 0.148 1$ ， $τ_{x y} = - 0.222$ （

Zienkiewicz et al.， 1992）。用两种不同的网格来测试这类单元所得应力解的精度，图10给出了对应的网格划分情况。求解过程中，AB和CD边上的分布力根据极坐标和直角坐标之间转换关系，求出直角坐标系下的应力，再施加到相应的边上。关于应力后处理，这里应用每一个点周边单元的应力平均值。

译

图10 使用8节点单元的两种网格划分

Fig.10 Different grids divided by 8-node elements

下载: 原图 | 高精图 | 低精图

从表9所示结果可知，Q8_L单元的应力与Q8单元对应的结果，精度总体上相当；畸变网格下有的应力精度甚至有所提高。

译

表9 点P应力解^1）

Table 9 Stress solutions for point P

Results	Mesh 1		Mesh 2
Results	Q8	Q8_L（m=n=4）	Q8	Q8_L（m=n=4）
$σ_{x}$	1.227 3（6.9%）	1.226 6（6.84%）	1.222 5（6.48%）	1.149 1（0.08%）
$σ_{y}$	-0.191 0（28.97%）	-0.190 0（28.29%）	-0.184 1（24.31%）	-0.179 3（21.07%）
$τ_{x y}$	-0.214 6（-3.42%）	-0.214 7（-3.38%）	-0.215 8（-2.88%）	-0.215 1（-3.20%）

1）括号内的数值表示相对于精确解的相对误差。

下载: CSV

下载: 表格图片

3 结论

本文所构造的Q8_Legendre单元在减少网格规模及抗畸变网格方面有很好的理论和实践价值。主要结论如下：

译

（1）从网格畸变会导致精度下降入手，通过增加以Legendre正交函数族为基础的泡状项函数来修正8节点等参元位移场，能够改进其性能。

译

（2）与Q8单元对比发现：常规网格下所提方法能够提高精度；随着 $m$ 、 $n$ 值增大提升越多，并逐步收敛于精确解；在畸变网格情况下Q8_L单元受畸变的影响小。

译

（3）相对于Q8单元，所得位移解的精度有一定的提升，约为10%~50%。和其他优秀单元（如PS）对比，位移解精度相当，甚至部分计算结果进一步提升。

译

（4）由于多了消除附加项的静力凝聚过程，Q8_L单元的时间成本比Q8单元大。但为达到同样精度，本文方法在离散结构时可缩小网格的规模，从而减小计算量。

译

参考文献

岑松，尚闫，周培蕾，等， 2017. 形状自由的高性能有限元方法研究的一些进展［J］. 工程力学， 34（3）：1-14. [百度学术]

龙志飞，岑松， 2001. 有限元法新论：原理·程序·进展［M］. 北京：中国水利水电出版社. [百度学术]

龙驭球，龙志飞，岑松， 2001. 有限元的几个问题和进展［C］//第十届全国结构工程学术会议论文集. 中国力学学会：57-74. [百度学术]

曾攀， 2009. 有限元基础教程［M］. 北京：高等教育出版社. [百度学术]

DING S Y， SHAO G J， HUANG Y Y， et al， 2021. The superconvergence gradient recovery method for linear finite element method with polygons［J］. Int J Numer Methods Eng， 122（16）： 4154-4171. [百度学术]

ESSONGUE S， COUÉGNAT G， MARTIN E， 2021. Finite element modelling of traction-free cracks： Benchmarking the augmented finite element method（AFEM）［J］. Fracture Mechanics， 253：107873. [百度学术]

HAN L， WANG J X ， LI H J， et al， 2020. A time-domain spectral element method with C¹ continuity for static and dynamic analysis of frame structures［J］. Structures， 28： 604-613. [百度学术]

HAU N N， KHANH C N， KHUONG D N， et al， 2021. A consecutive-interpolation polyhedral finite element method for solid structures［J］. Int J Numer Methods Eng， 122（20）：5692-5717. [百度学术]

HE D Z， LIU T， 2021. In-plane modal studies of arbitrary laminated triangular plates with elastic boundary constraints by the Chebyshev-Ritz approach［J］. Composite Structures， 271：114138. [百度学术]

HOUMAT A， 1997. Hierarchical finite element analysis of the vibration of membranes［J］. J Sound Vib， 201：465-472. [百度学术]

HOUMAT A， 2000. A triangular Fourier p-element for the analysis of membrane vibrations［J］. J Sound Vib， 230（1）：31-43. [百度学术]

HUANG Y Q， LI Q S， 2004. Four-node incompatible plane and axisymmetric elements with quadratic completeness in the physical space［J］. Int J Numer Methods Eng， 61（10）： 1603-1624. [百度学术]

LEE N S， BATHE K J， 1993. Effects of element distortions on the performance of isoparametric elements［J］. Int J Numer Methods Eng， 36（20）：3553-3576. [百度学术]

LIU G R， DAI K Y， NGUYEN T T， 2007. A smoothed finite element method for mechanics problems［J］. Computational Mechanics， 39（6）：859-877. [百度学术]

LONG Y Q， CEN S， LONG Z F， 2009. Advanced finite element method in structural engineering［M］. Beijing： Tsinghua University Press. [百度学术]

MACNEAL R H， 1989. Toward a defect-free four-noded membrane element ［J］. Finite Elem Anal Des， 5（1）： 31-37. [百度学术]

PIAN T H H， SUMIHARA K， 1984. Rational approach for assumed stress finite elements［J］. Int J Numer Methods Eng， 20（9）：1685-1695. [百度学术]

POZRIKIDIS C， 2014. Introduction to finite and spectral element methods using MATLAB［M］. Chapman and Hall/CRC. [百度学术]

PRZEMIENIECKI J S， 2009. Finite element structural analysis： New concepts［M］. AIAA （American Institute of Aeronautics & Astronautics）. [百度学术]

SOH A K， LONG Y Q， CEN S， 2000. Development of eight-node quadrilateral membrane elements using the area coordinates method［J］. Computational Mechanics， 25（4）：376-384. [百度学术]

SU Z， JIN G Y， YE T G， 2018. Vibration analysis of multiple-stepped functionally graded beams with general boundary conditions［J］. Composite Structures， 186（8）：315-323. [百度学术]

TATYANA S， ZHANG S Y， 2022. An interpolated Galerkin finite element method for the Poisson equation［J］. J Sci Comput， 92（2）：47. [百度学术]

WU C C， CHEUNG Y K， 1995. On optimization approaches of hybrid stress elements［J］. Finite Elem Anal Des， 21（1）：111-128. [百度学术]

ZIENKIEWICZ O C， ZHU J Z， 1992. The superconvergent patch recovery and a posteriori error estimates. Part 1： The recovery technique［J］. Int J Numer Methods Eng， 33（7）：1331-1364. [百度学术]

ZIENKIEWICZ O C， TAYLOR R L， 2005. The finite element method for solid and structural mechanics［M］. Oxford： Elsevier Butterworth-Heinemann. [百度学术]

摘要

Abstract

关键词

平面问题泡状函数位移修正模型高斯积分抗畸变

Keywords

plane problemBubble functiondisplacement correction modelGauss quadratureanti-distortion

references

岑松，尚闫，周培蕾，等， 2017. 形状自由的高性能有限元方法研究的一些进展［J］. 工程力学， 34（3）：1-14.

龙志飞，岑松， 2001. 有限元法新论：原理·程序·进展［M］. 北京：中国水利水电出版社.

龙驭球，龙志飞，岑松， 2001. 有限元的几个问题和进展［C］//第十届全国结构工程学术会议论文集. 中国力学学会：57-74.

曾攀， 2009. 有限元基础教程［M］. 北京：高等教育出版社.

HAN L， WANG J X ， LI H J， et al， 2020. A time-domain spectral element method with C1 continuity for static and dynamic analysis of frame structures［J］. Structures， 28： 604-613.

HAU N N， KHANH C N， KHUONG D N， et al， 2021. A consecutive-interpolation polyhedral finite element method for solid structures［J］. Int J Numer Methods Eng， 122（20）：5692-5717.

HE D Z， LIU T， 2021. In-plane modal studies of arbitrary laminated triangular plates with elastic boundary constraints by the Chebyshev-Ritz approach［J］. Composite Structures， 271：114138.

HOUMAT A， 1997. Hierarchical finite element analysis of the vibration of membranes［J］. J Sound Vib， 201：465-472.

HOUMAT A， 2000. A triangular Fourier p-element for the analysis of membrane vibrations［J］. J Sound Vib， 230（1）：31-43.

HUANG Y Q， LI Q S， 2004. Four-node incompatible plane and axisymmetric elements with quadratic completeness in the physical space［J］. Int J Numer Methods Eng， 61（10）： 1603-1624.

LEE N S， BATHE K J， 1993. Effects of element distortions on the performance of isoparametric elements［J］. Int J Numer Methods Eng， 36（20）：3553-3576.

LIU G R， DAI K Y， NGUYEN T T， 2007. A smoothed finite element method for mechanics problems［J］. Computational Mechanics， 39（6）：859-877.

LONG Y Q， CEN S， LONG Z F， 2009. Advanced finite element method in structural engineering［M］. Beijing： Tsinghua University Press.

MACNEAL R H， 1989. Toward a defect-free four-noded membrane element ［J］. Finite Elem Anal Des， 5（1）： 31-37.

PIAN T H H， SUMIHARA K， 1984. Rational approach for assumed stress finite elements［J］. Int J Numer Methods Eng， 20（9）：1685-1695.

POZRIKIDIS C， 2014. Introduction to finite and spectral element methods using MATLAB［M］. Chapman and Hall/CRC.

PRZEMIENIECKI J S， 2009. Finite element structural analysis： New concepts［M］. AIAA （American Institute of Aeronautics & Astronautics）.

SOH A K， LONG Y Q， CEN S， 2000. Development of eight-node quadrilateral membrane elements using the area coordinates method［J］. Computational Mechanics， 25（4）：376-384.

SU Z， JIN G Y， YE T G， 2018. Vibration analysis of multiple-stepped functionally graded beams with general boundary conditions［J］. Composite Structures， 186（8）：315-323.

TATYANA S， ZHANG S Y， 2022. An interpolated Galerkin finite element method for the Poisson equation［J］. J Sci Comput， 92（2）：47.

WU C C， CHEUNG Y K， 1995. On optimization approaches of hybrid stress elements［J］. Finite Elem Anal Des， 21（1）：111-128.

ZIENKIEWICZ O C， ZHU J Z， 1992. The superconvergent patch recovery and a posteriori error estimates. Part 1： The recovery technique［J］. Int J Numer Methods Eng， 33（7）：1331-1364.

ZIENKIEWICZ O C， TAYLOR R L， 2005. The finite element method for solid and structural mechanics［M］. Oxford： Elsevier Butterworth-Heinemann.

浏览量

191

下载量

CSCD

文章被引用时，请邮件提醒。

提交

工具集

关联资源

暂无数据