1. 从张量到向量为什么我们需要Voigt与Mandel标记法如果你处理过固体力学、材料科学或者晶体物理中的本构关系一定会对四阶弹性张量那令人头疼的81个分量感到窒息。在计算机里我们最终需要的是矩阵运算而不是一堆嵌套的索引。这就是Voigt标记法和Mandel标记法诞生的核心驱动力它们是将高阶对称张量尤其是二阶对称张量如应力、应变以及四阶对称张量如刚度、柔度压缩成向量和矩阵的“编码规则”是连接连续介质力学理论与有限元、材料计算等数值实现的关键桥梁。简单来说这两种方法都是为了“降维打击”。它们通过利用张量的对称性将6个独立的应力或应变分量排列成一个6×1的列向量将具有最多36个独立分量的四阶弹性张量排列成一个6×6的矩阵。这样一来复杂的张量双点积运算如广义胡克定律 σ C : ε就优雅地简化为了矩阵乘法运算{σ} [C] {ε}。这不仅仅是书写上的简化更是数值计算和程序实现的基石。没有这种标记法现代有限元分析软件的核心材料子程序将难以高效编写。那么Voigt和Mandel这两种看起来目标一致的标记法究竟有何不同为什么有了Voigt还需要Mandel这绝不仅仅是学术上的文字游戏而是关系到能量表示的一致性、数值计算的精度以及程序实现的便利性。对于从事计算力学、本构建模、材料数据分析和相关软件开发的朋友来说理解两者的区别与联系是写出正确、高效代码的必修课。接下来我将结合十多年的仿真和编程经验为你彻底拆解这两种标记法的“前世今生”与实战应用。2. 基石对称二阶张量的向量化表示在深入对比之前我们必须先统一“原材料”的处理方式。无论是Voigt还是Mandel它们首先都要处理应力和应变这类对称二阶张量。2.1 应力与应变张量的分量在三维空间中一个对称二阶张量例如柯西应力张量 σ有9个分量但由于对称性σ_ij σ_ji只有6个是独立的。我们通常按以下顺序列出它们 σ_11, σ_22, σ_33, σ_23, σ_13, σ_12。 应变张量 ε 同理其独立分量为ε_11, ε_22, ε_33, 2ε_23, 2ε_13, 2ε_12。 注意应变中的“2”因子这是为了在后续的功共轭关系中保持能量的一致性使得应力功密度满足 w σ : ε {σ}^T {ε}。这是一个关键细节很多初学者会在这里出错。2.2 Voigt标记法最直观的“打包”方式Voigt标记法是最古老、最直观的映射方式。它简单地将6个独立分量按顺序排列成一个列向量。对于应力Voigt向量定义为 {σ}_V [σ_11, σ_22, σ_33, σ_23, σ_13, σ_12]^T对于应变为了满足功共轭其Voigt向量定义为 {ε}_V [ε_11, ε_22, ε_33, 2ε_23, 2ε_13, 2ε_12]^T这样应力功密度就可以简洁地表示为 w σ : ε σ_11ε_11 σ_22ε_22 σ_33ε_33 2σ_23ε_23 2σ_13ε_13 2σ_12ε_12 {σ}_V · {ε}_V 这里“·”表示向量的点积。看Voigt标记法成功地将张量双点积转化为了向量点积。注意应变定义的“陷阱”。在不同的教科书或软件中你可能会看到两种应变Voigt向量定义一种是带因子2的工程剪应变如上所述另一种是不带因子2的张量剪应变γ_23, γ_13, γ_12。前者保证了能量形式的简洁w {σ}^T {ε}后者则与张量分量直接对应。在查阅任何资料或编写代码时第一件事就是确认它使用的是哪一种约定。ABAQUS、ANSYS等主流商用软件通常使用带因子2的工程剪应变约定。混淆二者是导致材料子程序计算结果错误的最常见原因之一。2.3 从Voigt到弹性矩阵一个隐藏的问题当我们用Voigt标记法表示四阶弹性刚度张量 C_ijkl 时会得到一个6×6的矩阵 [C]_V。广义胡克定律写作 {σ}_V [C]_V {ε}_V这里看起来一切完美。但是如果我们仔细审视这个矩阵乘法展开后的能量表达式会发现一个不协调之处。能量 w 作为标量应该满足 2w {ε}_V^T [C]_V {ε}_V 由于 [C]_V 是从四阶张量映射而来它本身应该是对称且正定的。然而如果我们直接使用上面的Voigt向量在能量表达式中剪应变分量带因子2的会引入额外的系数使得 [C]_V 矩阵中的某些元素与剪切相关的位置与原始四阶张量分量之间的关系变得“不纯粹”。具体来说矩阵 [C]_V 不再是四阶张量 C 在某种正交变换下的纯粹表示而是一种带有缩放因子的表示。这在数学上不够优雅在涉及张量不变量计算或坐标变换时会带来额外的复杂性。3. Mandel标记法为了能量正交性而生的优化Mandel标记法正是为了解决Voigt标记法在能量表示上的这一“不纯粹”问题而提出的。它的核心思想是对剪应力/剪应变分量进行归一化处理使得映射后的6维向量空间成为一个真正的、具有标准欧几里得内积的向量空间。3.1 Mandel向量的定义Mandel标记法对应力和应变向量采用了相同的归一化因子。其定义如下应力Mandel向量 {σ}_M [σ_11, σ_22, σ_33, √2 σ_23, √2 σ_13, √2 σ_12]^T应变Mandel向量 {ε}_M [ε_11, ε_22, ε_33, √2 ε_23, √2 ε_13, √2 ε_12]^T 请注意这里应变分量不再包含因子2而是统一使用 √2。这个 √2 因子从何而来它源于将一对对称的剪切分量如 σ_23 和 σ_32视为一个整体时在构建正交基过程中的归一化要求。可以证明在这种定义下应力功密度依然满足 w σ : ε {σ}_M · {ε}_M 并且更重要的是Mandel向量空间是一个“等距映射”的空间。这意味着原始张量空间中的Frobenius内积即双点积完全等同于Mandel向量空间中的标准点积没有任何缩放因子的扭曲。3.2 Mandel标记法的核心优势标准正交性这是Mandel标记法最强大的地方。在Mandel表示下6×6的弹性刚度矩阵 [C]_M 满足一个极其优美的性质 [C]_M [M] [C]_V [M]^{-1} 其中 [M] 是一个简单的对角缩放矩阵其对角线元素为 [1, 1, 1, √2, √2, √2]。更关键的是由于映射是等距的四阶弹性张量的所有内在对称性主要对称性 C_ijkl C_klij 和次要对标性 C_ijkl C_jikl C_ijlk在 [C]_M 中表现为矩阵的完全对称性。这使得坐标变换极其方便在Mandel表示下弹性矩阵的坐标变换例如从材料主轴转到全局坐标系可以表示为标准正交矩阵的相似变换。而在Voigt表示下变换矩阵不是正交的形式更复杂容易出错。本征问题更自然求解弹性张量的特征值和特征向量对应于材料的弹性模量和模式时在Mandel表示下就是一个标准的6×6实对称矩阵的特征值问题。在Voigt表示下由于矩阵不是标准内积下的对称矩阵需要谨慎处理。张量运算的保真度更高任何在原始张量空间中成立的恒等式在Mandel向量/矩阵空间中几乎保持相同的形式。这对于推导复杂的本构模型如塑性、损伤的数值算法特别有利。实操心得何时用Voigt何时用Mandel根据我的经验可以遵循一个简单的原则输入输出和存储用Voigt核心计算用Mandel。原因如下Voigt表示直观与大多数实验数据报告格式σ11, σ22, σ33, τ23, τ13, τ12和软件输入输出接口直接兼容。因此从文件读取数据或向用户显示结果时用Voigt。然而在内部进行材料状态更新、计算雅可比矩阵切线刚度、进行坐标变换或任何涉及张量不变量运算时转换为Mandel表示进行计算会简单、清晰且不易出错。许多先进的材料模型源代码内部都采用Mandel约定。4. 实战对比各向同性弹性矩阵的两种表示让我们以一个最简单的例子——各向同性线弹性材料——来直观感受两种标记法的差异。各向同性材料的刚度张量只有两个独立参数常用杨氏模量 E 和泊松比 ν 表示或用拉梅常数 λ 和 μ 表示。4.1 Voigt表示下的弹性矩阵 [C]_V使用拉梅常数 λ 和 μ剪切模量在Voigt标记法下6×6弹性矩阵为[λ2μ, λ, λ, 0, 0, 0] [λ, λ2μ, λ, 0, 0, 0] [λ, λ, λ2μ, 0, 0, 0] [0, 0, 0, μ, 0, 0] [0, 0, 0, 0, μ, 0] [0, 0, 0, 0, 0, μ]注意右下角3×3剪切块的对角线元素是 μ。如果我们用这个矩阵计算应力例如 σ_23 [C]_V(4,4) * ε_V(4) μ * (2ε_23)结果正是 2μ ε_23这与本构关系 σ_23 2μ ε_23 一致。这里可以看到Voigt表示中剪切分量隐含的因子2。4.2 Mandel表示下的弹性矩阵 [C]_M根据Mandel标记法的定义我们需要对Voigt矩阵进行一个缩放变换[C]_M [M] [C]_V [M]^{-1}。由于 [M] 是对角阵 diag(1,1,1,√2,√2,√2)这个变换实际上只影响与剪切相关的行和列。变换后Mandel表示的弹性矩阵为[λ2μ, λ, λ, 0, 0, 0] [λ, λ2μ, λ, 0, 0, 0] [λ, λ, λ2μ, 0, 0, 0] [0, 0, 0, 2μ, 0, 0] [0, 0, 0, 0, 2μ, 0] [0, 0, 0, 0, 0, 2μ]看在Mandel表示下剪切部分的对角线元素变成了 2μ。这个形式更加“干净”和“统一”。整个矩阵可以看作是单位矩阵的某种组合其特征值和特征向量的物理意义更加清晰。例如体积模量 K λ 2μ/3 相关的模式对应前三个分量的某种组合而剪切模量 μ 直接体现在后三个分量的系数上现在是2μ。4.3 坐标变换的复杂度对比假设我们有一个在材料坐标系例如纤维方向下定义的弹性矩阵 [C]_V‘现在需要将其旋转到全局坐标系下的 [C]_V。对于Voigt表示变换公式为 [C]_V [T]_V [C]_V‘ [T]_V^T 这里 [T]_V 是一个6×6的变换矩阵其元素是方向余弦的函数形式相当复杂因为它需要正确处理正应力与剪应力耦合、剪应力之间的耦合。推导和编程实现这个矩阵很容易出错。而对于Mandel表示变换公式在形式上完全相同 [C]_M [T]_M [C]_M‘ [T]_M^T 但关键在于[T]_M 是一个正交矩阵即 [T]_M^{-1} [T]_M^T。这意味着在Mandel空间中坐标变换就是标准的正交变换。其矩阵形式比Voigt的变换矩阵简洁得多并且具有更好的数值性质保持矩阵的范数和条件数。注意事项变换矩阵的编程实现。如果你正在编写一个需要处理材料取向如复合材料的有限元材料子程序强烈建议在内部使用Mandel表示进行坐标变换。你可以先实现Voigt到Mandel的转换函数然后在Mandel空间中进行正交变换最后再转回Voigt格式输出。这会大大减少代码调试的难度。网上有很多开源代码库如NumPy或Eigen的封装提供了现成的Mandel变换工具可以直接参考使用。5. 在有限元与材料子程序中的具体应用理论说得再好不如一行代码。在实际的有限元分析中材料子程序如ABAQUS的UMAT、ANSYS的USERMAT是应用这些标记法的主战场。5.1 UMAT中的应力应变传递约定以ABAQUS为例它传递给UMAT的应变增量数组STRAN和需要UMAT更新的应力数组STRESS默认使用的就是Voigt标记法并且是带因子2的工程剪应变约定。具体顺序通常是(ε11, ε22, ε33, γ12, γ13, γ23)其中γij 2εij。 应力顺序对应为(σ11, σ22, σ33, σ12, σ13, σ23)。因此当你从UMAT的输入参数中读取STRAN时心里要清楚它的后三个分量已经是工程剪应变了。如果你内部采用Mandel表示进行计算第一步就是进行转换ε_M [ε11, ε22, ε33, γ12/√2, γ13/√2, γ23/√2]^T。5.2 雅可比矩阵DDSDDE的计算UMAT要求你提供雅可比矩阵DDSDDE即 ∂Δσ / ∂Δε。这是切线刚度矩阵的Voigt表示。如果你在材料积分算法中在Mandel表示下推导并得到了雅可比矩阵J_M ∂{σ}_M / ∂{ε}_M那么你不能直接将其赋值给DDSDDE必须转换回Voigt表示。转换关系为DDSDDE(i,j) M_ii * J_M(i,j) * (M_jj)^{-1}其中M_ii是之前提到的对角缩放矩阵[M]的第i个对角元。对于i, j 1,2,3M_ii 1对于i, j 4,5,6M_ii √2。由于[M]是对角阵其逆就是每个对角元的倒数。所以具体操作是前3行前3列正应力-正应力部分直接复制J_M(1:3, 1:3)。前3行后3列正应力-剪应力部分DDSDDE(1:3, 4:6) J_M(1:3, 4:6) / √2。后3行前3列剪应力-正应力部分DDSDDE(4:6, 1:3) J_M(4:6, 1:3) / √2。后3行后3列剪应力-剪应力部分DDSDDE(4:6, 4:6) J_M(4:6, 4:6) / 2。忘记这个转换或者转换错误是导致材料子程序收敛缓慢甚至发散的一个非常隐蔽的原因。因为错误的雅可比矩阵会误导求解器的牛顿迭代方向。5.3 一个简单的Mandel工具函数示例为了方便你可以在代码中定义如下工具函数。这里用Python伪代码示意import numpy as np # 定义缩放因子 M_diag np.array([1., 1., 1., np.sqrt(2.), np.sqrt(2.), np.sqrt(2.)]) M_inv_diag 1.0 / M_diag def voigt_to_mandel(vec_v): 将Voigt向量转换为Mandel向量。假设输入vec_v是工程应变Voigt向量。 return M_diag * vec_v # 对于应变Voigt后三个分量是2*epsilon这里乘以M_diag后变为 sqrt(2)*2*epsilon需要根据定义调整。 # 更精确的如果输入是ABAQUS的STRAN (ε11, ε22, ε33, γ12, γ13, γ23) # 其中γ2ε那么要得到Mandel应变 [ε11, ε22, ε33, √2 ε12, √2 ε13, √2 ε23] # 转换应为ε_M [ε11, ε22, ε33, γ12/√2, γ13/√2, γ23/√2] def abaqus_stran_to_mandel(stran): mandel np.zeros(6) mandel[0:3] stran[0:3] mandel[3:6] stran[3:6] / np.sqrt(2.0) return mandel def mandel_to_abaqus_stress(sig_m): 将Mandel应力向量转换为ABAQUS Voigt应力向量。 stress_v np.zeros(6) stress_v[0:3] sig_m[0:3] stress_v[3:6] sig_m[3:6] * np.sqrt(2.0) # 注意ABAQUS应力顺序是(11,22,33,12,13,23) # 如果Mandel顺序是(11,22,33,23,13,12)则需要调整索引对应关系 return stress_v def mandel_jacobian_to_voigt_ddsdde(J_m): 将Mandel表示下的雅可比矩阵J_M转换为ABAQUS所需的DDSDDE (Voigt)。 D np.zeros((6,6)) # 正-正部分 D[0:3, 0:3] J_m[0:3, 0:3] # 正-剪部分 D[0:3, 3:6] J_m[0:3, 3:6] / np.sqrt(2.0) # 剪-正部分 D[3:6, 0:3] J_m[3:6, 0:3] / np.sqrt(2.0) # 剪-剪部分 D[3:6, 3:6] J_m[3:6, 3:6] / 2.0 return D常见问题顺序混淆。不同的软件、甚至不同的文献中Voigt向量的后三个剪切分量的排序可能不同。常见的有(12, 13, 23)和(23, 13, 12)。ABAQUS使用的是前者(12, 13, 23)。Mandel表示通常继承所基于的Voigt顺序。在编写代码或阅读数据时务必首先确认顺序。一个错误的顺序会导致完全错误的结果而且这种错误在简单测试中可能发现不了。6. 高级话题各向异性材料与张量不变量对于各向异性材料如单晶、复合材料或织构金属弹性张量具有更丰富的结构。此时Mandel标记法的优势更加明显。6.1 各向异性弹性矩阵的存储对于最一般的三斜晶系弹性刚度张量有21个独立常数。在Voigt表示下6×6矩阵[C]_V是对称的我们可以只存储上三角的21个元素。在Mandel表示下矩阵[C]_M同样是对称的存储方式相同。但由于Mandel表示是正交变换[C]_M的特征值就是弹性张量的主刚度其特征向量方向具有明确的物理意义对应于特定的应变模式。而在Voigt表示下直接对[C]_V进行特征值分解得到的“特征向量”并不是真正的张量特征模式因为它们存在于一个非标准内积的空间中。6.2 张量不变量的计算在塑性力学、损伤力学中我们经常需要计算应力的不变量如偏应力、等效应力von Mises应力、第三不变量等。这些不变量在张量形式下有清晰的定义。在Voigt表示下计算等效应力 σ_eq √(3/2 * s:s)其中s是偏应力张量。用Voigt向量表示偏应力{s}_V后这个双点积需要小心计算s:s s11^2 s22^2 s33^2 2*(s23^2 s13^2 s12^2)。注意剪切项前面的因子2。在Mandel表示下由于内积是标准的计算变得直接。将应力张量转换为Mandel向量 {σ}_M然后计算偏应力需要先计算静水压力部分。等效应力的平方可以表示为 (3/2) * {s}_M · {s}_M这里的点积就是标准的向量点积无需额外因子。这使得在Mandel空间内实现本构模型更新算法更加简洁和不易出错。6.3 应用于复杂本构模型对于涉及复杂屈服面、塑性流动和硬化法则的材料模型其数值积分算法如回映算法通常需要在应力空间中进行操作。在Mandel表示的6维欧几里得空间中屈服函数可以看作是一个超曲面塑性流动方向是该超曲面的梯度一个6维向量。所有的向量运算加法、点积、范数都是标准的这使得算法的推导和编程实现更加直观。例如在编写一个基于Mises屈服准则的塑性UMAT时如果你在Mandel表示下工作计算屈服函数值f ||s||_M - sqrt(2/3)*σ_y就非常简单其中||s||_M就是Mandel偏应力向量的模长。而塑性流动方向n ∂f/∂σ在Mandel表示下就是n_M s_M / ||s_M||这是一个单位向量。这些操作在Voigt表示下都需要格外注意各个分量的权重。7. 总结与选择建议经过以上长篇的拆解我们可以清晰地看到Voigt和Mandel标记法的定位Voigt标记法是接口标准和存储格式。它直观、历史悠久、被绝大多数商业软件、实验数据报表和教科书采用。当你需要与外部世界输入、输出、显示交换数据时使用Voigt。Mandel标记法是计算引擎和算法空间。它数学性质优美正交性、标准内积使坐标变换、特征值求解、张量运算和复杂本构模型的数值实现变得简洁、清晰且健壮。当你进行核心的材料状态更新、刚度矩阵变换或任何涉及张量代数的运算时在内部使用Mandel。我个人的实战经验是在开发新的材料模型或编写高性能的材料子程序时采用“Mandel内核Voigt外壳”的架构。即在子程序入口将输入的Voigt应变增量转换为Mandel应变增量。在Mandel空间中进行所有的本构积分、应力更新和雅可比矩阵计算。在子程序出口将更新后的Mandel应力转换回Voigt应力输出并将Mandel雅可比矩阵转换回Voigt格式的DDSDDE。这种方法结合了两种标记法的优点既保证了与外部环境的兼容性又享受了内部计算简洁可靠的福利。刚开始接触时你可能会觉得多了一层转换有点麻烦但一旦你经历过在Voigt表示下调试一个复杂的各向异性材料变换矩阵的噩梦你就会深深体会到Mandel标记法带来的清爽与安心。它就像一把精心打磨的瑞士军刀在纷繁的张量计算世界中为你提供了一个标准、统一的操作平台。
