Merge pull request #4409 from willzunker/mdr-rebase2
pair_style granular - MDR contact model
This commit is contained in:
Axel Kohlmeyer
@ -222,10 +222,10 @@ restart file, so that the operation of the fix continues in an
|
||||
uninterrupted fashion.
|
||||
|
||||
If the :code:`contacts` option is used, this fix generates a per-atom array
|
||||
with 8 columns as output, containing the contact information for owned
|
||||
with at least 8 columns as output, containing the contact information for owned
|
||||
particles (nlocal on each processor). All columns in this per-atom array will
|
||||
be zero if no contact has occurred. The values of these columns are listed in
|
||||
the following table:
|
||||
be zero if no contact has occurred. The first 8 values of these columns are
|
||||
listed in the following table.
|
||||
|
||||
+-------+----------------------------------------------------+----------------+
|
||||
| Index | Value | Units |
|
||||
@ -248,6 +248,14 @@ the following table:
|
||||
| 8 | Radius :math:`r` of atom | distance units |
|
||||
+-------+----------------------------------------------------+----------------+
|
||||
|
||||
If a granular submodel calculates additional contact information (e.g. the
|
||||
heat submodels calculate the amount of heat exchanged), these quantities
|
||||
are appended to the end of this array. First, any extra values from the
|
||||
normal submodel are appended followed by the damping, tangential, rolling,
|
||||
twisting, then heat models. See the descriptions of granular submodels in
|
||||
the :doc:`pair granular <pair_granular>` page for information on any extra
|
||||
quantities.
|
||||
|
||||
None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
|
||||
No parameter of this fix can be used with the *start/stop* keywords of the
|
||||
:doc:`run <run>` command. This fix is not invoked during :doc:`energy
|
||||
|
||||
@ -243,10 +243,10 @@ uninterrupted fashion.
|
||||
with a different region ID.
|
||||
|
||||
If the :code:`contacts` option is used, this fix generates a per-atom array
|
||||
with 8 columns as output, containing the contact information for owned
|
||||
with at least 8 columns as output, containing the contact information for owned
|
||||
particles (nlocal on each processor). All columns in this per-atom array will
|
||||
be zero if no contact has occurred. The values of these columns are listed in
|
||||
the following table:
|
||||
be zero if no contact has occurred. The first 8 values of these columns are
|
||||
listed in the following table.
|
||||
|
||||
+-------+----------------------------------------------------+----------------+
|
||||
| Index | Value | Units |
|
||||
@ -269,6 +269,14 @@ the following table:
|
||||
| 8 | Radius :math:`r` of atom | distance units |
|
||||
+-------+----------------------------------------------------+----------------+
|
||||
|
||||
If a granular submodel calculates additional contact information (e.g. the
|
||||
heat submodels calculate the amount of heat exchanged), these quantities
|
||||
are appended to the end of this array. First, any extra values from the
|
||||
normal submodel are appended followed by the damping, tangential, rolling,
|
||||
twisting, then heat models. See the descriptions of granular submodels in
|
||||
the :doc:`pair granular <pair_granular>` page for information on any extra
|
||||
quantities.
|
||||
|
||||
None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
|
||||
No parameter of this fix can be used with the *start/stop* keywords of the
|
||||
:doc:`run <run>` command. This fix is not invoked during :doc:`energy
|
||||
|
||||
@ -40,6 +40,9 @@ Examples
|
||||
pair_style granular
|
||||
pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat area 0.1
|
||||
|
||||
pair_style granular
|
||||
pair_coeff * * mdr 5e6 0.4 1.9e5 2.0 0.5 0.5 tangential linear_history 940.0 0.0 0.7 rolling sds 2.7e5 0.0 0.6 damping none
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
@ -82,6 +85,7 @@ and their required arguments are:
|
||||
3. *hertz/material* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`
|
||||
4. *dmt* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
|
||||
5. *jkr* : E, :math:`\eta_{n0}` (or :math:`e`), :math:`\nu`, :math:`\gamma`
|
||||
6. *mdr* : :math:`E`, :math:`\nu`, :math:`Y`, :math:`\Delta\gamma`, :math:`\psi_b`, :math:`e`
|
||||
|
||||
Here, :math:`k_n` is spring stiffness (with units that depend on model
|
||||
choice, see below); :math:`\eta_{n0}` is a damping prefactor (or, in its
|
||||
@ -162,6 +166,143 @@ initially will not experience force until they come into contact
|
||||
experience a tensile force up to :math:`3\pi\gamma R`, at which point they
|
||||
lose contact.
|
||||
|
||||
The *mdr* model is a mechanically-derived contact model designed to capture the
|
||||
contact response between adhesive elastic-plastic particles into large deformation.
|
||||
The theoretical foundations of the *mdr* model are detailed in the
|
||||
two-part series :ref:`Zunker and Kamrin Part I <Zunker2024I>` and
|
||||
:ref:`Zunker and Kamrin Part II <Zunker2024II>`. Further development
|
||||
and demonstrations of its application to industrially relevant powder
|
||||
compaction processes are presented in :ref:`Zunker et al. <Zunker2025>`.
|
||||
|
||||
The model requires the following inputs:
|
||||
|
||||
1. *Young's modulus* :math:`E > 0` : The Young's modulus is commonly reported
|
||||
for various powders.
|
||||
|
||||
2. *Poisson's ratio* :math:`0 \le \nu \le 0.5` : The Poisson's ratio is commonly
|
||||
reported for various powders.
|
||||
|
||||
3. *Yield stress* :math:`Y \ge 0` : The yield stress is often known for powders
|
||||
composed of materials such as metals but may be unreported for ductile organic
|
||||
materials, in which case it can be treated as a free parameter.
|
||||
|
||||
4. *Effective surface energy* :math:`\Delta\gamma \ge 0` : The effective surface
|
||||
energy for powder compaction applications is most easily determined through its
|
||||
relation to the more commonly reported critical stress intensity factor
|
||||
:math:`K_{Ic} = \sqrt{2\Delta\gamma E/(1-\nu^2)}`.
|
||||
|
||||
5. *Critical confinement ratio* :math:`0 \le \psi_b \le 1` : The critical confinement
|
||||
ratio is a tunable parameter that determines when the bulk elastic response is
|
||||
triggered. Lower values of :math:`\psi_b` delay the onset of the bulk elastic
|
||||
response.
|
||||
|
||||
6. *Coefficient of restitution* :math:`0 \le e \le 1` : The coefficient of
|
||||
restitution is a tunable parameter that controls damping in the normal direction.
|
||||
|
||||
.. note::
|
||||
|
||||
The values for :math:`E`, :math:`\nu`, :math:`Y`, and :math:`\Delta\gamma` (i.e.,
|
||||
:math:`K_{Ic}`) should be selected for zero porosity to reflect the intrinsic
|
||||
material property rather than the bulk powder property.
|
||||
|
||||
The *mdr* model produces a nonlinear force-displacement response, therefore the
|
||||
critical timestep :math:`\Delta t` depends on the inputs and level of
|
||||
deformation. As a conservative starting point the timestep can be assumed to be
|
||||
dictated by the bulk elastic response such that
|
||||
:math:`\Delta t = 0.35\sqrt{m/k_\textrm{bulk}}`, where :math:`m` is the mass of
|
||||
the smallest particle and :math:`k_\textrm{bulk} = \kappa R_\textrm{min}` is an
|
||||
effective stiffness related to the bulk elastic response.
|
||||
Here, :math:`\kappa = E/(3(1-2\nu))` is the bulk modulus and
|
||||
:math:`R_\textrm{min}` is the radius of the smallest particle.
|
||||
|
||||
.. note::
|
||||
|
||||
The *mdr* model requires some specific settings to function properly,
|
||||
please read the following text carefully to ensure all requirements are
|
||||
followed.
|
||||
|
||||
The *atom_style* must be set to *sphere 1* to enable dynamic particle
|
||||
radii. The *mdr* model is designed to respect the incompressibility of
|
||||
plastic deformation and inherently tracks free surface displacements
|
||||
induced by all particle contacts. In practice, this means that all particles
|
||||
begin with an initial radius, however as compaction occurs and plastic
|
||||
deformation is accumulated, a new enlarged apparent radius is defined to
|
||||
ensure that that volume change due to plastic deformation is not lost.
|
||||
This apparent radius is stored as the *atom radius* meaning it is used
|
||||
for subsequent neighbor list builds and contact detection checks. The
|
||||
advantage of this is that multi-neighbor dependent effects such as
|
||||
formation of secondary contacts caused by radial expansion are captured
|
||||
by the *mdr* model. Setting *atom_style sphere 1* ensures that updates to
|
||||
the particle radii are properly reflected throughout the simulation.
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
atom_style sphere 1
|
||||
|
||||
Newton's third law must be set to *off*. This ensures that the neighbor lists
|
||||
are constructed properly for the topological penalty algorithm used to screen
|
||||
for non-physical contacts occurring through obstructing particles, an issue
|
||||
prevalent under large deformation conditions. For more information on this
|
||||
algorithm see :ref:`Zunker et al. <Zunker2025>`.
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
newton off
|
||||
|
||||
The damping model must be set to *none*. The *mdr* model already has a built
|
||||
in damping model.
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_coeff * * mdr 5e6 0.4 1.9e5 2 0.5 0.5 damping none
|
||||
|
||||
The definition of multiple *mdr* models in the *pair_style* is currently not
|
||||
supported. Similarly, the *mdr* model cannot be combined with a different normal
|
||||
model in the *pair_style*. Physically this means that only one homogenous
|
||||
collection of particles governed by a single *mdr* model is allowed.
|
||||
|
||||
The *mdr* model currently only supports *fix wall/gran/region*, not
|
||||
*fix wall/gran*. If the *mdr* model is specified for the *pair_style*
|
||||
any *fix wall/gran/region* commands must also use the *mdr* model.
|
||||
Additionally, the following *mdr* inputs must match between the
|
||||
*pair_style* and *fix wall/gran/region* definitions: :math:`E`,
|
||||
:math:`\nu`, :math:`Y`, :math:`\psi_b`, and :math:`e`. The exception
|
||||
is :math:`\Delta\gamma`, which may vary, permitting different
|
||||
adhesive behaviors between particle-particle and particle-wall interactions.
|
||||
|
||||
.. note::
|
||||
|
||||
The *mdr* model has a number of custom *property/atom* and *pair/local* definitions that
|
||||
can be called in the input file. The useful properties for visualization
|
||||
and analysis are described below.
|
||||
|
||||
In addition to contact forces the *mdr* model also tracks the following
|
||||
quantities for each particle: elastic volume change, average normal
|
||||
stress components, total surface area involved in
|
||||
contact, and individual contact areas. In the input script, these quantities are
|
||||
initialized by calling *run 0* and can then be accessed using subsequent *compute*
|
||||
commands. The last *compute* command uses *pair/local p13* to calculate the pairwise
|
||||
contact areas for each active contact in the *group-ID*. Due to the use of an apparent
|
||||
radius in the *mdr* model, the keyword/arg pair *cutoff radius* must be specified for
|
||||
*pair/local* to properly detect existing contacts.
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
run 0
|
||||
compute ID group-ID property/atom d_Velas
|
||||
compute ID group-ID property/atom d_sigmaxx
|
||||
compute ID group-ID property/atom d_sigmayy
|
||||
compute ID group-ID property/atom d_sigmazz
|
||||
compute ID group-ID property/atom d_Acon1
|
||||
compute ID group-ID pair/local p13 cutoff radius
|
||||
|
||||
.. note::
|
||||
|
||||
The *mdr* model has two example input scripts within the
|
||||
*examples/granular* directory. The first is a die compaction
|
||||
simulation involving 200 particles named *in.tableting.200*.
|
||||
The second is a triaxial compaction simulation involving 12
|
||||
particles named *in.triaxial.compaction.12*.
|
||||
----------
|
||||
|
||||
In addition, the normal force is augmented by a damping term of the
|
||||
@ -674,7 +815,10 @@ supported are:
|
||||
2. *radius* : :math:`k_{s}`
|
||||
3. *area* : :math:`h_{s}`
|
||||
|
||||
If the *heat* keyword is not specified, the model defaults to *none*.
|
||||
If the *heat* keyword is not specified, the model defaults to *none*. All
|
||||
heat models calculate an additional pairwise quantity accessible by the
|
||||
single() function (described below) which is the heat conducted between the
|
||||
two particles.
|
||||
|
||||
For *heat* *radius*, the heat
|
||||
:math:`Q` conducted between two particles is given by
|
||||
@ -789,7 +933,7 @@ The single() function of these pair styles returns 0.0 for the energy
|
||||
of a pairwise interaction, since energy is not conserved in these
|
||||
dissipative potentials. It also returns only the normal component of
|
||||
the pairwise interaction force. However, the single() function also
|
||||
calculates 13 extra pairwise quantities. The first 3 are the
|
||||
calculates at least 13 extra pairwise quantities. The first 3 are the
|
||||
components of the tangential force between particles I and J, acting
|
||||
on particle I. The fourth is the magnitude of this tangential force.
|
||||
The next 3 (5-7) are the components of the rolling torque acting on
|
||||
@ -797,9 +941,17 @@ particle I. The next entry (8) is the magnitude of the rolling torque.
|
||||
The next entry (9) is the magnitude of the twisting torque acting
|
||||
about the vector connecting the two particle centers.
|
||||
The next 3 (10-12) are the components of the vector connecting
|
||||
the centers of the two particles (x_I - x_J). The last quantity (13)
|
||||
is the heat flow between the two particles, set to 0 if no heat model
|
||||
is active.
|
||||
the centers of the two particles (x_I - x_J). If a granular submodel
|
||||
calculates additional contact information (e.g. the heat submodels
|
||||
calculate the amount of heat exchanged), these quantities are appended
|
||||
to the end of this list. First, any extra values from the normal submodel
|
||||
are appended followed by the damping, tangential, rolling, twisting, then
|
||||
heat models. See the descriptions of specific granular submodels above
|
||||
for information on any extra quantities. If two or more models are
|
||||
defined by pair coefficients, the size of the array is set by the
|
||||
maximum number of extra quantities in a model but the order of quantities
|
||||
is determined by each model's specific set of submodels. Any unused
|
||||
quantities are zeroed.
|
||||
|
||||
These extra quantities can be accessed by the :doc:`compute pair/local <compute_pair_local>` command, as *p1*, *p2*, ...,
|
||||
*p12*\ .
|
||||
@ -870,10 +1022,32 @@ solids. Proc. R. Soc. Lond. A, 324(1558), 301-313.
|
||||
|
||||
.. _DMT1975:
|
||||
|
||||
**Derjaguin et al, 1975)** Derjaguin, B. V., Muller, V. M., & Toporov,
|
||||
**(Derjaguin et al, 1975)** Derjaguin, B. V., Muller, V. M., & Toporov,
|
||||
Y. P. (1975). Effect of contact deformations on the adhesion of
|
||||
particles. Journal of Colloid and interface science, 53(2), 314-326.
|
||||
|
||||
.. _Zunker2024I:
|
||||
|
||||
**(Zunker and Kamrin, 2024)** Zunker, W., & Kamrin, K. (2024).
|
||||
A mechanically-derived contact model for adhesive elastic-perfectly
|
||||
plastic particles, Part I: Utilizing the method of dimensionality
|
||||
reduction. Journal of the Mechanics and Physics of Solids, 183, 105492.
|
||||
|
||||
.. _Zunker2024II:
|
||||
|
||||
**(Zunker and Kamrin, 2024)** Zunker, W., & Kamrin, K. (2024).
|
||||
A mechanically-derived contact model for adhesive elastic-perfectly
|
||||
plastic particles, Part II: Contact under high compaction—modeling
|
||||
a bulk elastic response. Journal of the Mechanics and Physics of Solids,
|
||||
183, 105493.
|
||||
|
||||
.. _Zunker2025:
|
||||
|
||||
**(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S.,
|
||||
Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large
|
||||
deformation powder compaction: mechanically-derived contact model and
|
||||
screening of non-physical contacts. engrXiv.
|
||||
|
||||
.. _Luding2008:
|
||||
|
||||
**(Luding, 2008)** Luding, S. (2008). Cohesive, frictional powders:
|
||||
|
||||
@ -56,7 +56,7 @@ Syntax
|
||||
random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ternary(x,y,z),
|
||||
ramp(x,y), stagger(x,y), logfreq(x,y,z), logfreq2(x,y,z),
|
||||
logfreq3(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c),
|
||||
vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
|
||||
vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z), sign(x)
|
||||
group functions = count(group), mass(group), charge(group),
|
||||
xcm(group,dim), vcm(group,dim), fcm(group,dim),
|
||||
bound(group,dir), gyration(group), ke(group),
|
||||
@ -532,37 +532,37 @@ functions, special functions, feature functions, atom values, atom
|
||||
vectors, custom atom properties, compute references, fix references, and references to other
|
||||
variables.
|
||||
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Number | 0.2, 100, 1.0e20, -15.4, etc |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Constant | PI, version, on, off, true, false, yes, no |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Thermo keywords | vol, pe, ebond, etc |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Math operators | (), -x, x+y, x-y, x\*y, x/y, x\^y, x%y, x == y, x != y, x < y, x <= y, x > y, x >= y, x && y, x \|\| y, x \|\^ y, !x |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Math functions | sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ternary(x,y,z), ramp(x,y), stagger(x,y), logfreq(x,y,z), logfreq2(x,y,z), logfreq3(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z) |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Math functions | sqrt(x), exp(x), ln(x), log(x), abs(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ternary(x,y,z), ramp(x,y), stagger(x,y), logfreq(x,y,z), logfreq2(x,y,z), logfreq3(x,y,z), stride(x,y,z), stride2(x,y,z,a,b,c), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z), sign(x) |
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Group functions | count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), inertia(ID,dimdim), omega(ID,dim) |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Region functions | count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR) |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Special functions | sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), rsort(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label), is_timeout() |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Feature functions | is_available(category,feature), is_active(category,feature), is_defined(category,id) |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Atom values | id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i] |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Atom vectors | id, mass, type, mol, x, y, z, vx, vy, vz, fx, fy, fz, q |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Custom atom properties | i_name, d_name, i_name[i], d_name[i], i2_name[i], d2_name[i], i2_name[i][j], d_name[i][j] |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Compute references | c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i] |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Fix references | f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i] |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
| Other variables | v_name, v_name[i] |
|
||||
+------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
+------------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|
||||
|
||||
Most of the formula elements produce a scalar value. Some produce a
|
||||
global or per-atom vector of values. Global vectors can be produced
|
||||
@ -860,6 +860,9 @@ run, according to one of these formulas, where omega = 2 PI / period:
|
||||
|
||||
where dt = the timestep size.
|
||||
|
||||
The sign(x) function returns 1.0 if the value is greater than or equal
|
||||
to 0.0, and -1.0 otherwise.
|
||||
|
||||
The run begins on startstep. Startstep can span multiple runs, using
|
||||
the *start* keyword of the :doc:`run <run>` command. See the :doc:`run
|
||||
<run>` command for details of how to do this. Note that the
|
||||
|
||||
149
examples/granular/in.tableting.200
Normal file
149
examples/granular/in.tableting.200
Normal file
@ -0,0 +1,149 @@
|
||||
##################################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere 1
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 1.0e-3 bin
|
||||
neigh_modify every 10 delay 60 check no
|
||||
timestep 4e-6
|
||||
#processors 2 2 1
|
||||
|
||||
############################## SIMULATION BOUNDING BOX AND INSERT PARTICLES #################################
|
||||
|
||||
boundary f f f
|
||||
read_data spheres200.data
|
||||
|
||||
#################################### ADD DIE AND ATOM PARAMETERIZATION ######################################
|
||||
|
||||
variable atomRadius equal 0.44e-3*1.25
|
||||
variable atomDiameter equal 2*${atomRadius}
|
||||
variable atomDensity equal 1560
|
||||
variable atomMassAvg equal ${atomDensity}*4.0/3.0*PI*${atomRadius}^3.0
|
||||
variable dieRadius equal 4e-3
|
||||
variable dieHeight equal 1e-2
|
||||
|
||||
############################## PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ##################################
|
||||
|
||||
pair_style granular
|
||||
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
variable YoungsModulus equal 5e6
|
||||
variable YieldStress equal 1.9e5
|
||||
variable PoissonsRatio equal 0.4
|
||||
variable SurfaceEnergy equal 2
|
||||
variable SurfaceEnergyWall equal 0.0
|
||||
variable CoR equal 0.5
|
||||
variable psi_b equal 0.5
|
||||
|
||||
# linear_history = k_t, x_gammat, mu_s
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable kt_wall equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0
|
||||
variable mu_s equal 0.7
|
||||
variable mu_s_wall equal 0.1
|
||||
|
||||
# sds = mu_roll, k_roll, gamma_roll
|
||||
variable mu_roll equal 0.6
|
||||
variable k_roll equal 2.25*${mu_roll}*${mu_roll}*${YoungsModulus}*${atomRadius}
|
||||
variable gamma_roll equal 0.0
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} rolling sds ${k_roll} ${gamma_roll} ${mu_roll} damping none
|
||||
|
||||
######################################### ADD DIE AND PUNCH WALLS ############################################
|
||||
|
||||
variable disp_upper equal 0.0
|
||||
variable disp_lower equal 0.0
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergyWall} ${psi_b} ${CoR} tangential linear_history ${kt_wall} ${xgammat} ${mu_s_wall} rolling sds ${k_roll} ${gamma_roll} ${mu_roll} damping none"
|
||||
|
||||
variable dieHeight2 equal 2*${dieHeight}
|
||||
|
||||
region lowerPunch plane 0 0 0 0 0 1 side in units box move NULL NULL v_disp_lower units box
|
||||
region upperPunch plane 0 0 ${dieHeight} 0 0 -1 side in move NULL NULL v_disp_upper units box
|
||||
region die cylinder z 0 0 ${dieRadius} 0 ${dieHeight2} side in units box
|
||||
|
||||
fix lowerPunch all wall/gran/region ${wall_contact_string} region lowerPunch contacts
|
||||
fix upperPunch all wall/gran/region ${wall_contact_string} region upperPunch contacts
|
||||
fix die all wall/gran/region ${wall_contact_string} region die contacts
|
||||
|
||||
compute avgUpperPunchForce all reduce sum f_upperPunch[4]
|
||||
variable avgUpperPunchForce equal c_avgUpperPunchForce
|
||||
compute avgLowerPunchForce all reduce sum f_lowerPunch[4]
|
||||
variable avgLowerPunchForce equal c_avgLowerPunchForce
|
||||
|
||||
fix printFD all print 1 "${disp_upper} ${avgUpperPunchForce} ${avgLowerPunchForce}" file punch_force_disp_tableting200.csv screen no
|
||||
|
||||
##################################### INTEGRATION AND GRAVITY #################################################
|
||||
|
||||
fix 1 all nve/sphere
|
||||
fix grav all gravity 9.81 vector 0 0 -1
|
||||
|
||||
########################################### SCREEN OUTPUT ####################################################
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke vol v_disp_upper
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
compute ke all ke/atom
|
||||
variable output_rate equal round(1e-3/dt)
|
||||
|
||||
run 0
|
||||
|
||||
compute sigmaxx all property/atom d_sigmaxx
|
||||
compute sigmayy all property/atom d_sigmayy
|
||||
compute sigmazz all property/atom d_sigmazz
|
||||
compute Velas all property/atom d_Velas
|
||||
|
||||
compute sigmaxx_ave all reduce ave c_sigmaxx
|
||||
compute sigmayy_ave all reduce ave c_sigmayy
|
||||
compute sigmazz_ave all reduce ave c_sigmazz
|
||||
compute Velas_sum all reduce sum c_Velas
|
||||
|
||||
variable sxx_ave equal c_sigmaxx_ave
|
||||
variable syy_ave equal c_sigmayy_ave
|
||||
variable szz_ave equal c_sigmazz_ave
|
||||
variable Vparticles equal c_Velas_sum
|
||||
|
||||
fix log all print 1 "${sxx_ave} ${syy_ave} ${szz_ave} ${Vparticles}" file average_normal_stresses_tableting200.csv screen no
|
||||
dump dumpParticles all custom ${output_rate} tableting200.dump id type mass diameter x y z vx vy vz fx fy fz c_ke c_sigmaxx c_sigmayy c_sigmazz
|
||||
#dump dumpParticlesVTK all vtk ${output_rate} post/particles_*.vtk id x y z fx fy fz vx vy vz c_ke radius c_sigmaxx c_sigmayy c_sigmazz
|
||||
|
||||
############################################## RUN SIMULATION #################################################
|
||||
|
||||
variable upper_punch_stroke equal 0.6733*${dieHeight}
|
||||
variable vel_upper equal 0.25
|
||||
|
||||
variable settling_steps equal round(0.02/dt)
|
||||
variable compression_steps equal 2*round(${upper_punch_stroke}/${vel_upper}/dt)
|
||||
variable ejection_steps equal ${compression_steps}
|
||||
variable free_float_steps equal round(0.02/dt)
|
||||
|
||||
##### SETTLING #####
|
||||
|
||||
run ${settling_steps}
|
||||
|
||||
##### Compression & Release #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${compression_steps}/2)
|
||||
variable disp_upper equal -${upper_punch_stroke}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable short_release equal round(${compression_steps}*1.0)
|
||||
run ${short_release}
|
||||
|
||||
##### EJECTION #####
|
||||
|
||||
variable punch_frequency equal PI/2/(dt*${ejection_steps})
|
||||
variable disp_lower equal ${dieHeight}*sin(${punch_frequency}*elapsed*dt)
|
||||
variable disp_upper equal 0.9*v_disp_lower
|
||||
run ${ejection_steps}
|
||||
|
||||
##### FREE FLOAT #####
|
||||
|
||||
variable disp_lower equal ${dieHeight}
|
||||
variable disp_upper equal ${dieHeight}*0.9
|
||||
variable max_disp equal ${dieRadius}*0.75
|
||||
run ${free_float_steps}
|
||||
109
examples/granular/in.triaxial.compaction.12
Normal file
109
examples/granular/in.triaxial.compaction.12
Normal file
@ -0,0 +1,109 @@
|
||||
############################### SIMULATION SETTINGS ###################################################
|
||||
|
||||
atom_style sphere 1
|
||||
atom_modify map array
|
||||
comm_modify vel yes
|
||||
units si
|
||||
newton off
|
||||
neighbor 2 bin
|
||||
neigh_modify delay 0
|
||||
timestep 1e-6
|
||||
|
||||
##################### SIMULATION BOUNDING BOX, INSERT PARTICLES, AND INTEGRATION #######################
|
||||
|
||||
boundary f f f
|
||||
read_data spheres12.data
|
||||
fix integr all nve/sphere
|
||||
|
||||
# create pair group for contact area outputs
|
||||
group particles_1_12 id 1 12
|
||||
|
||||
########################### PARTICLE MATERIAL PROPERTIES AND FORCE MODEL ###############################
|
||||
|
||||
variable atomRadius equal 0.5
|
||||
|
||||
pair_style granular
|
||||
|
||||
# mdr = E, nu, Y, gamma, psi_b, CoR
|
||||
variable YoungsModulus equal 1e9
|
||||
variable PoissonsRatio equal 0.3
|
||||
variable YieldStress equal 50e6
|
||||
variable SurfaceEnergy equal 0.0
|
||||
variable psi_b equal 0.5
|
||||
variable CoR equal 0.5
|
||||
|
||||
# linear_history = k_t, x_gamma,t, mu_s
|
||||
variable kt equal 2/7*${YoungsModulus}*${atomRadius}
|
||||
variable xgammat equal 0.0
|
||||
variable mu_s equal 0.5
|
||||
|
||||
pair_coeff * * mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none
|
||||
|
||||
######################################### ADD IN PLANES ################################################
|
||||
|
||||
variable boxWidth equal 3
|
||||
variable halfBoxWidth equal ${boxWidth}/2
|
||||
|
||||
variable plane_disp equal 0.0
|
||||
variable plane_disp_neg equal 0.0
|
||||
|
||||
region plane_yz_pos plane ${halfBoxWidth} 0 0 -1 0 0 side in move v_plane_disp_neg NULL NULL units box
|
||||
region plane_yz_neg plane -${halfBoxWidth} 0 0 1 0 0 side in move v_plane_disp NULL NULL units box
|
||||
region plane_xz_pos plane 0 ${halfBoxWidth} 0 0 -1 0 side in move NULL v_plane_disp_neg NULL units box
|
||||
region plane_xz_neg plane 0 -${halfBoxWidth} 0 0 1 0 side in move NULL v_plane_disp NULL units box
|
||||
region plane_xy_pos plane 0 0 ${halfBoxWidth} 0 0 -1 side in move NULL NULL v_plane_disp_neg units box
|
||||
region plane_xy_neg plane 0 0 -${halfBoxWidth} 0 0 1 side in move NULL NULL v_plane_disp units box
|
||||
|
||||
variable wall_contact_string string "granular mdr ${YoungsModulus} ${PoissonsRatio} ${YieldStress} ${SurfaceEnergy} ${psi_b} ${CoR} tangential linear_history ${kt} ${xgammat} ${mu_s} damping none"
|
||||
|
||||
fix plane_yz_pos all wall/gran/region ${wall_contact_string} region plane_yz_pos contacts
|
||||
fix plane_yz_neg all wall/gran/region ${wall_contact_string} region plane_yz_neg contacts
|
||||
fix plane_xz_pos all wall/gran/region ${wall_contact_string} region plane_xz_pos contacts
|
||||
fix plane_xz_neg all wall/gran/region ${wall_contact_string} region plane_xz_neg contacts
|
||||
fix plane_xy_pos all wall/gran/region ${wall_contact_string} region plane_xy_pos contacts
|
||||
fix plane_xy_neg all wall/gran/region ${wall_contact_string} region plane_xy_neg contacts
|
||||
|
||||
compute plane_xy_neg_force all reduce sum f_plane_xy_neg[4]
|
||||
variable plane_xy_neg_force equal c_plane_xy_neg_force
|
||||
|
||||
compute plane_xz_neg_force all reduce sum f_plane_xz_neg[3]
|
||||
variable plane_xz_neg_force equal c_plane_xz_neg_force
|
||||
|
||||
compute plane_yz_neg_force all reduce sum f_plane_yz_neg[2]
|
||||
variable plane_yz_neg_force equal c_plane_yz_neg_force
|
||||
|
||||
fix print1 all print 1 "${plane_disp} ${plane_xy_neg_force} ${plane_xz_neg_force} ${plane_yz_neg_force}" file force_disp_triaxial12.csv screen no
|
||||
|
||||
######################################## SCREEN OUTPUT ####################################################
|
||||
|
||||
compute 1 all erotate/sphere
|
||||
thermo_style custom dt step atoms ke c_1 vol
|
||||
thermo 100
|
||||
thermo_modify lost ignore norm no
|
||||
|
||||
##################################### DEFINE WALL MOVEMENT #################################################
|
||||
|
||||
variable disp_max equal 0.499
|
||||
variable ddisp equal 0.00001
|
||||
variable compression_steps equal round(${disp_max}/${ddisp})
|
||||
variable output_rate equal round(${compression_steps}/100)
|
||||
|
||||
##################################### SET UP DUMP OUTPUTS ####################################################
|
||||
|
||||
dump dumpParticles all custom ${output_rate} triaxial_compaction_12.dump id type mass x y z vx vy vz fx fy fz radius
|
||||
#dump dmp all vtk ${output_rate} post/triaxial12particles_*.vtk id type mass x y z vx vy vz fx fy fz radius
|
||||
|
||||
#################################### COMPRESS THE PARTICLES ##################################################
|
||||
|
||||
run 0
|
||||
|
||||
# print out contact area evolution for particles 1 and 12
|
||||
compute Ac_1_12 particles_1_12 pair/local p13 cutoff radius
|
||||
compute Ac_1_12_sum particles_1_12 reduce sum c_Ac_1_12 inputs local
|
||||
variable Ac_1_12 equal c_Ac_1_12_sum
|
||||
fix logArea all print 100 "${plane_disp} ${Ac_1_12}" file pair_1_12_contact_area_triaxial12.csv screen no
|
||||
|
||||
variable plane_disp equal ${ddisp}*elapsed
|
||||
variable plane_disp_neg equal -${ddisp}*elapsed
|
||||
|
||||
run ${compression_steps}
|
||||
23
examples/granular/spheres12.data
Normal file
23
examples/granular/spheres12.data
Normal file
@ -0,0 +1,23 @@
|
||||
#LAMMPS data file created by matlab.
|
||||
12 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-10.0000000000 10.0000000000 xlo xhi
|
||||
-10.0000000000 10.0000000000 ylo yhi
|
||||
-10.0000000000 10.0000000000 zlo zhi
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 0.8000000000 1000.0000000000 0.0717535226 -0.2092222842 0.3662146798
|
||||
2 1 1.2000000000 1000.0000000000 -0.8233763986 -0.7426114800 -0.8263932264
|
||||
3 1 0.8000000000 1000.0000000000 -1.0685863278 -0.4494609702 0.2196698078
|
||||
4 1 0.8000000000 1000.0000000000 0.5829432471 -1.0098803839 -0.7607543861
|
||||
5 1 0.8000000000 1000.0000000000 -0.8658471132 0.6951192569 0.0107556658
|
||||
6 1 1.2000000000 1000.0000000000 0.3966456126 0.7215053869 -0.7540113087
|
||||
7 1 1.2000000000 1000.0000000000 0.7316242921 0.8996483982 0.6751483031
|
||||
8 1 1.0000000000 1000.0000000000 0.6267527768 -0.8419367233 0.6964197101
|
||||
9 1 0.8000000000 1000.0000000000 -0.0409043189 -0.1452314035 -1.0102948313
|
||||
10 1 0.8000000000 1000.0000000000 -0.9495107709 0.6760151650 -0.9220534482
|
||||
11 1 1.0000000000 1000.0000000000 -0.7488486472 0.2188003421 0.7892021020
|
||||
12 1 1.2000000000 1000.0000000000 0.8968590780 -0.2350366437 -0.2006719701
|
||||
211
examples/granular/spheres200.data
Normal file
211
examples/granular/spheres200.data
Normal file
@ -0,0 +1,211 @@
|
||||
#LAMMPS data file created by matlab.
|
||||
200 atoms
|
||||
|
||||
1 atom types
|
||||
|
||||
-0.005000 0.005000 xlo xhi
|
||||
-0.005000 0.005000 ylo yhi
|
||||
-0.001000 0.020000 zlo zhi
|
||||
|
||||
Atoms
|
||||
|
||||
1 1 0.001206 1560.000000 -0.000938 0.000556 0.000883
|
||||
2 1 0.000953 1560.000000 -0.002626 -0.000145 0.002778
|
||||
3 1 0.001035 1560.000000 -0.000434 0.000172 0.008458
|
||||
4 1 0.001225 1560.000000 -0.003126 -0.000604 0.004986
|
||||
5 1 0.001119 1560.000000 0.000772 0.002972 0.002568
|
||||
6 1 0.001243 1560.000000 -0.000363 0.001184 0.004927
|
||||
7 1 0.001173 1560.000000 0.000218 0.000243 0.005475
|
||||
8 1 0.000937 1560.000000 0.000033 0.000029 0.003141
|
||||
9 1 0.001055 1560.000000 -0.001660 0.001975 0.008611
|
||||
10 1 0.000938 1560.000000 -0.001818 0.002352 0.002534
|
||||
11 1 0.000990 1560.000000 0.001592 0.000435 0.004416
|
||||
12 1 0.000927 1560.000000 -0.001659 -0.000004 0.005901
|
||||
13 1 0.001272 1560.000000 0.002972 0.000553 0.007291
|
||||
14 1 0.001226 1560.000000 0.002090 0.000983 0.001406
|
||||
15 1 0.000957 1560.000000 0.002241 -0.001608 0.001304
|
||||
16 1 0.001020 1560.000000 -0.001944 0.001290 0.002030
|
||||
17 1 0.001289 1560.000000 -0.002256 -0.001173 0.003474
|
||||
18 1 0.000998 1560.000000 0.000771 0.002127 0.000906
|
||||
19 1 0.000927 1560.000000 0.000186 0.000567 0.001207
|
||||
20 1 0.001095 1560.000000 -0.000937 -0.003179 0.008173
|
||||
21 1 0.001006 1560.000000 -0.001736 0.000751 0.004618
|
||||
22 1 0.001037 1560.000000 0.000784 0.001844 0.002380
|
||||
23 1 0.001297 1560.000000 0.000234 -0.001597 0.008560
|
||||
24 1 0.001017 1560.000000 0.002454 -0.000505 0.001171
|
||||
25 1 0.001110 1560.000000 -0.000803 -0.000415 0.003714
|
||||
26 1 0.001192 1560.000000 0.002283 0.000648 0.003048
|
||||
27 1 0.000992 1560.000000 -0.000065 -0.000545 0.007062
|
||||
28 1 0.001116 1560.000000 0.002174 -0.001463 0.005830
|
||||
29 1 0.001258 1560.000000 0.001602 0.001853 0.007246
|
||||
30 1 0.001055 1560.000000 -0.001535 -0.002770 0.007196
|
||||
31 1 0.000958 1560.000000 -0.000438 -0.000260 0.004709
|
||||
32 1 0.001188 1560.000000 0.000339 -0.000355 0.009171
|
||||
33 1 0.001166 1560.000000 0.002513 -0.001215 0.004434
|
||||
34 1 0.000907 1560.000000 0.001905 -0.000373 0.004921
|
||||
35 1 0.001245 1560.000000 -0.000091 -0.002620 0.004150
|
||||
36 1 0.001302 1560.000000 0.003292 0.000184 0.005377
|
||||
37 1 0.001305 1560.000000 0.002099 0.001261 0.008939
|
||||
38 1 0.000988 1560.000000 0.003274 0.000136 0.003667
|
||||
39 1 0.000892 1560.000000 0.001798 -0.002104 0.008610
|
||||
40 1 0.001247 1560.000000 -0.003058 -0.000575 0.000948
|
||||
41 1 0.000900 1560.000000 -0.000258 -0.000469 0.001478
|
||||
42 1 0.000945 1560.000000 -0.001434 -0.001711 0.004610
|
||||
43 1 0.000977 1560.000000 -0.001410 0.002808 0.004963
|
||||
44 1 0.000930 1560.000000 -0.002110 -0.001362 0.006749
|
||||
45 1 0.000931 1560.000000 0.001256 -0.000876 0.000844
|
||||
46 1 0.000901 1560.000000 0.000899 -0.001189 0.005316
|
||||
47 1 0.000940 1560.000000 -0.002189 -0.000047 0.007240
|
||||
48 1 0.001217 1560.000000 -0.000108 -0.001333 0.002257
|
||||
49 1 0.001088 1560.000000 0.001364 -0.000594 0.002789
|
||||
50 1 0.001143 1560.000000 -0.000311 -0.001425 0.006092
|
||||
51 1 0.001054 1560.000000 0.002262 0.002312 0.004315
|
||||
52 1 0.001016 1560.000000 -0.000724 0.000741 0.003295
|
||||
53 1 0.001051 1560.000000 0.000527 -0.001987 0.003307
|
||||
54 1 0.000905 1560.000000 0.000827 0.001457 0.005868
|
||||
55 1 0.001195 1560.000000 -0.001176 -0.000645 0.000798
|
||||
56 1 0.001253 1560.000000 0.002583 -0.001847 0.003310
|
||||
57 1 0.000982 1560.000000 0.001551 -0.002803 0.005076
|
||||
58 1 0.000945 1560.000000 -0.000481 0.000354 0.007220
|
||||
59 1 0.001040 1560.000000 -0.002736 0.001076 0.008769
|
||||
60 1 0.000917 1560.000000 0.000826 -0.001887 0.006449
|
||||
61 1 0.000914 1560.000000 -0.001171 -0.001592 0.007266
|
||||
62 1 0.000959 1560.000000 0.000834 -0.002671 0.007105
|
||||
63 1 0.000990 1560.000000 -0.000251 -0.001327 0.004339
|
||||
64 1 0.001220 1560.000000 0.001384 0.002896 0.005874
|
||||
65 1 0.000949 1560.000000 -0.001340 -0.000608 0.007496
|
||||
66 1 0.001306 1560.000000 0.002187 0.002068 0.002629
|
||||
67 1 0.001206 1560.000000 0.000148 0.001506 0.008517
|
||||
68 1 0.001123 1560.000000 0.001288 -0.000303 0.006613
|
||||
69 1 0.001151 1560.000000 -0.000876 0.001549 0.001740
|
||||
70 1 0.001315 1560.000000 -0.001902 -0.002590 0.001344
|
||||
71 1 0.000927 1560.000000 0.002285 -0.000866 0.006900
|
||||
72 1 0.001279 1560.000000 -0.000165 0.002689 0.007449
|
||||
73 1 0.000910 1560.000000 0.001009 0.001054 0.005049
|
||||
74 1 0.001148 1560.000000 -0.002229 -0.001285 0.008736
|
||||
75 1 0.001067 1560.000000 -0.000261 -0.002945 0.002157
|
||||
76 1 0.000993 1560.000000 -0.001641 0.002272 0.007601
|
||||
77 1 0.001228 1560.000000 0.001939 -0.000214 0.008903
|
||||
78 1 0.001076 1560.000000 0.000767 0.001172 0.003556
|
||||
79 1 0.001105 1560.000000 -0.000561 0.002493 0.004214
|
||||
80 1 0.001195 1560.000000 0.002694 -0.000817 0.007949
|
||||
81 1 0.001239 1560.000000 -0.000968 -0.003145 0.006096
|
||||
82 1 0.001083 1560.000000 -0.000808 0.001813 0.006396
|
||||
83 1 0.000923 1560.000000 0.000632 -0.001437 0.001310
|
||||
84 1 0.000981 1560.000000 -0.001842 0.002774 0.006508
|
||||
85 1 0.000998 1560.000000 -0.002775 0.001616 0.001453
|
||||
86 1 0.000979 1560.000000 -0.002520 0.001715 0.007741
|
||||
87 1 0.001002 1560.000000 -0.001465 -0.001931 0.006048
|
||||
88 1 0.000958 1560.000000 0.003264 0.000707 0.001189
|
||||
89 1 0.001052 1560.000000 -0.001314 -0.000701 0.002721
|
||||
90 1 0.001096 1560.000000 0.001154 0.002129 0.004403
|
||||
91 1 0.001104 1560.000000 0.002118 0.001977 0.000794
|
||||
92 1 0.001263 1560.000000 -0.001499 -0.002764 0.003441
|
||||
93 1 0.001086 1560.000000 -0.001096 0.002514 0.001154
|
||||
94 1 0.000895 1560.000000 0.001130 0.000029 0.001045
|
||||
95 1 0.000964 1560.000000 0.000905 -0.003200 0.000542
|
||||
96 1 0.000898 1560.000000 -0.000868 0.003148 0.008306
|
||||
97 1 0.000907 1560.000000 -0.001406 0.001144 0.007862
|
||||
98 1 0.001176 1560.000000 0.001246 -0.001074 0.004327
|
||||
99 1 0.001148 1560.000000 0.001512 -0.002739 0.003346
|
||||
100 1 0.000922 1560.000000 0.001470 -0.000036 0.007695
|
||||
101 1 0.001031 1560.000000 -0.002751 0.000928 0.004124
|
||||
102 1 0.001030 1560.000000 -0.000177 -0.002370 0.005374
|
||||
103 1 0.000915 1560.000000 0.000824 0.000521 0.007070
|
||||
104 1 0.001085 1560.000000 -0.002281 -0.000023 0.009123
|
||||
105 1 0.001004 1560.000000 -0.000167 0.002610 0.008905
|
||||
106 1 0.001060 1560.000000 -0.000389 -0.002220 0.007688
|
||||
107 1 0.000920 1560.000000 -0.000483 0.003231 0.006505
|
||||
108 1 0.001122 1560.000000 0.001781 -0.001547 0.002237
|
||||
109 1 0.001172 1560.000000 -0.002650 0.000830 0.005429
|
||||
110 1 0.001137 1560.000000 -0.000030 -0.003246 0.001024
|
||||
111 1 0.001315 1560.000000 0.001470 -0.001735 0.007580
|
||||
112 1 0.001245 1560.000000 0.000481 -0.003067 0.006025
|
||||
113 1 0.000904 1560.000000 0.000632 -0.000184 0.002010
|
||||
114 1 0.000883 1560.000000 -0.001828 0.002191 0.003819
|
||||
115 1 0.000974 1560.000000 0.002167 0.001616 0.006226
|
||||
116 1 0.001150 1560.000000 0.000871 -0.002731 0.002136
|
||||
117 1 0.001312 1560.000000 -0.000326 -0.001971 0.001000
|
||||
118 1 0.000914 1560.000000 0.001020 0.000810 0.002086
|
||||
119 1 0.001136 1560.000000 -0.000101 -0.003277 0.007246
|
||||
120 1 0.000991 1560.000000 -0.001944 0.000576 0.003215
|
||||
121 1 0.001216 1560.000000 -0.000913 -0.001165 0.008857
|
||||
122 1 0.001045 1560.000000 -0.003110 0.001062 0.002973
|
||||
123 1 0.000918 1560.000000 0.000348 0.000365 0.004046
|
||||
124 1 0.001279 1560.000000 -0.000884 0.003087 0.002268
|
||||
125 1 0.001065 1560.000000 -0.002238 0.001309 0.006452
|
||||
126 1 0.001012 1560.000000 -0.002059 -0.001354 0.001935
|
||||
127 1 0.001142 1560.000000 -0.003011 0.000567 0.001739
|
||||
128 1 0.000921 1560.000000 0.001764 0.002804 0.008177
|
||||
129 1 0.001151 1560.000000 -0.003105 -0.000384 0.006602
|
||||
130 1 0.000967 1560.000000 0.000932 0.000588 0.008823
|
||||
131 1 0.000908 1560.000000 -0.001873 -0.001947 0.007825
|
||||
132 1 0.000923 1560.000000 -0.002993 0.000883 0.007425
|
||||
133 1 0.001171 1560.000000 0.003310 -0.000405 0.006558
|
||||
134 1 0.000977 1560.000000 -0.000098 -0.000180 0.000492
|
||||
135 1 0.000938 1560.000000 -0.000706 -0.000129 0.006085
|
||||
136 1 0.001008 1560.000000 -0.000256 0.002333 0.000550
|
||||
137 1 0.001073 1560.000000 0.000534 -0.000055 0.008080
|
||||
138 1 0.000890 1560.000000 0.000351 0.001695 0.007195
|
||||
139 1 0.000973 1560.000000 0.002593 0.001907 0.005394
|
||||
140 1 0.001176 1560.000000 -0.001862 -0.000534 0.004494
|
||||
141 1 0.001306 1560.000000 -0.000951 0.001053 0.009299
|
||||
142 1 0.001103 1560.000000 -0.001937 -0.002711 0.008485
|
||||
143 1 0.001262 1560.000000 -0.002947 -0.001470 0.007682
|
||||
144 1 0.000914 1560.000000 0.002047 0.000811 0.005504
|
||||
145 1 0.000954 1560.000000 0.001935 -0.002349 0.006632
|
||||
146 1 0.001003 1560.000000 0.000766 -0.002635 0.008483
|
||||
147 1 0.001137 1560.000000 0.000102 0.003195 0.004922
|
||||
148 1 0.001006 1560.000000 -0.001982 0.001014 0.000685
|
||||
149 1 0.001255 1560.000000 -0.000718 0.001939 0.003056
|
||||
150 1 0.001057 1560.000000 -0.001189 -0.001717 0.003045
|
||||
151 1 0.001228 1560.000000 0.001581 0.002926 0.003510
|
||||
152 1 0.001052 1560.000000 -0.002172 0.001949 0.004831
|
||||
153 1 0.000979 1560.000000 -0.001817 0.000291 0.002048
|
||||
154 1 0.001286 1560.000000 -0.002647 -0.001839 0.004620
|
||||
155 1 0.001085 1560.000000 -0.000081 0.000850 0.002139
|
||||
156 1 0.000990 1560.000000 -0.000081 0.002105 0.005587
|
||||
157 1 0.001043 1560.000000 0.001636 -0.000112 0.001860
|
||||
158 1 0.001309 1560.000000 0.003216 -0.000851 0.002791
|
||||
159 1 0.000913 1560.000000 0.000608 0.003148 0.006565
|
||||
160 1 0.000919 1560.000000 0.000536 -0.003106 0.003249
|
||||
161 1 0.000943 1560.000000 0.003145 -0.000528 0.008915
|
||||
162 1 0.000993 1560.000000 -0.002811 -0.000099 0.008110
|
||||
163 1 0.001125 1560.000000 0.001415 -0.002271 0.000643
|
||||
164 1 0.000919 1560.000000 -0.001406 0.000223 0.006781
|
||||
165 1 0.001040 1560.000000 0.000690 0.003193 0.008329
|
||||
166 1 0.001055 1560.000000 0.001075 0.002584 0.009093
|
||||
167 1 0.001176 1560.000000 0.000851 0.003176 0.000591
|
||||
168 1 0.001003 1560.000000 -0.001462 0.001511 0.005544
|
||||
169 1 0.001126 1560.000000 -0.000077 0.003324 0.001347
|
||||
170 1 0.001068 1560.000000 0.003110 0.000810 0.008495
|
||||
171 1 0.001011 1560.000000 -0.001661 0.000117 0.008201
|
||||
172 1 0.001066 1560.000000 -0.000359 -0.003279 0.009094
|
||||
173 1 0.001303 1560.000000 0.003066 0.001188 0.004082
|
||||
174 1 0.000983 1560.000000 0.000354 0.002261 0.003558
|
||||
175 1 0.001137 1560.000000 0.002860 -0.001571 0.009180
|
||||
176 1 0.001070 1560.000000 0.001246 -0.001279 0.009104
|
||||
177 1 0.000886 1560.000000 0.002271 -0.000316 0.003675
|
||||
178 1 0.000983 1560.000000 -0.001987 -0.002490 0.005377
|
||||
179 1 0.000939 1560.000000 0.000601 -0.000861 0.003477
|
||||
180 1 0.001177 1560.000000 0.001522 0.002902 0.001690
|
||||
181 1 0.001036 1560.000000 -0.001200 -0.002874 0.004750
|
||||
182 1 0.000898 1560.000000 -0.001705 -0.001140 0.005503
|
||||
183 1 0.001315 1560.000000 0.002732 0.001766 0.007885
|
||||
184 1 0.001318 1560.000000 -0.002909 -0.001610 0.005936
|
||||
185 1 0.001218 1560.000000 0.003213 0.000884 0.002316
|
||||
186 1 0.001234 1560.000000 -0.002394 -0.002298 0.002575
|
||||
187 1 0.001160 1560.000000 -0.003313 -0.000065 0.003625
|
||||
188 1 0.001022 1560.000000 -0.003096 -0.001048 0.002151
|
||||
189 1 0.000966 1560.000000 0.001891 -0.002093 0.004404
|
||||
190 1 0.001048 1560.000000 -0.002367 0.002338 0.000697
|
||||
191 1 0.000995 1560.000000 -0.001204 -0.001912 0.002030
|
||||
192 1 0.001136 1560.000000 -0.001152 -0.002402 0.009223
|
||||
193 1 0.001083 1560.000000 -0.002588 -0.001768 0.000753
|
||||
194 1 0.000946 1560.000000 -0.001338 -0.000741 0.006527
|
||||
195 1 0.000943 1560.000000 -0.000073 0.003254 0.003663
|
||||
196 1 0.001059 1560.000000 0.000087 0.000958 0.006388
|
||||
197 1 0.001131 1560.000000 0.001030 0.001019 0.000752
|
||||
198 1 0.001257 1560.000000 -0.001365 0.002946 0.009266
|
||||
199 1 0.000891 1560.000000 -0.000445 -0.000273 0.002382
|
||||
200 1 0.001055 1560.000000 0.001781 0.000748 0.006583
|
||||
2
src/.gitignore
vendored
2
src/.gitignore
vendored
@ -851,6 +851,8 @@
|
||||
/fix_ffl.h
|
||||
/fix_filter_corotate.cpp
|
||||
/fix_filter_corotate.h
|
||||
/fix_granular_mdr.cpp
|
||||
/fix_granular_mdr.h
|
||||
/fix_viscosity.cpp
|
||||
/fix_viscosity.h
|
||||
/fix_ehex.cpp
|
||||
|
||||
702
src/GRANULAR/fix_granular_mdr.cpp
Normal file
702
src/GRANULAR/fix_granular_mdr.cpp
Normal file
@ -0,0 +1,702 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors:
|
||||
William Zunker (MIT), Sachith Dunatunga (MIT),
|
||||
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
|
||||
----------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_granular_mdr.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "fix_wall_gran_region.h"
|
||||
#include "fix_neigh_history.h"
|
||||
#include "force.h"
|
||||
#include "granular_model.h"
|
||||
#include "gran_sub_mod_normal.h"
|
||||
#include "input.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neigh_list.h"
|
||||
#include "pair.h"
|
||||
#include "pair_granular.h"
|
||||
#include "region.h"
|
||||
#include "update.h"
|
||||
#include "variable.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
using namespace Granular_MDR_NS;
|
||||
using namespace FixConst;
|
||||
using MathConst::MY_PI;
|
||||
|
||||
static constexpr double EPSILON = 1e-16;
|
||||
static constexpr double OVERLAP_LIMIT = 0.75;
|
||||
|
||||
enum {COMM_1, COMM_2};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixGranularMDR::FixGranularMDR(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
comm_forward = 5;
|
||||
create_attribute = 1;
|
||||
|
||||
id_fix = nullptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixGranularMDR::~FixGranularMDR()
|
||||
{
|
||||
if (id_fix && modify->nfix)
|
||||
modify->delete_fix(id_fix);
|
||||
delete[] id_fix;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixGranularMDR::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::post_constructor()
|
||||
{
|
||||
int tmp1, tmp2;
|
||||
id_fix = utils::strdup("MDR_PARTICLE_HISTORY_VARIABLES");
|
||||
modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz ghost yes", id_fix));
|
||||
|
||||
index_Ro = atom->find_custom("Ro", tmp1, tmp2);
|
||||
index_Vcaps = atom->find_custom("Vcaps", tmp1, tmp2);
|
||||
index_Vgeo = atom->find_custom("Vgeo", tmp1, tmp2);
|
||||
index_Velas = atom->find_custom("Velas", tmp1, tmp2);
|
||||
index_eps_bar = atom->find_custom("eps_bar", tmp1, tmp2);
|
||||
index_dRnumerator = atom->find_custom("dRnumerator", tmp1, tmp2);
|
||||
index_dRdenominator = atom->find_custom("dRdenominator", tmp1, tmp2);
|
||||
index_Acon0 = atom->find_custom("Acon0", tmp1, tmp2);
|
||||
index_Acon1 = atom->find_custom("Acon1", tmp1, tmp2);
|
||||
index_Atot = atom->find_custom("Atot", tmp1, tmp2);
|
||||
index_Atot_sum = atom->find_custom("Atot_sum", tmp1, tmp2);
|
||||
index_ddelta_bar = atom->find_custom("ddelta_bar", tmp1, tmp2);
|
||||
index_psi = atom->find_custom("psi", tmp1, tmp2);
|
||||
index_history_setup_flag = atom->find_custom("history_setup_flag", tmp1, tmp2);
|
||||
index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2);
|
||||
index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2);
|
||||
index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::setup_pre_force(int /*vflag*/)
|
||||
{
|
||||
pair = dynamic_cast<PairGranular *>(force->pair_match("granular", 1));
|
||||
if (pair == nullptr)
|
||||
error->all(FLERR, "Must use pair granular with MDR model");
|
||||
|
||||
if (force->newton)
|
||||
error->all(FLERR, "MDR contact model requires Newton off");
|
||||
|
||||
// Confirm all MDR models are consistent
|
||||
|
||||
class GranularModel *pair_model, *fix_model;
|
||||
class GranularModel **models_list = pair->models_list;
|
||||
class GranSubModNormalMDR *norm_model = nullptr;
|
||||
for (int i = 0; i < pair->nmodels; i++) {
|
||||
pair_model = models_list[i];
|
||||
if (pair_model->normal_model->name == "mdr") {
|
||||
if (norm_model != nullptr)
|
||||
error->all(FLERR, "Cannot currently define multiple MDR normal models in the pairstyle");
|
||||
norm_model = dynamic_cast<GranSubModNormalMDR *>(pair_model->normal_model);
|
||||
} else {
|
||||
error->all(FLERR, "Cannot combine MDR normal model with a different normal model in the pairstyle");
|
||||
}
|
||||
}
|
||||
|
||||
if (norm_model == nullptr)
|
||||
error->all(FLERR, "Must specify MDR normal model with pair granular");
|
||||
psi_b_coeff = norm_model->psi_b;
|
||||
|
||||
fix_wall_list = modify->get_fix_by_style("wall/gran/region");
|
||||
class GranSubModNormalMDR *norm_model2;
|
||||
class FixWallGranRegion *fix;
|
||||
for (int i = 0; i < fix_wall_list.size(); i++) {
|
||||
if (!utils::strmatch(fix_wall_list[i]->style, "wall/gran/region"))
|
||||
error->all(FLERR, "MDR model currently only supports fix wall/gran/region, not fix wall/gran");
|
||||
|
||||
fix = dynamic_cast<FixWallGranRegion*>(fix_wall_list[i]);
|
||||
if (fix->model->normal_model->name != "mdr")
|
||||
error->all(FLERR, "Fix wall/gran/region must use an MDR normal model when using an MDR pair model");
|
||||
|
||||
norm_model2 = dynamic_cast<GranSubModNormalMDR *>(fix->model->normal_model);
|
||||
|
||||
if (norm_model->E != norm_model2->E)
|
||||
error->all(FLERR, "Young's modulus in pair style, {}, does not agree with value {} in fix gran/wall/region",
|
||||
norm_model->E, norm_model2->E);
|
||||
if (norm_model->nu != norm_model2->nu)
|
||||
error->all(FLERR, "Poisson's ratio in pair style, {}, does not agree with value {} in fix gran/wall/region",
|
||||
norm_model->nu, norm_model2->nu);
|
||||
if (norm_model->Y != norm_model2->Y)
|
||||
error->all(FLERR, "Yield stress in pair style, {}, does not agree with value {} in fix gran/wall/region",
|
||||
norm_model->Y, norm_model2->Y);
|
||||
if (norm_model->psi_b != norm_model2->psi_b)
|
||||
error->all(FLERR, "Bulk response trigger in pair style, {}, does not agree with value {} in fix gran/wall/region",
|
||||
norm_model->psi_b, norm_model2->psi_b);
|
||||
if (norm_model->CoR != norm_model2->CoR)
|
||||
error->all(FLERR, "Coefficient of restitution in pair style, {}, does not agree with value {} in fix gran/wall/region",
|
||||
norm_model->CoR, norm_model2->CoR);
|
||||
}
|
||||
|
||||
fix_history = dynamic_cast<FixNeighHistory *>(modify->get_fix_by_id("NEIGH_HISTORY_GRANULAR"));
|
||||
|
||||
pre_force(0);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::pre_force(int)
|
||||
{
|
||||
double *radius = atom->radius;
|
||||
double *Ro = atom->dvector[index_Ro];
|
||||
double *Vgeo = atom->dvector[index_Vgeo];
|
||||
double *Velas = atom->dvector[index_Velas];
|
||||
double *Vcaps = atom->dvector[index_Vcaps];
|
||||
double *eps_bar = atom->dvector[index_eps_bar];
|
||||
double *dRnumerator = atom->dvector[index_dRnumerator];
|
||||
double *dRdenominator = atom->dvector[index_dRdenominator];
|
||||
double *Acon0 = atom->dvector[index_Acon0];
|
||||
double *Acon1 = atom->dvector[index_Acon1];
|
||||
double *Atot = atom->dvector[index_Atot];
|
||||
double *Atot_sum = atom->dvector[index_Atot_sum];
|
||||
double *psi = atom->dvector[index_psi];
|
||||
double *ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
double *sigmaxx = atom->dvector[index_sigmaxx];
|
||||
double *sigmayy = atom->dvector[index_sigmayy];
|
||||
double *sigmazz = atom->dvector[index_sigmazz];
|
||||
double *history_setup_flag = atom->dvector[index_history_setup_flag];
|
||||
|
||||
int new_atom;
|
||||
int nlocal = atom->nlocal;
|
||||
int ntotal = nlocal + atom->nghost;
|
||||
for (int i = 0; i < ntotal; i++) {
|
||||
// initialize new atoms
|
||||
new_atom = 0;
|
||||
if (history_setup_flag[i] < EPSILON) {
|
||||
Ro[i] = radius[i];
|
||||
Acon0[i] = 0.0;
|
||||
Acon1[i] = 0.0;
|
||||
Vcaps[i] = 0.0;
|
||||
eps_bar[i] = 0.0;
|
||||
dRnumerator[i] = 0.0;
|
||||
dRdenominator[i] = 0.0;
|
||||
Atot_sum[i] = 0.0;
|
||||
ddelta_bar[i] = 0.0;
|
||||
sigmaxx[i] = 0.0;
|
||||
sigmayy[i] = 0.0;
|
||||
sigmazz[i] = 0.0;
|
||||
history_setup_flag[i] = 1.0;
|
||||
new_atom = 1;
|
||||
}
|
||||
|
||||
// update apparent radius
|
||||
|
||||
// will forward to ghosts
|
||||
if (i >= nlocal) continue;
|
||||
|
||||
// only update outside of setup (unless a new atom)
|
||||
if (update->setupflag && (!new_atom)) continue;
|
||||
|
||||
const double R = radius[i];
|
||||
const double Vo = 4.0 / 3.0 * MY_PI * pow(Ro[i], 3.0);
|
||||
const double Vgeoi = 4.0 / 3.0 * MY_PI * pow(R, 3.0) - Vcaps[i];
|
||||
|
||||
Vgeo[i] = MIN(Vgeoi, Vo);
|
||||
Velas[i] = Vo * (1.0 + eps_bar[i]);
|
||||
Atot[i] = 4.0 * MY_PI * pow(R, 2.0) + Atot_sum[i];
|
||||
psi[i] = (Atot[i] - Acon1[i]) / Atot[i];
|
||||
|
||||
if (psi_b_coeff < psi[i]) {
|
||||
const double dR = MAX(dRnumerator[i] / (dRdenominator[i] - 4.0 * MY_PI * pow(R, 2.0)), 0.0);
|
||||
if ((radius[i] + dR) < (1.5 * Ro[i]))
|
||||
radius[i] += dR;
|
||||
}
|
||||
Acon0[i] = Acon1[i];
|
||||
}
|
||||
|
||||
comm_stage = COMM_1;
|
||||
comm->forward_comm(this, 5);
|
||||
|
||||
// rezero temporary variables for all atoms, no need to communicate
|
||||
for (int i = 0; i < ntotal; i++) {
|
||||
ddelta_bar[i] = 0.0;
|
||||
if (!update->setupflag) {
|
||||
Vcaps[i] = 0.0;
|
||||
eps_bar[i] = 0.0;
|
||||
dRnumerator[i] = 0.0;
|
||||
dRdenominator[i] = 0.0;
|
||||
Acon1[i] = 0.0;
|
||||
Atot_sum[i] = 0.0;
|
||||
sigmaxx[i] = 0.0;
|
||||
sigmayy[i] = 0.0;
|
||||
sigmazz[i] = 0.0;
|
||||
}
|
||||
}
|
||||
if (!update->setupflag) {
|
||||
calculate_contact_penalty();
|
||||
mean_surf_disp();
|
||||
update_fix_gran_wall();
|
||||
}
|
||||
|
||||
comm_stage = COMM_2;
|
||||
comm->forward_comm(this, 1);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixGranularMDR::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,int * /*pbc*/)
|
||||
{
|
||||
double **dvector = atom->dvector;
|
||||
int m = 0;
|
||||
if (comm_stage == COMM_1) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
buf[m++] = dvector[index_Vgeo][j]; // 2
|
||||
buf[m++] = dvector[index_Velas][j]; // 3
|
||||
buf[m++] = dvector[index_Acon0][j]; // 8
|
||||
buf[m++] = dvector[index_Atot][j]; // 10
|
||||
buf[m++] = dvector[index_psi][j]; // 13
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
buf[m++] = dvector[index_ddelta_bar][j];
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
double **dvector = atom->dvector;
|
||||
int m = 0;
|
||||
int last = first + n;
|
||||
|
||||
if (comm_stage == COMM_1) {
|
||||
for (int i = first; i < last; i++) {
|
||||
dvector[index_Vgeo][i] = buf[m++]; // 2
|
||||
dvector[index_Velas][i] = buf[m++]; // 3
|
||||
dvector[index_Acon0][i] = buf[m++]; // 8
|
||||
dvector[index_Atot][i] = buf[m++]; // 10
|
||||
dvector[index_psi][i] = buf[m++]; // 13
|
||||
}
|
||||
} else {
|
||||
for (int i = first; i < last; i++) {
|
||||
dvector[index_ddelta_bar][i] = buf[m++];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
initialize setup flag to zero, called when atom is created
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::set_arrays(int i)
|
||||
{
|
||||
atom->dvector[index_history_setup_flag][i] = 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Screen for non-physical contacts occuring through obstructing particles.
|
||||
Assign non-zero penalties to these contacts to adjust force evaluation.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::calculate_contact_penalty()
|
||||
{
|
||||
NeighList * list = pair->list;
|
||||
const int size_history = pair->get_size_history();
|
||||
|
||||
int i, j, k, lv1, ii, jj, inum, jnum;
|
||||
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
int *touch, **firsttouch;
|
||||
double *history, *history_ij, *history_ik, *history_jk, *history_kj;
|
||||
double *allhistory, *allhistory_j, *allhistory_k, **firsthistory;
|
||||
|
||||
bool touchflag = false;
|
||||
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
firsttouch = fix_history->firstflag;
|
||||
firsthistory = fix_history->firstvalue;
|
||||
|
||||
// zero existing penalties
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
allhistory = firsthistory[i];
|
||||
jnum = numneigh[i];
|
||||
for (jj = 0; jj < jnum; jj++)
|
||||
(&allhistory[size_history * jj])[PENALTY] = 0.0;
|
||||
}
|
||||
|
||||
// contact penalty calculation
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
const double xtmp = x[i][0];
|
||||
const double ytmp = x[i][1];
|
||||
const double ztmp = x[i][2];
|
||||
allhistory = firsthistory[i];
|
||||
double radi = radius[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
double radj = radius[j];
|
||||
const double delx_ij = x[j][0] - xtmp;
|
||||
const double dely_ij = x[j][1] - ytmp;
|
||||
const double delz_ij = x[j][2] - ztmp;
|
||||
const double rsq_ij = delx_ij * delx_ij + dely_ij * dely_ij + delz_ij * delz_ij;
|
||||
const double r_ij = sqrt(rsq_ij);
|
||||
const double rinv_ij = 1.0 / r_ij;
|
||||
const double radsum_ij = radi + radj;
|
||||
const double deltan_ij = radsum_ij - r_ij;
|
||||
if (deltan_ij < 0.0) continue;
|
||||
for (int kk = jj + 1; kk < jnum; kk++) {
|
||||
k = jlist[kk];
|
||||
k &= NEIGHMASK;
|
||||
|
||||
const double delx_ik = x[k][0] - xtmp;
|
||||
const double dely_ik = x[k][1] - ytmp;
|
||||
const double delz_ik = x[k][2] - ztmp;
|
||||
const double rsq_ik = delx_ik * delx_ik + dely_ik * dely_ik + delz_ik *delz_ik;
|
||||
const double r_ik = sqrt(rsq_ik);
|
||||
const double rinv_ik = 1.0 / r_ik;
|
||||
const double radk = radius[k];
|
||||
const double radsum_ik = radi + radk;
|
||||
const double deltan_ik = radsum_ik - r_ik;
|
||||
|
||||
if (deltan_ik < 0.0) continue;
|
||||
|
||||
const double delx_jk = x[k][0] - x[j][0];
|
||||
const double dely_jk = x[k][1] - x[j][1];
|
||||
const double delz_jk = x[k][2] - x[j][2];
|
||||
const double rsq_jk = delx_jk * delx_jk + dely_jk * dely_jk + delz_jk *delz_jk;
|
||||
const double r_jk = sqrt(rsq_jk);
|
||||
const double rinv_jk = 1.0 / r_jk;
|
||||
const double radsum_jk = radj + radk;
|
||||
const double deltan_jk = radsum_jk - r_jk;
|
||||
|
||||
if (deltan_jk < 0.0) continue;
|
||||
|
||||
// pull ij history
|
||||
history_ij = &allhistory[size_history * jj];
|
||||
double * pij = &history_ij[PENALTY]; // penalty for contact i and j
|
||||
|
||||
// pull ik history
|
||||
history_ik = &allhistory[size_history * kk];
|
||||
double * pik = &history_ik[PENALTY]; // penalty for contact i and k
|
||||
|
||||
// Find pair of atoms with the smallest overlap, atoms a & b, 3rd atom c is central
|
||||
// if a & b are both local:
|
||||
// calculate ab penalty and add to the pab[0] history entry
|
||||
// if a is local & b is ghost or vice versa:
|
||||
// each processor has a-b in nlist and independently calculates + adds penalty
|
||||
// if a & b are both ghosts:
|
||||
// skip calculation since it's performed on other proc
|
||||
// This process requires newton off, or nlist may not include ab, ac, & bc
|
||||
|
||||
const double r_max = MAX(r_ij, MAX(r_ik, r_jk));
|
||||
if (r_ij == r_max) { // the central particle is k
|
||||
const double enx_ki = -delx_ik * rinv_ik;
|
||||
const double eny_ki = -dely_ik * rinv_ik;
|
||||
const double enz_ki = -delz_ik * rinv_ik;
|
||||
const double enx_kj = -delx_jk * rinv_jk;
|
||||
const double eny_kj = -dely_jk * rinv_jk;
|
||||
const double enz_kj = -delz_jk * rinv_jk;
|
||||
const double alpha = std::acos(enx_ki * enx_kj + eny_ki * eny_kj + enz_ki * enz_kj);
|
||||
pij[0] += 1.0 / (1.0 + std::exp(-50.0 * (alpha / MY_PI - 0.5)));
|
||||
} else if (r_ik == r_max) { // the central particle is j
|
||||
const double enx_ji = -delx_ij * rinv_ij;
|
||||
const double eny_ji = -dely_ij * rinv_ij;
|
||||
const double enz_ji = -delz_ij * rinv_ij;
|
||||
const double enx_jk = delx_jk * rinv_jk;
|
||||
const double eny_jk = dely_jk * rinv_jk;
|
||||
const double enz_jk = delz_jk * rinv_jk;
|
||||
const double alpha = std::acos(enx_ji * enx_jk + eny_ji * eny_jk + enz_ji * enz_jk);
|
||||
pik[0] += 1.0 / (1.0 + std::exp(-50.0 * (alpha / MY_PI - 0.5)));
|
||||
} else { // the central particle is i
|
||||
if (j < atom->nlocal || k < atom->nlocal) {
|
||||
const double enx_ij = delx_ij * rinv_ij;
|
||||
const double eny_ij = dely_ij * rinv_ij;
|
||||
const double enz_ij = delz_ij * rinv_ij;
|
||||
const double enx_ik = delx_ik * rinv_ik;
|
||||
const double eny_ik = dely_ik * rinv_ik;
|
||||
const double enz_ik = delz_ik * rinv_ik;
|
||||
const double alpha = std::acos(enx_ij * enx_ik + eny_ij * eny_ik + enz_ij * enz_ik);
|
||||
|
||||
// don't know who owns the contact, k may be in j's nlist or vice versa
|
||||
// need to search both to find owner
|
||||
double * pjk = nullptr;
|
||||
if (j < atom->nlocal) {
|
||||
int * const jklist = firstneigh[j];
|
||||
const int jknum = numneigh[j];
|
||||
for (int jk = 0; jk < jknum; jk++) {
|
||||
const int kneigh = jklist[jk] & NEIGHMASK;
|
||||
if (k == kneigh) {
|
||||
allhistory_j = firsthistory[j];
|
||||
history_jk = &allhistory_j[size_history * jk];
|
||||
pjk = &history_jk[PENALTY]; // penalty for contact j and k
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// check if j is in the neighbor list of k
|
||||
if (pjk == nullptr && k < atom->nlocal) {
|
||||
int * const kjlist = firstneigh[k];
|
||||
const int kjnum = numneigh[k];
|
||||
for (int kj = 0; kj < kjnum; kj++) {
|
||||
const int jneigh = kjlist[kj] & NEIGHMASK;
|
||||
if (j == jneigh) {
|
||||
allhistory_k = firsthistory[k];
|
||||
history_kj = &allhistory_k[size_history * kj];
|
||||
pjk = &history_kj[PENALTY]; // penalty for contact j and k
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (pjk == nullptr)
|
||||
error->one(FLERR, "Contact between a pair of particles was detected by the MDR model, however it is not reflected in the neighbor lists. To solve this issue either build the neighbor lists more frequently or increase their size (e.g. increase the skin distance).");
|
||||
|
||||
pjk[0] += 1.0 / (1.0 + std::exp(-50.0 * (alpha / MY_PI - 0.5)));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Calculate mean surface displacement increment for each particle
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::mean_surf_disp()
|
||||
{
|
||||
NeighList * list = pair->list;
|
||||
|
||||
const int size_history = pair->get_size_history();
|
||||
int i, j, k, ii, jj, inum, jnum, itype, jtype;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
int *touch, **firsttouch;
|
||||
double *history, *allhistory, **firsthistory;
|
||||
|
||||
bool touchflag = false;
|
||||
class GranularModel* model;
|
||||
class GranularModel** models_list = pair->models_list;
|
||||
int ** types_indices = pair->types_indices;
|
||||
|
||||
double **x = atom->x;
|
||||
int *type = atom->type;
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *Acon0 = atom->dvector[index_Acon0];
|
||||
double *ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
firsttouch = fix_history->firstflag;
|
||||
firsthistory = fix_history->firstvalue;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
itype = type[i];
|
||||
touch = firsttouch[i];
|
||||
allhistory = firsthistory[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
jtype = type[j];
|
||||
model = models_list[types_indices[itype][jtype]];
|
||||
|
||||
// Reset model and copy initial geometric data
|
||||
model->xi = x[i];
|
||||
model->xj = x[j];
|
||||
model->radi = radius[i];
|
||||
model->radj = radius[j];
|
||||
model->i = i;
|
||||
model->j = j;
|
||||
model->touch = touch[jj];
|
||||
touchflag = model->check_contact();
|
||||
|
||||
// is it necessary to clear the history here???
|
||||
if (!touchflag) {
|
||||
touch[jj] = 0;
|
||||
history = &allhistory[size_history * jj];
|
||||
for (k = 0; k < size_history; k++) history[k] = 0.0;
|
||||
continue;
|
||||
}
|
||||
|
||||
touch[jj] = 1;
|
||||
|
||||
history = &allhistory[size_history * jj];
|
||||
model->history = history;
|
||||
|
||||
const double delta = model->radsum - sqrt(model->rsq);
|
||||
|
||||
double deltamax = history[DELTA_MAX];
|
||||
double deltap0 = history[DELTAP_0];
|
||||
double deltap1 = history[DELTAP_1];
|
||||
|
||||
if (delta > deltamax) deltamax = delta;
|
||||
|
||||
double delta0old = history[DELTA_0];
|
||||
double delta1old = history[DELTA_1];
|
||||
|
||||
int i0;
|
||||
int i1;
|
||||
if (atom->tag[i] > atom->tag[j]) {
|
||||
i0 = i;
|
||||
i1 = j;
|
||||
} else {
|
||||
i0 = j;
|
||||
i1 = i;
|
||||
}
|
||||
|
||||
double R0 = radius[i0];
|
||||
double R1 = radius[i1];
|
||||
|
||||
double delta_geo0;
|
||||
double delta_geo1;
|
||||
double deltaOpt1 = deltamax * (deltamax - 2.0 * R1) / (2.0 * (deltamax - R0 - R1));
|
||||
double deltaOpt2 = deltamax * (deltamax - 2.0 * R0) / (2.0 * (deltamax - R0 - R1));
|
||||
(R0 < R1) ? delta_geo0 = MAX(deltaOpt1, deltaOpt2) : delta_geo0 = MIN(deltaOpt1, deltaOpt2);
|
||||
(R0 < R1) ? delta_geo1 = MIN(deltaOpt1, deltaOpt2) : delta_geo1 = MAX(deltaOpt1, deltaOpt2);
|
||||
|
||||
if (delta_geo0 / R0 > OVERLAP_LIMIT) {
|
||||
delta_geo0 = R0 * OVERLAP_LIMIT;
|
||||
delta_geo1 = deltamax - delta_geo0;
|
||||
} else if (delta_geo1 / R1 > OVERLAP_LIMIT) {
|
||||
delta_geo1 = R1 * OVERLAP_LIMIT;
|
||||
delta_geo0 = deltamax - delta_geo1;
|
||||
}
|
||||
|
||||
double deltap = deltap0 + deltap1;
|
||||
|
||||
double delta0 = delta_geo0 + (deltap0 - delta_geo0) / (deltap - deltamax) * (delta - deltamax);
|
||||
double delta1 = delta_geo1 + (deltap1 - delta_geo1) / (deltap - deltamax) * (delta - deltamax);
|
||||
|
||||
double ddel0 = delta0 - delta0old;
|
||||
double ddel1 = delta1 - delta1old;
|
||||
|
||||
if (Acon0[i0] != 0.0) {
|
||||
const double Ac_offset0 = history[AC_0];
|
||||
ddelta_bar[i0] += Ac_offset0 / Acon0[i0] * ddel0;
|
||||
}
|
||||
|
||||
if (Acon0[i1] != 0.0) {
|
||||
const double Ac_offset1 = history[AC_1];
|
||||
ddelta_bar[i1] += Ac_offset1 / Acon0[i1] * ddel1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Update instance of fix gran/wall
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixGranularMDR::update_fix_gran_wall()
|
||||
{
|
||||
int i, m, nc, iwall;
|
||||
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *Acon0 = atom->dvector[index_Acon0];
|
||||
double *ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
|
||||
for (int w = 0; w < fix_wall_list.size(); w++) {
|
||||
FixWallGranRegion *fix = dynamic_cast<FixWallGranRegion*>(fix_wall_list[w]);
|
||||
GranularModel *model = fix->model;
|
||||
const int size_history = model->size_history;
|
||||
Region *region = fix->region;
|
||||
|
||||
if (region->dynamic_check())
|
||||
region->prematch();
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
if (! region->match(x[i][0], x[i][1], x[i][2])) continue;
|
||||
|
||||
nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i] + model->pulloff_distance(radius[i], 0.0));
|
||||
|
||||
if (nc == 0) {
|
||||
fix->ncontact[i] = 0;
|
||||
continue;
|
||||
}
|
||||
if (nc == 1) {
|
||||
fix->c2r[0] = 0;
|
||||
iwall = region->contact[0].iwall;
|
||||
if (fix->ncontact[i] == 0) {
|
||||
fix->ncontact[i] = 1;
|
||||
fix->walls[i][0] = iwall;
|
||||
for (m = 0; m < size_history; m++) fix->history_many[i][0] [m] = 0.0;
|
||||
} else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0])
|
||||
fix->update_contacts(i, nc);
|
||||
} else
|
||||
fix->update_contacts(i, nc);
|
||||
|
||||
// process current contacts
|
||||
for (int ic = 0; ic < nc; ic++) {
|
||||
const double wij = 1.0;
|
||||
if (Acon0[i] != 0.0) {
|
||||
const double delta = radius[i] - region->contact[ic].r;
|
||||
const double delta_offset0 = fix->history_many[i][fix->c2r[ic]][0];
|
||||
const double ddelta = delta - delta_offset0;
|
||||
const double Ac_offset0 = fix->history_many[i][fix->c2r[ic]][18];
|
||||
ddelta_bar[i] += wij * Ac_offset0 / Acon0[i] * ddelta;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
107
src/GRANULAR/fix_granular_mdr.h
Normal file
107
src/GRANULAR/fix_granular_mdr.h
Normal file
@ -0,0 +1,107 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(GRANULAR/MDR,FixGranularMDR);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_GRANULAR_MDR_H
|
||||
#define LMP_FIX_GRANULAR_MDR_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
namespace Granular_MDR_NS {
|
||||
|
||||
enum HistoryIndex {
|
||||
DELTA_0 = 0, // apparent overlap
|
||||
DELTA_1,
|
||||
DELTAO_0, // displacement
|
||||
DELTAO_1,
|
||||
DELTA_MDR_0, // MDR apparent overlap
|
||||
DELTA_MDR_1,
|
||||
DELTA_BULK_0, // bulk displacement
|
||||
DELTA_BULK_1,
|
||||
DELTAMAX_MDR_0, // maximum MDR apparent overlap
|
||||
DELTAMAX_MDR_1,
|
||||
YFLAG_0, // yield flag
|
||||
YFLAG_1,
|
||||
DELTAY_0, // yield displacement
|
||||
DELTAY_1,
|
||||
CA_0, // contact area intercept
|
||||
CA_1,
|
||||
AADH_0, // adhesive contact radius
|
||||
AADH_1,
|
||||
AC_0, // contact area
|
||||
AC_1,
|
||||
EPS_BAR_0, // volume-averaged infinitesimal sor
|
||||
EPS_BAR_1,
|
||||
PENALTY, // contact penalty
|
||||
DELTA_MAX,
|
||||
DELTAP_0,
|
||||
DELTAP_1
|
||||
};
|
||||
|
||||
} // namespace Granular_MDR_NS
|
||||
|
||||
class FixGranularMDR : public Fix {
|
||||
public:
|
||||
FixGranularMDR(class LAMMPS *, int, char **);
|
||||
~FixGranularMDR() override;
|
||||
int setmask() override;
|
||||
void post_constructor() override;
|
||||
void setup_pre_force(int) override;
|
||||
void pre_force(int) override;
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
void set_arrays(int) override;
|
||||
|
||||
private:
|
||||
int comm_stage;
|
||||
char *id_fix;
|
||||
double psi_b_coeff;
|
||||
class PairGranular *pair;
|
||||
class FixNeighHistory *fix_history;
|
||||
std::vector<Fix *> fix_wall_list;
|
||||
|
||||
void mean_surf_disp();
|
||||
void calculate_contact_penalty();
|
||||
void update_fix_gran_wall();
|
||||
|
||||
int index_Ro; // initial radius
|
||||
int index_Vgeo; // geometric particle volume of apparent particle afterremoving spherical cap volume
|
||||
int index_Velas; // particle volume from linear elasticity
|
||||
int index_Vcaps; // spherical cap volume from intersection of apparentradius particle and contact planes
|
||||
int index_eps_bar; // volume-averaged infinitesimal strain tensor
|
||||
int index_dRnumerator; // summation of numerator terms in calculation of dR
|
||||
int index_dRdenominator; // summation of denominator terms in calculation of dR
|
||||
int index_Acon0; // total area involved in contacts: Acon^{n}
|
||||
int index_Acon1; // total area involved in contacts: Acon^{n+1}
|
||||
int index_Atot; // total particle area
|
||||
int index_Atot_sum; // running sum of contact area minus cap area
|
||||
int index_ddelta_bar; // change in mean surface displacement
|
||||
int index_psi; // ratio of free surface area to total surface area
|
||||
int index_sigmaxx; // xx-component of the stress tensor, not necessary forforce calculation
|
||||
int index_sigmayy; // yy-component of the stress tensor, not necessary forforce calculation
|
||||
int index_sigmazz; // zz-component of the stress tensor, not necessary forforce calculation
|
||||
int index_history_setup_flag; // flag to check if history variables have beeninitialized
|
||||
int index_contacts; // total contacts on particle
|
||||
int index_adhesive_length; // total length of adhesive contact on a particle
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
@ -200,7 +200,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"contacts") == 0) {
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = 8;
|
||||
size_peratom_cols = 8 + model->nsvector;
|
||||
peratom_freq = 1;
|
||||
iarg += 1;
|
||||
} else if (strcmp(arg[iarg],"temperature") == 0) {
|
||||
@ -365,7 +365,7 @@ void FixWallGran::setup(int vflag)
|
||||
|
||||
void FixWallGran::post_force(int /*vflag*/)
|
||||
{
|
||||
int i,j;
|
||||
int i,j,n;
|
||||
double dx,dy,dz,del1,del2,delxy,delr,rwall,meff;
|
||||
double *forces, *torquesi;
|
||||
double vwall[3];
|
||||
@ -437,7 +437,9 @@ void FixWallGran::post_force(int /*vflag*/)
|
||||
|
||||
rwall = 0.0;
|
||||
|
||||
model->calculate_svector = 0;
|
||||
if (peratom_flag) {
|
||||
model->calculate_svector = 1;
|
||||
clear_stored_contacts();
|
||||
}
|
||||
|
||||
@ -546,6 +548,9 @@ void FixWallGran::post_force(int /*vflag*/)
|
||||
array_atom[i][5] = x[i][1] - dy;
|
||||
array_atom[i][6] = x[i][2] - dz;
|
||||
array_atom[i][7] = radius[i];
|
||||
|
||||
for (n = 0; n < model->nsvector; n++)
|
||||
array_atom[i][8 + n] = model->svector[n];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -50,14 +50,14 @@ class FixWallGran : public Fix {
|
||||
int maxsize_restart() override;
|
||||
void reset_dt() override;
|
||||
|
||||
// for granular model choices
|
||||
class Granular_NS::GranularModel *model;
|
||||
|
||||
protected:
|
||||
int wallstyle, wiggle, wshear, axis;
|
||||
int nlevels_respa;
|
||||
bigint time_origin;
|
||||
|
||||
// for granular model choices
|
||||
class Granular_NS::GranularModel *model;
|
||||
|
||||
double lo, hi, cylradius;
|
||||
double amplitude, period, omega, vshear;
|
||||
double dt;
|
||||
@ -84,7 +84,7 @@ class FixWallGran : public Fix {
|
||||
|
||||
// store particle interactions
|
||||
|
||||
int store;
|
||||
int nsvector;
|
||||
|
||||
void clear_stored_contacts();
|
||||
};
|
||||
|
||||
@ -118,7 +118,7 @@ void FixWallGranRegion::init()
|
||||
|
||||
void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
{
|
||||
int i, m, nc, iwall;
|
||||
int i, n, m, nc, iwall;
|
||||
double *forces, *torquesi;
|
||||
double meff, vwall[3], w0[3] = {0.0, 0.0, 0.0};
|
||||
bool touchflag = false;
|
||||
@ -174,7 +174,11 @@ void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
region->set_velocity();
|
||||
}
|
||||
|
||||
if (peratom_flag) clear_stored_contacts();
|
||||
model->calculate_svector = 0;
|
||||
if (peratom_flag) {
|
||||
model->calculate_svector = 1;
|
||||
clear_stored_contacts();
|
||||
}
|
||||
|
||||
// Define constant wall properties (atom j)
|
||||
model->radj = 0.0;
|
||||
@ -228,6 +232,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
model->radi = radius[i];
|
||||
model->radj = region->contact[ic].radius;
|
||||
model->r = region->contact[ic].r;
|
||||
model->i = i;
|
||||
model->j = ic;
|
||||
|
||||
if (model->beyond_contact) model->touch = history_many[i][c2r[ic]][0];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
@ -280,6 +287,9 @@ void FixWallGranRegion::post_force(int /*vflag*/)
|
||||
array_atom[i][5] = x[i][1] - model->dx[1];
|
||||
array_atom[i][6] = x[i][2] - model->dx[2];
|
||||
array_atom[i][7] = radius[i];
|
||||
|
||||
for (n = 0; n < model->nsvector; n++)
|
||||
array_atom[i][8 + n] = model->svector[n];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -44,9 +44,8 @@ class FixWallGranRegion : public FixWallGran {
|
||||
int size_restart(int) override;
|
||||
int maxsize_restart() override;
|
||||
|
||||
private:
|
||||
class Region *region;
|
||||
int nregion;
|
||||
void update_contacts(int, int);
|
||||
|
||||
// shear history for multiple contacts per particle
|
||||
|
||||
@ -57,10 +56,11 @@ class FixWallGranRegion : public FixWallGran {
|
||||
int *c2r; // contact to region mapping
|
||||
// c2r[i] = index of Ith contact in
|
||||
// region-contact[] list of contacts
|
||||
private:
|
||||
|
||||
int nregion;
|
||||
int motion_resetflag; // used by restart to indicate that region
|
||||
// vel info is to be reset
|
||||
|
||||
void update_contacts(int, int);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -42,6 +42,7 @@ GranSubMod::GranSubMod(class GranularModel *gm, LAMMPS *lmp) : Pointers(lmp)
|
||||
beyond_contact = 0;
|
||||
num_coeffs = 0;
|
||||
contact_radius_flag = 0;
|
||||
nsvector = 0;
|
||||
|
||||
nondefault_history_transfer = 0;
|
||||
transfer_history_factor = nullptr;
|
||||
|
||||
@ -46,6 +46,8 @@ namespace Granular_NS {
|
||||
|
||||
GranularModel *gm;
|
||||
|
||||
int nsvector, index_svector;
|
||||
|
||||
protected:
|
||||
int allocated;
|
||||
|
||||
|
||||
@ -50,6 +50,7 @@ GranSubModHeatRadius::GranSubModHeatRadius(GranularModel *gm, LAMMPS *lmp) : Gra
|
||||
num_coeffs = 1;
|
||||
contact_radius_flag = 1;
|
||||
conductivity = 0.0;
|
||||
nsvector = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -65,7 +66,9 @@ void GranSubModHeatRadius::coeffs_to_local()
|
||||
|
||||
double GranSubModHeatRadius::calculate_heat()
|
||||
{
|
||||
return 2 * conductivity * gm->contact_radius * (gm->Tj - gm->Ti);
|
||||
double heat = 2 * conductivity * gm->contact_radius * (gm->Tj - gm->Ti);
|
||||
if (gm->calculate_svector) gm->svector[index_svector] = heat;
|
||||
return heat;
|
||||
}
|
||||
|
||||
|
||||
@ -78,6 +81,7 @@ GranSubModHeatArea::GranSubModHeatArea(GranularModel *gm, LAMMPS *lmp) : GranSub
|
||||
num_coeffs = 1;
|
||||
contact_radius_flag = 1;
|
||||
heat_transfer_coeff = 0.0;
|
||||
nsvector = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -93,5 +97,7 @@ void GranSubModHeatArea::coeffs_to_local()
|
||||
|
||||
double GranSubModHeatArea::calculate_heat()
|
||||
{
|
||||
return heat_transfer_coeff * MY_PI * gm->contact_radius * gm->contact_radius * (gm->Tj - gm->Ti);
|
||||
double heat = heat_transfer_coeff * MY_PI * gm->contact_radius * gm->contact_radius * (gm->Tj - gm->Ti);
|
||||
if (gm->calculate_svector) gm->svector[index_svector] = heat;
|
||||
return heat;
|
||||
}
|
||||
|
||||
@ -13,25 +13,70 @@
|
||||
|
||||
#include "gran_sub_mod_normal.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "error.h"
|
||||
#include "citeme.h"
|
||||
#include "fix_granular_mdr.h"
|
||||
#include "granular_model.h"
|
||||
#include "math_const.h"
|
||||
#include "modify.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <iomanip>
|
||||
#include <sstream>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace Granular_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
using MathConst::MY_2PI;
|
||||
using MathConst::MY_PI;
|
||||
|
||||
static constexpr double PI27SQ = 266.47931882941264802866; // 27*PI**2
|
||||
static constexpr double PISQ = 9.8696044010893579923; // PI^2
|
||||
static constexpr double PIINV = 0.318309886183790691216; // 1/PI
|
||||
static constexpr double PI27SQ = 266.479318829412648029; // 27*PI^2
|
||||
static constexpr double PITOFIVETHIRDS = 6.73880859569814116838; // PI^(5/3)
|
||||
static constexpr double CBRT2 = 1.25992104989487319067; // cbrt(2)
|
||||
static constexpr double SQRTHALFPI = 1.25331413731550012081; // sqrt(PI/2)
|
||||
static constexpr double CBRTHALFPI = 1.16244735150962652526; // cbrt(PI/2)
|
||||
static constexpr double FOURTHIRDS = 1.33333333333333333333; // 4/3
|
||||
static constexpr double THREEROOT3 = 5.19615242270663202362; // 3*sqrt(3)
|
||||
static constexpr double SIXROOT6 = 14.69693845669906728801; // 6*sqrt(6)
|
||||
static constexpr double INVROOT6 = 0.40824829046386307274; // 1/sqrt(6)
|
||||
static constexpr double FOURTHIRDS = (4.0 / 3.0); // 4/3
|
||||
static constexpr double JKRPREFIX = 1.2277228507842888; // cbrt(3*PI**2/16)
|
||||
|
||||
static constexpr int MDR_MAX_IT = 100; // Newton-Raphson for MDR
|
||||
static constexpr double MDR_EPSILON1 = 1e-10; // Newton-Raphson for MDR
|
||||
static constexpr double MDR_EPSILON2 = 1e-16; // Newton-Raphson for MDR
|
||||
static constexpr double MDR_EPSILON3 = 1e-20; // For precision checks
|
||||
static constexpr double MDR_OVERLAP_LIMIT = 0.75; // Maximum contact overlap for MDR
|
||||
|
||||
static const char cite_mdr[] =
|
||||
"MDR contact model command: (i) https://doi.org/10.1016/j.jmps.2023.105492 || (ii) https://doi.org/10.1016/j.jmps.2023.105493 || (iii) https://doi.org/10.31224/4289\n\n"
|
||||
"@Article{zunker2024mechanicallyI,\n"
|
||||
" author = {Zunker, William and Kamrin, Ken},\n"
|
||||
" title = {A mechanically-derived contact model for adhesive elastic-perfectly plastic particles,\n"
|
||||
" Part I: Utilizing the method of dimensionality reduction},\n"
|
||||
" journal = {Journal of the Mechanics and Physics of Solids},\n"
|
||||
" year = {2024},\n"
|
||||
" volume = {183},\n"
|
||||
" pages = {105492},\n"
|
||||
"}\n\n"
|
||||
"@Article{zunker2024mechanicallyII,\n"
|
||||
" author = {Zunker, William and Kamrin, Ken},\n"
|
||||
" title = {A mechanically-derived contact model for adhesive elastic-perfectly plastic particles,\n"
|
||||
" Part II: Contact under high compaction—modeling a bulk elastic response},\n"
|
||||
" journal = {Journal of the Mechanics and Physics of Solids},\n"
|
||||
" year = {2024},\n"
|
||||
" volume = {183},\n"
|
||||
" pages = {105493},\n"
|
||||
"}\n\n"
|
||||
"@Article{zunker2025experimentally,\n"
|
||||
" author = {Zunker, William and Dunatunga, Sachith and Thakur, Subhash and Tang, Pingjun and Kamrin, Ken},\n"
|
||||
" title = {Experimentally validated DEM for large deformation powder compaction:\n"
|
||||
" mechanically-derived contact model and screening of non-physical contacts},\n"
|
||||
" year = {2025},\n"
|
||||
" journal = {engrXiv},\n"
|
||||
"}\n\n";
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Default normal model
|
||||
------------------------------------------------------------------------- */
|
||||
@ -381,3 +426,574 @@ void GranSubModNormalJKR::set_fncrit()
|
||||
{
|
||||
Fncrit = fabs(Fne + 2.0 * F_pulloff);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
MDR contact model
|
||||
|
||||
Contributing authors:
|
||||
William Zunker (MIT), Sachith Dunatunga (MIT),
|
||||
Dan Bolintineanu (SNL), Joel Clemmer (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
|
||||
GranSubModNormal(gm, lmp)
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_mdr);
|
||||
|
||||
num_coeffs = 6;
|
||||
contact_radius_flag = 1;
|
||||
size_history = 26;
|
||||
nsvector = 1;
|
||||
fix_mdr_flag = 0;
|
||||
id_fix = nullptr;
|
||||
|
||||
nondefault_history_transfer = 1;
|
||||
transfer_history_factor = new double[size_history];
|
||||
for (int i = 0; i < size_history; i++) {
|
||||
transfer_history_factor[i] = +1;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
GranSubModNormalMDR::~GranSubModNormalMDR()
|
||||
{
|
||||
if (id_fix && modify->nfix)
|
||||
modify->delete_fix(id_fix);
|
||||
delete[] id_fix;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void GranSubModNormalMDR::coeffs_to_local()
|
||||
{
|
||||
E = coeffs[0]; // Young's modulus
|
||||
nu = coeffs[1]; // Poisson's ratio
|
||||
Y = coeffs[2]; // yield stress
|
||||
gamma = coeffs[3]; // effective surface energy
|
||||
psi_b = coeffs[4]; // bulk response trigger based on ratio of remaining free area: A_{free}/A_{total}
|
||||
CoR = coeffs[5]; // coefficent of restitution
|
||||
|
||||
if (E <= 0.0) error->all(FLERR, "Illegal MDR normal model, Young's modulus must be greater than 0");
|
||||
if (nu < 0.0 || nu > 0.5) error->all(FLERR, "Illegal MDR normal model, Poisson's ratio must be between 0 and 0.5");
|
||||
if (Y < 0.0) error->all(FLERR, "Illegal MDR normal model, yield stress must be greater than or equal to 0");
|
||||
if (gamma < 0.0) error->all(FLERR, "Illegal MDR normal model, effective surface energy must be greater than or equal to 0");
|
||||
if (psi_b < 0.0 || psi_b > 1.0) error->all(FLERR, "Illegal MDR normal model, psi_b must be between 0 and 1.0");
|
||||
if (CoR < 0.0 || CoR > 1.0) error->all(FLERR, "Illegal MDR normal model, coefficent of restitution must be between 0 and 1.0");
|
||||
|
||||
G = E / (2.0 * (1.0 + nu)); // shear modulus
|
||||
kappa = E / (3.0 * (1.0 - 2.0 * nu)); // bulk modulus
|
||||
Eeff = E / (1.0 - pow(nu, 2.0)); // composite plane strain modulus
|
||||
|
||||
// precomputing factors
|
||||
|
||||
Eeffinv = 1.0 / Eeff;
|
||||
Eeffsq = Eeff * Eeff;
|
||||
Eeffsqinv = Eeffinv * Eeffinv;
|
||||
|
||||
gammasq = gamma * gamma;
|
||||
gamma3 = gammasq * gamma;
|
||||
gamma4 = gammasq * gammasq;
|
||||
|
||||
warn_flag = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void GranSubModNormalMDR::init()
|
||||
{
|
||||
if (!fix_mdr_flag) {
|
||||
if (modify->get_fix_by_style("GRANULAR/MDR").size() == 0) {
|
||||
id_fix = utils::strdup("MDR");
|
||||
modify->add_fix(fmt::format("{} all GRANULAR/MDR", id_fix));
|
||||
}
|
||||
fix_mdr_flag = 1;
|
||||
}
|
||||
|
||||
// initialize particle history variables
|
||||
int tmp1, tmp2;
|
||||
index_Ro = atom->find_custom("Ro", tmp1, tmp2); // initial radius
|
||||
index_Vcaps = atom->find_custom("Vcaps", tmp1, tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
|
||||
index_Vgeo = atom->find_custom("Vgeo", tmp1, tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
|
||||
index_Velas = atom->find_custom("Velas", tmp1, tmp2); // particle volume from linear elasticity
|
||||
index_eps_bar = atom->find_custom("eps_bar", tmp1, tmp2); // volume-averaged infinitesimal strain tensor
|
||||
index_dRnumerator = atom->find_custom("dRnumerator", tmp1, tmp2); // summation of numerator terms in calculation of dR
|
||||
index_dRdenominator = atom->find_custom("dRdenominator", tmp1, tmp2); // summation of denominator terms in calculation of dR
|
||||
index_Acon0 = atom->find_custom("Acon0", tmp1, tmp2); // total area involved in contacts: Acon^{n}
|
||||
index_Acon1 = atom->find_custom("Acon1", tmp1, tmp2); // total area involved in contacts: Acon^{n+1}
|
||||
index_Atot = atom->find_custom("Atot", tmp1, tmp2); // total particle area
|
||||
index_Atot_sum = atom->find_custom("Atot_sum", tmp1, tmp2); // running sum of contact area minus cap area
|
||||
index_ddelta_bar = atom->find_custom("ddelta_bar", tmp1, tmp2); // change in mean surface displacement
|
||||
index_psi = atom->find_custom("psi", tmp1, tmp2); // ratio of free surface area to total surface area
|
||||
index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); // xx-component of the stress tensor, not necessary for force calculation
|
||||
index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); // yy-component of the stress tensor, not necessary for force calculation
|
||||
index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); // zz-component of the stress tensor, not necessary for force calculation
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double GranSubModNormalMDR::calculate_forces()
|
||||
{
|
||||
using namespace Granular_MDR_NS;
|
||||
// To understand the structure of the overall code it is important to consider
|
||||
// the following:
|
||||
//
|
||||
// The MDR contact model was developed by imagining individual particles being
|
||||
// squished between a number of rigid flats (references below). To allow
|
||||
// for many interacting particles, we extend the idea of isolated particles surrounded
|
||||
// by rigid flats. In particular, we imagine placing rigid flats at the overlaps
|
||||
// between particles. The force is calculated seperately on both sides
|
||||
// of the contact assuming interaction with a rigid flat. The two forces are then
|
||||
// averaged on either side of the contact to determine the final force. If the
|
||||
// contact is between a particle and wall then only one force evaluation is required.
|
||||
//
|
||||
// Zunker and Kamrin, 2024, Part I: https://doi.org/10.1016/j.jmps.2023.105492
|
||||
// Zunker and Kamrin, 2024, Part II: https://doi.org/10.1016/j.jmps.2023.105493
|
||||
// Zunker, Dunatunga, Thakur, Tang, and Kamrin, 2025:
|
||||
|
||||
double *Rinitial = atom->dvector[index_Ro];
|
||||
double *Vgeo = atom->dvector[index_Vgeo];
|
||||
double *Velas = atom->dvector[index_Velas];
|
||||
double *Vcaps = atom->dvector[index_Vcaps];
|
||||
double *eps_bar = atom->dvector[index_eps_bar];
|
||||
double *dRnumerator = atom->dvector[index_dRnumerator];
|
||||
double *dRdenominator = atom->dvector[index_dRdenominator];
|
||||
double *Acon0 = atom->dvector[index_Acon0];
|
||||
double *Acon1 = atom->dvector[index_Acon1];
|
||||
double *Atot = atom->dvector[index_Atot];
|
||||
double *Atot_sum = atom->dvector[index_Atot_sum];
|
||||
double *ddelta_bar = atom->dvector[index_ddelta_bar];
|
||||
double *psi = atom->dvector[index_psi];
|
||||
double *sigmaxx = atom->dvector[index_sigmaxx];
|
||||
double *sigmayy = atom->dvector[index_sigmayy];
|
||||
double *sigmazz = atom->dvector[index_sigmazz];
|
||||
|
||||
const int itag_true = atom->tag[gm->i]; // true i particle tag
|
||||
const int jtag_true = atom->tag[gm->j]; // true j particle tag
|
||||
const int i_true = gm->i; // true i particle index
|
||||
const int j_true = gm->j; // true j particle index
|
||||
const double radi_true = gm->radi; // true i particle initial radius
|
||||
const double radj_true = gm->radj; // true j particle initial radius
|
||||
|
||||
double F = 0.0; // average force
|
||||
double F0 = 0.0; // force on contact side 0
|
||||
double F1 = 0.0; // force on contact side 1
|
||||
double delta = gm->delta; // apparent overlap
|
||||
double Ac_avg = 0.0; // average contact area across both sides
|
||||
|
||||
double *history = & gm->history[history_index]; // load in all history variables
|
||||
int history_update = gm->history_update;
|
||||
|
||||
// Rigid flat placement scheme
|
||||
double *deltamax_offset = & history[DELTA_MAX];
|
||||
double *deltap_offset0 = & history[DELTAP_0];
|
||||
double *deltap_offset1 = & history[DELTAP_1];
|
||||
double deltap0 = *deltap_offset0;
|
||||
double deltap1 = *deltap_offset1;
|
||||
|
||||
// always update deltamax since gm->delta won't change until initial integrate
|
||||
// also need to define deltamax if an atom is created with an overlap
|
||||
double deltamaxi = *deltamax_offset;
|
||||
if (gm->delta >= *deltamax_offset) *deltamax_offset = gm->delta;
|
||||
double deltamax = *deltamax_offset;
|
||||
|
||||
|
||||
for (int contactSide = 0; contactSide < 2; contactSide++) {
|
||||
|
||||
double *delta_offset, *deltao_offset, *delta_MDR_offset, *delta_BULK_offset;
|
||||
double *deltamax_MDR_offset, *Yflag_offset, *deltaY_offset, *cA_offset, *aAdh_offset;
|
||||
double *Ac_offset, *eps_bar_offset, *penalty_offset, *deltap_offset;
|
||||
|
||||
if (gm->contact_type == PAIR) {
|
||||
// displacement partitioning only necessary for particle-particle contact
|
||||
|
||||
// itag and jtag persist after neighbor list builds, use tags to compare to match
|
||||
// contact history variables consistently across steps for a particle pair.
|
||||
if ((contactSide == 0 && itag_true > jtag_true) || (contactSide != 0 && itag_true < jtag_true)) {
|
||||
gm->i = i_true;
|
||||
gm->j = j_true;
|
||||
gm->radi = radi_true;
|
||||
gm->radj = radj_true;
|
||||
} else {
|
||||
gm->i = j_true;
|
||||
gm->j = i_true;
|
||||
gm->radi = radj_true;
|
||||
gm->radj = radi_true;
|
||||
}
|
||||
|
||||
// determine the two maximum experienced geometric overlaps on either side of rigid flat
|
||||
double delta_geo, delta_geo_alt;
|
||||
double denom = 1.0 / (2.0 * (deltamax - gm->radi - gm->radj));
|
||||
double delta_geoOpt1 = deltamax * (deltamax - 2.0 * gm->radj) * denom;
|
||||
double delta_geoOpt2 = deltamax * (deltamax - 2.0 * gm->radi) * denom;
|
||||
if (gm->radi < gm->radj) {
|
||||
delta_geo = MAX(delta_geoOpt1, delta_geoOpt2);
|
||||
delta_geo_alt = MIN(delta_geoOpt1,delta_geoOpt2);
|
||||
} else {
|
||||
delta_geo = MIN(delta_geoOpt1, delta_geoOpt2);
|
||||
delta_geo_alt = MAX(delta_geoOpt1, delta_geoOpt2);
|
||||
}
|
||||
|
||||
// cap displacement if exceeds the overlap limit, parition the remaining to the other side
|
||||
if (delta_geo / gm->radi > MDR_OVERLAP_LIMIT) {
|
||||
delta_geo = gm->radi * MDR_OVERLAP_LIMIT;
|
||||
} else if (delta_geo_alt / gm->radj > MDR_OVERLAP_LIMIT) {
|
||||
delta_geo = deltamax - gm->radj * MDR_OVERLAP_LIMIT;
|
||||
}
|
||||
|
||||
// determine final delta used for subsequent calculations
|
||||
double deltap = deltap0 + deltap1;
|
||||
if (contactSide == 0) {
|
||||
delta = delta_geo + (deltap0 - delta_geo) * (gm->delta - deltamax) / (deltap - deltamax);
|
||||
} else {
|
||||
delta = delta_geo + (deltap1 - delta_geo) * (gm->delta - deltamax) / (deltap - deltamax);
|
||||
}
|
||||
} else if (gm->contact_type != PAIR && contactSide != 0) {
|
||||
// contact with particle-wall requires only one evaluation
|
||||
break;
|
||||
}
|
||||
|
||||
delta_offset = &history[DELTA_0 + contactSide];
|
||||
deltao_offset = &history[DELTAO_0 + contactSide];
|
||||
delta_MDR_offset = &history[DELTA_MDR_0 + contactSide];
|
||||
delta_BULK_offset = &history[DELTA_BULK_0 + contactSide];
|
||||
deltamax_MDR_offset = &history[DELTAMAX_MDR_0 + contactSide];
|
||||
Yflag_offset = &history[YFLAG_0 + contactSide];
|
||||
deltaY_offset = &history[DELTAY_0 + contactSide];
|
||||
cA_offset = &history[CA_0 + contactSide];
|
||||
aAdh_offset = &history[AADH_0 + contactSide];
|
||||
Ac_offset = &history[AC_0 + contactSide];
|
||||
eps_bar_offset = &history[EPS_BAR_0 + contactSide];
|
||||
deltap_offset = &history[DELTAP_0 + contactSide];
|
||||
|
||||
// temporary i and j indices
|
||||
const int i = gm->i;
|
||||
const int j = gm->j;
|
||||
|
||||
// geometric property definitions
|
||||
const double Ro = Rinitial[i]; // initial radius
|
||||
const double R = gm->radi; // apparent radius
|
||||
|
||||
// kinematics
|
||||
const double ddelta = delta - *delta_offset;
|
||||
if (history_update) *delta_offset = delta;
|
||||
|
||||
const double deltao = delta - (R - Ro);
|
||||
const double ddeltao = deltao - *deltao_offset;
|
||||
if (history_update) *deltao_offset = deltao;
|
||||
|
||||
double ddelta_MDR, ddelta_BULK;
|
||||
if (psi[i] < psi_b) {
|
||||
// bulk response triggered, split displacement increment between MDR and BULK components
|
||||
ddelta_MDR = MIN(ddelta - ddelta_bar[i], delta - *delta_MDR_offset);
|
||||
ddelta_BULK = ddelta_bar[i];
|
||||
} else {
|
||||
// no bulk response, full displacement increment goes to the MDR component
|
||||
ddelta_BULK = 0.0;
|
||||
ddelta_MDR = ddelta;
|
||||
}
|
||||
|
||||
// calculate and update MDR/BULK displacements
|
||||
const double delta_MDR = *delta_MDR_offset + ddelta_MDR;
|
||||
if (history_update) *delta_MDR_offset = delta_MDR;
|
||||
const double delta_BULK = MAX(0.0, *delta_BULK_offset + ddelta_BULK);
|
||||
if (history_update) *delta_BULK_offset = delta_BULK;
|
||||
|
||||
if (delta_MDR > *deltamax_MDR_offset) *deltamax_MDR_offset = delta_MDR;
|
||||
const double deltamax_MDR = *deltamax_MDR_offset;
|
||||
|
||||
// average pressure along yield surface
|
||||
const double pY = Y * (1.75 * exp(-4.4 * deltamax_MDR / R) + 1.0);
|
||||
|
||||
if (*Yflag_offset == 0.0 && delta_MDR >= deltamax_MDR) {
|
||||
const double phertz = 4 * Eeff * sqrt(delta_MDR) / (3 * MY_PI * sqrt(R));
|
||||
if (!history_update && warn_flag && deltamaxi == 0 && phertz > pY) {
|
||||
error->warning(FLERR, "The newly inserted particles have pre-existing overlaps that "
|
||||
"have caused immediate plastic deformation. This could lead to "
|
||||
"non-physical results in the MDR model, as it handles some aspects "
|
||||
"related to plastic deformation incrementally.");
|
||||
warn_flag = 0;
|
||||
}
|
||||
if (history_update && phertz > pY) {
|
||||
*Yflag_offset = 1.0;
|
||||
*deltaY_offset = delta_MDR;
|
||||
*cA_offset = MY_PI * (pow(*deltaY_offset, 2) - *deltaY_offset * R);
|
||||
}
|
||||
}
|
||||
|
||||
// MDR force calculation
|
||||
double F_MDR;
|
||||
double A, Ainv; // height of elliptical indenter
|
||||
double B; // width of elliptical indenter
|
||||
double deltae1D; // transformed elastic displacement
|
||||
double deltaR; // displacement correction
|
||||
double amax, amaxsq; // maximum experienced contact radius
|
||||
const double cA = *cA_offset; // contact area intercept
|
||||
|
||||
if (*Yflag_offset == 0.0) {
|
||||
// elastic contact
|
||||
A = 4.0 * R;
|
||||
Ainv = 1.0 / A;
|
||||
B = 2.0 * R;
|
||||
deltae1D = delta_MDR;
|
||||
amax = sqrt(deltamax_MDR * R);
|
||||
} else {
|
||||
// plastic contact
|
||||
amax = sqrt(2.0 * deltamax_MDR * R - pow(deltamax_MDR, 2) + cA * PIINV);
|
||||
amaxsq = amax * amax;
|
||||
A = 4.0 * pY * Eeffinv * amax;
|
||||
Ainv = 1.0 / A;
|
||||
B = 2.0 * amax;
|
||||
|
||||
// maximum transformed elastic displacement
|
||||
const double deltae1Dmax = A * 0.5;
|
||||
|
||||
// force caused by full submersion of elliptical indenter to depth of A/2
|
||||
double Fmax = Eeff * (A * B * 0.25) * acos(1 - 2 * deltae1Dmax * Ainv);
|
||||
Fmax -= (2 - 4 * deltae1Dmax * Ainv) * sqrt(deltae1Dmax * Ainv - pow(deltae1Dmax * Ainv, 2));
|
||||
|
||||
// depth of particle center
|
||||
const double zR = R - (deltamax_MDR - deltae1Dmax);
|
||||
|
||||
deltaR = 2 * amaxsq * (-1 + nu) - (-1 + 2 * nu) * zR * (-zR + sqrt(amaxsq + pow(zR, 2)));
|
||||
deltaR *= Fmax / (MY_2PI * amaxsq * G * sqrt(amaxsq + pow(zR, 2)));
|
||||
|
||||
// transformed elastic displacement
|
||||
deltae1D = (delta_MDR - deltamax_MDR + deltae1Dmax + deltaR) / (1 + deltaR / deltae1Dmax);
|
||||
|
||||
// added for rigid flat placement
|
||||
if (history_update) *deltap_offset = deltamax_MDR - (deltae1Dmax + deltaR);
|
||||
}
|
||||
|
||||
double a_na;
|
||||
double a_fac = 0.99;
|
||||
(deltae1D >= 0.0) ? a_na = B * sqrt(A - deltae1D) * sqrt(deltae1D) * Ainv : a_na = 0.0;
|
||||
double aAdh = *aAdh_offset;
|
||||
if (aAdh > a_fac * amax) aAdh = a_fac * amax;
|
||||
|
||||
double Ainvsq = Ainv * Ainv;
|
||||
double Asq = A * A;
|
||||
double A3 = Asq * A;
|
||||
double A4 = Asq * Asq;
|
||||
|
||||
double Binv = 1.0 / B;
|
||||
double Bsq = B * B;
|
||||
double B4 = Bsq * Bsq;
|
||||
|
||||
if (gamma <= 0.0) {
|
||||
// non-adhesive contact
|
||||
|
||||
if (deltae1D <= 0.0) {
|
||||
F_MDR = 0.0;
|
||||
} else {
|
||||
F_MDR = calculate_nonadhesive_mdr_force(deltae1D, Ainv, Eeff, A, B);
|
||||
}
|
||||
|
||||
if (std::isnan(F_MDR)) {
|
||||
error->one(FLERR, "F_MDR is NaN, non-adhesive case");
|
||||
}
|
||||
|
||||
if (history_update) *aAdh_offset = a_na;
|
||||
} else {
|
||||
// adhesive contact
|
||||
double g_aAdh;
|
||||
|
||||
if (delta_MDR == deltamax_MDR || a_na >= aAdh) {
|
||||
// case 1: no tensile springs, purely compressive contact
|
||||
|
||||
if (deltae1D <= 0.0) {
|
||||
F_MDR = 0.0;
|
||||
} else {
|
||||
F_MDR = calculate_nonadhesive_mdr_force(deltae1D, Ainv, Eeff, A, B);
|
||||
}
|
||||
|
||||
if (std::isnan(F_MDR))
|
||||
error->one(FLERR, "F_MDR is NaN, case 1: no tensile springs");
|
||||
|
||||
if (history_update) *aAdh_offset = a_fac * a_na;
|
||||
} else {
|
||||
// case 2+3, tensile springs
|
||||
const double lmax = sqrt(MY_2PI * aAdh * gamma * Eeffinv);
|
||||
g_aAdh = A * 0.5 - A * Binv * sqrt(Bsq * 0.25 - pow(aAdh, 2));
|
||||
g_aAdh = round_up_negative_epsilon(g_aAdh);
|
||||
|
||||
double tmp = 27 * A4 * B4 * gamma * Eeffinv;
|
||||
tmp -= 2 * pow(B, 6) * gamma3 * PISQ * pow(Eeffinv, 3);
|
||||
tmp += sqrt(27) * Asq * B4 * sqrt(27 * A4 * Eeffsq * gammasq - 4 * Bsq * gamma4 * PISQ) * Eeffsqinv;
|
||||
tmp = cbrt(tmp);
|
||||
|
||||
double acrit = -Bsq * gamma * MY_PI * Ainvsq * Eeffinv;
|
||||
acrit += CBRT2 * B4 * gammasq * PITOFIVETHIRDS / (Asq * Eeffsq * tmp);
|
||||
acrit += CBRTHALFPI * tmp * Ainvsq;
|
||||
acrit /= 6;
|
||||
|
||||
if ((deltae1D + lmax - g_aAdh) >= 0.0) {
|
||||
// case 2: tensile springs do not exceed critical length --> deltae + lmax - g(aAdhes) >= 0
|
||||
const double deltaeAdh = g_aAdh;
|
||||
const double F_na = calculate_nonadhesive_mdr_force(deltaeAdh, Ainv, Eeff, A, B);
|
||||
const double F_Adhes = 2.0 * Eeff * (deltae1D - deltaeAdh) * aAdh;
|
||||
F_MDR = F_na + F_Adhes;
|
||||
if (std::isnan(F_MDR))
|
||||
error->one(FLERR, "F_MDR is NaN, case 2: tensile springs, but not exceeding critical length");
|
||||
} else {
|
||||
// case 3: tensile springs exceed critical length --> deltae + lmax - g(aAdhes) = 0
|
||||
|
||||
if (aAdh < acrit) {
|
||||
aAdh = 0.0;
|
||||
F_MDR = 0.0;
|
||||
} else {
|
||||
// newton-raphson to find aAdh
|
||||
double aAdh_tmp = aAdh;
|
||||
double fa, fa2, fa_tmp, dfda;
|
||||
for (int lv1 = 0; lv1 < MDR_MAX_IT; ++lv1) {
|
||||
fa_tmp = deltae1D - A * 0.5 + A * sqrt(Bsq * 0.25 - pow(aAdh_tmp, 2)) * Binv;
|
||||
fa = fa_tmp + sqrt(MY_2PI * aAdh_tmp * gamma * Eeffinv);
|
||||
if (abs(fa) < MDR_EPSILON1) {
|
||||
break;
|
||||
}
|
||||
dfda = -aAdh_tmp * A / (B * sqrt(-pow(aAdh_tmp, 2) + Bsq * 0.25));
|
||||
dfda += gamma * SQRTHALFPI / sqrt(aAdh_tmp * gamma * Eeff);
|
||||
aAdh_tmp = aAdh_tmp - fa / dfda;
|
||||
fa2 = fa_tmp + sqrt(MY_2PI * aAdh_tmp * gamma * Eeffinv);
|
||||
if (abs(fa - fa2) < MDR_EPSILON2) {
|
||||
break;
|
||||
}
|
||||
if (lv1 == MDR_MAX_IT - 1) {
|
||||
aAdh_tmp = 0.0;
|
||||
}
|
||||
}
|
||||
aAdh = aAdh_tmp;
|
||||
|
||||
g_aAdh = A * 0.5 - A * Binv * sqrt(Bsq * 0.25 - pow(aAdh, 2));
|
||||
g_aAdh = round_up_negative_epsilon(g_aAdh);
|
||||
|
||||
const double deltaeAdh = g_aAdh;
|
||||
const double F_na = calculate_nonadhesive_mdr_force(deltaeAdh, Ainv, Eeff, A, B);
|
||||
const double F_Adhes = 2.0 * Eeff * (deltae1D - deltaeAdh) * aAdh;
|
||||
F_MDR = F_na + F_Adhes;
|
||||
if (std::isnan(F_MDR))
|
||||
error->one(FLERR, "F_MDR is NaN, case 3: tensile springs exceed critical length");
|
||||
}
|
||||
if (history_update) *aAdh_offset = aAdh;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// contact penalty scheme
|
||||
penalty_offset = &history[PENALTY];
|
||||
double pij = *penalty_offset;
|
||||
const double wij = MAX(1.0 - pij, 0.0);
|
||||
|
||||
// area related calculations
|
||||
double Ac;
|
||||
(*Yflag_offset == 0.0) ? Ac = MY_PI * delta * R : Ac = MY_PI * (2.0 * delta * R - pow(delta, 2)) + cA;
|
||||
if (Ac < 0.0) Ac = 0.0;
|
||||
if (history_update) {
|
||||
Atot_sum[i] += wij * (Ac - MY_2PI * R * (deltamax_MDR + delta_BULK));
|
||||
Acon1[i] += wij * Ac;
|
||||
}
|
||||
Ac_avg += wij * Ac;
|
||||
|
||||
// bulk force calculation
|
||||
double F_BULK;
|
||||
(delta_BULK <= 0.0) ? F_BULK = 0.0 : F_BULK = (1.0 / Vgeo[i]) * Acon0[i] * delta_BULK * kappa * Ac;
|
||||
|
||||
// total force calculation
|
||||
(contactSide == 0) ? F0 = F_MDR + F_BULK : F1 = F_MDR + F_BULK;
|
||||
|
||||
if (history_update) {
|
||||
// mean surface displacement calculation
|
||||
*Ac_offset = wij * Ac;
|
||||
|
||||
// radius update scheme quantity calculation
|
||||
Vcaps[i] += MY_PI * THIRD * pow(delta, 2) * (3.0 * R - delta);
|
||||
}
|
||||
|
||||
const double Fntmp = wij * (F_MDR + F_BULK);
|
||||
const double fx = Fntmp * gm->nx[0];
|
||||
const double fy = Fntmp * gm->nx[1];
|
||||
const double fz = Fntmp * gm->nx[2];
|
||||
const double bx = -(Ro - deltao) * gm->nx[0];
|
||||
const double by = -(Ro - deltao) * gm->nx[1];
|
||||
const double bz = -(Ro - deltao) * gm->nx[2];
|
||||
const double eps_bar_contact = (fx * bx + fy * by + fz * bz) / (3 * kappa * Velas[i]);
|
||||
if (history_update) eps_bar[i] += eps_bar_contact;
|
||||
|
||||
double desp_bar_contact = eps_bar_contact - *eps_bar_offset;
|
||||
if (history_update && delta_MDR == deltamax_MDR && *Yflag_offset > 0.0 && F_MDR > 0.0) {
|
||||
const double Vo = FOURTHIRDS * MY_PI * pow(Ro, 3);
|
||||
dRnumerator[i] -= Vo * (eps_bar_contact - *eps_bar_offset);
|
||||
dRnumerator[i] -= wij * MY_PI * ddeltao * (2 * deltao * Ro - pow(deltao, 2) + pow(R, 2) - pow(Ro, 2));
|
||||
dRdenominator[i] += wij * 2.0 * MY_PI * R * (deltao + R - Ro);
|
||||
}
|
||||
|
||||
if (history_update) {
|
||||
*eps_bar_offset = eps_bar_contact;
|
||||
sigmaxx[i] += fx * bx / Velas[i];
|
||||
sigmayy[i] += fy * by / Velas[i];
|
||||
sigmazz[i] += fz * bz / Velas[i];
|
||||
}
|
||||
}
|
||||
|
||||
// save contact area
|
||||
if (gm->calculate_svector) gm->svector[index_svector] = Ac_avg * 0.5;
|
||||
|
||||
gm->i = i_true;
|
||||
gm->j = j_true;
|
||||
gm->radi = radi_true;
|
||||
gm->radj = radj_true;
|
||||
|
||||
double *penalty_offset = &history[PENALTY];
|
||||
const double pij = *penalty_offset;
|
||||
const double wij = MAX(1.0 - pij, 0.0);
|
||||
|
||||
// assign final force
|
||||
if (gm->contact_type != PAIR) {
|
||||
F = wij * F0;
|
||||
} else {
|
||||
F = wij * (F0 + F1) * 0.5;
|
||||
}
|
||||
|
||||
// calculate damping force
|
||||
if (F > 0.0) {
|
||||
double Eeff2;
|
||||
double Reff2;
|
||||
if (gm->contact_type == PAIR) {
|
||||
Eeff2 = E / (2.0 * (1.0 - pow(nu, 2)));
|
||||
Reff2 = 1.0 / ((1.0 / gm->radi + 1.0 / gm->radj));
|
||||
} else {
|
||||
Eeff2 = E / (1.0 - pow(nu, 2));
|
||||
Reff2 = gm->radi;
|
||||
}
|
||||
const double kn = Eeff2 * Reff2;
|
||||
const double beta = -log(CoR) / sqrt(pow(log(CoR), 2) + PISQ);
|
||||
const double damp_prefactor = beta * sqrt(gm->meff * kn);
|
||||
const double F_DAMP = -damp_prefactor * gm->vnnr;
|
||||
|
||||
F += wij * F_DAMP;
|
||||
}
|
||||
|
||||
return F;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double GranSubModNormalMDR::calculate_nonadhesive_mdr_force(double delta, double Ainv, double Eeff, double A, double B)
|
||||
{
|
||||
double F_na = acos(1.0 - 2.0 * delta * Ainv);
|
||||
F_na -= (2 - 4 * delta * Ainv) * sqrt(delta * Ainv - pow(delta * Ainv, 2));
|
||||
F_na *= 0.25 * Eeff * A * B;
|
||||
|
||||
return F_na;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
round values within (-EPSILON,0.0) due to machine precision errors to zero
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double GranSubModNormalMDR::round_up_negative_epsilon(double value)
|
||||
{
|
||||
if (value < 0.0 && value > -MDR_EPSILON3) value = 0.0;
|
||||
return value;
|
||||
}
|
||||
|
||||
@ -19,6 +19,7 @@ GranSubModStyle(hertz,GranSubModNormalHertz,NORMAL);
|
||||
GranSubModStyle(hertz/material,GranSubModNormalHertzMaterial,NORMAL);
|
||||
GranSubModStyle(dmt,GranSubModNormalDMT,NORMAL);
|
||||
GranSubModStyle(jkr,GranSubModNormalJKR,NORMAL);
|
||||
GranSubModStyle(mdr,GranSubModNormalMDR,NORMAL);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
@ -133,6 +134,35 @@ namespace Granular_NS {
|
||||
int mixed_coefficients;
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
class GranSubModNormalMDR : public GranSubModNormal {
|
||||
public:
|
||||
GranSubModNormalMDR(class GranularModel *, class LAMMPS *);
|
||||
~GranSubModNormalMDR() override;
|
||||
void coeffs_to_local() override;
|
||||
void init() override;
|
||||
double calculate_forces() override;
|
||||
double E, nu, Y, gamma, CoR, psi_b; // specified coeffs
|
||||
|
||||
protected:
|
||||
double G, kappa, Eeff; // derived coeffs
|
||||
double Eeffsq, Eeffinv, Eeffsqinv;
|
||||
double gammasq, gamma3, gamma4;
|
||||
|
||||
int warn_flag;
|
||||
|
||||
int index_Ro, index_Vgeo, index_Velas, index_Vcaps, index_eps_bar, index_dRnumerator;
|
||||
int index_dRdenominator, index_Acon0, index_Acon1, index_Atot, index_Atot_sum, index_ddelta_bar;
|
||||
int index_psi, index_sigmaxx, index_sigmayy, index_sigmazz, index_contacts, index_adhesive_length;
|
||||
int fix_mdr_flag;
|
||||
|
||||
char *id_fix;
|
||||
|
||||
inline double calculate_nonadhesive_mdr_force(double, double, double, double, double);
|
||||
inline double round_up_negative_epsilon(double);
|
||||
};
|
||||
|
||||
} // namespace Granular_NS
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
|
||||
@ -27,6 +27,7 @@
|
||||
#include "force.h"
|
||||
#include "gran_sub_mod.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
|
||||
#include "style_gran_sub_mod.h" // IWYU pragma: keep
|
||||
|
||||
@ -64,6 +65,10 @@ GranularModel::GranularModel(LAMMPS *lmp) : Pointers(lmp)
|
||||
twisting_model = nullptr;
|
||||
heat_model = nullptr;
|
||||
|
||||
calculate_svector = 0;
|
||||
nsvector = 0;
|
||||
svector = nullptr;
|
||||
|
||||
for (int i = 0; i < NSUBMODELS; i++) sub_models[i] = nullptr;
|
||||
transfer_history_factor = nullptr;
|
||||
|
||||
@ -100,6 +105,7 @@ GranularModel::~GranularModel()
|
||||
delete[] gran_sub_mod_class;
|
||||
delete[] gran_sub_mod_names;
|
||||
delete[] gran_sub_mod_types;
|
||||
delete[] svector;
|
||||
|
||||
for (int i = 0; i < NSUBMODELS; i++) delete sub_models[i];
|
||||
}
|
||||
@ -243,7 +249,12 @@ void GranularModel::init()
|
||||
|
||||
// Must have valid normal, damping, and tangential models
|
||||
if (normal_model->name == "none") error->all(FLERR, "Must specify normal granular model");
|
||||
if (normal_model->name == "mdr") {
|
||||
if (damping_model->name != "none")
|
||||
error->all(FLERR, "MDR require 'none' damping model. To damp, specify a coefficient of restitution < 1.");
|
||||
} else {
|
||||
if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model");
|
||||
}
|
||||
if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model");
|
||||
|
||||
// Twisting, rolling, and heat are optional
|
||||
@ -293,6 +304,21 @@ void GranularModel::init()
|
||||
}
|
||||
|
||||
for (int i = 0; i < NSUBMODELS; i++) sub_models[i]->init();
|
||||
|
||||
nsvector = 0;
|
||||
int index_svector = 0;
|
||||
for (int i = 0; i < NSUBMODELS; i++) {
|
||||
if (sub_models[i]->nsvector != 0) {
|
||||
sub_models[i]->index_svector = index_svector;
|
||||
nsvector += sub_models[i]->nsvector;
|
||||
index_svector += sub_models[i]->nsvector;
|
||||
}
|
||||
}
|
||||
|
||||
if (nsvector != 0) {
|
||||
delete[] svector;
|
||||
svector = new double[nsvector];
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -493,10 +519,9 @@ void GranularModel::calculate_forces()
|
||||
if (contact_type == PAIR) sub3(torquesj, tortwist, torquesj);
|
||||
}
|
||||
|
||||
if (heat_defined) {
|
||||
if (heat_defined)
|
||||
dq = heat_model->calculate_heat();
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute pull-off distance (beyond contact) for a given radius and atom type
|
||||
|
||||
@ -74,6 +74,9 @@ class GranularModel : protected Pointers {
|
||||
int beyond_contact, limit_damping, history_update;
|
||||
ContactType contact_type;
|
||||
|
||||
// Particle identifiers
|
||||
int i, j, itype, jtype;
|
||||
|
||||
// History variables
|
||||
int size_history, nondefault_history_transfer;
|
||||
double *transfer_history_factor;
|
||||
@ -93,6 +96,10 @@ class GranularModel : protected Pointers {
|
||||
double magtwist;
|
||||
bool touch;
|
||||
|
||||
// Extra output
|
||||
int calculate_svector, nsvector;
|
||||
double *svector;
|
||||
|
||||
protected:
|
||||
int rolling_defined, twisting_defined, heat_defined; // Flag optional sub models
|
||||
int classic_model; // Flag original pair/gran calculations
|
||||
|
||||
@ -199,6 +199,11 @@ void PairGranular::compute(int eflag, int vflag)
|
||||
model->xj = x[j];
|
||||
model->radi = radius[i];
|
||||
model->radj = radius[j];
|
||||
model->i = i;
|
||||
model->j = j;
|
||||
model->itype = itype;
|
||||
model->jtype = jtype;
|
||||
|
||||
if (use_history) model->touch = touch[jj];
|
||||
|
||||
touchflag = model->check_contact();
|
||||
@ -412,8 +417,10 @@ void PairGranular::init_style()
|
||||
error->all(FLERR,"Heat conduction in pair granular requires atom style with heatflow property");
|
||||
}
|
||||
|
||||
// allocate history and initialize models
|
||||
// allocate history and aggregate model information
|
||||
class GranularModel* model;
|
||||
double nsvector_total;
|
||||
extra_svector = 0;
|
||||
int size_max[NSUBMODELS] = {0};
|
||||
for (int n = 0; n < nmodels; n++) {
|
||||
model = models_list[n];
|
||||
@ -424,13 +431,23 @@ void PairGranular::init_style()
|
||||
}
|
||||
if (model->size_history != 0) use_history = 1;
|
||||
|
||||
for (int i = 0; i < NSUBMODELS; i++)
|
||||
nsvector_total = 0;
|
||||
for (int i = 0; i < NSUBMODELS; i++) {
|
||||
nsvector_total += model->sub_models[i]->nsvector;
|
||||
if (model->sub_models[i]->size_history > size_max[i])
|
||||
size_max[i] = model->sub_models[i]->size_history;
|
||||
}
|
||||
extra_svector = MAX(extra_svector, nsvector_total);
|
||||
|
||||
if (model->nondefault_history_transfer) nondefault_history_transfer = 1;
|
||||
}
|
||||
|
||||
if (extra_svector != 0) {
|
||||
single_extra = 12 + extra_svector;
|
||||
delete[] svector;
|
||||
svector = new double[single_extra];
|
||||
}
|
||||
|
||||
size_history = 0;
|
||||
if (use_history) {
|
||||
for (int i = 0; i < NSUBMODELS; i++) size_history += size_max[i];
|
||||
@ -711,6 +728,10 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
||||
model->xj = x[j];
|
||||
model->radi = radius[i];
|
||||
model->radj = radius[j];
|
||||
model->i = i;
|
||||
model->j = j;
|
||||
model->itype = itype;
|
||||
model->jtype = jtype;
|
||||
model->history_update = 0; // Don't update history
|
||||
|
||||
// If history is needed
|
||||
@ -765,7 +786,9 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
||||
model->omegaj = omega[j];
|
||||
model->history = history;
|
||||
|
||||
model->calculate_svector = 1;
|
||||
model->calculate_forces();
|
||||
model->calculate_svector = 0;
|
||||
|
||||
// apply forces & torques
|
||||
// Calculate normal component, normalized by r
|
||||
@ -785,6 +808,14 @@ double PairGranular::single(int i, int j, int itype, int jtype,
|
||||
svector[10] = model->dx[1];
|
||||
svector[11] = model->dx[2];
|
||||
|
||||
// add submodel-specific quantities
|
||||
for (int n = 0; n < model->nsvector; n++)
|
||||
svector[12 + n] = model->svector[n];
|
||||
|
||||
// zero any values unused by this specific model
|
||||
for (int n = 12 + model->nsvector; n < single_extra; n++)
|
||||
svector[n] = 0.0;
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
||||
@ -46,6 +46,12 @@ class PairGranular : public Pair {
|
||||
double memory_usage() override;
|
||||
double atom2cut(int) override;
|
||||
double radii2cut(double, double) override;
|
||||
int get_size_history() const { return size_history; }
|
||||
|
||||
// granular models
|
||||
class Granular_NS::GranularModel** models_list;
|
||||
int **types_indices;
|
||||
int nmodels, maxmodels;
|
||||
|
||||
protected:
|
||||
int freeze_group_bit;
|
||||
@ -73,14 +79,11 @@ class PairGranular : public Pair {
|
||||
int size_history;
|
||||
int heat_flag;
|
||||
|
||||
// granular models
|
||||
int nmodels, maxmodels;
|
||||
class Granular_NS::GranularModel **models_list;
|
||||
int **types_indices;
|
||||
|
||||
// optional user-specified global cutoff, per-type user-specified cutoffs
|
||||
double **cutoff_type;
|
||||
double cutoff_global;
|
||||
|
||||
int extra_svector;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -77,7 +77,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
|
||||
SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
|
||||
RANDOM,NORMAL,CEIL,FLOOR,ROUND,TERNARY,
|
||||
RAMP,STAGGER,LOGFREQ,LOGFREQ2,LOGFREQ3,STRIDE,STRIDE2,
|
||||
VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,
|
||||
VDISPLACE,SWIGGLE,CWIGGLE,SIGN,GMASK,RMASK,
|
||||
GRMASK,IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,IS_FILE,EXTRACT_SETTING,
|
||||
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY};
|
||||
|
||||
@ -2594,7 +2594,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
|
||||
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),ternary(x,y,z),
|
||||
ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
|
||||
logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),swiggle(x,y,z),
|
||||
cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
|
||||
cwiggle(x,y,z),sign(x),gmask(x),rmask(x),grmask(x,y)
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
double Variable::collapse_tree(Tree *tree)
|
||||
@ -3144,6 +3144,14 @@ double Variable::collapse_tree(Tree *tree)
|
||||
return tree->value;
|
||||
}
|
||||
|
||||
if (tree->type == SIGN) {
|
||||
arg1 = collapse_tree(tree->first);
|
||||
if (tree->first->type != VALUE) return 0.0;
|
||||
tree->type = VALUE;
|
||||
tree->value = (arg1 >= 0.0) ? 1.0 : -1.0; // sign(arg1);
|
||||
return tree->value;
|
||||
}
|
||||
|
||||
// mask functions do not become a single collapsed value
|
||||
|
||||
if (tree->type == GMASK) return 0.0;
|
||||
@ -3162,7 +3170,7 @@ double Variable::collapse_tree(Tree *tree)
|
||||
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),ternary(x,y,z),
|
||||
ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
|
||||
logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),
|
||||
swiggle(x,y,z),cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
|
||||
swiggle(x,y,z),cwiggle(x,y,z),sign(x),gmask(x),rmask(x),grmask(x,y)
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
double Variable::eval_tree(Tree *tree, int i)
|
||||
@ -3478,6 +3486,9 @@ double Variable::eval_tree(Tree *tree, int i)
|
||||
return arg;
|
||||
}
|
||||
|
||||
if (tree->type == SIGN)
|
||||
return (eval_tree(tree->first,i) >= 0.0) ? 1.0 : -1.0; // sign(eval_tree(tree->first,i));
|
||||
|
||||
if (tree->type == GMASK) {
|
||||
if (atom->mask[i] & tree->ivalue) return 1.0;
|
||||
else return 0.0;
|
||||
@ -3651,7 +3662,7 @@ tagint Variable::int_between_brackets(char *&ptr, int varallow)
|
||||
atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),ternary(),
|
||||
ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
|
||||
logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),
|
||||
swiggle(x,y,z),cwiggle(x,y,z)
|
||||
swiggle(x,y,z),cwiggle(x,y,z),sign(x)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int Variable::math_function(char *word, char *contents, Tree **tree, Tree **treestack,
|
||||
@ -3669,7 +3680,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree
|
||||
strcmp(word,"logfreq") != 0 && strcmp(word,"logfreq2") != 0 &&
|
||||
strcmp(word,"logfreq3") != 0 && strcmp(word,"stride") != 0 &&
|
||||
strcmp(word,"stride2") != 0 && strcmp(word,"vdisplace") != 0 &&
|
||||
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0)
|
||||
strcmp(word,"swiggle") != 0 && strcmp(word,"cwiggle") != 0 && strcmp(word,"sign") != 0)
|
||||
return 0;
|
||||
|
||||
// parse contents for comma-separated args
|
||||
@ -4064,6 +4075,11 @@ int Variable::math_function(char *word, char *contents, Tree **tree, Tree **tree
|
||||
double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
|
||||
argstack[nargstack++] = value;
|
||||
}
|
||||
} else if (strcmp(word,"sign") == 0) {
|
||||
if (narg != 1)
|
||||
print_var_error(FLERR,"Invalid math function in variable formula",ivar);
|
||||
if (tree) newtree->type = SIGN;
|
||||
else argstack[nargstack++] = (value1 >= 0.0) ? 1.0 : -1.0; // sign(value1);
|
||||
}
|
||||
|
||||
// delete stored args
|
||||
|
||||
Reference in New Issue
Block a user