您可以点此 下载链接 下载实验报告的 PDF 扫描版本，或者继续阅读本页面，其由 OCR 技术生成但可能存在识别误差的电子文档，建议以扫描版为准，感谢您的查阅。
一、实验预习
1. 实验目的
（1）以MgO晶体为计算对象，学习在Materials Studio（CASTEP模块）中通过对晶体施加沿X轴（$[100]$方向）的单轴拉伸应变，计算各应变下体系的应力张量，获得完整的单轴应力-应变曲线；
（2）从线性弹性阶段提取MgO的杨氏模量 $E$（沿$[100]$方向）和泊松比 $\nu$，并识别弹性极限点（弹性-非线性转变点），理解DFT所得应力-应变曲线与宏观实验曲线的对应关系；
（3）理解DFT计算应力的理论基础（Hellmann-Feynman定理与应力定理），以及单轴拉伸与弹性常数张量的关系。
2. 实验原理
2.1 应力的DFT计算——应力定理
在DFT框架下，应力张量由总能对应变的变分导出（Nielsen-Martin应力定理）：
$$ \sigma_{\alpha\beta} = -\frac{1}{\Omega}\frac{\partial E_{\text{total}}}{\partial \varepsilon_{\alpha\beta}} $$
其中 $\Omega$ 为晶胞体积，$\varepsilon_{\alpha\beta}$ 为应变张量分量。CASTEP在每次自洽计算后自动输出应力张量的笛卡尔分量（单位GPa）：
$$ \boldsymbol{\sigma} = \begin{pmatrix}\sigma_{xx}&\sigma_{xy}&\sigma_{xz}\\\sigma_{yx}&\sigma_{yy}&\sigma_{yz}\\\sigma_{zx}&\sigma_{zy}&\sigma_{zz}\end{pmatrix} $$
2.2 单轴拉伸应变的施加方式
对平衡结构沿X轴施加单轴拉伸工程应变 $\varepsilon_{xx} = \varepsilon$，晶胞变换矩阵为：
$$ \mathbf{h} = \mathbf{h}_0 \cdot (\mathbf{I} + \boldsymbol{\varepsilon}), \quad \boldsymbol{\varepsilon} = \begin{pmatrix}\varepsilon & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{pmatrix} $$
即将晶格矢量 $a$ 拉伸为 $a(1+\varepsilon)$，而 $b, c$ 方向不施加外部约束，允许弛豫（但在此实验中设置为 $y, z$ 方向是刚性锁定的，仅使用经验泊松比释放应力）。
2.3 线弹性本构关系
在线弹性范围内，单轴应力满足广义胡克定律：
$$ \sigma_{xx} = E \cdot \varepsilon_{xx} $$
实验中通过线性拟合应力-应变数据提取杨氏模量 $E$：
$$ \sigma = E\varepsilon + b $$
截距 $b$ 来源于平衡结构的残余应力（理想情况下 $b \approx 0$）。
对于立方晶体，$[100]$ 方向的杨氏模量与弹性常数的关系为：
$$ \frac{1}{E_{[100]}} = S_{11} = \frac{C_{11}+C_{12}}{(C_{11}-C_{12})(C_{11}+2C_{12})} $$
泊松比由横向收缩与纵向伸长之比定义：
$$ \nu = -\frac{\varepsilon_{yy}}{\varepsilon_{xx}} $$
对立方晶体，理论值为：
$$ \nu = -\frac{S_{12}}{S_{11}} = \frac{C_{12}}{C_{11}+C_{12}} $$
2.4 应力-应变曲线的物理分区
完整应力-应变曲线分为以下阶段：
| 阶段 | 特征 | DFT计算范围 |
|---|---|---|
| 线弹性区 | $\sigma \propto \varepsilon$，可逆 | $\varepsilon \lesssim \varepsilon_y$ |
| 弹性极限点 | 线性偏离起点 $(\varepsilon_y, \sigma_y)$ | 识别拐点 |
| 非线性弹性区 | 仍可逆但非线性 | $\varepsilon_y < \varepsilon < \varepsilon_{\text{fracture}}$ |
| 塑性/断裂区 | 不可逆，宏观实验才能观测 | DFT理想晶体通常直接断裂 |
DFT计算的是理想完整晶体的应力-应变响应，不包含位错、晶界等缺陷效应，因此得到的是理论强度（理论极限值），通常远高于实验测量值。
二、实验过程
1. 实验过程记录
（1）平衡结构准备
以实验一BM方程拟合所得 $a_0 = 4.285$ Å为初始参数，在CASTEP中对MgO进行充分几何优化（GeomOpt），使力收敛至 $< 0.01$ eV/Å，应力收敛至 $< 0.1$ GPa，得到充分弛豫的平衡结构（实验截图中Balanced MgO结构，显示5个晶胞约束，约束条件为 $1\,1\,1\,0\,0\,0$，即固定晶胞形状和体积，仅优化原子位置）。
计算参数：GGA-PBE，截断能370 eV，K点 $3\times3\times3$，Gaussian展宽0.1 eV，Ultrasoft赝势。
（2）单轴应变序列生成
在Materials Studio中，对平衡结构沿X轴（$a$方向）施加一系列拉伸应变，应变值 $\varepsilon$ 从0.005到0.100，步长0.005，共20个应变点。每个应变结构的生成方式为：修改晶格参数，令 $a' = a_0(1+\varepsilon)$，$b' = c' \approx a_0 \times [1 - \nu \cdot \epsilon]$，保持角度90°不变，对每个结构设置晶胞约束（固定晶胞），在CASTEP中提交GeomOpt计算。
（3）应力数据提取
对每个计算完成的晶体结构，从.castep输出文件中提取应力张量的笛卡尔分量（GPa）。实验截图中显示了0.0655 Å应变点的应力张量输出：
Cartesian components (GPa):
x y z
x 10.124799 -0.000015 0.000070
y -0.000015 1.199875 -0.000029
z 0.000070 -0.000029 -4.1749取 $\sigma_{xx}$ 值作为单轴应力数据点，整理全部20个应变点的 $(\varepsilon_i, \sigma_{xx,i})$ 数据，由Python脚本（stress_strain_analysis.py）自动完成拟合、弹性区识别和图形输出。
三、分析讨论
1. 实验结果与分析
MgO单轴拉伸应力-应变曲线：
曲线参数汇总：
$$ a_0 = 4.28513\ \text{Å},\quad \nu = 0.18 $$
线性拟合结果（弹性范围 $\varepsilon \in [0, 0.025]$）：
$$ \sigma = 202.11\,\varepsilon + 0.0727, \quad R^2 = 0.99801596 $$
由此提取杨氏模量：
$$ \boxed{E = 202.11\ \text{GPa}} $$
截距 $b = 0.0727$ GPa接近零，说明弛豫充分，平衡结构残余应力极小。
弹性极限点由分析脚本自动识别：
$$ \boxed{\varepsilon_y = 0.030,\quad \sigma_y = 5.81\ \text{GPa}} $$
当 $\varepsilon > 0.030$ 时，应力偏离线性，进入非线性弹性区（图中蓝色圆点转为空心圆点标注），曲线斜率逐渐降低，应力继续增大但增速减缓。
这与GGA对MgO弹性常数的系统性低估有关，同时也与本次仅对 $b, c$ 方向进行的经验泊松比约束设置有关，若完全放开横向弛豫则 $\nu$ 将更接近理论值。
与实验值对比：
| 参数 | 本次DFT计算 | 实验值 |
|---|---|---|
| $E_{[100]}$ (GPa) | 202.11 | 248–310 |
| $\varepsilon_y$（理论极限） | 0.030 | $\sim 0.001$（实验多晶） |
杨氏模量低于实验值约25–35%，原因在于：GGA-PBE系统性低估结合能，导致能量曲面曲率（即弹性常数）偏低；实验测量的多晶MgO包含织构效应，$[100]$ 单晶方向模量本身也因晶向各向异性而不同于多晶均值。
弹性极限应变（理论值0.030）远高于实验多晶断裂应变（$\sim 0.001$），这是DFT计算理想完整晶体的本征结果——真实材料的断裂由位错滑移、晶界开裂等缺陷机制主导，DFT无法描述这些过程，给出的是不含缺陷的理论强度上限。
2. 问题提出与讨论
问题：DFT给出的应力-应变曲线与宏观实验曲线（附图所示）有何本质差异，弹性极限如何界定？
宏观实验应力-应变曲线（见附图）完整包含线弹性区（Linear Elastic Region）、比例极限（Proportional Limit）、屈服应力（Yield Stress）、完全塑性屈服（Perfectly Plastic Yielding）、应变硬化（Strain Hardening）、颈缩（Necking）直至断裂（Fracture）全过程，并区分工程应力/应变曲线与真实应力/应变曲线（True Stress/Strain）。
DFT计算的理想晶体应力-应变曲线在物理上对应的是完整单晶、零温、无缺陷条件下的响应，具体差异为：
$$ \underbrace{\sigma_y^{\text{DFT}}}_{\text{理论极限}} \gg \underbrace{\sigma_y^{\text{exp}}}_{\text{实验值}} $$
DFT曲线不存在位错介导的屈服平台和应变硬化，超过弹性极限后直接进入非线性响应直至键断裂，因此本质上只有弹性区和断裂前的非线性区两个阶段。DFT弹性极限的界定方式是：以线性拟合的残差为判据，当残差超过阈值（本实验取中位数残差的3.5倍）时认为线性关系不再成立，对应的最后一个线性点即为弹性极限 $(\varepsilon_y, \sigma_y)$。
此外，从实验截图的应力张量输出可以看出，在施加 $x$ 方向单轴应变后，$y$ 和 $z$ 方向出现了非零正应力分量（$\sigma_{yy} \approx 1.20$ GPa，$\sigma_{zz} \approx -4.17$ GPa），这正是泊松效应的直接体现——纵向拉伸引起横向收缩趋势，由于晶胞在 $y, z$ 方向被约束，该收缩趋势转化为横向应力。若要获得严格意义上的单轴应力状态，应允许 $b, c$ 两个方向自由弛豫，使 $\sigma_{yy} = \sigma_{zz} = 0$。
作者：GARFIELDTOM
邮箱：coolerxde@gt.ac.cn