复杂三维曲梁结构的无闭锁等几何分析算法研究

夏 阳,廖 科

(大连理工大学运载工程与力学学部汽车工程学院,工业装备结构分析国家重点实验室,辽宁,大连 116024)

摘 要:梁结构在工程中应用广泛,梁结构的仿真分析是计算力学的一个重要研究内容。该文研究了复杂三维曲梁结构的等几何分析方法,首次应用拟协调有限元中的多套函数技术,使用降阶基函数逼近梁内应变项,解决复杂三维曲梁结构仿真中的闭锁问题。利用全局坐标系列式方法,避免了单元刚度阵组装时的复杂坐标变换过程,提高计算效率。使用多片 NURBS(非均匀有理 B样条)数据表示复杂三维梁结构,可精确描述曲梁结构的几何形状,与有限元方法等仿真技术相比避免了网格生成过程,减少了几何误差。数值结果表明该文算法可有效解决闭锁问题,适于复杂三维曲梁结构的仿真分析。

关键词:等几何分析;三维曲梁;闭锁问题;多套函数;拟协调有限元

梁结构在工程中应用广泛,在机械、土木行业的大规模承载结构和电子元件等微结构设计中发挥着重要作用。梁结构采用有限元方法进行仿真分析时,需要将几何模型剖分为网格模型,对曲梁结构采用以直代曲的方式生成有限元模型。这个过程中引入了几何误差,降低了计算精度。Hughes等[1-3]提出了等几何分析方法,采用几何模型直接进行仿真分析,避免了引入几何误差,提高计算精度,融合计算机辅助设计和仿真分析。等几何分析方法可以准确地描述曲梁结构的几何特征,将曲梁结构的曲率、挠率等几何本征特性融入列式。Zhang等[4]用等几何分析处理了三维空间曲梁的静力学问题。Luu和 Cottrell等[5-6]采用等几何分析,基于Timoshenko梁理论分析了自由曲面梁结构的动力学特性。Cazzani等[7]利用基于梁模型的等几何分析方法,研究建筑中拱形结构的受力状态和设计方法。这些研究表明了等几何分析在处理空间曲梁结构时具有优势。

闭锁问题[8]是梁结构有限元仿真分析中的一个经典问题。等几何梁结构分析中也存在闭锁问题。Bouclier等[9]利用 B-bar方法处理了厚梁的闭锁问题。Beirao和 Marino等[10-11]结合配点列式方法处理梁结构等几何分析的闭锁问题。Greco等[12]利用混合列式方法处理了平面 Kirchhoff梁中的膜闭锁问题。Adam等[13]改进了等几何分析中的积分策略,给出了缓解闭锁问题的积分方法。

本文研究了三维曲梁结构的无闭锁等几何分析列式方法。将拟协调有限元中的多套函数列式技术[14―16]推广到等几何分析列式中,对三维梁等几何单元的剪切和膜应变项采用低阶函数逼近(Orderreduction,OR),有效克服闭锁问题。与B-bar等方法相比,本文列式简便有效,充分考虑的曲梁结构的曲率、挠率等所有几何特性,并采用全局列式方法,适于处理复杂框架组合梁结构。算例结果表明本文列式可有效解决闭锁问题,在处理复杂框架梁结构时具有较好的精度。

1 NURBS简介

非均匀有理B样条(non-uniform rational B-splines,NURBS)是 CAD模型中描述几何采用的数学表达方式。等几何分析采用NURBS作为单元中几何和物理场的逼近基函数,可以直接利用 CAD几何数据进行仿真分析,实现CAD设计和CAE仿真模拟数据格式的统一。本节对NURBS进行简要介绍,更详细的内容请参考专著[17―18]

一条P次曲线可采用NURBS表示为:

其中:

式中:Pi为第i个控制点;表示第iP次NURBS基函数;表示第iP次B样条基函数;wi为第i个控制点对应的权值。对于给定节点向量:P阶B样条基函数定义为:

p=0时:

p≥1时,递推形式如下:

2 三维曲梁理论

2.1 坐标系

三维曲梁列式需要利用全局坐标系和基于Frenet标架定义的局部坐标系。全局坐标系如图 1中的(x,y,z)坐标系所示。Frenet标架(t,n,b)定义局部坐标系为:

图1 Frenet坐标系(t,n,b)和全局坐标系(x,y,z)
Fig.1 Frenet coordinate system(t,n,b)and global coordinate system(x,y,z)

式中:表示单位切向量;表示单位主法线向量;表示单位副法线向量;为曲线的参数表示;ξ为参数坐标;表示对参数坐标求导。由微分几何理论,曲线的Frenet-Serret方程为:

式中:表示曲线的曲率;表示挠率。

任意向量如果用Frenet标架表示为:

该向量对弧长微分可表示为如下形式,其中s表示弧长:

2.2 本构方程

列式中采用的本构方程如式(7)所示:

其中:

式中:FtFnFb表示Frenet标架下等几何单元端部的轴向力;MtMnMb表示Frenet标架下等几何单元端部的扭转力矩;E表示弹性模量;G表示剪切模量;A表示梁的截面积;k表示剪切校正因子;ItInIb表示Frenet标架下梁横截面上的惯性矩。假设梁截面存在两个对称轴,则弹性矩阵为式(7)所示的对角阵形式。

2.3 Frenet坐标系下应变方程

按照 Timoshenko梁理论[19―20],在局部坐标系下应变-位移方程为:

其中:ξt是切向应变;γnγbtntb平面内的剪切应变;χnχb表示tbtn平面内的曲率;ηt表示nb横截面上的转角。

2.4 全局坐标系下的应变方程

在式(8)中,应变由局部坐标系下的位移分量表示。在单元刚度阵组装到整体刚度阵时,需要进行局部和全局坐标系位移分量的转换,增大计算量。本节我们将局部坐标系下的应变方程推广到全局坐标系下。应用该全局坐标系下几何方程列式,可以对单元刚度阵直接进行组装,减少计算量。局部坐标系下位移和转角变量表示为:

式中,表示局部坐标系下的物理量。通过旋转矩阵将其转换成全局坐标系下的位移和转角变量:

转换公式如下所示:

将式(9)和式(10)代入式(8),并利用式(6),得到全局坐标系下应变方程为:

其中:uxuyuzxθyθzθ分别表示三维曲梁六个自由度的节点位移,可以采用NURBS基函数进行插值表示为:

式中:是第i个控制点对应的位移和转角在全局坐标系下分量;是NURBS基函数。

2.5 刚度阵组装方法

如图2所示典型多片梁结构。现有等几何分析算法,如参考文献[21]所描述算法首先计算Frenet局部坐标系下的单片刚度阵列式,总刚度矩阵的组装过程需要对单片刚度阵进行转换,组装过程如下。

假设以第一片梁的局部坐标系为整体坐标系。对于第一片梁的单元刚度阵,可以直接组装到总刚度矩阵:

第二片梁的单元刚度矩阵,需要变换到第一片梁所在的局部坐标下,即整体坐标系下,然后再组装到总刚度矩阵:

第三片梁的单元刚度矩阵的处理方法和第二片量类似,需要先变换到第一片梁所在的局部坐标下,最后组装到总刚度矩阵:

特别地,对于复杂的多片梁结构,总刚度矩阵的组装需要大量的坐标变换的计算。

式中:Ti,1表示第i片梁的单元刚度矩阵变换到第1片梁所在的局部坐标下的方向变换矩阵。本文采用2.4节式(11)所示的全局坐标列式,可直接进行单片刚度矩阵的组装。如图2所示多片梁结构,按照本文算法总刚度矩阵的组装过程如下:

图2 典型多片梁结构
Fig.2 Multi-patch beam structure

式中:表示第i片梁在全局坐标系下的单元刚度矩阵。所以,相比传统等几何分析算法下总刚度矩阵的组装方式,本文算法在全局坐标下直接计算各个片的刚度阵,总刚度矩阵的组装不需要频繁的进行坐标变换,大大提升了复杂多片梁结构的计算效率。

3 多套函数逼近的无剪切闭锁算法

3.1 多套函数列式理论

取参数域上的等几何单元为例,设表示经典的等几何分析列式中采用的应变逼近式,带有上标“-”的表示本文的多套函数列式理论采用的应变逼近式,右上角“h”表示其为离散形式。经典等几何分析中应变的逼近式由位移逼近式直接求导得到。在多套函数列式中,采用独立于位移逼近式的B样条基函数进行逼近:

式中:表示逼近式的广义自由度;是多套函数列式中采用的B样条基函数。该基函数相对位移逼近式使用的B样条基函数阶次较低,在3.3节将讨论该基函数的阶次选择策略。

采用最小二乘法,使逼近ψh

将式(18)代入式(19)得:

表示逼近残差的平方和:

其取最小值时,对每个的微分为零:

整理得到:

为保证所求有解,需要插值基函数满足线性无关。本文选取B样条基函数作为插值基函数,可以自然满足上述条件。

由加权积分方法也可以得到式(23)。将B样条基函数作为权函数,代入式(18)则得到:

通过式(23)可以解出广义自由度。将求得的代回式(18)即可得到应变矩阵:

即可求出单元刚度阵:

3.2 三维曲梁应变项处理策略

对于平面曲梁,我们对膜应变项tξ和剪切应变项γn分别进行如3.1节所述处理,很好地抑制了平面曲梁中的闭锁问题[19]。本节探讨三维曲梁列式时的应变项处理策略。图3所示为曲梁模型,弹性模量为200 GPa,泊松比0.3,旋转半径R=5 m,梁横截面半径,受集中载荷节点向量对应的权重NURBS模型控制顶点的坐标为

图3 曲梁模型图
Fig.3 Curved beam

应变项处理策略分为以下情况,情况 1:处理tξnγbγbχ四项;情况2:处理tξnγbγ三项;情况3:处理nγbγ两项。以上三种情况分别表示为情况1、情况2、情况3。采用这三种策略,计算图3模型在集中载荷作用下的B端位移。采用两个高斯积分点进行计算。这三种应变项处理策略对应结果与等几何分析的对比见图4。图中纵坐标表示仿真计算解与解析解的相对误差,横坐标表示NURBS单元的个数。

由以上测试如图4所示,对于三维曲梁应变项的处理策略,推荐使用情况2对应的处理tξnγ、bγ三项应变项。该方案具有较快的收敛率和较高的计算精度,同时也具有相对较低的计算成本。

图4 三种不同应变项处理策略的结果
Fig.4 Results of three different strain term treatment strategies

3.3 三维曲梁列式降阶策略

下文举例来说明三维曲梁列式中的降阶策略。首先考虑节点向量的选择方式。对于阶次p=3的基函数,选用如下节点向量:

对基函数进行一次降阶,p=2,那么节点向量可以为:

或者:

对应节点向量式(28)、式(29)的基函数在节点0.5处的连续性不相同。但我们的计算表明,采用这两种节点向量得到的结果是一样的。如需详细了解重控制点问题的研究可参考文献[22]。因此,为减少计算量,采用的策略是选用式(29)的节点向量进行降阶处理。

基函数降阶策略如图5所示。本文测试了以下三种形式。

情况 1:两端单元阶次为P=1,中间单元阶次P=0;情况2:两端前两个单元阶次为P=1,中间单元阶次P=0;情况3:所有单元阶次为P=1。

图5 基函数降阶策略
Fig.5 Reduction strategy of base function

为了保证降阶后基函数具有足够阶次,算法具有足够的稳定性,对于边界单元,进行降阶处理应保证阶次至少为1,即p≥1,内部单元p≥0。

这三种策略对3.2所述算例的计算结果见图6。

图6 三种不同降阶策略的结果
Fig.6 Results of three different order reduction strategies

计算中采用2个高斯积分点的积分策略。基函数阶次为2,图6(b)中采用8个等几何单元。由图6所示结果,对于三维曲梁建议采用情况1对应的两端单元为P=1,中间单元P=0的降阶策略。例如,对于p=2的B样条基函数有以下节点向量:

在该片内,第一个单元(0-0.5)对应节点向量:

第二个单元(0.5-0.6)对应节点向量:

第三个单元(06-1)对应节点向量:

对于第一个单元降阶以后,p=1,节点向量为:

对于第二个单元降阶以后,p=0,节点向量为:

对于第三个单元降阶以后,p=1,节点向量为:

4 数值算例

本文约定所有算例均为三维模型,梁的横截面均为圆形,NURBS样条基函数的阶次都为P=2。

4.1 多片曲梁

如图7所示为一个由两片NURBS模型描述的三维曲梁,采用2个高斯积分点积分策略,共4个单元;弹性模量 200 GPa,泊松比 0.3,长度L=10m ,高度H=10m,截面圆形半径r是变化量。添加集中载荷Fy=-50 N。则应用本文提出的OR方法与等几何方法对比结果分析如图8所示。图中的纵坐标表示计算解与参考解的相对误差。

图7 两片三维曲梁
Fig.7 Two-patch 3-D curved beam

由图8,采用4个等几何单元,当长细比变大时,等几何分析结果出现严重的闭锁现象。本文方法可有效地抑制多片三维曲梁的闭锁现象,始终得到高精度解。

图8 不同长细比下多片曲梁的处理结果
Fig.8 Treatment results of multi-patch curved beam

4.2 半圆形曲梁

图9所示为一个两端固定的半圆形曲梁,采用3个高斯积分点积分策略,共16个单元;弹性模量200 GPa,泊松比0.3,旋转半径R=3 m。截面半径r为变化量,受集中载荷参考数值解,等几何数值解和本文提出的降阶方法得到的计算结果见表1。两种方法的对比结果分析如图10所示。

图9 半圆形曲梁
Fig.9 Semicircle curved beam

图 10横坐标表示直梁长度和半径的比值,即长细比。由图10,等几何分析结果出现严重的闭锁现象。本文方法可有效地抑制三维半圆形曲梁的闭锁现象,始终得到高精度解。

图10 半圆形曲梁的处理结果
Fig.10 Treatment results of semicircle curved beams

表1 参考解与数值分析结果
Table 1 Reference solutions and numerical results

4.3 螺旋曲梁

图11 螺旋曲梁
Fig.11 Helix curved beam

如图 11所示为一端固定的螺旋型三维曲梁,NURBS基函数阶次P=3,采用三个高斯积分点积分策略,共 24个单元;弹性模量 200 GPa,泊松比0.3,旋转半径R=3 m,高度H=2 m,截面半径r为变化量,受集中载荷F=-1000 N。则应用OR方法与等几何方法对比结果分析如图12所示。

图12 螺旋曲梁的处理结果
Fig.12 Treatment results of helix curved beams

由图12,当长细比变大时,等几何分析结果出现较大的误差。本文方法可有效地抑制三维螺旋梁的闭锁现象,始终得到高精度解。

4.4 多片梁结构-车身

如图13所示为一多片NURBS车身模型,梁结构材料参数为:弹性模量200 GPa,泊松比0.3,车身长车身高车身宽梁截面半径模拟车身弯曲刚度试验,受集中载荷Fy=-5000 N。计算图13模型在集中载荷作用下前侧作用点的位移。参考计算解为 7.68286×10-5m,由采用极密网格的有限元方法计算得到。则应用OR方法与等几何方法对比结果分析如图14所示。

图13 车身
Fig.13 Car body

图14 车身的处理结果
Fig.14 Treatment results of car body

由图14,当单元个数较少的时候等几何分析结果出现严重的闭锁现象。本文方法可有效地抑制多片梁结构的闭锁现象,计算精度得到有效的提高。

5 结论

本文首次采用拟协调有限元中的多套函数技术,使用降阶逼近应变策略,运用于复杂三维曲梁结构仿真中。本文采用全局坐标系下的应变-几何方程,避免了全局-局部坐标系间的频繁变换,能有效地提升计算效率。本文对多套函数技术处理三维复杂梁结构的应变项处理方法、降阶策略进行了探讨,给出了最优化策略,使三维曲梁等几何算法的计算精度大大提高。实例证明,本方法能有效的克服复杂三维梁结构的闭锁问题。本方法算法理论清晰简洁,易于编程实现,可以推广到板壳结构的等几何分析中,用于构造高性能的等几何板壳分析算法。

参考文献:

[1]Hughes T J R,Cottrell J A,Bazilevs Y.Isogeometric analysis: CAD,finite elements,NURBS,exact geometry and mesh refinement [J].Computer Methods in Applied Mechanics and Engineering,2005,194(39/40/41):4135-4195.

[2]Bazilevs Y,Da Veiga LB,Cottrell JA,et al.Isogeometric analysis: approximation,stability and error estimates for h-refined meshes [J].Mathematical Models & Methods in Applied Sciences,2006,16(7): 1031-1090.

[3]吴紫俊,黄正东,左兵权,等.等几何分析研究概述[J].机械工程学报,2015(5): 114-129.Wu Zijun,Huang Zhengdong,Zuo Bingquan,et al.Perspectives on isogeometric analysis [J].Journal of Mechanical Engineering,2015(5): 114-129.(in Chinese)

[4]Zhang G,Alberdi R,Khandelwal K.Analysis of three-dimensional curved beams using isogeometric approach [J].Engineering Structures,2016,117: 560-574.

[5]Luu A T,Kim N I,Lee J.Isogeometric vibration analysis of free-form Timoshenko curved beams [J].Meccanica,2015,50(1): 169-187.

[6]Cottrell J A,Reali A,Bazilevs Y,Hughes T J R.Isogeometric analysis of structural vibrations [J].Computer Methods in Applied Mechanics and Engineering,2006,195(41/42/43): 5257-5296.

[7]Cazzani A,Malagu M,Turco E.Isogeometric analysis: A powerful numerical tool for the elastic analysis of historical masonry arches [J].Continuum Mechanics and Thermodynamics,2016,28(1/2): 139-156.

[8]王勖成.有限单元法[M].北京: 清华大学出版社,2003: 354-364.Wang Xucheng.Finite element method [M].Beijing:Tsinghua University Press,2003: 354-364.(in Chinese)[9]Bouclier R,Elguedj T,Combescure A.Locking free isogeometric formulations of curved thick beams [J].Computer Methods in Applied Mechanics and Engineering,2012,245: 144-162.

[10]Beirao da Veiga L,Lovadina C,Reali A.Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods [J].Computer Methods In Applied Mechanics And Engineering,2012,241/242/243/244: 38-51.

[11]Marino E.Locking-free isogeometric collocation formulation for three-dimensional geometrically exact shear-deformable beams with arbitrary initial curvature[J].Computer Methods in Applied Mechanics And Engineering,2017,324: 546-572.

[12]Greco L,Cuomo M,Contrafatto L,et al.An efficient blended mixed B-spline formulation for removing membrane locking in plane curved Kirchhoff rods [J].Computer Methods in Applied Mechanics and Engineering,2017,324: 476-511.

[13]Adam C,Bouabdallah S,Zarroug M,et al.Improved numerical integration for locking treatment in isogeometric structural elements,Part I: Beams [J].Computer Methods in Applied Mechanics and Engineering,2014,279: 1-28.

[14]唐立民,陈万吉,刘迎曦.有限元分析中的拟协调元[J].大连工学院学报,1980,19(2): 19-35.Tang Limin,Chen Wanji,Liu Yingxi.Quasi-conforming elements for finite element analysis [J].Journal of Dalian Institute of Technology,1980,19(2): 19-35.(in Chinese)

[15]张鸿庆,王鸣.多套函数有限元逼近与拟协调板元[J].应用数学和力学,1985,6(1): 41-52.Zhang Hongqing,Wang Ming.Finite element approximations with multiple sets of functions and quasi-conforming elements for plate bending problems[J].Applied Mathematics and Mechanics,1985,6(1):41-52.(in Chinese)

[16]胡清元,夏阳,胡平,等.假设位移拟协调平面单元应变离散算法研究[J].工程力学,2016,33(9): 30-39.Hu Qingyuan,Xia Yang,Hu Ping,et al.Research on the strain discretization algorithmin assumed displacement quasi-conforming plane finite element method [J].Engineering Mechanics,2016,33(9): 30-39.(in Chinese)

[17]Piegl L A,Tiller W.The NURBS book [M].2nd ed.,Germany: Springer Verlag,1997.

[18]过斌,葛建立,杨国来,等.三维实体结构 NURBS等几何分析[J].工程力学,2015,32(9): 42-48.Guo Bin,Ge Jianli,Yang Guolai,et al.Nurbs-based isogeometric analysis of three-dimensional solid structures [J].Engineering Mechanics,2015,32(9): 42-48.(in Chinese)

[19]Hu P,Hu Q,Xia Y.Order reduction method for locking free isogeometric analysis of Timoshenko beams [J].Computer Methods in Applied Mechanics and Engineering,2016,308: 1-22.

[20]曾森,陈少峰,曲婷,等.大位移小转角空间曲梁的弹性力学方程[J].工程力学,2010,27(12): 14-20.Zeng Sen,Chen ShaoFeng,Qu Ting,et al.Elasticity equations for spatial curved beams with large displacement and small rotation [J].Engineering Mechanics,2010,27(12): 14-20.(in Chinese)

[21]Hu Q,Xia Y,Zou R,Hu P.A global formulation for complex rod structures in isogeometric analysis [J].International Journal of Mechanical Sciences,2016,115/116: 736-745.

[22]张勇,林皋,胡志强.等几何分析方法中重控制点问题的研究与应用[J].工程力学,2013,30(2): 1-7.Zhang Yong,Lin Gao,Hu Zhiqiang.Repeated control points issue in isogeometric analysis and its application[J].Engineering Mechanics,2013,30(2): 1-7.(in Chinese)

LOCKING-FREE ISOGEOMETRIC ANALYSIS OF COMPLEX THREE-DIMENSIONAL BEAM STRUCTURES

XIA Yang,LIAO Ke
(State Key Laboratory of Structural Analysis for Industrial Equipment,School of Automotive Engineering,Dalian University of Technology,Dalian 116024,China)

Abstract:The beam structure is widely used in engineering.The numerical simulation of beam structures is an important topic in computational mechanics.In this paper,the locking-free isogeometric analysis of complex three-dimensional beam structures is investigated.The technique of multiple sets of approximation functions originated from quasi-conforming finite element method is first applied to the isogeometric analysis of three-dimensional beam structures to solve the locking problem.Order-reduced approximation functions are applied to simulate the strains of beams.Global formulation of beam strains is applied,and the stiffness matrices of beam elements and patches can be combined without transformation between local and global coordinate systems.The beam structure is described by multi-patch non-uniform rational B-spline functions.The geometry is exactly described,and the geometrical error introduced by finite element mesh can be avoided.The numerical experiments prove that the proposed algorithm can effectively avoid the locking problem in Timoshenko beam formulation,and is suitable for the analysis of complex three-dimensional beam structures.

Key words:isogeometric analysis; three-dimensional beam; locking free; multiple sets of approximation functions; quasi-conforming finite element method

作者简介:廖科(1993―),男,四川人,硕士生,主要从事等几何分析研究(E-mail: liaokk11@dlut.edu.cn).

通讯作者:夏阳(1987―),男,河南人,讲师,博士,主要从事等几何分析和增材制造工艺力学研究(E-mail: yangxia@dlut.edu.cn).

基金项目:中国自然科学基金项目(11702056,61572021);中央高校基本科研业务费专项资金项目(DUT17JC32).

收稿日期:2017-08-15;修改日期:2018-01-16

文章编号:1000-4750(2018)11-0017-09

doi:10.6052/j.issn.1000-4750.2017.08.0602

文献标志码:A

中图分类号:O241.8;TU323.3