Fix ehex and how to set the energy flux into the reservoir

Dear All,

I’m trying to calculate thermal conductivity of a nanofluid using the direct method. I changed the in.ehex KAPPA exmaple to fit my problem. The input is:

fix hot all ehex 1 0.075 region hot
fix cold all ehex 1 -0.075 region cold

I tried 0.075 and 0.15 as provided in examples of the fix ehex doc page.

I get this error “Fix ehex kinetic energy went negative (src/RIGID/fix_ehex.cpp:269)”

My question is, how do I set the energy flux to solve this error?
How should I choose the value of the parameter “energy flux into the
reservoir” using the real units? How can I use trial and error to get it right?

Best regards,
Yunes Salman

Input file:

``````###############################################
# El Mehdi Achhal, Hicham Jabraoui, Soukaina Zeroual, Hamid Loulijat, Abdellatif Hasnaoui, Said Ouaskit, Modeling and simulations of nanofluids using classical molecular dynamics: Particle size and temperature effects on thermal conductivity, Advanced Powder Technology, Volume 29, Issue 10, 2018
###############################################

# ------------------------ INITIALIZATION ----------------------------
units real
dimension	3
boundary	p	p	p
atom_style	atomic
#bond_style harmonic
#angle_style harmonic
neighbor 2.0 bin

variable latparam equal 5.72

variable    T equal 86
variable    V equal vol
variable    dt equal 1
variable    p equal 20000     # correlation length
variable    s equal 1         # sample interval
variable    d equal \$p*\$s     # dump interval
variable    kB equal 1.3806504e-23    # [J/K] Boltzmann
variable    ev2J equal 1.60217e-19
variable    kCal2J equal 4186.0/6.02214e23
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    convert equal \${kCal2J}*\${kCal2J}/\${fs2s}/\${A2m}
variable dZ equal 3.35
variable dX equal 30

# ----------------------- ATOM DEFINITION ----------------------------
lattice fcc \${latparam}
region   box block  -34.32 34.32 -34.32 34.32 -34.32 34.32 units box
create_box  2 box

region liquid_1 sphere 0 0 0 6 side out units box
region liquid_0 block  -34.32 34.32 -34.32 34.32 -34.32 34.32 units box
region liquid  intersect 2 liquid_1 liquid_0 units box

#liquid
create_atoms 2 region liquid #ratio 0.997974244 74637

region particle_1 sphere 0 0 0 6 units box
#create_atoms 1 region particle_1 #ratio 0.363636364 74637
create_atoms 1 random 14 489571 particle_1

group solid type 1
group liquid type 2

mass       1     63.546  #Cu
mass       2     39.948 #Ar

# ------------------------ FORCE FIELDS ------------------------------
pair_style	hybrid/overlay lj/cut 8.5 eam
#pair_style	eam
pair_coeff 1 1 eam Cu_u3.eam

pair_coeff  2 2 lj/cut  0.24068 3.405 8.5
pair_coeff  1 2 lj/cut 1.499 2.872 8.5

#pair_style lj/cut 8.8
#pair_coef type_i type_j epsilon sigma rcut
#pair_coeff  2 2  0.24068 3.405 8.8
#pair_coeff  1 1 9.4584 2.3377 8.8
#pair_coeff  1 2 1.499 2.871 8.8

info system
print "Finished INFO"

# ---------------------Minimization---------------------------
reset_timestep	0
min_style  cg
minimize  1.0e-13 1.0e-12 1000 10000

print "Finished Minimizing"

# ------------------------- Heat Regions -----------------------------

variable T equal   86 	# target temperature
variable thi equal 116 	# temperature of hot region
variable tlo equal 56 	# temperature of cold region

variable Lhot equal -34.32
variable Rhot equal  -28.32
variable Lcold equal 28.32
variable Rcold equal 34.32

region hot block \${Lhot} \${Rhot} INF INF INF INF units box	# Define hot region X1,X2,Y1,Y2,Z1,Z2 (INF => PBC)
region cold block \${Lcold} \${Rcold} INF INF INF INF units box	# Define cold region

compute         Thot all temp/region hot
compute         Tcold all temp/region cold

# ---------------------1st equilibration run-----------------------------
reset_timestep	0
timestep \${dt}

velocity all create 86.0 54664

fix 1 solid rigid/nvt group 1 solid temp 86.0 86.0 \$(100.0*dt) reinit no
fix 2 liquid nvt temp 86.0 86.0 \$(100.0*dt) tchain 1

compute 1 solid temp
compute 2 liquid temp

#Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp density c_1 c_2

dump 	1 all custom 1000 dump.cunp.lammpstrj id type x y z mass

run 2000000

unfix 1
unfix 2

# ---------------------2nd equilibration run-----------------------------

fix		1 all nve
fix             hot all ehex 1 100.0 region hot
fix             cold all ehex 1 -100.0 region cold

thermo_style    custom step temp c_Thot c_Tcold
thermo	        1000
run             1000000

# ------------Thermal_Conductivity_Calculation----------------

compute		ke all ke/atom
variable	temp atom c_ke/1.5

compute         layers all chunk/atom bin/1d z lower 0.05 units reduced
fix		2 all ave/chunk 10 100 1000 layers v_temp file profile.ehex

variable        tdiff equal f_2[11][3]-f_2[1][3]
fix             ave all ave/time 1 1 1000 v_tdiff ave running start 13000
thermo_style    custom step temp c_Thot c_Tcold v_tdiff f_ave

run             1000000
``````

Obviously, you must use a smaller amount of energy.

How much depends on the details of the method. You need to study the corresponding publication(s) describing it and you can use `compute ke/atom` and `compute reduce/region` to determine the current amount of kinetic energy in those regions.

Please also note that your script needs additional modifications for real units, since many conversion constants and parameters are 1.0 in reduced units, but not in other settings.
Again, you need to refer to the corresponding text book(s) or publication(s) describing the methods and properties how to correctly determine and convert them.

1 Like