Dear lammps users,
I am modeling a 2D flow. There are 100,000 fluid atoms which are confined by two solid atoms wall perpendicular to y axis (each solid wall consist four layers and there are 6,000 solid atoms for per wall) flow in x direction. Firstly, all fluid and solid atoms are initialised velocities corresponding the given temperature 70K, and run this for 2 ns (my tilmestep is 1.0e-5 ns). Then all fluid atoms are added a body force (1.0e-12 N) in x direction and langevin thermostat is used to control the temperature for all fluid atoms equal to 120K. For solid wall atoms, the temperature is still 70K. The following 3 ns is used for the system to reach a steady state. After that, 10 ns is adopted for temperature sampling used ‘compute temp/chunk’ and ‘compute temp/profile’ command respectively. The result is attached. For the to different temperature sampling methods, there are around 444 fluid atoms in each bin. The temperature for bins near the solid wall deviate from 120K and reach to 132K. I will be grateful for any review. My input file is :
units nano
variable T_wall equal 70
variable T_fluid equal 120
variable spring_constant equal 70000 #70.0e3
dimension 2
boundary p f p
atom_style atomic
neighbor 1.0 bin
neigh_modify delay 1 check yes
#-------------------------model----------------------------
lattice sq 0.4
region box block 0 1200 0 68.0001 -0.25 0.25
create_box 2 box
#-------------------------region----------------------------
region r_bot_wall block INF INF INF 4.01 INF INF
region r_top_wall block INF INF 63.99 INF INF INF
region r_flow block INF INF 4.01 63.99 INF INF
region r_flow_q block INF INF 8 60 INF INF
#--------------------------create atoms----------------------
create_atoms 1 region box
delete_atoms region r_flow
group g_bot_wall region r_bot_wall
group g_top_wall region r_top_wall
group g_boundary union g_bot_wall g_top_wall
create_atoms 2 random 100000 987654 r_flow_q
group g_flow region r_flow
mass 1 0.0001793448 #179.3448/1.0e6
mass 2 0.000066424 #66.424/1.0e6
#-----------------------------LJ potential------------------
pair_style lj/cut 0.85125
pair_coeff 1 1 83.5 0.2475 0.85125
pair_coeff 2 2 1.657 0.3405 0.85125
pair_coeff 1 2 0.9613 0.2978 0.85125
#------------------------initialization velocity------------
velocity g_boundary create ${T_wall} 987654 dist gaussian
velocity g_flow create ${T_wall} 987654 dist gaussian
#-------------------------spring model for wall atoms----------------------------------------------
fix spring_boundary_wall g_boundary spring/self ${spring_constant}
#-------------------------energy minimization-------------------------------------
minimize 1.0e-6 1.0e-9 50000 20000
min_style cg
#-------------------------thermostat for boundary---------------------------------
compute th_com_wall g_bot_wall temp
fix th_fix_wall g_bot_wall langevin {T_wall} {T_wall} 0.001 699483 #1.0e-3
fix_modify th_fix_wall temp th_com_wall
compute th_com_wall_0 g_top_wall temp
fix th_fix_wall_0 g_top_wall langevin {T_wall} {T_wall} 0.001 699483 #1.0e-3
fix_modify th_fix_wall_0 temp th_com_wall_0
#-------------------------thermostat for liquid------------------------------------
compute th_com_flow g_flow temp
fix th_scale_flow g_flow temp/rescale 1 {T_wall} {T_wall} 0.1 1
fix_modify th_scale_flow temp th_com_flow
fix 1 g_flow nve
fix 2D all enforce2d
fix spring_recenter_1 g_boundary recenter INIT INIT INIT
dump e1 all atom 10000 edump.flow.lammpstrj
thermo 10000
timestep 0.00001 #1e-5
run 200000
#------------------------------------fluid definition-------------------------------
unfix th_scale_flow
unfix 1
unfix 2D
unfix spring_recenter_1
fix flow_definition g_flow addforce 1.0 0.0 0.0
compute fd_com g_flow temp/profile 1 0 0 xy 1 1
fix fd_temp g_flow langevin {T_fluid} {T_fluid} 0.001 699483
fix_modify fd_temp temp fd_com
fix fd_time_1 g_flow nve
fix 2D all enforce2d
fix spring_recenter_2 g_boundary recenter INIT INIT INIT
thermo 10000
thermo_modify temp fd_com
timestep 0.00001 #1e-5
run 300000
#------------------------------------sampling---------------------------------------
compute 1ps_com_20_4 g_flow chunk/atom bin/2d x 0.0 32.0 y 0.0 1.6 units box
compute 1ps_temp_20_4 g_flow temp/chunk 1ps_com_20_4 temp com yes
fix 1ps_rusult_20_4 g_flow ave/chunk 2 500000 1600000 1ps_com_20_4 vx vy density/number temp bias 1ps_temp_20_4 file p_6.out format %20.16g
compute s_com g_flow temp/profile 1 0 0 xy 15 17 out bin
fix s_re g_flow ave/time 2 500000 1600000 c_s_com[1] c_s_com[2] file cc_t.out mode vector format %20.16g
#------------------------------------equilibrium------------------------------------
thermo 10000
timestep 0.00001
run 2130000
