文章目录
1.计算平衡晶格常数
- 晶格是晶体结构的数学表示，晶格中的每个格点代表一个基元
- 平衡晶格常数对应的晶格具有最小的结合能
- 找到结合能曲线的最低点，对应的就是最小的结合能 E ( a 0 ) E(a_{0}) E(a0) 和平衡晶格常数 a 0 a_{0} a0
1.Ar结合能的单点计算
-
units lj # LJ单位制，无量纲 boundary p p p atom_style atomic lattice fcc 1.00 # LJ单位制下为约化密度 region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box mass 1 1.0 pair_style lj/cut 4.0 pair_coeff 1 1 1.0 1.0 4.0 variable P equal pe # 定义变量 variable L equal (count(all)/1.00)^(1/3) # 变量L为晶格常数 run 0 print "Lattice constant:$L" print "Cohensive Energy of Ar:$P" - Ar 的LJ势参数： σ = 0.3405 n m ε = 1.654 × 1 0 − 21 J = 0.010324 e V \sigma = 0.3405 \; nm\quad\quad\varepsilon =1.654×10^{-21}\; J=0.010324\; eV σ=0.3405nmε=1.654×10−21J=0.010324eV Anikeenko, A. V., G. G. Malenkov, and Yu I. Naberukhin. “Visualization of vortex movements in a molecular dynamics model of liquid argon.” Doklady Physical Chemistry. Vol. 472. Pleiades Publishing, 2017.
- 加入循环控制
units lj boundary p p p atom_style atomic label loop_i # 用指定ID来标记脚本中该命令所在的这一行 variable i loop 50 variable x equal 1.02+0.002*$i lattice fcc $x region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box mass 1 1.0 pair_style lj/cut 4.0 pair_coeff 1 1 1.0 1.0 4.0 variable P equal pe variable L equal (count(all)/1.00)^(1/3) run 0 print "Cohensive Energy of Ar a = $L E = $P" clear # 删除所有原子，将所有设置恢复到默认值，并释放由lammps分配的所有内存 next i # 将值列表中的下一个值赋给变量 jump SELF loop_i # 跳转到指定脚本的指定位置 # file = SELF时，将重新打开当前的输入脚本 - Python进行结果分析
import numpy as np import pandas as pd from matplotlib import pyplot as plt # 读取数据 (之前存储的csv文件) filename = "Ar_fcc.csv" DataFrame = pd.read_csv(filename) lat = np.array(DataFrame['lat']) ene = np.array(DataFrame['ene']) # 用二次函数拟合 a,b,c = np.polyfit(lat,ene,2) mesh = np.linspace(np.min(lat),np.max(lat),1000) opt_lat = -b/(2*a) Emin = np.polyval([a,b,c],opt_lat) print("a,b,c =",a,b,c) print("Emin =",Emin,"opt_lat =",opt_lat) # 数据可视化 plt.scatter(lat,ene,10) plt.plot(mesh,a*mesh**2+b*mesh+c) plt.xlabel('Lattice parameter') plt.ylabel('Energy') plt.show() # 模拟结果 a0 = 1.54 σ - 与理论值对比
FCC结构 Ar 气体的理论平衡晶格常数： ∂ u ∂ r = − 2 ε [ n A n ( σ n r n + 1 ) − m A m ( σ m r m + 1 ) ] = 0 \frac{\partial u}{\partial r} =-2\varepsilon [nA_{n}(\frac{\sigma ^{n}}{r^{n+1}})-mA_{m}(\frac{\sigma ^{m}}{r^{m+1}})]=0 ∂r∂u=−2ε[nAn(rn+1σn)−mAm(rm+1σm)]=0 平衡距离： r 0 = ( 2 A 12 A 6 ) 1 / 6 σ = 1.09 σ r_{0}=(\frac{2A_{12}}{A_{6}})^{1/6}\sigma=1.09\sigma r0=(A62A12)1/6σ=1.09σ 平衡晶格常数： a = 2 r 0 = 1.54 σ a=\sqrt{2}r_{0}=1.54\sigma a=2r0=1.54σ
2.金属Al的平衡晶格常数
- 晶格常数循环控制法：循环产生一系列连续变化的晶格常数，计算各值下的能量，最低值对应的晶格常数即为平衡晶格常数
units metal boundary p p p atom_style atomic variable i loop 30 variable x equal 3.90+0.01*$i lattice fcc $x region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box pair_style eam/fs # eam/fs势进行计算 pair_coeff * * Al_mm.eam.fs Al variable n eqaul count(all) variable P equal pe/$n run 0 print "Cohensive Energy of Al a = $x E = $P" clear next i jump SELF - 晶格结构弛豫法：使用共轭梯度法优化晶体结构
units metal dimension 3 boundary p p p atom_style atomic lattice fcc 4.00 region box block 0 1 0 1 0 1 units lattice create_box 1 box create_atoms 1 box pair_style eam/fs # eam/fs势函数 pair_coeff * * Al_mm.eam.fs Al neighbor 2.0 bin # 近邻列表 neigh_modify delay 10 check yes compute eng all pe/atom # 定义计算和输出 compute eatoms all reduce sum c_eng dump 1 all atom 1 relax.lammpstrj # 能量最小化 reset_timestep 0 fix 1 all box/relax iso 0.0 vmax 0.001 # 外部压强为0，体积变化最大为0.1% thermo 10 thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms min_style cg # 共轭梯度算法结构优化，最速下降法(sd) minimize 1e-25 1e-25 5000 10000 variable natoms equal count(all) variable teng equal "c_eatoms" variable a equal lx variable ecoh equal "v_teng/v_natoms" print "Total energy (eV)= ${teng};" print "Number of atoms = ${natoms};" print "Lattice constant (Angstroms) = ${a};" print "Cohensive energy (eV/atom) = ${ecoh};" - 如何预测其他元素的最稳定构型？
列举常见的晶格类型，通过改变晶格常数获得该晶体类型下的最低能量，比较不同类型晶体的最低能量即为最稳定构型 - 晶格常数的拟合区间不宜过小，过小会放大截断误差，影响计算精度
- 区间也不宜过大，过大会有非谐效应的出现，使实际曲线偏离二次型较大
2.计算体积模量
- 体积模量：是衡量物质可压缩性的指标。对于体积模量Bulk Modulus，有时也称体变模量。假设，在 P 0 P_{0} P0 的情况下体积为 V 0 V_{0} V0，若压强变化 d P dP dP，则体积变化 d V dV dV
- 体积模量 B 为： B = − V ( ∂ P ∂ V ) T B=-V(\frac{\partial P}{\partial V} )_{T} B=−V(∂V∂P)T
- 在 T=0 的情况下，压强 P 为 P = − ∂ U ∂ V P=-\frac{\partial U}{\partial V} P=−∂V∂U 其中 U 为体系的总能 B = − V ( ∂ P ∂ V ) T = V ∂ 2 U ∂ V 2 B=-V(\frac{\partial P}{\partial V} )_{T}=V\frac{\partial^{2}U}{\partial V^{2}} B=−V(∂V∂P)T=V∂V2∂2U 定义单原子的能量和体积 u = U N v = V N u=\frac{U}{N}\quad\quad v=\frac{V}{N} u=NUv=NV B = v ∂ 2 u ∂ v 2 B=v\frac{\partial ^{2}u}{\partial v^{2}} B=v∂v2∂2u
1.LJ体系的体积模量
u = 2 ε [ ∑ i ≠ j ( 1 M i j ) 12 ( σ r ) 12 − ∑ i ≠ j ( 1 M i j ) 6 ( σ r ) 6 ] = 2 ε [ A 12 ( σ r ) 12 − A 6 ( σ r ) 6 ] u=2\varepsilon[\sum_{i\ne j}(\frac{1}{M_{ij}})^{12}(\frac{\sigma}{r})^{12}-\sum_{i\ne j}(\frac{1}{M_{ij}})^{6}(\frac{\sigma}{r})^{6}] = 2\varepsilon [A_{12}(\frac{\sigma}{r})^{12}-A_{6}(\frac{\sigma}{r})^{6}] u=2ε[i=j∑(Mij1)12(rσ)12−i=j∑(Mij1)6(rσ)6]=2ε[A12(rσ)12−A6(rσ)6] A 12 = ∑ j 1 M i j 12 = 12.13 A 6 = ∑ j 1 M i j 6 = 14.45 A_{12}=\sum_{j}\frac{1}{M_{ij}^{12}}=12.13\quad\quad A_{6}=\sum_{j}\frac{1}{M_{ij}^{6}}=14.45 A12=j∑Mij121=12.13A6=j∑Mij61=14.45 根据链式求导法则： ∂ u ∂ r = ∂ u ∂ v ⋅ ∂ v ∂ r \frac{\partial u}{\partial r} =\frac{\partial u}{\partial v} \cdot \frac{\partial v}{\partial r} ∂r∂u=∂v∂u⋅∂r∂v 对于 FCC 晶体： v = a 3 4 = r 3 2 ( ∂ v ∂ r ) 2 = 9 r 4 2 v=\frac{a^{3}}{4}=\frac{r^{3}}{\sqrt{2}}\quad\quad (\frac{\partial v}{\partial r} )^{2}=\frac{9r^{4}}{2} v=4a3=2r3(∂r∂v)2=29r4 体积模量 B： B = v ∂ 2 u ∂ v 2 = v ( ∂ v ∂ r ) 2 ⋅ ∂ 2 u ∂ r 2 = 2 9 r ⋅ ∂ 2 u ∂ r 2 ∣ r = r 0 B=v\frac{\partial^{2} u}{\partial v^{2}} =\frac{v}{(\frac{\partial v}{\partial r} )^{2}}\cdot \frac{\partial^{2} u}{\partial r^{2}}=\frac{\sqrt{2}}{9r}\cdot\frac{\partial^{2} u}{\partial r^{2}}|r=r_{0} B=v∂v2∂2u=(∂r∂v)2v⋅∂r2∂2u=9r2⋅∂r2∂2u∣r=r0
- LJ体积模量的lammps脚本
结果误差来源：E-V 数据存在一定的非谐效应units lj boundary p p p atom_style atomic label LOOP variable i loop 40 # 晶格常数控制体积变化 variable x equal 1.05+0.002*$i lattice fcc $x region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box pair_style lj/cut 8.50 # 使用lj/cut势函数 pair_coeff * * 1.0 1.0 mass 1 1.0 thermo_style custom etotal variable P equal pe # 定义计算变量，P为单原子势能 variable v eqaul (1.0/$x) # 单原子体积 dump 1 all custom 10 trj.lammpstrj id type x y z vx vy vz run 0 print "fcc lattice rho = $v E = $P" clear next i jump SELF LOOP - 拟合物态方程求体弹模量：
拟合Birch-Murnaghan方程： E ( V ) = E 0 + B 0 V B 0 ′ ( ( V 0 / V ) B 0 ′ B 0 ′ − 1 + 1 ) − B 0 V 0 B 0 ′ − 1 E(V)=E_{0}+\frac{B_{0}V}{B'_{0}}(\frac{(V_{0}/V)^{B'_{0}}}{B'_{0}-1}+1)-\frac{B_{0}V_{0}}{B'_{0}-1} E(V)=E0+B0′B0V(B0′−1(V0/V)B0′+1)−B0′−1B0V0 待拟合参数： E 0 E_{0} E0：平衡状态结合能， V 0 V_{0} V0：平衡状态体积， B 0 B_{0} B0：平衡状态体积模量， B 0 ′ B'_{0} B0′：初始值取 3.5 - Python代码拟合Birch-Murnaghan方程
import pandas as pd import numpy as np # 读取数据 (之前存储的csv文件): df = pd.read_csv('data.csv') v = np.array(df['v']) e = np.array(df['e']) # 先做一次多项式拟合，根据多项式拟合的结果设置初始猜测解x0 a,b,c = np.polyfit(v,e,2) v0 = -b/(2*a) e0 = a*v0**2 + b*v0 + c b0 = 2*a*v0 bP = 3.5 x0 = [e0,b0,bP,v0] def Murnaghan(parameters.vol): E0 = parameters[0] B0 = parameters[1] BP = parameters[2] V0 = parameters[3] E = E0 + B0*vol/BP*(((V0/vol)**BP)/(BP-1)+1) - V0*B0/(BP-1.) return E # 定义误差 (优化目标) def residual(pars,y,x): err = y - Murnaghan(pars,x) return err # 进行拟合 from scipy.optimize import leastsq murnpars,ier = leastsq(residual,x0,args=(e,v)) # 输出结果 print('Bulk Modulus:' + str(murnpars[1])) print('lattice constant:',murnpars[3]**(1/3)) # 拟合效果可视化 from matplotlib import pyplot as plt v_mesh = np.linspace(np.min(v),np.max(v),1000) plt.scatter(v,e,10) plt.plot(v_mesh,Murnaghan(murnpars,v_mesh)) plt.ylabel('Energy') plt.xlabel('Volume') plt.show() - Birch-Murnaghan 方程因为更接近于物理实际，所有得到的结果与解析解最为接近；体积模量的拟合，应优先使用Brich-Murnaghan方程
2.计算Cu的体积模量
- Cu 的LJ势参数： σ = 0.2338 n m ε = 65.58 × 1 0 − 21 J = 0.40933 e V \sigma = 0.2338\; nm\quad\quad\varepsilon =65.58 × 10^{-21}\;J=0.40933\;eV σ=0.2338nmε=65.58×10−21J=0.40933eV Lv, Jizu, et al. “The molecular dynamic simulation on impact and friction characters of nanofluids with many nanoparticles system.” Nanoscale research letters 6 (2011): 1-8.
- 根据LJ势参数计算
原子间平衡距离： r 0 = 1.09 σ = 2.55 A ˚ r_{0}=1.09\sigma =2.55\;\AA r0=1.09σ=2.55A˚ 结合能： u 0 = 2 ε [ A 12 ( σ r 0 ) 12 − A 6 ( σ r 0 ) 6 ] = − 8.6069 ε = − 3.52 e V u_{0}=2\varepsilon[A_{12}(\frac{\sigma}{r_{0}})^{12}-A_{6}(\frac{\sigma}{r_{0}})^{6}]=-8.6069\varepsilon =-3.52\;eV u0=2ε[A12(r0σ)12−A6(r0σ)6]=−8.6069ε=−3.52eV 实验中 Cu-Cu 键长：2.56 Å，Cu 的结合能：-3.49 eV - LJ势计算Cu的体积模量
units metal boundary p p p atom_style atomic # 设置loop variable i loop 40 variable x equal 3.40+0.01*$i lattice fcc $x region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box mass 1 64 # 使用Cu的LJ势参数 pair_style lj/cut 10.0 pair_coeff 1 1 0.40933 2.338 variable v equal ($x)^3 variable n equal count(all) variable P equal pe # 运行 run 0 print "Cohensive Energy of Cu v = $v x = $x E = $P" clear next i jump SELF
计算值与实验值差距很大，原因：LJ 势函数不能准确描述金属原子相互作用
- EAM势计算Cu体积模量
units metal boundary p p p atom_style atomic # 设置loop variable i loop 40 variable x equal 3.40+0.01*$i lattice fcc $x region box block 0 1 0 1 0 1 create_box 1 box create_atoms 1 box mass 1 64 # 使用EAM势 pair_style eam pair_coeff * * Cu_u3.eam variable v equal ($x)^3 variable n equal count(all) variable P equal pe # 运行 run 0 print "Cohensive Energy of Cu v = $v x = $x E = $P" clear next i jump SELF
计算体积模量：136.93 GPa，实验值范围：130 ~ 145 GPa
3.材料拉伸压缩模拟
- 应力和应变之间的关系称为该材料的应力-应变曲线
应力： σ = F A 0 \sigma=\frac{F}{A_{0}} σ=A0F 应变： ε = L − L 0 L 0 = Δ L L 0 \varepsilon =\frac{L-L_{0}}{L_{0}}=\frac{\Delta L}{L_{0}} ε=L0L−L0=L0ΔL
1.拉伸压缩基础命令
- 建模命令
read_data、lattice、region、create_atoms - 边界条件设置
deform拉伸，盒子等比例拉伸：拉伸方向为p
velocity拉伸，固定两端，设置速度场：拉伸方向为s - 弛豫过程
拉伸方向尽量在npt(不适用s边界条件)下弛豫，消除内应力 - 应变公式
variable tmp equal "ly" #定义tmp变量为ly variable L0 equal ${tmp} #将第一次tmp取值存入L0 variable strain equal "(ly-v_L0)/v_L0" #定义应变变量 variable p1 equal "v_strain" - 应变速率
variable srate equal 1.0e10 variable srate1 equal "v_srate/1.0e12" #速率单位转换 - 应力计算
variable p2 equal "-pyy/10000" #单位转换 # 或 compute peratom all stress/atom NULL #计算所有原子的单原子应力，定义为变量peratom compute p all reduce sum c_peratom[2] #将所有原子的Y方向的应力加和，定义为p variable p2 equal c_p/vol #平均应力=应力和/体积 - 应力与应变数据保存到文件
fix def1 all print 100 "${p1}${p2}" file strain_stress.txt screen no #每隔100步输出所有原子的应力应变，保存为指定文件，屏幕不输出
2.单层石墨烯拉伸
- lammps脚本
units real # real单位制 timestep 1 variable fpress equal 0.000101325 # fpress和fenergy用于单位转换，atm -> GPa variable fenergy equal 0.043 # Kcal/mole -> eV # 单层石墨烯 dimension 3 boundary p p f atom_style charge # 指定原子属性包括种类、坐标和电荷 read_data grapheneLarge.dat # 导入石墨烯模型 pair_style reax/c NULL checkqeq no # 使用反应力场势函数，lammps需编译USER-REAXC包 pair_coeff * * ffield.reax C # 温度设置 variable temp equal 300 variable seed equal 1356 velocity all create ${temp} ${seed} dist uniform # 正态分布产生温度为300K的初始温度 # 平衡 fix fnpt all npt temp ${temp} ${temp} 10 x 0 0 500 y 0 0 500 # 所有原子温度为${temp}，温度阻尼为10，控制x和y方向上压力为0，压力阻尼为500，z方向不施加压力控制 thermo 100 run 1000 # 在NPT系综下弛豫1000步 unfix fnpt # 输出，假设石墨烯厚度为 3.35 A variable tmp equal lx variable lx0 equal ${tmp} # lx0 初始x方向长度 variable tmp equal ly # 对于equal的变量，每次引用该变量时，其值都会被重新计算 variable ly0 equal ${tmp} # ly0 初始y方向长度 variable Eavg equal etotal/atoms*${fenergy} # Eavg 原子平均能量 variable pe equal pe/atoms*${fenergy} variable ke equal ke/atoms*${fenergy} variable strainx equal (lx-${lx0})/${lx0} variable strainy equal (ly-${ly0})/${ly0} variable stressx equal -pxx*(lz/3.35)*${fpress} variable stressy equal -pyy*(lz/3.35)*${fpress} thermo_style custom step time temp etotal press v_Eavg v_pe v_ke v_strainx v_stressx v_strainy v_stressy # 变形 fix boxdeform all deform 1 x scale 2 remap x # fix deform命令改变盒子的体积或形状，每跑1步变形1次，X方向最终形变为初始的2倍长度 ，原子位置随box同比例形变 fix fnpt all npt temp ${temp} ${temp} 10 y 0 0 500 thermo 100 run 10000 - 模拟结果：极限拉伸强度为165 GPa，线弹性区最大至应变20%，拉伸到应变37%发生断裂，杨氏模量： E = σ / ε = 0.814 T P a E=\sigma/\varepsilon =0.814\;TPa E=σ/ε=0.814TPa
- 实验测量杨氏模量：0.9 ~ 1.1 TPa
Lee, Changgu, et al. “Measurement of the elastic properties and intrinsic strength of monolayer graphene.” science 321.5887 (2008): 385-388.
3.Fe三向拉伸
- lammps脚本
# Fe三向拉伸 #------------模型设置------------- units metal #单位制 dimension 3 #3维 boundary p p p #周期性边界条件 atom_style atomic #原子类型 neighbor 2 bin #邻居列表 timestep 0.005 #时间步长 #-----------建立模型-------------- lattice bcc 2.8665 #晶格 region box block 0 20 0 20 0 20 #建立box区域 create_box 1 box #设置原子类型数 create_atoms 1 box #box区域生成类型1原子 #----------力场设置--------------- pair_style eam/fs #eam势 pair_coeff * * Fe.eam.fs Fe #----------能量最小化------------ thermo 100 #每隔100步输出 thermo_style custom step temp press pe #输出步数、温度、压力、势能 dump 1 all custom 100 Fe_mini.xyz id type x y z #每隔100步存储原子ID 类型 xyz坐标 min_style cg #cg法进行能量最小化，共轭梯度法 minimize 1e-5 1e-15 5000 5000 #能量容差，力容差，最大迭代次数，最大评估次数 undump 1 #取消转存 #----------速度初始化------------- velocity all create 300 876848 dist gaussian #初始温度300K，高斯分布 #----------npt弛豫----------------- fix 1 all npt temp 300 300 0.5 iso 0 0 5 #恒温恒压，300K dump 1 all custom 1000 Fe_npt.xyz id type x y z #每隔1000步存储原子ID 类型 xyz坐标 run 10000 #运行步数 unfix 1 undump 1 #取消命令 #-------------应力应变设置--------------- variable tmpx equal "lx" #x方向 variable Lx0 equal ${tmpx} variable strainx equal "(lx - v_Lx0)/v_Lx0" variable stressx equal "-pxx/10000" variable tmpy equal "ly" #y方向 variable Ly0 equal ${tmpy} variable strainy equal "(ly - v_Ly0)/v_Ly0" variable stressy equal "-pyy/10000" variable tmpz equal "lz" #z方向 variable Lz0 equal ${tmpz} variable strainz equal "(lz - v_Lz0)/v_Lz0" variable stressz equal "-pzz/10000" #----------------应力应变保存文件------------------------ fix def3 all print 100 "${strainx} ${stressx} ${strainy} ${stressy} ${strainz} ${stressz}" & file stress_strain.dat screen no #每100步保存文件，不输出屏幕 #----------------三向拉伸---------------------------- fix 1 all nve #nve系综 fix 2 all temp/rescale 10 300 300 10 1.0 #速度标定法，每10步重新缩放，温度差超过10的时候进行缩放，分数=1.0表示温度超出后重置为期望值 fix 3 all deform 100 x erate 0.005 y erate 0.005 z erate 0.005 remap x units box #盒子单位，remap x将原子坐标重新映射到盒子 dump 1 all custom 100 Fe_tension.xyz id type x y z #每100步存储 run 50000
4.高熵合金FeNiCrCoCu拉伸
- lammps脚本
#初始参数设置 units metal #单位 dimension 3 #维度 boundary p p p #边界条件 atom_style atomic #原子类型 timestep 0.001 #时间步长 neighbor 2.0 bin #邻居列表 #模型建立 lattice fcc 3.6149 #Cu晶格 region box block 0 10 0 20 0 10 #定义box区域 create_box 5 box #原子类型数 create_atoms 1 box #生成类型1原子 #按比例将类型1原子替换为其他原子 set type 1 type/ratio 2 0.2 87393 set type 1 type/ratio 3 0.5 87393 set type 1 type/ratio 4 0.5 87393 set type 3 type/ratio 5 0.5 87393 #设置摩尔质量 mass 1 55.845 #Fe mass 2 58.6934 #Ni mass 3 51.9961 #Cr mass 4 58.933194 #Co mass 5 63.546 #Cu #保存模型data文件 write_data al1.xyz #设置原子间相互作用势 pair_style eam/alloy #eam势 pair_coeff * * FeNiCrCoCu-heafixed.setfl Fe Ni Cr Co Cu #输出热力学信息 thermo 100 thermo_style custom step temp pe ke press #每100步输出步数、温度、动能、势能、压力 #能量最小化 dump 1 all atom 10 mini.xyz #每10步将原子相关信息保存为文件 min_style cg #使用cg方法进行能量最小化 minimize 1e-15 1e-15 5000 5000 #能量容差、力容差、最大迭代次数、最大评估次数 undump 1 #取消 reset_timestep 0 #设置当前时间步为第0步 #温度初始化 velocity all create 300 88989 #初始温度300K #npt系综下弛豫 dump 1 all custom 1000 npt.xyz id type x y z #每1000步保存原子信息 fix 1 all npt temp 300 300 0.1 iso 0 0 1 #恒温恒压 run 50000 undump 1 unfix 1 #取消 reset_timestep 0 #应力应变输出 variable tmp equal "ly" #y方向 variable L0 equal ${tmp} variable strain equal "(ly-v_L0)/v_L0" variable p1 equal "v_strain" variable p2 equal "-pyy/10000" fix def1 all print 100 "${p1} ${p2}" file strain_stress.txt screen no #每100步保存应力应变，不输出屏幕 #定义拉伸应变速率 variable srate equal 1.0e10 variable srate1 equal "v_srate / 1.0e12" #保存拉伸轨迹文件 shell "mkdir tension" #执行shell命令，新建tension文件夹 dump 1 all custom 100 tension/dump.tensile_*.xyz type x y z #‘*’表示每100步保存一次 #nvt系综下使用deform单轴拉伸 fix 1 all nvt temp 298 298 0.1 fix 2 all deform 10 y erate ${srate1} remap x units box run 25000
5.铝和石墨烯复合物压缩变形
- lammps脚本
#模型基本设置 units metal #单位 dimension 3 #维度 boundary p p p #边界条件 timestep 0.001 #时间步长 neighbor 0.3 bin #邻居列表 neigh_modify delay 0 #box设置 region box block 0 50 0 50 0 40 units box #设置box区域，使用盒子单位 create_box 2 box #box区域有两种原子 #生成Al原子 lattice fcc 4.0495 #Al晶格 region al block 0 50 0 50 0 40 units box #设置al区域，使用盒子单位 create_atoms 1 region al #在al区域生成类型1原子 #删除中间一层Al原子 region del block 0 50 0 50 18 23 units box #设置待删除区域，使用盒子单位 delete_atoms region del compress yes #compress yes删除区域后，对原子重新排序 #自定义石墨烯晶格 lattice custom 2.4768 a1 1.0 0.0 0.0 a2 0.0 1.732 0.0 a3 0.0 0.0 1.3727 & basis 0.0 0.33333 0.0 & basis 0.0 0.66667 0.0 & basis 0.5 0.16667 0.0 & basis 0.5 0.83333 0.0 & #在中间层插入石墨烯原子 region graphene block 0 50 0 50 18 23 units box #设置石墨烯区域，使用盒子单位 create_atoms 2 region graphene #在石墨烯区域生成类型2原子 #设置原子质量 mass 1 26.9815385 #Al mass 2 12.011 #C #设置力场参数 pair_style hybrid lj/cut 10 eam/fs airebo 3.0 #使用混合势 pair_coeff 1 2 lj/cut 0.0234943 3.1482 #Al-C之间 pair_coeff * * eam/fs Al1.eam.fs Al NULL #Al-Al之间 pair_coeff * * airebo CH.airebo NULL C #C-C之间 #能量最小化 dump 1 all atom 1 1_mini.lammpstrj #每隔1步保存原子相关信息 min_style cg #以cg法进行能量最小化 minimize 1e-15 1e-15 10000 10000 #能量容差、力容差、最大迭代次数、最大评估次数 undump 1 #取消输出 #温度初始化 velocity all create 300 8877423 #初始温度300K #npt系综下进行弛豫 fix 1 all npt temp 300 300 0.01 iso 0 0 1 #恒温恒压 dump 1 all atom 1000 2_npt.lammpstrj #每1000步将原子相关信息保存为文件 run 50000 unfix 1 #取消 undump 1 reset_timestep 0 #设置当前步为第0步 #应力应变设置 variable tmp equal "lz" #z方向 variable L0 equal ${tmp} variable strain equal "(v_L0-lz)/v_L0" #应变 variable stressz equal "pzz/10000" #应力 variable srate equal 1.0e10 #应变速率 variable srate1 equal "-1*v_srate/1.0e12" #应力应变保存文件 fix def3 all print 100 "${strain} ${stressz}" file stress_strain.dat screen no #每100步保存，不输出屏幕 #热力学输出 thermo 100 thermo_style custom step temp lx ly lz vol pxx pyy pzz #每100步输出步数、温度、xyz长度、体积、xyz方向压力 #压缩设置 fix 1 all npt temp 300 300 0.01 x 0 0 0.1 y 0 0 0.1 #恒温，xy方向压缩 fix 2 all deform 1 z erate ${srate1} units box remap x #使用盒子单位，z方向压缩，remap x映射到盒子 #输出轨迹文件 dump 1 all atom 100 3_deform.lammpstrj #每100步保存原子相关信息 run 30000
4.其他力学性质计算
1.弹性常数
- 弹性常数(elastic constants)表征材料弹性的物理量，弹性形变时应力与应变满足胡克定律。联系各向异性介质中应力和应变关系的广义弹性张量共有81个分量，但这81个分量并不独立，可以根据晶体的对称性简化为21个独立分量 σ i j = C i j k l ε k l \sigma_{ij}=C_{ijkl}\varepsilon _{kl} σij=Cijklεkl c i j m n = c j i m n = c i j n m = c j i n m = c m n i j = c n m i j = c m n j i = c n m j i c_{ijmn}=c_{jimn}=c_{ijnm}=c_{jinm}=c_{mnij}=c_{nmij}=c_{mnji}=c_{nmji} cijmn=cjimn=cijnm=cjinm=cmnij=cnmij=cmnji=cnmji
- 以FCC的晶体Cu为例，根据对称性和独立弹性常数个数的关系可知，具有立方对称性的材料有3个独立弹性常数：
C
=
(
C
11
C
12
C
12
0
0
0
C
12
C
11
C
12
0
0
0
C
12
C
12
C
11
0
0
0
0
0
0
C
44
0
0
0
0
0
0
C
44
0
0
0
0
0
0
C
44
)
C=\begin{pmatrix}C_{11} &C_{12}&C_{12}&0&0&0 \\ C_{12}&C_{11}&C_{12}&0&0&0 \\ C_{12}&C_{12}&C_{11}&0&0&0 \\ 0&0&0&C_{44}&0&0 \\ 0&0&0&0&C_{44}&0 \\ 0&0&0&0&0&C_{44} \end{pmatrix}
C=
C11C12C12000C12C11C12000C12C12C11000000C44000000C44000000C44
体系的弹性内能： E e l a s / V = 1 2 ∑ i j C i j ε i ε j E^{elas}/V=\frac{1}{2}\sum_{ij}C_{ij}\varepsilon _{i}\varepsilon _{j} Eelas/V=21ij∑Cijεiεj 可构造相应的弹性形变来计算弹性常数 - 获得弹性常数，构造出刚度张量，就可以计算材料的各种弹性性质
刚度张量： C = ( C 11 C 12 C 12 0 0 0 C 12 C 11 C 12 0 0 0 C 12 C 12 C 11 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 ) C=\begin{pmatrix}C_{11} &C_{12}&C_{12}&0&0&0 \\ C_{12}&C_{11}&C_{12}&0&0&0 \\ C_{12}&C_{12}&C_{11}&0&0&0 \\ 0&0&0&C_{44}&0&0 \\ 0&0&0&0&C_{44}&0 \\ 0&0&0&0&0&C_{44} \end{pmatrix} C= C11C12C12000C12C11C12000C12C12C11000000C44000000C44000000C44 柔度张量 [ S ] = [ C ] − 1 [S]=[C]^{-1} [S]=[C]−1
2.杨氏模量和泊松比
- 弹性材料承受正向应力时会产生正向应变，定义正向应力与正向应变的比值为材料的杨氏模量(Young’s modulus)： Y = σ ε = F / A 0 Δ L / L 0 = F L 0 A 0 Δ L Y=\frac{\sigma}{\varepsilon}=\frac{F/A_{0}}{\Delta L/L_{0}}=\frac{FL_{0}}{A_{0}\Delta L} Y=εσ=ΔL/L0F/A0=A0ΔLFL0 弹性材料承受单轴拉伸时，唯一的应力为 σ x x \sigma_{xx} σxx ε x x = S 11 σ x x \varepsilon_{xx}=S_{11}\sigma_{xx} εxx=S11σxx Y = σ x x ε x x = 1 S 11 = C 11 2 + C 11 C 12 − 2 C 12 2 C 11 + C 12 Y=\frac{\sigma_{xx}}{\varepsilon_{xx}}=\frac{1}{S_{11}}=\frac{C_{11}^{2}+C_{11}C_{12}-2C_{12}^{2}}{C_{11}+C_{12}} Y=εxxσxx=S111=C11+C12C112+C11C12−2C122
- 材料受拉伸或压缩力时，其横向变形量与纵向变形量的比值称为泊松比，是一个无量纲的物理量 v = − ε l a t e r a l ε a x i a l v=-\frac{\varepsilon_{lateral}}{\varepsilon_{axial}} v=−εaxialεlateral ε x x = S 11 σ x x ε y y = S 21 σ x x \varepsilon_{xx}=S_{11}\sigma_{xx}\quad\quad\varepsilon_{yy}=S_{21}\sigma_{xx} εxx=S11σxxεyy=S21σxx v = ε y y ε x x = S 21 S 11 = − C 12 C 11 + C 12 v=\frac{\varepsilon_{yy}}{\varepsilon_{xx}}=\frac{S_{21}}{S_{11}}=\frac{-C_{12}}{C_{11}+C_{12}} v=εxxεyy=S11S21=C11+C12−C12
扫一扫
1638
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
