（一般力学与力学基础专业论文）接触碰撞问题的算法研究.pdf

南京航空航天大学硕士学位论文 i 摘 要 接触碰撞问题在工程系统中广泛存在，其计算涉及到几何非线性、材料非线性和边 界条件非线性，对其进行分析和研究具有重要意义。接触碰撞问题的有限元分析是解决 结构接触碰撞问题的有效方法。本文针对结构接触碰撞问题，在以下方面进行了研究： 1. 建立了适于接触碰撞问题动态分析的非线性有限元格式， 并对其求解过程中涉及 到的有限元空间离散模型、非线性几何关系方程、材料本构模型、时间积分格式等领域 进行了分析。 2. 提出了一种新的接触局部搜索算法。 基于矢量运算的思想， 首先在全局搜索中确 定接触测试对，然后在局部搜索中采用再次筛选的方法，最后对剩余的接触测试对进行 投影计算。 3. 发展了直接以接触点为防御节点的接触力计算方法。 较之使用最为广泛的罚函数 法而言， 接触约束条件在该方法中能得到充分的满足； 同时， 该方法未引入新的未知量， 因此有限元求解方程仍是解耦的，且可用显式时间积分法进行分析，提高了计算效率。 关键词：接触碰撞，有限元，非线性，接触搜索，接触力 接触碰撞问题的算法研究 ii abstract contact-impact widely exists in engineering. the behavior of contact-impact computa tion presents geometric nonlinearity, material nonlinearity and boundary condition nonli nearity. so the analysis and research of contact-impact phenomena is very important. to solve the contact-impact problems, it is an efficient way to adopt the numerical si mulation with the help of finite element method. doing some helpful contribution to t his field is right the aim of the present paper. the main work can be shown as follo wings: 1. accurate equations are formulated for nonlinear contact-impact finite element an alysis. during the process of solving equations, the following issues are discussed: th e spacial discrete model of finite element, nonlinear geometrical equations, nonlinear material constitute model, time integration algorithm, etc. 2. a new local contact search algorithm is given. in the new algorithm, before cal culating the projection, the contact test pairs are lessened again with the technique of vector calculation. 3. the defence node algorithm has been developed. here, we consider the contac t point as defence node directly. compared with the penalty function method, which i s the most widely used method, this method can meet the contact constraint condition s completely. at the same time, this method dont introduce a new unknown variable, so the finite element equations are still decoupled, and also can be analyzed using e xplicit time integration method, which can improve the computational efficiency. key words: contact-impact, finite element, nonlinearity, contact search, contact force 南京航空航天大学硕士学位论文 v 图 清 单 图 2.1 接触前后系统的构形.5 图 2.2 典型的厚壳单元 .6 图 2.3 弹塑性连续体中已屈服点上的应力增量变化.15 图 2.4 弹塑性连续体中初始屈服点上的应力增量变化.16 图 2.5 接触碰撞系统的计算流程图 .21 图 3.1 典型的立方格结构 .28 图 3.2 单元面、面域和立方格 .30 图 3.3 主击节点与单元面的矢量关系 .32 图 4.1 防御节点的建立 .36 图 4.2 主击节点 s 与靶单元面 m 在点 c 处发生接触 .37 接触碰撞问题的算法研究 vi 表 清 单 表 2.1 高斯点弹塑性状态判断 .15 南京航空航天大学硕士学位论文 vii 注 释 表 ( )e m 单元质量矩阵 ( )e r 单元节点内力矩阵 ( )e q 单元节点外力矩阵 ( )e c f 单元接触力矩阵 tl 全 lagrange 格式 l 线性应变 nl 非线性应变 s 应力 u 位移 d 弹性矩阵 e 弹性模量 0 y 材料初始屈服应力 单位体积材料密度 h 线性应变强化参数 y 屈服应力 2 j 应力偏量的第二不变量 cell n 立方格总数 摩擦系数 承诺书 本人郑重声明：所呈交的学位论文，是本人在导师指导下，独立进行 研究工作所取得的成果。尽我所知，除文中已经注明引用的内容外，本学 位论文的研究成果不包含任何他人享有著作权的内容。对本论文所涉及的 研究工作做出贡献的其它个人和集体，均已在文中以明确方式标明。 本人授权南京航空航天大学可以有权保留送交论文的复印件， 允许论文 被查阅和借阅,可以将学位论文的全部或部分内容编入有关数据库进行检 索，可以采用影印、缩印或其它复制手段保存论文。 作者签名： 日 期： 南京航空航天大学硕士学位论文 1 第一章 绪论 1.1 引言 接触碰撞问题存在于许多工程领域，如空间伸展机构在伸展过程中与障碍物发生的 接触碰撞1，对接飞船的舱门捕获和锁钩对接2，高速铁道机车的轮轨碰撞3，机械手 在抓取重物时的接触碰撞4，汽车结构耐撞性仿真分析5，金属模具的成形过程6等。 在冲击载荷作用下，结构可能发生大变形，材料可能进入塑性状态，原来没有任何联系 的两个部件或一个部件的两个部分可能接触到一起。因此，这是一个涉及几何非线性、 材料非线性、边界条件非线性的难题7。正是因为该问题的普遍性和复杂性，对该问题 的研究也显得极为重要。 有限元法作为一种高效的工程数值分析方法，在接触碰撞问题的分析过程中得到越 来越广泛的应用8,9。 1.2 接触碰撞有限元分析研究现状 从应用数学角度来看，有限单元法基本思想的提出，可以追溯到courant在1943年的 工作， 他第一次尝试应用定义在三角形区域上的分片连续函数和最小位能原理相结合来 求解st.venna扭转问题10。一些应用数学家、物理学家和工程师由于各种原因都涉足过 有限单元的概念。但是直到1960年以后，随着计算机的广泛应用，有限单元法的发展速 度才显著加快。 现代有限单元法的第一个成功的尝试， 是将刚架位移法推广应用于弹性力学平面问 题，这是turner、clough等人在分析飞机结构时于1956年得到的成果11。他们首次给出 了用三角形单元求得平面应力问题的正确格式。 三角形单元的单元特性是由弹性理论方 程来确定的，采用的是直接刚度法。他们的研究工作打开了利用计算机求解复杂平面弹 性问题的新局面。1960年，clough进一步处理了平面弹性问题，并第一次提出了“有限 单元法”的名称，使人们开始认识了有限单元法的功效12。 接触碰撞有限元分析是一个复杂的系统工程，它涉及到以下几个研究问题：(1) 连 续介质问题的变分原理；(2) 材料本构模型；(3) 有限元单元模型；(4) 时间积分算法； (5) 接触算法。其中所采用的时间积分算法和接触算法尤为关键。 对于有限元系统所形成的二阶常微分方程组，通常有隐式和显式两种时间积分算 法。尽管动态分析中隐式算法是无条件稳定的，但在大规模动态冲击接触分析中却很少 采用， 这是因为分析的模型具有很强的非线性特性， 而且撞击载荷的有效频率成分很高， 为保证隐式算法迭代过程的收敛性和具有足够的精度，时间步长必须取的足够小，因而 接触碰撞问题的算法研究 2 使计算量十分大。另外，在大型程序分析中，隐式算法占据计算机的内存很大。相反， 显式算法正好克服上述方法在冲击接触分析中存在的不足， 也就是可以解决隐式算法收 敛性小的问题，当质量矩阵采用对角阵时，无需建立和求解联立方程组，所占据的内存 空间小，计算速度快，而且稳定性准则能自动限制时间步长的大小，保证了时间积分的 精度。因此，在接触碰撞分析中广泛采用显式算法13,14。 hertz首先研究了静态接触问题15。他假定接触面是平面，而且跟接触面相接触的 结构已经有一个预先设定的初始应变。hughes等在接触碰撞有限元领域做出了有益的 探索，他们通过接触面节点对之间的几何关系判断接触穿越情况，这是一个很粗糙的模 型，由于接触面之间的接触点往往不在构成接触面的节点上，仅仅通过节点对的几何位 置不能准确的判断接触情况，目前己经没有人使用15-17。lawrence livermore国家试验 室的hallquist在开发著名的有限元动态分析程序dyna时，采用节点面元 (node-segment)模型来处理接触问题18-20。hallquist的理论被很多动态接触/碰撞分析有 限元理论采用，如nike,hemp,toody等21-23。hallquist的接触算法是后来的主/从滑 移面接触算法的基础。主/从滑移面接触算法是广泛使用的接触算法。主/从滑移面算法 预先将潜在的接触面对分为主面和从面， 然后将主面上的面元与从面上的节点之间的几 何关系作为穿越判据。 主/从滑移面接触搜索算法可以分为接触搜索算法和接触力算法两 部分，其中接触搜索算法一般来说可以分为两步：全局搜索算法和局部搜索算法。全局 搜索就是要粗略地确定可能发生接触的节点和单元，即找出潜在接触块（在此，我们称 相互接触的一个节点与一个单元构成接触块24）；局部搜索就是要针对已找到的潜在接 触块，精确地定位节点与单元之间的接触状态、接触位置及贯入量。目前有影响的全局 接触搜索算法主要包括主从面算法、单曲面算法、级域算法、位码算法及nbs算法； 局部搜索算法主要有基于“点面”的算法、基于“小球”的算法、基于“光滑曲面（曲 线）”的算法。 接触力算法目前主要有拉氏乘子法和罚函数法。在拉氏乘子法中，要求精确满足接 触界面无穿透的约束条件。 因此， 拉氏乘子法是精确方法， 但由于它引入了新的未知量， 增加了方程组的未知量数目，在系数矩阵的对角线上出现零元素，使方程组不再是解耦 的，从而使问题变的更难求解25 。由于此方法与显式算法不相容，因此一般不使用该 方法计算接触力。 罚函数法不增加系统未知量的数目， 通过引入的罚参数与界面穿透量的乘积作为接 触力，使无穿透的约束条件近似得到满足。该方法简单易用，但当所引入的罚参数特别 大时，将使方程组呈病态。同时，罚参数的选取也与单元的刚度有关，但没有一个明确 的规则，要求使用者具有一定的经验25。 为了克服罚函数法和拉氏乘子法各自的缺点， 出现了改变拉氏乘子法非零项的摄动 南京航空航天大学硕士学位论文 3 拉氏法(perturbed lagrangian method) 26、拉氏乘子法和罚函数法相结合的增广拉氏法 (augmented lagrangian method)27。这两种方法主要用于静力分析或使用隐式算法的动 力分析程序中。为了在显式有限元程序中有效地进行接触力计算，zhong提出了防御节 点法(defence node algorithm)26, carpenter等提出了前增量拉氏乘子法(forward increment lagrange multiplier)28, pollock等提出了前增量位移中心差分方法(forward increment displacement-central difference)25。这几种方法中，防御节点法不需要迭代， 效率较高；后两种解决了拉氏乘子法与显式算法的相容问题，但计算时间较长。 最早的接触碰撞有限元分析系统有美国lawrence livemore 国家实验室开发的 hemp22和sandia国家实验室开发的toody23。 随着接触碰撞有限元理论的成熟， 陆续 出现了一大批能够进行接触/碰撞分析的有限元系统， 其中比较著名的通用有限元分析商 业软件有msc.dytran, dyna, marc等。类似的程序基本思路是： 把实体结构(网格划分为二维单元或三维单元)相互接触的二个表面称为主表面(其 上的节点称为主节点)和从表面(其上的节点称为从节点)，在每一计算的时间步长内，检 查从(主)节点是否穿越主(从)表面，如果没有穿越则对该节点不作处理；如有穿越，则 对该节点与相应的表面进行处理，其方式多种多样：(1) 引入一个较大的交界面接触力， 其大小与穿越深度成正比，称为罚函数，其物理意义相当于在从节点与被穿越的主表面 之间放置一法向弹簧，以限制穿越；(2) 减小时间步长，重复上次的计算，使该节点不 穿越对应的表面，而且对其中刚到达主表面的从节点施加接触条件，对所有已经与主表 面接触的节点施加约束条件以保持其接触。此外，检查节点与对应接触表面之间是否存 在受拉交界面，如有拉力存在，则用释放条件使该节点脱离对应接触表面。关于施加的 接触条件或接触力有很多种定义方法，一般来说，其中或多或少引入了一些经验系数以 保证其算法稳定和结果精度。 在这些程序中，存在着几个重要问题：(1) 在一些大规模的接触碰撞问题中，计算 时间过长29，其中接触搜索占总计算时间的比例最高可达608030。因此，接触搜 索算法的计算效率至关重要；(2) 计算精度不理想。只能定性地得到变形趋势，定量的 变形结果与实验结果还有较大出入31。 1.3 本文的主要工作 本文针对接触碰撞有限元分析现状，主要在以下几个方面进行了研究。 (1) 建立了超参壳元的非线性动力分析的tl计算格式，并对其求解过程进行了深入 的分析与研究。 (2) 分析了现有接触搜索算法，然后对基于点面算法的局部搜索算法做了一定的改 进，改进后的局部搜索算法在进行精确的接触判断前，再次对接触测试对进行筛选，有 接触碰撞问题的算法研究 4 效地提高了局部搜索的精度与效率。 (3) 推导了直接以接触点作为防御节点的接触力计算公式。 南京航空航天大学硕士学位论文 5 第二章 接触碰撞问题的有限元理论及算法 本文研究的接触碰撞问题指的是基于连续介质力学的两个或者多个连续体之间具 有共同边界， 因相互作用而产生的结构变形或者破坏的现象。 计算响应时间为0秒到0.1 秒之间。跟一般的非线性动态有限元不同，接触碰撞有限元法涉及到更多的研究领域。 2.1 接触系统的约束条件 考虑两个弹性体的接触问题，如图2.1所示。假定有两个物体a和b在 1 t 时刻占据 空间 a v和 b v（图中虚线） ，当时间变为 2 t时，系统的构形发生了改变，两个物体发生了 接触（图中实线） 。我们可以用下面的关系来描述该系统32，其中 a s和 b s分别指a物 体和b物体的表面积， a f和 b f分别指a物体和b物体边界上的接触力。 图 2.1 接触前后系统的构形 a. 无贯入条件 () abb vvs= (2.1) b. 两物体间存在共同边界 ab ss (2.2) c. 在边界上存在接触力 0 ab ff+= (2.3) d. 法向接触力不能处于拉伸状态，即 0fn ,a b= (2.4) e. 切向接触力必须符合库仑摩擦定律 tn ff 是 否 说明高斯点此前已经屈服。现检查是否有 ttt ee ss 说明高斯点此前并未屈服，现检查是否有 0t ey s 是 否 是 否 说明高斯点此前已经屈 服且应力仍在增加，因此 所有的应力增量 () t e rs，必须退回到 屈服面上， 如图 2.3 所示。 令缩减系数1r= 说明高斯点正在弹 性卸载，直接转到 步骤 e 说明高斯点已经屈服，大 于材料屈服值的应力部分 () t e rs， 必须退回到屈 服面上，如图 2.4 所示。 令缩减系数 0t ey ttt ee sbc r acss = 说明高斯点仍为弹 性， 直接转到步骤 e 注： tt e s 是tt时刻按弹性性能求得应力的等效应力； 0 y 材料工作强化开始前的初始单向屈服 应力。 图 2.3 弹塑性连续体中已屈服点上的应力增量变化 接触碰撞问题的算法研究 16 图 2.4 弹塑性连续体中初始屈服点上的应力增量变化 对于已经屈服的高斯点，应力分为两部分，即落在屈服面以内的部分和落在屈服面 以外的部分，其计算公式分别为() ()1 t tt e r +ss和() t e rs。落在屈服面以外的部 分必须通过一定手段拉回到屈服面上来，本文采用切向预测径向返回的方法35,39。 ? 切向预测 高斯点进入塑性状态后，本构关系为 ep =sd (2.67) ep d为弹塑性矩阵，可用下式计算 t dd ep t d h = + d d dd d a (2.68) 其中 d =dda (2.69) 这里，d为弹性矩阵，a为流动矢量，按下式计算 2 3 222 2 t xyzxyyzxz j =a (2.70) 其中， , xyz 是应力偏量，表达式如下 南京航空航天大学硕士学位论文 17 () () () 1 3 1 3 1 3 xxxyz yyxyz zzxyz =+ =+ =+ (2.71) 将式(2.68)代入式(2.67)中，得到 tt dddd tt dd hh = = + d dd d sdd d ad a t d t d h = + = a d d d a da d d d (2.72) 在计算s时，式(2.72)中的d 用() t e rs代换。 这样可得到t时刻应力的预测值 ts () ()1 t ttt e r =+ssss () t tt ed =+ssd d (2.73) ? 径向返回 因为在切向预测的计算过程中，使用的 ep d是起点切线刚度，又由于加载曲线是外 凸的，所以s是在加载曲面的切线方向，并且 ts 总是在加载曲线之外，但是屈服准则 要求应力 ts 只能在加载曲面之上或者之内， 所以接着采用径向返回的方法以求得满足屈 服条件的 ts ，具体做法是令 tt r=ss (2.74) 这里，r是比例因子，由下式计算 0 t yp t h r s + = (2.75) 其中， ts 是应力 ts 的等效应力，由式(2.65)计算， t p 是t时刻的等效塑性应变，h是 线性应变强化参数，分别由下面两式计算 () t ttt pp r =+ (2.76) 1 t t e h ee = (2.77) 这里，() t 是() t 的等效应变， t e是弹塑性特性曲线的斜率，e为弹性模量。 当然，若超出屈服面的应力部分() t e rs比较大，且应力点在大曲率区附近，则通 过上述方法还是不能精确预测屈服面上的最终点d， 图2.5描绘了这种情况下弹性应力 退回到屈服面的过程，即将() t e rs分成若干段，分段将超出屈服面的部分拉回到屈服 面上，显然，划分的段数越多，精确度也就越高，但同时计算费用也就越高，为均衡这 接触碰撞问题的算法研究 18 两方面，文献39给出了一种折中的方案，即将超过的应力部分() t e rs划分成m部分， m为小于 () 0 81 t e y rs + (2.78) 的、最接近的整数值，其中( ) t e s是() t e s的等效应力， 0 y 是材料工作强化开始前的 初始单向屈服应力。 e. 对于弹性高斯点直接得出本步的应力 对于处于弹性加载或弹性卸载的高斯点，直接将 ttt e =+sss作为本步应力。 2.6 时间历程响应的直接积分法 有限元动力平衡方程(2.24)的求解方法可以分为两类：直接积分法和模态法。模态 法在时间积分之前进行了基的变换， 即从有限元坐标基变换为广义特征空间 2 =km 的特征向量基。从数学上看，由n个特征向量所张成的空间和由有限元n个节点位移所 张成的空间是相同的，所以两种分析所得到的解必定是相同的。模态法最终用duhamel 积分求解通过基变换后形成的n个互不耦合的方程，而直接积分法则直接对方程进行逐 步数值积分。模态法应用于低频长时间的动力响应问题，如地震、结构振动等；而直接 积分法适用于高频、短时间的动力问题，如冲击问题。直接积分法在解决非线性系统、 时变系统具有很强的优势，本文详细讨论了这种方法。 通常的直接积分法基于两个概念： 一是只是在相隔t的离散点上而不是在任一时刻 上满足运动方程，即方程(2.24)中的惯性力、内力、外载和接触力的平衡是在求解区间 上的一些离散时刻上获得； 二是可以假定位移、 速度和加速度在每一时间区间t内的变 化形式。 当然， 这些变化形式决定了动力平衡方程解的精度、 稳定性和求解过程的费用。 基于求解tt+时刻的位移、速度和加速度所用的是时刻tt+还是时刻t上的平衡 条件，直接积分法分为隐式积分法和显式积分法两大类。隐式积分法在每一个时间步都 需要迭代求解，对于线性问题，它允许较大的时间步长，是无条件稳定的，而对于非线 性问题，为保证迭代过程的收敛性与精度，时间步长须取的足够小，因而工作量很大， 一般在碰撞问题中不采用该方法。 显式积分法结构简单，易于处理非线性问题，同时它也避免了求解大规模代数方程 组，程序实现时算法简单，占用内存小。显式积分法是条件稳定的积分算法，时间步长 不能超过一定的限度。由于冲击问题一般是非线性的，要求解的动态响应整体时间也不 会太长，所以接触碰撞问题大多采用显式时间积分方法。 中心差分法是最典型的显式求解方法，广泛应用于冲击动力学问题中，下面对该方 法进行详细描述。 南京航空航天大学硕士学位论文 19 2.6.1 基于中心差分法的位移计算公式 在中心差分法中，加速度和速度可以用位移来表示35： () () 2 1 2 1 2 tttttt ttttt uuuu t uuu t + + =+ =+ ? ? (2.79) 其中t为均匀的时间步长；, ttttt uuu + 分别为t时刻及其前、后t时刻的节点位移。 将上式代入式(2.24)中，得到tt+时刻的节点位移 ttu+ () 12 2 tttttttt c t + =+ umqfrm umu (2.80) 当m为对角阵时，上式变为解耦的5n（n为节点总数）个方程 () 2 1 2 tttttttt iiciiiiii i utqfrm umu m + =+ ()1,2,5in=?(2.81) 在一般情况下，运动的初始条件都是给出0t =时刻的初始位移 0u 和初始速度 0u ? ， 而不会给出更前一个t时刻的位移 tu ， 所以无法直接按上式进行第一步的后推， 为此， 又利用中心差分式 () () 00 2 0 1 2 1 2 tt tt uuuu t uuu t =+ =+ ? ? (2.82) 把上式代入0t =时的动力方程 0000 c +=+?m urqf (2.83) 得 () 1200000 1 22 2 t c tt =+ ?umqfrm umu (2.84) 这个方程右端项均为已知，可以解得 t u。 如上，第一步先由 00 ,?uu求得 t u，之后就可以由 0 , t uu按式(2.80)推出 t u，依次 类推，可解得各时刻的节点位移。 中心差分法是条件稳定算法，即利用它求解具体问题时，时间步长必须小于由该问 题求解性质所决定的某个临界时间步长 cr t，否则算法将不稳定40。中心差分法的稳定 性条件是 max 2 n cr t tt = (2.85) 其中 max 是有限元系统的最大频率， n t是最小振动周期。理论上可以利用一般矩阵特征 值问题的求解方法得到 n t，实际上只需要求解单个单元的最小振动周期即可，又因为系 统的最小振动周期 n t总大于或等于单个单元的最小振动周期 e n t，因此，通过研究单个 单元，就可简单的估计出临界时间步长，而求解偏安全。在工程应用中，通常采用下面 接触碰撞问题的算法研究 20 的近似表达式来估计时间步长41 cr l t c (2.86) 其中l是单元的特征长度，是任意两个节点的最小长度；是取决于单元类型的系数；c 是声速，可通过下式计算得到 () 2 1 e c v = (2.87) 其中, ,ev分别是单元的弹性模量，密度和泊松比。 通过式(2.86)和式(2.87)我们可以看出，临界时间步长与单元的尺寸和材料系数有 关， 当各单元材料系数相同时， 中心差分法的时间步长由网格中的最小尺寸的单元决定， 它的尺寸越小，将使时间步长越小，从而使计算费用越高。 2.6.2 质量矩阵的计算 按照式(2.28)计算得到的单元质量矩阵称为一致质量矩阵，在中心差分法中，一般 不用一致质量矩阵，而使用集中质量矩阵。集中质量矩阵的构造方法很多，针对超参壳 元，cheung42建议首先生成一致质量矩阵，然后，如果对角项为位移项，则将本行的位 移项相加得到新对角项，本行（列）的其它项置0；如果对角项为转角项，则将本行的 转动项相加得到新对角项，本行（列）的其它项置0。为了计算方便，本文采用的方法 是43： a. 求出单元的总质量 e m； b. 求出单元一致质量矩阵中的对角元素 ii m； c. 求出对角元素和 0ii sm=； d. 调整对角项 0 e iiii mm ms=； e. 其它非对角元素清零。 这样可以生成对角质量矩阵且保证单元总质量不变，具体形式如下： ( ) 110 220 0 e e e e kk m ms m ms m ms = m ? ()5kn= (2.88) 2.7 计算流程图 图2.5给出了接触碰撞问题的计算流程图 南京航空航天大学硕士学位论文 21 开始 输入给定的几何尺寸、边界条件、材料性能数据；给出；选择时间步长 00 , u u ? t 建立壳的节点坐标系；局部坐标系；计算应变矩阵 构造质量矩阵；计算外载对应的等效节点外力 ； 根据初始应力计算节点内力 计算 t u 0t = 依据上两步位移；计算该时间步位移t t+ u 节点位移输出 m q b ? l tt= 计算该时间步单元内力；外载对应的等效 节点外力；接触搜索及接触力计算 ttt ttt + = = uu uu ttt= + , ttt u u 结束 r y n 图 2.5 接触碰撞系统的计算流程图 接触碰撞问题的算法研究 22 2.8 本章小结 首先描述了接触系统的约束条件，然后针对超参壳元进行了有限元的离散模型分 析，利用虚功原理建立了tl格式的接触系统运动平衡方程。讨论了载荷、初始与边界 条件的处理方法，推导了超参壳元在大变形条件下的green应变的表达式，给出了同时 考虑材料非线性与几何非线性时的应力计算过程， 最后给出了接触碰撞系统的计算流程 图。 南京航空航天大学硕士学位论文 23 第三章 接触搜索算法 在第一章中，对目前的接触搜索算法做了一个简要的综述，在本章中，将对现有的 接触搜索算法进行详细的分析和研究，给出这些方法的优点和缺点，然后在此基础上， 给出一种新的、有效的接触搜索算法。 3.1 接触搜索算法概述 所谓接触搜索， 就是要确定在整个系统中有哪些部位发生了接触或者有哪些原已接 触的部位发生了滑移或脱离。在有限元法中，接触一般是针对节点与单元表面间的相互 关系而言的。因此，接触搜索就是要确定系统中各节点与系统中所有单元间的关系，即 是否有接触、贯入或脱开。对于一大型结构，单元数量及时间积分的步数均很大，因此 接触搜索算法的效率和精度对结构分析至关重要。 在确定确切的接触状态前， 我们将接触对和潜在接触对统称为“接触测试对”（在此， 我们称相互接触的一个节点与一个单元构成接触对26， 具有接触可能的一个节点与一个 单元，或者距离很近的一个节点与一个单元称为潜在接触对）。 接触搜索的任务首先是找出系统中存在哪些接触测试对， 然后针对找出的每个接触 测试对确定其接触状态，即接触、贯入或脱开。对于接触或贯入的情况，接触测试对已 经转变为实际的接触对，此时还要找出接触点及贯入量。 上述寻找接触测试对的过程，一般被称为全局搜索；确定接触测试对的接触状态及 找出接触点与贯入量的过程，一般被称为局部搜索30。在早期的接触搜索算法中，一般 将全局搜索与局部搜索结合在一起进行的，不进行区分，如广泛使用的主从面算法19 就是如此。但是，近年来，随着问题规模的不断扩大，对二者进行区分的必要性越来越 明显。 利用全局搜索，就可以在做精确的局部搜索前，对单元面进行筛选，从中排除在一 次或几次求解循环（求解循环在隐式算法中是指一次迭代运算；在显式算法中是指一个 增量步的运算）过程中不可能发生接触的单元，只留下可能发生接触的测试对。这样， 在局部搜索过程中，就可集中力量对保留下来的接触测试对进行更精确的计算，从而精 确确定节点与单元面的相互关系。通过全局搜索与局部搜索的结合，就可达到既减少搜 索时间，提高搜索效率，又提高搜索精度的目的。 根据上述特点，全局搜索过程有时也被称为预搜索过程，局部搜索过程有时也被称 为精确搜索过程44。 在建立新的搜索算法时，全局搜索算法的关键是运算量，即计算效率；局部搜索算 接触碰撞问题的算法研究 24 法的关键是计算的精度与稳定性。 3.2 已有方法的研究 正如前文所述，接触搜索算法分为全局搜索算法与局部搜索算法，下面分别对这两 种算法的目前研究情况作一详细的分析。 3.2.1 全局搜索算法 a. 强力搜索方法 这是最容易被理解，也是在小规模的接触问题中常采用的一种方法。实现过程非常 简单，即对系统中所有单元进行循环，对于某一特定单元，对系统中所有节点循环，依 次判断每个节点到单元面间的距离，根据这个距离找出所有接触测试对，然后调用局部 搜索算法对每个接触测试对进行分析，计算出实际的接触状态，这种搜索方法称为“强 力搜索算法”(brute force search method)24,30。若假定系统的单元数与节点数均为n，则 上述搜索过程包含两个对所有单元和所有节点的算法，运算量即为o(n2)。对于单元数 达数十万的大型结构，接触搜索的运算量随单元数目的增加呈二次方关系增加，所以这 种算法不实用。 b. 主从面算法 主从面算法(master-slave algorithm)19是最早用来求解大规模冲击接触问题的接触 搜索算法，也是目前商用显式有限元软件使用最为广泛的接触搜索算法之一，在 dyna3d和abaqus均有该算法。该算法最初是由hallquist19提出的，以后又有很多 学者对该算法进行了发展。在该算法中，相互接触的两个边界面被分别定为主面(master surface)和从面(slave surface)，主面上的节点被定义为主节点(master node)，从面上的 节点被定义为从节点(slave node)。在该算法中，从面上的节点不允许穿透主面，主面 上的节点可以穿透从面，因此，应用时只需搜索与主面接触的从节点，即找出由从节点 与主面上的单元面所形成的接触对。 该算法包含三个步骤：首先是针对每一个从节点，找出与其距离最近的主节点；其 次在包含此主节点的所有单元中找出离从节点最近的主单元面，形成接触测试对；最后 在此测试对中精确定位节点与单元面的关系。 尽管主从面算法被许多有限元软件采纳，但这种算法存在明显的缺点：首先，它不 能用于相对初始构形有严重变形的接触表面， 因为在该算法中认为接触只发生在从节点 与包含离从节点最近的主节点的单元面之间； 其次， 该算法只能处理两个表面间的接触， 且在有限元前处理程序中就要明确说明主面与从面，因此，它不能处理如单一曲面在发 生屈曲变形时所产生的自身接触。 南京航空航天大学硕士学位论文 25 c. 单曲面算法 单曲面算法(single surface algorithm)30主要是为了解决主从面算法不能处理单一 曲面自身接触的情况出现的。它最早由benson和hallquist30提出，这是目前dyna3d 进行接触搜索的主要算法。 在这种算法中，不需要提前人为地指定主面与从面，而将所有结构（不管是否属于 同一部件）的表面看成是一个表面。搜索时对结构中的所有节点与单元面进行循环。为 了更迅速地找出系统中的所有接触测试对，该算法使用了子域排序算法(bucked sort algorithm)来寻找接触测试对。像多数排序算法一样，排序是在一维方向上进行的，通 过三层嵌套排序来完成三维空间的搜索过程。 单曲面算法尽管可以处理严重变形的表面间接触， 但由于它的搜索是通过三个嵌套 的排序算法实现的，因此，该算法编程复杂，运算量大，同时，可并行性差。 d. 级域算法 级域算法（hierarchy-territory algorithm）是由zhong24提出的。级域法是针对多 物体接触提出的，它使用了两个基本概念：级（hierarchy）和域（territory）。在该算 法中， 将整个结构分成不同级别的若干级别， 即接触体 （contact body） 、 接触面 （contact surface） 、 接触片 （contact segment） 、 接触边 （contact edge） 、 接触点 （contact node） 。 在搜索时，先在较高级别的两个域之间进行，当发现两个处于同一级别的域不存在公共 部分时，就不必进行下一个级别间的接触检测，如两个域间存在公共部分，则在下一级 的两个域间进行搜索检测。 在这种算法中，实际构成的接触测试对只占系统中节点与单元面全部组合的一部 分，避免了许多不必要的操作，计算效率高。但是，这种方法编程和使用并不方便，因 为它涉及不同级别接触元素的定义和管理，在接触搜索过程中也要求逐级进行，使得数 据的管理很复杂，同时它的可并行性也很低。 e. 位码算法 位码算法（position-code algorithm）是oldedburg和nilsson32于1994年提出的， 它的基本思想是将三维空间划分成若干立方格，根据立方格的位置，为每个立方格赋予 一个编码， 在搜索时， 对每一个单元， 根据其位置可计算出所有与此单元相交的立方格， 然后通过折半搜索，找出这些立方格中的节点，这样，这些立方格中的节点均与该单元 构成接触测试对。 位码算法避免了三维嵌套搜索，减小了编程及接触搜索的计算量，对大规模计算很 有好处。 f. nbs算法 nbs算法（no binary search）是1998年由munjiza和andrews45针对离散粒子的 接触碰撞问题的算法研究 26 接触系统提出的搜索算法。 该算法适合于离散粒子的大规模动态接触系统， 但这种算法只能用于尺寸相同的圆 球的接触系统，而有限元的接触系统中单元形状各异，该方法不能直接纳入有限元计算 中。王福军46将位码算法的基本思想融入nbs算法中，另外提出了新的搜索算法，新算 法计算量小，可适用于任何形状的单元，包括实体单元、壳单元、板单元等，是一种比 较有前景的全局搜索方法。 3.2.2 局部搜索算法 a. 基于点面算法的局部搜索算法 基于点面算法的方案是最早用于局部接触搜索的算法， 也是到目前为止应用最为广 泛的局部搜索算法。 在这类算法中，一般将给定的有限元节点称为“从节点”（slave node），与其相对的单元面称为“主面”（master segment）， 在这类方法中，接触点的 寻找变成了求解点到曲面的最小距离问题，即求解偏微分方程组 (,) (,)0 (,) (,)0 cchcc cchcc u wu w u u wu w w = = x xx x xx (3.1) 所得到的(,) cc u wx，实际是节点 h x在单元面( , )x u w上的投影点，即节点 h x到单元面 ( , )x u w的最近点。当节点到单元面的距离为0或负值时，表示节点已与单元面接触，此 时称这个投影点为接触点。得到(,) cc u wx后，可通过计算接触点处的法向矢量 , , (,)(,) (,)(,) uccwcc uccwcc xu wxu w xu wxu w = n (3.2) 得到主击节点到单元面的距离 n g (,) nhcc gu w=nxx (3.3) 直接求解偏微分方程组(3.1)是件很困难的事47，因为，第一，建立单元面的曲面方 程非常耗时；第二，采用newton-raphson迭代算法求解时，算法不稳定，经常不收敛， 初值很难选。为此许多学者47-49采用了一些特殊处理方法来解决这个问题，从某种程度 上解决了迭代算法的不稳定等问题 b. 基于小球算法的局部搜索算法 a) 小球算法 小球算法（pinball algorithm）是belytschko50于1991年提出的。该算法将接触单元 近似看成是由等体积的小圆球构成，单元中心为球心。若两个单位中心距离小于两个球 半径之和，即判为接触。 b) 圆管算法 圆管算法是whirley51于1994年在小球算法的基础上提出的。用“球”来构造单元的 南京航空航天大学硕士学位论文 27 角，用“管”来构造单元的边。这样，两个单元在边界处的接触就比较好判断了，计算精 度比小球算法要高很多，该算法已经在dyna3d中成功应用。 c) 球形排序算法 球形排序算法（