📚 课程目录第2章应力应变理论基础

2. 应力应变理论基础

第一性原理计算弹性常数的全部数学基础，都建立在连续介质力学对“应力”与“应变”的张量描述之上。理解这套符号体系，是后续把能量-应变关系、应力-应变关系转化为可计算量的前提。本节从定义出发，逐步建立应力张量、应变张量、广义胡克定律，以及弹性刚度矩阵与柔度矩阵的完整表述，并说明它们各自的对称性来源。

一、应力张量：一点处九个分量的完整描述

在连续介质中，物体内部某一点的受力状态不能用一个标量或矢量完整描述，而必须用二阶张量。取一个以该点为中心、各面法线分别平行于坐标轴 x、y、z 的微小立方体，每个面上都有一个面力矢量。规定：正应力以拉为正、压为负；剪应力分量 σ_ij 表示作用在法线为 i 方向的面上、指向 j 方向的应力分量。于是九个分量构成应力张量：

        ⎡ σ_xx  σ_xy  σ_xz ⎤
σ_ij =  ⎢ σ_yx  σ_yy  σ_yz ⎥
        ⎣ σ_zx  σ_zy  σ_zz ⎦

由微元体绕各轴的力矩平衡（角动量守恒），可证明 σ_ij = σ_ji，即应力张量是对称张量。这一对称性把独立的应力分量从 9 个减少到 6 个：σ_xx、σ_yy、σ_zz、σ_yz、σ_xz、σ_xy。在 Voigt 记法中，常把这 6 个分量排成列向量 σ = (σ₁, σ₂, σ₃, σ₄, σ₅, σ₆)，对应关系为 1→xx、2→yy、3→zz、4→yz、5→xz、6→xy。这一约定在后续刚度矩阵的指标缩并中会反复使用，必须与计算软件的输出约定保持一致。

二、应变张量：几何变形的对称描述

应变描述的是物体内各点相对位移的几何变化。设位移场为 u(x, y, z)，小变形假设下，应变张量定义为位移梯度的对称部分：

ε_ij = ½ ( ∂u_i/∂x_j + ∂u_j/∂x_i )

展开后，对角分量 ε_xx、ε_yy、ε_zz 为正应变，表示沿坐标轴方向的相对伸长；非对角分量 ε_xy、ε_yz、ε_xz 为工程剪应变的一半，即张量剪应变。由于定义中已显式取对称，应变张量天然满足 ε_ij = ε_ji，同样只有 6 个独立分量。这里有一个极易出错的细节：工程剪应变 γ_xy = 2ε_xy，在把第一性原理计算得到的应变施加到晶胞时，若把 γ 直接当作张量分量使用，会导致剪应变被放大一倍，进而使算出的弹性常数出现系统性偏差。因此，在构造形变矩阵时必须明确区分“张量剪应变”和“工程剪应变”。

三、广义胡克定律与刚度、柔度矩阵

在弹性范围内，应力与应变成线性关系，这就是广义胡克定律。用张量形式写作：

σ_ij = C_ijkl · ε_kl        （刚度形式）
ε_ij = S_ijkl · σ_kl        （柔度形式）

其中 C_ijkl 为四阶弹性刚度张量，S_ijkl 为四阶弹性柔度张量，二者互为逆张量。由于 σ 和 ε 各自只有 6 个独立分量，四阶张量可缩并为 6×6 的矩阵。采用 Voigt 记法后，广义胡克定律写成紧凑的矩阵形式：

σ = C · ε        ε = S · σ        S = C⁻¹

以立方晶系为例，其刚度矩阵只有三个独立常数 C₁₁、C₁₂、C₄₄，矩阵形式为：

      ⎡ C11  C12  C12   0    0    0  ⎤
      ⎢ C12  C11  C12   0    0    0  ⎥
C =   ⎢ C12  C12  C11   0    0    0  ⎥
      ⎢  0    0    0   C44   0    0  ⎥
      ⎢  0    0    0    0   C44   0  ⎥
      ⎣  0    0    0    0    0   C44 ⎦

柔度矩阵 S 与 C 结构相同，只是独立常数变为 S₁₁、S₁₂、S₄₄。工程上常用的弹性模量（杨氏模量 E、剪切模量 G、泊松比 ν）都可以由 S 矩阵分量换算得到，例如立方晶系沿 [100] 方向的杨氏模量 E = 1/S₁₁，剪切模量 G = 1/S₄₄。这说明：第一性原理计算直接得到的是 C 矩阵，而工程可读的模量需要先求逆得到 S，再做换算。

四、对称性的三个层次

弹性张量的对称性并非人为约定，而是有明确的物理与数学来源，可分为三个层次。第一层是应力张量和应变张量自身的对称性，即 σ_ij = σ_ji、ε_ij = ε_ji，这来自角动量守恒和应变定义的对称化，使 9 个分量降为 6 个。第二层是弹性张量的“次对称性”：C_ijkl = C_jikl = C_ijlk，即前两个指标和后两个指标分别可交换，这直接继承自 σ 和 ε 的对称性，把 81 个分量压缩到 36 个。第三层是“主对称性”：C_ijkl = C_klij，即指标对 (ij) 与 (kl) 可整体交换。这一条来自应变能密度函数的存在性——弹性变形是保守过程，存在应变能 U = ½ C_ijkl ε_ij ε_kl，对其求二阶偏导可得 ∂²U/∂ε_ij∂ε_kl = ∂²U/∂ε_kl∂ε_ij，从而证明主对称性。主对称性把 36 个分量进一步压缩到 21 个独立常数，这是三斜晶系（最低对称性）的上限。

在此基础上，晶体的点群对称性会进一步减少独立常数数目：单斜晶系 13 个，正交晶系 9 个，四方晶系 6 个（部分点群为 7 个），立方晶系 3 个，六方晶系 5 个，各向同性体只有 2 个。这一递减规律是判断计算是否合理的强约束：如果对立方晶体算出的 C 矩阵中出现了本应为零的 C₁₆ 或 C₄₅ 分量，且数值不可忽略，通常意味着形变施加方式、原子弛豫设置或对称性识别环节存在问题。

五、从张量到可计算量：Voigt 缩并的实操要点

把上述理论落到第一性原理计算中，核心操作是：对晶胞施加一组独立的应变模式，计算对应的应力或能量变化，再拟合出 C 矩阵。以立方晶系为例，常用三种形变模式分别提取三个独立常数。第一种是体积形变（静水应变）：ε₁ = ε₂ = ε₃ = δ，其余为零，此时应力响应给出 C₁₁ + 2C₁₂。第二种是沿 [100] 的单轴应变：ε₁ = δ，ε₂ = ε₃ = 0，配合第一种模式即可分离出 C₁₁ 和 C₁₂。第三种是纯剪切：ε₆ = γ（即 ε_xy = γ/2），直接给出 C₄₄。每种模式取若干应变幅度（例如 δ = ±0.005、±0.010、±0.015），对应力-应变数据做线性拟合，斜率即为相应的弹性常数组合。

在施加剪切应变时，务必确认代码中写入的是张量剪应变 ε_xy 还是工程剪应变 γ_xy。若形变矩阵按 γ 构造，而拟合公式按 ε 推导，C₄₄ 会被系统性低估为真实值的一半。建议在脚本中显式注释所用约定，并用一个已知弹性常数的参考体系做一次自检。

六、结果合理性判断

得到 C 矩阵后，不能只看数值大小，还要做三项检查。其一，对称性检查：对照该晶系的标准矩阵形式，确认本应为零的分量是否接近零（通常要求绝对值小于 C₁₁ 的 1%～2%）。其二，力学稳定性检查：C 矩阵必须正定，即所有特征值为正；对立方晶系，等价条件为 C₁₁ > 0、C₄₄ > 0、C₁₁ > |C₁₂|、C₁₁ + 2C₁₂ > 0。其三，物理量检查：由 S = C⁻¹ 换算出的杨氏模量、剪切模量和泊松比应落在该材料大类的合理区间内，例如金属的泊松比通常在 0.25～0.45 之间，若算出负值或超过 0.5，说明拟合或形变设置存在错误。这三项检查构成了从数学表述走向可信计算结果的基本闭环。