从VASP弹性常数到连续介质模型
将Si实际弹性张量导出为SI单位，再计算不同厚度悬臂梁的挠度和频率。
阅读本篇目录
VASP给出的弹性常数可以作为更大尺度力学模型的材料输入。本文从集群实算的Si弹性张量出发，导出SI单位材料文件，再用悬臂梁公式计算一个小例子，帮助理解原子尺度参数如何传递到器件尺度。
1. 先得到可用的弹性张量#
沿用Si弹性计算教程的有限应变结果：PBE PAW64、ENCUT=550 eV、EDIFF=1E-6；各应变下允许内部原子弛豫，自由原子力要求0.02 eV/Å。该原子尺度计算经Slurm完成，使用NVHPC25.3/CUDA12.8构建的vasp/6.6.1-p100。
本例输入为内部弛豫后的立方晶体C11、C12、C44，以及同一模型的质量密度。先检查弹性矩阵正定，再求逆得到柔度S，沿[100]方向的杨氏模量为1/S11。
得到E[100]=122.473313 GPa，密度2280.906 kg/m³。这里采用特定晶向的单晶模量，不能用多晶平均值直接替换后仍称为同一个模型。
2. 导出材料参数#
附件已有本例的ELASTIC.json。完成自己弹性教程中的提取后，可以换成自己的同格式结果，再运行：
module load ase/3.29.0
python export_continuum.py
输出CONTINUUM-MATERIAL.json包含6×6矩阵、密度、晶轴和约定：
| 字段 | 本例约定 |
|---|---|
| 弹性矩阵 | Pa；原GPa数值乘10⁹ |
| Voigt顺序 | xx、yy、zz、yz、xz、xy |
| 应变向量 | εxx、εyy、εzz、2εyz、2εxz、2εxy |
| 应力向量 | σxx、σyy、σzz、σyz、σxz、σxy |
| 密度 | kg/m³ |
| 晶向 | 梁长沿Si [100] |
不同有限元程序的剪切项顺序、工程剪切应变与张量剪切应变可能不同。导入前按目标程序调整，不能只复制36个数字。如果器件坐标与晶轴不同，还需旋转四阶弹性张量；本例没有做任意晶向旋转。
3. 计算三个厚度的悬臂梁#
模型长L=1 mm、宽b=100 μm，端部横向载荷F=10 μN，厚度t分别为10、20、40 μm。采用均匀Euler–Bernoulli梁、小变形、固定端无柔度的解析近似：
A = b t
I = b t³ / 12
端部挠度 δ = F L³ / (3 E I)
根部最大弯曲应力 σ = 6 F L / (b t²)
第一阶弯曲频率 f1 = 1.875104² / (2π L²) × sqrt(E I / (ρ A))
脚本同时生成CONTINUUM-DATA.json，数值如下：
| 厚度 / μm | 端部挠度 / μm | 最大弯曲应力 / MPa | 第一阶频率 / kHz |
|---|---|---|---|
| 10 | 3.266018 | 6.000 | 11.837 |
| 20 | 0.408252 | 1.500 | 23.674 |
| 40 | 0.051032 | 0.375 | 47.349 |
在相同端部载荷下，厚度从10增加到20 μm，弯曲刚度增加8倍，挠度降为1/8；第一阶频率增加到2倍。这是给定梁模型的解析结果，其中材料参数来自真实DFT计算。
4. 换成自己的模型时检查什么#
先确认晶向、单位、剪切约定和密度，再确认结构尺度是否适合连续介质近似。本例的δ/L都小于1%，使用小挠度模型；厚梁、复杂夹持、残余应力、表面效应、有限温度或大变形需要更合适的模型。
本文没有进行有限元网格求解，也没有验证真实加工器件的强度或寿命。计算出的最大应力不是材料破坏强度，不能据此自动给出安全载荷。相场、压电或半导体器件还需要自由能、梯度系数、耦合张量等额外参数，单靠弹性常数无法补齐。