用于有限Cosserat微极弹性理论的SO(3)旋转场的结构保持离散化方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：A structure-preserving discretisation of SO(3)-rotation fields for finite Cosserat micropolar elasticity

　　•有限应变Cosserat弹性理论简述。•Eringen公式与其位错形式之间的新关联概述。•旋转张量离散化与客观性问题的简要回顾。•一种将MITC理念扩展到非线性三维系统的新型结构保持插值方法。•用于验证该方法并展示其作用的全面数值示例。

引言
经典连续介质力学将材料视为由无数点构成的连续体，每个点都具有三个平动自由度[1]。根据角动量守恒定律，柯西应力张量是对称的，即σ=σT，且材料的响应仅取决于变形映射φ:V?R3→Vφ?R3，这一映射通过其一阶梯度F=Dφ来描述[2]。虽然这一框架能够准确描述许多宏观现象，但在应变场变化剧烈的情况下则失效[3]。在应变梯度较高或存在显著微观结构相互作用的区域，如泡沫、复合材料、颗粒介质和超材料中，内部旋转和偶应力变得重要，此时经典理论已不再适用[4]。

在Cosserat微极模型中，自然地考虑了材料点的内在旋转和偶应力[5]。该理论通过为每个材料点引入独立的旋转自由度，丰富了经典连续介质力学，从而能够体现微观结构对宏观响应的影响[6]。在此框架下，连续体中的每一点都由位置向量φ以及一个表示局部微旋转的旋转张量Rˉ:V?R3→SO(3)来描述[7]，该微旋转与宏观旋转R=polarF相互独立。Fran?ois和Eugène Cosserat在1909年的著作《可变形物体理论》中提出的初衷，是构建一种统一的可变形物体场论，使得力的平衡和力矩的平衡都能从欧几里得变换下不变的单一变分原理中得出[8]。因此，Cosserat微极模型可以通过最小化作用量函数A(φ,Rˉ)→minw.r.t.{φ,Rˉ}来方便地描述[7]。

如今，三维Cosserat微极连续介质模型被广泛用于那些由微旋转、偶应力及固有长度尺度决定材料行为的场合[9]。它能有效描述具有弯曲主导型微观结构的构制超材料和手性超材料，能够捕捉与尺寸相关的刚度以及扭转耦合效应[10]、[11]，在地球科学中也被用于模拟地震过程中的颗粒聚集体和断块相互作用[12]。在纳米尺度上，同一框架还可解释诸如中心对称磁性材料中的磁机械耦合现象，即柔顺磁性[13]。对于薄结构，Cosserat壳层理论可将三维模型简化为具有旋转特性的二维模型，从而实现自然的混合维度耦合[14]。在这方面，最近的Γ收敛结果阐明了三维能量与简化后能量之间的关系，并给出了高效的建模和参数识别所需的显式曲率表达式[15]、[16]，进一步的分析还表明通过域重构，该理论也可应用于非定向流形[17]。对于动态Cosserat板模型，人们已经确定了适当的初始条件和边界条件，以确保在薄微结构介质中的瞬态波作用下模型的合理性[18]。不过，该理论最常用的特殊形式是一维Cosserat杆模型，它能够对细长梁进行几何上精确的描述。实际上，Simo–Reissner梁[19]、[20]、[21]就是它的特例，在相同的运动学假设下包含了剪切变形和有限旋转。过去几十年来，Cosserat杆模型在工程和应用物理学领域发挥着核心作用，应用范围涵盖软体机器人中的实时仿真与控制[22]、[23]，生物力学和生物物理学中DNA、胶原蛋白等生物纤维的建模[24]、[25]、[26]，以及海上和能源领域中钻杆、升管和系泊绳的非线性动力学分析[27]、[28]、[29]，此外还有许多其他发展和应用，例如可参考[30]、[31]、[32]、[33]。由于该模型的广泛应用，其在计算领域也受到了高度重视。这里的主要挑战在于，在大角度旋转的情况下保持客观性和几何一致性，同时避免出现膜锁定或剪切锁定之类的伪现象。为了解决这一问题，人们基于李群方法以及SO(3)上的旋转插值，提出了几何上精确的公式，以确保框架的不变性[34]、[35]、[36]，而锁定效应则通过混合公式、选择性积分[37]，以及采用不同试函数和测试空间针对特定能量贡献的Petrov–Galerkin方案得到缓解[38]、[39]。需要指出的是，锁定现象[40]、[41]、[42]可被视为任何维度下Cosserat连续介质所固有的结构保持问题[14]、[43]、[44]。

几何上精确的Cosserat公式利用旋转张量Rˉ:V→SO(3)来描述材料点的独立有限旋转。这些独立旋转的效应由两个张量来衡量[45]。应变由Cosserat应变张量RˉTF?1来表示，它由Biot型伸长张量Uˉ=RˉTF和二阶单位张量1组成。曲率则可以用经典的Eringen形式中的wryness张量Kˉ=axl(RˉTRˉ,i)?ei来表示[46]，或者用位错形式下的微位错张量Aˉ=RˉTCurlRˉ来表示[47]。本质上，Cosserat应变张量的耦合机制存在于所有衍生模型中，正是它导致了相应梁理论和壳层理论中众所周知的锁定现象。但更一般而言，这是一个结构保持问题，可以表述为如何构建一个离散的Rˉh，使其恰好满足极分解式Fh=RhUh所给出的polarFh=Rh的条件。为便于理解，假设某个材料参数μc→+∞，那么由其控制的耦合能量μc‖RˉhTFh?1‖2必须保持有限值。这就意味着‖RˉhTFh?1‖2必须以至少与μc趋于无穷大的相同速率趋近于零。在标准的离散化方法中，变形梯度被定义为Fh=Dφh，其中变形映射φ=φiei是通过经典的C0(V)-连续拉格朗日元素φh∈R3?CGp(V)来逐行近似的，即φi∈CGp(V)?H1(V)。显然，由于p≥1，且CGp(T)=Pp(T)在元素层面是一个完整的多项式空间，其梯度已经包含了常数空间R3×3?D[R3?CGp(T)] [48]。然而，任何离散的Rˉh都必须属于H1[V;SO(3)]的一个共形子空间，因为旋转场在SO(3)中取值，而SO(3)是一个非线性的李群流形，而非向量空间。换句话说，无法逼近单位矩阵1的根本原因在于RˉhT和Fh的离散化之间存在不匹配，这使得无法逐点重建RˉhTDφh=1。要理解这一点，只需应用极分解式，并要求RˉhTFh=RˉhTRhUh=1。既然我们已经知道离散变形梯度Fh能够生成一个常数伸长张量Uh=1，那么正是Rh=polarFh与Rˉh在弱意义上的不相容性导致了我们所称的锁定现象。再次强调，由于SO(3)不是一个代数空间，因此不存在针对它的标准多项式插值方法。所以，所有的插值都依赖于某种从有限元空间到目标空间的映射。实现这种映射的方法有多种，比如利用欧拉–罗德里格斯闭式公式通过指数映射[49]、[50]、[51]、[52]，使用欧拉旋转矩阵[53]、[54]，四元数[55]、[56]、[57]、[58]、[59]，测地线元素[36]、[60]、[61]、[62]、[63]、[64]，或是基于投影的元素[65]、[66]。然而，上述方法单独使用都无法精确匹配φh∈R3?CGp(V)时的polarFh。原因很简单：要么是函数的规律性不足，要么是插值方法有问题，或者两者都有问题。简而言之，由于φ∈R3?H1(V)且Dφ∈R3?H(curl,V)，因此Rˉ∈H1[V;SO(3)]?polar[R3?H(curl,V)]?R=polarF。同时，Rˉh∈Xh(V)?Yh(V)?Rh=polarFh，这里的Xh(V)和Yh(V)是某些共形子空间。在几何线性化模型中，耦合效应可简化为μc‖skewDu?A‖2≈μc‖(1/2)curlu?a‖2，此处u∈R3?CGp(V)是位移场，A∈so(3)?CGp(V)是斜对称线性化旋转矩阵，a=axlA∈R3?CGp(V)是其轴向量，而显然curlu?R3?CGp(V) [43]。最简单的类似情况是当t→0时Reissner–Mindlin板模型中的项(1/t2)‖?w??‖2，此时挠度w∈CGp(A)，旋转量?∈R2?CGp(A)，因为?w?R2?CGp(A) [41]。在所有情况下，根本的难题都是高阶规则的旋转空间与低阶规则的梯度空间之间的不兼容性。从密度角度分析，这种不兼容性问题在无限维度情形下，或者在具有足够丰富结构的有限元空间中，如高阶空间或宏观空间，可能不会出现[67]、[68]。不过，这一观点在实际有限元离散化中使用的离散子空间中通常并不成立。因此，本研究的主题就是重新建立精确的匹配关系。我们提出的解决方案拓展了混合插值张量分量（MITC）方法[69]、[70]、[71]、[72]以及壳体中的Regge插值方法[73]的原有理念，将其应用到更高维度和非线性有限元空间中。具体而言，我们在离散作用量函数中提出如下构造方式：即Rˉh由测地线有限元GEq(V)给出，但对于耦合项，则先通过Πcp?1将其插值到Nédélec空间NIIp?1(V) [48]、[74]、[75]、[76]中，然后再投影回旋转的李群SO(3)中，从而得到一种低阶规则的基于投影的插值方法。由于Fh=Dφh∈R3?NIIp?1(V)?R3?H(curl,V)，其极分解后的结果位于polar[R3?NIIp?1(V)]中，这样就实现了精确匹配。通过使用测地线有限元来离散化微旋转，真实的变动率以及曲率值得以保留，同时插值过程也保持了客观性。鉴于该方法在几何连续介质力学和几何有限元空间方面的理论基础（ΓEΩMETPIA），我们将该方法命名为几何结构保持插值法（Γ-SPIN）。还需指出的是，对于Cosserat杆模型，我们的方法与运用Petrov–Galerkin方法来减轻剪切锁定的技术存在关联[38]、[59]，因为后者通过特定的测试函数选择，同样可以实现消除虚假耦合的投影效果。

本文的结构如下。首先介绍Cosserat模型，并探讨其有限弹性和偶应力理论的极限情况。接着，阐述我们的结构保持方法，包括对SO(3)场离散化以及稳定离散配对选项的讨论。随后，给出数值示例并分析其结果。最后，总结研究结论并展望未来。

本研究全程采用以下符号表示。特殊情况会在具体语境中明确说明。
•我们用{e1,e2,e3}表示笛卡尔坐标基。
•对于重复出现的指标φiei，默认遵循爱因斯坦求和约定。
•向量用加粗的小写字母表示，如v、ξ∈Rd。
•二阶张量用加粗的大写字母表示，如T∈Rd×d。
•四阶张量则用黑板格式表示，如?∈Rd×d×d……
•尖括号用于表示任意维度的标量积，如〈v,u〉=viui，〈T,F〉=TijFij。
•矩阵乘积用于表示高阶张量与低阶张量之间的所有部分缩并，如Tv=Tijvjei，?T=CijklTklei?ej。
•二阶单位张量定义为1=ei?ei，因此1v=v。
•物理域的体积、表面积和曲线分别用V、A和s表示。
•旋转矩阵的空间由SO(d)表示，即SO(d)={R∈Rd×d∣detR=1,R?1=RT}，它与极投影算子polar:Rd×d→SO(d)相关联。
•我们将斜对称二阶张量的常数空间定义为so(d)={T∈Rd×d∣T=?TT}。
•空间so(d)与运算符skewT=(1/2)(T?TT)∈so(d)、Antiv=v×1∈so(3)及其逆运算axl(Antiv)=v相关联。
•勒贝格空间以及相关的希尔伯特函数空间在域V?R3上用常规符号L2(V)、H1(V)表示。
•空间之间的张量积是一种表示多个空间实例的方式，如R3?H1(V)=[H1(V)]3=H1(V;R3)=H1(V)×H1(V)×H1(V)，也就是说，若φiei=φ∈R3?H1(V)，则可知φi∈H1(V)。
•纳布拉算子定义为?=ei?i。
•左梯度用?表示，即?λ=??λ。
•对于向量和高阶张量，右梯度用D表示，即Dv=v??。
•我们定义矢量散度为divv=〈?,v〉。
•张量散度则用DivT=T?表示，即沿行进行单次缩并操作。
•矢量旋度用curlv=?×v表示。
•对于张量，旋度运算符表示为CurlT=?T×?，也是沿行操作。

片段内容
非线性微极模型
在下文中，我们将运用虚功原理简要推导出准静态的、几何非线性的Cosserat微极模型。目前，微极模型主要有两种著名的表述形式，即经典的Eringen形式[46]和位错形式[47]。这两种形式的等价性源于SO(3)矩阵场的梯度完全由其旋度所控制[77]。这两种形式对应的等效材料参数以及相应的转换公式如下。

客观且保持结构的离散化方法
接下来，我们将介绍一种能够保持客观性的作用量函数离散化方法，即便在耦合项中存在投影效应也是如此。我们首先讨论旋转的离散化方法以及客观性保障措施，涉及欧拉–罗德里格斯公式、欧拉矩阵、测地线插值以及基于投影的插值方法。为便于全面介绍这些新提出的理念，我们先回顾一些关于SO(3)张量插值的常见概念。

数值示例
在下文中，我们将通过若干基准测试来验证和量化我们的结构保持方法。为简化实现过程，这里我们采用欧拉矩阵，以避免需要进行乘法更新，同时也无需求解测地线问题，因为后者需要通过复杂的纯算法逐元素求解最小化问题[62]。为了在使用欧拉矩阵的情况下仍保持测试的有效性，我们将基准测试的范围限制在单轴旋转下的适度变化情况。因此，云台锁定问题得到了解决。

**结论与展望**
在本研究中，我们提出了用于有限应变Cosserat微极模型中保持结构特性的几何结构保持插值方法（Γ-SPIN）。该方法主张使用测地线有限元来进行插值运算Rˉh。这样的选择能够在存在刚体运动叠加的情况下保证结果的客观性，并且能够准确测量曲率变化率。当Cosserat耦合模量趋于无穷大时，即μc→+∞，该模型便趋近于其对应的耦合应力模型。

**CRediT作者贡献声明**
Lucca Schek：写作——审阅与编辑，写作——初稿撰写，可视化，验证，研究分析，形式化分析。
Peter Lewintan：写作——审阅与编辑，写作——初稿撰写，研究。
Wolfgang Müller：写作——审阅与编辑，监督指导，资金筹集。
Ingo Muench：写作——审阅与编辑。
Andreas Zilian：写作——审阅与编辑，方法论设计。
Stéphane P.A. Bordas：写作——审阅与编辑，监督指导。
Patrizio Neff：写作——审阅与编辑。

**利益冲突声明**
作者们声明存在以下可能被视为潜在利益冲突的财务利益或个人关系：Adam Sky表示获得了欧盟委员会的财政支持；Lucca Schek表示获得了德国研究基金会的财政支持；Patrizio Neff也表示获得了德国研究基金会的财政支持。其他作者则声明没有已知的财务利益或个人关系可能构成利益冲突。

**致谢**
Adam Sky感谢欧洲委员会在“玛丽·居里行动”博士后奖学金计划中的支持，具体项目为AGE2M，资助编号为101152236。Lucca Schek感谢德国研究基金会在该项目下的支持，项目编号为525235558。Patrizio Neff感谢德国研究基金会“2256优先计划——用于预测工程结构和材料中复杂现象的变分方法”项目中的支持，项目编号为440935806。

Lucca Schek | Peter Lewintan | Wolfgang Müller | Ingo Muench | Andreas Zilian | Stéphane P.A. Bordas | Patrizio Neff | Adam Sky