1. 课程导入与弹性常数概述

在材料力学与固体物理中,弹性常数是描述材料在外力作用下发生可逆形变能力的一组核心参数。它把宏观的“软硬”“刚柔”转化为可计算、可比较、可传递的物理量。对工程结构而言,弹性常数直接决定构件在载荷下的位移、固有频率、屈曲临界值与热应力水平;对材料研究而言,它又是判断相稳定性、键合强度和各向异性程度的重要依据。本课程的目标,是建立一条从弹性理论基础到第一性原理计算、再到结果判读的完整链路,使学习者能够独立完成常见晶体弹性常数的计算与验证。

1.1 弹性常数的定义与物理意义

弹性常数最一般的定义来自广义胡克定律。在弹性极限内,应力张量 σij 与应变张量 εkl 之间近似满足线性关系:

σ_ij = C_ijkl · ε_kl        (爱因斯坦求和约定)
ε_ij = S_ijkl · σ_kl
C_ijkl = C_klij = C_jikl = C_ijlk
S = C⁻¹

其中 Cijkl 称为弹性刚度常数(单位 GPa 或 N/m²),Sijkl 称为弹性柔度常数。四阶张量原本有 81 个分量,但应力张量与应变张量的对称性使 Cijkl = Cjikl = Cijlk,能量二阶导的对称性又给出 Cijkl = Cklij,于是独立分量最多降到 21 个。采用 Voigt 记法把 (11,22,33,23,13,12) 映射为 (1,2,3,4,5,6),即可写成 6×6 的 Cαβ 矩阵。

晶体对称性会进一步削减独立分量:三斜晶系 21 个,单斜 13 个,正交 9 个,三方与四方各 6~7 个,六方 5 个,立方 3 个(C11、C12、C44),各向同性体只有 2 个(如体积模量 K 与剪切模量 G)。因此,判断一个晶体属于哪个点群,往往等价于判断其弹性矩阵中哪些元素为零、哪些元素相等。这也是第一性原理计算中“先定对称性、再定独立分量”的根本原因。

弹性常数的物理意义可以从三个层面理解。第一,它是能量对形变的二阶导数:Cαβ = (1/V)·∂²E/∂εα∂εβ,反映晶格抵抗形变的能力。第二,它决定声速与晶格动力学行为,例如立方晶体的 [100] 纵波声速 vL = √(C11/ρ),[100] 横波声速 vT = √(C44/ρ)。第三,它通过 Voigt–Reuss–Hill 平均与多晶弹性模量、泊松比、德拜温度、硬度趋势等工程指标建立联系。

1.2 工程应用场景

弹性常数在工程中并非抽象参数,而是直接进入设计公式。例如:

  • 结构刚度与轻量化:比刚度 E/ρ 决定航空、汽车构件的减重潜力,需要准确的 Cij 来推算各向异性板壳的弯曲刚度。
  • 热障涂层与薄膜:涂层与基体的弹性失配会引发界面应力,Cij 是热应力与剥落寿命评估的输入。
  • 相变与稳定性:立方相向四方相转变时,若 C11−C12 趋近于零,说明剪切失稳临近,可用于筛选形状记忆合金与铁电材料。
  • 高温合金与陶瓷:弹性常数随温度下降的速率,是判断高温承载能力与热震抗性的关键。
  • 地质与高压物理:地幔矿物在高压下的 Cij 决定地震波速,是解释地球内部结构的基础数据。

传统上这些数据依赖超声脉冲回波、共振法或纳米压痕实验。但实验面临单晶样品难制备、微区与薄膜测量误差大、高温高压条件苛刻、含缺陷样品数据分散等问题。这正是第一性原理计算能够补位的场景。

1.3 第一性原理计算的角色

第一性原理(ab initio)方法从电子结构出发,不依赖经验势参数,通过求解 Kohn–Sham 方程获得体系总能量与应力。计算弹性常数的核心思路是:对晶胞施加一系列微小应变,拟合总能量或应力对应变的响应。常用两条路径:

  1. 能量–应变法:施加应变 ε,计算总能 E(ε),用二次多项式拟合 E(ε) = E₀ + ½V₀Cε²,二阶系数即弹性常数。适合小应变、高精度场景。
  2. 应力–应变法:施加应变后直接读取维里应力 σ,由 σ = Cε 线性拟合得到 C。一次计算可同时获得多个分量,效率更高。

以立方晶体为例,只需三组独立形变即可解出 C11、C12、C44:

体积形变:ε₁=ε₂=ε₃=δ,得 K = (C₁₁+2C₁₂)/3
四方剪切:ε₁=ε₂=δ,ε₃=-2δ,得 (C₁₁-C₁₂)
单斜剪切:ε₄=ε₅=ε₆=δ,得 C₄₄

实际计算中,应变幅度通常取 ±0.5%、±1.0%、±1.5%、±2.0%,需保证最大应变仍在弹性范围内;平面波截断能、k 点网格、电子自洽收敛阈值必须经过收敛测试,否则拟合出的 Cij 会出现系统性偏差。此外,离子位置需在每次施加应变后弛豫,而晶胞形状固定,这是“应力–应变法”与“能量–应变法”共同遵循的原则。

1.4 课程路线图

本课程共 30 章,整体按“理论—方法—实操—判读—拓展”推进。前段建立弹性理论与张量基础,中段进入 DFT 计算流程,包括结构优化、应变施加、能量/应力提取与拟合,后段覆盖各晶系实例、多晶模量平均、温度与压力效应、常见报错与结果验证。下图给出课程主线与关键环节的依赖关系。

弹性常数第一性原理计算 · 课程主线 ① 理论基础 应力/应变张量 Voigt 记法与对称性 ② DFT 基础 Kohn–Sham 方程 赝势 / 截断能 / k 点 ③ 计算流程 结构优化 → 施加应变 能量/应力提取 → 拟合 ④ 结果判读 稳定性判据 多晶模量平均 核心计算环节(以立方晶体为例) 结构优化 晶格常数 a₀ 原子位置弛豫 施加应变 δ = ±0.5% ~ ±2% 固定晶胞形状 能量/应力 E(ε) 或 σ(ε) 维里应力张量 拟合求解 二次/线性拟合 C₁₁, C₁₂, C₄₄ 验证 Born 判据 文献对比 输出与应用 弹性矩阵 C_αβ 各向异性数据 多晶模量 K, G, E, ν 声速与德拜温度 热力学性质 相稳定性判据 材料筛选

需要特别强调的是,第一性原理计算给出的弹性常数是 0 K、无缺陷、理想单晶的结果。要与室温实验对比,还需考虑零点振动、热膨胀与声子贡献,这部分将在后续温度效应章节展开。对工程应用而言,0 K 数据常用于趋势筛选与相对比较,而绝对数值需结合准简谐近似或分子动力学修正。

学习建议:在进入具体计算前,先手算一遍立方晶体的 C11、C12、C44 与 K、G、E、ν 的换算关系,并记住 Born 稳定性判据:立方晶体需满足 C11−C12 > 0、C11+2C12 > 0、C44 > 0。这三条不等式是判断计算结果是否物理合理的底线。