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
#read_data argoncopper1per.data
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