首页 > 其他分享> > Lammps例子：弹性常数的计算（转摘）

Lammps例子：弹性常数的计算（转摘）

2022-07-17 05:00:21 作者：互联网

转摘自知乎：https://zhuanlan.zhihu.com/p/359014805?utm_medium=social&utm_oi=869273158199443456。这个解释很清晰，帮助理解。

一、广义胡克定律

材料在线弹性范围内，应力 $\sigma$ 与应变 $\varepsilon$ 成正。由于一些剪切条件的等价性和对称性，对于对称性最差的三斜晶系，独立的弹性常数可减至21个。可用下面的矩阵表示：

$\begin{bmatrix}\sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_4 \\ \sigma_5 \\ \sigma_6 \end{bmatrix}=\begin{bmatrix} c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\ c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\ c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\ c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\ c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\ c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\ \end{bmatrix} \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5 \\ \varepsilon_6 \end{bmatrix}\\$

由于对称性中只需要考虑一个半三角就行。

二、Voigt标记

$\varepsilon = \begin{pmatrix} e_1 & \frac{1}{2}e_6 & \frac{1}{2}e_5 \\ \frac{1}{2}e_6 & e_2 & \frac{1}{2}e_4 \\ \frac{1}{2}e_5 & \frac{1}{2}e_4 & e_3 \end{pmatrix} = (e_1 \quad e_2 \quad e_3\quad e_4 \quad e_5 \quad e_6)\\$

三、应力应变法求解弹性常数

Lammps例子中的文件有 in.elastic.lmp、init.mod、potential.mod、displace.mod 以及对应的势文件。下面挑点有用的讲。主要是跟着 in.elastic.lmp 的流程理一下。

（1） init.mod

应变量up（最好测试一下微应变与最后得到的弹性常数的无关性，换几个值试一下）

单位转换不多说

定义盒子和填充原子可以用Lammps的建模方法，也可以用read_data，但是对于data文件里面是正交盒子的，需要在data文件里面添加 0 0 0 xy xz yz 的tilt倾角值，然后在read_data前面加个命令 box tilt large 不然在后面施加切应变时Lammps无法对你的正交盒子操作。

Lammps的定义盒子的方式挺好，后面可以讲讲。

（2） potential.mod

这个没啥好说的（不知道反应力场行不行，没试过）

（3） in.elastic.lmp

主流程，用来不断调用 displace.mod

先让整体静水压为零，然后弛豫，记录平衡体系的压强和盒子尺寸；然后定义计算弹性常数的variable（公式里面是去掉残余压强项，压强/应变*压强和弹性常数的系数）。

开始施加应变模式，此处直接转到下面一小段，后面没啥好说的了。

（4）displace.mod

dir决定了想到得到的弹性常数的列或行（对称性，故无区别），根据施加的是正应变还是切应变，确定了对应的len0，在前面定义的弹性常数variable中需要，如下：

if "${dir} == 1" then &
   "variable len0 equal ${lx0}" 
if "${dir} == 2" then &
   "variable len0 equal ${ly0}" 
if "${dir} == 3" then &
   "variable len0 equal ${lz0}" 
if "${dir} == 4" then &
   "variable len0 equal ${lz0}" 
if "${dir} == 5" then &
   "variable len0 equal ${lz0}" 
if "${dir} == 6" then &
   "variable len0 equal ${ly0}"

由于Lammps定义盒子是固定x方向（在我看来是这样的），然后用 xy、xz、yz 三个倾角来定义另外两个方向与x方向的关系，很直观。所以在施加切应变的时候，以dir=4为例，对应原长是lz0（x方向不动，z方向转动）。

重新读取平衡状态，定义potential，然后分别做正、负应变（正负重复实验减少误差）

最关键的是下面对施加的应变模式的理解，使用了change_box命令，改变盒子尺寸，然后remap原子坐标，然后弛豫获取压强计算弹性常数，思路就是这样。

对于deformation那段命令的理解：

variable delta   equal -${up}*${len0}
variable deltaxy equal -${up}*xy
variable deltaxz equal -${up}*xz
variable deltayz equal -${up}*yz
if "${dir} == 1" then &
   "change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box"
if "${dir} == 2" then &
   "change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box"
if "${dir} == 3" then &
   "change_box all z delta 0 ${delta} remap units box"
if "${dir} == 4" then &
   "change_box all yz delta ${delta} remap units box"
if "${dir} == 5" then &
   "change_box all xz delta ${delta} remap units box"
if "${dir} == 6" then &
   "change_box all xy delta ${delta} remap units box"

所有的delta就是施加的应变模式，这里是把每次取一个应变，如x方向：

$\begin{bmatrix}\sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_4 \\ \sigma_5 \\ \sigma_6 \end{bmatrix}=\begin{bmatrix} c_{11} & c_{12} & c_{13} & c_{14} & c_{15} & c_{16}\\ c_{21} & c_{22} & c_{23} & c_{24} & c_{25} & c_{26}\\ c_{31} & c_{32} & c_{33} & c_{34} & c_{35} & c_{36}\\ c_{41} & c_{42} & c_{43} & c_{44} & c_{45} & c_{46}\\ c_{51} & c_{52} & c_{53} & c_{54} & c_{55} & c_{56}\\ c_{61} & c_{62} & c_{63} & c_{64} & c_{65} & c_{66}\\ \end{bmatrix} \begin{bmatrix} \delta \\ 0 \\ 0\\ 0 \\0 \\0 \end{bmatrix}=\begin{bmatrix}c_{11}\\c_{21}\\c_{31}\\c_{41}\\c_{51}\\c_{61}\end{bmatrix}\delta\\$

应力应变比值就是对应的弹性常数，然后正负各1次，求平均值，6个应变量各取1次，共12次计算就得到全部的弹性常数。当然知道晶体的对称性也可以稍微改下直接求需要的元素，不过没啥意义哈。

之前困扰我的是x和y方面施加应变的时候，为啥有个 deltaxy 这些。后来请教了人特别好的戴老师，思考后才发现是因为考虑到晶胞如果不是正交的，那正应变也会导致切应变，就是得把正应变分解下。

标签：box,应变,Lammps,equal,转摘,delta,variable,常数,dir
来源： https://www.cnblogs.com/sysu/p/16485820.html