REG: SUPERCRTICAL BASED NANOFLUID (Cu - scCO2)

Dear all,
I am trying to find out the thermal conductivity of a supercritical based nanofluid( base fluid CO2 with copper nanoparticle) by Green Kubo formalism. I validated the code with the available NIST data for supercritical CO2 and extend the code by simply putting Cu parameters and potentials. I have taken the large correlation length (p =4000) and very small sample length (s=2). I have also checked that the system reached the equilibrium before I start the post processing. I have taken the Green kubo code from the lammps directory for the post processing (please refer to the code shown below). The system has 9600 CO2 molecules with 1205 Cu atoms (3nm diameter). The simulation domain is cube of side 100 Angstrom. I have also made regions around the nanoparticle with bin size 10 A (radially), to check for the density distribution around the nanoparticle whether it is forming nanolayer or not.

Problem I am facing:

  1. I am getting very high value of thermal conductivity when i put nanoparticle in the system (200 times approximately) at 310 K
  2. The density distribution shows that near the nanoparticle region , there is no layer formation and also the density far away from the nanoparticle is not coming as bulk CO2 density.
  3. Is there a standard bin size to make for regions to check whether nanolayer formed has more density than the bulk fluid?
  4. Also the peak in the rdf plot is more in sc-nanofluid (scCO2-Cu) when compared to the rdf plot of supercritical co2 only. Is this okay?
  5. I am also getting this error : Fix in variable not computed at compatible time. I have found the solution of this in the LAMMPS forum but could not able to apply to my code and still getting the same error.

Please suggest me where I am wrong and what changes should I incorporate in the code attached to get the correct physics.

#CODE

units real
atom_style full

variable T equal 310
variable V equal vol
variable dt equal 0.1
variable p equal 4000 # correlation length
variable s equal 2 # sample interval
variable d equal $p*$s # dump interval

convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
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}

LJ parameters for CO2 and Cu

read_data data.9600sco2
pair_style lj/cut/coul/long 10 10
pair_coeff 1 1 0.164254 3.0
pair_coeff 2 2 0.0573209 2.7918
pair_coeff 1 2 0.097 2.9
pair_coeff 1 3 3.7788 2.6835
pair_coeff 2 3 2.23427 2.5485
pair_coeff 3 3 9.45 2.337
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4
bond_coeff 1 469 1.16
angle_coeff 1 147.5 180.0

Cube of 100 A side

neighbor 2.0 bin
neigh_modify every 10 delay 20 check yes

#Grouping

#1=O
#2=C
group co2 type 1 2
group np type 3

Regions made for density distribution (bin size 10)

region 0 sphere 50.0 50.0 50.0 15 side out units box
region 1 sphere 50.0 50.0 50.0 25 side in units box
region 01 intersect 2 0 1

region 2 sphere 50.0 50.0 50.0 25 side out units box
region 3 sphere 50.0 50.0 50.0 35 side in units box
region 23 intersect 2 2 3

region 4 sphere 50.0 50.0 50.0 35 side out units box
region 5 sphere 50.0 50.0 50.0 45 side in units box
region 45 intersect 2 4 5

region 6 sphere 50.0 50.0 50.0 45 side out units box
region 7 sphere 50.0 50.0 50.0 55 side in units box
region 67 intersect 2 6 7

region 8 sphere 50.0 50.0 50.0 55 side out units box
region 9 sphere 50.0 50.0 50.0 65 side in units box
region 89 intersect 2 8 9

region 10 sphere 50.0 50.0 50.0 65 side out units box
region 11 sphere 50.0 50.0 50.0 75 side in units box
region 1011 intersect 2 10 11

region 12 sphere 50.0 50.0 50.0 75 side out units box
region 13 sphere 50.0 50.0 50.0 85 side in units box
region 1213 intersect 2 12 13

region 14 sphere 50.0 50.0 50.0 85 side out units box
region 15 sphere 50.0 50.0 50.0 95 side in units box
region 1415 intersect 2 14 15

region 16 sphere 50.0 50.0 50.0 95 side out units box
region 17 block 0.0 100.0 0.0 100.0 0.0 100.0
region 1617 intersect 2 16 17

min_style cg
dump 1 all custom 1000 water2_k_MSD_RDF2.lammpstrj id mol type x y z ix iy iz
thermo 1000
thermo_style custom step etotal temp vol press density
minimize 1e-5 1e-7 5000000 10000000

equilibration and thermalization

neighbor 3.0 bin
neigh_modify delay 0 every 1 check yes
neigh_modify exclude molecule co2
velocity all create $T 102486 mom yes rot yes dist gaussian

NVE

fix NVE all nve
thermo_style custom step etotal temp vol press density
thermo 1000
timestep ${dt}
run 1000000

unfix NVE

NVT

fix NVT all nvt temp $T $T 100

variable a1 equal count(co2,01)
variable a2 equal count(co2,23)
variable a3 equal count(co2,45)
variable a4 equal count(co2,67)
variable a5 equal count(co2,89)
variable a6 equal count(co2,1011)
variable a7 equal count(co2,1213)
variable a8 equal count(co2,1415)
variable a9 equal count(co2,1617)

thermo_style custom step etotal temp vol press density v_a1 v_a2 v_a3 v_a4 v_a5 v_a6 v_a7 v_a8 v_a9
timestep ${dt}
thermo 1000
dump 2 all custom 1000 cd2000_conductivity.lammpstrj id mol type x y z vx vy vz ix iy iz
run 1000000

compute msd co2 msd com yes
fix msd co2 ave/time 2 100 1000 c_msd[4] file test.msd

compute rdf co2 rdf 1000 1 1 1 2 2 2
fix rdf co2 ave/time 2 100 1000 c_rdf file test.rdf mode vector

thermo 1000

timestep ${dt}
thermo_style custom step etotal temp vol press density
run 1000000

compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom NULL virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix JJ all ave/correlate $s $p d & c_flux[1] c_flux[2] c_flux[3] type auto file J1Jt.dat ave running variable scale equal {convert}/${kB}/$T/$T/$V*s*{dt}
variable k11 equal trap(f_JJ[3]){scale} variable k22 equal trap(f_JJ[4])*{scale}
variable k33 equal trap(f_JJ[5])
{scale} thermo 1000 thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 timestep {dt}
run 1000000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print “average conductivity: $k[W/mK] @ T K, {ndens} /A^3”