electrostatics

Submitted by pain08 on Mon, 06/01/2020 - 17:35

hello everyone i hope ur all doing well i am a new ligggths user i need help i want to add an electrostatic field to attract granuls or push them if possible with positive and negative voltage to make a granular separation i am trying to modify the chute wear scrit by adding efield command but it's not working getting the following error thank you all

atom_style hybrid charge granular
atom_modify map array
boundary f f f
newton off

communicate single vel yes

units si

region domain block -0.5 0.1 -0.2 0.2 -0.4 0.15 units box
create_box 1 domain

neighbor 0.002 bin
neigh_modify delay 0

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global k_finnie peratomtypepair 1 1.0

#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *

timestep 0.00001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix kick external-field efield 1.0 1.0 1.0

#the chute
fix cad all mesh/surface/stress file meshes/simple_chute.stl type 1 wear finnie
fix inface all mesh/surface file meshes/insertion_face.stl type 1
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes cad

#distributions for insertion
fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2500 radius constant 0.0015
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2500 radius constant 0.0025
fix pdd1 all particledistribution/discrete 32452843 2 pts1 0.3 pts2 0.7

#region and insertion
group nve_group region domain
region bc cylinder z 0.0 0.0 0.015 0.05 0.12 units box

#particle insertion
fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
nparticles 6000 massrate 0.1 insert_every 1000 overlapcheck yes all_in no vel constant 0.0 0.0 -1.0 &
insertion_face inface

#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere

#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no

#insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 200 post/chute_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstress all mesh/gran/VTK 200 post/mesh_*.vtk stress wear cad

#insert particles
run 100000 upto
unfix ins

ERROR: Could not find fix group ID (../modify.cpp:782)

pain08 | Tue, 06/02/2020 - 17:10

thank you that solved my problm unfortunatly i have another error
ERROR: Fix efield requires atom attribute q (../fix_efield.cpp:134)
if u could help me with that i'd appriciatee it thank you again for ur help

tjleps | Wed, 06/03/2020 - 23:46

You need to use "atom_style hybrid" charge granular with "fix_efield"

pain08 | Thu, 06/04/2020 - 21:02

thank you so much i made a mistake by using hybride sphere granular and by changing it to charge granular it worked great i just hav a little question the program is working fine but i don't see the granuls changing direction because of the electric field that's why i hav some doubts i modifyied many variables to make some particles go right and others go left like for the electrostatic separation by appliying like 20 kv on a side and -20 kv the other side for example if u have and idea on what i should modify in the script it would mean the world to me and thanks for ur help :

#Simple chute wear test

atom_style hybrid charge granular
atom_modify map array
boundary f f f
newton off

communicate single vel yes

units si

region domain block -0.5 0.1 -0.2 0.2 -0.4 0.15 units box
create_box 1 domain

neighbor 0.002 bin
neigh_modify delay 0

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global k_finnie peratomtypepair 1 1.0

#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *

timestep 0.00001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix efield all efield 0 200000 0
fix efield all efield 0 -20000 0

#the chute
fix cad all mesh/surface/stress file meshes/simple_chute.stl type 1 wear finnie
fix inface all mesh/surface file meshes/insertion_face.stl type 1
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes cad

#distributions for insertion
fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2500 radius constant 0.0015
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2500 radius constant 0.0025
fix pdd1 all particledistribution/discrete 32452843 2 pts1 0.3 pts2 0.7

#region and insertion
group nve_group region domain
region bc cylinder z 0.0 0.0 0.015 0.05 0.12 units box

#particle insertion
fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
nparticles 1000 massrate 0.1 insert_every 1000 overlapcheck yes all_in no vel constant 0.0 0.0 -1.0 &
insertion_face inface

#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere

#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no

#insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 200 post/chute_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstress all mesh/gran/VTK 200 post/mesh_*.vtk stress wear cad

#insert particles
run 100000 upto
unfix ins

tjleps | Thu, 06/04/2020 - 21:44

You need to set the charge for the atoms, right now they're all uncharged so they won't be affected by the efield. Also fix_efield sets an electric field, not a potential, you're just overriding one field with the other.

pain08 | Thu, 06/04/2020 - 22:07

wish commands and lines should i change sir can u help me correct it by charging particles and create a potential please thank you for ur precious help

tjleps | Thu, 06/04/2020 - 23:12

I recommend you go through a couple unedited tutorial scripts, looking up each command in the documentation, until you fully grasp what the structure of a LIGGGHTS script is, it will help a lot when you start making your own scripts. In the meantime, use the set command:

set atom all charge 0.000001
will put one microcoulomb on every atom in your system.

fix_efield 0 50000 0
will create the equivalent of a 20kv potential across your simulation box since V=E*d for a constant electric field.

pain08 | Fri, 06/05/2020 - 01:33

thank you verry much for ur help sir i hope u could guide me with ligggths and if u have publications that i can site email them to me i hope we can work on a project someday meanwhile here my email painrikosu8@gmail.com and thanks again for sharing ur knowledge

tjleps | Wed, 06/10/2020 - 02:01

No problem, it can be a little bit of a steep learning curve. I don't have any publications that would make sense to cite I don't think, but I'm glad I could help.

pain08 | Sun, 07/05/2020 - 16:46

hello tjleps i hope ur doing great i have a little problm i hope u can help me with i have trouble charging particles the set command u told me about only charges atom or molecules i need it to charge particles like granuls the adapt command is not working properly either can you please help me i am kinda stuck on it thank you

tjleps | Sun, 07/05/2020 - 19:44

If you could attach some code and specify more precisely, exactly what you're trying to have happen, it would help with trying to figure out a solution.

pain08 | Mon, 07/06/2020 - 20:03

i am trying to do something like this example https://www.cfdem.com/electrostatics with two different granular particles types that are mixed and two electric fields so that type one goes left and type two go right if u could help me it would be great and thank you so much ur help is priceless my friend

#falling particles and way flow coupling

atom_style granular
atom_modify map array
boundary f f f
newton off

communicate single vel yes

units si

region domain block -0.155 0.155 -0.015 0.015 -0.005 0.235 units box
create_box 1 domain

neighbor 0.002 bin
neigh_modify delay 0

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5

#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *

timestep 0.00001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

#distributions for insertion
fix pts all particletemplate/sphere 15485863 atom_type 1 density constant 1000 radius constant 0.004

fix pdd all particledistribution/discrete 15485867 1 pts 1.0

region factory block -0.14 -0.10 -0.006 0.006 0.215 0.225 units box
lattice sc 0.0001
create_atoms 1 region factory

fix ins all insert/rate/region seed 32452843 distributiontemplate pdd &
nparticles 1000 massrate 1.66666667 insert_every 1000 &
overlapcheck yes vel constant 0. 0. -1. region factory ntry_mc 10000

# Geometry
fix wal11 all wall/gran model hertz tangential history primitive type 1 zplane 0.0
fix wall2 all wall/gran model hertz tangential history primitive type 1 xplane -0.15
fix wall3 all wall/gran model hertz tangential history primitive type 1 xplane 0.15
fix wall4 all wall/gran model hertz tangential history primitive type 1 yplane -0.01
fix wall5 all wall/gran model hertz tangential history primitive type 1 yplane 0.01

#apply efield and atom charges

set atom all charge 0.000001
fix efield all efield 0 50000 0

#apply nve integration to all particles that are inserted as single particles
fix integr all nve/sphere

# Thermo settings
thermo_style custom step atoms ke cpu
thermo 10000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

# Run 1 step to check timestep and initialize
fix ctg all check/timestep/gran 1 0.01 0.01
run 1
unfix ctg
# Dump particle positions
dump dmp all custom/vtk 100 post/fall*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
### Execution and further settings
# Run remainder
run 200000 upto

tjleps | Mon, 07/06/2020 - 21:47

I'm not sure what you mean by "charges atom or molecules i need it to charge particles like granuls" Atoms in LIGGGHTS are generic to all particles including grains, owing to it's molecular dynamics roots. The set command as you have it should work to give your grains a charge.

I can try to run your code to see if I can figure out what's going on, but it could be a few days.

pain08 | Mon, 07/20/2020 - 17:41

hi tjleps hope ur doing great i wanted to ask you about the dump files i cannot get atom charge or electric charge after adding q to dump i get zeros hope u can help thanks
dump dmp all custom/vtk 100 post/fall*.vtk id type type x y z ix iy iz vx vy vz fx fy fz q omegax omegay omegaz radius

tjleps | Mon, 07/20/2020 - 22:10

The reason is that the set command only applies to atoms after they are inserted, so you need to set the charge after you've done your insertion commands. I'm not sure how, or if it's possible, to set the charge before an insertion. If you move your set command to after the insertion step it should work fine.

Also your lattice is extremely small (2.5% of your radius) in your example code, resulting in millions of overlapping particles.

pain08 | Thu, 07/23/2020 - 18:39

hi tjleps hope ur doing great i also wanted to ask you about the dump files i cannot get atom charge or electric charge after adding q to dump i get zeros hope u can help thanks

zihengw | Thu, 09/08/2022 - 19:16

Hi, @tjleps and @pain08. Your discussion about the electrostatics module is very helpful. I think my need is very similar to yours @pain08, which is to simulate the triboelectric between granular materials. I am encountering the following two problems, maybe you have already developed some insights that you are willing to share.

1. The first problem is similar to what @pain08 was facing near the end of this post. I was trying to modify the chute_wear example problem by introducing an electric field (fix electricfield all efield 23.34 0.0 0.0), changing the atom style to hybrid (atom_style hybrid charge granular), and attempting to assign charge to individual granular particle (set atom 1 charge -29.9792). Please feel free to check out the input file below. As you may know, for the chute_wear simulation, the granular particles were generated using fix insert/stream command. The objects generated by this command seem to be granular style not atom style. So when I ran the attached input file, the set atom 1 charge -29.9792 command was upset because it couldn't find any atoms. The exact error message was ERROR: Set command with no atoms existing (../set.cpp:108) . Just curious if anyone has successfully assigned charges to the granular particles generated by fix insert/stream, fix insert/pack or other similar fix command. Because I'd love to leverage the capability of fix insert commands on controlling particle size distribution while the particles can be charged. However, I don't have much hope as the description of the set command only mentioned supporting read_data, read_restart or create_atoms commands.

Just a side note, I've tried create_atoms command to create an atom style particle. With the set command, the particle did respond to the electric field. My workaround right now is to generate particles (with controlled distribution) and write them into a txt file; then, use read_data command to finish the particle generation process for the simulation.


#Electrostatic powder spreader simulation

atom_style hybrid charge granular
atom_modify map array
boundary f f f # fixed boundries -> particles will be deleted if leaving the simulation box
newton off # default

communicate single vel yes # default

units cgs

#Definition of boundaries

variable xmin equal -20.5
variable xmax equal 20.5

variable ymin equal -15.5
variable ymax equal 15.5

variable zmin equal -1.5
variable zmax equal 13.5

#Definition of the time step

variable dt equal 1e-5 #Each iteration step represents the defined seconds of the real time.

#Definition of material properties

variable nmaterials equal 2 #Define how many types of materials in the simulation.
# e.g., #1 -> particles and/or parts of the spreader; #2-#n -> other parts of the sprader

#Young Modulus

variable youngmodulus1 equal 1.9e+8 #barye = 1.0e-6 bars
variable youngmodulus2 equal 2.7e+6 #barye = 1.0e-6 bars

#When Youngs Modulus is > 1.e10 use the command 'hard_particles yes'

#hard_particles yes

#Poission ratio

variable poission1 equal 0.28
variable poission2 equal 0.39

#Coefficient of restitution

variable CoR11 equal 0.6
variable CoR12 equal 0.6
variable CoR21 equal 0.6
variable CoR22 equal 0.6

#Sliding friction coefficient

variable sf11 equal 0.3
variable sf12 equal 0.36
variable sf21 equal 0.36
variable sf22 equal 0.0

#Rolling friction coefficient

#variable rf11 equal 0.8
#variable rf12 equal 0.8
#variable rf21 equal 0.8
#variable rf22 equal 0.8

#Particle sizes and properties

variable nradii equal 2
variable radius1 equal 0.025 #cm
variable radius2 equal 0.01 #cm
variable frac1 equal 0.5
variable frac2 equal 0.5

variable density equal 7.85 #g/cm³

variable charges atom mass/100

#filling parameters

variable filltime equal 0.005 #seconds
variable fillmass equal 1 #gram
variable fillmassrate equal ${fillmass}/${filltime} #g/s
variable fillsteps equal ${filltime}/${dt} #Convert time to simulation steps

#settle time

variable settletime equal 1 #seconds
variable settlesteps equal ${settletime}/${dt} #Convert time to simulation steps

#deposition time

#variable depositiontime equal #seconds
#variable depositionsteps equal ${depositiontime}/${dt} #Convert time to simulation steps

#Definition of simulation domain

region domain block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} units box
create_box 2 domain
#create_atoms 1 single 0.0 0.0 3.0

neighbor 0.01 bin #Param controls how frequent the neighbor lists are updated (0.1 cm by default)
neigh_modify delay 0 #Default

#Definition of the contact model

pair_style gran model hertz tangential history #Contact model
pair_coeff * * #Default

timestep ${dt}

fix integrator all nve/sphere

fix ts_check all check/timestep/gran 100 0.1 0.1 warn yes

fix gravi all gravity 980.665 vector 0.0 0.0 -1.0 #Gravity (cm/s²)

fix electricfield all efield 23.34 0.0 0.0 #Electric field (statvolt/cm)

#Definition of mterial properties required for contact

fix m1 all property/global youngsModulus peratomtype ${youngmodulus1} ${youngmodulus2}
fix m2 all property/global poissonsRatio peratomtype ${poission1} ${poission2}
fix m3 all property/global coefficientRestitution peratomtypepair ${nmaterials} ${CoR11} ${CoR12} ${CoR21} ${CoR22}
fix m4 all property/global coefficientFriction peratomtypepair ${nmaterials} ${sf11} ${sf12} ${sf21} ${sf22}
#fix m5 all property/global coefficientRollingFriction peratomtypepair ${nmaterials} ${rf11} ${rf12} ${rf21} ${rf22}
#fix m6 all property/atom q scalar charge

#Input the geometries

fix pan all mesh/surface file mesh/powder_pan.stl type 1 scale 0.1 #Scale 0.1 because the stl was exported in unit of mm
fix plate all mesh/surface file mesh/build_plate.stl type 1 scale 0.1
#fix electrode all mesh/surface file mesh/electrode.stl type scale 0.1

fix walls all wall/gran model hertz tangential history mesh n_meshes 2 meshes pan plate #Define the particle-wall interactions

#Distributions for particle insertion

fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant ${density} radius constant ${radius1}
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant ${density} radius constant ${radius2}
fix pdd1 all particledistribution/discrete 32452843 ${nradii} pts1 ${frac1} pts2 ${frac2}

fix inface all mesh/surface/planar file mesh/insertion_face.stl type 1 scale 0.1 #Define the insertion plane for particles

#particle insertion
fix ins all insert/stream seed 86028157 distributiontemplate pdd1 &
mass ${fillmass} massrate ${fillmassrate} overlapcheck yes all_in no vel constant 0.0 0.0 -1 &
insertion_face inface extrude_length 0.2

set atom 1 charge -29.9792 #Assign charges to all particles (statcoulombs)

#Simulation output and post-processing

shell mkdir post

#Definition of the dumptime

variable dumptime equal 5e-4 #Every dumptime save 1 image
variable dumpstep equal ${dumptime}/${dt} #Convert to iteration steps

run 1 #Insert some particles first so that dump is not empty

dump dumpparticles all custom/vtk ${dumpstep} post/particles_*.vtk id type x y z vx vy vz fx fy fz q radius mass
dump dumpplate all mesh/stl ${dumpstep} post/Plate*.stl plate
dump dumppan all mesh/stl ${dumpstep} post/Pan*.stl pan
#dump dumpelectode all mesh/stl ${dumpstep} post/Electrode*.stl electrode

#Run the simulation filling steps

run ${fillsteps}

#unfix ins #Stop the insertion process

#Run the simulation settling steps

run ${settlesteps}

2. I haven't had much time to explore the second problem yet but I think I will eventually run into it at some points. Wondering if LIGGGHTS has the physics to capture the charge exchange between particles or between particle and surrounding geometries?

I will for sure follow-up on this if I found anything. At the meantime, it is very much appreciated if anyone can share their ideas.

mschramm | Sun, 09/11/2022 - 19:45

For your first problem, I was able to get it working without the set command by adding
I think you only need the q one but I don't have the time to go into the code and find out...
fix m9 all property/atom q scalar yes no no -29.9792
fix m0 all property/atom charge scalar yes no no -29.9792

and adjusting your insert statement
again I think you only need the q one.

#particle insertion
fix ins all insert/stream seed 86028157 distributiontemplate pdd1 &
mass ${fillmass} massrate ${fillmassrate} overlapcheck yes all_in no vel constant 0.0 0.0 -1 &
set_property q -29.9792 set_property charge -29.9792 insertion_face inface extrude_length 0.2

I do not know why you do not get the right VTK output as it is pulling from atom->q but there is no additional force variable just for q as the potential is added directly to f.

q is a purely atom quantity so if you want to make charge a global variable, you will need to add this to the contact physics and extend the skin of each neighbor (for binning) to include the force travel distance.
Or write an extension of the efield class that tracks the centroid of geometries. I may keep thinking about this but there is no ready to go version of this.

zihengw | Mon, 09/12/2022 - 23:19

@mschramm, thanks for your response.

Your approach is very interesting. I've tested it. It didn't raise any errors. But as you mentioned, the q values in the VTK output are all zeros. Of course, the particles didn't respond to the applied electric field at all. (I've tried to increase both the field strength and the charge. No particle changed its trajectory.) I went through the documentation again. The wordings that describe fix property/atom and fix insert/stream set_property gave me a sense that they are just creating a new attribution for the particles. Despite the fact that we named this attribute q or charge, it is not the inherent attribute q or charge which will be accessed by the applied electric field during the electrostatic force calculation.

On the other hand, using the create_atoms plus set command changes the q values in the VTK nicely. And these atoms can be moved by an applied electric field.

I will keep posting my progress. Thanks again for your insights.

mschramm | Tue, 09/13/2022 - 05:43

Hello,
I didn't like my answer so I just coded it myself... If you have access to the source files and can build liggghts again, add these lines in fix_efield.cpp
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
int const * type = atom->type;
// Dirty hack, set q of atoms
int max_type = atom->get_properties()->max_type();
double const * q_val = static_cast(modify->find_fix_property("atomCharge","property/global","peratomtype",max_type,0,"electricField"))->get_values();
fprintf(screen, "max_type = %d\n", max_type);
fprintf(screen, "q = [%f, %f]\n", q_val[0], q_val[1]);
for (int i = 0; i < nlocal; i++) {
q[i] = q_val[type[i]-1];
}
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------

after the
int nlocal = atom->nlocal;
line in the function post_force (should be around line 200).
Then all you need to add to your script is the line
fix m7 all property/global atomCharge peratomtype -29.9792 0.0

I may do some re-edits to tidy things up and make it not so dirty, but that should do what you want (and it is MPI enabled...).

zihengw | Tue, 09/13/2022 - 21:17

@mschramm, I got the following error message when I tried to rebuild LIGGGHTS after adding the lines into fix_efield.cpp. Did you encounter this?

-------------------------------------------------------
../fix_efield.cpp: In member function ‘virtual void LAMMPS_NS::FixEfield::post_force(int)’:
../fix_efield.cpp:204:37: error: expected ‘<’ before ‘(’ token
204 | double const * q_val = static_cast(modify->find_fix_property("atomCharge","property/global"
,"peratomtype",max_type,0,"electricField"))->get_values();
|
../fix_efield.cpp:204:37: error: expected type-specifier before ‘(’ token
../fix_efield.cpp:204:37: error: expected ‘>’ before ‘(’ token
--------------------------------------------------------

Also, if I understand correctly, you are trying to set the atom charge as a global variable. Could you help me see what is the benefit of that? Realistically, I may eventually want to be able to set the atom charge as property/atom. The initial charge of a granular particle depends mostly on its surface area, material property and electric field strength.

mschramm | Thu, 09/15/2022 - 01:22

As to the second part of your post. I see charge as a specific atom type of property and I would treat it as I would with Young's Modulus.
Now, if you wish to have the size of the particle effect its charge, I would change
property/global atomCharge peratomtype ...
to
property/global atomChargeAreaDensity ...
then in the fix_efield.cpp file, I would set q[i] = pi*q_val[type[i]-1]*atom->radius[i]*atom->radius[i]