其他分享
首页 > 其他分享> > Lammps例子:弹性常数的计算(转摘)

Lammps例子:弹性常数的计算(转摘)

作者:互联网

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

 

一、广义胡克定律

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

[公式]

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

二、Voigt标记

[公式]

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

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方向:

[公式]

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

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

 

 

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