一种基于表面参数化与残余柔顺性的柔性体接触动力学凸优化方法

《Applied Mathematical Modelling》：A convex optimization method for contact dynamics between flexible bodies with surface parameterization and residual compliance

　　•在锥互补框架中引入残差合规性校正。•采用边界首次展平方法对潜在接触面进行参数化处理。•几何结构、模态场以及残差合规性数据均离线存储。•接触区域可根据需求被识别并离散化，且对网格的依赖性较低。

引言
在工程系统中，涉及可变形体的接触现象极为普遍，应用范围涵盖多体机构、承重设备、精密装配以及机电装置等，例如间隙接头[[3],[4],[5]]、滚动与滑动部件[[6],[7],[8]]以及过盈配合操作[9,10]。与刚体接触不同，接触引发的局部高频结构响应会显著改变接触区域内的有效合规性，进而改变接触力的空间分布[11]。此外，接触过程中的间歇性啮合与分离，以及粘着状态与滑动状态之间的摩擦模式转换，还会带来明显的非光滑性[12]，从而导致主动接触位置随时间变化。因此，准确预测接触力是柔性体接触动力学建模与计算中的核心难题，而这又依赖于可靠的接触检测机制以及对接触引发局部变形的精准描述。

接触力计算的数值方法通常分为基于合规力的方法和基于互补性的方法两类[13]。基于合规力的方法通过对单边约束进行正则化处理，可轻松融入标准的动力学公式中[[14],[15],[16]]。然而，正则化参数需要在穿透误差与数值刚度之间取得平衡：参数值过小可能导致明显的相互穿透现象，而参数值过大则可能引发数值条件不佳、收敛性下降，同时还限制了允许的时间步长[17,18]。相比之下，基于互补性的非光滑公式[12,19]直接施加不可穿透约束，并能在不同接触模式转换时保持约束条件的一致性。二维摩擦接触问题往往可以表述为线性互补问题（LCP）[[20],[21],[22]]。但在三维情况下，库仑允许集变为锥形，原有的LCP公式不再适用[23]。通过多面体近似，虽能以牺牲摩擦锥近似为代价恢复LCP公式[[24],[25],[26]]，而精确锥形公式则能够保留真实的库仑摩擦锥，通常表现为非线性互补问题（NCP）[[27],[28],[29]]或锥互补问题（CCP）[[30],[31],[32]]。为提升计算效率，Tasora和Anitescu[30]提出了一种简化方法，将其转化为更易处理的二阶锥线性互补问题（SOCLCP）。后续研究进一步优化并扩展了这一公式[[33],[34],[35]]，同时开发出了高效的锥互补问题求解器[[36],[37],[38],[39]]以及非光滑时间积分方案[40,41]。

要实现可靠的接触力预测，还需要对结构灵活性进行准确且高效的表征。有限元模型能够精确描述接触引发的局部变形，已被广泛应用于梁-梁[42]、壳-壳[16]、实体-实体[43]以及刚体-柔性体[44,45]接触问题。然而，在潜在接触区域进行高密度离散化会导致系统规模过大。固定界面组件模态综合方法虽然能够减少内部自由度[[46],[47],[48]]，但对于范围较大或不断移动的接触区域，仍可能需要保留大量界面坐标，从而限制了降阶效果。为解决这一问题，人们提出了多种策略，包括模态集切换[49]、任意拉格朗日-欧拉公式[50,51]，以及在浮动参考框架公式体系下的各类降阶方法，如界面模态[52]、移动模态[53]、全局模态[54]以及基于载荷位置参数的模态[55]。自由界面组件模态综合方法[56]则提供了一种新途径，无需在整个潜在接触区域内预先定义界面坐标。只需少量自由界面法向模态即可表示整体动态响应，而残差合规性则可为因模态截断而产生的局部响应提供准静态修正，从而无需反复重构模态基函数即可更新界面上的载荷。不过，该方法中的接触处理是基于惩罚项实现的，而对于面积较大的潜在接触区域，直接存储残差合规性矩阵会占用大量内存。

除了接触力公式构建与局部变形建模之外，准确的几何表示对于接触区域的识别以及接触集的生成也至关重要。基于面片的方法能与柔性体建模框架自然结合[57,58]，但其分段平面表面以及接触点定位仍依赖于底层网格。而基于样条函数或计算机辅助设计（CAD）的表示方法，如等几何分析中所使用的那种[59,60]，以及隐式距离场方法[61,62]，虽然能实现更平滑的几何形状，但分别依赖于现有的平滑描述或距离场的实时更新。对于离散网格而言，网格到平滑表面的重建精度与稳定性取决于参数化的合理性，这样的参数化还能支持可调整的接触点集以及接触区域网格[64]。不过，直接从离散表面网格构建此类参数化仍然具有很大难度。元素级参数化与切平面映射通常仅能在局部邻域内重建平滑几何结构[[65],[66],[67]]。计算保角映射法则将表面块映射到平面域上。与最小二乘保角映射（LSCM）[69]以及基于角度的展平方法（ABF++）[70]相比，边界首次展平方法（BFF）[71,72]能够在保持较高保角精度的同时实现直接的边界控制。由此得到的规则矩形域便于进行多项式近似[73]，为接触区域的表示以及所需表面和界面数据的参数域拟合提供了合适的基础。

尽管在柔性体接触建模方面已取得显著进展，但仍存在三个相关挑战。首先，完整的有限元模型虽然能够描述局部接触变形，但会导致系统规模过大。模态降阶虽能减少动态自由度数量，但若仅保留少数模态，则可能遗漏局部响应；而若包含过多高阶模态，则又会增加计算成本并提升数值刚度。残差合规性虽可用于描述因模态截断而产生的准静态效应，但将其纳入锥互补公式的研究相对较少。若能实现这种整合，就能在保证不可穿透性和库仑摩擦特性的同时，避免人为引入接触刚度，从而体现截断模态的局部灵活性。其次，面积较大的潜在接触区域会产生庞大的界面数据集，尤其是那些包含点对点耦合信息的残差合规性矩阵，直接离散存储会带来很大负担。因此，需要一种紧凑的参数域表示方法来整理和压缩接触几何结构、模态场以及界面残差合规性数据，而从离散表面网格构建稳定且实用的参数化方法仍然具有挑战性。第三，一些传统方法会在整个潜在接触区域内，于原始有限元表面网格上预先定义界面位移自由度。当该区域范围较广时，就需要引入大量变量，而且表面表示与接触点分布仍受底层网格的限制。

为应对这些挑战，本文提出的方法整合了以下几项创新技术。首先，将因模态截断而产生的残差合规性纳入锥互补公式中，从而将保留下来的模态所对应的整体动力学行为与局部准静态响应联系起来。在此框架下，无需人为引入接触刚度，即可实现不可穿透性和库仑摩擦约束，且在适当条件下，该公式还可转化为凸优化问题。其次，采用边界首次展平方法在统一的参数域中表示接触几何结构、模态场以及界面残差合规性。切比雪夫拟合方法能为这些数据提供紧凑的离线表示形式，并便于在线查询；此外，通过互反性处理及系数阈值筛选，还可进一步减少残差合规性数据的存储量。第三，在接触区域构建以及局部接触信息整合过程中，采用参数化的平滑表面表示方法。通过对拟合后的数据进行在线评估，可方便地在已识别的区域内生成或调整接触点及局部网格，从而降低对原始有限元网格的依赖。再加上所采用的自由界面降阶表示方法，更有助于避免在整个潜在接触区域内预先定义大量的界面位移自由度。

本文的其余部分结构如下：第2节介绍了结合残差合规性修正的柔性体接触锥互补公式，推导了时间离散化的系统动力学方程，并概述了整体求解流程。在此基础上，第3节阐述了残差合规性矩阵的计算方法，以及如何通过边界首次展平映射实现这些矩阵的离线拟合与压缩存储。第4节则探讨了如何利用拟合后的几何信息进行接触检测，包括接触区域的配对生成及相关接触量的计算。在具备上述各项组件之后，第5节通过典型算例展示了所提框架的准确性与有效性。最后，第6节对全文进行总结，并指出未来的研究方向。

章节节选
具有残差合规性修正的锥互补接触动力学
为构建本文所采用的接触动力学模型，首先介绍了基于真实摩擦锥的经典锥互补公式，随后引入残差合规性修正机制，用以补偿局部接触变形，同时仅保留少量模态用于描述整体动态响应。结合时间离散化的系统动力学方程，即可得到每个时间步长的接触力学公式以及整体的求解流程。

残差合规性的计算与参数化存储
为高效地表示和存储界面残差合规性数据，首先需计算离散形式的残差合规性矩阵。接着运用边界首次展平方法对潜在接触面进行参数化处理，建立统一的参数域表示形式。在此基础上，对该参数域中的残差合规性矩阵进行多项式拟合，再通过系数压缩处理，得到可用于在线查询的矩阵函数表示形式。

参数化表面上的接触检测与局部接触信息获取
基于第3节中介绍的参数化表示方法与拟合存储方案，本节阐述了如何为锥互补公式构建接触区域表示形式以及获取局部接触信息。这些局部量是通过逆映射、主从投影关系、局部坐标变换，以及在对检测到的接触区域进行整合后获得的。

数值算例
本章通过六个算例来评估所提方法的准确性及数值稳定性。除非另有说明，所有算例的时间步长均设为1×10??秒。第5.1至5.4节中的所有ABAQUS参考模型均采用四节点线性四面体单元（C3D4）进行离散化，各算例对应的网格尺寸及网格统计信息也在文中给出。对于涉及节点力对比的算例，还给出了归一化绝对误差ei(Q)以及区域……

结论
本文提出了一种用于处理柔性多体系统非光滑摩擦接触动力学的计算方法。在该方法中，不可穿透性和库仑摩擦约束以锥互补形式加以体现，同时引入界面残差合规性矩阵来修正局部接触变形。通过这种方式，整体动态响应得以以基于模态的降阶形式呈现，而接触力与局部变形之间的耦合关系则……