3D Poiseuille flow of MD by using lammps

Hi,

I am using lammps to simulate the 3D Poiseuille flow. My model is a cubic and I use the fix addforce command to give each atom a force which is along with x direction. Unlike the other 3D poiseuille flow, there are two walls used to confine the fluid zone in y and z direction respectively.

However, when my code is run about 500,000 timesteps and there is an error “ERROR: Lost atoms: original 20480 current 19873 (…/thermo.cpp:442)”. I found that some fluid atoms will pass through the solid wall by outputting the number density of fluid atom using commands “compute mychunk1 flow chunk/atom bin/1d z 1.635 0.545 units box” and “fix result1 flow ave/chunk 2 2000001 5000000 mychunk1 density/number temp bias mytemp file result1.out format %20.16g”.

My input file is shown as follows. I will be grateful if you can help me fix this problem. Looking forward for your reply.

units nano

variable T_fluid equal 300

variable T_wall equal 300

dimension 3

boundary p f f

atom_style atomic

neighbor 2.0 bin

neigh_modify delay 1 check yes

#geometry create

lattice fcc 0.545

region box block 0 20 0 16 0 16

create_box 2 box

create_atoms 1 box

mass 1 6.6335209/1e5

mass 2 6.6335209/1e5

#define groups

region 1 block 0 20 3 13 3 13

group flow region 1

group boundary subtract all flow

set group boundary type 2

#initial velocities:

velocity boundary create ${T_wall} 99999 dist gaussian mom yes

fix wall_force boundary spring/self 70e3

velocity flow create ${T_fluid} 99999 dist gaussian mom yes

#L-J potential

pair_style lj/cut 0.85125

pair_coeff 1 1 1.657 0.3405 0.85125

pair_coeff 1 2 1.657 0.3405 0.85125

pair_coeff 2 2 1.657 0.3405 0.85125

#control the temperature of wall

compute temp_wall boundary temp

fix wall_1 boundary nve

fix wall_2 boundary temp/rescale 50 {T_wall} {T_wall} 0.1 1

#definition of flow

fix definition_flow flow addforce 10.0 0.0 0.0

#control the temperature of fluid

compute mytemp flow temp/com

fix flow_1 flow nve

fix flow_2 flow temp/rescale 50 {T_fluid} {T_fluid} 0.1 1

fix_modify flow_2 temp mytemp

#output results

compute mychunk1 flow chunk/atom bin/1d z 1.635 0.545 units box # the zone of 0.5453=1.635<=z<=0.54513=7.085 is the fluid region.

compute mychunk2 flow chunk/atom bin/2d x 0.0 0.545 z 1.635 0.545 units box # the zone of 0.5453=1.635<=z<=0.54513=7.085 is the fluid region.

fix result1 flow ave/chunk 2 2000001 5000000 mychunk1 density/number temp bias mytemp file result1.out format %20.16g

fix result2 flow ave/chunk 2 2000001 5000000 mychunk2 density/number temp bias mytemp file result2.out format %20.16g

#run

thermo 20000

thermo_modify temp mytemp

timestep 1e-6

run 5000000

Best wishes,

Qiangqiang Sun

in.flow (1.79 KB)

1 Like