清华大学计算固体力学第十次课件接触-碰撞

1、非线性有限元第10-碰撞 计算固体力学1第10章 接触-碰撞引言接触界面方程 摩擦模型 弱形式 有限元离散 显式方法 21 引言 在制造加工过程的仿真中包括接触和碰撞。例如，在薄金属板的成型中，模具和工件之间接触面的模拟；汽车碰撞和挤压成型的模拟，都需要有接触滑移界面。接触问题的失效形式是相互作用的接触面发生破坏。 对于接触-碰撞，将展示Lagrangian网格的控制方程和有限元程序。在接触中，物体的控制方程与前面介绍的方程是一致的，但在接触界面上，需要增加动力学和运动学的条件。 关键条件是不可侵彻性条件：即两个物体不能互相侵入的条件，不可侵彻性不能表示为一个简单的方程，所以，发展了几种简单的

2、方法。其中两个基于最近点映射的形式： 1）显式动态问题适用的率形式； 2）主要适用于隐式方法的平衡解答。 此外，经典的Coulomb摩擦模型和界面本构模型。31 引言 给出控制方程的弱形式，处理接触界面约束的4种方法：1 Lagrange乘子法，2 罚方法，3 增广的Lagrangian法，4 摄动的Lagrangian法。 由Lagrange乘子法，在接触问题的离散化中，在接触界面上乘子必须是近似的。乘子必须满足法向面力是压力的约束。 在罚方法中，面力不等式源于Heaviside分步函数，该函数被嵌入在罚力之中。4 1 引言 接触-碰撞问题是属于最困难的非线性问题之一，因为在接触-碰撞问题中

3、的响应是不平滑的。 当发生瞬时接触时，垂直于接触界面的速度是瞬时不连续的。对于Coulomb摩擦模型，当出现粘性滑移行为时，沿着界面的切向速度是不连续的。 接触-碰撞问题的这些特性给离散方程的时间积分带来了明显的困难，削弱了Newton算法的功能。 因此，选择适当的方法和算法是至关重要的，并且在获得强健的求解程序中，规则化的技术是非常有用的。52 接触界面方程模拟接触-碰撞问题的标记 接触界面包含两个物体表面的交界。 接触界面包括两个物体处于接触的两个物理表面，它们是重合的，在数值计算中，两个表面一般不重合，分为主控和从属表面。接触界面是时间的函数，确定它是接触-碰撞问题解答的重要部分。 在通

4、用有限元软件中，接触算法能够处理多个物体的相互作用，然而多个物体的接触包含成对物体的相互作用。因此，从考虑两个物体的问题入手。62 接触界面方程 在主控接触表面的每一点建立局部坐标系统，可以构造相切于主控物体表面的单位矢量：物体A的法线给出为 在接触界面上有即两个物体的法线方向相反。以局部分量的形式表示速度场 在三维问题中希腊字母下角标的取值范围为2；当问题是二维时，接触表面成为一条线，取值为1。 72 接触界面方程 物体由标准场方程控制：质量、动量和能量的守恒，应变度量，以及本构方程。接触增加了条件： 1 在界面上，两个物体不可相互侵入和面力满足动量守恒； 2 横跨接触界面的法向面力不能为拉

5、力。 按照要求分类： 1 对于位移和速度的要求作为运动学条件， 2 对于面力的要求作为动力学条件。 8不可侵彻性条件2 接触界面方程一对物体的不可侵彻性条件可以表示为交集为零 例如，如果物体在旋转中，对于P 点接触Q 点是可能的，而一个不同的相对运动可能导致P 点与S 点接触。结论是，除了以一般的形式，找不到其它的方程表示P 点没有侵入物体A 的事实。 两个物体不允许重叠，这可以视为一个协调条件。对于大位移问题，不可侵彻性条件是高度非线性的，并且一般不能以位移的形式表示为一个代数方程或者微分方程。其困难源于在一个任意运动中，不可能预先估计到两个物体的哪些点将发生接触。9不可侵彻性条件运动学2

6、接触界面方程 由于以位移的形式表示交集为零的公式是不可能的，所以，在接触过程的每一阶段中以率形式或者增量形式表示不可侵彻性方程是很方便的。其率形式应用到物体A和B上发生接触的部分，即是位于接触表面上的那些点 两个物体的相互侵彻速率 利用 点乘 得到上两式 并且利用法线是正交于与平面相切的单位矢量的事实。 10不可侵彻性条件2 接触界面方程 对于在接触表面上的任意点，不可侵彻性条件限制了相互侵彻速率成为负值，表示当两个物体发生接触时，它们或者必须保持接触 或者脱离 。 对于接触区域上的所有点当满足上式时，精确满足不可侵彻性条件。然而，公式与交集为零不是等价的。在大多数数值方法中，仅在瞬时时刻注意

7、到该式，对于接近分离而没有接触的点，相互侵彻是可能的。该式仅适用于处于接触的成对点。11不可侵彻性条件2 接触界面方程 上面强化公式将不连续性引入速度时间历史中。在接触之前，法向速度是不相等的，而在随后发生碰撞，法向速度分量必须满足公式。在时间上的这些不连续性使得离散方程的时间积分变得很复杂。 间隙率是相互侵彻率的负数。 相对切向速度给出为12动力学面力条件2 接触界面方程 横跨接触界面，面力必须服从动量平衡。由于界面上没有质量，这就要求两个物体上面力的合力为零： 在法线方向上，不考虑在接触表面之间的任何粘性，法向面力不能是拉力。其条件表示为 这个条件要求 为正数，物体B上的面力在A的单位法线

8、上的投影， 它指向物体B。对应于物体A和B，注意到上面的表达式是不对称的。为了定义法向面力，选择其中一个物体的法向，并且物体法向面力的符号将取决于选择的这个法向。 13定义切向面力为面力条件2 接触界面方程切向面力投影到主控接触表面上的面力合力，由动量平衡要求当应用接触的无摩擦模型时，切向面力为零 在前面建立接触界面方程中，尽管选择了其中一个物体为主控物体，当两个接触表面是重合时，且满足公式 ，对于物体，这些方程是对称的。因此，选择哪个物体作为主控物体是没有关系的。但是，当两个表面不重合时，如在大多数数值求解中，则主控物体的选择会改变结果。14单一接触条件2 接触界面方程由法向面力不是拉力的条

9、件 由不可侵彻性条件得到单一接触条件 当物体发生接触并且保持接触时， 而当接触停止 ，并且法向面力消失 ，乘积总是为零。 这个方程也可以表示为接触力的法向分量不工作的事实。 15相互侵彻度量2 接触界面方程 在物体B上的点P侵入到物体A的内部，定义为至物体A的表面上任意点的最小距离。在用坐标 表示点P到物体A表面上的任意点之间的距离给出为 图示物体B上的点P已经侵入物体A。寻找相互侵彻的度量，它表示为 在物体A上的点Q是最接近于物体B上的点P：它是点P在A上的正交映射。16相互侵彻度量2 接触界面方程相互侵彻量 为上式的最小值，并且考虑到仅当P在物体A内部时才是非零的。通过检验法线到物体A在

10、上的投影，可以检验后面的条件：当投影是负值时，点P是在物体A的内部，因此有相互侵彻，否则P不在A的内部，没有相互侵彻。所以，相互侵彻的定义是 令相互侵彻度量对坐标求导数取值为零，求出 最小值。172 接触界面方程 考虑发生部分侵彻的两个表面。主控物体是9节点等参单元，所以表面A 的3个节点是二次映射定义： 从属物体B的表面为一条水平线，给出为例子10.1在例子中的相互侵彻已经被夸大了。注意到沿着界面有 对于在表面B上的点P，找到相互侵彻。求点Q 正交投影的最小值 182 接触界面方程例子10.1取最小化给出为数值求解上式的根为r0.2451，因此Q点位置为， 对于在表面B上的点P，找到相互侵彻

11、。求点Q 正交投影的最小值 192 接触界面方程 当两个物体是不光滑或者不是局部地凸状时，这种定义相互侵彻的方法将会遇到困难。如在图示情况下， 的最小值是不唯一的：这里有两个点为P的正交投影。在这种情况下，难以建立一种方法，能够唯一地定义相互侵彻的度量。通过一个有转折表面的侵彻，说明正交映射点求解的不唯一性。 相互侵彻度量202 接触界面方程 在ABAQUS/Standard和/Explicit中的接触模拟功能具有明显的差异。 在/Standard中的接触模拟或者是基于表面（surface）或者是基于接触单元（contact element）。因此，必须在模型的各个部件上创建可能发生接触的表面

12、。然后，必须判断哪一对表面可能发生彼此接触，称之为接触对。最后，必须定义控制各接触面之间相互作用的本构模型。这些接触面相互作用的定义包括如摩擦行为等。 在/Explicit中的接触模拟可以利用通用（“自动”）接触算法或者接触对算法。通常定义一个接触模拟只需简单地指定所采用的接触算法和将会发生接触作用的表面。在某些情况下，当默认的接触设置不满足需要时，可以指定接触模拟的其它形式；例如，考虑摩擦的相互作用力学模型。21 在电子工业中，为了评估产品的耐久性，仿真分析正代替跌落试验。电子产品跌落模拟 通用接触算法(general contact)22Global ModelLCD Submodel电子

13、产品跌落模拟 非线性有限元的应用232 接触界面方程 在/Standard和/Explicit中的力学接触算法具有本质的区别，体现在如何定义接触条件：1 /Standard在施加接触约束时采用严格主从权重：约束从属表面节点不能侵入主控表面。而主控表面上的节点原则上可以侵入从属表面。/Explicit包括这个公式，但是典型地默认应用平衡主从权重。242 /St和/Ex都提供了有限滑动接触公式。在/St中的二维有限滑动公式要求主控表面是光滑的；而在/Ex的主控表面是由面元构成的，除非是解析刚性表面，其主控表面可以是光滑的。3 /St和/Ex都提供了小滑移接触公式。在/St中的小滑移公式根据从属节点

14、的当前位置向主控节点传递载荷。/Ex总是通过固定点传递载荷。4 /Ex在接触逻辑中可以考虑壳和膜的当前厚度和中面偏移，/St不能。 /Ex通用接触算法的许多优势在/St中不具备。由于上述的差异，所以在/St分析中定义的接触不能导入/Ex中，反之亦然。2 接触界面方程253 摩擦模型 一般将切向面力的模型称之为摩擦模型。基本上有三种形式的摩擦模型：1 Coulomb摩擦模型，它是基于经典摩擦理论的模型。2 界面本构方程，由方程给出切向力，类似于材料本构方程。3 粗糙润滑模型，模拟界面的物理特性，常用于微观尺度。摩擦能量消耗，本构问题磨损材料损耗，本构问题263 摩擦模型 在Coulomb摩擦模型

15、应用于连续体中，它应用在接触界面的每一点，如果A和B是在x处接触，则Coulomb摩擦模型 是一个变量，由动量方程的解答确定。两个物体在一点处接触的条件意味着法向力 ，因此，两个表达式的右端项恒正。 条件(a)作为粘着条件：当在一点处的切向面力是小于临界值时，不允许相对切向运动，即两个物体为粘着的。 条件(b)对应于滑动条件，该方程的第二部分表示切向摩擦的方向必须是与相对切向速度的方向相反。273 摩擦模型Coulomb摩擦模型 Coulomb摩擦更类似于刚塑性材料本构。如果将切向速度理解为应变，将切向面力理解为应力，则上式中的a)式可以理解为屈服函数。当屈服准则不满足时，切向速度将为零；一旦

16、满足了屈服准则，则切向速度沿着由b)式确定的方向。这些特征平行于刚塑性材料模型。 283 摩擦模型界面本构方程 界面行为的本构模型源于凸凹不平导致的表面粗糙度，发生在微观尺度上即使是最光滑的表面。在滑动中，摩擦是由粗糙部分的相互作用生成的。最初的滑动引起了这些粗糙部分的弹性变形，所以，真正的粘着条件不会自然产生，是所观察到行为的理想化。在滑动过程中，伴随着粗糙部分的弹性变形是粗糙表面的研磨。粗糙的弹性变形是可逆的，而研磨是不可逆的，因此，自然地将初始滑动归属于弹性特性，而后面的滑动归属于塑性特性，属于非关联塑性。29摩擦滑移屈服表面 Mohr-Coulomb本构模型 滑移方向（塑性流动）是水平

17、的（沿Q 的方向）而不是垂直屈服面。这是非关联塑性流动的例子。对于连续体和多轴应力应变状态的行为，M-C准则具有普适性。它应用于模拟颗粒状的土壤和岩石。 M-C准则是基于这样的概念，即当任意面上的切应力和平均法向应力达到临界组合时在材料中发生屈服 c是内聚力，通过 定义内摩擦角。 3 摩擦模型303 摩擦模型界面本构方程对于Coulomb类型性能的屈服函数，对应于Coulomb摩擦条件： 在二维情况下，屈服函数包括以斜率为 的两条斜线。 在三维情况下，屈服函数是一个圆锥。 314 弱形式 对于Lagrangian网格，建立动量方程和接触界面条件的弱形式。当接触表面是作为L格式处理时，这一形式也

18、适用于ALE网格。为了简单，从无摩擦接触开始，后面介绍切向面力的处理。将下面的公式限制在面力或者位移边界的情况，进而描述相应的所有面力或者速度分量。接触表面既不是面力也不是位移边界。物体A的全部边界为 注意到 在接触问题中，奇异性发生在边界，如赫兹接触(1881)。然而与断裂力学不同的是，奇异性没有显示任何的工程意义，因为表面粗糙度的存在抵消了在应力中出现的更接近奇异性的行为。 32 自然变分原理是对原物理问题的微分方程和边界条件建立对应的泛函，再使泛函取驻值得到问题的解答，但是其未知场函数需要满足一定的附加条件。而广义变分原理（或称约束变分方程）不需要事先满足附加条件，采用lagrange乘

19、子法和罚函数法将附加条件引入泛函，重新构造一个修正泛函，将问题转化为求修正泛函的驻值。称为无附加条件的变分原理。 对于罚函数方法，将罚参数取正值，对修正泛函得到的近似解只是近似地满足附加条件，罚参数值越大，附加条件的满足程度就越好。而在实际计算中，罚函数只能取有限值，所以利用罚函数求解只能得到近似解。4 弱形式334 弱形式Lagrange乘子弱形式强加接触约束的通常方法是借助于Lagrange乘子。 令Lagrange乘子试函数为 相应的变分函数为 弱形式为其中侵彻率变分为 这个弱形式不等式是等价于动量方程、面力边界条件、内部连续条件（广义的动量平衡）和接触界面条件：不可侵彻性、法向面力的动

20、量平衡和无摩擦条件。Lagrange乘子场仅要求为C-1连续，因为它的导数并不出现在弱形式中。要求法向界面力是压力，这是对Lagrange乘子在试空间的一种限制。 虚功率为344 弱形式Lagrange乘子弱形式 比较Hu-Washizu变分原理，借助于Lagrange乘子，上面的方法是在弱形式中附加约束的标准方法。与Hu-W形式的唯一区别是约束为一个不等式。利用高斯积分关系，通过积分证明 从弱不等式可以推导出强形式，但是，必须考虑关于变分函数的不等式。一旦变分函数是不受约束时，对于与变分函数相乘的项的符号则没有限制，并且由密度原理该项必须为零。 从上式中的前两个积分得到 即在物体A和B上，满

21、足动量方程和自然边界条件。 Lagrange乘子强形式35从上式联立消去 得到关于法向面力的动量平衡条件： 4 弱形式Lagrange乘子强形式 在接触表面被积函数的所有项中，除最后一项外，变分函数是任意的，因此获得了等式， 由上式被积函数最后一项的变分函数为负值 推论 必然是非正的，即弱不等式表示不可侵彻条件 364 弱形式 在罚方法中，以沿接触表面施加不可侵彻性约束作为罚法向面力。对比于Lagrange乘子法，罚方法允许一些相互侵彻。然而，它更容易编程并且广泛应用。考虑两种形式的罚方法：1 罚数正比于相互侵彻率 的平方； 2 罚数为相互侵彻及其率的任意函数。 在非线性问题的应用中，第二种方

22、法是更有用的，因为速度相关罚数允许更多的相互侵彻。 在罚方法中，变分和试函数与在Lagrange乘子法中的项完全相同。对于罚方法，弱形式到强形式的等价性可以表述如下如果 其中 率相关的罚方法（penalty） 侵彻脱离 374 弱形式率相关的罚方法（penalty） 罚参数 可以是空间坐标的函数。相应于罚方法的弱形式不是一个不等式；由式中出现的Heaviside分步函数引入了接触问题的非连续性质。这种弱形式并不意味着不可侵彻条件，在罚方法中，它仅仅近似地得到满足。类似L方法，证明弱形式包含着强形式，得到 不像L乘子弱形式，罚弱形式没有强制横跨接触界面速度的连续性；事实上，横跨界面的速度将是不连

23、续的。可以从公式中得到不连续的量级，给出 在相关法向速度分量中的不连续反比于罚参数，随罚增加而减小。 384 弱形式接触约束方法 对于任何节点和单元之间的缝隙，g为正，不存在接触；g为负值或者为零，存在接触的节点对（单元对），施加约束条件。 通过简单地用乘子乘以上式给出的缝隙条件，给出Lagrange乘子的方法。接触缝隙条件 视为作用在每个节点上的一个“力”以防止侵彻。 39 在罚函数方法中，最终缝隙值将不是零，而是一个小量，取决于所选择的参数值。罚函数方法的优点是可以根据判别,能够给出一个可接受结果的实际需要的参数值，进行适当的调整。实际上，在一个复杂的问题中，这不是简单的工作，特别是对于涉

24、及梁、板、或壳单元与实体单元之间的接触问题。 4 弱形式罚函数方法 一种能够避免Lagrange乘子求解困难的方法是罚函数。罚方法是在接触节点之间附加一个刚度为beta的非线性弹簧。 404 弱形式摄动的Lagrangian弱形式(Perturbed)其中 上式被积函数的第二项为Lagrange乘子弱形式的摄动；由于 是大数， 很小。 增广的Lagrangian弱形式(Augmented) 是求解过程待定的正参数，被积函数第二项为侵彻率的高阶项。 通过在弱形式中附加上强化切向面力连续性的一项，所有上面的公式可以修改以便处理界面摩擦。 415 有限元离散 对于接触-碰撞问题的所有方法（罚方法、L

25、agrange乘子法等），在弱形式的表述中包括标准虚功率和接触界面变分的贡献合成。当在无接触状态时，可以精确地离散标准虚功率，因此，这里关注各种接触界面弱形式的离散化。 应用Lagrangian网格，包括UL和TL格式。在TL格式中，必须以变形表面的形式施加接触界面条件。离散化也适用于ALE格式，只要在接触表面的节点为L节点。它们不能直接应用到Eulerian格式，因为假设已经有了描述接触表面的参考坐标系。这个坐标系不能在一个E网格中定义。在L网格中，接触表面对应于网格边界的一个子集。 425 有限元离散 Lagrange乘子法横跨接触界面，两个物体的速度不一定必须是连续的，可以C0连续；相互

26、侵彻条件源于弱形式的离散化。注意到速度场的近似也定义了位移场的近似。试函数的空间定义如下 由于是L网格，以材料坐标表示速度场的有限元近似，也可写成单元坐标的形式，因为这两组坐标是等价的。速度场为 如果物体A和B的节点编号不同，两个速度场可写成为一个表达式：435 有限元离散 Lagrange乘子法在接触表面上Lagrange乘子场是C-1形状函数 变分函数为 弱形式为 r为残数，以节点速度的形式可以表示相互侵彻率： 可以消除在物体A和B的节点之间的区别，将上式表示为445 有限元离散 Lagrange乘子法利用形状函数，给出法向分量和速度的乘积的近似和侵彻率 将离散弱形式写成为 因为是不等式，

27、所以控制方程的提取必须小心。对于那些没在接触界面上的节点，可以直接地从第一个求和项中提取方程，由于节点速度的变分为任意的，得到标准的节点运动方程：在第一个求和项中的余下部分重新以接触界面的局部坐标系写成，并组合第二个求和项，给出 455 有限元离散 Lagrange乘子法 由于切向节点速度是没有约束的，对于节点速度的系数，弱不等式服从一个等式 上面的运动方程和接触界面条件，加之应变位移方程和本构方程，对于半离散模型组成了方程的完备系统。对于没有在接触界面上的节点运动方程是与没有约束的情况一样的。 以矩阵形式写出运动方程，通过组合这种形式与内部、外部和惯性的功率的矩阵形式，给出 考虑到变分项的任

28、意性，得到运动方程和相互侵彻条件 (为作用在每个节点上一个防止侵彻的“力”) 465 有限元离散 Lagrange乘子网格 构造Lagrange乘子网格具有一定的难度。一般说来，两个接触物体的节点是不重合的，如左图所示。因此，有必要建立一种方法处理不相邻的节点。 一种可能性表示在左图中，选择Lagrange乘子场中的节点为主控物体的接触节点。当一个物体比另一个的网格更加细划时，这种简单的方法是无效的。Lagrange乘子的粗网格则导致相互侵彻。 另一种方法是无论在物体A或B上出现一个节点，则放置Lagrange乘子节点，如右图所示。这种方法的不足之处是当物体A和B上的节点接近时，一些L乘子单元

29、非常小。这可能导致方程的病态条件。在三维情况下，这种方法是不可行的。475 有限元离散 Lagrange乘子网格 对于一般性的应用，Lagrange乘子必须单独构造网格，这种网格独立于其它任何网格，但是，至少细划到二者之中较为细划的那个网格程度。 Lagrange乘子法的主要缺点是对于L乘子网格的需要。像在简单二维例子中看到的，这可能引入了复杂性。在三维问题中，这个工作是更为复杂的。当接触界面变化时，网格必须随着时间变化。在罚方法中，没有必要建立附加的网格。 与罚方法比较，Lagrange乘子法的优点是没有用户设定的参数，并且当节点相邻时，接触约束几乎精确得到满足。当节点不相邻时，可能会稍微的

30、违背不可侵彻性，但是不会像罚方法那么明显。然而，对于高速碰撞，L乘子法常常导致非常不平顺的结果，因此，L乘子法更适合于静态和低速接触问题。 485 有限元离散 例10.2 一维接触-碰撞有限元方程接触条件 当两个节点发生接触时，节点1的速度必须小于或等于节点2的速度。如果两个节点速度相等，节点保持接触状态。运动方程为： 左端最后一项是在节点1和2之间由接触引起的节点力。这些力大小相等且方向相反，当Lagrange乘子为零时消失。除了在发生接触的节点，运动方程与无约束的有限元方程是一致的。 采用Lagrange乘子法495 有限元离散 对于取单位面积的对角化质量矩阵，方程可以写成为其中 ，第四个

31、方程为不可侵彻性条件 其中，kI 为单元 I 的刚度。在无接触的装配刚度矩阵中，即左上角的33矩阵是奇异的，而当应用附加的接触界面条件，完整的44矩阵成为非奇异的。通过组合G矩阵（E10.2.2）与装配刚度得到关于小位移弹性静力学的方程： 例10.2 一维接触-碰撞有限元方程505 有限元离散 例10.3 二维问题 采用4节点四边形单元模拟发生接触的两个物体。以物体A上边界坐标的形式（单元坐标投影到接触线上）写出在接触线上的速度场为： 单位法线为其中515 有限元离散 例10.3 二维问题由公式（10.5.17）给出矩阵：G矩阵（10.5.19）给出为 由一个线性场近似Lagrange乘子：

32、上面的矩阵类似于杆的一致质量矩阵：在节点1处的接触导致在节点2处的力，反之亦然。接触节点力是严格地沿着y方向；由于矩阵的奇数列为零，所有接触节点力的x方向分量为零，接触界面是无摩擦的。52 在碰撞接触中，罚方法发挥了类似的作用。在Lagrange乘子法中，在接触界面处发生碰撞的时刻，速度在时间上是不连续的，如图所示。这些不连续性在物体中传播，并且导致不可忽视的振荡。罚方法可以考虑作为接触界面条件的规则化：它平顺了不连续的速度，并且保持了动量守恒。通过允许两个物体的部分重叠，它仅放松了一个条件，即不可侵彻性条件。这是为了平滑结果所付出的一个小的代价。 5 有限元离散 在规则化的过程中，使得难以处理的解答，即由于解答引起的不连续性或奇异性，经过人为