基于边界面法的完整实体应力分析理论与应用

摘 要:提出基于边界面法(Boundary Face Method,BFM)的完整实体应力分析方法. 在该分析中,避免对结构作几何上的简化,结构的所有局部细节都按实际形状尺寸作为三维实体处理. 以边界积分方程为理论基础的BFM是完整实体应力分析的自然选择. 在该方法中,边界积分和场变量插值都在实体边界曲面的参数空间里实现. 高斯积分点的几何数据,如坐标、雅可比和外法向量都直接由曲面算得,而不是通过单元插值近似获得,从而避免几何误差. 该方法的实现直接基于边界表征的CAD模型,可做到与CAD软件的无缝连接. 线弹性问题的应用实例表明,该方法可以简单有效地模拟具有细小特征的复杂结构,并且计算结果的应力精度比边界元法(Boundary Element Method,BEM)和有限元法(Finite Element Method,FEM)高.

关键词:完整实体应力分析; 边界积分方程; 边界面法; CAD模型

中图分类号:O241.4; TB115 文献标志码:

Theory and application of stress analysis on

complete solids based on boundary face method

ZHANG Jianming

(State Key Lab. of Advanced Design & Manufacturing for Vehicle Body,

Hunan Univ., Changsha 410082, China)

Abstract: The method of stress analysis on complete solids based on Boundary Face Method(BFM) is proposed. Without any geometric simplification of a structure, the analysis takes into account all local and small features of the structure according to their real shapes and sizes. BFM which is based on boundary integral equations is a natural choice for implementation of the complete solid stress analysis. In the BFM, both boundary integration and interpolation of field variables are performed in the parametric space of each boundary face. The geometric data at Gaussian integration points, such as the coordinates, the Jacobians and the out normals are calculated directly from the faces rather than from element interpolation, and thus the geometric errors can be avoided. The method can seamlessly interact with CAD software, because it is implemented by making direct use of a CAD model through its boundary representation data. The applications on linear elastic problem indicate that the proposed method has more accurate precision of stress and is more effective and easier to use than the traditional Boundary Element Method(BEM) and Finite Element Method(FEM) for simulation on complex structure with small features.

Key words: complete solid stress analysis; boundary integral equation; boundary face method; CAD model

0 引 言

目前,结构应力分析以有限元法(Finite Element Method,FEM)为主体.FEM需要在结构的几何模型上另外生成1个离散的网格模型.[1]网格单元用于进行能量积分和变量插值,还用于插值逼近原始的CAD实体几何.对于复杂且含有细小特征的实体,离散成具有形状良好的实体单元比较困难.此外,FEM的实现基于所求问题控制方程和边界条件的等效积分“弱”形式,要求其试函数(一般为插值函数)至少具有1阶连续性.[1]在求解力学问题时,应力精度总比位移精度低1阶,导致应力结果精度不高,但在实际问题中更关注应力值,如应力集中部位及最大值.

图1 建筑外用电梯

标准节建筑外用电梯标准节见

图1.对于这样的刚架结构,FEM一般采用梁单元模拟,离散后的计算模型与实际结构在几何和拓扑上都存在很大差别,以致于变形结果相对较精确而应力精度很差,更不能得到焊缝应力的精确值.焊缝通常是应力集中点,也是实际中大多数结构破坏的源发点,因此设计工程师更关心焊缝应力的大小.骨骼组织高分辨率有限元分析模型见图2.

图 2 骨骼组织高分辨率有限元分析模型

① 该文已被《计算力学学报》录用,待发表,作者为覃先云、张见明和庄超.目前,流行的高分辨率有限元分析采用细小的8节点六面体单元逼近骨骼的几何体,然而无论单元划分有多细,其表面都不光滑,必然会产生虚假的应力集中.随着计算机硬件的快速发展,计算机的计算能力不断提高,可直接在结构的几何模型上进行应力分析,本文称这种分析为完整实体应力分析.完整实体应力分析避免使用抽象的一维、二维单元,避免对结构作几何上的简化;所有结构(包括细长和薄型结构)从整体到局部细节,如机械结构中的倒角、焊缝、退刀槽和小圆孔等,都按照实际形状尺寸作三维实体处理.

利用HUGHES等[2]和BAZILEVS等[3]提出的等几何分析方法可达到完整实体应力分析的目标.但该方法的实现仍基于FEM的基本理论,需对整个实体域在三维参数空间中离散.[2,3]对于任意复杂的实体,这种参数化离散相当困难,至今该方法只运用在由非裁剪曲面组成的简单实体上.[2,3]本文利用由ZHANG等[4]、QIN等[5]、覃先云[6]及文献“基于参数曲面三维势问题的边界面法”①提出的边界面法(Boundary Face Method, BFM)实现完整实体应力分析的目标.BFM以成熟的边界积分方程理论为基础,是对边界元法(Boundary Element Method, BEM)[7,10]和边界点法(Boundary Node Method, BNM)[11,12]的继承和发展.BFM只需实体边界的参数曲面信息,边界积分和变量插值直接在曲面的二维参数空间中实现;具有与实体CAD造型系统无缝连接的天然优势(这是因为在所有实体造型系统中都有现成的边界表征数据结构);另外,BFM与BEM一样,只需分析实体的边界离散,具有精度高,便于处理无限域问题、边界及界面问题和奇异性问题等优点.BEM在工程数值分析领域有大量的应用[7,10],为BFM的工程应用打下良好基础.因此,BFM或其他边界类型方法[11,13,14]是实现完整实体应力分析的合适选择之一.

1 理论基础――BFM

BFM是在传统BEM和BNM基础上发展起来的新数值方法.该方法继承以边界积分方程为基础的边界类型数值方法[7,10]的优势,同时克服传统方法在分析模型等方面存在的缺点.利用BFM,数值求解的应力与位移具有相同的精度.在分析结构应力时,BFM比FEM具有优势,因为FEM计算结果的应力精度比位移精度低1阶.[1]本文对BFM只作简单介绍,详细内容可参考文献[4,6]和文献①.

在传统的BEM或BNM中,将CAD几何模型离散为分析模型后,CAD模型的原始几何信息基本丢失.这2种方法所用的分析几何数据是基于网格单元通过Lagrange或Hermite插值方法近似得到的.[7,12]基于网格单元的几何插值引起的几何误差,会从根本上导致计算精度问题[2],甚至对有些计算起决定性影响.另外,分析模型和CAD几何模型的分离使设计和分析成为2个互相独立的过程.在自适应网格细分过程中,需反复与CAD系统交互,而每个阶段的交互很复杂.为克服上述缺点,在BFM中将分析模型与CAD几何模型融为一体,使结构分析与几何造型集成于统一框架.

BFM的实现直接基于边界表征的实体造型数据结构,实体的边界曲面都以参数形式表达.在其实现中,不论是对边界的数值积分还是对场变量的插值都在边界曲面的二维参数空间中实现.[4,6]该方法需将每个参数曲面离散成若干个参数曲面单元;每个单元定义在所在曲面的二维参数空间上,而非三维物理空间上;这种曲面单元相当于分片曲面(surface patch),保持曲面的原始几何信息;在数值积分过程中,被积函数的几何变量,如高斯积分点的坐标、雅可比和外法向量直接通过曲面单元中的曲面参数变量计算获得,而不通过分段多项式插值近似得到.[4,6]组成圆环体曲面在参数空间的离散见图3.曲面参数单元映射到物理空间后仍保持几何精确性:每个单元在物理空间中是光滑的,且单元间也是光滑连续过渡的;这使得求解域的整个表面光滑,避免几何误差,有利于提高数值计算精度.

(a)圆环体曲面在参数空间中的

离散 (b)参数单元映射到物理

空间

图 3 组成圆环体曲面在参数空间的离散

在曲面参数空间中离散,不但使参数网格保持精确的几何信息,而且可有效离散较复杂的曲面模拟复杂结构.复合材料的代表性体元模型见图4.该结构为ZHANG等[15,16]在进行纳米复合材料研究时的代表模型;所有复杂纳米体元的边界直接在其对应的曲面参数空间中离散,精确、有效地模拟这种复杂结构.对于这种复杂结构,在FEM中很难得到合适的体网格离散.同样,BEM也需高质量的边界单元模拟复杂边界.

图 4 复合材料的代表性体元模型

在曲面参数空间里,可以利用移动最小二乘法(Moving Least,Square method, MLS)[17]、非均匀有理B样条(NURBS)[18]和分段多项式(见文献[4,6]和文献①)等方法对未知变量进行插值逼近.如果利用MLS等无网格法中的变量逼近方法时,BFM具有像BNM和杂交边界点法(HdBNM)等无网格法[13, 17]的特点,那么曲面单元仅作为背景积分单元.如果利用分段多项式插值,那么BFM具有基于网格插值BEM的特点,可更有效地处理具有任意裁剪和细小特征结构的实体问题.在该法的程序实现框架中,建立统一的实现不同插值方法的数据结构,从而根据不同的结构采用合适的插值方式,使BFM的实现灵活自由.如图5(a)所示的肘形弯管[4],其表面的边界变量利用MLS逼近,图中的曲面参数单元作为背景积分单元,每个背景单元里的中心点为MLS插值点.在图5(b)中,利用分段多项式插值具有细小孔平面上的未知变量,在曲面的边界上采用非连续单元,而在内部利用连续单元.[5,6]这种单元的分布方式有利于灵活离散实体边界,可采用不同密度的网格有效地离散复杂结构.[5,6]在离散的边界积分方程[7,8]中,不要求插值形函数的连续性,因此在BFM采用非连续单元可满足理论要求,而在FEM中不可行.

(a) 利用MLS逼近

边界变量(b) 利用分段多项式

插值边界变量图 5 不同曲面的边界变量插值方案

直接在曲面的参数空间内进行边界积分和变量插值,直接利用CAD造型系统中参数曲面的几何信息是BFM与当今主流CAE软件中分析方法的重要不同之处.在FEM和传统BEM中,函数的插值和数值积分(能量积分或边界积分)都在单元内实现,且必须依赖于单元.在本文方法中需要的分析计算的几何变量直接来自CAD几何模型,因而与实体造型系统自然地融为一体,实现完整实体应力分析.

2 应用实例与讨论

以分析圆环体和三通体的线弹性问题为例,评估基于BFM进行完整实体应力分析方法的可行性及优势.通过分析圆环体,比较BFM与BEM的计算精度及效率.分析三通体的混合边界值问题,比较BFM与FEM的计算精度和模拟复杂结构的能力.

为估计误差和研究算法的收敛性,定义L2范数相对误差Иe=1v

max[KF(]1NNi=1(v

(e)i-v

(n)i)2[KF)](1)И[WTBX]式中:v

max为N个样本点的最大值;v

(e)i和v

(n)i分别为精确解(或参考解)和数值解.

2.1 圆环体Dirichlet边界值问题

以大圆半径为10,截面圆半径为3,大圆圆心在原点的圆环体Dirichlet问题为例,分析比较BFM与传统BEM的求解精度.材料的弹性模量假定为

1.0,泊松比为0.25,不计重力.为便于数值计算结果与精确解进行比较,考虑下面的位移解析场:Иux=-2x2+3y2+3z2

uy=3x2-2y2+3z2

uz=3x2+3y2-2z2(2)ИЦ据弹性力学的理论导出该位移场的应力解И[JB({]σ

xx=-16(3x+y+z)

σ

yy=-16(x+3y+z)

σ

zz=-16(x+y+3z)[JB)], [JB({]σ

xy=24(x+y)

σ

xz=24(x+z)

σ

yz=24(y+z)[JB)][JY](3)ИХ直鹩BFM和BEM进行分析计算,整个边界表面的位移按式(2)施加,最后将得到的应力结果与式(3)比较.

在同一程序框架中,分别用BFM和BEM分析表1中的9种离散方案.每次给圆环体的表面施加对于式(2)的本质边界条件.利用线性三角形单元插值边界变量将BFM单元定义在曲面的参数空间中,而将BEM单元定义在三维物理空间.图6为用于BFM分析的1种离散方案,其曲面单元和节点数分别为528和306.利用式(1)计算各种情况下分析结果的节点应力分量的相对误差,其结果见表1,其中:符号Tx,Ty和Tz分别为沿坐标x,y和z方向的节点应力分量的相对误差;t

Mat为获得方程系数矩阵所需的时间.

表 1 BFM与BEM的计算结果比较单元数量/个节点数量/个BFMBEMTx/%Ty/%Tz/%t

Mat/sTx/%Ty/%Tz/%t

Mat/s1121542023125287761 3402 1362 920711001271893024367381 1531 5625.0544.4803.9143.0312.1391.6811.1540.8470.691[BH]4.7304.2233.4922.8071.9901.5761.1040.8090.691[BH]9.9709.1356.3015.4043.342.5961.6751.1750.924[BH]

图 6 圆环体边界离散

由表1可知,随着单元数量的增加,这2种方法的计算结果精度都有所提高;但对于相同的离散方案,BFM的相对误差总小于BEM,特别是在网格较稀疏的情况下表现得更为明显.对于每种离散方案,2种方法获得方程系数矩阵所需的时间相当,说明BFM直接从参数曲面计算几何变量与BEM利用单元插值近似几何的计算效率相当.2种方法计算结果中沿x方向节点应力分量的相对误差见图7,可知在网格相对稀疏时BEM数值结果对网格密度的敏感性大,而BFM相对较小,说明BFM数值结果比BEM的稳定.图 7 2种方法计算结果中沿x方向节点

应力分量的相对误差

BFM具有上述优点,是因为BFM所需的几何信息直接通过参数曲面获得,而BEM则用多项式在单元内进行插值近似获得.这样,BFM的分析模型避免几何误差,而BEM在网格较稀疏时具有明显的几何误差.几何误差对分析精度有着重要影响,甚至起决定性影响.

2.2 三通体混合边界值问题

利用BFM分析具有细小倒圆角的三通体的混图 8 三通体几何结构和边界条件合边界值问题,并与FEM比较.该实体的CAD几何结构和边界条件见图8.在面faceA上约束沿坐标x,y和z方向的位移,在面faceB上施加大小为500

N/mm2均布拉力.材料的弹性模量为100000MPa,泊松比为0.25,不计重力.

利用BFM和FEM根据不同的离散方案计算该问题,并比较11个给定计算点的von Mises应力精度.计算点沿圆周方向均匀分布在圆柱面上,见图8.利用MSC Patran和MSC Nastran完成有限元分析.[19]每种离散方案所用的单元和节点数量见表2,BFM和FEM分别采用线性曲面单元和线性四面体单元.其中的2种离散方案为:图9(a)为BFM边界离散模型,其边界单元和节点数分别为2 149个和

1 647个;图9(b)为有限元法体离散模型,其体单元和节点数分别为111 026个和2 661个.

表 2 计算点von Mises应力相对误差比较计算方法[]单元数量/个节点数量/个相对误差/%BFMFEM1 175 7834.941 403 1 0133.732 149 1 6472.064 199 1 1507.3811 026 2 6615.8436 110 7 8594.3977 289(2次)116 9910

(a) 边界面法(b) 有限元法图 9 三通体的离散模型

为便于分析计算结果精度,以利用精细离散的有限元模型数值计算结果作为参考解(假定的精确解),所用的2次10节点四面体单元和节点数分别

为77 289个和116 991个.利用式(1)估算每种情况下给定计算点的von Mises应力相对误差,结果见表2,可见BFM用较少的单元可获得比FEM更好的计算结果.BFM用1 403个单元时计算点的应力相对误差只为2.46%,而FEM要用36 110个单元(误差为2.42%)才能获得同等的精度.

利用不同方法模拟结果的计算点von Mises应力值见图10,并与参考值进行比较,其中横坐标参数u∈[0,2π]表示圆柱参数曲面的参数坐标值,沿圆柱的圆周方向.由图10可知,BFM用较少的单元可获得比FEM更好的计算结果.更重要的是,在应力较大的地方(u=1.0～2.0)FEM计算结果偏离参考值较大,而BFM利用较少节点的计算结果偏离参考值相对较小,说明在利用相同阶次插值函数时,BFM比FEM在模拟应力集中方面更具有优势.

图 10 von Mises应力值比较

BFM具有较好的精度,主要是因为该方法与BEM一样,其计算结果的位移与应力具有相同的精度.[7,10]在FEM中,首先计算出节点位移,然后根据位移求出应力,这导致应力精度比位移精度低.[1]另一主要原因是BFM的计算几何信息直接来源于CAD几何模型,避免几何误差,而FEM用插值获得近似模型,具有几何误差.图11(a)和11(b)分别为图9中BFM和FEM模型中倒圆角处的局部网格.可知,图11(a)曲面网格较精确地模拟细小圆角,而图11(b)利用体网格近似模型具有明显的几何误差.BFM的网格定义在曲面的参数空间中,可很方便、精确地模拟圆角等细小几何特征.FEM的网格定义在三维物理空间中,需要较精密的网格才能模拟细小几何特征.用精密的网格会明显增加有限元的计算规模,通常将细小特征忽略,而细小特征往往是应力集中的地方.在该实例中的倒圆角处就是应力集中部位,不能忽略.

(a) 边界面法(b) 有限元法图 11 倒圆角处局部网格

由于BFM具有较精确模拟细小特征的优势,对具有细小特征的结构可直接基于三维实体弹性力学理论进行结构应力分析.某焊接刚架结构及局部焊接特征见图12,在FEM中采用抽象的单元模拟刚架结构和简化焊接连接部位,这导致应力分析结果不可靠.BFM可较精确地模拟焊接连接部位及刚架的CAD结构,结合该方法基于实体弹性理论分析该类刚架结构具有重要意义.该焊接应力分析问题为多域问题并且计算规模较大,需进一步与多域算法、快速算法[14,16, 20,21]相结合来实现完整实体应力分析.图 12 某焊接刚架结构及局部焊接特征

3 结 论

提出基于BFM的完整实体应力分析方法.在该分析方法中,所有结构(包括细长和薄型结构)从整体到局部细节都按照实际形状尺寸作三维实体处理.利用以边界积分方程为理论基础的BFM实现完整实体应力分析.在BFM中只需实体边界曲面的参数离散,边界积分和变量插值都直接在曲面参数空间里实现.该方法的实现直接基于边界表征的实体造型数据结构,可做到与CAD软件的无缝连接.

应用实例说明本文方法不但比传统BEM的计算精度高,而且计算结果对网格密度的敏感性小;用较少曲面单元就能获得需用较多体单元FEM计算结果的同等应力精度;另外,可有效、精确地模拟复杂结构,在结构应力分析中具有独特的优势.

后续工作是将BFM与多域算法、快速算法结合,实现与CAD实体造型系统无缝连接,以求解任意复杂几何形状和任意材料构成的大规模复杂工程问题.

参考文献:

[1] 王勖成, 邵敏. 有限单元基本原理和数值方法[M]. 2版. 北京: 清华大学出版社, 1997: 3,5.

[2] HUGHES T J R, COTTRELL J A, BAZILEVS Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[J]. Comput Methods Appl Mech Eng, 2005, 194(39): 4135,4195.

[3] BAZILEVS Y, CALO V M, COTTRELL J A, et al. Isogeometric analysis using T,splines[J]. Comput Methods Appl Mech Eng, 2010, 199(5): 229,263.

[4] ZHANG Jianming, QIN Xianyun, HAN Xu, et al. A boundary face method for potential problems in three dimensions[J]. Int J Numer Meth Eng, 2009, 80(3): 320,337.

[5] QIN Xianyun, ZHANG Jianming, LI Guangyao, et al. An element implementation of the boundary face method for 3D potential problems[J/OL]. Eng Anal Bound Elem, 2010: 1,10 [2010,03,29]. /science?_ob=MImg&_imagekey=B6V2N,50F36BB,1,1F&_cdi=5707&_user=2932517&_pii=S0955799710001517&_orig=search&_coverDate=07%2F01%2F2010&_sk=999999999&view=c&wchp=dGLbVtz,zSkWb&md5=0679158eb5eed5b183320c79dc91cb1a&ie=/sdarticle.pdf.

[6] 覃先云. 基于参数曲面边界面法的研究[D]. 长沙: 湖南大学, 2010.

[7] 布瑞比亚 C A. 边界单元法理论和工程应用[M]. 龙述尧, 刘腾喜, 蔡松柏, 译. 长沙: 国防业工出版社, 1988: 1,6.

[8] 龙述尧. 边界单元法概论[M]. 香港: 中国科学文化出版社, 2002: 1,4.

[9] 黎在良, 王乘. 高等边界元法[M]. 北京: 科学出版社, 2008: 21,55.

[10] 姜弘道. 弹性力学问题的边界元法[M]. 北京: 中国水利水电出版社, 2008: 79,109.

[11] MUKHERJEE Y X, MUKHERJEE S. The boundary node method for potential problems[J]. Int J Numer Meth Eng, 1997, 40(5): 797,815.

[12] CHATI M K, MUKHERJEE S. The boundary node method for three,dimensional problems in potential theory[J]. Int J Numer Meth Eng, 2000, 47(9): 1523,1547.

[13] ZHANG Jianming, YAO Zhenhan, LI Hong. A hybrid boundary node method[J]. Int J Numer Meth Eng, 2002, 53(1): 751,763.

[14] ZHANG Jianming, TANAKA Masataka. Adaptive spatial decomposition in fast multipole method[J]. J Comput Phys, 2007, 226(1): 17,28.

[15] ZHANG Jianming, TANAKA Masataka. Fast HdBNM for large,scale thermal analysis of CNT,reinforced composites[J]. Comput Mech, 2008,

41(6): 777,787.

[16] ZHANG Jianming, TANAKA Masataka, MATSUMOTO Toshiro. A simplified approach for heat conduction analysis of CNT,based nano,composites[J]. Comput Methods Appl Mech Eng, 2004, 193(52): 5597,5609.

[17] ZHANG Jianming, TANAKA Masataka, MATSUMOTO Toshiro. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method[J]. Int J Numer Meth Eng, 2004, 59(9): 1147,1166.

[18] SHAW A, ROY D. NURBS,based parametric mesh,free methods[J]. Comput Methods Appl Mech Eng, 2008, 197(17): 1541,1567.

[19] MSC.Software. MSC Patran & MSC Nastran使用指南[K]. 2002: 694,699.

[20] SMAJIC J, ANDJELIC Z, BEBENDOF M. Fast BEM for eddy,current problems using H,matrices and adaptive cross approximation[J]. IEEE Transactions on Magnetics, 2007, 43(4): 1269,1272.

[21] OSTROWSKIL J, ANDJELIC Z, BEBENDOF M, et al. Fast BEM,solution of laplace problems with H,H,matrices and

ACA[J]. IEEE Transactions on Magnetics, 2006, 42(2): 627,630.