LAMMPS讲解67-二维粗糙度铜块对磨模拟

带有粗糙度的摩擦需要一些很小新的设置。这里举例说明，如何施加额定的法向载荷和切向载荷，局部控温和载荷监控。具体方法在in文件的注释中。这里要说明的是，磨损过程中温度会升高很快，这是正常的，因为剧烈的摩擦产生大量的热不能及时传递出去，采用较慢的摩擦速度会降低温度升高这一现象。这里采用了铜块。单质铜是比较软的，所以摩擦过程中会出现焊接现象。改用一硬一软的材料应该就会出现实际中的磨损现象。粗糙度建模分为两步首先在LAMMPS中或者其他软件建立一个金属块，然后通过python代码切除粗糙度。data文件中的模型是一个长方体的金属材料。粗糙度在z方向上生成。如果要生成上下两块则可以分别生成然后在Excel中将两块组合起来。

该方法的原理是使用Weierstrass-Mandelbrot (W-M)函数生成坐标，然后基于生成 的坐标删除多余的原子。

import math

import matplotlib.pyplot as plt

#写出原始data文件的路径和修改后data的路径

data = open("E:\\work\\0_wenku\\磨损\\1.data", 'r')

data_modify = open("E:\\work\\0_wenku\\磨损\\3.data", 'w')

line = data.readline()

line = data.readline()

line = data.readline().split()

natom = int(line[0])

line = data.readline()

line = data.readline()

line = data.readline().split()

xlo = float(line[0])

xhi = float(line[1])

line = data.readline().split()

ylo = float(line[0])

yhi = float(line[1])

line = data.readline().split()

zlo = float(line[0])

zhi = float(line[1])

line = data.readline()

line = data.readline()

line = data.readline()

line = data.readline()

line = data.readline()

line = data.readline()

line = data.readline()

line = data.readline()

data_modify.write("LAMMPS data of ice slab with rough\n")

data_modify.write("\n")

data_modify.write(str(natom) + " atoms\n")

data_modify.write("\n")

data_modify.write("1 atom types\n")

data_modify.write("\n")

data_modify.write(str(xlo) + " " + str(xhi) + " xlo xhi\n")

data_modify.write(str(ylo) + " " + str(yhi) + " ylo yhi\n")

data_modify.write(str(zlo) + " " + str(zhi) + " zlo zhi\n")

data_modify.write("\n")

data_modify.write("Atoms\n")

data_modify.write("\n")

#调整Gx，Dx和rx可获得不同形状的粗糙度

Gx=500.0

Dx=1.555

rx=1.45

x=0

jx=1

z=0

nx = []

ny = []

ix = 0.1

inter=18.075 #金属块x方向长度除以20

num = 24； #num*inter要超过金属块在x方向上的长度

while ix < num: #

    x = ix*inter

    z = 0.0

    jx=1

    ix += 1

    while jx < 100:

        z_temp = (rx**(Dx*jx-2*jx))*(math.cos(2*(math.pi)*(rx**jx)*x))

        z += z_temp

        jx += 1

    z = (Gx**(Dx-1))*z+50 #50是粗糙度的z方向上的偏移量，调整以得到合适的坐标

    nx.append(x)

    ny.append(z)

ny = list(reversed(ny))

#将粗糙度画出来看是否合理

plt.plot(nx,ny)

plt.xlim(0,361.5)

plt.ylim(0,361.5)

plt.show()  

natom_modify = 0

done = 0

while not done:

    line = data.readline()

    if line != '\n':

        pos = line.split()

        x = float(pos[2])

        i = int((x)//inter)

        z0 = ((ny[i+1] - ny[i])/(nx[i+1] - nx[i]))*(x - nx[i]) + ny[i]

        if float(pos[4]) >= z0:

            data_modify.write(line)

            natom_modify += 1

    else:

        done = 1

data.close()

data_modify.close()

#最后要用打印出的原子数替换修改后的data文件的原子数

print(str(natom_modify))

本例中的data文件就是采用上面代码生成的，以下是模拟的in文件：

dimension    3

boundary     p p s

units        metal

atom_style   atomic

pair_style   eam

read_data system.data

pair_coeff * * Cu_u3.eam

#上下两块都会分为三层，分别是固定层，温控层和自由层，通过region命令将三层分别定义为不同组

region upfix block INF INF INF INF 377 INF units box

group upfix region upfix

region downfix block INF INF INF INF INF 1.9 units box

group downfix region downfix

region upt block INF INF INF INF 357 377 units box

group upt region upt

region downt block INF INF INF INF 1.9 21.9 units box

group downt region downt

#为了后面统计受力取消下面一块金属固定层之间的相互作用

neighbor 3.0 bin

neigh_modify every 1 delay 0 check yes exclude group downfix downfix

group up type 2

group down type 1

group upfree subtract up upfix upt

group downfree subtract down downfix downt

group t union upt downt

group free union upfree downfree

timestep 0.002

thermo 100

#分别监控温控层和自由层的温度

compute myTfree free temp

compute myTt t temp

#下面一块金属的固定层的受力都来自上方，因此计算该固定层的切向力和法向力就可以得到切向摩擦力和载荷的大小

compute fx downfix reduce sum fx

compute fz downfix reduce sum fz

#监控输出温度和载荷

thermo_style custom step c_myTfree c_myTt c_fx c_fz

thermo_modify flush yes

velocity t create 300 5645354 temp myTt loop local dist gaussian rot yes

velocity free create 300 565454 temp myTfree loop local dist gaussian rot yes

dump mydump1 all atom 1000 atom1.lammpstrj

#对自由层和温控层进行弛豫

fix fxnvt1 t nvt temp 300.0 300.0 0.2

fix fxnvt2 free nvt temp 300.0 300.0 0.2

run 100000

#取消自由层上的温度控制

unfix fxnvt2

fix fxnve free nve

#结合使用一下三个fix施加法向载荷使两块金属接触

fix 1 upfix setforce 0.0 0.0 NULL

fix 2 upfix aveforce 0.0 0.0 -0.01

fix 3 upfix nve

run 100000

#施加切向载荷开始剪切

unfix 1

unfix 2

fix 1 upfix setforce NULL 0.0 NULL

fix 2 upfix aveforce 0.1 0.0 -0.01

run 1000000

结果如下：

温度变化的基本规律是正确的，但是这些尖峰是不对的。这可能是因为温控层也在移动，没有减去宏观速度就进温控的原因。载荷的变化是符合预期的，刚开始法向和切向载荷都为零。然后上下两块开始接触固体块出现弹性振荡，两个载荷也出现振荡。当开始剪切的时候切向载荷会迅速升高，因为粗糙度之间开始磨损。待粗糙度磨平后，载荷开始下降并维持在基本不变的范围内。可以看到，由于原子数没有足够大，载荷会出现明显的涨落。这是我们应该对载荷取平均来获得平均载荷用以计算摩擦系数。为了改进温度的输出，采用一下改进的in文件（未经测试，欢迎加微信baolu_yao，和我讨论）

dimension    3

boundary     p p s

units        metal

atom_style   atomic

pair_style   eam

read_data system.data

pair_coeff * * Cu_u3.eam

#上下两块都会分为三层，分别是固定层，温控层和自由层，通过region命令将三层分别定义为不同组

region upfix block INF INF INF INF 377 INF units box

group upfix region upfix

region downfix block INF INF INF INF INF 1.9 units box

group downfix region downfix

region upt block INF INF INF INF 357 377 units box

group upt region upt

region downt block INF INF INF INF 1.9 21.9 units box

group downt region downt

#为了后面统计受力取消下面一块金属固定层之间的相互作用

neighbor 3.0 bin

neigh_modify every 1 delay 0 check yes exclude group downfix downfix

group up type 2

group down type 1

group upfree subtract up upfix upt

group downfree subtract down downfix downt

group t union upt downt

group free union upfree downfree

timestep 0.002

thermo 100

#分别监控温控层和自由层的温度

compute myTupfree upfree temp/partial 0 1 0

compute myTdownfree downfree temp/partial 0 1 0

compute myTupt upt temp/partial 0 1 0

compute myTdownt downt temp/partial 0 1 0

#下面一块金属的固定层的受力都来自上方，因此计算该固定层的切向力和法向力就可以得到切向摩擦力和载荷的大小

compute fx downfix reduce sum fx

compute fz downfix reduce sum fz

#监控输出温度和载荷

thermo_style custom step c_myTupfree c_myTdownfree c_myTupt c_myTdownt c_fx c_fz

thermo_modify flush yes

velocity upt create 300 5643454 temp myTupt loop local dist gaussian rot yes

velocity downt create 300 564654354 temp myTdownt loop local dist gaussian rot yes

velocity upfree create 300 567854 temp myTupfree loop local dist gaussian rot yes

velocity downfree create 300 590454 temp myTdownfree loop local dist gaussian rot yes

dump mydump1 all atom 1000 atom1.lammpstrj

#对自由层和温控层进行弛豫

fix fxnvt1 t nvt temp 300.0 300.0 0.2

fix fxnvt2 free nvt temp 300.0 300.0 0.2

run 100000

#取消自由层上的温度控制

unfix fxnvt1

unfix fxnvt2

fix fxnvt1 upt nvt temp 300.0 300.0 0.2

fix_modify fxnvt1 temp myTupt

fix fxnvt2 downt nvt temp 300.0 300.0 0.2

fix_modify fxnvt2 temp myTdownt

fix fxnve free nve

#结合使用一下三个fix施加法向载荷使两块金属接触

fix 1 upfix setforce 0.0 0.0 NULL

fix 2 upfix aveforce 0.0 0.0 -0.01

fix 3 upfix nve

run 100000

#施加切向载荷开始剪切

unfix 1

unfix 2

fix 1 upfix setforce NULL 0.0 NULL

fix 2 upfix aveforce 0.1 0.0 -0.01

run 1000000

感谢鲍路瑶老师的分享，内容来自于鲍老师分享出来的资料

如有需要添加微信：lmp_zhushou  进入微信群，帮助他人，共建社区

获取完整版lammps讲义可以加微信lmp_zhushou或加入QQ群994359511