基于图论与螺旋理论的机器人效率分析

机器效率通过图论和螺旋理论确定并应用于机器人技术

摘要

提出了一种基于图论与螺旋理论确定复杂机构（包括并联机器人）效率的新方法。该方法可用于任意类型关节以任意方式连接的机构。图论用于描述拓扑结构，螺旋理论用于描述关节几何特征。该方法是戴维斯方法的扩展，其中使用主动耦合来建模功率输入与输出运动链。通过拆分棱柱副，获得更高的精度，从而得到更真实的模型。文中以滑块曲柄机构和3‐UPU空间操作器为例，详细介绍了棱柱关节、旋转关节和万向节的建模方法。

1. 引言

通常，在机器人技术领域，效率被认为不太重要，相关文献也较为匮乏。造成这一现象的一个原因是典型的串联机器人所需的功率相对较低。然而，并联机器人较高的负载能力正开始改变这种认知。克雷格 [1]估计，平均约有25%的电机扭矩大小用于克服关节摩擦。并联机器人通常具有更多的关节，且由于某些应用涉及大量能量，效率可能起到关键作用。此外，并联机器人中普遍存在过约束[2] ，这会进一步增加摩擦。

准确的损耗模型可用于在机器人制造前评估其效率。因此，效率可作为选择不同并联机器人结构时的决策标准。对于给定的结构，效率计算能够检测哪些关节消耗更多能量，从而以最低成本优化能耗。

摩擦模型可分为两部分：摩擦作用大小和摩擦作用的几何支撑。作为摩擦建模的先决条件，摩擦作用大小在文献[3]中已有充分论述。而几何因素及其对摩擦作用大小的影响尚未得到应有的关注，如第3节所示。

许多研究并联机器人效率的作者承认摩擦是一个重要因素，但在评估模型[4,5] 或考虑传动中的摩擦[6–8] 时忽略其影响。这对串联机器人可能足够，但对并联机构未必适用。杜邦 [6]指出，机器人部件中的法向力随关节位置、速度和加速度而变化。然而，在他的棱柱关节模型中，杜邦 [6]将损耗归因于螺杆与螺母螺纹之间的摩擦。滑动副之间的摩擦被忽略。杜邦 [6],还指出，使用简单的综合摩擦模型来处理每个关节所有可能的摩擦效应并不一定正确。

蒂施勒等人 [9]使用公式(3)对多环连杆机构中的摩擦进行建模。使用了运动图 [10,11] ，并通过向图中添加路径来内化外部作用。这种选择不必要地增加了变量数量。蒂施勒等人 [9]将摩擦作用的大小视为独立量，通过将传递作用于联轴器的作用大小与摩擦系数线性相乘得到。该过程迭代直至解收敛。似乎蒂施勒等人 [9]未考虑摩擦作用可能具有高度非线性，如第3节所述。显然，在该摩擦模型中忽略了联轴器传递的扭矩。

弗拉先科 [12]开发并实现了机械系统的分布式仿真方法。据称该结果是一种精确的非迭代算法，适用于具有任意关节类型和任意拓扑结构的机构，包括分支和运动学闭环。提供了两个示例：一个串联机器人和一个汽车悬架。由于未详细说明摩擦模型，弗拉先科 [12]很可能也忽略了第3节中讨论的摩擦模式。

法尔哈特等人 [13] 研究了建模过程产生的非线性微分方程组，即刚性微分方程的解。所使用的摩擦模型由同一研究小组在先前的一篇论文中提出 [7] ，且仅在执行器中包含摩擦，类似于杜邦的 [6] 模型。

本文提出了一种计算包括并联机器人在内的复杂机械效率的通用算法。该方法具有很强的通用性，可适用于任意运动链和任意类型的关节。采用图论描述拓扑结构，利用螺旋理论建模关节的几何特征。选用库仑摩擦模型，因其简单且具有非线性特性，但也可采用其他摩擦模型。

在下一节中，将回顾戴维斯方程和螺旋符号。在第3节中，将基于作用螺旋定义摩擦模式与模型，并讨论对棱柱、旋转关节和万向节建模时的考虑因素。第4节给出了第一个示例，主要目标有两个：引入符号表示，并提供对方法的理解。在第 5节 ，主要示例被提出。在第 5.5节中 做了一些评述 ，并在第 6节中得出结论。

2. 戴维斯方程

本文的表述紧密遵循[14] 。熟悉戴维斯方程的读者可考虑跳过本节。

将基尔霍夫定律应用于多体系统的基础是通过图来表示耦合网络，该图称为耦合图，其中每个连杆（刚体）由一个节点表示，连杆之间的每个直接耦合由一条边表示。因此，边变量建模了通过耦合传递的作用力（through variable）以及耦合允许的运动（ across variable）。这些变量在两个方面不同于其电学对应量：第一，每个运动和每个作用力在几何上都是一个螺旋，需要 d 个坐标， 其中 1 ≤ d ≤ 6；第二，一个耦合器最多可传递/允许 d 个独立作用力/运动。对于由连杆和耦合器构成的网络，Davies [15] 提供了相当于基尔霍夫电压定律和电流定律的方程。Davies [10,11] 还给出了另外两个基于虚功率的方程。这四个方程在表1中重现，它们提供了连杆对可能经历的运动大小（一阶运动学）之间的关系，这些运动源于欠约束，以及可能存在的作用力大小之间的关系，这些作用力源于过约束。在第4.2、4.3、5.3.1和5.3.2节中，详细举例说明了公式(1)和(2)的应用。

方程(1)的解类似于仅包含互连电压源的电气网络的网络方程的解。同样，方程(2)的解类似于仅包含互连电流源的电气网络的网络方程的解。据作者们所知，机械网络方程(1–4)首次引入了电阻力的机械类比，用于建模参考文献[14]中由摩擦引起的功率损耗。

2.1. 螺旋表示法与功率

本文再次采用Davies引入的符号[10,11,16,17] ，并做了一些扩展。所考虑的所有运动螺旋和作用螺旋所属的螺旋系统的最小阶数为维度d（1 ≤ d ≤ 6）；k 和 分别为独立割集和回路（环路）的数量；C 和 F 分别为约束自由度和自由度的总量，即所有耦合器的约束自由度c和自由度f之和。矩阵ˆ M N 的每一列包含一个耦合器特征的f个独立单位运动螺旋中某个运动螺旋的运动螺旋坐标，该f个独立单位运动螺旋张成运动螺旋的f‐系统。同样，矩阵ˆ A N 的每一列包含c个独立单位作用螺旋中某个作用螺旋的作用螺旋坐标。向量ψ和分别包含运动螺旋和作用螺旋的大小。

一个力向量 Q ={ U, V, W }和一个在原点处的力矩矢量 P ={ R, S, T }表示一个作用螺旋 $ a 。一个角速度向量 p ={ r, s, t} 和一个位于连杆上瞬时处于原点处点的线速度向量 q ={ u, v , w} 表示一个运动螺旋 $ m 。作用螺旋 $ a ={ R, S, T ; U, V, W} 被称为按轴坐标顺序书写，而运动螺旋 $ m ={ r, s, t; u, v , w} 则称为按射线坐标顺序书写。螺旋的前一部分与第二部分之间的关系称为螺杆螺距 h，其单位为长度。纯力和纯角速度的螺杆螺距为零，h = 0 。纯扭矩和纯线速度的螺杆螺距为无限大，h → ∞ 。扭矩与力的组合称为力系（wrench）。同样，线速度与角速度的组合称为旋量（twist）。如果螺旋螺距是有限的，则瞬时螺旋轴（ISA） 的方向由运动螺旋中的 p 的单位向量 ˆ p 给出，或由作用螺旋中 Q 的单位向量 ˆ Q 给出，且该瞬时螺旋轴的位置可由这两个向量确定。然而，如果螺距是无限的，则其方向由 ˆ q 或 ˆ P 给出，位置则无法确定。具有有限螺距的作用螺旋的瞬时螺旋轴通常称为作用线。

有时将螺旋表示为一个大小乘以一个归一化的螺旋是方便的，即$= ˆ $ ，其中ϱ是大小，ˆ $ 是归一化螺旋（一种纯几何实体）。根据需要，会添加上标a（表示作用）和m（表示运动）。大小ϱ是一个实数，对于作用螺旋，其量纲为力或扭矩（仅当作用螺旋为纯扭矩时，h→ ∞ ）；对于运动螺旋，其量纲为角速度或平移速度（仅当运动为纯平移时，h→ ∞ ）。通常，与某一根正则轴对齐的特定螺旋的大小用坐标标签表示。例如，通过耦合器C传递且与x轴平行的纯力的大小用U C 表示，而不是 a C。这种做法有助于快速识别螺旋的几何支撑。

在[11],中解释了引入这种符号的原因之一，即作用螺旋 $ a 对运动螺旋 $ m 所消耗的功率 P 可以方便地记为内积：
$$P=$ a •$ m= r R+ s S+ t T+ u U+ v V+ w W. (5)$$
如果能量瞬时进入联轴器网络，则该产品将为正标量；如果能量离开，则为负值。

相同的表达式，使用亨特对普吕克线坐标[18]的改进，为：
$$L 1 P 2∗+ M 1 Q 2∗+ N 1 R 2∗+ L 2 P 1∗+ M 2 Q 1∗+ N 2 R 1∗$$
其中，下标1和2可不加区分地表示作用或运动。

3. 摩擦模式与模型

在刚体假设下，可以认为任何力都集中在其相应作用线上的一点。但在实际中，这种集中力并不存在，所有通过机械方式施加到物体上的外力都是分布在有限接触面积上的。摩擦的考虑需要在这两个概念之间进行协调。出于教学目的，本节假设为平面几何。该假设意味着作用螺旋属于第四类特殊三系螺旋系统，即平面静力学螺旋 [18] 。

当外力 Q作用于在导轨内以速度 q移动的滑块时，如图1和图2所示，可能会产生摩擦力 Q F。该摩擦力取决于速度 q、法向反力 Q N,材料以及接触表面。在图1 a中，外力 Q和法向反力 Q N的作用线通过接触区域中心。反应力 Q N在接触区域上按照某种分布规律分布；通常对于刚体接触，此规律被假定为均匀的。

在图1 b中， Q的作用线不通过接触区域中心，但仍与该区域相交。因此，与图1 a中的情况相比，法向反力分布在较小的面积上。对于依赖于接触区域的摩擦模型，必须考虑这一事实。

另一方面，如果摩擦模型不依赖于接触区域，则图1所示的两种情况被视为等效。在这两种情况下，摩擦力 Q F的作用线位于接触区域内，并与法向反力 Q N的作用线相交。此时，摩擦力 Q F的大小不依赖于外力 Q的作用位置，只要该外力的作用线与接触区域相交即可。

图1中的力分布有效，直到外力 Q的作用线移出接触区域。当外力 Q的作用线到达接触区域的边缘时，有效接触区域变为一个点。在刚体假设下，这不会对摩擦力模型产生显著影响，除非该模型依赖于接触区域。然而，当外力 Q的作用线位置超出接触区域边缘时， 将对摩擦力模型产生更严重的影响。

在图2 a中，外力 Q的作用线通过接触区域边缘的右侧。如果忽略导轨与滑块之间的间隙，或者仅允许少量柔性时，法向反力将分布在导轨上下两侧的两个区域。通常假设力分布遵循三角形分布规律，因此法向反力 Q N ′和 Q N ′′的位置是已知的。法向反力 Q N ′和 Q N ′′ 的大小可根据外力 Q、滑块长度以及外力 Q在滑块中心产生的力矩来计算。一旦确定了法向反力 Q N ′和 Q N ′′ ，便可基于相同的摩擦模型计算摩擦力 Q F ′和 Q F′ ′。

本工作所采用的网络方法的一个优点是，系统中传递的作用以螺旋形式表示，即力向量以及该力在原点处产生的力矩向量。通过坐标变换，可以在任意感兴趣点计算该力矩矢量。

如果导轨与滑块之间的间隙不可忽略，滑块将发生微小旋转，并形成两个接触点，如图2 b所示。假设间隙很小，因此旋转产生的唯一显著影响是形成了接触点，那么图2 a和图2 b所述情况之间唯一明显的区别在于法向反力 Q N ′和 Q N ′′ 的作用位置。例如， 若采用库仑摩擦模型，只需将摩擦系数简单地乘以3/2，即可使两种方法得到的结果相容[19] 。

已确定两种工作模式。第一种模式如图1所示，只有一个法向反力，其作用线与外力的作用线相同。第二种模式如图2所示。在此模式中，存在两个法向反力，其作用线的位置与外力的作用位置无关，但大小取决于外力的作用位置以及滑块宽度。

这两种模式之间的区分基于外力的作用位置。有时，这种区分可以预先确定。但在大多数情况下，两种模式的可能性会使建模过程变得不必要地复杂。如果可能，应进行简化以避免这些复杂性。一种常见的简化因素是接触区域的尺寸相对较小。另一种可能性是将单个联轴器分解为两个或多个并联作用的耦合器（例如，参见第3.2节中的图4）。

3.1. 接触表面几何

图1和图2所示滑块的横截面几何形状也会影响法向反力分布，进而影响法向反力和摩擦力的大小。通常，文献中提供的摩擦系数仅适用于平面。图3展示了外加力 Q通过圆柱轴线时的法向反力分布情况。销钉与衬套的啮合角用2 β表示。

通常的做法是根据法向反力分布规律调整摩擦系数。对于圆柱表面，Baranov [19]推荐使用两种分布规律：新副采用常数分布； 旧副因已磨损并相互贴合，采用余弦分布。本文假设完全接触、2 β= πrad、旧副；因此，当接触表面为圆柱面时，适当的摩擦系数可通过将平面情况的摩擦系数乘以 4 π 得到。

3.2. 棱柱和圆柱副关节模型

Davies [11]使用被动圆柱副等多种类型来说明公式(1)至(4)的应用。在静力学分析中，这些副被视为能够传递四个独立作用力： 两个力和两个扭矩。只要这些力的作用线与副轴线相交，其作用位置可视为任意，因为在求解公式(2)或(3)时，扭矩会补偿由于作用位置不同而产生的差异。

将这一推理扩展到图1和图2的平面情况，法向反力由一个平行于 Q N（或 Q N ′和 Q N ′′ ）但位置任意的单个力以及一个扭矩来表示。如果已知等效力和扭矩，则可以判断工作模式，并相应地施加正确的法向反力。然而，在网络问题中，这些等效作用依赖于摩擦力，而摩擦力又依赖于相应的法向反力。

为了解决这一难题，可将滑块拆分为两个小型并联滑块，如图4所示。每个滑块能够传递一个独立的力。这些力对应于图2中的法向反力，或等效于图1中的单一法向反力。

使用分体式滑块时无需确定工作模式。然而，摩擦力的垂直位置取决于相应法向反力的方向（比较图4 a和b中的 Q F ′ ）。如果滑块高度较小，则可忽略此因素；否则，摩擦力作用位置取决于对应法向反力大小的符号，即联轴器传递的力的符号。因此，公式 (2)和(3)变为非线性。

如果系统相对较小，仅包含少数符号函数，则有可能像第4节中那样找到符号解。然而，对于较大的问题，只能得到数值解。第 5节讨论了数值解。

将滑块拆分为两部分具有第二个优点。通常，棱柱副或圆柱副的结构设计使得分体滑块之间的距离不是常数分布，如图5所示。与其考虑一个能够传递两个垂直力且其位置难以确定的单个联轴器，不如更实用地考虑两个并联联轴器，每个仅传递一个易于确定位置的垂直力。

然而，在考虑运动副传递的其他作用力时应谨慎。例如，对于移动副，绕运动副轴线的扭矩可视为由任意一个并联耦合器传递， 但不能同时由两个耦合器传递。如果两个耦合器都能传递扭矩，则可能会在由并联耦合器形成的闭合回路中锁定一个不确定扭矩 [16,17] 。该扭矩将构成一种新的过约束，这既不必要（因为它并不能改进对系统的描述），也不理想（因为它需要引入额外的主变量）。

第5节给出了使用移动副的一个示例。

3.3. 转动关节模型

阻碍转动关节运动的作用力是一个平行于关节轴的纯扭矩。通常，这种关节可以由任何回转曲面构成。本文假设采用圆柱销结构。在此几何条件下，形成两个接触表面：一个为圆柱面，另一个则较不明显，为平面副。关节传递的力包括作用于这些表面上的法向力， 当关节运动时，会产生与中心销半径成正比的摩擦扭矩。

如第 3.1节所 述 ，需对作用于圆柱表面的法向力进行摩擦系数的调整。

对于平面副表面，分析稍微复杂一些，但对于旧的销钉而言，情况会简化，且无需调整。 有兴趣的读者可参考[20,21] 。

关节传递的扭矩与关节轴正交，也会在接触表面产生力分布，并导致部分摩擦扭矩。然而，如果中心销钉较短，则这种效应可以忽略不计（详见[21] ）。

旋转关节可以 以类似于棱柱关节和圆柱关节的方式分解（见第3.2节）。优势 is 以避免未知因素的影响力的位置。

3.4. 万向或虎克关节

万向节或胡克铰传递三个独立力和一个独立扭矩[22] 。当其运动时，在存在摩擦的情况下，会产生两个附加的摩擦扭矩。由于这些扭矩是正交的，因此可以按照第3.3节中讨论的方式进行建模。

4. 滑块曲柄机构

对于效率分析，图6所示的滑块曲柄机构可被视为由四个刚性连杆组成，这些连杆通过三个转动关节和一个圆柱副或棱柱关节连接。此处，该后者副被视为棱柱副，因为运动链中其余构件完全限制了D副（见图6）的旋转。

在刚性连杆假设下。将联轴器 D视为棱柱副会降低维度 d ，且这更多地是一种便利性而非必要性而非必需而非一种必要性。假设为理想连杆几何，因此忽略与过约束相关的z轴方向的力以及平行于z = 0平面的扭矩作用。为了简化示例，仅对耦合器D中的摩擦进行建模。

通常情况下，像耦合器D这样的棱柱副的位置是任意的。然而，从摩擦考虑出发，通过仔细放置此类副可以实现简化。一个沿x轴方向的外力和一个平行于x轴的摩擦力作用在耦合器D上。摩擦力的实际位置取决于D上垂直反作用力的方向。目前，忽略摩擦力的垂直偏移。因此，将耦合器D与C重合放置较为方便，这样耦合器C和D中的垂直力在原点产生的扭矩可以相互补偿。此外，耦合器D的行为类似于图1中的滑块，且仅能发生第3节中讨论的第一种操作模式。

联轴器B、C和D的位置可以用θ、曲柄角度以及杆长a1和a2表示为：
$$B x= a 1 cos θ$$
$$B y = a 1 sin θ$$
$$Cx = D x= B x+ √ a 2 2 − a 1 2+ B x 2$$
$$= a 1 cos θ+ √ a 2 2 − a 1 2+ a 1 2 cos 2 θ. (6)$$

4.1. 螺旋系统与摩擦模型

图6所示耦合网络的运动螺旋根据亨特的[18]分类属于第五类特殊三系螺旋系统。通过将z轴定向为与标记为A、B和C的耦合器的瞬时螺旋轴平行，x轴定向为与耦合器D的瞬时螺旋轴平行，运动可由三个运动螺旋{ t ; u , v }张成。与各耦合器相关的运动螺旋详见 表2。

图6所示耦合网络中的作用螺旋，可由{T; U, V}张成，属于第四类特殊三系螺旋系统，即平面静力学螺旋[18] 。这些作用螺旋在 表3中进行了描述。注意，在表3中，t A是施加于该网络的外力扭矩的大小，u D是施加于该网络的外力的大小。u D D 是

表 2 运动张成各系统的运动允许由图联轴器在图 6中

联标签轴，器类型	平面位置	旋转或平移	速度（方向）	大小	单位运动螺旋坐标 按射线坐标顺序
	x	y	角速度	原点处点的速度
			t	u	v
A，旋转关节	0	0	旋转 ( z)	tA	1 0 0
B，旋转关节	Bx	By	旋转 ( z)	t B	1 By −Bx
C，旋转关节	Cx	0	旋转 ( z)	t C	1 0 −Cx
D，棱柱关节	Cx	0	Transl ( x)	u D	0 1 0

表 3 作用跨越各系统的作用可由图 6中的耦合器传递的图 6 。

联轴器 标签，类型	平面位置	力或扭矩 （方向）	大小	在…处的力矩 原点	单位作用螺旋坐标 按轴坐标顺序
	x	y	T	U	V
			扭矩( z )
A，旋转关节	0	0	力( x )	UA	0 1 0
			力 ( y)	VA	0 0 1
B，旋转关节	Bx	By	力 ( x)	U B	−By 1 0
			力 ( y)	V B	Bx 0 1
C，旋转关节	C x	0	力 ( x)	U C	0 1 0
			力 ( y)	V C	Cx 0 1
D，棱柱关节	Cx	0	扭矩 ( z)	T D	1 0 0
			力 ( x)	U D	0 1 0
			力 ( y)	V D	Cx 0 1
			力 ( x)	U D D	0 1 0

摩擦力的大小由以下给出
$$U D D= −μD sign( V D) V D sign( u D) (7)$$
其中，V D 是通过耦合器D传递的垂直力（法向力）的大小，μD是经典意义上的无量纲摩擦系数，u D 是耦合器D的平移速度大小（滑块相对于固定坐标系的速度）。式(7)中的负号是由于螺旋$a UD D 和$a uD 的方向均从左向右所致。式(7)实际上是库仑摩擦模型， 但也可以采用更复杂的摩擦模型。

4.2. 运动分析

图6所示单环运动链机构的运动分析是一个简单的问题，因为网络单元运动矩阵 ˆ M N 等于由表2得到的直接耦合 ˆ M D 的单位运动矩阵。
$$Mˆ N= Mˆ D=[ tA t B t C u D t 1 1 1 0 u 0 B y 0 1 v 0 −Bx −Cx . (8)$$
其中 B x 、B y 和 C x 由公式(6) 给出。ˆ M N 和 ˆ M D 的零空间，也称为核，是
$$Null(M ˆ N)= ⎡⎢⎣ tA Cx − Bx t B −Cx t C Bx u D B y Cx ⎤⎥⎦ 1 C x −Bx . (9)$$
方程(1)的解与Null ( ˆ M N) 的唯一向量成正比。表2中列出的任何变量均可选为主变量。任意选择t A，其运动螺旋矩阵每列包含 一个螺旋，表示为
$$M= Mˆ D diag(Null( Mˆ N)) tA (10)$$
其中，利用 Eqs.(8)和(9) ，可得
$$M=[ $ m A $ m B $ m C $ m D t Cx − B x −Cx B x 0 u 0 −B y Cx 0 B y Cx v 0 B x Cx −B x Cx t A C x − B x . (11)$$

4.3. 作用分析

图6中示意性显示的耦合网络的拓扑结构由图7的耦合图表示。通过将边A、B和C作为分支，可直接观察得到G C 的基本割集矩阵。
$$Q=[ A B C D A 1 0 0 B 0 1 0 −1 C 0 0 1 . (12)$$
在这种情况下，任何支路选择都会得到相同的割集矩阵。有关基本割集矩阵的获取方法，请参考文献 [10,11,15,23–25] 中的任意一篇。

作用图的基本割集矩阵通过 G A的列复制来获得 ，如下
$$Q A= ⎡⎣ A TA A UA A VA B U B B V B C U C C V C D T D D U D D V D D U D D A 1 1 1 0 0 0 0 −1 −1 −1 −1 B 0 0 0 1 1 0 0 −1 −1 −1 −1 C 0 0 0 0 0 1 1 −1 −1 −1 −1 ⎤⎦ (13)$$
其中，Q 的每一列被复制的次数等于相应联轴器传递的独立作用力的数量（见表3）。或者，可以参考文献[11],绘制作用图G A， 然后通过观察组装矩阵Q A。

直接耦合器的作用矩阵根据表3组装如下
$$A D=[TA UA VA U B V B U C V C T D U D V D U D D T 1 0 0 −By Bx 0 Cx 1 0 Cx 0 U 0 1 0 1 0 1 0 0 1 0 1 V 0 0 1 0 1 0 1 0 0 1 . (14)$$
网络单元作用矩阵由公式(13)和(14)组装而成
$$[Aˆ N]d k ×C= ⎡⎢⎢⎢⎣ [Aˆ D]d×C[ Q 1] C×C [Aˆ D]d×C[ Q 2] C×C … [Aˆ D]d×C[ Q k] C×C ⎤⎥⎥⎥⎦ d k ×C (15)$$
其中 [ Q i] C×C = diag ([ Q A] i) 是对角矩阵 ( i= 1, 2 , . . . , k) 。[ Q i ] C × C 的对角元素是 [ Q A ] k × C 的第i行的元素。因此
$$Aˆ N= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ TA UA VA U B V B U C V C T D U D V D U D D T 1 0 0 0 0 0 0 −1 0 −Cx 0 U 0 1 0 0 0 0 0 0 −1 0 −1 V 0 0 1 0 0 0 0 0 0 −1 0 T 0 0 0 −B y Bx 0 0 −1 0 −Cx 0 U 0 0 0 1 0 0 0 0 −1 0 −1 V 0 0 0 0 1 0 0 0 0 −1 0 T 0 0 0 0 0 0 Cx −1 0 −Cx 0 U 0 0 0 0 0 1 0 0 −1 0 −1 V 0 0 0 0 0 0 1 0 0 −1 0 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ . (16)$$
式(16) 仅包含拓扑和几何信息；不包含元件特性。表达元件特性的方程称为本构方程 [26] 。例如，该系统的单一本构方程是 式 (7)。

用于对网络中作用关系进行建模的增广作用矩阵 augmented action matrix D A 通过将方程Eq.(7)以矩阵形式代入方程Eq.(16)得到
$$D A = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ TA U A VA U B V B U C V C T D U D V D U D D T 1 0 0 0 0 0 0 −1 0 −Cx 0 U 0 1 0 0 0 0 0 0 −1 0 −1 V 0 0 1 0 0 0 0 0 0 −1 0 T 0 0 0 −B y B x 0 0 −1 0 −Cx 0 U 0 0 0 1 0 0 0 0 −1 0 −1 V 0 0 0 0 1 0 0 0 0 −1 0 T 0 0 0 0 0 0 Cx −1 0 −Cx 0 U 0 0 0 0 0 1 0 0 −1 0 −1 V 0 0 0 0 0 0 1 0 0 −1 0 c.e. 0 0 0 0 0 0 0 0 0 a μ D 1 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (17)$$
其中 a = sign( V D) sign( u D) , 并且标签 c.e. 表示最后一行来自一个构成本方程。

该矩阵用于替代方程(2)中的 ˆ A N ，且其零空间为
$$Null( DA )= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ TA By Cx UA Bx −Cx VA By U B Bx −Cx V B By T C Bx −Cx T C By T D 0 U D Bx −Cx + a μD By V D By U D D −a μD By ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ 1 By Cx (18)$$
其中，再次， a= sign( V D) sign( u D) 。

作用幅值向量与Null ( D A) 中的唯一向量成正比。因此，在任何情况下，由于耦合器D的位置被设置为与C重合，T D为零。

如果耦合器D位于其他位置，则T D可能不为零，且V D将具有不同的值，需要对式(7)进行修改。

任何幅值变量都可以被选为主变量，但T D除外，因为它在方程(18)中为零。选择t A时，作用幅值的向量由以下给出
$$= Null( DA ) TA (19)$$
其中 Null( D A)由方程(18)给出。

可以使用另一个变量，而不是 T D ，作为主变量，只需相应地调整方程 (19) 即可。例如，如果希望使用 V D 作为主变量，则将 t A 替换为 C x V D，且向量 Null( D A) C x V D 的最后一个元素将仅为 V D。

每列包含一个作用螺旋的作用矩阵 ，该矩阵可通过以下方式计算： ，可以被计算为
$$A= ˆ A D diag( ) = ⎡⎣ $ a TA $ a UA $ a VA $ a U B $ a V B $ a U C $ a V C $ a T D $ a U D $ a V D $ a U D D T By Cx 0 0 By ( Cx −Bx ) Bx By 0 By Cx 0 0 By Cx 0 U 0 Bx −Cx 0 Bx −Cx 0 Bx −Cx 0 0 Bx −Cx + a μD By 0 −a μD By V 0 0 By 0

5. 自由度机构： 3‐UPU

本节致力于对具有并联运动链的复杂机械进行摩擦分析，旨在确定其效率。所选作为示例的机械是一种称为3‐UPU的并联机器人。
3‐UPU [27]是一种三自由度空间并联操作手，由三个支链组成，支链通过万向节或胡克铰连接到基座和平台，如图8所示。标签 与图9和图10相关，并将在后文进一步说明。每个支链包含一个驱动棱柱副，并连接两个万向节。假设各支链相同，且它们与基座和平 台的连接呈对称分布（万向节中心形成两个等边三角形，一个位于基座，另一个位于平台，与图8所示三角形平行）。所提出的方法原 则上可应用于非对称机构，但采用对称结构更为方便，因为它可以实现更简洁的描述。此外，系统被假定为准静态，因此忽略惯性效 应和控制误差。

5.1. 逆运动学

评估3‐UPU的整体效率需要计算在执行任务期间平台的运动和作用。通常，任务由携带负载的平台所期望的运动规定。由于3‐UPU 仅能进行平移，因此该任务的描述等同于平台上任意一点运动的描述。注意，在奇异条件下可能发生平台旋转[28], ，但在正常操作 下可避免奇异性。平台中心的点标记为点M，并在点M处对平台施加外部载荷（见图11和图13）。

表示连杆如何连接以构成整个3‐UPU的耦合图如图9所示。在运动链中包含了一个称为阿苏尔虚拟链[29,30] , 的虚拟耦合，也标 记为耦合M，它直接耦合基座和平 台。该额外耦合的作用有两个：一是允许使用与平台相关的变量作为主运动变量（逆运动学）；二 是将作用于平台的外部作用（逆动力学）即外部载荷内化。耦合M在图9中由虚线表示，具有六个自由度且 一个约束度 。因此，联轴器 M是一种主动式联轴器 ，从这个意义上说，功率可能被消耗在它上面或由它提供 。

5.1.1. 联轴器位置和方向

运动和作用分析所需的附加细节包括耦合器的位置以及轴的方位，所有这些都相对于一个参考坐标系明确写出。

参考坐标系的原点位于与基座相连的万向节中心所构成的等边三角形中心处，z轴垂直于基座平面，x轴指向一个万向节中心（见 图11和图13）。该特定万向节属于所谓的第一支链。由于所有支链结构相同，仅位置不同，因此本文仅描述此第一支链，其示意图如 图10所示。在图10中，标记为A的万向类型联轴器将标记为0的基座连接到连杆1。连杆1通过棱柱关节与连杆2相连，该棱柱关节由联 轴器B和C建模（更多细节见第5.2节）。最后，连杆2通过另一个标记为D的万向节连接到连杆3，即平台。

点 A位于万向节中心，且标签相同。连杆 2在点 B处连接到连杆 1 ，且 A与 B之间的距离为 常数分布 。对于点 C和 D，存在类似情况。

各支链中与基座相连的万向节中心（例如，图10中的点A）相对于参考坐标系的位置由位置向量给出
$$e 1 =[ x cos α y sin α z 0] r b (35)$$
其中 r b 是从基座中心（原点）到任意与基座相连的万向节中心的距离，α=0 , 2 π/ 3, −2 π/ 3 rad 分别对应第一、第二和第 三支链。由于平台旋转被限制，每个万向节第一轴的方向由单位向量 ˆ e 1 表示 e 1 ，参见图10。

向量点 从基座的万向节中心到每个支链的平台的万向节中心由
$$e 2 =[ x Mx +( r t − r b) cos α y M y+ ( r t − r b) sin α z Mz ] (36)$$
其中，r t 是平台中心到与平台相连的任意万向节中心的距离；M x 、M y 和 M z 是平台中心（点M）相对于参考坐标系的坐标； r b 和 α 定义见公式(35)。值得注意的是，由公式(36)给出的向量依赖于当前平台位置，而公式(35)中的向量为常数分布。各支链棱 柱关节的方向（轴线）由 e 2 中的单位向量ê2 给出。

每个万向节的第二轴的方向由单位向量 ˆ e 3 给出，该向量与 ˆ e 1 和 ê2 正交。联轴器B和C传递平行于 ˆ e 3 和 ê4 的力；ê4 与 ê2 和 ˆ e 3 正交。万向节传递围绕单位向量 ˆ e 5 的扭矩，该向量与 ˆ e 1 和 ˆ e 3 正交。

五个单位向量 ˆ e i ( i = 1, 2 ,… 5) 取决于由 α 给出的相应支链的位置。对于第一支链，α= 0 ，图10 中显示了 ˆ e i。

基座上的万向节中心位置（即点A、E和I）由公式(35)给出，为 e 1，并取合适的角度α。此外，平台上的万向节中心位置（即点 D、H和L）由公式(35)和(36)给出，为 e1 + e 2。忽略棱柱关节的径向长度，连杆1和2之间的两个接触点由
$$B= e 1 + l B eˆ 2 C= e 1 + e 2 − l C eˆ 2 (37)$$
其中 l B是常数分布 dist点A和点B之间的距离，以及点C和点D之间的距离。<方程>(37)</方程>也被用于确定 ne 点 F， G， J，和 K使用 a适用于值的 α。

5.2. 棱柱关节

如第3.2节所述，摩擦考虑要求将棱柱关节拆分为两个联轴器，以简化传递力的布置。因此，每个支链中的棱柱关节被拆分为两个 联轴器B和C（第二和第三支链分别为F‐G和J‐K，见图9）。联轴器B和C必须共同表现为一个棱柱副；即必须具有一个自由度和五个约 束自由度。挑战在于选择两种联轴器类型，当它们并联连接时，仅允许一个平移运动，且不产生任何过约束。这可以通过使用一个圆 形和一个方形环形关节来实现。

圆形环形关节允许三个旋转，其中两个是无穷小的，以及一个平移。它可以被想象成手指上的戒指：可以自由地绕手指轴线旋转， 也可以绕垂直于手指的任意方向旋转，但幅度很小。因此，该关节传递两个力。

表 4 运动张成由系统的运动允许通过耦合器在 图 10中

联轴器大小	运动	方向
A pA 1	旋转	ˆ e1
pA 2	旋转	ˆ e3
B q B	平移	ˆ e2
p B 1	旋转	ˆ e3
p B 2	旋转	ˆ e4
C q C	平移	ˆ e2
p C 1	旋转	ˆ e2
p C 2	旋转	ˆ e3
p C 3	旋转	ˆ e4
D p D 1	旋转	ˆ e1
p D 2	旋转	ˆ e 3

方形环形关节允许两个无穷小的旋转和一个平移。因此，它传递两个力和一个扭矩。将并联连接的两个环形关节的轴线对齐后， 仅允许平移，该组合结构表现为棱柱关节。

需要一个方形环形关节来防止绕关节轴的旋转。然而，不建议使用两个方形环形件，因为在这种情况下，并联子网络中需要锁定 扭矩，这会使作用分析变得不必要的繁琐。

方形环位于B点，圆形环位于C点。将环的类型分配给每个点完全是任意的。

5.3. 分析

第4节中描述的分析是借助计算机代数系统完全通过符号方式完成的。然而，3‐UPU的需求远远超出了现有软件和硬件的能力。 因此，编写了一个用C语言实现的软件程序，用于模拟3‐UPU执行任务的过程。单位运动和直接耦合的作用矩阵 ˆ M D 和 ˆ A D、作 用图 G A 的基本割集矩阵 Q A，以及运动图 G M 的基本回路矩阵 B M 均使用计算机代数系统生成，并被硬编码到软件中。在接下来 的章节中，将描述这些矩阵的获取、摩擦模型以及软件实现。

5.3.1. 运动分析

第一支链耦合器允许的运动与表4中所示的运动线性相关。直接耦合 ˆ M D 的单元运动矩阵的前33列是根据表4中的数据结合公式 (35)到(37)给出的耦合位置组装而成。最后6列与耦合M相关，由以下给出
$$Mˆ D M = ⎡⎢⎢⎢⎢⎢⎣ r M s M t M u M v M w M r 1 0 0 0 0 0 s 0 1 0 0 0 0 t 0 0 1 0 0 0 u 0 −M z M y 1 0 0 v M z 0 −M x 0 1 0 w −M y M x 0 0 0 1 ⎤⎥⎥⎥⎥⎥⎦ .$$
通过检查图9中的耦合图 G C，并以边C、G、H、K、L和M作为弦，构造出耦合图 G C 的基本回路矩阵 B。运动图 G M 的基本 回路矩阵 B M 通过复制 B 的列得到。由矩阵 ˆ M D 和 B M 共得到36个线性方程，包含39个变量。施加点M的速度大小 u M 、v M 和 w M 后，求解该线性方程组以获得其余的运动幅值。

5.3.2. 作用分析

表5 描述了3‐UPU第一支链中耦合器传递的作用力系统所涵盖的各个作用。第一支链的执行器产生大小为 Q C A 的力。摩擦力的 大小分别为 Q B F 和 Q C F ，摩擦扭矩的大小分别为 P A F 1 、P A F 2 、P D F 1 和 P D F 2 ，这些在第5.3.3节中有详细描述。表5 中其余 的十三个作用与第一支链每个耦合器的约束自由度相关。

联轴器A和D能够传递三个独立的力。通常，用于表示这些力的螺旋会与正则轴对齐。然而，所应用的摩擦系数取决于接触表面的 几何和材料特性。因此，构成联轴器A和D所传递的力的力的方向分别垂直于圆柱面和平面。

表 5 作用涵盖由系统的作用所能被传递的耦合器在图 10 中。

联轴器大小	作用	方向
A PA	传递的扭矩	ˆ e5
PA F 1	第一摩擦扭矩	ˆ e1
PA F 2	第二摩擦转矩	ˆ e3
QA 1	传递的力	ˆ e1
QA 2	传递的力	ˆ e3
QA 3	传递的力	ˆ e5
B P B	传递的扭矩	ˆ e2
Q B 1	传递的第一个力	ˆ e3
Q B 2	传递的第二个力	ˆ e4
Q B F	摩擦力	ˆ e2
C Q C 1	传递的第一个力	ˆ e3
Q C 2	传递的第二个力	ˆ e4
Q C A	驱动力	ˆ e2
Q C F	摩擦力	ˆ e2
D P D	传递的扭矩	ˆ e5
P D F 1	第一摩擦扭矩	ˆ e1
P D F 2	第二摩擦转矩	ˆ e3
Q D 1	传递的力	ˆ e1
Q D 2	传递的力	ˆ e3
Q D 3	传递的力	ˆ e 5

直接耦合的单位动作矩阵ˆ A D 的前60列根据表5中的数据以及由公式(35)–(37)给出的耦合位置组装而成。最后一列与耦合M相 关，由以下给出
$$Aˆ D M = ⎡⎢⎢⎢⎢⎢⎣ W M R My S −Mx T 0 U 0 V 0 W 1 ⎤⎥⎥⎥⎥⎥⎦ .$$
通过检查图9中的边A、B、D、E、F、I和J作为分支，可得到耦合图 G C 的基本割集矩阵 Q。动作图 G A 的基本割集矩阵 Q A 可 通过复制 Q 的列获得。由矩阵 ˆ A D 和 Q A 可得到总共42个线性方程和61个变量。

5.3.3. 摩擦模型

本研究的重点是摩擦损失的网络效应。目标是建立一种通用方法，能够轻松适应任何摩擦模型。因此，尽管局部摩擦模型对最终 结果有深远影响，但其重要性在此被弱化，因为所提出的方法可以简单地适配到所有摩擦模型。同时，为了保持示例尽可能简单并维 持模型的合理性，假设采用库仑摩擦。

摩擦扭矩被分配给每个万向节旋转轴。这些扭矩的大小取决于关节传递的作用力以及相对速度的方向。考虑到中心销钉的长度较 小，所传递的作用力是将中心销钉压向衬套的力。接触表面可以是圆柱形或平面，具体取决于力的方向。对于平面，摩擦扭矩 P F 由 以下公式给出
$$P F= μ r Q N (38)$$
其中 r是中心销钉半径 ， Q N是法向力的大小[21] 。对于圆柱面 ，方程 (38)右边需乘以 4/ π[21] 。

联轴器A为万向类型，其第一轴与ˆ e 1 对齐；因此，力Q A 作用于平面，而力Q A 2 和Q A 3 作用于圆柱面。因此，第一轴的摩 擦扭矩的大小可估计为
$$P A F 1 = −μ A r A sign(p A 1 )( | Q A 1 | + 4 π | Q A 2 | + 4 π | Q A 3 | ) (39)$$
其中 r A是 轴的半径。

耦合A的第二个轴与ˆ e 3 对齐；因此Q A 1和Q A 3作用于圆柱面，而Q A2 作用于平面（见表5）。第二个轴的摩擦扭矩大小由下式给出
$$PA F 2= −μA rA sign(pA 2)( 4 π | QA 1 | + | QA 2 | + 4 π | QA 3 | ). (40)$$
公式(39)和(40)也可用于联轴器 D ，只需将索引 A替换为 D 。

第一个分支的棱柱关节上的摩擦力被分解为耦合器B和C。耦合器B上的摩擦力大小由以下公式给出
$$Q B F= −μB sign( q B)( | Q B 1 | + | Q B 2 | ) (41)$$
并且方程(41)也用于耦合器 C，将索引 B替换为 C 。耦合器 C（圆形环）中的摩擦扭矩未被建模，因为其旋转被耦合器 B（方形环）所阻止。

共从公式(39)到(41)得到18个本构方程。由于使用了绝对值算子，这些方程是非线性的（符号函数仅用于已知的运动变量）。结 合从矩阵Q A和ˆ A N, 得到的42个网络线性方程，构成一个60个方程61个变量的系统。

5.3.4. 软件实现

使用C语言编写了一个软件程序，专门用于模拟3‐UPU在执行任务期间的行为，重点在于功率损耗计算。在每个离散时间点，该 软件执行的步骤由算法1描述。

算法 1 3‐UPU效率仿真流程。
1. 计算点M的坐标及其线速度；
2. 组装网络单元运动矩阵 ˆ M N；
3. 求解公式(1) 得 到 ψ；
4. 组装网络单元作用矩阵 ˆ A N；
5. 组装本构方程（见第5.3.3节）并将其附加到 ˆ A N；
6. 求解公式(2) 得到 ；
7. 计算瞬时功率。

在 Algorithm 1中，第3步 通过LU分解执行。第6步 更为复杂 ，且需要求解一个非线性方程组的解 。此处使用的特定解法与在其他 研究中使用的方法不同 ，如[13]所示。

求解包含绝对值算子的线性方程组的传统方法是假设所有变量符号的可能组合，然后求解该线性方程组。当解与假设矛盾时，该 解将被舍弃。应预期会出现多个解。

对于3‐UPU作用分析，需要确定30个变量的符号：Q A 1 , Q A 2 , Q A 2 , Q B 1 , Q B 2 , Q C 1 , Q C 2 , Q D 1 , Q D 2 , Q D 2 , Q E 1 ,Q E 2 , Q E 3 , Q F 1 , Q F 2 , Q G 1 , Q G 2 , Q H 1 , Q H 2 , Q H 3 , Q I 1 , Q I 2 , Q I 3 , Q J 1 , Q J 2 , Q K 1 , Q K 2 , Q L 1 , Q L 2 ,和Q L 3 。采用传统方法，需要 求解60个方程的线性系统2^30次，这在实践中不可行。此外，可以合理预期获得的多个解中会有多于一个有效，但只能采用其中一个。

因此，本文采用了较为非传统的同伦延拓法。

同伦延拓法涉及通过从一个已知问题的解出发，并随着该已知问题被同伦变形（即逐渐变形）为给定问题，数值地求得该给定问 题的解[31] 。

在这种情况下，对于第一个时刻，求解过程首先在忽略摩擦的情况下确定作用幅值，并通过逐步增加摩擦系数来继续，直到找到 解。忽略摩擦会将问题转化为具有平凡解的线性问题。在同伦过程中，在获知30个变量的实际值之前，需要知道它们的符号。这些符 号取自已知问题（忽略摩擦）结束时对应变量的值，并存储在变量中；符号变量只能取{ −1, 0,+1}中的值。如果在过程中有任何变 量改变其符号，则更新符号变量并重新启动该过程。对于下一个时刻，已知问题是先前已完成的问题，即位置逐渐变形到当前位置。

该过程耗时，但肯定比尝试寻找230种可能的解要快。此外，还得到了一个合理的解。在下一节中，将讨论使用该软件获得的结 果。

5.4. 仿真结果

执行特定任务用于评估3‐UPU的效率。第二个任务用于评估软件性能和精度。所采用的两个任务是：对负载进行提升以及通过椭 圆路径进行平移。

5.4.1. 模拟的3-UPU物理描述

两个任务均使用具有以下尺寸的 3‐UPU进行了仿真
r b = 1. 0m r t = 0. 5ml B = l F = l J = 0. 6m l C = l G = l K = 0. 6m
万向节的中心销半径为 r i = 15 mm ，i = A, D, E, H, I 和 L。摩擦系数采用 μi= 0. 11 ，i = A . . . , L，如 Shigley 和 Mischke [32] 对润滑钢对 钢接触所建议的那样。

5.4.2. 用于效率测定3‐UPU效率的任务是在1秒内以恒定速度提升恒定有效载荷。作用在点M的总载荷为W M = 10 0 0 0 N，且为常数分布。点M的参数方 程为

$$M x = 0. 1mM y = 0. 2m M z = 0. 5+ 0. 4 t[ m] (42)$$
其中 t是以秒为单位的时间。对方程(42)关于时间求导，得到
$$u M= 0v M = 0w M = 0. 4 m/ s.$$
选择该轨迹是因为它简单，所做功易于计算，从而验证过程直接明了；更重要的是，在执行过程中无相对运动方向变化。因此， 公式(39)–(41)中的符号函数将保持常数分布。

在图 11中，使用仿真数据绘制了 的初始和最终构型，重点展示了每个耦合M的位置 。

采用的时间步长为 5毫秒。如预期所示，负载处消耗的功率 P M, 为常数分布，如图12 所示。在图12 中，还显示了每个执行器提供的功率以及 总功率。

对总功率进行积分，得到执行器供应的能量为 6 146.83 焦耳。传递给负载（通过耦合M离开网络）的能量为 [4 0 0 0]J，该数值对 应于重力势能的增加量。由这两个数值的比值得出的机械效率计算为 η= 0. 6517 。

注意，η= 0. 6517是与此特定任务相关的效率。在不同任务的情况下，效率可能会有所不同。注意，在图12 中，输入功率是变 化的，但输出是恒定的，因此效率在任务执行过程中是变化的。例如，从相同位置开始，用一半的时间提升相同负载经过一半的距离， 其效率被发现为 η= 0. 7372 。

对于原始任务，棱柱耦合器占损耗的98.55%。在万向关节中，较小的摩擦扭矩与低速度结合导致相对较低的功率损耗。

The effects pro duced by a reversa l of运动是 discusse d in th下一节。

5.4.3. 椭圆路径

第5.4.2节中描述的任务适用于效率测定。然而，运动反向的影响——这是验证该方法正确性的一个重要特征——在该任务中并不存在。因此，为 了评估此类

效果 a第二个任务被使用。该任务包含使负载沿由参数方程描述的椭圆路径移动。
$$M x = 0. 2 cos ( 2 π t) M y = 0. 2 sin ( 2 π t) M z = 0. 2 sin ( 2 π t)+ 0. 7 (43)$$
其中是以米为单位，t 以秒为单位。对时间求导后，方程 (43) 得到
$$u M = −0. 4 π sin ( 2 π t) v M = 0. 4 π cos ( 2 π t) w M = 0. 4 π cos ( 2 π t) (44)$$
单位为米每秒。3‐UPU 的初始构型和椭圆路径如图13所示。在图14中，显示了执行器提供的总功率P CA + P GA + P KA 以及负载消耗 的功率P M。在任务开始时，图14中P M的符号为负值，表示功率正从3‐UPU流出，换句话说，负载正在被提升。

运动方向反转会导致摩擦作用发生突然变化，需要急剧调整功率输入。这些急剧调整在功率图中表现为不连续性或趋势变化。第 一次明显的不连续性出现在 t ≈ 66.43ms，此时基座与第一支链（图9和图10中的连杆1和2）之间的夹角最小。因此，p A 2 和 p D 2 的符号发生变化，导致摩擦扭矩方向反转，并需要执行器力突然增加（Q K A 增加了 4.13%）。在此之前的两个不连续性分别发生在 t ≈ 43.54ms（由 p E 1 和 p H 1 引起）和 t ≈ 47.57ms（由 p I 2 和 p L 2 引起），但由于用于绘制图14的平滑算法，它们几乎不可 见。

例如，当t ≈ 0.33275s时，执行器B‐C的运动方向发生反转，从而引起趋势上的细微变化。

通过对比忽略摩擦的情况来考察执行器驱动力的变化，不仅发现力的幅值更高，而且还出现了不连续性和尖点。图15显示了第一 支链执行器施加的力 Q CA 在两种情况下的表现：考虑摩擦和忽略摩擦。在无摩擦情况下，执行器只需支撑负载重量，且平滑轨迹要 求执行器提供平滑的力。当考虑摩擦时，在 t ≈ 0.33275s 处趋势变化的原因明显表现为驱动力的急剧不连续性。图15还显示了在机器 其他部分由于运动方向反转而引起的小幅不连续性和尖点。

图 14和 15 中所示的行为正是具有非线性摩擦的复杂系统所预期的表现。

5.4.4. 与其他方法的比较

如引言中所述，通常研究人员在摩擦模型中未包含传递转矩，忽略了其对实际接触力幅值[6,7,9,12,13]的影响。为了评估忽略传 递转矩所带来的影响，将各支链中的环形联轴器替换为单个棱柱副联轴器。对于这些棱柱关节，考虑了三个位置：每个支链的两个端 点和中点。位置上的差异通过联轴器传递的扭矩进行补偿，且传递力在所有三个位置必须相同。

摩擦力的大小是通过使用 方程(41)估算得出的，因此，由棱柱副联轴器传递的扭矩被忽略。

在执行第一项任务（见第5.4.2节）且考虑条件相同的情况下，所得效率为 η=0 . 9904 ，与棱柱关节的位置无关。该值接近于1， 因为仅根据传递力计算出的接触力远小于实际值。

正如预期的那样，传递的力在不同的关节位置下具有相同的值，这解释了为什么效率在所有三种情况中都相同。

5.5. 备注

作用在圆形环面上的摩擦力（联轴器C、G和K）可视为与法向力成正比（其大小分别为 Q C 1 和 Q C 2 ），而非这些力大小之和， 这是由于轴对称性所致。因此，摩擦力应由
$$Q C F = −μC sign( q C) √ Q C 1 2+ Q C 2 2 (45)$$
代替而不是由方程(41)给出。使用 方程(45)会引入二次方程到系统中。然而，这些方程比包含绝对值算子的方程更容易求解。

由于同样的原因，作用在万向关节上的摩擦扭矩的估算公式 和 E q s. ( 39) 和 ( 40) 可以进行修改。值得注意的是，这些修改后对整体效率没有观 察到显著影响。

图Fig. 9中的耦合图有 275种可能的生成树。在Sections 5.3.1和 5.3.2中使用了同一棵任意选定的生成树。

6. 结论

本文提出了一种对戴维斯方法的改进，该方法允许功率通过端口或耗散元件进入和离开联轴器网络。

用于模拟摩擦的耗散元件，类似于电气电阻，被包含在联轴器网络中。

有源元件类似于通常向网络提供功率的电流源，但也可以从网络吸收功率。这些元件被包含在耦合网络中，用于建模有时作为输 入、有时作为输出的执行器和外部负载。

该方法在考虑功率损耗时，简化了复杂机械的建模。

在复杂机械中，效率取决于具体的任务。适用于所有情况的总效率概念是无意义的。如果要将效率作为对不同机械进行排序的准 则，则必须将任务作为准则的一部分加以规定（见第5.4节）。

非线性方程组采用同伦延拓法求解。该方法在3‐UPU作用分析的求解中表现出良好的性能。在最坏情况下，仅需六次迭代。数值 解所需时间在可接受范围内：评估200个时间步长约需2分钟（见第5.4节）。显然，计算所需时间高度依赖于运行程序所用的硬件。

本文提到的时间是使用HP G61笔记本电脑获得的。

棱柱副和圆柱副中的摩擦可以通过将关节拆分为两个滑块来直接建模（见第3.2节）。这些滑块可以建模为两个圆形环面用于圆柱 副，或一个圆形环面与一个方形环的组合用于棱柱副（见第5.2节）。如果有关节结构的信息，也可以建立更详细的模型。实现驱动棱 柱副的两种常见方式是使用液压油缸和螺旋。Bonchis 等人报道了液压油缸的摩擦模型[33] ，而杜邦则研究了螺旋的摩擦模型[6] 。

在3‐UPU仿真任务中，棱柱耦合器导致的损耗分别占直线路径的98.55%和椭圆路径的97.71%。试图用滚动轴承替代万向节中的 滑动轴承似乎益处不大（见第5.4节）。所提出方法的一个重要特点是，即使在复杂机构的设计阶段，它也能指出运动链中哪一部分需 要更多关注。

在未来的[19,21]工作中，将探讨其他关节类型（如球面副、平面副以及球杆球组合）的建模。文献中提供的摩擦系数通常适用 于平面情况。对于圆柱和球面几何结构，使用这些摩擦系数需要仔细考虑，相关内容将在未来进行探索。此外，还可研究包含此类传 动的传动系统和串联机器人的效率。