在Port-Hamilton框架下，对Cosserat杆动力学采用混合公式化方法及结构保持离散化处理

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Mixed formulation and structure-preserving discretization of Cosserat rod dynamics in a port-Hamiltonian framework

　　菲利普·L·基农（Philipp L. Kinon）| 西蒙·R·欧格斯特尔（Simon R. Eugster）| 彼得·贝奇（Peter Betsch）
德国卡尔斯鲁厄理工学院（KIT）

**摘要**
本文提出了一种基于能量的建模框架，用于研究在较大位移和旋转作用下的空间Cosserat杆的非线性动力学。该混合公式具有独立的位移、速度和应力变量，并且不含目标函数和锁死现象。有限旋转通过一种称为“导向器公式”（director formulation）来表示，这种公式可以避免奇异性，并生成一个恒定的质量矩阵。由此形成的是一个无限维的非线性Port-Hamiltonian（PH）系统，该系统由偏微分-代数方程控制，并具有二次能量泛函。通过使用应力-应变关系的时间微分形式，可以施加运动约束（如不可伸长性或剪切刚性）。采用保持结构特征的有限元离散化方法，可以将系统转换为有限维的PH结构，从而便于设计能量-动量一致性的积分方案。耗散材料行为（通过广义Maxwell模型）和非标准驱动方式（如气动腔或肌腱）可以自然地融入该框架。通过选定的数值示例可以看出，该框架为涉及有限旋转的计算力学中的能量-动量一致性公式提供了一种新的方法。

**1. 引言**
由刚体和柔性组件组成的柔性多体系统（[3]）在软体机器人[4]或缆索动力学[5]等领域有着广泛的应用。这些系统本质上是复杂的，受到几何非线性和复杂材料行为的影响，因此需要准确且稳健的仿真方法[6]。
为此，本文提出了一种用于Cosserat杆动力学的新型建模框架[7]。在此框架中，控制方程被重新表述为一个Port-Hamiltonian（PH）系统[8]。这种方法天然地采用了独立的位移、速度和应力变量，有助于在空间和时间上进行结构保持的离散化。

**Cosserat杆**
Cosserat杆也被称为Simo-Reissner梁[9]、[10]或几何精确梁[11]、[12]。这些模型已经变得非常具有影响力[11]、[12]、[13]、[14]、[15]、[16]、[17]、[18]、[19]、[20]、[21]、[22]，因为它们避免了线性或更高阶的近似。最近，这些模型在软体机器人[23]、[24]、[25]、[26]中也越来越多地被使用。
该模型的复杂性主要源于其几何非线性，这种非线性来自于有限旋转，其参数化方法有多种[3]。参数化旋转自由度的候选方式包括但不限于欧拉角、（投影的）单位四元数或相对旋转向量。本文将考虑一种无奇异性的参数化方法，这种方法完全避免了使用旋转坐标，并生成一个恒定且对角线的质量矩阵——即导向器框架。由于其优越的特性，Cosserat杆的导向器公式在有限元方法的设计中已成为广泛应用的选择。相应的有限元实现方法要么依赖于代数节点约束的施加[11]、[15]、[16]、[27]、[28]，要么引入节点旋转自由度[12]、[29]、[30]、[31]、[32]。关于旋转参数选择的更多细节，请参见相关文献[33]、[34]、[35]、[36]。

**Port-Hamiltonian系统**
Port-Hamiltonian（PH）框架为互联动力系统的建模、仿真和控制提供了一种统一的方法[8]，非常适合柔性多体系统。PH描述考虑了能量耗散和与环境的相互作用，因此可以被视为经典哈密顿动力学的扩展。PH模型最初是为有限维系统提出的[37]，后来被扩展到无限维系统[38]，并在多个领域进行了广泛研究[39]。在结构动力学的背景下，已经为几何精确的弦[40]、[41]、[42]、板[43]、[44]以及受限于小或中等变形的梁[45]、[46]、[47]、[48]、[49]开发了PH公式。非线性基于应变的梁模型也在相关研究中得到探讨[50]、[51]、[52]。这些公式完全基于固定于物体的坐标系，与Hodges对Cosserat杆的描述密切相关[53]。PH描述与结构力学中常用的描述方法之间的等价性已在文献[54]中得到证明，相关的PH表示方法在文献[55]中有详细说明。

**保持结构特征的离散化**
PH公式天然具有能量一致性，因此非常适合设计保持结构特征的离散化方法[6]，这些方法能确保在离散层面上遵循基本的功率平衡。这可以避免数值耗散和能量不一致性，特别是在长期仿真和控制导向问题中非常重要[56]、[57]。同时保持能量和动量的数值方案通常被称为EM方案。尽管[58]是第一个用于空间Cosserat杆的EM方法，但它违反了客观性这一重要性质[13]。第一个保留客观性的空间Cosserat杆EM方法归功于[12]（另见[11]、[28]），该方法利用了上述的导向器框架。

在无限维PH系统的背景下，空间离散化同样通过有限元（FE）方法实现[59]，从而得到有限维的PH系统（参见[40]、[41]、[45]、[47]、[60]、[61]在结构动力学中的应用）。需要注意的是，PH方法中典型使用的能量变量本质上会导致混合FE公式，从而也包括独立的应力和/或应变场。同样，时间离散化对于实现具有能量一致性的精确仿真至关重要。用于PH系统的能量一致积分方案包括离散梯度公式[41]、[62]、[63]、[64]和基于投影的方法[65]。

**1.1. 研究空白**
据作者所知，针对Cosserat杆动力学的混合公式和FE方法仍然相对较少。尽管对该梁公式和PH框架进行了大量研究，但它们的集成工作却受到有限关注。现有的基于应变的PH公式[50]、[52]、[55]需要重建绝对位置和方向，这在PH社区之外的研究中也是如此[66]、[67]。此外，[50]中的李群设置引入了额外的几何复杂性，而[50]、[52]、[55]中的离散化方法也没有使用有限元。相比之下，计算力学社区提出的公式主要限于静态场景[21]、[68]、[69]、[70]、[71]、[72]。
在细长机械结构的仿真中，另一个挑战是横向剪切锁定现象。这种数值硬化效应在离散化过程中会出现，导致不合理的变形。虽然选择性减少积分[12]、[27]是一种流行的解决方案，但其使用混合方法的处理较少被讨论（见例外情况[21]、[61]、[68]、[71]、[73]在PH背景下）。
线性本构关系可以自然地通过基于Hellinger-Reissner变分原理的顺应性方程纳入PH公式[21]、[68]、[70]、[71]。所谓的速度-应力公式在PH社区中很流行[40]、[47]、[74]，它们考虑了时间微分形式的顺应性方程。然而，这些工作并未明确指出与Hellinger-Reissner原理的联系。对于Cosserat杆而言，这样的公式还有助于轻松防止某些变形模式，从而无需进一步修改公式即可模拟剪切刚性或不可伸长的梁[21]。然而，将这些公式应用于Cosserat杆的动力学，尤其是与PH公式结合使用，仍然是一个未解决的问题。

最后，用广义Maxwell模型[75]、[76]、[77]描述粘弹性以及使用气动腔和肌腱[23]、[24]、[24]、[25]、[26]进行驱动在连续体机器人技术中非常重要。将其纳入PH框架仍然是一个未解决的问题。

**1.2. 贡献**
先前的研究[42]、[61]表明，混合PH公式可以推进非线性结构动力学的结构保持方法的发展。在这些研究中，特定的PH框架[62]、[78]、[79]、[80]也被证明在无限维环境中适用。在本研究中，我们采用了这种公式来处理空间Cosserat杆，解决了与有限旋转相关的特定挑战，并提出了以下改进：
(C1) 我们提出了一种用于空间Cosserat杆的混合无限维PH公式，该公式自然地包含了独立的位移、速度和应力变量，同时确保了能量-动量的一致性。所采用的有限旋转导向器参数化方法生成了恒定且对角线的质量矩阵。
(C2) 我们展示了使用混合有限元素和二阶时间步进方法进行的结构保持空间和时间离散化，这得益于混合PH公式。我们的方法确保了功率和角动量平衡的精确离散表示。
(C3) 我们证明了由PH框架引导的FE方法可以减轻剪切锁定现象，并直接包括剪切刚性和不可伸长梁的公式作为特殊情况。这可以在时间微分形式的顺应性方程背景下轻松实现。
(C4) 我们将给定的PH模型扩展到考虑粘弹性材料行为（通过广义Maxwell模型）和气动腔及肌腱驱动。这些扩展与所提出的PH公式非常契合。

**1.3. 符号说明**
我们遵循以下符号规则：标量f、g、n、m、s∈R用斜体字表示；向量及其d元组等效于一列矩阵，例如a、b、c∈R↑d×1，用粗体字表示。所有表示物理量的三维向量都直接用坐标表示，即相对于固定的惯性正交基的R3中的三元组。因此，该惯性基的基向量具有规范坐标表示：e1=(1,0,0)、e2=(0,1,0)和e3=(0,0,1)。使用爱因斯坦求和约定，即对于拉丁字母指标i、j、k∈{1,2,3}和希腊字母指标α、β∈{1,2}，三维向量可以写作a=∑i=13aiei=aiei=(a1,a2,a3)。为了方便表示，使用简短括号符号?[a?,b?]??(a,b)。方块符号↑表示占位符。标量积用a?b=aibi表示，R3中的叉积表示为?c=a×b=?ijkajbkei=a?b，其中我们使用了Levi-Civita符号??ijk，表示R3与斜对称矩阵空间之间的映射so(3)?R3×3。Kronecker delta δij在i=j时为1，在i≠j时为0。矩阵A∈Rn×m用粗体表示，涉及导数的微分矩阵运算符用花体大写字母A表示。矩阵转置用A?表示，运算符A的形式伴随表示为A?。矩阵-向量乘积和运算符的作用分别表示为Ab和Ab。单位矩阵用I表示，全零矩阵用0表示，其维度从上下文中明确。函数f的偏导数表示为?f/?x。我们进一步引入了?和?2等符号。在一维域Ω上的内积表示为?(a,b)Ω?∫Ωa?bds或?(f,g)Ω?∫Ωfgds，以考虑正确的维度。

**1.4. 结构安排**
本文的其余部分结构如下：第2节介绍了与空间Cosserat杆动力学相关的问题。第3节将模型重新表述为一阶系统，并展示了其PH结构，从而自然引出了功率平衡方程。借助PH公式，这种功率平衡方程可以通过结构保持的离散化方法转换为离散形式（见第4节）。第5节讨论了将公式扩展到粘弹性和两种驱动机制，这些在连续体机器人应用中非常重要。最后，第6节提供了数值示例，验证了所提方法的性能。这些示例在第7节的结论中进行了总结。不直接影响主要论点的推导和细节包含在附录A（导向器速度和虚拟导向器位移）、附录B（虚拟功原理）、附录C（形式斜伴随结构算子证明）、附录D（总角动量平衡）和附录E（空间离散化后的系统矩阵）中。

**2. Cosserat杆模型**
**2.1. 运动学**
以下概念基于文献[12]、[27]中的模型。考虑一个一维材料梁域Ω=[0,L]?R，其中L∈R，梁上每个材料点的运动通过φ(s,t)∈R3描述。这里的材料坐标s∈Ω是无应力初始配置Ω0={φ(s,0)=sLe3|s∈Ω}的弧长参数化，t∈[0,tend]?R表示时间。这种几何精确的运动学描述使得每个横截面在运动过程中保持平面且刚性，其方向通过旋转矩阵(1)R(s,t)=[d1(s,t)d2(s,t)d3(s,t)]=dk(s,t)ek?描述，其中ek∈SO(3)是特殊正交群。R的列形成了局部正交导向器di(s,t)∈R3（i=1,2,3）。参见图1了解运动学的示意图。假设导向器d3(s,t)垂直于每个横截面，而d1(s,t)和d2(s,t)垂直于d3(s,t)，即d3(s,t)=d1(s,t)×d2(s,t)。为了简化问题，我们进一步假设d1(s,t)和d2(s,t)在任何时候都与截面的主轴惯性重合。在下文中，省略了对空间和时间的显式依赖性，只有在必要时才会提及。滥用符号，我们也将R=R(d)解释为垂直结合的指向量d=(d1,d2,d3)∈R9的函数。此外，我们称(φ,d)为位移量，因为它们唯一地表征了梁配置上每个材料点的位移。由于R(d)∈SO(3)，因此有det(R(d))=+1且R(d)?R(d)=R(d)R(d)?=I。因此，有6个正交条件：(2)0=g(d)∈R6?12di?dj?δij=0，对于(i,j)∈{(1,1),(2,2),(3,3),(1,2),(2,3),(3,1)}。我们使用客观的材料应变度量Γ=R?φ,s来描述剪切变形和膨胀，以及K?=R?R,s来描述曲率和扭矩，它们可以等价地表示为(3)?Γ=(dk?φ,s)ek=d1?φ,sd2?φ,sd3?φ,s,K=12?ijk(dk?dj,s)ei=12d3?d2,s?d2?d3,sd1?d3,s?d3?d1,sd2?d1,s?d1?d2,s，其中我们转换到了K?的轴向量。由于Γ(d,φ,s)和K(d,d,s)都是双线性的，可以使用两种不同的表示方法来紧凑地重写它们，即(4)Γ=R(d)?φ,s=JN(φ,s)?d,K=L(d,s)d=?L(d)d,s。旋转矩阵(1)以及新引入的矩阵(5)都是它们参数的线性函数。注意，由于(2)，L(d)d=0。

为了简化分析，我们进一步假设d1(s,t)和d2(s,t)在任何时候都与截面的主轴惯性重合。在接下来的讨论中，省略了对空间和时间的显式依赖性，只有在必要时才会提到。滥用符号，我们也将R(R(d))解释为指向量d=(d1,d2,d3)∈R^9的垂直组合的函数。此外，我们称(φ,d)为位移量，因为它们唯一地表征了梁上每个材料点的位移。由于R(d)属于SO(3)，因此有det(R(d))=+1且R(d)^T R(d)=R(d)R(d)^T = I。因此，存在6个正交条件：(2)0 = g(d)∈R^6 ? 12di^T dj ? δij = 0，对于(i,j) ∈ {(1,1), (2,2), (3,3), (1,2), (2,3), (3,1)}。

我们使用客观的材料应变度量Γ=R^T φ, s来描述剪切变形和膨胀，以及K?=R^T R, s来描述曲率和扭矩，这些可以等价地表示为：(3)?Γ=(dk^T φ, s)ek=d1^T φ, s d2^T φ, s d3^T φ, s, K=1/2 ?ijk(dk^T dj, s)ei=1/2 d3^T d2, s?d2^T d3, s d1^T d3, s?d3^T d1, s d2^T d1, s?d1^T d2, s。这里我们转换到了K?的轴向量。由于Γ(d,φ, s)和K(d, d, s)都是双线性的，可以使用两种不同的表示方法来紧凑地重写它们，即(4)Γ=R(d)^T φ, s=JN(φ, s)^T d, K=L(d, s)d=?L(d)d, s。

图1. Cosserat杆配置的示意图。

注释2.1
我们假设初始配置是直的且无应力的，这作为参考配置。因此，截面法线d3(s, t=0)与中心线切向量φ, s(s, t=0)对齐，这意味着均匀的参考应变约为Γ0(s) ≈ Γ(s, t=0) = (0,0,1)且Γ0(s) ≈ K(s, t=0) = 0。对于预弯曲的梁，参考应变会有所不同，但这不影响本工作中开发的理论框架的普遍性。此外，Γ描述了中心线（即切向量）的变化率，而K描述了沿梁移动时的方向变化率。K是从斜对称矩阵K?导出的，这可以追溯到对正交条件(2)关于s求导。

2.2. 构成建模
从基于标量值应变能函数W=W?(Γ, K)的超弹性材料模型的一般情况出发，相应的功 conjugated 应变量可以推导为：(6)N=?Γ W?, M=?K W?。这里，用大写字母表示的N=Niei是材料法向力和剪力，M=Miei是材料弯矩和扭矩矩。空间应力量n和m可以通过以下方式与它们的材料对应量联系起来：(7)n=R(d)N=Nidi, m=R(d)M=Midi。

在接下来的分析中，我们假设线性本构律，使得W具有二次型。这是文献中的常见假设，因为通常预计局部只有小的应变。因此我们写出：(8)W?(Γ, K)=1/2(Γ?Γ0)^T CΓ(Γ?Γ0)+1/2(K?K0)^T CK(K?K0)，其中常数弹性矩阵CΓ=diag(ks1, ks2, ke)∈R^3×3和CK=diag(kb1, kb2, kt)∈R^3×3，以及Γ0和K0根据注释2.1。

二次超弹性应变能函数也可以用应力表示为互补能量函数W?(N, M)，即：(9)W?(N, M)=1/2N^T CNN+1/2M^T CMM，其中 Compliance 矩阵CN=diag(ks1?1, ks2?1, ke?1)∈R^3×3, CM=diag(kb1?1, kb2?1, kt?1)∈R^3×3。注意，如果某些刚度参数趋于无穷大（例如，ks1, ks2→∞），这些矩阵可能是奇异的，例如在众所周知的Kirchhoff梁模型中，CN=diag(0,0, ke?1)。类似于(6)，应变量可以通过微分得到：(10)Γ?Γ0=?N W?=CNN, K?K0=?M W?=CMM。

在超弹性应力-应变关系无法求解的情况下，可能需要采用Hu–Washizu类型的方法[41], [61]。

2.3. 运动方程
控制Cosserat杆动力学的偏微分方程可以从虚功原理推导出来，见附录B，即：(11)ρAφ?=n, s+n?, Iρω?+ω×Iρω=m, s，分别表示空间形式的线性和角动量平衡[7]，[9]。这里，我们引入了每单位长度均匀分布的质量密度ρ≥0和截面面积A>0，ω(s, t)∈R^3是空间角速度，Iρ∈R^3×3是惯量张量的空间表示，n?(s, t), m?(s, t)∈R^3是每单位长度的外部分布力和力矩。

与(11)不同，我们目标是使用上面介绍的指向量来表述。考虑到附录B中的推导，可以得到Cosserat杆动力学的等效表示形式：(12a)ρAφ?=n, s+n?, (12b)Mρd?=?JNN?JMM?Gd^T λ+Tm?, (12c)g(d)=0，其中JN已在(5)中引入。此外，我们引入了一个常数和对角质量矩阵Mρ=diag(Mρ11I, Mρ22I, 0)∈R^9×9，以及算子矩阵JM，其及其形式伴随矩阵JK=JM?的定义见(13)。在(12b)中，我们还引入了T∈R^9×3，定义为(14)T(d)=?1/2d?1d?2d?3，其中T(d)^T=1/2[d?1d?2d?3]，它连接角速度和指向量速度，见附录A中的(85)。最后，约束梯度矩阵(15)和包含所有拉格朗日乘数λ∈R^6的向量在(12b)中编码了约束力。这些力确保满足正交条件(12c)，这些条件等同于(2)。注意，N和M是通过(4), (6)根据φ和d计算得出的。相关的初始边界值问题通过适当的初始和边界条件来完成。

纯粹基于位移的公式(12)在能量-动量一致离散化的背景下容易锁定且处理起来更繁琐。因此，我们在下一节中采用混合的Port-Hamiltonian（PH）描述。

3. 混合公式作为Port-Hamiltonian系统
3.1. 一阶系统的局部形式
我们的目标是将当前问题重写为一组带有限制条件的一阶偏微分方程。为此，我们引入了独立的应力和速度变量，并通过时间微分形式的合规性方程将它们联系起来，同时也考虑了速度级别的正交性约束。让我们引入独立的中心线速度vφ=φ?和与指向量相关的速度vd=d?，以及位移和速度向量q, v∈R^12，使得(16)q=φd, v=vφvd，且q?=v。利用新引入的速度，我们可以将平衡方程(12a), (12b)重写为(17)ρAv?φ=?s(RN)+n?, Mρv?d=?JNN?JMM?Gd^T λ+Tm?。接下来，对合规性方程(10)和(4)进行时间微分，得到应力的演化方程(18)CNN?=R(d)^T vφ, s+JN(φ, s)^T vd, CM?=JK(d, d, s)vd。注意我们在(18)中使用了算子矩阵JK(d, d, s)，它在(13)中定义。最后，使用指向量作为独立量需要明确执行6个正交条件(12c)，我们也将它们以时间微分形式考虑，即(19)G(q)v=0, Gd(d)vφvd=0，其中约束矩阵Gd(d)在(15)中定义。最终，结合适当的初始和边界条件，方程(16), (17), (18), (19)构成了Cosserat杆动力学的混合PH公式。为了简洁起见，我们在下一节中引入(20)来紧凑地重写(17), (18)。

3.2. Port-Hamiltonian表示
当前控制的偏微分-代数方程的PH表示包括方程(16), (17), (18), (19)。使用定义(20)，我们可以将它们重写为(21)，其中所谓的描述符矩阵E=diag(I, Mo, C, 0)是常数，由于6个代数约束和奇异的质量矩阵，其秩至少差9。如果假设某些刚度趋于无穷大（例如，描述一个剪切刚性梁），则秩差会相应增加，因为合规性矩阵C变得奇异。简而言之，可以用[78], [79], [80]2将(21)重写为(22)Ex?=J(x)z(x)+B(x)uΩ, yΩ=B(x)^T z(x)，其中状态向量x∈R^3^6包括主要未知数(23)x=(q, v, σ, λ)。(22)中引入的其他新量可以通过与(21)比较来识别。此外，J(x)是一个微分算子矩阵，称为结构算子矩阵，它涉及偏导数。进一步，引入了代价函数z，其形式为z(x)=Qx，其中Q=diag(0, I, I, I)。在(22)中，我们附加了分布式输出关系，得到yΩ=(vφ, ω)，即中心线速度和空间角速度ω。根据定义，这种分布式输出与分布式输入力和扭矩uΩ=(n?, m?)共轭。

此时，需要证明上述方程具有PH结构。这要求(i)代价函数z正确地与系统的底层哈密顿量关联，以及(ii)结构算子J是形式上的斜对称的，即J(x)=?J(x)^?。这两个条件确保PH系统具有底层的Stokes–Dirac结构，因此满足功率平衡方程。条件(i)要求(24)E^T z(x)=δxH，其中δxH是适当哈密顿量H的变分导数，见[78], [79], [80]。根据E和z(x)的定义，对于由(25)给出的哈密顿量H(x)=∫Ωh(x)ds=1/2∫Ωx^T Q^T Ex ds=1/2∫Ωv^T Mov+σ^T Cσ ds，这是系统的总机械能，即动能和由于变形储存在梁中的互补能之和。因此，δxH=?xh=Q^T Ex，且(24)成立，因为Q^T E=E^T Q。注意哈密顿量不依赖于q和λ，因为它们不影响系统的能量。这些属性通过矩阵Q和E中的相应零块反映出来。

关于条件(ii)，J在(21)中的形式斜对称性可以看作是矩阵对微分算子斜对称性的延伸。简而言之，这意味着计算z, JzΩ只能在分部积分后得到边界项。对于不涉及空间导数的元素（如I和G），这通过检查斜对称性可以直接得出。然而，对于微分算子Jσ和Jv，需要更多的注意。最终，根据附录C中的方法得到结果(26)?z, JzΩ=vφ^T n+ω^T m?Ω?y?Ω?u?Ω。

3.3. 弱形式
问题(29a)ρAδvφ, φ??vφΩ=0, (29b)Mρδvd, d??vdΩ=0, (29c)δφ, ρAv?φ?n?Ω+δφ, s, RNΩ?δφ, nN?Ω=0, (29d)δd, Mρv?d+JNN+Gd^T λ?Tm?Ω+JKδd, MΩ?δd, TmN?Ω=0, (29e)δN, CNN??R^T vφ, s?JN^T vdΩ=0, (29f)δM, CMM??JKvdΩ=0, (29g)δλ, GdvdΩ=0，对所有可接受的和足够光滑的测试函数都成立。它可以从微分形式(21)直接得出，并作为以下FE离散化的基础。将(21)两边预先乘以适当的测试函数(Moδv, δq, δσ, δλ)，在空间域Ω上积分并应用分部积分，得到上述弱关系。这种弱形式也可以与Hellinger–Reissner原理联系起来，见附录B.2。

注释3.2
如果考虑在可接受测试函数空间中施加的Dirichlet边界条件，弱形式(29)是有效的。然而，如果以弱形式施加Dirichlet边界条件，则必须稍作修改，例如，通过引入额外的拉格朗日乘数并考虑规定的边界值。通常，对于PH系统施加非均匀混合边界条件需要额外的注意[40], [81], [82], [83]。

注释3.3
值得注意的是，所提出的公式还满足总角动量的基本平衡，这可以通过弱形式(29)证明，见附录D。

3.4. 客观性
让我们讨论一下在叠加刚体运动的情况下，或者等效地，在观察者变化的情况下，所提出的公式（具有弱形式(29)的行为。为此，考虑c(t)∈R^3和Φ(t)∈SO(3)，它们规定了梁的叠加平移和旋转。因此，运动量变为(30)φ#=c+Φφ, di#=Φdi, vφ#=c?+Φvφ+Φ?φ, vd, i#=Φvd, i+Φ?d，对于i∈{1,2,3}，这里表示考虑了叠加刚体运动的转换量。外部分布式载荷经过变换后变为（31）n?#=Λn?，m?#=Λm?。材料应变（3）在这种变换下保持不变，参见[27，第1781页]和[12，第1689页]。相应地，可以声明（32）Γ(d#,?sφ#)=Γ(d,?sφ)，K(d#,?sd#)=K(d,?sd）。利用Λ??Λ+Λ?Λ?=0，还可以证明相应的应变率Γ?和K?也不变，这样在时间微分形式（18）中的顺应性方程就可以被变换后的量完全满足。这确保了应力测量及其时间导数的客观性。此外，速度水平上的正交性约束（19）不会受到这种叠加的刚体运动的影响。最后，将运动学关系（16）和平衡方程（17）分别与Λ和diag(Λ,Λ,Λ)相乘，会得到用变换后的量表示的相同关系。这个证明也可以通过使用适当变换的测试函数的弱形式（29）来完成。因此，当前对于几何精确梁的混合公式是客观的。值得注意的是，这一性质不会被后续章节中提出的离散化方法所改变，参见备注4.1.4。

保持结构的离散化
4.1. 在空间中的混合有限元离散化
我们采用与[61]中的2D情况相同的混合有限元（FE）方法来保持空间的离散化。Cosserat杆动力学的PH系统（21）可以通过其弱形式（29）方便地进行空间离散化。
离散化假设
我们根据通用标准程序[84]将域离散为ne个有限元，即Ω=∪e=1neΩe。根据假设函数的选择，这对应于FE网格上的nn个节点。位移和速度量（q,v）用p阶的C0-连续拉格朗日多项式假设函数来近似，而应力类型量σ则使用p?1阶的不连续拉格朗日多项式假设函数来近似。相应的空间表示为（33）Vh={qh∈H1(Ω,R12)|qh|Ωe∈Pp(Ωe,R12)?Ωe∈Ω,qh|?ΩD=qD}，V0h={qh∈H1(Ω,R12)|qh|Ωe∈Pp(Ωe,R12)?Ωe∈Ω,qh|?ΩD=0}，Sh={σh∈L2(Ω,R6)|σh|Ωe∈Pp?1(Ωe,R6)?Ωe∈Ω}，其中p∈{1,2}，表示近似。此外，H1表示平方可积的C0-连续函数的Sobolev空间，L2表示平方可积函数的空间，Pp表示p阶的多项式函数空间。通过这样的选择，我们有nn=nep+1。
拉格朗日乘数是按照配置法计算的，即在节点处取离散值，在其他地方为零，它们的近似空间表示为Lh。这可以通过delta函数[85，第1.11章]正式实现。相应地，方向量的正交性约束（19）被放宽到FE网格的节点点上。这便于使用拉格朗日多项式来处理（q,v），并且在文献[12]、[27]、[28]中是一种常见的方法。
关于节点之间的插值，“方向量框架的正交性的缺失可以被视为离散化误差，如果增加元素的数量，这种误差会减小”[15]。正如后面所使用的，使用二次假设函数进行插值得到的插值与考虑正确流形的SE(3)有限元相比几乎没有区别，见[32]中的图2。
然后，使用Bubnov–Galerkin方法的离散弱形式可以表述如下。寻找（qh,vh,σh,λh）∈Vh×Vh×Sh×Lh，使得（29）对所有（δqh,δvh,δσh,δλh）∈V0h×V0h×Sh×Lh都成立。为此，我们用假设矩阵Φ,Ψ,Ξ来表示近似量（34）qh(s,t)=Φ(s)q?(t)∈Vh，vh(s,t)=Φ(s)v?(t)∈Vh，σh(s,t)=Ψ(s)σ?(t)∈Sh，λh(s,t)=Ξ(s)λ?(t)∈Lh，其中q?,v?∈R12(nep+1)，σ?∈R6nep，λ?∈R6(nep+1)包含了所有有限元和相应节点的未知数。此外，对于相应测试函数的近似，考虑与假设相同的基。在这个上下文中，上述假设矩阵包含了多项式假设函数，例如，qh(s,t)=Φ(s)q?(t)=∑i=1nnNi(s)q?i(t)。更多细节见[61]中的附录A.2（针对平面情况）。我们进一步遵循常见的FE程序，例如[84]，其中还包括了考虑相应自由度的全局系统矩阵和向量的组装。此外，我们都在元素级别进行所有计算，并使用等参概念将所有计算转换为标准参考元素，其中ξ∈Ω?=[?1,1]。
有限维的PH系统
通过用假设空间（34）近似弱形式（29）来获得离散系统动力学。考虑到测试函数节点值的任意性，我们得到了PH微分-代数方程（PH-DAE）（35）。
离散系统矩阵可以在附录E中找到。类似于无限维版本（22），我们可以在PH框架[78]、[79]、[80]中重新构造半离散方程组（36）E?x??=J?(x?)z?(x?)+B?(x?)u?，y?=B?(x?)?z?(x?)，其中我们还加入了离散的输出方程y?=(y?Ω,y??)，包含了分布式的节点和边界项。这些量在物理上有速度的含义，它们是与节点分布式和边界输入力和力矩的功率共轭的。更多细节见附录E。这里的有限维版本的状态（23）包含了所有有限元素的自由度，表示为（37）x?=(q?,v?,σ?,λ?)。
将假设（34）代入（25）可以得到近似的哈密顿量（38）H(xh)=H?(x?)=1/2v??M?ov?+1/2σ??C?σ?=1/2x??Q??E?x?，其中Q?=diag(0,I,I,I)，更多的系统矩阵可以在附录E中找到。
PH结构和功率平衡
我们可以再次检查手头方程的PH结构。首先，根据无限维设置中的条件（i），可以声明（39）E??z?(x?)=?H?(x?)对于半离散系统是成立的，因为离散哈密顿量（38）和PH-DAE（35）中E?和z?(x?)=Q?x?的定义。其次，结构算子J的形式上斜对称的性质得到了保持，即结构矩阵J?(x?)=?J?(x?)?在（35）中是斜对称的。注意J?的状态依赖性突出了当前问题的几何非线性。这对应于无限维设置中的条件（ii）。
如前所述，PH公式的一个核心性质是它满足形式化的功率平衡方程（28）。由于（36）是一个离散的PH-DAE系统，它满足半离散功率平衡方程（40）H??=?H?(x?)?x??=z?(x?)?E?x??=z?(x?)?J?(x?)z?(x?)+B?(x?)u?=y??u?，其中我们使用了（36）、（39）和J?(x?)的斜对称性。这个方程类似于（28）的离散版本，表明在离散化过程中保持了被动性和无损耗性。此外，总应变等于每个分支中的应变，即 (56)????(φ,d)=?∞=?el(i)+α(i)，对于每个 i，其中 α(i) 是相应分支中阻尼元件的粘性应变，而 ??∞ 是纯弹性分支中的应变。一个常见的假设是粘性应力 (57)σvisc(i)=V(i)α?(i) 通过一个对称的、正半定的粘性矩阵 V(i) 与粘性应变率线性相关。下载：下载高分辨率图像（80KB）下载：下载全尺寸图像图2. 广义麦克斯韦方法的流变学模型示意图。为了考虑 PH 模型中的这种线性粘弹性行为，我们首先为每个分支引入额外的状态变量 σel(i)=σvisc(i)。接下来，我们根据 (55) 将额外的应力贡献包括在动量平衡 (21) 中。最后，我们根据 (56) 中的两个等式的时间导数，考虑演化方程 (58)?C∞σ?∞=??(φ,d)，?C(i)σ?el(i)=??(φ,d)?V(i)?1σel(i)，i=1,…,m，对于长期弹性应力 σ∞ 和每个单独的弹性应力贡献 σel(i)。我们还考虑了本构关系 (57) 以及 ?C∞σ∞=?∞ 和 ?C(i)σel(i)=?el(i)。注意，为了与纯弹性梁模型一致，其弹性矩阵为 CΓ 和 CK，请参见方程 (8)，所有分支的弹性矩阵必须加起来等于这些矩阵。

注 5.1 一种常见的方法 [75]，[88] 是假设粘性矩阵 V(i) 与弹性矩阵具有相同的结构。为此，可以引入动态粘度 ηG，ηE≥0 分别用于剪切和扭转以及拉伸和弯曲，使得 (59)τG=ηGG，τE=ηEE 是相应的特征松弛时间。这里，E，G≥0 分别表示杨氏模量和剪切模量。这最终导致粘性矩阵 (60)，其中 CΓ(i) 和 CK(i) 与第 i 个粘性分支的顺应矩阵 C(i) 一致。

注 5.2 如果考虑 m=1 的耗散麦克斯韦分支，Kelvin–Voigt 散热模型（例如在 [19]，[88] 中使用）包含在当前框架中。在这个分支中，弹性元件的顺应性必须为零，以便它表现得像刚性的一样。

耗散 PH 模型。耗散系统可以再次从 PH 视角来理解。为此，我们将所有麦克斯韦分支的应力连接起来，得到新的 PH 系统的状态和余态 (61)x=(φ,v,σ∞,σvis,λ)，z(x)=(0,v,σ∞,σvis,λ)。以及哈密顿量 (62)H(x)=∫Ω12v?Mov+12σ∞?C∞σ∞+12σvis?Cvisσvisds，包括之前的动能以及与广义麦克斯韦模型中的弹性元件相关的所有补充应变能项，其中 Cvis=diag(C(1)，…,C(m))。同样，我们定义 V?1=diag((V(1))?1,…,(V(m))?1)，Jvvis=Jv…Jv 和 Jσvis=Jσ?…Jσ??，每个量都有 m 个副本。因此，我们将动力学表示为 (63)。简而言之，这再次符合 PH 框架 [78]，[79]，[80] (64)Ex?=(J(x)?R)z(x)+B(x)uΩ,yΩ=B(x)?z，其中 R=diag(0,0,0,V?1,0) 是一个常数、对称的、正半定的矩阵。

6 正半定性质可以通过块对角结构以及 (V(i))?1 的对称性和正半定性来验证，见注 5.1。这个 PH 系统引起的功率平衡方程扩展了纯弹性模型 (28)，即 (65)H?=y?Ω?u?Ω+∫ΩyΩ?uΩds?D。在这里，我们引入了由于粘性变形而产生的耗散功率（也称为耗散函数 [87]）(66)D=∫Ωz?Rzds=∫Ωσvis?V?1σvisds=∫Ω∑i=1m(σel(i))?(V(i))?1σel(i)ds≥0。因此，功率平衡方程 (65) 立即导致了耗散不等式 (67)H?≤y?Ω?u?Ω+∫ΩyΩ?uΩds，证明了耗散无限维 PH 系统的被动性。

注 5.3 如 [42] 所示，对于几何上精确的弦，我们也可以通过功率保持互连来获得模型 (63)。在这种情况下，所有粘性麦克斯韦分支被解释为纯耗散的 PH 模型，而弹性分支则被视为其自身的模型，类似于 (21)。整个模型是通过被动子系统的回转器互连获得的。

5.2. 扩展到气动和肌腱驱动在软体机器人操纵器的背景下，两种先进的驱动方法是由充当所谓肌腱的弹性缆绳 [23]，[24] 和气动腔室 [24]，[25] 提供的。这两种方法都可以以统一的方式无缝集成到当前的梁模型中。我们遵循 [24]，[25] 中的方法，也参见 [26]，并包括作为分布式驱动项的额外应力贡献。

因此，我们通过引入额外项来扩展运动方程 (21)，即 (68)N→N+R(d)?nu，M→M+R(d)?mu，其中 nu，mu∈R3 是来自输入驱动的额外内部力和力矩。对于肌腱和压力腔室驱动，它们分别具有形式 (69)nu=∑k=1Nτktk，mu=∑k=1Nρk×τktk，其中 k 是驱动元件的索引，即肌腱或压力腔室，以及 (70)τk=pk(t)Ak≤0 对于压力驱动，Tk(t)≥0 对于肌腱驱动，Ak 是加压腔室的截面面积，而 Tk 是肌腱中的拉力。pk(t) 和 Tk(t) 都被视为已知输入量，它们将与分布式驱动项相关联。此外，在 (69) 中，tk 是相应的单位方向向量，ρk(s,t)=ρkαdα(s,t)，对于 α=1,2，描述了从梁的中心线到第 k 个驱动元件的质心线的向量。这里，ρkα 是固定的材料坐标。相应地，rk=φ+ρk。两种驱动方法之间的唯一区别是肌腱在肌腱的切线方向上产生牵引力，而压力腔室在截面法线方向上产生牵引力。这是通过单位方向向量 tk 实现的，即 (71)tk=d3 对于压力驱动，rk,s||rk,s|| 对于肌腱驱动。

对于肌腱驱动的情况，见图 3 的示意图。在强形式 (21) 中，这些项导致额外的分布式输入形式 Bu(x)τ，其中 τ=(τ1,…,τN)∈RN 是额外的输入量。为了考虑 (17)，(68)，(69)，我们引入了分布式端口矩阵 (72)Bu(x)=(tk,s)k=1N?(JNR?tk+JM(R?ρ?ktk))k=1N∈R6×N，其中矩阵的每一列对应于一个驱动元件的影响。注意，端口矩阵的状态依赖性不仅来源于操作符，还来源于 (71) 中引入的单位方向向量 tk。同样，这也导致 PH 动力学的微分形式 (22) 中的额外输入项，最终扩展到 (73)。

类似的，这些额外贡献也通过额外的项 δφ,s，nuΩ 和 JKδd,R?muΩ 在 (29c)，(29d) 中得到考虑。

下载：下载高分辨率图像（144KB）下载：下载全尺寸图像图 3. Cosserat 棒的驱动示意图，其中显示了第 k 个驱动元件和相关的运动学量。

6. 数值示例让我们将提出的框架应用于代表性的数值示例。使用牛顿方法在每个时间步长解决非线性离散时间方程 (42)，对残差的范数有一个绝对容差 ??。空间积分通过三点高斯求积法近似，对于多项式假设函数我们选择 p=2。位移和速度的狄利克雷边界条件直接纳入假设函数中，见注 3.2，并且没有考虑重力效应。计算是使用有限元代码 moofeKIT 执行的，生成的数据发表在 [1] 中并存档在 [2] 下。第一个示例用于评估所用离散化方法的结构保持能力。此外，还显示了与分析解的良好一致性。在第二个基准测试中，我们的结果与文献中报告的数据非常吻合，并根据第 4.3 节验证了运动学一致性属性。这个示例还处理了第 5.1 节中的粘弹性扩展，并强调了当前框架对于不可伸展的 Kirchhoff 梁的适用性。第三个示例在准静态场景中研究了公式化，表明当前公式化与分析解相匹配并且避免了锁定。最后，我们使用我们的公式来模拟由气动腔室驱动的软体机器人臂的运动，如第 5.2 节所述。

表 1. 飞行意大利面：物理和数值参数。htendne??LρAMρ11=Mρ22kb1=kb2=ktks1=ks2=ke空单元0.1151010?1110110103104下载：下载高分辨率图像（481KB）下载：下载全尺寸图像图 4. 飞行意大利面：以方位角和仰角 (55,15) 观察的皮划艇划行运动的快照，以及时间 tn 的颜色图。

6.1. 飞行意大利面这个第一个基准问题最初在 [58]，[89] 中提出，并被 [90]，[91]，[92]，[93]，[94]，[95] 等多次讨论。这个自由飞行问题的主要目的是验证所提出的梁格式的有益属性，并将结果与文献中的数据进行比较。我们采用 [91]，[92] 中选择的参数，见表 1。最初在倾斜配置下静止，φh(s=0,t=0)=(6,0,0) 和 φh(s=L,t=0)=(0,0,8)，梁在 s=L 时受到外部力和扭矩的作用，这两者都作用了一段时间，即 n(s=L,t)=ν(t) 和 m(s=L,t)=μ(t)，其中 (74)ν(t)=f(t)11000,μ(t)=f(t)0112,f(t)=80t，对于 2.5≤t≤5，400?80t，对于 2.55。在 t=5 时，系统关闭，因为两个边界输入都变为零，梁继续自由飞行，应保持总能量、线性和角动量。所得到的运动与 [91]，[93]，[94] 中的描述非常吻合，见图 4 中的快照。此外，从图 5 可以看出总能量的守恒。在每个时间步长，外部载荷所做的功 Wextn=hu?n+1/2?y?n+1/2 完全匹配加载阶段的能量增加 H?(x?n+1)?H?(x?n)。相应地，在所有时间，能量平衡违反 ?ΔEn?H?(x?n+1)?H?(x?n)?Wextn 保持在机器精度范围内，确认了离散时间功率平衡 (44)。加载阶段后的总线性和角动量守恒可以从图 6 中看出，确认了附录 D 中的结果。此外，对于 t≥5，p=(50,0,0) 的总线性动量与分析解完全吻合，在机器精度范围内。此外，负载函数 (74) 以及选择的边界条件允许在运动期间对梁的质心位置进行解析求解 [90]，[92]，由 (75)r=r1(t)04,r1(t)=3+215t3，对于 2.5≤t≤5，436?5t+2t2?215t3，对于 2.55。这与我们方法获得的模拟质心 ?rh(t)?1L∫Ωφh(s,t)ds 完全匹配，见图 7。在加载阶段仍然存在一些微小的偏差 ?Δri?ei?rh?ri，但在 t=5 之后，所有三条曲线都在计算机精度范围内。此外，进行了收敛性研究，以调查自由端中心线位置和速度的时间离散化精度。图 8 中的结果显示，随着时间分辨率的提高，相对误差 (76) 的减小了两个数量级。这与中点型离散化的预期行为一致。这里，评估时间选择为 t?=5。此外，h=10?3 获得的精细离散化参考解。在这种情况下，所有模拟的牛顿容差都降低到了 10?8。

时间积分方法的另一个显著优点是它的稳定性和鲁棒性。理论上，时间步长不受积分器稳定性的限制，而只受牛顿方法由于在过大步长时初始猜测不准确而导致的恶化影响。实际上，该方法必须适当地解决加载函数轮廓 (74)。值得注意的是，选择 h=2.5 对所提出的方法来说没有困难，而其他模拟该示例的工作通常报告的时间步长为 h=0.01 或 0.02。然而，鉴于前面的收敛性分析，对于更大的时间步长，结果的准确性确实应该受到质疑。根据[19]，我们假设一个不可伸长的Kirchhoff梁，并施加相关的剪切刚度和不可伸长性的运动约束。为此，我们在合规方程(18)中将CN设为0，从而确保Γ?=0。此外，请注意，由于所采用的能量-动量时间积分方法（42）的稳健性，我们可以选择相对较大的时间步长h=10^-3（见第一个示例中的段落）。与[19]中的原始设置相比，其中最大的时间步长为h=4?10^-5，即使对于h=5?10^-2，我们也能获得结果。梁最初处于静止状态，并沿着e1方向放置，然后在s=L处受到外力和扭矩的作用，这两者都作用了一小段时间，即n(s=L,t)=ν(t)和m(s=L,t)=μ(t)，其中(77)ν(t)=f(t)0，μ(t)=f(t)0.2500，f(t)=1-√(1-π^2t)，对于0≤t≤0.05，以及f(t)=1/2-√(1-2πt^2)，对于t>0.05。
下载：下载高分辨率图像（228KB）
下载：下载全尺寸图像

图9. 悬臂的非线性振动：截面坐标系中的尖端速度 Vitip。我们的方法()和参考解([19])的结果一致。与参考解相比，符号的变化和分量的不同顺序是由于导向框架方向的不同。
下载：下载高分辨率图像（199KB）
下载：下载全尺寸图像

图10. 悬臂的非线性振动：惯性坐标系φitip中的尖端位置。我们的方法()和参考解([19])的结果一致。与参考解相比，符号的变化和分量的不同顺序是由于导向框架方向的不同。
下载：下载高分辨率图像（134KB）
下载：下载全尺寸图像

图11. 悬臂的非线性振动：正交性约束 gih，其中i∈{...}。

下载：下载高分辨率图像（203KB）
下载：下载全尺寸图像

图12. 悬臂的非线性振动：从位移量计算出的剪切和膨胀应变分量的L2范数，即||Γi||Ω，其中i∈{...}。从应力量计算出的曲率和扭转应变测量的L2范数，以及从位移量计算出的L2范数，表示为||ΔKj||Ω，对于j∈{...}。

我们可以观察到材料尖端速度?Vitip?ei?R(dh(s=L,t=tn))?vφh(s=L,t=tn)（见图9）和空间尖端中心线位置?φitip?ei?φh(s=L,t=tn)（见图10）与[19]的结果有很好的一致性。此外，根据(48)可以验证，在整个仿真过程中，中心FE节点的导向的正交性约束（2）?gih?gi(dh(s=L/2,t=tn))在数值上为零（见图11）。这些结果也适用于其他节点。此外，还可以验证运动学的准确性。从混合变量导出的膨胀和剪切应变分量完全相同。通过(4)从位移计算出的这些应变也与施加的运动学约束保持一致（见变分恒等式(54)）。曲率和扭转分量同样被一致地近似。图12中显示的结果是用空间域上的L2范数表示的，即?||Γi||Ω?(∫0LΓi(qhn,?sqhn)2ds)^(1/2)和?||ΔKj||Ω?(∫0L((Kj)n?Kj(qhn,?sqhn))2ds)^(1/2)。

最后，我们考虑了第5.1节中介绍的具有m=1粘性分支和τE=τG=τ=0.08作为松弛时间的粘弹性广义Maxwell模型。刚度参数在两个分支上一致分配，纯弹性分支占1/4，粘性分支占3/4。请注意，在这个具有粘性效应的仿真中，我们仍然保持上述的运动学假设，因此只有弯曲和扭转模式被阻尼，且VΓ^-1=0。相应地，图13扩展了第6.1节中关于飞意大利面问题的发现，以验证功率平衡方程的一致离散时间表示。总能量随时间减少，耗散能量?D?n?hz?(x?n+1/2)?Rz?(x?n+1/2)严格为正。因此，图13右侧的子图显示，离散时间能量平衡(44)的违反?ΔEn?H?(x?n+1)?H?(x?n)?Wextn?D?(x?n)在数值上为零。

下载：下载高分辨率图像（207KB）
下载：下载全尺寸图像

图13. 悬臂的非线性振动：总能量H?(xn)带有耗散（ ）和无耗散（ ）。每个时间步长的严格正耗散能量和能量平衡违反

6.3. 准静态悬臂问题
这个数值实验展示了所提出的公式在准静态使用案例中的适用性，取自[21]、[86]。我们展示了一个初始时为直线且固定的梁，其长度L=2π，初始配置与e1轴对齐。其他参数详见表3。我们考虑了一个在自由端s=L处具有点力和力矩的载荷情况，即n(s=L,t)=ν(t)和m(s=L,t)=μ(t)，其中(78)ν(t)=t0?P0，μ(t)=t0√2.5P，其中t∈[0,1]。在这个背景下，时间t被视为载荷因子（更多解释见[61]）。相应地，惯性参数被设置为零以实现准静态仿真，这些仿真仅发生在e1-e2平面上。我们区分了两种情况：一种是剪切变形且可伸长的（“无约束”），另一种是剪切刚性且不可伸长的（“有约束”）。对于后者，可以给出欧拉弹性体的解析解[86]。两种约束和无约束设置的结果分别显示在图14中，并与参考结果进行了比较。这些结果不仅证实了正确的实现，还证实了没有锁定效应。

表3. 准静态悬臂问题：物理和数值参数

htendne??LρAMρ11=Mρ22
kb1=kb2
ktks1=ks2
keP10^-2
1
8
10^-1
2
π
0
0
2
0.5
1
5
10
kb1/L2

下载：下载高分辨率图像（258KB）
下载：下载全尺寸图像

图14. 准静态悬臂问题：对于tn∈{0,0.1,0.2,0.4,1}，我们的结果（,）与[21]()中报告的无约束设置的结果以及[86]中的解析解（）一致。尖端中心线沿虚线（,）移动。

6.4. 软体机器臂的动态操控
在这个最后一个例子中，我们模拟了一个在底部固定的软体机器臂，该臂通过气压腔室驱动（见第5.2节）。悬臂梁的横截面是圆形的，直径d=0.03，我们假设材料参数对应于Silicone DragonSkin? 30，即杨氏模量E=6?10^5和剪切模量G=2?10^5 [25]。这对应于泊松比ν=0.35。此外，横截面的面积惯性矩I1=I2=1/2IT=πd^4/64。梁通过三个气压腔室驱动，其质心线平行于梁的中心线，距离为rk=6.5?10^-2，对于所有k都成立。相对于全局e1轴，这些腔室最初位于角度αk∈{π/6, 5π/6, 9π/6}，从而产生截面坐标（ρk1,ρk2）=(rkcosαk,rksinαk)对于k=1,2,3。我们规定气压力（79）pk(t)Ak=1/2f(t)[1+cos(?(t)?αk)]，其中时间依赖函数（80）f(t)为：
f(t)=1/2fmax(1?cos(πt-t1))，对于0≤t≤t1，
fmax,t1其中，我们设置fmax=?50，t1=0.5，t2=3.5，t=4。这些函数在图15中可视化。从e3轴开始，软体机器臂首先在e1-e3平面上进行弯曲运动，直到t=t1。然后是旋转运动，尖端位置?φtip?φh(s=l,t=tn)大致保持在一个恒定高度，见图16中的右侧描绘。在t=t2时，臂完成一次旋转并再次弯曲直到t=t。图16中的快照仅为了清晰起见跟踪了圆形运动阶段。

表4. 软体机器臂的动态操控：物理和数值参数

htendne??lρamρ11=mρ22
kb1=kb2
ktks1=ks2
ke0.05
4
10
10^-1
1
0.175
5
10
8
π
d
2/4
ρieigi

下载：下载高分辨率图像（268kb）
下载：下载全尺寸图像

图15. 软体机器臂的动态操控：相位角函数、振幅函数和相应的气压力。尖端位置分量φitip=ei?φtip，其中i∈{...}。

下载：下载高分辨率图像（289kb）
下载：下载全尺寸图像

图16. 软体机器臂的动态操控：以(55,15)的方位角和仰角视角拍摄的中心线快照。尖端位置φtip沿着圆形路径移动（）。对于修改后的输入函数（81），尖端沿心形路径移动（）。

这种方法可以进一步扩展，使尖端在e1-e2平面上跟随更复杂的路径，例如心形路径。为此，对于?∈[0,2π]，心形路径的坐标为x(?)=sin(?)和y(?)=cos(?)?cos(2?)。相对于原点的相应距离由r(?)=x(?)^2+y(?)^2给出，并可以标准化为r?(?)=r(?)/max(?∈[0,2π]r(?)。通过?=2πtt和??(?)=atan2(y(?)，我们可以最终在(79)中指定振幅和相位角函数（81）为：f(t)=fmaxr?^(2πtt)，?(t)=??^(2πtt)。

soft robotic arm的尖端随后呈现出心形轨迹，见图16。通过包括肌腱和气压腔室，或者将梁分成多个具有独立驱动输入力的部分，可以实现更复杂的变形模式。

7. 结论与展望
在这项工作中，我们提出了一种用于空间cosserat杆动力学的混合公式和结构保持离散化程序。所提出的框架是客观的、无锁定的，并且由于其端口哈密顿结构而具有许多有利特性。通过合规方程，也可以在当前框架内无缝模拟剪切刚性和不可伸长梁的公式。采用的ph形式自然导致了一种结构保持的混合fe方法。在这种设置下，得到的ph微分-代数方程允许基于纯中间点近似的特别简单的时间积分策略。该方案在能量和动量上是一致的，提供二阶精度，并且非常稳健。另一个结果是关键运动关系的准确表示。合规方程和正交性约束都已对时间进行了微分，以揭示ph结构。所提出的时间积分方案准确积分了这些关系。为了方便起见，我们将自己限制在初始时为直线且无应力、横截面几何形状和材料参数恒定的梁问题上。将所提出的框架扩展到更一般的情况将仅仅增加更多的代数工作。虽然所提出的公式具有完整的几何非线性，即允许在三维中的大位移和旋转，但我们始终假设线性本构关系。当前框架可以沿着[41]、[42]的路线扩展到一般的超弹性，这需要包含离散梯度时间积分方法，这与ph结构[41]、[62]、[63]、[64]相吻合。控制导向的驱动机制和粘弹性材料行为已经直接纳入ph框架中，使其非常适合有趣的应用领域[4]、[5]。未来的研究方向之一是将我们的模型纳入新颖的控制策略中，如通过互连或能量整形方法[56]、[57]、[80]。通过采用无奇点的导向参数化，规避了有限旋转数值方法的陷阱。尽管概念上的简单性，导向公式并没有在有限元素内部正确再现旋转，因为一般情况下元素内部的导向正交性条件被违反了。将来，可以使用投影单元四元数[14]、[17]、[20]、[22]来获得cosserat杆动力学的流形一致离散ph公式，但代价是额外的非线性和配置依赖的质量矩阵。

credit作者贡献声明
philipp l. kinon：撰写——原始草稿、可视化、验证、软件、方法论、资金获取、形式分析、数据整理、概念化。
simon r. eugster：撰写——审阅与编辑、验证、监督、资源、方法论、形式分析。
peter betsch：撰写——审阅与编辑、验证、监督、资源、项目管理、方法论、资金获取、形式分析、概念化。 其中，我们设置fmax=?50，t1=0.5，t2=3.5，T=4。这些函数在图15中可视化。从e3轴开始，软体机器臂首先在e1-e3平面上进行弯曲运动，直到t=t1。然后是旋转运动，尖端位置?φtip?φh(s=L,t=tn)大致保持在一个恒定高度，见图16中的右侧描绘。在t=t2时，臂完成一次旋转并再次弯曲直到t=T。图16中的快照仅为了清晰起见跟踪了圆形运动阶段。 表4. 软体机器臂的动态操控：物理和数值参数 htendne??lρamρ11=Mρ22 kb1=kb2 ktks1=ks2 ke0.05 4 10 10^-1 1 0.175 5 10 8 π d 2 4 ρieigi 下载：下载高分辨率图像（268kb） 下载：下载全尺寸图像 图15. 软体机器臂的动态操控：相位角函数、振幅函数和相应的气压力。尖端位置分量φitip=ei?φtip，其中i∈{...}。 下载：下载高分辨率图像（289kb） 下载：下载全尺寸图像 图16. 软体机器臂的动态操控：以(55,15)的方位角和仰角视角拍摄的中心线快照。尖端位置φtip沿着圆形路径移动（）。对于修改后的输入函数（81），尖端沿心形路径移动（）。 这种方法可以进一步扩展，使尖端在e1-e2平面上跟随更复杂的路径，例如心形路径。为此，对于?∈[0,2π]，心形路径的坐标为x(?)=sin(?)和y(?)=cos(?)?cos(2?)。相对于原点的相应距离由r(?)=x(?)^2+y(?)^2给出，并可以标准化为r?(?)=r(?)/max(?∈[0,2π]r(?)。通过?=2πtT和??(?)=atan2(y(?)，我们可以最终在(79)中指定振幅和相位角函数（81）为：f(t)=fmaxr?^(2πtT)，?(t)=??^(2πtT)。 soft robotic arm的尖端随后呈现出心形轨迹，见图16。通过包括肌腱和气压腔室，或者将梁分成多个具有独立驱动输入力的部分，可以实现更复杂的变形模式。 7. 结论与展望 在这项工作中，我们提出了一种用于空间cosserat杆动力学的混合公式和结构保持离散化程序。所提出的框架是客观的、无锁定的，并且由于其端口哈密顿结构而具有许多有利特性。通过合规方程，也可以在当前框架内无缝模拟剪切刚性和不可伸长梁的公式。采用的ph形式自然导致了一种结构保持的混合fe方法。在这种设置下，得到的ph微分-代数方程允许基于纯中间点近似的特别简单的时间积分策略。该方案在能量和动量上是一致的，提供二阶精度，并且非常稳健。另一个结果是关键运动关系的准确表示。合规方程和正交性约束都已对时间进行了微分，以揭示ph结构。所提出的时间积分方案准确积分了这些关系。为了方便起见，我们将自己限制在初始时为直线且无应力、横截面几何形状和材料参数恒定的梁问题上。将所提出的框架扩展到更一般的情况将仅仅增加更多的代数工作。虽然所提出的公式具有完整的几何非线性，即允许在三维中的大位移和旋转，但我们始终假设线性本构关系。当前框架可以沿着[41]、[42]的路线扩展到一般的超弹性，这需要包含离散梯度时间积分方法，这与ph结构[41]、[62]、[63]、[64]相吻合。控制导向的驱动机制和粘弹性材料行为已经直接纳入ph框架中，使其非常适合有趣的应用领域[4]、[5]。未来的研究方向之一是将我们的模型纳入新颖的控制策略中，如通过互连或能量整形方法[56]、[57]、[80]。通过采用无奇点的导向参数化，规避了有限旋转数值方法的陷阱。尽管概念上的简单性，导向公式并没有在有限元素内部正确再现旋转，因为一般情况下元素内部的导向正交性条件被违反了。将来，可以使用投影单元四元数[14]、[17]、[20]、[22]来获得cosserat杆动力学的流形一致离散ph公式，但代价是额外的非线性和配置依赖的质量矩阵。 credit作者贡献声明 philipp l. kinon：撰写——原始草稿、可视化、验证、软件、方法论、资金获取、形式分析、数据整理、概念化。 simon r. eugster：撰写——审阅与编辑、验证、监督、资源、方法论、形式分析。 peter>其中，我们设置fmax=?50，t1=0.5，t2=3.5，t=4。这些函数在图15中可视化。从e3轴开始，软体机器臂首先在e1-e3平面上进行弯曲运动，直到t=t1。然后是旋转运动，尖端位置?φtip?φh(s=l,t=tn)大致保持在一个恒定高度，见图16中的右侧描绘。在t=t2时，臂完成一次旋转并再次弯曲直到t=t。图16中的快照仅为了清晰起见跟踪了圆形运动阶段。

表4. 软体机器臂的动态操控：物理和数值参数

htendne??lρamρ11=mρ22
kb1=kb2
ktks1=ks2
ke0.05
4
10
10^-1
1
0.175
5
10
8
π
d
2/4
ρieigi

下载：下载高分辨率图像（268kb）
下载：下载全尺寸图像

图15. 软体机器臂的动态操控：相位角函数、振幅函数和相应的气压力。尖端位置分量φitip=ei?φtip，其中i∈{...}。

下载：下载高分辨率图像（289kb）
下载：下载全尺寸图像

图16. 软体机器臂的动态操控：以(55,15)的方位角和仰角视角拍摄的中心线快照。尖端位置φtip沿着圆形路径移动（）。对于修改后的输入函数（81），尖端沿心形路径移动（）。

这种方法可以进一步扩展，使尖端在e1-e2平面上跟随更复杂的路径，例如心形路径。为此，对于?∈[0,2π]，心形路径的坐标为x(?)=sin(?)和y(?)=cos(?)?cos(2?)。相对于原点的相应距离由r(?)=x(?)^2+y(?)^2给出，并可以标准化为r?(?)=r(?)/max(?∈[0,2π]r(?)。通过?=2πtt和??(?)=atan2(y(?)，我们可以最终在(79)中指定振幅和相位角函数（81）为：f(t)=fmaxr?^(2πtt)，?(t)=??^(2πtt)。

soft robotic arm的尖端随后呈现出心形轨迹，见图16。通过包括肌腱和气压腔室，或者将梁分成多个具有独立驱动输入力的部分，可以实现更复杂的变形模式。

7. 结论与展望
在这项工作中，我们提出了一种用于空间cosserat杆动力学的混合公式和结构保持离散化程序。所提出的框架是客观的、无锁定的，并且由于其端口哈密顿结构而具有许多有利特性。通过合规方程，也可以在当前框架内无缝模拟剪切刚性和不可伸长梁的公式。采用的ph形式自然导致了一种结构保持的混合fe方法。在这种设置下，得到的ph微分-代数方程允许基于纯中间点近似的特别简单的时间积分策略。该方案在能量和动量上是一致的，提供二阶精度，并且非常稳健。另一个结果是关键运动关系的准确表示。合规方程和正交性约束都已对时间进行了微分，以揭示ph结构。所提出的时间积分方案准确积分了这些关系。为了方便起见，我们将自己限制在初始时为直线且无应力、横截面几何形状和材料参数恒定的梁问题上。将所提出的框架扩展到更一般的情况将仅仅增加更多的代数工作。虽然所提出的公式具有完整的几何非线性，即允许在三维中的大位移和旋转，但我们始终假设线性本构关系。当前框架可以沿着[41]、[42]的路线扩展到一般的超弹性，这需要包含离散梯度时间积分方法，这与ph结构[41]、[62]、[63]、[64]相吻合。控制导向的驱动机制和粘弹性材料行为已经直接纳入ph框架中，使其非常适合有趣的应用领域[4]、[5]。未来的研究方向之一是将我们的模型纳入新颖的控制策略中，如通过互连或能量整形方法[56]、[57]、[80]。通过采用无奇点的导向参数化，规避了有限旋转数值方法的陷阱。尽管概念上的简单性，导向公式并没有在有限元素内部正确再现旋转，因为一般情况下元素内部的导向正交性条件被违反了。将来，可以使用投影单元四元数[14]、[17]、[20]、[22]来获得cosserat杆动力学的流形一致离散ph公式，但代价是额外的非线性和配置依赖的质量矩阵。

credit作者贡献声明
philipp l. kinon：撰写——原始草稿、可视化、验证、软件、方法论、资金获取、形式分析、数据整理、概念化。
simon r. eugster：撰写——审阅与编辑、验证、监督、资源、方法论、形式分析。
peter betsch：撰写——审阅与编辑、验证、监督、资源、项目管理、方法论、资金获取、形式分析、概念化。>