problem with evaporation on a flat surface

Hello Lammps users and developers,
I am trying to simulate the evaporation of liquid Ar on a flat copper
surface at T=140K. The solid surface contains 9 layers of copper
atoms. The first layer of Copper atoms stayed as a boundary wall to
keep the volume of the system constant (Group S1); the following
inside two layers were set as heat source from which heat flux was
generated (Group S2); the last six layers were set as solid walls
through which heat conducted to argon fluids (Group S3). The group
"SOLID" and "liquid", "Vapor" and "LV" contain all solid atoms,liquid
atoms, vapor atoms and liquid+vapor atoms, respectively,
The temperature history during equlibration period is O.K and all the
system is in T=90K but in Evaporation period the temperature of Ar
liquid (T_LV) increases (more than 150 K) and get higher than the
temperature of solid surface (temperature of group 3) which is not
realistic. Would you please help me what is wrong in my input script
and simulation set up?

units metal
dimension 3
boundary p sf p
atom_style atomic
neighbor 2 bin
neigh_modify delay 5

region box block 0 45 0 600 0 45 units box ### simulation box
create_box 2 box

lattice fcc 3.61 ### lattice for solid region
region 1 block 0 45 0 1 0 45 units box ### fix solid atoms region
region 2 block 0 45 1 4.9 0 45 units box ### The region above
the first layer (second region)
region 3 block 0 45 4.9 15 0 45 units box ### last layers of solid wall
create_atoms 1 region 1
create_atoms 1 region 2
create_atoms 1 region 3
lattice none

lattice fcc 5.76 ### lattice for liquid region
region 4 block 0 45 15 30 0 45 units box ### liquid atoms region
create_atoms 2 region 4
lattice none

lattice fcc 32.86 ### lattice for vapor region
region 5 block 0 45 30 600 0 45 units box ### vapor atoms region
create_atoms 2 region 5
lattice none

group S1 region 1
group S2 region 2
group S3 region 3
group liquid region 4
group vapor region 5
group LV union vapor liquid ### liquid+vapor atoms

group solid union S1 S2 S3
group whole1 subtract all S1 ### all atoms-wall_fix
group whole2 subtract all S2 ### all atoms - S2

mass 1 63.5463 ### mass of solid atoms
mass 2 39.948 ### mass of liquid and vapor atoms
set group solid type 1
set group LV type 2

pair_style lj/cut 12
pair_coeff 1 1 0.409598855 2.34 12 ### S-S (free layer with itself)
pair_coeff 2 2 0.0104233 3.40 12 ### L-L, L-V and V-V
pair_coeff 1 2 0.065007161 2.8689 12 ### L-S and V-S(free layers)

compute T3 S3 temp
compute_modify T3 dynamic yes

compute T_LV LV temp
compute_modify T_LV dynamic yes

min_style cg
minimize 1.0e-6 1.0e-8 10000 10000
reset_timestep 0

velocity S1 set 0.0 0.0 0.0 units box
velocity whole1 create 90 482748 units box
fix 1 S1 setforce 0.0 0.0 0.0
fix 2 all nve
fix 3 solid langevin 90 90 100.0 699483
fix 4 LV langevin 90 90 100.0 699483
fix 6 all wall/reflect yhi EDGE units box

thermo 1000
thermo_style custom step c_T3 c_T_LV
thermo_modify line one lost error norm no
dump 1 all atom 500 Cu-flat.vmd

timestep 0.001
run 100000
############# Equlibration Period#################
unfix 3
unfix 4
fix 5 S2 langevin 90 90 100.0 699483
run 200000
############# Evaporation Period#################
unfix 2
unfix 5
fix 8 S2 nvt temp 140 140 0.1
fix 9 whole2 nve
run 2000000

Hello Lammps users and developers,
I am trying to simulate the evaporation of liquid Ar on a flat copper
surface at T=140K. The solid surface contains 9 layers of copper
atoms. The first layer of Copper atoms stayed as a boundary wall to
keep the volume of the system constant (Group S1); the following
inside two layers were set as heat source from which heat flux was
generated (Group S2); the last six layers were set as solid walls
through which heat conducted to argon fluids (Group S3). The group
"SOLID" and "liquid", "Vapor" and "LV" contain all solid atoms,liquid
atoms, vapor atoms and liquid+vapor atoms, respectively,
The temperature history during equlibration period is O.K and all the
system is in T=90K but in Evaporation period the temperature of Ar
liquid (T_LV) increases (more than 150 K) and get higher than the
temperature of solid surface (temperature of group 3) which is not
realistic. Would you please help me what is wrong in my input script
and simulation set up?

just visualize your trajectory and it will be obvious.

your system is falling apart already in the first "run"
statement. just looking at the temperature is not enough.

axel.

Hi Hossein

i checked your script. the system collapses.
during the 1st run the fluid temperature increases over 1000K and in the last run
the region 2 act as a heat sink and the fluid temperature decreases to about 400K.

first of all you should fix the atoms in region 1. for that purpose you should set their initial velocity
to zero via velocity command and then be care that they not to be included in any of the fixes.
it makes their positions remain constant because they never will be time integrated.
for confidence you can also fix two layers in right and left in a similar way.

before the evaporation apply the Nose-Hoover thermostat on both solid and fluid phases for some steps.
(for finding appropriate temperature according to your lattice-constants just take a look at the
argon phase-diagram attached)

and finally increase the solid temperature. let the solid in nvt ensemble whilethe fluid ensemble is nve.

best,
Fahim

phaseDiagram_argon.png

Dear Alex and Fahim,
I really appreciate your helps. You are right, The input script isn't
work properly. I think that I made a mistake during the simplifying my
main input script because it was a bit complicated so I decided to
delete unnecessary parts and upload a simpler one to LAMMPS user
forum. I will send the corrected version of it ASAP.
Thanks agian

I sent an incorrect version of my input scripts.

Dear Steve, Alex,Fahim and lammps users,
I tried to simplify my input script again, please see the revised
version of it. As I mentioned in my first email, Before starting the
evaporation, everything is O.K but during the Evaporation the
temperature of Ar liquid (T_LV) increases (more than 150K) and get
higher than the temperature of solid surface (T=140K) which I couldn't
find the reason of it. I also have another problem with pressure
calculation. The calculation of pressure during the equilibration
period is correct and my results have good agreement with experimental
data (please see the attached files) while during the evaporation it
seems that the value and trend of pressure is not correct. Please see
the attached files.
I guess (i) maybe the upper boundary condition (wall/reflect) produce
additional energy in the system but I am not sure, (ii) something
wrong is in integration part of the script.

Thank you very much for your kind help,
Best Regards
Hossein

Input Script:

variable rcut equal 12.5
units metal
dimension 3
boundary p sf p
atom_style atomic
neighbor 2 bin
neigh_modify delay 5

region box block 0 72 0 600 0 72 units box ### simulation box
create_box 2 box

lattice fcc 3.61 ###
lattice for solid region
region 1 block 0 72 0 1 0 72 units box ### fix solid atoms region (S1)
region 2 block 0 72 1 4.9 0 72 units box ### The region above
the first layer (second region i.e, S2)
region 3 block 0 72 4.9 15 0 72 units box ### last layers of
solid wall (S3)
create_atoms 1 region 1
create_atoms 1 region 2
create_atoms 1 region 3
lattice none

lattice fcc 5.76 ###
lattice for liquid region
region 4 block 0 72 15 30 0 72 units box ### liquid atoms region
create_atoms 2 region 4
lattice none

lattice fcc 32.86 ###
lattice for vapor region
region 5 block 0 72 30 600 0 72 units box ### vapor atoms region
create_atoms 2 region 5
lattice none

group S1 region 1
group S2 region 2
group S3 region 3
group liquid region 4
group vapor region 5
group LV union vapor liquid ### liquid+vapor atoms

group solid union S1 S2 S3
group whole1 subtract all S1 ### all atoms-wall_fix
group whole2 subtract all S2 ### all atoms - S2

mass 1 63.5463 ### mass of solid atoms
mass 2 39.948 ### mass of liquid and vapor atoms
set group solid type 1
set group LV type 2

pair_style lj/cut ${rcut}

pair_coeff 1 1 0.409598855 2.34 \{rcut\} pair\_coeff 2 2 0\.0104233 3\.40 {rcut}
pair_coeff 1 2 0.065007161 2.8689 ${rcut}

compute T3 S3 temp
compute_modify T3 dynamic yes

compute T_LV LV temp
compute_modify T_LV dynamic yes

compute peratomV vapor stress/atom
compute pV vapor reduce sum c_peratomV[1] c_peratomV[2] c_peratomV[3]
variable P_V equal -(c_pV[1]+c_pV[2]+c_pV[3])/(3*3032640) ###
3032640 is the volume of LV region

min_style cg
minimize 1.0e-6 1.0e-8 10000 10000
reset_timestep 0

group SS union S2 S3
velocity S1 set 0.0 0.0 0.0 units box
velocity whole1 create 90 482748 units box
fix 1 S1 setforce 0.0 0.0 0.0
fix 2 all nve
fix 3 SS langevin 90 90 100.0 699483
fix 4 LV langevin 90 90 100.0 699483
fix 6 all wall/reflect yhi EDGE units box

thermo 1000
thermo_style custom step c_T3 c_T_LV v_P_V
thermo_modify line one lost error norm no
dump 1 all atom 500 Cu-flat.vmd

timestep 0.001
run 100000
############# Equlibration Period#################
unfix 3
unfix 4
fix 5 S2 langevin 90 90 100.0 699483
run 200000
############# Evaporation Period#################
unfix 2
unfix 5
fix 8 S2 nvt temp 140 140 0.1
fix 9 whole2 nve
run 2000000

equlibrium pressure.bmp (2.2 MB)

pressure.bmp (1.5 MB)

temperature.bmp (1.33 MB)

A complicated input script is like a program. The
only way to debug it is to isolate one option
at a time and test that it works like you expect.
No one else is likely to do that for you.

Steve

Steve:
Thank you very much for your answer. I think higher temperature of
Argon atoms with respect to solid wall is due to the uncorrect
boundary condition at top of system (fix wall/reflect) that I used in
my problem. Probably, the "fix wall/reflect" add energy to the system
which cause the temperature liquid get higher. Is the collision of
atoms with this boundary elastic?
I guess after collision of atoms the energy of them increases. Am i right?
Thanks again for your helps in advance,
Best Regards,
Hossein

Steve:
Thank you very much for your answer. I think higher temperature of
Argon atoms with respect to solid wall is due to the uncorrect
boundary condition at top of system (fix wall/reflect) that I used in
my problem. Probably, the "fix wall/reflect" add energy to the system
which cause the temperature liquid get higher. Is the collision of
atoms with this boundary elastic?

you thank steve for his advice, but you don't follow it,
but rather re-iterate what you stated before. if you didn't
get an answer before, why should you get one now??

I guess after collision of atoms the energy of them increases. Am i right?

very few problems are solved by guessing,
but rather by systematic exclusion of possible
sources. (or paraphrasing from sir arthur conan doyle:
the way to find the solution, my dear watson,
is to eliminate the impossible. whatever is left,
however improbable, has to be the solution).

in this case, your guess has no foundation.
where please should the additional energy
be coming from? all the wall does is to
change the sign of the velocity at the point
where a particle hits the wall.
1/2*m*v**2 = 1/2*m*(-v)**2

Thanks again for your helps in advance,

you have to give yourself a hand here.
take your input and then uncomment
everything and the add the various features
and gimmicks step by step and monitor
when things start to go south. then you
will know what is triggering the problem
and you are half the way to understand
what is the cause and then you'll have
the solution.

short of offering an attractive "finder's fee",
i doubt that you'll find somebody to do
what is ultimately your job.

cheers,
    axel.

Hi

Yes the wall will flip the sign of the velocity i.e the colission is elastic. However I’m thinking the wall could act as a momentum mirror and genreate an shock wave or expansion wave that propagate out of the wall. however for generate a steady wave all velocity -up should be added to the particles. This phenomena may cause the atoms sorrounding the wall to have a higher kinetic energy than the ones far away from the wall

PS: I havent read all the messages in this thread, so i’m only guessing what could happen …Shame on me if im wrong or just saying crazy things.

Oscar G.

Alex:
Thank you for your quick response. In my last email, I asked a new
question about "fix wall/reflect". I didn't re-iterate my first
question.
Best Regards,
Hossein

Alex and Oscar:
I really appreciate your helps and guidance.
Best Regards,
Hossein

Axel is correct that fix wall/reflect should not add kinetic
energy to the system. On a per-particle basis, energy
is conserved in the wall collision.

Steve

Steve:
Thanks a lot for your answer.