Dear all,
I am doing a heat conduction between two parallel silver walls separated by liquid argon. I thermostat one silver wall side to be low temperature, and the other silver wall side to be high temperature, so a steady heat flux is generated in the system. In my case, the heat flux at x direction should be a constant, and in other directions should be zero. But my heat flux results are like this:
time step Jx Jy Jz
2995000 6.50608e-05 -1.64195e-05 -3.25935e-05
2996000 -1.1548e-06 3.25386e-05 -1.9882e-05
2997000 -1.26967e-05 7.93451e-05 -2.38833e-05
2998000 6.46115e-05 -4.37606e-06 5.89027e-05
2999000 -2.30768e-05 4.15145e-05 -3.38428e-05
3000000 -4.3106e-05 -3.81433e-06 2.88722e-05
Below is my script:
units metal
boundary f p p
atom_style full
read_restart 110-90.restart.1000000
pair_style hybrid eam lj/cut 10.0
pair_coeff 2 2 lj/cut 0.0103 3.405
pair_coeff 1 1 eam Ag_u3.eam
pair_coeff 1 2 lj/cut 0.06013 2.9895
region hol block -30 -29 INF INF INF INF units box
region col block 57 58 INF INF INF INF units box
group ho-layer region hol
group co-layer region col
group argon type 2
group hotwall1 id <= 648
group coldwall1 subtract all argon hotwall1
group hotwall subtract hotwall1 ho-layer
group coldwall subtract coldwall1 co-layer
group myall subtract all ho-layer co-layer
timestep 0.004
fix 3 hotwall langevin 110 110 0.01 699483
fix 4 coldwall langevin 90 90 0.01 890765
fix 5 myall nve
compute pe myall pe/atom pair
compute ke myall ke/atom
compute stress myall stress/atom pair
compute flux myall heat/flux ke pe stress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix 90 myall ave/time 1 1000 1000 v_Jx v_Jy v_Jz file flux.date
thermo 10000
thermo_style custom step v_Jx v_Jy v_Jz
run 2000000
Is there anything I did wrong about heat flux calculation?
Thanks!
ziyuan