接触-碰撞有限元计算算法综述

接触-碰撞有限元计算算法综述

1相关研究现状实际设计中存在着变形体之间的接触冲突（或动态接触）。接触-碰撞数值算法的研究大致起始于20世纪70年代,迄今已有很多学者开展了相关研究,发展了许多优秀算法。Bourago等基于拉格朗日框架的显式有限元是处理变形体间接触-碰撞最常用的数值方法,因此本文主要针对显式拉格朗日有限元相关的接触-碰撞算法进行总结分析。2接触界面的色散2.1处理工具中界面接触的描述变形体间接触-碰撞问题的有限元理论与基本数值方法可参考文献接触-碰撞算法可处理任意多个物体间的相互作用,多物体间的接触可归结为物体间的两两作用。不失一般性地,本文以两个物体间的接触为例给出其数学描述。如图1所示,相互接触的两个变形体A和B,t时刻的构形分别标识为Ω接触-碰撞系统除应满足通常的连续体控制方程外以及接触界面上还应满足一定的运动学条件与动力学条件即两物体间的间隙要求大于或等于0,0值对应于接触状态,g如果不考虑接触面的焊接或粘接作用,界面法向作用力p2.2接触-碰撞问题的拉格朗日乘子法离散方程与不含接触的问题相比,接触-碰撞系统仅需在接触界面上额外引入约束条件(1~3)。因此,可方便地建立接触-碰撞系统的控制方程弱形式式中σ对应于接触界面约束的罚函数法弱形式。若p即为拉格朗日乘子法的弱形式。对方程(4)进行有限元离散,可得半离散方程:式中M为质量矩阵,a中心差分法是显式有限元最常用的积分方法,此时方程(7)的时域离散可概括为式中Δt采用显式积分求解时,方程(8)右端项中唯一的未知量为t接触-碰撞问题中,接触面通常不能事前确定,并且随时间动态变化,因此接触力计算可进一步分解为接触搜索与接触约束施加两个子问题。接触搜索的目的是确定系统中哪些部位发生了接触或者哪些原已接触的部位发生了分离或滑移,接触约束施加则根据系统当前接触状态计算出界面接触力以满足接触约束条件。2.3接触界面离散模型有限元方法中,相互接触的变形体(接触体)离散为有限单元,接触界面即为变形体外表面的单元面集合或节点集合。这样,接触体间的作用转化为离散的节点间、节点与单元面或单元面间的相互作用。为便于问题描述,本文采用Hallquist等根据接触界面离散方式的不同,界面模型可分为三类。(1)点-点NTN(node-to-node)模型(2)点-面NTS(node-to-segment)模型(3)面-面STS(segment-to-segment)模型接触界面离散方式在很大程度上决定了接触界面搜索方式和接触力的计算方法。鉴于NTS模型的广泛应用,本文将针对基于该类模型发展的相关算法进行讨论。3接触搜索的实现如上文所述,接触界面的确定是接触-碰撞问题数值计算中最耗时的部分。如采用最直接的暴力搜索方法确定接触对,其计算量将为O(N为加速接触对确定,接触搜索通常分全局搜索和局部搜索两个步骤进行。全局搜索利用某些简单的准则快速排除掉系统中不可能发生接触的部位,或挑选出最可能发生接触的从点与主片对,即确定出接触测试对。局部搜索则从接触测试对中精确确定出真实发生接触的从点与主片对,即接触对。一般而言,全局搜索决定了接触计算效率;局部搜索主要影响接触计算精度,同时对计算效率也有一定影响。典型的全局搜索算法有主从面算法3.1整个搜索方法3.1.1确定接触部位主从面算法主从面算法的计算复杂度仍为O(N主从面算法对于含有少量接触,且事前可大致确定出接触部位的问题非常有效。但该算法只能处理两个表面间的接触,不能用于单一曲面接触,也不能用于含材料断裂破坏的侵蚀接触分析,甚至不能用于相对初始构形有严重变形的接触问题3.1.2节点分组优化桶排序法是目前应用最广泛的全局搜索算法,该类算法的基本思想是桶排序法通过将节点分组(同一桶中的节点视为一组)极大缩小了全局接触搜索范围,从而大幅提高计算效率。该算法不需要预先指定接触主从面,因此具有良好的适用性,可以用于多体接触对三维空间问题,经典桶排序法Heinstein等3.1.3接触域和接触搜索级域算法层域是某层级实体占据的空间,为避免每个迭代步都进行全局搜索,可将此区域适当扩展而形成扩展域。当从点位于某主点、线或片的接触域内时,认为该从点与相应的主点、线或片接触。点、线与片接触域的定义如图4所示。片的接触域由其边、外法线方向及某种特征厚度确定,线的接触域由其外法线方向及其所在的片厚度确定,点的接触域则由其所在线的外法线方向及相应片厚度确定。接触搜索时,先在较高级的层域间进行,若两个域存在重叠,则继续进行下一级层域的相交检查,直至确定出接触测试对;若两域不存在重叠,则两接触实体间不存在接触关系,接触搜索停止。级域算法利用了接触体的层级结构,在较高级实体中排除了不必要的接触检查,理论上可大幅提高接触搜索效率,这对于同级域重叠部分较小的接触问题非常有效。但对于同级实体存在较多重叠的接触问题,如单曲面接触,计算效率不如位置码算法3.1.4基于位置码的接触分析Oldenburg等定义了主片域后,位于主片域内从点的三维空间搜索利用位置码转换到一维空间中进行:接触界面空间分割成若干桶;根据节点所在的桶,对每一个节点赋予一个位置码;根据位置码与节点坐标将三维空间中的节点排序到一维数组(称为位置码向量);对每一个主片,确定出主片域占据的桶,利用二分查找算法在位置码向量中找出这些桶中的节点(潜在接触从点)。同时,位置码算法采用层级概念定义了三层接触实体,即接触面、接触片和接触点,每一个接触面定义一个位置码向量。位置码算法将三维空间的节点排序和搜索映射到一维数组中进行,并利用二分查找法加速接触从点的搜索速度,具有很高的计算效率(计算复杂度为O(Nlog位置码算法的效率有一定的网格方向依赖性,这是由桶的编码方式决定的。Diekmann等3.1.5创建接触片几何信息树形算法是新近发展的一类全局搜索方法。该方法的基本思想是,采用树形数据结构,如二叉树和四叉树,存储接触片几何信息,利用树的快速遍历方法确定出有可能参与接触的从点与主片或者从片与主片对。针对冲压成形中板料与模具间的动态接触,Bruneel等与经典桶排序算法3.1.6接触搜索算法的发展Papadopoulos等Mahadevaiah近年,针对离散粒子接触系统,发展了大量优秀的接触搜索算法,虽然不能直接用于有限元接触计算,但为有限元接触搜索算法的设计提供了重要借鉴。王福军等3.2主要使用的算法与全局搜索相比,局部搜索算法的研究相对较少,目前主要采用的方法有点面算法和小球类算法两大类。另有少量搜索算法主要是针对两类算法的不足而提出的改进方案,如内外算法、曲面构造法及静动态检测法等。3.2.1点面搜索算法的不足点面算法假定从点n式中c式中m=(c式中t为从点n式中n若l<0,表明从点对主片有侵入,则根据接触约束算法在从点与主片间施加接触约束力。当接触界面光滑且网格质量较好时,点面搜索算法非常有效,这也是该方法受到广泛采用的原因。但算法存在如下不足。(1)经典点面算法(2)以主片外法线方向作为从点的投影方向,存在投影点不唯一或无投影点等奇异性问题(3)接触面网格质量较差时,如细长比过大,可能会发生接触搜索遗漏,如图9所示。(4)由于该方法侵入检查是针对从点与主片的关系进行的,因此不能处理如图10所示的边-面和面-面接触情况。3.2.2从点部分植物的侵入量计算面向金属成形模拟,针对经典点面算法内外算法中,以从点的外法线方向作为从点对主片的投影方向(同时也是接触力施加方向)。从点外法向定义为与该节点相连的接触片法线方向的平均值,如图11所示。从点相对于主片一条边的内或外状态,通过从点与该边构成的三角形在主片上的投影来确定,如图12所示,若Δ如果从点位于主片内部,则利用式(18)计算从点对主片的侵入量,式中x为投影点,利用式(17)可以解析给出;x内外算法的优势在于,接触点可以唯一确定,避免了从点对主片的多重接触,并且接触点计算不需要迭代。因此,该算法具有很好的健壮性与较高的计算效率。由于内外算法假定接触点处主从面的法向共轴(但方向相反),当从点位于两个低阶单元的交界线上时,尤其是当网格不够精细时,内外算法给出的接触方向可能会与实际存在较大偏差3.2.3接触力的计算精度点面算法存在投影奇异性的根本原因是低阶有限元近似使得接触界面仅C构造曲面法不仅可以克服接触盲区,也有利于提高接触力计算精度,消除接触滑移时接触力的非物理震荡。关于构造曲面法,已有很多学者开展相关研究,提出了一些光滑化算法,如三次样条法FFS方法利用接触片节点坐标构造参数化光滑曲面片。如图13所示,任意一个由四条边构成的接触片,参数化光滑曲面可表示为式中x为曲面的整体坐标,u和w为曲面的局部坐标,cFFS算法构造的曲面具有C3.2.4移动点接触检查经典点面局部搜索算法(1)动态检测。根据从点和主片的当前位置及下一时间步的预测速度,利用从点和主片接触时刻的共面关系,确定出接触时刻与接触点。三维情况下,这种接触检查简化为判断一个移动三角形和一个移动点之间的接触问题(对于四边形主片,可以利用主片中心点与四个边将主片分解为四个三角形),如图14所示。求得接触时刻与接触点后,进一步判断接触时刻是否满足步长稳定性要求,及接触点是否位于主片内部,如满足则发生接触。(2)静态检测。静态检测基于预测构形进行,从点回退方向根据接触面的凹凸性质分别处理。对于凹面,接触点为从点在主片上的投影点,或者是主片顶点,从点对主片的侵彻方向如图15所示;对于凸面,侵彻方向认为沿着上一迭代步主片的法向,即图16所示的方向m^3.2.5壳单元间的接触检查与分裂材料小球算法(pinballalgorithm)是Belytschko等引入小球概念后,单元间的侵入检查变得极为简单。假定主面单元与从面单元等效小球的半径分别为R侵入深度g由式(21)计算得到,式中d=C小球算法可用于实体单元或薄壳单元的任意类型接触分析。但用于薄壳单元时,若壳元厚度相对于其特征尺寸很小,则小球算法精度会非常差;若壳单元间存在初始接触,小球算法可能会失效。为此,Belytschko等与小球算法类似,分裂小球算法以小球代替每一个参与接触的表面单元,但小球的半径要足够大以完全覆盖相应单元,该小球称为父球。父球间若有重叠区则表示两单元间可能存在接触,此时父球分裂为下一级的子球进一步作重叠区检查,直至最后一级子球的半径与壳厚度相当。四边形壳元的小球层级如图18所示,接触检查与小球逐级分裂的一维示意如图19所示。小球算法和分裂小球算法以统一的方式自动处理边-边、边-面及面-体等接触类型,算法便于向量化和并行计算;单元间的接触检查与侵入量计算简单高效。小球算法和分裂小球算法适用于碰撞和侵彻等问题的模拟。由于算法未考虑接触界面间的摩擦效应,不能应用于滑移和摩擦为关键因素的接触分析中。另外,当接触界面单元特征尺寸差别较大时,接触力的计算精度较差,不准确的几何近似和小球体积不变假设限制了该算法的应用4接触限制算法接触约束算法,也称为接触力算法,基本方法有拉氏乘子法和罚函数法4.1接触力的计算罚函数法的基本原理是,如果从点对主片没有侵入则不做处理;如果从点侵入主片,则在从点与主片间引入法向接触力,将从点推回到变形后的接触主片上,以满足不可侵入条件。罚函数法不增加系统方程的自由度数,接触力计算无需联立求解方程组,程序实现简单且计算效率高,广泛应用于显式有限元程序中,如LS-Dyna通用程序一般会对罚因子进行规则化处理,以便用户选取相对合理的罚参数。以LS-Dyna式中l为从点对主片的侵入深度式中n规则化后的接触罚参数f式中ue7884.2显式有限元接触力的计算拉氏乘子法中,以接触对间的接触力作为乘子来限制接触体间的相互侵入,不可侵入条件可精确满足。但该方法引入了新的未知变量,增加了方程组的自由度数;接触力计算需要联立方程求解,这与显式有限元计算不相容。为了在显式有限元中应用乘子法进行接触力计算,必须对经典的拉氏乘子法做必要的格式修正。Carpenter等防御节点法的基本思想是,将从点与主片之间的接触作用转化为从点与防御节点之间的接触,防御节点携带了主片的运动学与动力学信息,如位移、速度及节点力等。防御节点位于接触点,表征主片的运动,因此其速度和加速度可以利用主片的节点速度和加速度描述,式中φ为插值形函数。防御节点的运动为F为除接触力以外的力,f为从点与防御节点间的接触力。计算f时,施加不可侵入条件。得到f后,接触力就可以分布到主点上。5接触界面的全局动态性接触-碰撞问题通常采用显式有限元进行计算。显式虽然避免了内存的庞大需求,但积分稳定性要求时间步长必须足够小,这导致接触-碰撞问题求解通常需要大量迭代,对于大规模或复杂几何构形问题非常耗时,因此并行计算是接触-碰撞数值分析的必然选择。区域分解(spatialdomaindecompositionmethod)是实现接触-碰撞问题大规模并行计算的最有效方式。区域分解的基本思想与不含接触的动力学问题或准静态接触问题不同,接触-碰撞问题的接触域通常具有全局性与动态性,计算节点之间的通信关系在系统求解之前不能确定。如图20所示的一个简单接触问题,开始时刻下方物体上表面的左前部位与上面物体发生接触,而在某个时刻下方物体上表面的右后部位与上面物体发生接触。接触界面的全局性和动态性对接触-碰撞问题的区域分割和计算节点间的通信关系设计带来极大挑战。(1)对有限元计算(主要是内力计算)高效的某种区域分解在接触计算中往往不能取得好的并行性能,通常需要针对接触而引入另外的区域分割。(2)动态接触的全局性导致计算节点间存在大量的数据交互,最大限度地减小通信量才能保证算法的扩展性与并行效率。(3)接触对动态变化,局部接触搜索和约束施加在不同时刻可能由不同的计算节点完成,极易造成各节点机负载的不均衡。基于区域分解的并行接触算法研究始于20世纪90年代。Malone等上述算法仅适用于少量CPU核的集群环境,有良好扩展性、适用于大型并行系统(数千CPU核)的并行接触算法研究还较少。目前该方面的研究主要限于美国Sandia国家实验室基于区域分解实现接触-碰撞问题的粗粒度并行有单一区域分解法和双重区域分解法两种并行方案。双重区域分解法中,有限元计算采用一种子域剖分(称为主分区),接触计算采用另外一套分区(称为次分区);而单一区域分解法,接触区域不作单独的子区分割,接触的并行计算基于主分区进行。通常,双重区域分解法可获得更好的负载均衡,但节点机间的通信复杂;单一分区法通信拓扑简单,但易导致负载不均衡。5.1静态区域剖分法设计单一区域分解方法中,接触计算的并行基于主分区进行。接触算法并行的关键在于,每个计算节点可有效得到本地计算所需的全部信息。因此,该类算法通常引入影像区(ghost)或接触域(contactdomain)存储本地接触计算需要但驻留在其他节点机上的信息。基于单一区域分解的并行接触算法设计相对简单,本文以Malone等(1)不考虑接触界面,对求解域作区域剖分,分配给各计算节点。所有计算节点作如下相同的计算。(2)对任意计算节点PA,构建覆盖该计算节点单元几何空间的包围盒(boundingbox),并传递到其他计算节点,同时接收其他计算节点的包围盒信息。(3)利用体积检查(volumechecking)判断PA包围盒与其他包围盒是否有重叠。若无重叠,则两个计算节点的单元间无潜在接触;若有重叠部分,则将与PA包围盒有重叠的计算节点上的单元信息收集到PA的接触域,同时将PA单元传递给对方计算节点的接触域。(4)对每个计算节点利用串行接触算法进行接触计算(全局搜索、局部接触和约束施加)。基于单一静态区域剖分的并行接触算法概念简单,程序实现容易。由于引入影像区存储本地计算所需的全部信息,不需要每步计算都进行数据通信。但由于区域剖分时未考虑接触计算的影响,容易导致计算节点间严重的负载不均衡,算法并行扩展性差5.2分区区域分割方法双重区域分解方法中,单元计算和接触计算分别基于主分区和次分区进行。主分区通常采用图剖分法