Dear LAMMPS users and developers
I have a query relating to damping while simulating poly disperse (in size) granular systems.
In ‘pair_style granular’, we have 4 types of damping: velocity, mass_velocity, viscoelastic, tsuji.
(pair_style granular command — LAMMPS documentation)
I will concentrate on the first two (velocity, mass_velocity) in this query.
As per documentation, we (the user) has to specify \eta_{n0}.
In case of damping velocity, \eta_{n0} specified by the user should be in terms of mass/time units, and damping force is evaluated as:
F_{damp} = \eta_{n0} \times v
(v=normal or tangential velocity, as the case may be.)
In case of damping velocity, \eta_{n0} specified by the user should be in terms of 1/time units, and damping force is evaluated as:
F_{damp} = \eta_{n0} \times meff \times v
where meff = effective mass of colliding particles meff=(mi*mj/(mi+mj))
(see lines 363-374 here: lammps/pair_granular.cpp at a4ceda9706e56920f9168ef5987c87d7a327244d · lammps/lammps · GitHub)
Evaluation of \eta_{n0} in both of the above cases is as follows:
Following Silbert(2001),
if we are to use damping velocity, then \eta_{n0}=\frac{2 \ln(e_n) \sqrt{m_{eff}k}}{\sqrt{\pi^2 + \ln^2(e_n)}} (units: mass/time) [Eqn. 1]
and, if we are using damping velocity, then \eta_{n0}=\frac{2 \ln(e_n) \sqrt{k/m_{eff}}}{\sqrt{\pi^2 + \ln^2(e_n)}} (units: 1/time) [Eqn. 2]
(Some more details about [Eqn. 1, 2] are given towards the end in the footnote.)
With the above expressions of \eta_{n0}, we shall have \eta_{n0} in appropriate units required/expected by LAMMPS as given in the documentation.
(see also: [Update documentation] pair_style granular command damping mass velocity · Issue #3016 · lammps/lammps · GitHub)
My concern is that the term meff appears in the expression used to evaluate \eta_{n0}.
Hence, we have an issue that these damping styles will work only for a mono-disperse system.
In a poly-disperse system (eg. 10000 particle system whose dia is uniformly distributed in a range, say 0.8 to 1.2), the meff of the colliding particles is not known apriori before the simulation starts.
Particles of any two sizes can possibly collide. In such a scenario, we do not know meff of colliding particles and as a consequence, corresponding \eta_{n0} cannot be evaluated and specified in the input script by the user.
Kindly confirm if I am correct so far.
Please correct me if I am missing something in here. I’d be very grateful.
If my understanding detailed above is indeed correct, the code needs to be modified to take care of poly-disperse collisions.
This can be achieved in the following manner:
- meff term should be pulled out of \eta_{n0} in ‘damping velocity’ and/or ‘damping mass_velocity’.
In such a case
F_{damp} = \eta_{n0} \times sqrt(meff) \times v
where, user specified \eta_{n0} = \frac{2 \ln(e_n) \sqrt{k}}{\sqrt{\pi^2 + \ln^2(e_n)}}
(units: mass^0.5/time)
Or an entirely new damping can be included if we wish to ensure that the old input scripts of other users are not affected by these changes in newer versions of LAMMPS.
If that is the case, we can in fact go one step further and ask the user to just specify normal(and tangential) coefficient of restitution (e_n and e_t) since this is more physical.
Then, F_{damp} can be internally calculated as
F_{damp} = \frac{2 \ln(e_n) \sqrt{m_{eff}k}}{\sqrt{\pi^2 + \ln^2(e_n)}} \times v
We have to ensure that the above changes are made for:
- both the normal and tangential damping cases for the case of particle-particle collision. (src/GRANULAR/pair_granular.cpp file)
- both the normal and tangential damping cases for the case of particle-wall collision. (src/GRANULAR/fix_wall_gran.cpp)
Footnote:
Q: Why should we believe your expressions for \eta_{n0} [Eqn. 1,2].
A: I have obtained this by solving the (linear spring,damper,mass) system (Silbert. 2001).
Furthermore, I have checked and confirmed LAMMPS behaviour by running a few test simulations. I have considered the case of two particle collision in which one particle is initially at rest (v=0) and other impacts it with velocity (v=1). We can obtain final velocities of the particles and evaluate e=(v2-v1)_{final}/(v2-v1)_{initial} = (v2-v1)_{final}/1 = v2-v1
Case 1: Damping mass velocity.
Choose any kn value.
Choose any e value.
Choose dia of the colliding particles
Choose density of the particles.
Calculate the masses of two particles and meff.
Calculate \eta_{n0} value as per [Eqn. 1]
Specify this \eta_{n0} value in the input script.
Run the simulation.
Obtain final velocities of the particles.
Calculate e=(v2-v1)_{final}. This should be the same as that of our chosen e, thus confirming the correctness of [Eqn. 1] to evaluate \eta_{n0} for damping velocity.
The same can be carried out with damping mass velocity.
Comment/uncomment appropriate line in the input script.
Use [Eqn. 2] to evaluate \eta_{n0} for this case.
I am copying below three files: (new user-can’t upload files)
- in.dampingcheck.txt: input file
- grains.txt : read_data file
- etan0_calculator.py to calculate {eta_n0} value.
#input script begin
units si
dimension 2
atom_style sphere
boundary p p p
newton off
comm_modify vel yes
atom_modify map yes
read_data grains.txt
neighbor 0.2 bin #0.002 bin
neigh_modify delay 0
variable knpp equal 100000
variable ktpp equal (2/7)*${knpp}
variable xgammat equal 1 #tgt_damp = xgammat * normal_damp
variable mupp equal 0.5 #coeff. of friction between particle-particle
#gammapp = eta_n0; can use etan0_calculator.py to get this value.
pair_style granular
#comment/uncomment one of the two pair_coeff commands. use appropriate gammanpp value.
#variable gammanpp equal 64.025 #64.025-->en=0.5 damping velocity
#pair_coeff * * hooke ${knpp} ${gammanpp} tangential linear_history ${ktpp} ${xgammat} ${mupp} damping velocity #limit_damping
variable gammanpp equal 290.0136 #290.0136-->en=0.5 damping mass_velocity
pair_coeff * * hooke ${knpp} ${gammanpp} tangential linear_history ${ktpp} ${xgammat} ${mupp} damping mass_velocity #limit_damping
fix 1 all nve/sphere #disc
group id1 id 1
group id2 id 2
group id3 id 3
compute 1 id1 property/atom vy
compute v1 id1 reduce sum c_1
compute 2 id2 property/atom vy
compute v2 id2 reduce sum c_2
thermo 10000
thermo_style custom step time c_v1 c_v2
variable dumpfreq equal 1000000
dump depositions all custom ${dumpfreq} test2.dump.* id x y z vx vy vz type diameter radius fx fy fz mass omegaz
fix 2dfix all enforce2d
timestep 0.000001
run 5000000
print "Simulation Complete."
#input script end
#read_data file grains.txt
atom config
2 atoms
10 atom types
0.0 10.0 xlo xhi
0.0 20.0 ylo yhi
-0.5 0.5 zlo zhi
Atoms
1 6 1.0 1.0 3.0 11.2 0.0
2 6 0.9 1.0 3.0 9.1 0.0
Velocities
1 -0.0 -1.0 0.0 0 0 0
2 0.0 0.0 -0.0 0 0 0
#end of read_data file
#.py script to calculate eta_n0 values
#etan0calculator.py
import numpy as np
pi = np.pi
dia1 = 1.0
dia2 = 0.9 #1.0
density = 1.0
en = 0.5 # normal coeff. of restitution
kn = 10**5 # hookean spring
m1 = density * (4/3) * pi * ((dia1/2)**3)
m2 = density * (4/3) * pi * ((dia2/2)**3)
meff = m1*m2/(m1+m2)
eta_n0_vel = -(2 * np.log(en) * np.sqrt(meff*kn))/np.sqrt(pi**2 + np.log(en)**2)
eta_n0_massvel = -(2 * np.log(en) * np.sqrt(kn/meff))/np.sqrt(pi**2 + np.log(en)**2)
print('eta_n0_vel = ', eta_n0_vel)
print('eta_n0_massvel = ', eta_n0_massvel)
#end of .py script
Thank you.
Vikas K
IIT Kanpur
India.