基于混合Hu-Washizu型公式的Port-Hamiltonian表述方法及有限弹性的结构保持离散化

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Port-Hamiltonian formulation and structure-preserving discretization of finite elasticity based on a mixed Hu-Washizu-type formulation

　　莫里茨·希勒 | 彼得·贝奇 | 马龙·弗兰克
德国卡尔斯鲁厄理工学院（KIT）力学研究所

**摘要**
我们提出了一种基于Port-Hamiltonian公式和结构保持离散化的有限弹性理论。能量泛函（或哈密顿量）基于存储能量的多凸表示，并产生了三个应变型场，这些场在Port-Hamiltonian公式中充当能量变量。我们展示了如何使用Livens变分原理的Hu-Washizu型扩展来（i）推导连续的Port-Hamiltonian公式以及（ii）进行结构保持的空间离散化。特别是，我们证明了基础混合公式的空间有限元离散化会产生一个离散的Port-Hamiltonian系统。此外，基础连续公式的时间离散化产生了一个新的能量-动量一致框架，该框架可以容纳替代的有限元公式。这个新框架特别适用于处理准不可压缩材料的行为。通过数值示例评估了新设计的能量-动量方案的数值性能和稳定性。

**1. 引言**
与多物理系统相关的演化方程的Port-Hamiltonian（pH）公式为建模、仿真和控制提供了一个有吸引力的框架[1]，[2]。基于能量的pH描述通过端口考虑了耗散和与环境的相互作用，因此可以被视为经典哈密顿系统向开放系统的扩展。
有限维pH系统最初在[3]中引入，它们基于Dirac结构，该结构将系统属性（如无源性和能量一致性）编码到系统方程中。pH建模范式向无限维系统的扩展基于Stokes-Dirac结构，可以追溯到[4]。从那时起，分布式pH系统在多物理系统建模、控制技术分析和结构保持离散化方法等领域受到了广泛研究，如综合文献综述[5]中所分类的。
设计结构保持离散化技术是保持pH系统内在结构和属性（如无源性和能量一致性）的自然目标。基于Galerkin的数值方法用于pH系统的结构保持离散化最近才变得流行[6]。特别是，有限元方法已被应用于离散化粘弹性弦[7]、梁[8]、[9]、[10]、[11]、板[12]、波动方程[13]、多孔弹性[14]和热弹性[15]的pH公式。值得注意的是，pH方法中能量变量的典型使用本质上产生了混合有限元公式。
有趣的是，非线性弹性动力学领域直到最近才在pH框架内得到关注。在[16]中提出的所谓速度-应力公式的pH结构导致了一种线性隐式能量一致时间步进方案的设计。同样，在[17]中开发的线性隐式能量-动量（EM）方案基于依赖于位移、速度和应力场近似的混合有限元描述。这两项工作[16]、[17]都与Hellinger-Reissner变分公式密切相关，并且仅限于St. Venant-Kirchhoff材料。在[18]中提出的pH公式基于位移、线性动量密度和Green-Lagrangian应变场的近似，并考虑了一般超弹性本构律。虽然[18]中的半离散运动方程保留了pH结构，但在时间离散化中使用了St?rmer-Verlet方案或中点规则，从而通常会损害离散模型中的能量平衡定律。
上述在pH公式框架内的发展补充了早期关于非线性弹性动力学EM方案设计的工作。实际上，对数值稳定时间步进方案的探索导致了在开创性工作[19]中提出的第一个非线性弹性动力学的EM方案。该方法是在纯位移基有限元的框架内开发的，并依赖于本构层面的隐式中点规则的修改。然后，这种修改被推广到任意超弹性材料，应用了离散导数的概念[20]。本质上，在本构层面应用离散导数产生了EM方法的一个显著特征[21]。应当注意的是，EM方案也可以通过时间上的Galerkin方法来设计[22]。特别是在[23]中设计了更高阶的EM方案。在混合有限元的框架内也开发了EM方案，这些方案基于超弹性材料的多凸表示[24]。
在目前的工作中，我们提出了一种新的基于pH的途径来设计有限变形弹性动力学的EM方案。类似于[24]，我们从超弹性材料的多凸表示开始。所提出的基于pH的方法使得能够容纳替代的混合有限元公式，从而产生相关的EM方法。通过将[25]、[26]中开发的混合有限元家族纳入当前方法，可以证明这些混合有限元在包括准不可压缩材料行为在内的有限变形问题中表现良好。
本文的其余部分结构如下：第2节提供了有限应变弹性动力学的混合变分原理，从中将推导出连续的pH公式。完全混合公式的相应弱形式为第3节讨论的能量-动量一致时间离散化提供了起点。第4节表明，完全混合弱形式的空间离散化产生了适合离散pH系统框架的半离散运动方程。从基础完全混合描述的时间离散弱形式出发，第5节纳入了替代的空间离散化，从而得到相应的EM方案。第6节提出了数值研究，第7节得出了结论。

**2. 非线性弹性动力学的Port-Hamiltonian公式**
本节介绍了所提出的非线性弹性动力学pH公式的理论框架。我们首先回顾了连续介质力学和张量叉积的基础知识。随后引入了多凸能量密度函数，并使用变分原理推导出控制方程的强形式。局部场方程被重新表述为Port-Hamiltonian形式，并简要讨论了替代的简化阶pH公式。我们继续推导弱形式，并通过验证能量和角动量的平衡定律来结束。

**2.1. 连续介质力学量和张量叉积**
我们考虑一个在参考配置Ω0?R3中的可变形体，如图1所示，该可变形体由材料点X∈R3描述。映射x(t)=φ(X,t)将点X映射到时间t时体在空间配置Ω中的对应点x=x(t)∈R3，即Ω(t)=φ(Ω0,t)?R3，其中I=[t0,tN]是感兴趣的时间区间。域的边界?Ω0被分解为具有预设位置φˉ(t)的Dirichlet边界?φΩ0和具有预设牵引力tˉ=tˉ(t)的Neumann边界?σΩ0。此外，还有预设在体上的体力gˉ=gˉ(t)∈R3和分布式体控制力bˉ=bˉ(t)∈R3。

**下载：** 下载高分辨率图像（138KB）
**下载：** 下载全尺寸图像
图1. 可变形体Ω的配置和运动学

为了将材料线元素dX映射到它们的空间对应物dx，我们引入了变形梯度Fφ，这是一个二阶张量，表示为
(1)
dx=FφdX,
Fφ=?φ(X,t),
其中?表示相对于材料坐标的梯度。注意，符号Fφ明确表示变形梯度依赖于变形映射。这种符号在后续内容中用于区分依赖于变形的量和混合公式中的相关量。
在本文中，我们使用了两个二阶张量A、B∈R3×3之间的所谓张量叉积
(2)
(A××B)ij=?iαβ?jabAαaBβb,
其中重复的指标a、b、α、β、i、j∈{1, 2, 3}表示根据爱因斯坦约定进行求和，?ijk表示三阶置换张量。张量叉积可以用来重写二阶张量A的余因子和行列式
(3)
cof(A)=1/2A××A,
det(A)=1/6(A××A):A。
这种重新表述由于张量叉积的适当性质而允许一些直接的计算。

**2.2. 多凸能量密度函数和运动方程**
多凸性是广泛接受的超弹性材料在大变形范围内的概念。因此，我们假设所考虑的材料行为可以通过多凸能量密度函数来描述
(8)
W=W^(Fφ)=W?(Fφ, cof(Fφ), det(Fφ)),
其中W?被构造为其参数的凸函数。为了获得一个客观的表示，W被重新表示为
(9)
W=W^(Cφ)=W?(Cφ, Gφ, Jφ),
其中运动学量Cφ、Gφ和Jφ在(6)和(7)中定义。
为了推导弹性体运动的控制方程，我们采用了基于变分泛函的Hu-Washizu型扩展的Livens原理，定义为
(10)
S=∫t0tN∫Ω0[1/2ρ0v·v?W?(C,G,J)+π·(φ˙?v)?ΛC:(Cφ?C)?ΛG:(Gφ?G)?λ(Jφ?J)]dVdt+∫t0tN[∫Ω0(gˉ+bˉ)·φdV+∫?σΩ0tˉ·φdA]dt,
其中ρ0=ρ0(X)是参考配置中的质量密度。这里，位移场φ∈Vφ在函数空间中寻找
(11)
Vφ={φ:Ω0→R3|φi∈H1(Ω0)∧φ=φˉ?φΩ0},
其中H1(Ω0)是Sobolev空间。同样，速度场v∈Vv，其中函数空间Vv定义为
(12)
Vv={v:Ω0→R3|vi∈L2(Ω0)},
其中L2(Ω0)是Ω0上的平方可积函数空间。运动学变量C∈VC、G∈VG和J∈VJ被引入作为表示右Cauchy-Green张量、其余因子和变形梯度行列式的独立场。它们与依赖于变形的量Cφ、Gφ和Jφ的等价性通过张量值拉格朗日乘数ΛC∈VC、ΛG∈VG和标量值拉格朗日乘数λ∈VJ来强制。混合变量的函数空间及其拉格朗日乘数定义为
(13)
VC=VG={A:Ω0→S|Aij∈L2(Ω0)},
VJ={J:Ω0→R+?{0}|J∈L2(Ω0)},
其中S是对称的3×3矩阵空间。注意，在泛函S中，速度场v∈Vv被引入为一个独立场，它通过向量值拉格朗日乘数π∈Vv与变形场φ˙的时间导数相关联。
强加泛函S的稳定性产生了相关的Euler-Lagrange方程。因此，我们得到π=ρ0v，使得乘数π对应于线性动量密度。除此之外，剩余的稳定条件为
(14)
φ˙=v,
ρ0v,
v˙=Div(FφS)+gˉ+bˉ,
ΛC=?CW?(C,G,J),
ΛG=?GW?(C,G,J),
λ=?JW?(C,G,J),
C=Cφ,
G=Gφ,
J=Jφ}
在Ω0×I中，其中第二个Piola-Kirchhoff应力张量由
(15)
S=2ΛC+2ΛG××Cφ+λJφGφ给出。
注意，在(14)中，由于在离散化过程中我们更倾向于使用速度场，因此选择了速度场而不是线性动量密度场。
除了局部场方程(14)之外，泛函S的稳定性还意味着边界条件
(16)
FφSN=tˉ?σΩ0,
其中N表示边界?σΩ0上的单位外法向量。注意，在强加泛函S的稳定性时，已经考虑了所有X∈Ω0在时间t0和tN处的变形场的通常端点条件。

**注释1**
泛函S在(10)中可以被视为Livens变分原理（见[29]）和Hu-Washizu原理（见[30]中的“广义原理”）的组合，适用于多凸弹性。除了拉格朗日乘数π和{ΛC, ΛG, λ}之外，还可以使用另一个向量值拉格朗日乘数来强制在?φΩ0上实现Dirichlet边界条件φ=φˉ（参见[30]）。这与我们当前的方法不同，在当前方法中，通过选择φ∈Vφ来考虑Dirichlet边界条件。

**2.3. 局部场方程的Port-Hamiltonian形式**
接下来，我们将局部场方程转换为Port-Hamiltonian形式。为了识别Port-Hamiltonian公式相关的能量变量，我们引入了弹性体的总能量，也称为哈密顿泛函。因此，我们定义
(17)
H(x)=∫Ω0(1/2ρ0v·v+W?(C,G,J)?gˉ·φ)dV。
注意，体力gˉ的势能已包含在哈密顿量中，而作用在弹性体边界上的外部牵引力tˉ以及分布式控制体力bˉ被视为Port-Hamiltonian公式中的输入。公式(17)的特定形式产生了能量变量
(18)
x=[φvCGJ]。
现在，哈密顿量的泛函导数可以写为
(19)
δxH(x)=[?gˉρ0v?CW??GW??JW?]。
现在的目标是将当前问题的局部场方程转换为pH形式，即
(20)
Ex˙=J(x)z(x)+Bu,y=B*z(x),
其中
(21)
ETz(x)=δxH(x)
（见[31]，[32]）。特别是，我们需要从局部场方程中识别微分算子J(x)以及矩阵E。注意，B是一个可以解释作用在系统上的分布式控制的微分算子。共位和共轭输出y在(20)的第二个方程中定义，其中B*表示B的伴随矩阵。向量值函数z(x)在方程(21)中根据哈密顿量的变分导数定义。首先我们将局部场方程(14)重写为：
(22)
φ˙=vρ0
v˙=Div(FφS)+gˉ+bˉ
ΛC=?CW?(C,G,J)
ΛG=?GW?(C,G,J)
λ=?JW?(C,G,J)
C˙=ddt
Cφ=2sym(FφT?v)
G˙=ddt
Gφ=2Cφ××sym(FφT?v)
J˙=ddt
Jφ=Jφ?1
Gφ:sym(FφT?v)
}
在Ω0×I中，
注意第二个Piola-Kirchhoff应力张量S已在(15)中定义，因此是根据拉格朗日乘数{ΛC, ΛG, λ}的本构关系得出的。(22)中的最后三个方程包含了pH形式(20)所需的速率形式下的运动学关系。现在通过将(22)中的拉格朗日乘数表达式代入(15)，可以直接将其转换为pH形式(20)。因此，我们得到了微分算子J(x)，其表达式为：
(23)
J(x)=[0I
000
?I0
Div(2Fφ□)
Div(2FφCφ××□)
Div(Jφ?1FφGφ□)
0
2sym(FφT?(□))
0
000
0
Cφ××sym(FφT?(□))
0
000
Jφ?1Gφ:sym(FφT?(□))
0
000，
以及矩阵E和B，它们的形式为：
(24)
E=ET=[I
000
0
0
0
0
0
0
0
0
I
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
对时间导数求导并应用链式法则得到：
$$
\frac{d}{dt}U = \int_{\Omega} \int_{0}^{\partial C} w_C \hat{h} \, dV + \int_{\Omega} \int_{0}^{\partial G} w_G \hat{h} \, dV + \int_{\Omega} \int_{0}^{\partial J} w_J \hat{h} \, dV = \left\{ C^\prime \right\} V \int_{\Omega} \int_{0}^{N_C} \partial C \hat{h} \, dV + \left\{ G^\prime \right\} V \int_{\Omega} \int_{0}^{N_C} \partial G \hat{h} \, dV + J^\prime \int_{\Omega} \int_{0}^{N_J} \partial J \hat{h} \, dV
$$
因此，离散存储能量函数 $U$ 在 (76) 中定义的偏导数分别为：
$$
\frac{\partial C^U}{\partial t} = \int_{\Omega} \int_{0}^{N_C} \partial C \hat{h} \, dV, \quad \frac{\partial G^U}{\partial t} = \int_{\Omega} \int_{0}^{N_C} \partial G \hat{h} \, dV, \quad \frac{\partial J^U}{\partial t} = \int_{\Omega} \int_{0}^{N_J} \partial J \hat{h} \, dV
$$
为了避免繁琐的符号表示，$\frac{\partial C^U}$ 和 $\frac{\partial G^U}$ 分别代表 $\frac{\partial U}{\partial C^V}$ 和 $\frac{\partial U}{\partial G^V}$。考虑到对测试函数和试验函数的上述近似，弱形式 (Wequation.30) 可以写为：
$$
\int_{\Omega} w_\pi h \cdot (f^{\prime}h - v_h) \, dV = 0 \quad \int_{\Omega} w \phi h \cdot \rho_0 v^{\prime}h + 2 \text{sym}(F\phi hT \nabla w\phi h) : (\Lambda_C h + \Lambda_G h) \times C\phi h + \lambda h^2 J\phi hG\phi h) \, dV = \int_{\Omega} w \phi h \cdot (g^\prime + b^\prime) \, dV + \int_{\partial \Omega} w \phi h \cdot t^\prime \, dA \quad \int_{\Omega} w_C h : (\Lambda_C h - \partial C \hat{h}) \, dV = 0 \quad \int_{\Omega} w_G h : (\Lambda_G h - \partial G \hat{h}) \, dV = 0 \quad \int_{\Omega} w_J h \left( \lambda h - \partial J \hat{h} \right) \, dV = 0 \quad \int_{\Omega} w \Lambda_C h : (C^\prime h - 2 \text{sym}(F\phi hT \nabla v_h) \, dV = 0 \quad \int_{\Omega} w \Lambda_G h : (G^\prime h - 2 C\phi h \times \text{sym}(F\phi hT \nabla v_h) \, dV = 0 \quad \int_{\Omega} w \lambda h \left( J^\prime h - J\phi h - 1 G\phi h : \text{sym}(F\phi hT \nabla v_h) \, dV = 0
$$
将特定的有限元近似值代入 (79)，我们得到以下方程组：
$$
\begin{align*}
f^{\prime} - v^ &= 0 \\
w^T f^{\prime} &= M_f \rho v^^{\prime} + \Phi(f^) \{ \Lambda^C \} V + \Psi(f^) \{ \Lambda^G \} V + \Xi(f^) \lambda^ - g - b - t &= 0 \\
&= M_f \left\{ \Lambda^C \} V \int_{\Omega} N_C \partial C \hat{h} \, dV + \left\{ G^\prime \} V \int_{\Omega} N_C \partial G \hat{h} \, dV + J^\prime \int_{\Omega} N_J \partial J \hat{h} \, dV \\
\end{align*}
$$
其中，$M_f$, $M_g$, $M_j$, $\Phi(f^)$, $\Psi(f^)$, $\Xi(f^)$, $g$, $b$, $t$ 分别定义为：
$$
\begin{align*}
M_f &= \int_{\Omega} \rho_0 N_C \phi_t N_C \phi \, dV, \quad M_g &= \int_{\Omega} \int_{0}^{N_C} N_C \phi_t N_J \, dV, \quad M_j &= \int_{\Omega} \int_{0}^{N_J} N_J \phi_t N_J \, dV, \\
\Phi(f^) &= \int_{\Omega} 2 B \phi^T N_C \phi \, dV, \quad \Psi(f^) &= \int_{\Omega} 2 B \phi^T D(C\phi h) N_C \phi \, dV, \\
\Xi(f^) &= \int_{\Omega} B \phi^T J\phi h - 1 \{ G\phi h \} N_J \, dV, \quad g &= \int_{\Omega} \int_{0}^{N_C} N_C \phi_t g^\prime \, dV, \\
b &= \int_{\Omega} \int_{\partial \Omega} N_C \phi_t t^\prime \, dA
\end{align*}
$$
在这些定义中，$B_{\phi^} \in \mathbb{R}^{6 \times 3}$ 是通过以下公式定义的经典节点算子矩阵：
$$
S(f^) := A \left\{ V \right\} T B_{\phi^} w^ \phi
$$
为了完整性，节点算子矩阵包含在附录 A.1 中。矩阵 $D(\Box) \in \mathbb{R}^{6 \times 6}$ 对应于 Voigt 符号的张量叉积，满足以下关系：
$$
Q = \{ O \} V^T D(Q) \{ P \} V
$$
其中 $O$ 是一个对称的类似应变的二阶张量，$P$ 是一个对称的类似应力的二阶张量。$D(\Box)$ 的显式表示可以在附录 A.2.4 中找到。

**空间离散化表述为离散的 Port-Hamiltonian 系统**
接下来我们展示从 (80) 导出的空间离散化表述可以转化为离散的 pH 系统。为此，我们考虑离散测试函数的节点值的任意性，并在元素级别对拉格朗日乘数 $\{\Lambda^C\}V$, $\{\Lambda^G\}V$ 和 $\lambda^$ 进行静态凝聚。相应地，(84) 可以重写为：
$$
\begin{align*}
\{\Lambda^C\}V &= M_c^{-1} \frac{\partial C^U}{\partial t}, \quad \{\Lambda^G\}V &= M_c^{-1} \frac{\partial G^U}{\partial t}, \quad \lambda^ = M_j^{-1} \frac{\partial J^U}{\partial t}, \\
(80) &= \phi^{\prime} = v^ \quad M_f \rho v^^{\prime} = - \Phi(f^) M_c^{-1} \frac{\partial C^U}{\partial t} - \Psi(f^) M_c^{-1} \frac{\partial G^U}{\partial t} - \Xi(f^) M_j^{-1} \frac{\partial J^U}{\partial t} + g + b + t \quad M_c \{ C^\prime \} V = \Phi^T(f^) v^ \quad M_g \{ G^\prime \} V = \Psi^T(f^) v^ \quad M_j \j^ = \Xi^T(f^) v^ \\
\end{align*}
$$
注意，在元素级别，$\{C^\prime\}V$, $\{G^\prime\}V$ 和 $J^\prime$ 可以使用 (85) 的最后三个方程用节点量 $\phi^$ 和 $v^$ 表示。这意味着在时间离散化之后，通过元素级别的静态凝聚可以消除与应变类型场 $Ch$, $Gh$ 和 $Jh$ 相对应的自由度（见第 5.1.2 节）。然而，在进行时间离散化之前，我们首先展示 (85) 可以形成一个离散的 pH 系统。为此，我们将 (85) 重写为：
$$
\begin{align*}
[I_{000000} M_f \rho_{000000} M_c \{000000} M_c \{000000} M_j] \[ \phi^{\prime} v^^{\prime} \{ C^\prime\}V \{ G^\prime\}V J^\prime \] = [0 I_{0000} - I_0 - \Phi(f^) - \Psi(f^) - \Xi(f^) 0 \Phi^T(f^) 0000 \Psi^T(f^) 0000 \Xi^T(f^) 0000 ] \[ - g v^ M_c^{-1} \frac{\partial C^U}{\partial t} M_c^{-1} \frac{\partial G^U}{\partial t} M_j^{-1} \frac{\partial J^U}{\partial t} + [00 II_{0000000} ] \\
\end{align*}
$$
上述常微分方程组可以重写为：
$$
\begin{align*}
E^x^{\prime} = J(f^) z(x^) + B_u^,
\end{align*}
$$
其中
$$
E^z(x^) = \nabla H(x^)
$$
类似于 (20) 和 (21) 中的连续 pH 系统，方程 (87) 和 (88) 构成一个离散的 pH 系统 [31], [32]。可以很容易地观察到结构矩阵 $J$ 是斜对称的。此外，离散哈密顿量 $H(x^)$ 是通过插入之前处理的空间离散近似从哈密顿量 (17) 得出的。因此，
$$
\begin{align*}
H(x^) := H_h = \int_{\Omega} (12 \rho_0 v_h \cdot v_h + W^\hat(C,h,G,h) - g^\prime \cdot \phi^) \, dV = 12 v^T M_f \rho v^ + U(\{ C^ \}V,\{ G^ \}V,J^) - \phi^T g,
\end{align*}
$$
其中离散存储能量 $U$ 已在 (76) 中引入。现在，离散哈密顿量的导数可以计算为：
$$
\begin{align*}
\nabla H(x^) = [- g M_f \rho v^ \frac{\partial C^U}{\partial t} \frac{\partial G^U}{\partial t} \frac{\partial J^U}.
\end{align*}
$$
注意，导数 $\frac{\partial C^U}$, $\frac{\partial G^U}$ 和 $\frac{\partial J^U}$ 已在 (78) 中定义。现在可以很容易地验证 (88) 的有效性。

**注释 2**
在元素级别上的离散 pH 系统会产生一个相应的系统级别的离散 pH 系统，因为 (i) 离散 pH 系统之间的功率保持互连再次产生一个离散的 pH 系统，以及 (ii) 有限元组装过程可以被视为元素贡献的功率保持互连。

**注释 3**
半离散运动方程的 pH 形式保证了空间离散化表述的能量一致性。为了看到这一点，考虑：
$$
\begin{align*}
\frac{d}{dt} H(x^) = \nabla H(x^) \cdot x^^{\prime} = (E^z(x^)) \cdot x^^{\prime} = z(x^) \cdot E^x^{\prime} = z(x^) \cdot (J(f^)z(x^) + B_u^) = z(x^) \cdot B_u^ = y^ \cdot u^,
\end{align*}
$$
其中 $y^ = B^T z(x^)$ 是离散输出。因此，如果系统没有受到外部控制力 $u^$ 的作用，则总能量是守恒的。

**注释 4**
为了使 (87) 和 (88) 形成一个离散的 pH 系统，满足 (88) 是至关重要的。应该注意的是，在空间离散化过程中，弱形式中需要拉格朗日乘数的本构方程的存在。然而，有一个例外：在应变类型混合场 $C$, $G$ 和 $J$ 的常数近似情况下，逐点执行拉格朗日乘数的本构方程是足够的。在这种情况下，$N_C = I_6$ 和 $N_J = 1$，导致 $M_c = V_e I_6$ 和 $M_j = V_e$，其中 $V_e$ 是参考配置中元素的体积。

**5. 空间和时间离散化**
在本节中，第 3 节中介绍的时间离散弱形式 (53) 进一步在空间上离散化。尽管我们将处理其他空间离散化方法，但我们总是从完全混合公式对应的时间离散弱形式 (53) 开始。

**5.1. 基于应变类型场 {C, G, J} 的完全混合元素**
首先，我们应用第 4 节中设计的完全混合有限元公式。因此，时间离散弱形式 (53) 产生：
$$
\begin{align*}
\int_{\Omega} w_\pi h \cdot (f^{n+1}h - f^n \Delta_t n - v^{n+1/2}h) \, dV = 0 \\
\int_{\Omega} w \phi h \cdot \rho_0 v^{n+1}h - v^{n+1} \Delta_t n + 2 \text{sym}(F\phi^{n+1/2}hT \nabla w\phi^) : (\Lambda_C^{n+1}h + \Lambda_G^{n+1}h) \times C\phi^{n+1/2}h + \lambda_h^2 J\phi^{n+1/2}hG\phi^{n+1/2}h) \, dV = \int_{\Omega} w \phi h \cdot (g^\prime^{n+1/2} + b^\prime^{n+1/2) \, dV + \int_{\partial \Omega} w \phi h \cdot t^\prime^{n+1/2}dA \\
\int_{\Omega} w_C h : (\Lambda_C^{n+1}h - D_C \hat{h}) \, dV = 0 \\
\int_{\Omega} w_G h : (\Lambda_G^{n+1}h - D_G \hat{h}) \, dV = 0 \\
\int_{\Omega} w_J h [ \lambda^{n+1}h - D_J \hat{h} \, dV = 0 \\
\int_{\Omega} w \Lambda_C h : (C^{n+1}h - C_n \Delta_t n - 2 \text{sym}(F\phi^{n+1/2}hT \nabla v^{n+1/2}h)] \, dV = 0 \\
\int_{\Omega} w \Lambda_G h : (G^{n+1}h - G_n \Delta_t n - 2 C\phi^{n+1/2}h \times \text{sym}(F\phi^{n+1/2}hT \nabla v^{n+1/2}h)] \, dV = 0 \\
\int_{\Omega} w \lambda_h [ \lambda^{n+1}h - J_n \Delta_t n - (J\phi^{n+1/2}h) - 1 G\phi^{n+1/2}h : \text{sym}(F\phi^{n+1/2}hT \nabla v^{n+1/2}h)] \, dV = 0
\end{align*}
$$
注意，离散导数 $D_ViW^\hat(V_n,V_{n+1}$，$V_i \in \{Ch,G_h,J_h\}$ 已在第 3.5.1 节中引入。

**5.1.1. 完全混合离散 Port-Hamiltonian 公式**
由于完全混合公式的时间和空间离散化是可交换的，因此从 (92) 导出的方案也可以写成：
$$
\begin{align*}
E(x^{n+1} - x^n) = \Delta_t n(J(f^{n+1/2}z^\bar(x^n,x^n+1) + B_u^{n+1/2),
\end{align*}
$$
其中
$$
\begin{align*}
E^{n+1}z^\bar(x^n,x^n+1) = \nabla^\bar H(x^n,x^n+1).
\end{align*}
$$
这里，矩阵 $E$, $J(f^)$ 和 $B$ 可以再次从离散 pH 系统 (86) 中推断出来。此外，离散哈密顿量 (89) 的离散导数由下式给出：
$$
\begin{在两种情况下，对所使用的材料模型关于参考配置进行线性化处理后，可以得到准不可压缩的材料行为，其有效的泊松比分别为ν=0.495和ν=0.499。表2. 来自[36]的Mooney-Rivlin模型参数。表3. 来自[25]的改进Mooney-Rivlin模型参数，有效泊松比为ν=0.495。表4. 受[25]启发的改进Mooney-Rivlin模型参数，有效泊松比为ν=0.499。现在可以轻松计算出完全混合公式（第5.1节）、尺寸缩减混合公式（第5.2节）和尺寸缩减位移公式（第5.3节）中所需的存储能量函数的离散导数。对于Mooney-Rivlin模型（111），我们得到（113）DCW?=aI，DGW?=bI，DJW?=c2(Jn+1?1)2?dln(Jn+1)?(c2(Jn?1)2?dln(Jn))Jn+1?Jn。类似地，改进的Mooney-Rivlin模型（112）的离散导数表示为（114）DCW?=αtr(Cn+1/2)I，DGW?=βtr(Gn+1/2)I，DJW?=?γ(ln(Jn+1)?ln(Jn))+?1(Jn+12?2?Jn2?2+Jn+1?2?2?Jn?2?2)Jn+1?Jn。请注意，如示例1.6.1中所述，所考虑的两种存储能量函数都是可分离的。

这项初步研究的目的是评估新开发的公式和积分器与标准中点积分器相比的鲁棒性、长期稳定性和离散守恒特性。我们考虑了一个L形体的初始边界值问题，参考了[19]、[24]、[37]的工作。几何形状和初始设置以及边界条件如图2所示。特别是，施加的外部牵引力为tˉ1(t)=?tˉ2(t)=f(t)，其中f(t)的定义如下：f(t)={t<2.5s时为5，2.5s5s时为0。L形域被划分为117个有限元素，模拟时间设置为T=100s。图3显示了运动快照以及相应的冯·米塞斯应力分布，清楚地显示了大的旋转、变形和应力。图2展示了L形体的2D几何投影（左）和边界条件的3D示意图（右）。图4显示了使用所研究的积分器得到的总能量。注意，在t=5s之后，系统上没有外力作用，因此总能量应该是一个常数。

本例的目的是研究准不可压缩材料行为下的数值鲁棒性和离散守恒特性。我们主要应用了尺寸缩减混合公式（第5.2节），即RM-Q1Q0-EM和RM-Q2Q1-EM，因为底层混合有限元公式已被证明适用于准不可压缩有限弹性（参见[25]、[26]）。这些示例来源于[38]、[39]。对于纯弯曲梁，也可参考[17]。几何形状和初始设置以及边界条件如图7所示。为了获得与[17]、[39]中相似的变形，我们选择了一个密度（见表3），该密度是[17]、[39]中使用的密度的四分之一。模拟使用了RM-Q1Q0-EM公式处理弯曲梁，RM-Q2Q1-EM公式处理扭转梁。在弯曲梁的情况下，总时间为T=20s，时间步长为Δtn=0.05s，牛顿容忍度为?=5e-5；对于扭转梁，参数设置为T=4s，?=5e-6，时间步长分别为Δtn=0.01s（对于表3中的材料参数）和Δtn=0.002s（对于表4中的材料参数）。图8展示了RM-Q1Q0-EM公式在t={0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6}s时的变形快照和冯·米塞斯应力分布σvM [Pa]。图10显示了表3中参数对应的弯曲（左）和扭转（右）梁的总能量。尽管时间步长相对较粗，但RM-Q1Q0-EM公式在弯曲梁和扭转梁的模拟过程中都保持了能量的守恒。

我们提出了一个新的有限弹性的Port-Hamiltonian公式及其在空间和时间上的结构保持离散化方法。我们展示了Livens变分原理与适应于多凸有限弹性的Hu-Washizu型扩展的结合可以用来（i）推导出连续的pH公式，并（ii）构建一个完全混合的弱形式，这为在空间和时间上的结构保持离散化提供了基础。我们发现，完全混合弱形式的直接空间离散化导致了半离散的运动方程，这些方程可以表示为离散的pH系统。此外，完全混合弱形式的直接时间离散化提供了一个框架，可以容纳不同的有限元公式，并导致结构保持的能量-动量方案。这样，就可以选择适合处理准不可压缩有限弹性的替代有限元并将其扩展到动态范围。特别是对于[25]、[26]中开发的混合有限元公式，这一点得到了证明。底层pH描述的一个独特特点是应变-位移关系以速率形式表示。在这方面，当前方法与[40]和[39]中的先前发展有一些相似之处。当前基于pH的方法的主要优点是它以自然的方式促进了能量-动量一致方法的设计。这已经在使用中点规则和应用于存储能量函数的离散导数的概念的第二阶精确时间步进方案的背景下得到了证明。本质上，应变-位移关系的速率形式使得为整个混合有限元公式家族设计能量-动量一致的方法成为可能。相应的尺寸缩减程序已经为混合有限元素（见[25]、[26]的第5.2节）和纯位移元素（见第5.3节）展示。这与之前的相关工作[24]不同，后者需要特定的级联形式的应变型变量，从而排除了使用其他先前开发的混合有限元进行有限弹性分析的可能性。另一方面，当前基于pH的方法的多功能性是以需要在每个时间步长从存储中检索应变型历史变量为代价的。特别是，在混合有限元素的静态凝聚过程中以及对于完全缩减的位移元素，会出现应变型历史变量。事实上，这一特定特征将基于pH的方法与之前设计的EM方法区分开来。未来的研究可以探讨将当前方法扩展到多物理场以及非弹性材料模型的可能性，进一步探索Port-Hamiltonian框架的可能性。