Temperature profile for 2D flow

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

T.jpg

This is just fluid mechanics feedback - A velocity profile v_x (or it’s gradient dv_x/dy) and a few stress profiles may help you with this.