多边形应力杂交单元的接触算法研究*

杨 锋, 郭 然

(昆明理工大学 建筑工程学院, 昆明 650000)

摘要: 针对单元较小的情况下,现有的非均匀模型难以得到精确应力场和应力集中现象的问题,采用直接约束法,提出了一种基于多边形应力杂交单元的优化接触算法多边形应力杂交单元在应力函数的构造以及积分区域划分上的优势,使其能适应复杂的模型边界与材料边界,更易划分网格根据上述理论研究成果编制了完整的计算程序,算例结果发现该方法能够得到粉末压制过程的宏观非线性力学响应、高精度的应力场以及明显的应力集中现象,为复杂优化问题的求解提供了有效手段

关 键 词: 多边形应力杂交单元; 直接约束法; 接触算法

引 言

接触问题在实际工程中广泛存在,滑动接触或滚动摩擦会引起疲劳等应力集中现象,因此接触过程的分析改进方法一直是国内外学者的研究重点如对岩质边坡稳定性的研究[1],有助于为坡体崩塌的防治提供理论基础;对凿岩爆破过程中接触、裂纹发育情况的研究[2],有助于优化地下爆破效果;对粉末冶金过程的模拟[3],可研究该过程中压制力、温度对孔隙率、密度等属性的影响

1882年,Hertz提出了理想弹性体之间相互无摩擦接触的基本假设与解析解的公式,推动了接触分析的研究此后,Mindlin[4]针对两个弹性体产生滑动接触的问题进行研究,对Hertz理论给出补充;Muskhelishvili[5]将复应力函数引入平弹性问题,给出了求解接触问题的复变函数求解法,不过仅限于弹性近似情况;Rahman等[6]提出了一种基于解析几何的直接迭代法[7],但该方法计算量大,计算机难以实现;Shi[8]提出了非连续变形分析方法(discontinuous deformation analysis),但该方法的单元应力仅由中心点应力值描述,无法得到精确的应力场和应力集中现象此外,罚函数法由于添加罚因子[9]增加了系统矩阵的带宽,传统的Lagrange乘子法具有计算规模过大的缺点,而离散单元法(discrete element method)将颗粒假设为球形等具有特定形状的刚性对象,这就忽略了物体自身的变形事实上接触过程中物体之间始终存在复杂的动态接触关系,需要充分考虑其相互作用及变形情况,很显然这些对于接触问题的计算分析方法是有缺陷的

普通位移有限元法可以根据实际结构建立有限元网格,进行实际模拟计算,但是网格必须非常细密,计算量巨大因此,1964年,Pian(卞学璜)[10]基于修正余能原理首次提出了应力杂交元的概念,现已成为杂交元法中发展最早、应用最多的一类这是一种细分的网格生成技术,以克服传统位移有限元法无法解决单元边数大于4时位移场插值困难的缺点

鉴于现有研究方法的不完善,本文推出一种基于应力杂交元的多边形应力平衡单元,即多边形应力杂交元法(polygonal hybrid stress element method,以下简称PHSEM)由于单元边数可变,PHSEM能够灵活适应复杂的模型边界与材料边界,更容易划分网格将该方法用于接触问题的分析中,能够使用较少的单元,获得高阶应力场以及接触引起的应力集中现象

1 多边形应力杂交单元法的基本理论

在许多工程问题中,由于分析模型的材料组分种类众多且空间分布复杂,常常导致常规单元的网格划分困难为了克服这种困难,本文提出一种多边形应力杂交单元,分析模型的网格划分由边数不定的多边形构成,如图1所示这样,可以根据材料的空间分布及结构的特点,方便地进行网格划分

图1 多边形应力杂交单元划分示意图图2 单个多边形应力杂交单元
Fig. 1 Schematic diagram of polygonal hybrid stress elements Fig. 2 A polygonal hybrid stress element

针对一个典型的多边形应力杂交单元,单元构成如图2所示,单元余能采用卞学璜提出的应力杂交单元的修正余能泛函[10]

(1)

在此基础上,假设多边形单元内部满足平衡的应力场[11]

σ=

(2)

通过节点的广义位移插值来假定位移场:

u=Ld

(3)

根据修正余能驻值原理,得到求解广义位移的方程组:

(4)

其中K为刚度矩阵,具体推导过程见文献[11]

由于接触问题中材料和边界的复杂性, 构造单元时会导致多边形应力杂交单元的边数 n>4而位移有限元法在处理n>4的单元时,会遇到难以构造n次插值多项式的问题,进而无法保证位移场在单元内部以及边界上协调一致,并且易导致刚度矩阵的奇异且普通位移有限元法往往需要很密的网格来描述应力场,这就需要很丰富的网格划分经验,同时也增加了计算量而PHSEM在处理n>4的单元时,位移插值仅在单元边界上处理完成以满足边界位移的协调连续,使得单元边数的增加不会阻碍插值函数的构造,能够适应复杂的材料与模型边界而自由划分网格,更适用于复杂的接触问题

2 弹性接触问题的基本理论

鉴于求解接触问题的传统数值方法的不足,本文采用直接约束法,建立主从节点关系进行接触搜索与判断

2.1 直接约束法的基本原理和求解过程

直接约束法的基本思想是定义一个物体为接触体,另一个物体为目标体,接触力发生在接触体节点与目标体表面之间[12]为了将从接触体(e2)和主接触体(e1)的单元结合在一起,需要在单元刚度矩阵中消除从节点自由度,这样能够较好地引导接触节点在目标表面上的位移来满足接触条件

图3 基于主从自由度的接触关系
Fig. 3 The contact relationship based on the master/slave freedom

以图3所示的两个四节点单元为例,其节点位移向量为

δe=Ui Vi Uj Vj Uk Vk Ul Vl Um Vm

Un Vn Up Vp Uq VqT

(5)

e2上的节点m(从节点)与e1ij(主节点)边发生接触,则m点位移可通过ij点位移线性插值得到[13]

(6)

其中U表示水平位移,V表示竖直位移im长度为l1jm长度为l2,2l=l1+l2由于UmVm并不独立,引入新的节点位移向量:

(δ′)e=Ui Vi Uj Vj Uk Vk Ul Vl

Un Vn Up Vp Uq VqT

(7)

δe=e,其中T为16×14主从自由度约束矩阵:

(8)

将其引入刚度阵K,则

K′=TTKT

(9)

其中K′为对称正定矩阵,这便是引入主从自由度后刚度阵的组装主从自由度法即是在刚度矩阵组装时将以上主从约束关系引入,从而达到把主从节点“耦合”在一起的目的[13-14]其关键是实时跟踪可能发生接触的几何体,用该方法判断接触主要通过空间的几何构形,这样在检测到接触发生时便可直接通过施加动态约束来传递相应的自由度和节点力,在求解涉及大变形的接触问题方面有很大优势[13]且该方法对接触描述的精度高,不增加系统的自由度,更适用于静态非线性问题的求解,不需要定义特殊的界面单元,也不涉及复杂的接触条件变化,因而具有很强的适应性

2.2 接触搜索与判断

数值计算接触过程中难以精确描述节点恰好在一个接触段,因此用接触段上的接触距离容限来解决这个问题,如图4,其中dist表示接触容限,B表示偏移系数(bias factor)一般将穿透归于正常的接触,如图5

图4 接触面附近的接触距离容限 图5 穿透
Fig. 4 The contact distance tolerance near the contact surface Fig. 5 Penetration

图6 法向-切向接触条件 图7 主从节点对
Fig. 6 The normal-tangential contact Fig. 7 The master-slave pairs

实际的数值模拟中由于各节点的力平衡方程不能精确满足,可能在发生分离时节点上仍有一个小的正反力存在,因此需设置一个门槛值以避免发生不必要的分离

如图6,以e2(从接触体)为受力体,则切向作用力Pτ的方向应与e2相对于e1(主接触体)的运动方向相反,设接触体摩擦因数为μ则当e2有正向相对运动趋势(静摩擦状态)或发生正向相对滑动时满足Pτμ|Pn|;有反向相对运动趋势(静摩擦状态)或发生反向相对滑动时满足Pτ≤-μ|Pn|;无相对运动趋势(静摩擦状态)或不发生相对滑动时满足-μ|Pn|≤Pτμ

综合以上接触条件,建立的主从节点关系如图7应用以上算法,采用FORTRAN进行编程,流程图如图8

(a) 分离判断
(a) Judgment of separation

(b) 滑动判断 (c) 接触搜索
(b) Judgment of sliding (c) Judgment of contact
图8 接触算法流程图
Fig. 8 The flow chart of the contact algorithm

3 算例和讨论

3.1 算例1

有一平弹性冲头,其上作用压力为fy=9.8 MPa=9.8 N/mm2,已知弹性模量EA=EB=206 GPa=2.06×105 N/mm2,Poisson比νA=νB=0.3,其余尺寸如图9所示[16]

图9 平弹性基础与平弹性冲头之间的接触对
Fig. 9 Contact pairs between the flat elastic foundation and the flat elastic punch

如图10、11,采用PHSEM划分共469个8节点单元,1 523个节点,冲头边缘处最小单元边长为5 mm;应用MARC划分共31 948个8节点单元,32 286个节点,冲头边缘处最小单元边长为1 mm显然PHSEM所用单元数远远少于普通有限元方法

图10 PHSEM网格划分示意图 图11 MARC网格划分示意图
Fig. 10 Schematic diagram of mesh generation with PHSEMFig. 11 Schematic diagram of mesh generation with MARC

图12 PHSEM接触应力图 图13 MARC接触应力图
Fig. 12 Contact stress diagram with PHSEM Fig. 13 Contact stress diagram with MARC

为了解释图中的颜色,读者可以参考本文的电子网页版本,后同

该模型为对称结构,现只观察右半侧模型压力分布情况从图12、13中可以发现PHSEM得到的整体应力分布对称性好,应力集中现象明显,可以认为与MARC的计算结果一致图14中的压力分布说明两种方法的接触应力大小和分布拟合较好,且符合客观规律,即接触压力在中间最小,保持平稳向边缘过渡,接近边缘处时开始急剧增大,理论上来说弹性情况中边缘处应力具有奇异性,网格划分越密压力越大,最终可趋于无限大图12为本文所用PHSEM的计算结果,与图13的MARC计算结果为对比,并与文献[16]中的边界元法计算结果(此处篇幅所限,仅引用)参照对比现通过与MARC两种方法计算得到的结果与边界元法计算结果[16]对比可知,PHSEM的计算结果是有效且正确的,且能明显反映出在使用少量单元得到高阶应力场和明显应力集中现象上的优越性

3.2 算例2

先将粉末压制过程简化,如图15所示,取3个不规则块体验证准确性,其中A受正压之后与B、C接触,弹性模量EA=EB=206 GPa=2.06×105 N/mm2,Poisson比νA=νB=0.3

如图16、17,采用PHSEM划分共130个单元,471个节点;MARC划分了31 745个单元,共32 615个节点PHSEM仅在两个可能接触区域进行了网格加密操作,剩余部分仅用一个单元,这便是PHSEM能够根据复杂的接触情况灵活形成单元的优点,从而大大减少了计算量

图14 平弹性基础与平弹性冲头压力分布
Fig. 14 The pressure distribution between the flat elastic foundation and the flat elastic punch

图15 不规则块体之间的接触对
Fig. 15 Contact pairs between irregular blocks

图16 PHSEM网格划分示意图图17 MARC网格划分示意图
Fig. 16 Schematic diagram of mesh generation with PHSEMFig. 17 Schematic diagram of mesh generation with MARC

该模型简化了粉末压制过程,研究其中3个颗粒互相挤压的情况,同样属于弹性问题因而接触区域存在应力奇异性图18~21分别表示PHSEM接触应力图以及MARC接触应力图,采用σy作为观察指标从图18、20中可以发现PHSEM得到的整体应力分布明确并符合客观规律,且块体棱角处可观察到明显的应力集中情况实际上单元内应力函数是连续的,但此处PHSEM在画图时取值点较分散,因此图20中应力等值线不够光滑,且单元间应力场不连续,而PHSEM并未进行应力磨平操作此外,该模型中PHSEM对接触中心区域进行网格加密操作以观察其应力集中情况,而对远离核心观察区域的部分仅用一个单元描述,继续加密能够得到更好的应力场,但是也会提高计算成本从图19、21中可以发现虽然MARC也出现了明显的应力集中情况,但接触核心区的应力界限不规则现通过与MARC的结果对比可知,本研究的计算结果是有效且正确的,为模拟粉末压制过程提供了基础

图18 PHSEM接触应力图(左) 图19 MARC接触应力图(左)
Fig. 18 Contact stress contour with PHSEM(left)Fig. 19 Contact stress contour with MARC (left)

图20 PHSEM接触应力图(右) 图21 MARC接触应力图(右)
Fig. 20 Contact stress contour with PHSEM(right)Fig. 21 Contact stress contour with MARC (right)

3.3 多体粉末压制

粉末成形是粉末冶金工艺的基本工序之一,其研究重点包括压制后的硬度、变形率的计算、致密化过程及各点的相对密度等现采用YF06喷雾制粒料的粉末电镜图放大100倍[17],取20 mm×20 mm正方形区域内的颗粒作为研究对象(图22)压制方式为双向压制,每一侧单步给定位移0.02 mm,弹性模量EA=206 GPa=2.06×105 N/mm2,Poisson比νA=0.3

如图23,每个粉体颗粒为一个多边形应力杂交单元,共24个在可能接触区域可以自由加密节点,使应力集中现象更明显,从而提高计算精度,充分体现了该单元在自由划分网格上的优势

该模型模拟粉末压制过程,分析了24个颗粒互相接触的情况,PHSEM得到的整体应力分布图(图24),图中各颗粒受挤压而发生明显的变形图25中横轴表示压制过程中的载荷步数,图26中横轴表示压制反力随压制过程的增加情况;由于该模型为二维平面问题,因此纵轴变形率=加载后面积/原面积

图25中加载前期由于均匀加载且颗粒间孔隙较大,变形率(变形后/变形前)变化大致为线性;图26中压制反力随系统变形而增大,后期反力增大变快是由于颗粒间孔隙减少,系统发生等量变形需要更大的压力;总体接触关系搜索正常,接触区域应力集中现象明显,结果符合客观规律,继续增加每个颗粒的节点数能够继续提高计算精度,不过也会增大计算量以上结果说明PHSEM研究粉末压制过程行为是有效的,相对于普通有限元方法网格及其细密的情况有相当大的优越性

图22 粉末扫描电镜图图23 PHSEM网格划分示意图
Fig. 22 The scanning electron microscope Fig. 23 Schematic diagram of mesh image of the powder generation with PHSEM

图24 PHSEM接触变形-σy应力图
Fig. 24 Contact deformation -σy stress contour(PHSEM)

图25 加载前期接触变形率-载荷步变化关系图26 压制反力-变形率变化关系
Fig. 25 The relation between the contact deformation rate and Fig. 26 The relation between the compressive the incremental step in the early loading period counterforce and the deformation rate

4 结 论

本研究引入了多边形杂交单元替代普通单元以模拟接触关系,建立了非连续体接触的高效模型将计算结果与边界元法的结果[16]进行对比,结果显示符合客观规律,PHSEM得到的结果是有效且正确的定性分析了多边形杂交单元法相对于普通单元在接触分析的精度、以及该接触算法在接触搜索的计算效率上的特点,证明了PHSEM具有网格划分灵活,能够使用较少的单元获得高阶应力场以及应力集中现象的优势

参考文献

[1] 王曼灵, 王爱涛, 李阳, 等. 基于FEMLIP的全风化边坡失稳破坏全过程的数值模拟研究[J]. 应用数学和力学, 2019, 40(3): 269-281.(WANG Manling, WANG Aitao, LI Yang, et al. Numerical simulation of the whole instability and destruction process for fully weathered slopes based on the FEMLIP[J]. Applied Mathematics and Mechanics, 2019, 40(3): 269-281.(in Chinese))

[2] 江成, 甯尤军, 武鑫. 地应力条件下的凿岩爆破数值模拟[J]. 工程爆破, 2017, 23(1): 16-20.(JIANG Cheng, NING Youjun, WU Xin. Numerical simulation of drilling blasting under earth stress[J]. Engineering Blasting, 2017, 3(1): 16-20.(in Chinese))

[3] TAN Liming, HE Guoai, LIU Feng, et al. Effects of temperature and pressure of hot isostatic pressing on the grain structure of powder metallurgy superalloy[J]. Materials, 2018, 11(2): 328-338.

[4] MINDLIN R D. Compliance of elastic bodies in contact[J]. Journal of Applied Mechanics, 1949, 16(9): 259-268.

[5] MUSKHELISHVILI N I. Some Basic Problems of the Mathematical Theory of Elasticity: Fundamental Equations Plane Theory of Elasticity Torsion and Bending[M]. Springer, 1977.

[6] RAHMAN M U, ROWLANDS R E, COOK R D, et al. An iterative procedure for finite-element stress analysis of frictional contact problems[J]. Computers & Structures, 1984, 18(6): 947-954.

[7] 陈文胜, 郑宏. 岩石块体三维接触判断的侵入边法[J]. 岩石力学与工程学报, 2004, 23(4): 565-571.(CHEN Wensheng, ZHENG Hong. Detection of 3D rock block contacts by penetration edges[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(4): 565-571.(in Chinese))

[8] SHI Genhua. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block system[J]. Engineering Computations, 1992, 9: 157-168.

[9] UNDERHILL W R C, DOKAINISH M A, ORAVAS G ☞. A method for contact problems using virtual elements[J]. Computer Methods in Applied Mechanics & Engineering, 1997, 143(3): 229-247.

[10] PIAN T H H. Derivation of element stiffness matrices by assumed stress distributions[J]. AIAA Journal, 1964, 2(3): 1333-1336.

[11] 郭然. 颗粒增强复合材料界面脱层和热机疲劳的数值模拟[D]. 博士学位论文. 北京: 清华大学, 2003.(GUO Ran. Numerical simulation of interfacial delamination and thermal-mechanical fatigue of particulate reinforced composites[D]. PhD Thesis. Beijing: Tsinghua University, 2003.(in Chinese))

[12] 赵建才, 郝旺身, 姚振强. 基于直接约束法的密封条与车门之间的接触特性研究[J]. 弹性体, 2005, 15(3): 17-20.(ZHAO Jiancai, HAO Wangshen, YAO Zhenqiang. Study on contact characteristics between seal and door based on direct restriction[J]. Elastomer, 2005, 15(3): 17-20.(in Chinese))

[13] 杜丽惠, 刘先行, 刘宁. 评价高边坡断层的直接约束方法[J]. 清华大学学报(自然科学版), 2005, 45(12): 1613-1616.(DU Lihui, LIU Xianxing, LIU Ning. Direct constraint procedure forevaluating faults on steep slopes[J]. Journal of Tsinghua University(Natural Science Edition), 2005, 45(12): 1613-1616.(in Chinese))

[14] FARAHANI K, MOFID M, VAFAI A. A solution method for general contact-impact problems[J]. Computer Methods in Applied Mechanics & Engineering, 2000, 187(1): 69-77.

[15] 李小凯, 郑宏. 基于线性互补的非连续变形分析[J]. 岩土力学, 2014, 35(6): 1777-1794.(LI Xiaokai, ZHENG Hong. Discontinuous deformation analysis based on linear complementarity theory[J]. Rock and Soil Mechanics, 2014, 35(6): 1777-1794.(in Chinese))

[16] 蒲军平. 二维含缺陷弹性体移动和滚动接触的边界元法[D]. 博士学位论文. 北京: 清华大学, 2000.(PU Junping. Boundary element method for moving and rolling contact of 2D elastic bodies with defects[D]. PhD Thesis. Beijing: Tsinghua University, 2000.(in Chinese))

[17] 周书助, 彭卫珍, 杜亨全. YF06纳米硬质合金粉末压制性能的研究[J]. 硬质合金, 2004, 21(3): 138-141.(ZHOU Shuzhu, PENG Weizhen, DU Hengquan. The study of the pressing-pepformance of the peg-added and spraryed nanometer cemented carbide powderm[J]. Cemented Carbide, 2004, 21(3): 138-141.(in Chinese))

Study on Contact Algorithms for the Polygonal Hybrid Stress Element Method

YANG Feng, GUO Ran

(College of Architecture Engineering, Kunming University of Science and Technology, Kunming 650000, P.R.China)

Abstract: In the case of small elements, it is difficult to obtain accurate stress fields and stress concentration with existing non-uniform models. An optimization algorithm for the theory of the PHSEM (polygonal hybrid stress element method) was proposed with the direct restriction method. The advantages of the PHSEM in the construction of stress functions and the division of integral regions make it more suitable for complex model boundaries and material boundaries and easier to realize element meshing. The complete computing program was made according to the theoretical analysis. The results show that, the PHSEM can obtain macroscopic nonlinear mechanical responses, highly accurate stress fields and obvious stress concentration phenomena in the powder compaction process, which provides an effective means for the solution of the complex optimization problems.

Key words: polygonal hybrid stress element; direct restriction; contact algorithm

文章编号:1000-0887(2019)10-1059-12

ⓒ 应用数学和力学编委会,ISSN 1000-0887

*收稿日期: 2019-01-25; 修订日期:2019-03-08

基金项目: 国家自然科学基金(11572142)

作者简介: 杨锋(1994—),男,硕士(E-mail: yangfeng329@163.com);郭然(1968—),男,教授,博士,博士生导师(通讯作者. E-mail: prof.guo@189.cn).

中图分类号: O343.3

文献标志码:A

DOI: 10.21656/1000-0887.400046

Foundation item:The National Natural Science Foundation of China(11572142)

引用本文/Cite this paper:杨锋, 郭然. 多边形应力杂交单元的接触算法研究[J]. 应用数学和力学, 2019,40(10): 1059-1070.YANG Feng, GUO Ran. Study on contact algorithms for the polygonal hybrid stress element method[J]. Applied Mathematics and Mechanics, 2019, 40(10): 1059-1070.