diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 5c46ad892e..556d683c72 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -17,7 +17,7 @@ Syntax * zmin = minimal number of neighbors for reproducing kernels * zero or more keyword/value pairs may be appended to args * keyword = *thermal* or *interface/reconstruct* or *surface/detection* or - *shift* or *rho/sum* or *density* or *speed/sound* + *shift* or *rho/sum* or *density* or *self/mass* or *speed/sound* .. parsed-literal:: @@ -29,6 +29,7 @@ Syntax *limit/splash* = threshold for splash particles *shift* values = none, turns on velocity shifting *rho/sum* values = none, uses the kernel to compute the density of particles + *self/mass* values = none, a particle uses its own mass in a rho summation *density* values = *rho01*, ... *rho0N* (density) *speed/sound* values = *cs0*, ... *csN* (velocity) @@ -106,24 +107,38 @@ threshold for this classification is set by the numerical value of By default, RHEO integrates particles' densities using a mass diffusion equation. Alternatively, one can update densities every timestep by performing -a kernel summation of the masses of neighboring particles by specifying the *rho/sum* keyword. +a kernel summation of the masses of neighboring particles by specifying the *rho/sum* +keyword. -The *density* is used to specify the equilbrium density of each of the N +The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*. +Typically, the density :math:`\rho` of a particle is calculated as the sum + +.. math:: + \rho_i = \Sum_{j} W_{ij} M_j + +where the summation is over neighbors, :math:`W_{ij}` is the kernel, and :math:`M_j` +is the mass of particle :math:`j`. The *self/mass* keyword augments this expression +by replacing :math:`M_j` with :math:`M_i`. This may be useful in simulations of +multiple fluid phases with large differences in density, :ref:`(Hu) `. + +The *density* keyword is used to specify the equilbrium density of each of the N particle types. It must be followed by N numerical values specifying each type's equilibrium density *rho0*. -The *density* is used to specify the speed of sound of each of the N particle -types. It must be followed by N numerical values specifying each type's speed -of sound *cs*. +The *speed/sound* keyword is used to specify the speed of sound of each of the +N particle types. It must be followed by N numerical values specifying each +type's speed of sound *cs*. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options +No information about this fix is written to :doc:`binary restart files `. +None of the :doc:`fix_modify ` options are relevant to this fix. No global or per-atom quantities are stored by this fix for access by various :doc:`output commands `. No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. Restrictions """""""""""" @@ -138,7 +153,8 @@ set to all. Only one instance of fix rheo may be defined and it must be defined prior to all other RHEO fixes. This fix is part of the RHEO package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. Related commands """""""""""""""" @@ -156,6 +172,10 @@ Default ---------- -.. _howto-howto_rheo_palermo: +.. _howto_rheo_palermo: **(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation. + +.. _fix_rheo_hu: + +**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006). diff --git a/doc/src/fix_rheo_viscosity.rst b/doc/src/fix_rheo_viscosity.rst index 8c403f8d0b..8175178787 100644 --- a/doc/src/fix_rheo_viscosity.rst +++ b/doc/src/fix_rheo_viscosity.rst @@ -14,20 +14,26 @@ Syntax * rheo/viscosity = style name of this fix command * one or more types and viscosity styles must be appended * types = lists of types (see below) -* vstyle = *constant* +* vstyle = *constant* or *power* .. parsed-literal:: *constant* args = *eta* *eta* = viscosity + *power* args = *eta*, *gd0*, *K*, *n* + *eta* = viscosity + *gd0* = critical strain rate + *K* = consistency index + *n* = power-law exponent + Examples """""""" .. code-block:: LAMMPS fix 1 all rheo/viscosity * constant 1.0 - fix 1 all rheo/viscosity 1 constant 1.0 2 constant 2.0 + fix 1 all rheo/viscosity 1 constant 1.0 2 power 0.1 5e-4 0.001 0.5 Description """"""""""" @@ -47,18 +53,38 @@ means all types from 1 to n (inclusive). A trailing asterisk means all types from m to :math:`N` (inclusive). A middle asterisk means all types from m to n (inclusive). -The *types* definition is followed by the viscosity style, *vstyle*. Currently, -the only option is *constant*. Style *constant* simply applies a constant value -of the viscosity *eta* to each particle of the assigned type. +The *types* definition is followed by the viscosity style, *vstyle*. Two +options are available, *constant* and *power*. Style *constant* simply +applies a constant value of the viscosity *eta* to each particle of the +assigned type. Style *power* is a Hershchel-Bulkley constitutive equation +for the stress :math:`\tau` + +.. math:: + + \tau = \left(\frac{\tau_0}{\dot{\gamma}} + K \dot{\gamma}^{n - 1}\right) \dot{\gamma}, \tau \ge \tau_0 + +where :math:`\dot{\gamma}` is the strain rate and :math:`tau_0` is the critical +yield stress, below which :math:`\dot{\gamma} = 0.0`. To avoid divergences, this +expression is regularized by defining a critical strain rate *gd0*. If the local +strain rate on a particle falls below this limit, a constant viscosity of *eta* +is assigned. This implies a value of + +.. math:: + \tau_0 = \eta * \dot{\gamma}_0 - K \dot{\gamma}_0^N + +as further discussed in :ref:`(Palermo) `. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options +No information about this fix is written to :doc:`binary restart files `. +None of the :doc:`fix_modify ` options are relevant to this fix. No global or per-atom quantities are stored by this fix for access by various :doc:`output commands `. No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. Restrictions """""""""""" @@ -69,7 +95,8 @@ conjuction with :doc:`fix rheo `. The fix group must be set to all. Only one instance of fix rheo/viscosity can be defined. This fix is part of the RHEO package. It is only enabled if -LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. Related commands """""""""""""""" @@ -82,3 +109,9 @@ Default """"""" none + +---------- + +.. _howto_rheo_palermo: + +**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation. diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index 8322ad4ad4..c7270897a4 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -35,7 +35,9 @@ using namespace LAMMPS_NS; ComputeRHEORhoSum::ComputeRHEORhoSum(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute RHEO/rho command"); + if (narg != 4) error->all(FLERR,"Illegal compute RHEO/rho command"); + + self_mass_flag = utils::bnumeric(FLERR, arg[3], false, lmp); comm_forward = 1; comm_reverse = 1; @@ -117,9 +119,15 @@ void ComputeRHEORhoSum::compute_peratom() rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { w = compute_kernel->calc_w(i, j, delx, dely, delz, sqrt(rsq)); - rho[i] += w * mass[type[i]]; - if (newton || j < nlocal) { - rho[j] += w * mass[type[j]]; + + if (self_mass_flag) { + rho[i] += w * mass[type[i]]; + if (newton || j < nlocal) + rho[j] += w * mass[type[j]]; + } else { + rho[i] += w * mass[type[j]]; + if (newton || j < nlocal) + rho[j] += w * mass[type[i]]; } } } diff --git a/src/RHEO/compute_rheo_rho_sum.h b/src/RHEO/compute_rheo_rho_sum.h index 6ec2547b95..491e61ea81 100644 --- a/src/RHEO/compute_rheo_rho_sum.h +++ b/src/RHEO/compute_rheo_rho_sum.h @@ -39,6 +39,7 @@ class ComputeRHEORhoSum : public Compute { class FixRHEO *fix_rheo; private: + int self_mass_flag; double cut, cutsq; class NeighList *list; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 98382b07b5..a355aee5d8 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -60,6 +60,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : interface_flag = 0; surface_flag = 0; oxidation_flag = 0; + self_mass_flag = 0; int i; int n = atom->ntypes; @@ -124,6 +125,8 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : interface_flag = 1; } else if (strcmp(arg[iarg], "rho/sum") == 0) { rhosum_flag = 1; + } else if (strcmp(arg[iarg], "self/mass")) { + self_mass_flag = 1; } else if (strcmp(arg[iarg], "density") == 0) { if (iarg + n >= narg) error->all(FLERR, "Illegal rho0 option in fix rheo"); for (i = 1; i <= n; i++) @@ -142,6 +145,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : iarg += 1; } + if (self_mass_flag && !rhosum_flag) + error->all(FLERR, "Cannot use self/mass setting without rho/sum"); + if (lmp->citeme) lmp->citeme->add(cite_rheo); } @@ -178,7 +184,7 @@ void FixRHEO::post_constructor() if (rhosum_flag) { compute_rhosum = dynamic_cast(modify->add_compute( - "rheo_rhosum all RHEO/RHO/SUM")); + fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag))); compute_rhosum->fix_rheo = this; } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 51e0962000..89a0591824 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -41,6 +41,7 @@ class FixRHEO : public Fix { // Model parameters double h, cut; double *rho0, *csq; + int self_mass_flag; int zmin_kernel, zmin_surface, zmin_splash; int kernel_style, surface_style; double divr_surface; diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index 8539f04277..9f19d4f4da 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -141,6 +141,10 @@ void FixRHEOOxidation::setup_pre_force(int /*vflag*/) void FixRHEOOxidation::pre_force(int /*vflag*/) { + int *status = atom->status; + for (int i = 0; i < atom->nlocal; i++) + if (num_bond[i] != 0) + status[i] |= STATUS_NO_SHIFT; } /* ---------------------------------------------------------------------- */ @@ -215,7 +219,7 @@ void FixRHEOOxidation::post_integrate() if (bflag) continue; // Add bonds to owned atoms - // If newton bond, add to both, otherwise add to whichever has a smaller tag + // If newton bond off, add to both, otherwise add to whichever has a smaller tag if (!newton_bond || tagi < tagj) { if (num_bond[i] == atom->bond_per_atom) diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 5d4134a461..24a6a73021 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -674,7 +674,7 @@ void FixRHEOThermal::create_bonds() if (rsq > cutsq_bond) continue; // Add bonds to owned atoms - // If newton bond, add to both, otherwise add to whichever has a smaller tag + // If newton bond off, add to both, otherwise add to whichever has a smaller tag if (i < nlocal && (!newton_bond || tag[i] < tag[j])) { if (num_bond[i] == atom->bond_per_atom) error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal"); diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index abaa55ca70..2d5d58d494 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -56,6 +56,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : memory->create(gd0, n + 1, "rheo:gd0"); memory->create(K, n + 1, "rheo:K"); memory->create(npow, n + 1, "rheo:npow"); + memory->create(tau0, n + 1, "rheo:tau0"); for (i = 1; i <= n; i++) viscosity_style[i] = NONE; int iarg = 3; @@ -76,7 +77,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : viscosity_style[i] = CONSTANT; eta[i] = eta_one; } - } else if (strcmp(arg[iarg], "power") == 0) { + } else if (strcmp(arg[iarg + 1], "power") == 0) { if (iarg + 5 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/viscosity power", error); comm_forward = 1;