diff --git a/doc/src/fix_wall_gran.txt b/doc/src/fix_wall_gran.txt index 198b22aa3f..b517d48cca 100644 --- a/doc/src/fix_wall_gran.txt +++ b/doc/src/fix_wall_gran.txt @@ -59,31 +59,32 @@ close enough to touch it. The nature of the wall/particle interactions are determined by the {fstyle} setting. It can be any of the styles defined by the -"pair_style gran/*"_pair_gran.html or the more general "pair_style granular"_pair_granular.html" -commands. Currently the options are {hooke}, {hooke/history}, or {hertz/history} for the former, -and {granular} with all the possible options of the associated {pair_coeff} command -for the latter. The equation for the -force between the wall and particles touching it is the same as the -corresponding equation on the "pair_style gran/*"_pair_gran.html -and "pair_style_granular"_pair_granular.html doc -pages, in the limit of one of the two particles going to infinite -radius and mass (flat wall). Specifically, delta = radius - r = -overlap of particle with wall, m_eff = mass of particle, and the -effective radius of contact = RiRj/Ri+Rj is set to the radius -of the particle. +"pair_style gran/*"_pair_gran.html or the more general "pair_style +granular"_pair_granular.html" commands. Currently the options are +{hooke}, {hooke/history}, or {hertz/history} for the former, and +{granular} with all the possible options of the associated +{pair_coeff} command for the latter. The equation for the force +between the wall and particles touching it is the same as the +corresponding equation on the "pair_style gran/*"_pair_gran.html and +"pair_style_granular"_pair_granular.html doc pages, in the limit of +one of the two particles going to infinite radius and mass (flat +wall). Specifically, delta = radius - r = overlap of particle with +wall, m_eff = mass of particle, and the effective radius of contact = +RiRj/Ri+Rj is set to the radius of the particle. The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag} have the same meaning and units as those specified with the -"pair_style gran/*"_pair_gran.html commands. This means a NULL can -be used for either {Kt} or {gamma_t} as described on that page. If a +"pair_style gran/*"_pair_gran.html commands. This means a NULL can be +used for either {Kt} or {gamma_t} as described on that page. If a NULL is used for {Kt}, then a default value is used where {Kt} = 2/7 {Kn}. If a NULL is used for {gamma_t}, then a default value is used where {gamma_t} = 1/2 {gamma_n}. -All the model choices for cohesion, tangential friction, rolling friction -and twisting friction supported by the "pair_style granular"_pair_granular.html -through its {pair_coeff} command are also supported for walls. These are discussed -in greater detail on the doc page for "pair_style granular"_pair_granular.html. +All the model choices for cohesion, tangential friction, rolling +friction and twisting friction supported by the "pair_style +granular"_pair_granular.html through its {pair_coeff} command are also +supported for walls. These are discussed in greater detail on the doc +page for "pair_style granular"_pair_granular.html. Note that you can choose a different force styles and/or different values for the wall/particle coefficients than for particle/particle @@ -121,14 +122,14 @@ Optionally, the wall can be moving, if the {wiggle} or {shear} keywords are appended. Both keywords cannot be used together. For the {wiggle} keyword, the wall oscillates sinusoidally, similar to -the oscillations of particles which can be specified by the -"fix move"_fix_move.html command. This is useful in packing -simulations of granular particles. The arguments to the {wiggle} -keyword specify a dimension for the motion, as well as it's -{amplitude} and {period}. Note that if the dimension is in the plane -of the wall, this is effectively a shearing motion. If the dimension -is perpendicular to the wall, it is more of a shaking motion. A -{zcylinder} wall can only be wiggled in the z dimension. +the oscillations of particles which can be specified by the "fix +move"_fix_move.html command. This is useful in packing simulations of +granular particles. The arguments to the {wiggle} keyword specify a +dimension for the motion, as well as it's {amplitude} and {period}. +Note that if the dimension is in the plane of the wall, this is +effectively a shearing motion. If the dimension is perpendicular to +the wall, it is more of a shaking motion. A {zcylinder} wall can only +be wiggled in the z dimension. Each timestep, the position of a wiggled wall in the appropriate {dim} is set according to this equation: diff --git a/doc/src/fix_wall_gran_region.txt b/doc/src/fix_wall_gran_region.txt index d54b3dc009..8a39d6b642 100644 --- a/doc/src/fix_wall_gran_region.txt +++ b/doc/src/fix_wall_gran_region.txt @@ -48,8 +48,8 @@ Here are snapshots of example models using this command. Corresponding input scripts can be found in examples/granregion. Click on the images to see a bigger picture. Movies of these simulations are "here on the Movies -page"_http://lammps.sandia.gov/movies.html#granregion of the -LAMMPS web site. +page"_http://lammps.sandia.gov/movies.html#granregion of the LAMMPS +web site. :image(JPG/gran_funnel_small.jpg,JPG/gran_funnel.png) :image(JPG/gran_mixer_small.jpg,JPG/gran_mixer.png) @@ -129,15 +129,16 @@ to make the two faces differ by epsilon in their position. The nature of the wall/particle interactions are determined by the {fstyle} setting. It can be any of the styles defined by the -"pair_style gran/*"_pair_gran.html or the more general "pair_style granular"_pair_granular.html" -commands. Currently the options are {hooke}, {hooke/history}, or {hertz/history} for the former, -and {granular} with all the possible options of the associated {pair_coeff} command -for the latter. The equation for the -force between the wall and particles touching it is the same as the -corresponding equation on the "pair_style gran/*"_pair_gran.html -and "pair_style_granular"_pair_granular.html doc -pages, but the effective radius is calculated using the radius of the -particle and the radius of curvature of the wall at the contact point. +"pair_style gran/*"_pair_gran.html or the more general "pair_style +granular"_pair_granular.html" commands. Currently the options are +{hooke}, {hooke/history}, or {hertz/history} for the former, and +{granular} with all the possible options of the associated +{pair_coeff} command for the latter. The equation for the force +between the wall and particles touching it is the same as the +corresponding equation on the "pair_style gran/*"_pair_gran.html and +"pair_style_granular"_pair_granular.html doc pages, but the effective +radius is calculated using the radius of the particle and the radius +of curvature of the wall at the contact point. Specifically, delta = radius - r = overlap of particle with wall, m_eff = mass of particle, and RiRj/Ri+Rj is the effective radius, with @@ -150,17 +151,17 @@ particle. The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag} have the same meaning and units as those specified with the -"pair_style gran/*"_pair_gran.html commands. This means a NULL can -be used for either {Kt} or {gamma_t} as described on that page. If a +"pair_style gran/*"_pair_gran.html commands. This means a NULL can be +used for either {Kt} or {gamma_t} as described on that page. If a NULL is used for {Kt}, then a default value is used where {Kt} = 2/7 {Kn}. If a NULL is used for {gamma_t}, then a default value is used where {gamma_t} = 1/2 {gamma_n}. - -All the model choices for cohesion, tangential friction, rolling friction -and twisting friction supported by the "pair_style granular"_pair_granular.html -through its {pair_coeff} command are also supported for walls. These are discussed -in greater detail on the doc page for "pair_style granular"_pair_granular.html. +All the model choices for cohesion, tangential friction, rolling +friction and twisting friction supported by the "pair_style +granular"_pair_granular.html through its {pair_coeff} command are also +supported for walls. These are discussed in greater detail on the doc +page for "pair_style granular"_pair_granular.html. Note that you can choose a different force styles and/or different values for the 6 wall/particle coefficients than for particle/particle @@ -169,9 +170,9 @@ material. [Restart, fix_modify, output, run start/stop, minimize info:] -Similar to "fix wall/gran"_fix_wall_gran.html command, this fix -writes the shear friction state of atoms interacting with the wall to -"binary restart files"_restart.html, so that a simulation can continue +Similar to "fix wall/gran"_fix_wall_gran.html command, this fix writes +the shear friction state of atoms interacting with the wall to "binary +restart files"_restart.html, so that a simulation can continue correctly if granular potentials with shear "history" effects are being used. This fix also includes info about a moving region in the restart file. See the "read_restart"_read_restart.html command for @@ -185,14 +186,14 @@ So you must re-define your region and if it is a moving region, define its motion attributes in a way that is consistent with the simulation that wrote the restart file. In particular, if you want to change the region motion attributes (e.g. its velocity), then you should ensure -the position/orientation of the region at the initial restart -timestep is the same as it was on the timestep the restart file was -written. If this is not possible, you may need to ignore info in the -restart file by defining a new fix wall/gran/region command in your -restart script, e.g. with a different fix ID. Or if you want to keep -the shear history info but discard the region motion information, you -can use the same fix ID for fix wall/gran/region, but assign it a -region with a different region ID. +the position/orientation of the region at the initial restart timestep +is the same as it was on the timestep the restart file was written. +If this is not possible, you may need to ignore info in the restart +file by defining a new fix wall/gran/region command in your restart +script, e.g. with a different fix ID. Or if you want to keep the +shear history info but discard the region motion information, you can +use the same fix ID for fix wall/gran/region, but assign it a region +with a different region ID. None of the "fix_modify"_fix_modify.html options are relevant to this fix. No global or per-atom quantities are stored by this fix for diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index 73c2bbdd3b..e4b9bb3250 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -19,8 +19,7 @@ pair_style granular command :h3 pair_style granular cutoff :pre -cutoff = global cutoff (optional). See discussion below. :l -:ule +cutoff = global cutoff (optional). See discussion below. :l [Examples:] @@ -44,29 +43,39 @@ pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sd [Description:] -The {granular} styles support a variety of options for the normal, tangential, rolling and twisting -forces resulting from contact between two granular particles. This expands on the options offered -by the "pair gran/*"_pair_gran.html pair styles. The total computed forces and torques are the -sum of various models selected for the normal, tangential, rolling and twisting modes of motion. +The {granular} styles support a variety of options for the normal, +tangential, rolling and twisting forces resulting from contact between +two granular particles. This expands on the options offered by the +"pair gran/*"_pair_gran.html pair styles. The total computed forces +and torques are the sum of various models selected for the normal, +tangential, rolling and twisting modes of motion. -All model choices and parameters are entered in the "pair_coeff"_pair_coeff.html command, as described below. -Unlike e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not global, but can be set to different values for -different combinations of particle types, as determined by the "pair_coeff"_pair_coeff.html command. -If the contact model choice is the same for two particle types, the mixing for the cross-coefficients can be carried out -automatically. This is shown in the second example, where model choices are the same for type 1 - type 1 as for type 2 - type2 -interactions, but coefficients are different. In this case, the coefficients for type 2 - type interactions can be -determined from mixing rules discussed below. -For additional flexibility, coefficients as well as model forms can vary between particle types, -as shown in the third example: -type 1- type 1 interactions are based on a Hertzian normal contact model and 2-2 interactions are based on a DMT cohesive model (see below). -In that example, 1-1 and 2-2 interactions have different model forms, in which case -mixing of coefficients cannot be determined, so 1-2 interactions must be explicitly defined via the -{pair_coeff 1 2} command, otherwise an error would result. +All model choices and parameters are entered in the +"pair_coeff"_pair_coeff.html command, as described below. Unlike +e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not +global, but can be set to different values for different combinations +of particle types, as determined by the "pair_coeff"_pair_coeff.html +command. If the contact model choice is the same for two particle +types, the mixing for the cross-coefficients can be carried out +automatically. This is shown in the second example, where model +choices are the same for type 1 - type 1 as for type 2 - type2 +interactions, but coefficients are different. In this case, the +coefficients for type 2 - type interactions can be determined from +mixing rules discussed below. For additional flexibility, +coefficients as well as model forms can vary between particle types, +as shown in the third example: type 1- type 1 interactions are based +on a Hertzian normal contact model and 2-2 interactions are based on a +DMT cohesive model (see below). In that example, 1-1 and 2-2 +interactions have different model forms, in which case mixing of +coefficients cannot be determined, so 1-2 interactions must be +explicitly defined via the {pair_coeff 1 2} command, otherwise an +error would result. :line -The first required keyword for the {pair_coeff} command is the normal contact model. Currently supported options -for normal contact models and their required arguments are: +The first required keyword for the {pair_coeff} command is the normal +contact model. Currently supported options for normal contact models +and their required arguments are: {hooke} : \(k_n\), \(\eta_\{n0\}\) (or \(e\)) {hertz} : \(k_n\), \(\eta_\{n0\}\) (or \(e\)) @@ -74,137 +83,165 @@ for normal contact models and their required arguments are: {dmt} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\) {jkr} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\) :ol -Here, \(k_n\) is spring stiffness (with units that depend on model choice, see below); -\(\eta_\{n0\}\) is a damping prefactor (or, in its place a coefficient of restitution -\(e\), depending on the choice of damping mode, see below); E is Young's modulus -in units of {force}/{length}^2, i.e. {pressure}; \(\nu\) is Poisson's ratio -and \(\gamma\) is a surface energy density, in units of {energy}/{length}^2. +Here, \(k_n\) is spring stiffness (with units that depend on model +choice, see below); \(\eta_\{n0\}\) is a damping prefactor (or, in its +place a coefficient of restitution \(e\), depending on the choice of +damping mode, see below); E is Young's modulus in units of +{force}/{length}^2, i.e. {pressure}; \(\nu\) is Poisson's ratio and +\(\gamma\) is a surface energy density, in units of +{energy}/{length}^2. + +For the {hooke} model, the normal, elastic component of force acting +on particle {i} due to contact with particle {j} is given by: -For the {hooke} model, the normal, elastic component of force acting on particle {i} due to -contact with particle {j} is given by: \begin\{equation\} \mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\} \end\{equation\} -Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle overlap, -\(R_i, R_j\) are the particle radii, -\(\mathbf\{r\}_\{ij\} = \mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the -two particle centers (note the i-j ordering so that \(F_\{ne\}\) is positive for repulsion), -and \(\mathbf\{n\} = \frac\{\mathbf\{r\}_\{ij\}\}\{\|\mathbf\{r\}_\{ij\}\|\}\). -Therefore, for {hooke}, the units of the spring constant \(k_n\) are {force}/{distance}, -or equivalently {mass}/{time^2}. +Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle +overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} = +\mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two +particle centers (note the i-j ordering so that \(F_\{ne\}\) is +positive for repulsion), and \(\mathbf\{n\} = +\frac\{\mathbf\{r\}_\{ij\}\}\{\|\mathbf\{r\}_\{ij\}\|\}\). Therefore, +for {hooke}, the units of the spring constant \(k_n\) are +{force}/{distance}, or equivalently {mass}/{time^2}. For the {hertz} model, the normal component of force is given by: + \begin\{equation\} \mathbf\{F\}_\{ne, Hertz\} = k_N R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\} \end\{equation\} -Here, \(R_\{eff\} = \frac\{R_i R_j\}\{R_i + R_j\}\) is the effective radius, denoted for simplicity as {R} from here on. -For {hertz}, the units of the spring constant \(k_n\) are {force}/{length}^2, or equivalently -{pressure}. +Here, \(R_\{eff\} = \frac\{R_i R_j\}\{R_i + R_j\}\) is the effective +radius, denoted for simplicity as {R} from here on. For {hertz}, the +units of the spring constant \(k_n\) are {force}/{length}^2, or +equivalently {pressure}. For the {hertz/material} model, the force is given by: + \begin\{equation\} \mathbf\{F\}_\{ne, Hertz/material\} = \frac\{4\}\{3\} E_\{eff\} R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\} \end\{equation\} -Here, \(E_\{eff\} = E = \left(\frac\{1-\nu_i^2\}\{E_i\} + \frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}\) -is the effective Young's modulus, -with \(\nu_i, \nu_j \) the Poisson ratios of the particles of types {i} and {j}. Note that -if the elastic modulus and the shear modulus of the two particles are the same, the {hertz/material} -model is equivalent to the {hertz} model with \(k_N = 4/3 E_\{eff\}\) +Here, \(E_\{eff\} = E = \left(\frac\{1-\nu_i^2\}\{E_i\} + +\frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}\) is the effective Young's +modulus, with \(\nu_i, \nu_j \) the Poisson ratios of the particles of +types {i} and {j}. Note that if the elastic modulus and the shear +modulus of the two particles are the same, the {hertz/material} model +is equivalent to the {hertz} model with \(k_N = 4/3 E_\{eff\}\) + +The {dmt} model corresponds to the +"(Derjaguin-Muller-Toporov)"_#DMT1975 cohesive model, where the force +is simply Hertz with an additional attractive cohesion term: -The {dmt} model corresponds to the "(Derjaguin-Muller-Toporov)"_#DMT1975 cohesive model, -where the force is simply Hertz with an additional attractive cohesion term: \begin\{equation\} \mathbf\{F\}_\{ne, dmt\} = \left(\frac\{4\}\{3\} E R^\{1/2\}\delta_\{ij\}^\{3/2\} - 4\pi\gamma R\right)\mathbf\{n\} \end\{equation\} -The {jkr} model is the "(Johnson-Kendall-Roberts)"_#JKR1971 model, where the force is computed as: +The {jkr} model is the "(Johnson-Kendall-Roberts)"_#JKR1971 model, +where the force is computed as: + \begin\{equation\} \label\{eq:force_jkr\} \mathbf\{F\}_\{ne, jkr\} = \left(\frac\{4Ea^3\}\{3R\} - 2\pi a^2\sqrt\{\frac\{4\gamma E\}\{\pi a\}\}\right)\mathbf\{n\} \end\{equation\} -Here, {a} is the radius of the contact zone, related to the overlap \(\delta\) according to: +Here, {a} is the radius of the contact zone, related to the overlap +\(\delta\) according to: + \begin\{equation\} \delta = a^2/R - 2\sqrt\{\pi \gamma a/E\} \end\{equation\} -LAMMPS internally inverts the equation above to solve for {a} in terms of \(\delta\), then solves for -the force in the previous equation. Additionally, note that the JKR model allows for a tensile force beyond -contact (i.e. for \(\delta < 0\)), up to a maximum of \(3\pi\gamma R\) (also known as -the 'pull-off' force). -Note that this is a hysteretic effect, where particles that are not contacting initially -will not experience force until they come into contact \(\delta \geq 0\); as they move apart -and (\(\delta < 0\)), they experience a tensile force up to \(3\pi\gamma R\), -at which point they lose contact. +LAMMPS internally inverts the equation above to solve for {a} in terms +of \(\delta\), then solves for the force in the previous +equation. Additionally, note that the JKR model allows for a tensile +force beyond contact (i.e. for \(\delta < 0\)), up to a maximum of +\(3\pi\gamma R\) (also known as the 'pull-off' force). Note that this +is a hysteretic effect, where particles that are not contacting +initially will not experience force until they come into contact +\(\delta \geq 0\); as they move apart and (\(\delta < 0\)), they +experience a tensile force up to \(3\pi\gamma R\), at which point they +lose contact. :line -In addition, the normal force is augmented by a damping term of the following -general form: +In addition, the normal force is augmented by a damping term of the +following general form: \begin\{equation\} \mathbf\{F\}_\{n,damp\} = -\eta_n \mathbf\{v\}_\{n,rel\} \end\{equation\} -Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i) \cdot \mathbf\{n\}\) -is the component of relative velocity along \(\mathbf\{n\}\). +Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i) +\cdot \mathbf\{n\}\) is the component of relative velocity along +\(\mathbf\{n\}\). -The optional {damping} keyword to the {pair_coeff} command followed by a keyword -determines the model form of the damping factor \(\eta_n\), and the interpretation -of the \(\eta_\{n0\}\) or \(e\) coefficients specified as part of the normal contact -model settings. The {damping} keyword and corresponding -model form selection may be appended anywhere in the {pair coeff} command. -Note that the choice of damping model affects both the -normal and tangential damping (and depending on other settings, potentially also the twisting damping). -The options for the damping model currently supported are: +The optional {damping} keyword to the {pair_coeff} command followed by +a keyword determines the model form of the damping factor \(\eta_n\), +and the interpretation of the \(\eta_\{n0\}\) or \(e\) coefficients +specified as part of the normal contact model settings. The {damping} +keyword and corresponding model form selection may be appended +anywhere in the {pair coeff} command. Note that the choice of damping +model affects both the normal and tangential damping (and depending on +other settings, potentially also the twisting damping). The options +for the damping model currently supported are: {velocity} {viscoelastic} {tsuji} :ol -If the {damping} keyword is not specified, the {viscoelastic} model is used by default. +If the {damping} keyword is not specified, the {viscoelastic} model is +used by default. -For {damping velocity}, the normal damping is simply equal to the user-specified damping -coefficient in the {normal} model: +For {damping velocity}, the normal damping is simply equal to the +user-specified damping coefficient in the {normal} model: \begin\{equation\} \eta_n = \eta_\{n0\}\ \end\{equation\} -Here, \(\gamma_n\) is the damping coefficient specified for the normal contact model, in units of {mass}/{time}, +Here, \(\gamma_n\) is the damping coefficient specified for the normal +contact model, in units of {mass}/{time}, + +The {damping viscoelastic} model is based on the viscoelastic +treatment of "(Brilliantov et al)"_#Brill1996, where the normal +damping is given by: -The {damping viscoelastic} model is based on the viscoelastic treatment of "(Brilliantov et al)"_#Brill1996, -where the normal damping is given by: \begin\{equation\} \eta_n = \eta_\{n0\}\ a m_\{eff\} \end\{equation\} -Here, \(m_\{eff\} = m_i m_j/(m_i + m_j)\) is the effective mass, {a} is the contact radius, given by \(a =\sqrt\{R\delta\}\) -for all models except {jkr}, for which it is given implicitly according to \(delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\). -In this case, \eta_\{n0\}\ is in units of 1/({time}*{distance}). +Here, \(m_\{eff\} = m_i m_j/(m_i + m_j)\) is the effective mass, {a} +is the contact radius, given by \(a =\sqrt\{R\delta\}\) for all models +except {jkr}, for which it is given implicitly according to \(delta = +a^2/R - 2\sqrt\{\pi \gamma a/E\}\). In this case, \eta_\{n0\}\ is in +units of 1/({time}*{distance}). -The {tsuji} model is based on the work of "(Tsuji et al)"_#Tsuji1992. Here, the -damping coefficient specified as part of the normal model is interpreted -as a restitution coefficient \(e\). The damping constant \(\eta_n\) is given by: +The {tsuji} model is based on the work of "(Tsuji et +al)"_#Tsuji1992. Here, the damping coefficient specified as part of +the normal model is interpreted as a restitution coefficient +\(e\). The damping constant \(\eta_n\) is given by: \begin\{equation\} \eta_n = \alpha (m_\{eff\}k_n)^\{1/2\} \end\{equation\} -For normal contact models based on material parameters, \(k_n = 4/3Ea\). -The parameter \(\alpha\) is related to the restitution coefficient {e} according to: +For normal contact models based on material parameters, \(k_n = +4/3Ea\). The parameter \(\alpha\) is related to the restitution +coefficient {e} according to: \begin\{equation\} \alpha = 1.2728-4.2783e+11.087e^2-22.348e^3+27.467e^4-18.022e^5+4.8218e^6 \end\{equation\} -The dimensionless coefficient of restitution \(e\) specified as part of the normal contact model -parameters should be between 0 and 1, but no error check is performed on this. +The dimensionless coefficient of restitution \(e\) specified as part +of the normal contact model parameters should be between 0 and 1, but +no error check is performed on this. -The total normal force is computed as the sum of the elastic and damping components: +The total normal force is computed as the sum of the elastic and +damping components: \begin\{equation\} \mathbf\{F\}_n = \mathbf\{F\}_\{ne\} + \mathbf\{F\}_\{n,damp\} @@ -212,24 +249,24 @@ The total normal force is computed as the sum of the elastic and damping compone :line -The {pair_coeff} command also requires specification -of the tangential contact model. The required keyword {tangential} is expected, followed by the model -choice and associated parameters. Currently supported tangential model choices and their -expected parameters are as follows: +The {pair_coeff} command also requires specification of the tangential +contact model. The required keyword {tangential} is expected, followed +by the model choice and associated parameters. Currently supported +tangential model choices and their expected parameters are as follows: {linear_nohistory} : \(x_\{\gamma,t\}\), \(\mu_s\) {linear_history} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\) {mindlin} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) {mindlin_rescale} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) :ol -Here, \(x_\{\gamma,t\}\) is a dimensionless multiplier for the normal damping \(\eta_n\) -that determines the magnitude of the -tangential damping, \(\mu_t\) is the tangential (or sliding) friction +Here, \(x_\{\gamma,t\}\) is a dimensionless multiplier for the normal +damping \(\eta_n\) that determines the magnitude of the tangential +damping, \(\mu_t\) is the tangential (or sliding) friction coefficient, and \(k_t\) is the tangential stiffness coefficient. -For {tangential linear_nohistory}, a simple velocity-dependent Coulomb friction criterion is used, -which mimics the behavior -of the {pair gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by: +For {tangential linear_nohistory}, a simple velocity-dependent Coulomb +friction criterion is used, which mimics the behavior of the {pair +gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by: \begin\{equation\} \mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|\mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\} @@ -241,41 +278,52 @@ The tangential damping force \(\mathbf\{F\}_\mathrm\{t,damp\}\) is given by: \mathbf\{F\}_\mathrm\{t,damp\} = -\eta_t \mathbf\{v\}_\{t,rel\} \end\{equation\} -The tangential damping prefactor \(\eta_t\) is calculated by scaling the normal damping \(\eta_n\) (see above): +The tangential damping prefactor \(\eta_t\) is calculated by scaling +the normal damping \(\eta_n\) (see above): + \begin\{equation\} \eta_t = -x_\{\gamma,t\} \eta_n \end\{equation\} -The normal damping prefactor \(\eta_n\) is determined by the choice of the {damping} keyword, as discussed above. -Thus, the {damping} keyword also affects the tangential damping. -The parameter \(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the literature use -\(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009, "Tsuji et al"_#Tsuji1992, "Silbert et al"_#Silbert2001). -The relative tangential velocity at the point of contact is given by -\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf\{n\}\), -where \(\mathbf\{v\}_\{t\} = \mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), -\(\mathbf\{v\}_r = \mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force is -\(\mathbf\{t\} = \mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\). +The normal damping prefactor \(\eta_n\) is determined by the choice of +the {damping} keyword, as discussed above. Thus, the {damping} +keyword also affects the tangential damping. The parameter +\(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the +literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009, +"Tsuji et al"_#Tsuji1992, "Silbert et al"_#Silbert2001). The relative +tangential velocity at the point of contact is given by +\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i + +R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} = +\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r = +\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force +is \(\mathbf\{t\} = +\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\). The normal force value \(F_\{n0\}\) used to compute the critical force depends on the form of the contact model. For non-cohesive models -({hertz}, {hertz/material}, {hooke}), it is given by the magnitude of the normal force: +({hertz}, {hertz/material}, {hooke}), it is given by the magnitude of +the normal force: \begin\{equation\} F_\{n0\} = \|\mathbf\{F\}_n\| \end\{equation\} -For cohesive models such as {jkr} and {dmt}, the critical force is adjusted so that the critical tangential -force approaches \(\mu_t F_\{pulloff\}\), see "Marshall"_#Marshall2009, equation 43, and "Thornton"_#Thornton1991. -For both models, \(F_\{n0\}\) takes the form: +For cohesive models such as {jkr} and {dmt}, the critical force is +adjusted so that the critical tangential force approaches \(\mu_t +F_\{pulloff\}\), see "Marshall"_#Marshall2009, equation 43, and +"Thornton"_#Thornton1991. For both models, \(F_\{n0\}\) takes the +form: \begin\{equation\} F_\{n0\} = \|\mathbf\{F\}_ne + 2 F_\{pulloff\}\| \end\{equation\} -Where \(F_\{pulloff\} = 3\pi \gamma R \) for {jkr}, and \(F_\{pulloff\} = 4\pi \gamma R \) for {dmt}. +Where \(F_\{pulloff\} = 3\pi \gamma R \) for {jkr}, and +\(F_\{pulloff\} = 4\pi \gamma R \) for {dmt}. -The remaining tangential options all use accumulated tangential displacement (i.e. contact history). This -is discussed below in the context of the {linear_history} option, but the same treatment of the +The remaining tangential options all use accumulated tangential +displacement (i.e. contact history). This is discussed below in the +context of the {linear_history} option, but the same treatment of the accumulated displacement applies to the other options as well. For {tangential linear_history}, the tangential force is given by: @@ -284,49 +332,55 @@ For {tangential linear_history}, the tangential force is given by: \mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t\mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\} \end\{equation\} -Here, \(\mathbf\{\xi\}\) is the tangential displacement accumulated during the entire -duration of the contact: +Here, \(\mathbf\{\xi\}\) is the tangential displacement accumulated +during the entire duration of the contact: \begin\{equation\} \mathbf\{\xi\} = \int_\{t0\}^t \mathbf\{v\}_\{t,rel\}(\tau) \mathrm\{d\}\tau \end\{equation\} -This accumulated tangential displacement must be adjusted to account for changes -in the frame of reference -of the contacting pair of particles during contact. This occurs due to the overall motion of the contacting particles -in a rigid-body-like fashion during the duration of the contact. There are two modes of motion -that are relevant: the 'tumbling' rotation of the contacting pair, which changes the orientation of the -plane in which tangential displacement occurs; and 'spinning' rotation of the contacting pair -about the vector connecting their centers of mass (\(\mathbf\{n\}\)). -Corrections due to the former mode of motion are -made by rotating the accumulated displacement into the plane that is tangential -to the contact vector at each step, -or equivalently removing any component of the tangential displacement -that lies along \(\mathbf\{n\}\), and rescaling to preserve the magnitude. -This follows the discussion in "Luding"_#Luding2008, see equation 17 and -relevant discussion in that work: +This accumulated tangential displacement must be adjusted to account +for changes in the frame of reference of the contacting pair of +particles during contact. This occurs due to the overall motion of the +contacting particles in a rigid-body-like fashion during the duration +of the contact. There are two modes of motion that are relevant: the +'tumbling' rotation of the contacting pair, which changes the +orientation of the plane in which tangential displacement occurs; and +'spinning' rotation of the contacting pair about the vector connecting +their centers of mass (\(\mathbf\{n\}\)). Corrections due to the +former mode of motion are made by rotating the accumulated +displacement into the plane that is tangential to the contact vector +at each step, or equivalently removing any component of the tangential +displacement that lies along \(\mathbf\{n\}\), and rescaling to +preserve the magnitude. This follows the discussion in +"Luding"_#Luding2008, see equation 17 and relevant discussion in that +work: \begin\{equation\} \mathbf\{\xi\} = \left(\mathbf\{\xi'\} - (\mathbf\{n\} \cdot \mathbf\{\xi'\})\mathbf\{n\}\right) \frac\{\|\mathbf\{\xi'\}\|\}\{\|\mathbf\{\xi'\}\| - \mathbf\{n\}\cdot\mathbf\{\xi'\}\} \label\{eq:rotate_displacements\} \end\{equation\} -Here, \(\mathbf\{\xi'\}\) is the accumulated displacement prior to the current time step and -\(\mathbf\{\xi\}\) is the corrected displacement. Corrections to the displacement -due to the second mode of motion described above (rotations about \(\mathbf\{n\}\)) -are not currently implemented, but are expected to be minor for most simulations. +Here, \(\mathbf\{\xi'\}\) is the accumulated displacement prior to the +current time step and \(\mathbf\{\xi\}\) is the corrected +displacement. Corrections to the displacement due to the second mode +of motion described above (rotations about \(\mathbf\{n\}\)) are not +currently implemented, but are expected to be minor for most +simulations. -Furthermore, when the tangential force exceeds the critical force, -the tangential displacement is re-scaled to match the value for the critical force (see "Luding"_#Luding2008, -equation 20 and related discussion): +Furthermore, when the tangential force exceeds the critical force, the +tangential displacement is re-scaled to match the value for the +critical force (see "Luding"_#Luding2008, equation 20 and related +discussion): \begin\{equation\} \mathbf\{\xi\} = -\frac\{1\}\{k_t\}\left(\mu_t F_\{n0\}\mathbf\{t\} + \mathbf\{F\}_\{t,damp\}\right) \end\{equation\} -The tangential force is added to the total normal force (elastic plus damping) to produce the total force -on the particle. The tangential force also acts at the contact point (defined as the center of the overlap region) -to induce a torque on each particle according to: +The tangential force is added to the total normal force (elastic plus +damping) to produce the total force on the particle. The tangential +force also acts at the contact point (defined as the center of the +overlap region) to induce a torque on each particle according to: \begin\{equation\} \mathbf\{\tau\}_i = -(R_i - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t @@ -343,72 +397,90 @@ option by an additional factor of {a}, the radius of the contact region. The tan \mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\} \end\{equation\} -Here, {a} is the radius of the contact region, given by \(a = \delta R\) for all normal contact models, -except for {jkr}, where it is given implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), -see discussion above. To match the Mindlin solution, one should set \(k_t = 8G\), where -\(G\) is the shear modulus, related to Young's modulus \(E\) by \(G = E/(2(1+\nu))\), where \(\nu\) -is Poisson's ratio. This can also be achieved by specifying {NULL} for \(k_t\), in which case -a normal contact model that specifies material parameters \(E\) and \(\nu\) is required (e.g. {hertz/material}, -{dmt} or {jkr}). In this case, mixing of the shear modulus for different particle types {i} and {j} is done according -to: +Here, {a} is the radius of the contact region, given by \(a = \delta +R\) for all normal contact models, except for {jkr}, where it is given +implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see +discussion above. To match the Mindlin solution, one should set \(k_t += 8G\), where \(G\) is the shear modulus, related to Young's modulus +\(E\) by \(G = E/(2(1+\nu))\), where \(\nu\) is Poisson's ratio. This +can also be achieved by specifying {NULL} for \(k_t\), in which case a +normal contact model that specifies material parameters \(E\) and +\(\nu\) is required (e.g. {hertz/material}, {dmt} or {jkr}). In this +case, mixing of the shear modulus for different particle types {i} and +{j} is done according to: + \begin\{equation\} 1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j \end\{equation\} -The {mindlin_rescale} option uses the same form as {mindlin}, but the magnitude of the tangential -displacement is re-scaled as the contact unloads, i.e. if \(a < a_\{t_\{n-1\}\}\): +The {mindlin_rescale} option uses the same form as {mindlin}, but the +magnitude of the tangential displacement is re-scaled as the contact +unloads, i.e. if \(a < a_\{t_\{n-1\}\}\): + \begin\{equation\} \mathbf\{\xi\} = \mathbf\{\xi_\{t_\{n-1\}\}\} \frac\{a\}\{a_\{t_\{n-1\}\}\} \end\{equation\} -Here, \(t_\{n-1\}\) indicates the value at the previous time step. This rescaling -accounts for the fact that a decrease in the contact area upon unloading leads to the contact -being unable to support the previous tangential loading, and spurious energy is created -without the rescaling above ("Walton"_#WaltonPC ). See also discussion in "Thornton et al, 2013"_#Thornton2013 -, particularly equation 18(b) of that work and associated discussion. +Here, \(t_\{n-1\}\) indicates the value at the previous time +step. This rescaling accounts for the fact that a decrease in the +contact area upon unloading leads to the contact being unable to +support the previous tangential loading, and spurious energy is +created without the rescaling above ("Walton"_#WaltonPC ). See also +discussion in "Thornton et al, 2013"_#Thornton2013 , particularly +equation 18(b) of that work and associated discussion. :line -The optional {rolling} keyword enables rolling friction, which resists pure rolling -motion of particles. The options currently supported are: +The optional {rolling} keyword enables rolling friction, which resists +pure rolling motion of particles. The options currently supported are: {none} {sds} : \(k_\{roll\}\), \(\gamma_\{roll\}\), \(\mu_\{roll\}\) :ol If the {rolling} keyword is not specified, the model defaults to {none}. -For {rolling sds}, rolling friction is computed via a spring-dashpot-slider, using a -'pseudo-force' formulation, as detailed by "Luding"_#Luding2008. Unlike the formulation -in "Marshall"_#Marshall2009, this allows for the required adjustment of -rolling displacement due to changes in the frame of reference of the contacting pair. -The rolling pseudo-force is computed analogously to the tangential force: +For {rolling sds}, rolling friction is computed via a +spring-dashpot-slider, using a 'pseudo-force' formulation, as detailed +by "Luding"_#Luding2008. Unlike the formulation in +"Marshall"_#Marshall2009, this allows for the required adjustment of +rolling displacement due to changes in the frame of reference of the +contacting pair. The rolling pseudo-force is computed analogously to +the tangential force: \begin\{equation\} \mathbf\{F\}_\{roll,0\} = k_\{roll\} \mathbf\{\xi\}_\{roll\} - \gamma_\{roll\} \mathbf\{v\}_\{roll\} \end\{equation\} -Here, \(\mathbf\{v\}_\{roll\} = -R(\mathbf\{\Omega\}_i - \mathbf\{\Omega\}_j) \times \mathbf\{n\}\) is the -relative rolling velocity, as given in "Wang et al"_#Wang2015 and "Luding"_#Luding2008. This differs -from the expressions given by "Kuhn and Bagi"_#Kuhn2004 and used in "Marshall"_#Marshall2009; -see "Wang et al"_#Wang2015 for details. The rolling displacement is given by: +Here, \(\mathbf\{v\}_\{roll\} = -R(\mathbf\{\Omega\}_i - +\mathbf\{\Omega\}_j) \times \mathbf\{n\}\) is the relative rolling +velocity, as given in "Wang et al"_#Wang2015 and +"Luding"_#Luding2008. This differs from the expressions given by "Kuhn +and Bagi"_#Kuhn2004 and used in "Marshall"_#Marshall2009; see "Wang et +al"_#Wang2015 for details. The rolling displacement is given by: \begin\{equation\} \mathbf\{\xi\}_\{roll\} = \int_\{t_0\}^t \mathbf\{v\}_\{roll\} (\tau) \mathrm\{d\} \tau \end\{equation\} -A Coulomb friction criterion truncates the rolling pseudo-force if it exceeds a critical value: +A Coulomb friction criterion truncates the rolling pseudo-force if it +exceeds a critical value: + \begin\{equation\} \mathbf\{F\}_\{roll\} = min(\mu_\{roll\} F_\{n,0\}, \|\mathbf\{F\}_\{roll,0\}\|)\mathbf\{k\} \end\{equation\} -Here, \(\mathbf\{k\} = \mathbf\{v\}_\{roll\}/\|\mathbf\{v\}_\{roll\}\|\) is the direction of the pseudo-force. -As with tangential displacement, the rolling displacement is rescaled when the critical -force is exceeded, so that the spring length corresponds the critical force. Additionally, the -displacement is adjusted to account for rotations of the frame of reference of the two -contacting particles in a manner analogous to the tangential displacement. +Here, \(\mathbf\{k\} = +\mathbf\{v\}_\{roll\}/\|\mathbf\{v\}_\{roll\}\|\) is the direction of +the pseudo-force. As with tangential displacement, the rolling +displacement is rescaled when the critical force is exceeded, so that +the spring length corresponds the critical force. Additionally, the +displacement is adjusted to account for rotations of the frame of +reference of the two contacting particles in a manner analogous to the +tangential displacement. -The rolling pseudo-force does not contribute to the total force on either particle (hence 'pseudo'), -but acts only to induce an equal and opposite torque on each particle, according to: +The rolling pseudo-force does not contribute to the total force on +either particle (hence 'pseudo'), but acts only to induce an equal and +opposite torque on each particle, according to: \begin\{equation\} \tau_\{roll,i\} = R_\{eff\} \mathbf\{n\} \times \mathbf\{F\}_\{roll\} @@ -420,9 +492,10 @@ but acts only to induce an equal and opposite torque on each particle, according :line -The optional {twisting} keyword enables twisting friction, which resists -rotation of two contacting particles about the vector \(\mathbf\{n\}\) that connects their -centers. The options currently supported are: +The optional {twisting} keyword enables twisting friction, which +resists rotation of two contacting particles about the vector +\(\mathbf\{n\}\) that connects their centers. The options currently +supported are: {none} {sds} : \(k_\{twist\}\), \(\gamma_\{twist\}\), \(\mu_\{twist\}\) @@ -430,36 +503,42 @@ centers. The options currently supported are: If the {twisting} keyword is not specified, the model defaults to {none}. -For both {twisting sds} and {twisting marshall}, a history-dependent spring-dashpot-slider is used to compute the twisting -torque. Because twisting displacement is a scalar, there is no need to adjust for changes -in the frame of reference due to rotations of the particle pair. The formulation in -"Marshall"_#Marshall2009 therefore provides the most straightforward treatment: +For both {twisting sds} and {twisting marshall}, a history-dependent +spring-dashpot-slider is used to compute the twisting torque. Because +twisting displacement is a scalar, there is no need to adjust for +changes in the frame of reference due to rotations of the particle +pair. The formulation in "Marshall"_#Marshall2009 therefore provides +the most straightforward treatment: \begin\{equation\} \tau_\{twist,0\} = -k_\{twist\}\xi_\{twist\} - \gamma_\{twist\}\Omega_\{twist\} \end\{equation\} -Here \(\xi_\{twist\} = \int_\{t_0\}^t \Omega_\{twist\} (\tau) \mathrm\{d\}\tau\) is the twisting -angular displacement, and \(\Omega_\{twist\} = (\mathbf\{\Omega\}_i - \mathbf\{\Omega\}_j) \cdot \mathbf\{n\}\) -is the relative twisting angular velocity. The torque is then truncated according to: +Here \(\xi_\{twist\} = \int_\{t_0\}^t \Omega_\{twist\} (\tau) +\mathrm\{d\}\tau\) is the twisting angular displacement, and +\(\Omega_\{twist\} = (\mathbf\{\Omega\}_i - \mathbf\{\Omega\}_j) \cdot +\mathbf\{n\}\) is the relative twisting angular velocity. The torque +is then truncated according to: \begin\{equation\} \tau_\{twist\} = min(\mu_\{twist\} F_\{n,0\}, \tau_\{twist,0\}) \end\{equation\} -Similar to the sliding and rolling displacement, the angular displacement is -rescaled so that it corresponds to the critical value if the twisting torque -exceeds this critical value: +Similar to the sliding and rolling displacement, the angular +displacement is rescaled so that it corresponds to the critical value +if the twisting torque exceeds this critical value: \begin\{equation\} \xi_\{twist\} = \frac\{1\}\{k_\{twist\}\} (\mu_\{twist\} F_\{n,0\}sgn(\Omega_\{twist\}) - \gamma_\{twist\}\Omega_\{twist\}) \end\{equation\} -For {twisting sds}, the coefficients \(k_\{twist\}, \gamma_\{twist\}\) and \(\mu_\{twist\}\) are -simply the user input parameters that follow the {twisting sds} keywords in the {pair_coeff} command. +For {twisting sds}, the coefficients \(k_\{twist\}, \gamma_\{twist\}\) +and \(\mu_\{twist\}\) are simply the user input parameters that follow +the {twisting sds} keywords in the {pair_coeff} command. -For {twisting_marshall}, the coefficients are expressed in terms of sliding friction coefficients, -as discussed in "Marshall"_#Marshall2009 (see equations 32 and 33 of that work): +For {twisting_marshall}, the coefficients are expressed in terms of +sliding friction coefficients, as discussed in +"Marshall"_#Marshall2009 (see equations 32 and 33 of that work): \begin\{equation\} k_\{twist\} = 0.5k_ta^2 @@ -485,19 +564,25 @@ Finally, the twisting torque on each particle is given by: :line -LAMMPS automatically sets pairwise cutoff values for {pair_style granular} based on particle radii (and in the case -of {jkr} pull-off distances). In the vast majority of situations, this is adequate. -However, a cutoff value can optionally be appended to the {pair_style granular} command to specify -a global cutoff (i.e. a cutoff for all atom types). Additionally, the optional {cutoff} keyword -can be passed to the {pair_coeff} command, followed by a cutoff value. -This will set a pairwise cutoff for the atom types in the {pair_coeff} command. -These options may be useful in some rare cases where the automatic cutoff determination is not sufficient, e.g. -if particle diameters are being modified via the {fix adapt} command. In that case, the global cutoff -specified as part of the {pair_style granular} command is applied to all atom types, unless it is -overridden for a given atom type combination by the {cutoff} value specified in the {pair coeff} command. -If {cutoff} is only specified in the {pair coeff} command and no global -cutoff is appended to the {pair_style granular} command, then LAMMPS will use that cutoff for the specified -atom type combination, and automatically set pairwise cutoffs for the remaining atom types. +LAMMPS automatically sets pairwise cutoff values for {pair_style +granular} based on particle radii (and in the case of {jkr} pull-off +distances). In the vast majority of situations, this is adequate. +However, a cutoff value can optionally be appended to the {pair_style +granular} command to specify a global cutoff (i.e. a cutoff for all +atom types). Additionally, the optional {cutoff} keyword can be passed +to the {pair_coeff} command, followed by a cutoff value. This will +set a pairwise cutoff for the atom types in the {pair_coeff} command. +These options may be useful in some rare cases where the automatic +cutoff determination is not sufficient, e.g. if particle diameters +are being modified via the {fix adapt} command. In that case, the +global cutoff specified as part of the {pair_style granular} command +is applied to all atom types, unless it is overridden for a given atom +type combination by the {cutoff} value specified in the {pair coeff} +command. If {cutoff} is only specified in the {pair coeff} command +and no global cutoff is appended to the {pair_style granular} command, +then LAMMPS will use that cutoff for the specified atom type +combination, and automatically set pairwise cutoffs for the remaining +atom types. :line @@ -529,20 +614,21 @@ The "pair_modify"_pair_modify.html mix, shift, table, and tail options are not relevant for granular pair styles. Mixing of coefficients is carried out using geometric averaging for -most quantities, e.g. if friction coefficient for type 1-type 1 interactions -is set to \(\mu_1\), and friction coefficient for type 2-type 2 interactions -is set to \(\mu_2\), the friction coefficient for type1-type2 interactions -is computed as \(\sqrt\{\mu_1\mu_2\}\) (unless explicitly specified to -a different value by a {pair_coeff 1 2 ...} command. The exception to this is -elastic modulus, only applicable to {hertz/material}, {dmt} and {jkr} normal -contact models. In that case, the effective elastic modulus is computed as: +most quantities, e.g. if friction coefficient for type 1-type 1 +interactions is set to \(\mu_1\), and friction coefficient for type +2-type 2 interactions is set to \(\mu_2\), the friction coefficient +for type1-type2 interactions is computed as \(\sqrt\{\mu_1\mu_2\}\) +(unless explicitly specified to a different value by a {pair_coeff 1 2 +...} command. The exception to this is elastic modulus, only +applicable to {hertz/material}, {dmt} and {jkr} normal contact +models. In that case, the effective elastic modulus is computed as: \begin\{equation\} E_\{eff,ij\} = \left(\frac\{1-\nu_i^2\}\{E_i\} + \frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\} \end\{equation\} -If the {i-j} coefficients \(E_\{ij\}\) and \(\nu_\{ij\}\) are explicitly specified, -the effective modulus is computed as: +If the {i-j} coefficients \(E_\{ij\}\) and \(\nu_\{ij\}\) are +explicitly specified, the effective modulus is computed as: \begin\{equation\} E_\{eff,ij\} = \left(\frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\} + \frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\}\right)^\{-1\} @@ -610,57 +696,70 @@ compute depend on atom velocities. See the [Default:] -For the {pair_coeff} settings: {damping viscoelastic}, {rolling none}, {twisting none} +For the {pair_coeff} settings: {damping viscoelastic}, {rolling none}, +{twisting none}. [References:] :link(Brill1996) -[(Brilliantov et al, 1996)] Brilliantov, N. V., Spahn, F., Hertzsch, J. M., & Poschel, T. (1996). -Model for collisions in granular gases. Physical review E, 53(5), 5382. +[(Brilliantov et al, 1996)] Brilliantov, N. V., Spahn, F., Hertzsch, +J. M., & Poschel, T. (1996). Model for collisions in granular +gases. Physical review E, 53(5), 5382. -:link(Tsuji1992) -[(Tsuji et al, 1992)] Tsuji, Y., Tanaka, T., & Ishida, T. (1992). Lagrangian numerical simulation of plug flow of -cohesionless particles in a horizontal pipe. Powder technology, 71(3), 239-250. +:link(Tsuji1992) +[(Tsuji et al, 1992)] Tsuji, Y., Tanaka, T., & Ishida, +T. (1992). Lagrangian numerical simulation of plug flow of +cohesionless particles in a horizontal pipe. Powder technology, 71(3), +239-250. :link(JKR1971) -[(Johnson et al, 1971)] Johnson, K. L., Kendall, K., & Roberts, A. D. (1971). -Surface energy and the contact of elastic solids. Proc. R. Soc. Lond. A, 324(1558), 301-313. +[(Johnson et al, 1971)] Johnson, K. L., Kendall, K., & Roberts, +A. D. (1971). Surface energy and the contact of elastic +solids. Proc. R. Soc. Lond. A, 324(1558), 301-313. :link(DMT1975) -[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. +[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. :link(Luding2008) -[(Luding, 2008)] Luding, S. (2008). Cohesive, frictional powders: contact models for tension. Granular matter, 10(4), 235. +[(Luding, 2008)] Luding, S. (2008). Cohesive, frictional powders: +contact models for tension. Granular matter, 10(4), 235. :link(Marshall2009) -[(Marshall, 2009)] Marshall, J. S. (2009). Discrete-element modeling of particulate aerosol flows. -Journal of Computational Physics, 228(5), 1541-1561. +[(Marshall, 2009)] Marshall, J. S. (2009). Discrete-element modeling +of particulate aerosol flows. Journal of Computational Physics, +228(5), 1541-1561. :link(Silbert2001) -[(Silbert, 2001)] Silbert, L. E., Ertas, D., Grest, G. S., Halsey, T. C., Levine, D., & Plimpton, S. J. (2001). -Granular flow down an inclined plane: Bagnold scaling and rheology. Physical Review E, 64(5), 051302. +[(Silbert, 2001)] Silbert, L. E., Ertas, D., Grest, G. S., Halsey, +T. C., Levine, D., & Plimpton, S. J. (2001). Granular flow down an +inclined plane: Bagnold scaling and rheology. Physical Review E, +64(5), 051302. :link(Kuhn2004) -[(Kuhn and Bagi, 2005)] Kuhn, M. R., & Bagi, K. (2004). Contact rolling and deformation in granular media. -International journal of solids and structures, 41(21), 5793-5820. +[(Kuhn and Bagi, 2005)] Kuhn, M. R., & Bagi, K. (2004). Contact +rolling and deformation in granular media. International journal of +solids and structures, 41(21), 5793-5820. :link(Wang2015) -[(Wang et al, 2015)] Wang, Y., Alonso-Marroquin, F., & Guo, W. W. (2015). -Rolling and sliding in 3-D discrete element models. Particuology, 23, 49-55. +[(Wang et al, 2015)] Wang, Y., Alonso-Marroquin, F., & Guo, +W. W. (2015). Rolling and sliding in 3-D discrete element +models. Particuology, 23, 49-55. :link(Thornton1991) -[(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the presence of adhesion. -J. Phys. D: Appl. Phys. 24 1942 +[(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the +presence of adhesion. J. Phys. D: Appl. Phys. 24 1942 :link(Mindlin1949) -[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies in contact. -J. Appl. Mech., ASME 16, 259-268. +[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies +in contact. J. Appl. Mech., ASME 16, 259-268. :link(Thornton2013) -[(Thornton et al, 2013)] Thornton, C., Cummins, S. J., & Cleary, P. W. (2013). -An investigation of the comparative behaviour of alternative contact force models -during inelastic collisions. Powder Technology, 233, 30-46. +[(Thornton et al, 2013)] Thornton, C., Cummins, S. J., & Cleary, +P. W. (2013). An investigation of the comparative behaviour of +alternative contact force models during inelastic collisions. Powder +Technology, 233, 30-46. :link(WaltonPC) [(Otis R. Walton)] Walton, O.R., Personal Communication diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index c63a2c0ba8..9925c37e4b 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -124,21 +124,27 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : roll_model = twist_model = NONE; while (iarg < narg) { if (strcmp(arg[iarg], "hooke") == 0) { - if (iarg + 2 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hooke option"); + if (iarg + 2 >= narg) + error->all(FLERR,"Illegal fix wall/gran command, " + "not enough parameters provided for Hooke option"); normal_model = NORMAL_HOOKE; normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping iarg += 3; } else if (strcmp(arg[iarg], "hertz") == 0) { int num_coeffs = 2; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + if (iarg + num_coeffs >= narg) + error->all(FLERR,"Illegal fix wall/gran command, " + "not enough parameters provided for Hertz option"); normal_model = NORMAL_HERTZ; normal_coeffs[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "hertz/material") == 0) { int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + if (iarg + num_coeffs >= narg) + error->all(FLERR,"Illegal fix wall/gran command, " + "not enough parameters provided for Hertz option"); normal_model = HERTZ_MATERIAL; Emod = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping @@ -147,7 +153,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : normal_coeffs[2] = poiss; iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "dmt") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal fix wall/gran command, not enough parameters provided for Hertz option"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal fix wall/gran command, " + "not enough parameters provided for Hertz option"); normal_model = DMT; Emod = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //damping @@ -157,7 +165,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } else if (strcmp(arg[iarg], "jkr") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for JKR option"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal wall/gran command, " + "not enough parameters provided for JKR option"); beyond_contact = 1; normal_model = JKR; Emod = force->numeric(FLERR,arg[iarg+1]); //E @@ -168,7 +178,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : normal_coeffs[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } else if (strcmp(arg[iarg], "damping") == 0) { - if (iarg+1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters provided for damping model"); + if (iarg+1 >= narg) + error->all(FLERR, "Illegal wall/gran command, " + "not enough parameters provided for damping model"); if (strcmp(arg[iarg+1], "velocity") == 0) { damping_model = VELOCITY; iarg += 1; @@ -178,58 +190,80 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg+1], "tsuji") == 0) { damping_model = TSUJI; iarg += 1; - } else error->all(FLERR, "Illegal wall/gran command, unrecognized damping model"); + } else error->all(FLERR, "Illegal wall/gran command, " + "unrecognized damping model"); iarg += 1; } else if (strcmp(arg[iarg], "tangential") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword"); + if (iarg + 1 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "must specify tangential model after tangential keyword"); if (strcmp(arg[iarg+1], "linear_nohistory") == 0) { - if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (iarg + 3 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); tangential_model = TANGENTIAL_NOHISTORY; tangential_coeffs[0] = 0; - tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); //gammat - tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff. + // gammat and friction coeff + tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+2]); + tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) || (strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strcmp(arg[iarg+1], "linear_history") == 0) tangential_model = TANGENTIAL_HISTORY; - else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model = TANGENTIAL_MINDLIN; - else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model = TANGENTIAL_MINDLIN_RESCALE; - if ((tangential_model == TANGENTIAL_MINDLIN || tangential_model == TANGENTIAL_MINDLIN_RESCALE) && + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); + if (strcmp(arg[iarg+1], "linear_history") == 0) + tangential_model = TANGENTIAL_HISTORY; + else if (strcmp(arg[iarg+1], "mindlin") == 0) + tangential_model = TANGENTIAL_MINDLIN; + else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) + tangential_model = TANGENTIAL_MINDLIN_RESCALE; + if ((tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_RESCALE) && (strcmp(arg[iarg+2], "NULL") == 0)) { if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) { - error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties"); + error->all(FLERR, "NULL setting for Mindlin tangential " + "stiffness requires a normal contact model " + "that specifies material properties"); } tangential_coeffs[0] = 4*(2-poiss)*(1+poiss)/Emod; } else { tangential_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt } tangential_history = 1; - tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + // gammat and friction coeff + tangential_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); + tangential_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else { - error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); + error->all(FLERR, "Illegal pair_coeff command, " + "tangential model not recognized"); } } else if (strcmp(arg[iarg], "rolling") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters"); + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal wall/gran command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0) { roll_model = ROLL_NONE; iarg += 2; } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for rolling model"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal wall/gran command, " + "not enough parameters provided for rolling model"); roll_model = ROLL_SDS; roll_history = 1; - roll_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kR - roll_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. + // kR, gammaR, rolling friction coeff + roll_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); + roll_coeffs[1] = force->numeric(FLERR,arg[iarg+3]); + roll_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else { - error->all(FLERR, "Illegal wall/gran command, rolling friction model not recognized"); + error->all(FLERR, "Illegal wall/gran command, " + "rolling friction model not recognized"); } } else if (strcmp(arg[iarg], "twisting") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Illegal wall/gran command, not enough parameters"); + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal wall/gran command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0) { twist_model = TWIST_NONE; iarg += 2; @@ -238,7 +272,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : twist_history = 1; iarg += 2; } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal wall/gran command, not enough parameters provided for twist model"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal wall/gran command, " + "not enough parameters provided for twist model"); twist_model = TWIST_SDS; twist_history = 1; twist_coeffs[0] = force->numeric(FLERR,arg[iarg+2]); //kt @@ -246,7 +282,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : twist_coeffs[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. iarg += 5; } else { - error->all(FLERR, "Illegal wall/gran command, twisting friction model not recognized"); + error->all(FLERR, "Illegal wall/gran command, " + "twisting friction model not recognized"); } } else if (strcmp(arg[iarg], "xplane") == 0 || strcmp(arg[iarg], "yplane") == 0 || @@ -472,8 +509,6 @@ void FixWallGran::init() 27.467*pow(cor,4)-18.022*pow(cor,5)+ 4.8218*pow(cor,6); } - - } /* ---------------------------------------------------------------------- */ @@ -616,7 +651,8 @@ void FixWallGran::post_force(int /*vflag*/) else { if (pairstyle == GRANULAR && normal_model == JKR && use_history) { if ((history_one[i][0] == 0) && (rsq > radius[i]*radius[i])) { - // Particles have not contacted yet, and are outside of contact distance + // Particles have not contacted yet, + // and are outside of contact distance for (j = 0; j < size_history; j++) history_one[i][j] = 0.0; continue; @@ -824,7 +860,8 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, history[1] += vtr2*dt; history[2] += vtr3*dt; } - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); // rotate shear displacements @@ -954,7 +991,8 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, history[1] += vtr2*dt; history[2] += vtr3*dt; } - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); // rotate history displacements @@ -1015,372 +1053,382 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, torque[2] -= radius*tor3; } +/* ---------------------------------------------------------------------- */ void FixWallGran::granular(double rsq, double dx, double dy, double dz, - double *vwall, double rwall, double *v, - double *f, double *omega, double *torque, - double radius, double meff, double *history, - double *contact) + double *vwall, double rwall, double *v, + double *f, double *omega, double *torque, + double radius, double meff, double *history, + double *contact) { - double fx,fy,fz,nx,ny,nz; - double radsum,r,rinv; - double Reff, delta, dR, dR2; + double fx,fy,fz,nx,ny,nz; + double radsum,r,rinv; + double Reff, delta, dR, dR2; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; + double wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; - double knfac, damp_normal, damp_normal_prefactor; - double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; - double fs, fs1, fs2, fs3; + double knfac, damp_normal, damp_normal_prefactor; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; + double fs, fs1, fs2, fs3; - double tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; + double tor1,tor2,tor3; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; - //For JKR - double R2, coh, F_pulloff, a, a2, E; - double t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3; + // for JKR + double R2, coh, F_pulloff, a, a2, E; + double t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3; - //Rolling - double k_roll, damp_roll; - double torroll1, torroll2, torroll3; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; + // rolling + double k_roll, damp_roll; + double torroll1, torroll2, torroll3; + double rollmag, rolldotn, scalefac; + double fr, fr1, fr2, fr3; - //Twisting - double k_twist, damp_twist, mu_twist; - double signtwist, magtwist, magtortwist, Mtcrit; - double tortwist1, tortwist2, tortwist3; + // twisting + double k_twist, damp_twist, mu_twist; + double signtwist, magtwist, magtortwist, Mtcrit; + double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht; - r = sqrt(rsq); - radsum = rwall + radius; + r = sqrt(rsq); + radsum = rwall + radius; - E = normal_coeffs[0]; + E = normal_coeffs[0]; - radsum = radius + rwall; - if (rwall == 0) Reff = radius; - else Reff = radius*rwall/(radius+rwall); + radsum = radius + rwall; + if (rwall == 0) Reff = radius; + else Reff = radius*rwall/(radius+rwall); - rinv = 1.0/r; + rinv = 1.0/r; - nx = dx*rinv; - ny = dy*rinv; - nz = dz*rinv; + nx = dx*rinv; + ny = dy*rinv; + nz = dz*rinv; - // relative translational velocity + // relative translational velocity - vr1 = v[0] - vwall[0]; - vr2 = v[1] - vwall[1]; - vr3 = v[2] - vwall[2]; + vr1 = v[0] - vwall[0]; + vr2 = v[1] - vwall[1]; + vr3 = v[2] - vwall[2]; - // normal component + // normal component - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; + vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n + vn1 = nx*vnnr; + vn2 = ny*vnnr; + vn3 = nz*vnnr; - delta = radsum - r; - dR = delta*Reff; - if (normal_model == JKR) { - history[0] = 1.0; - E *= THREEQUARTERS; - R2=Reff*Reff; - coh = normal_coeffs[3]; - dR2 = dR*dR; - t0 = coh*coh*R2*R2*E; - t1 = PI27SQ*t0; - t2 = 8*dR*dR2*E*E*E; - t3 = 4*dR2*E; - sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues - t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); - t5 = t3/t4 + t4/E; - sqrt2 = MAX(0, 2*dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); - a = INVROOT6*(t6 + sqrt(sqrt3)); - a2 = a*a; - knfac = normal_coeffs[0]*a; - Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); - } - else{ - knfac = E; //Hooke - a = sqrt(dR); - if (normal_model != HOOKE) { - Fne *= a; - knfac *= a; - } - Fne = knfac*delta; - if (normal_model == DMT) - Fne -= 4*MY_PI*normal_coeffs[3]*Reff; - } + delta = radsum - r; + dR = delta*Reff; + if (normal_model == JKR) { + history[0] = 1.0; + E *= THREEQUARTERS; + R2=Reff*Reff; + coh = normal_coeffs[3]; + dR2 = dR*dR; + t0 = coh*coh*R2*R2*E; + t1 = PI27SQ*t0; + t2 = 8*dR*dR2*E*E*E; + t3 = 4*dR2*E; + sqrt1 = MAX(0, t0*(t1+2*t2)); // in case sqrt(0) < 0 due to precision issues + t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); + t5 = t3/t4 + t4/E; + sqrt2 = MAX(0, 2*dR + t5); + t6 = sqrt(sqrt2); + sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = normal_coeffs[0]*a; + Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); + } + else{ + knfac = E; //Hooke + a = sqrt(dR); + if (normal_model != HOOKE) { + Fne *= a; + knfac *= a; + } + Fne = knfac*delta; + if (normal_model == DMT) + Fne -= 4*MY_PI*normal_coeffs[3]*Reff; + } - if (damping_model == VELOCITY) { - damp_normal = 1; - } - else if (damping_model == VISCOELASTIC) { - damp_normal = a*meff; - } - else if (damping_model == TSUJI) { - damp_normal = sqrt(meff*knfac); - } + if (damping_model == VELOCITY) { + damp_normal = 1; + } + else if (damping_model == VISCOELASTIC) { + damp_normal = a*meff; + } + else if (damping_model == TSUJI) { + damp_normal = sqrt(meff*knfac); + } - damp_normal_prefactor = normal_coeffs[1]*damp_normal; - Fdamp = -damp_normal_prefactor*vnnr; + damp_normal_prefactor = normal_coeffs[1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; - Fntot = Fne + Fdamp; + Fntot = Fne + Fdamp; - //**************************************** - //Tangential force, including history effects - //**************************************** + //**************************************** + // tangential force, including history effects + //**************************************** - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; + // tangential component + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; - // relative rotational velocity - wr1 = radius*omega[0] * rinv; - wr2 = radius*omega[1] * rinv; - wr3 = radius*omega[2] * rinv; + // relative rotational velocity + wr1 = radius*omega[0] * rinv; + wr2 = radius*omega[1] * rinv; + wr3 = radius*omega[2] * rinv; - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); + // relative tangential velocities + vtr1 = vt1 - (nz*wr2-ny*wr3); + vtr2 = vt2 - (nx*wr3-nz*wr1); + vtr3 = vt3 - (ny*wr1-nx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); - if (normal_model == JKR) { - F_pulloff = 3*M_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } - else if (normal_model == DMT) { - F_pulloff = 4*M_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } - else{ - Fncrit = fabs(Fntot); - } + if (normal_model == JKR) { + F_pulloff = 3*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } + else if (normal_model == DMT) { + F_pulloff = 4*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } + else{ + Fncrit = fabs(Fntot); + } - //------------------------------ - //Tangential forces - //------------------------------ - k_tangential = tangential_coeffs[0]; - damp_tangential = tangential_coeffs[1]*damp_normal_prefactor; + //------------------------------ + // tangential forces + //------------------------------ - int thist0 = tangential_history_index; - int thist1 = thist0 + 1; - int thist2 = thist1 + 1; + k_tangential = tangential_coeffs[0]; + damp_tangential = tangential_coeffs[1]*damp_normal_prefactor; - if (tangential_history) { - if (tangential_model == TANGENTIAL_MINDLIN) { - k_tangential *= a; - } - else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { - k_tangential *= a; - if (a < history[3]) { //On unloading, rescale the shear displacements - double factor = a/history[thist2+1]; - history[thist0] *= factor; - history[thist1] *= factor; - history[thist2] *= factor; - } - } - shrmag = sqrt(history[thist0]*history[thist0] + history[thist1]*history[thist1] + - history[thist2]*history[thist2]); + int thist0 = tangential_history_index; + int thist1 = thist0 + 1; + int thist2 = thist1 + 1; - // Rotate and update displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (history_update) { - rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0) { - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - history[thist0] -= rsht*nx; - history[thist1] -= rsht*ny; - history[thist2] -= rsht*nz; - //Also rescale to preserve magnitude - history[thist0] *= scalefac; - history[thist1] *= scalefac; - history[thist2] *= scalefac; - } - //Update history - history[thist0] += vtr1*dt; - history[thist1] += vtr2*dt; - history[thist2] += vtr3*dt; - } + if (tangential_history) { + if (tangential_model == TANGENTIAL_MINDLIN) { + k_tangential *= a; + } + else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + k_tangential *= a; + if (a < history[3]) { //On unloading, rescale the shear displacements + double factor = a/history[thist2+1]; + history[thist0] *= factor; + history[thist1] *= factor; + history[thist2] *= factor; + } + } + shrmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); - // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; + // rotate and update displacements. + // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 + if (history_update) { + rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz; + if (fabs(rsht) < EPSILON) rsht = 0; + if (rsht > 0) { + // if rhst == shrmag, contacting pair has rotated 90 deg in one step, + // in which case you deserve a crash! + scalefac = shrmag/(shrmag - rsht); + history[thist0] -= rsht*nx; + history[thist1] -= rsht*ny; + history[thist2] -= rsht*nz; + // also rescale to preserve magnitude + history[thist0] *= scalefac; + history[thist1] *= scalefac; + history[thist2] *= scalefac; + } + // update history + history[thist0] += vtr1*dt; + history[thist1] += vtr2*dt; + history[thist2] += vtr3*dt; + } - // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[2] * Fncrit; - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); - history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); - history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - } - else{ //Classic pair gran/hooke (no history) - fs = meff*damp_tangential*vrel; - if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; - else Ft = 0.0; - fs1 = -Ft*vtr1; - fs2 = -Ft*vtr2; - fs3 = -Ft*vtr3; - } + // tangential forces = history + tangential velocity damping + fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; - //**************************************** - // Rolling resistance - //**************************************** + // rescale frictional displacements and forces if needed + Fscrit = tangential_coeffs[2] * Fncrit; + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + if (fs > Fscrit) { + if (shrmag != 0.0) { + history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + + damp_tangential*vtr1); + history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs + + damp_tangential*vtr2); + history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs + + damp_tangential*vtr3); + fs1 *= Fscrit/fs; + fs2 *= Fscrit/fs; + fs3 *= Fscrit/fs; + } else fs1 = fs2 = fs3 = 0.0; + } + } else { // classic pair gran/hooke (no history) + fs = meff*damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + } - if (roll_model != ROLL_NONE) { - relrot1 = omega[0]; - relrot2 = omega[1]; - relrot3 = omega[2]; + //**************************************** + // rolling resistance + //**************************************** - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; + if (roll_model != ROLL_NONE) { + relrot1 = omega[0]; + relrot2 = omega[1]; + relrot3 = omega[2]; - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; + // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // This is different from the Marshall papers, + // which use the Bagi/Kuhn formulation + // for rolling velocity (see Wang et al for why the latter is wrong) + vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; + vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; + vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - // Rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + // rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); - if (history_update) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - //Also rescale to preserve magnitude - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; - } - history[rhist0] += vrl1*dt; - history[rhist1] += vrl2*dt; - history[rhist2] += vrl3*dt; - } + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - k_roll = roll_coeffs[0]; - damp_roll = roll_coeffs[1]; - fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; - fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; - fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; + if (history_update) { + if (fabs(rolldotn) < EPSILON) rolldotn = 0; + if (rolldotn > 0) { // rotate into tangential plane + scalefac = rollmag/(rollmag - rolldotn); + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; + } + history[rhist0] += vrl1*dt; + history[rhist1] += vrl2*dt; + history[rhist2] += vrl3*dt; + } - // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[2] * Fncrit; + k_roll = roll_coeffs[0]; + damp_roll = roll_coeffs[1]; + fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; + fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; + fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); - history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); - history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - } + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[2] * Fncrit; - //**************************************** - // Twisting torque, including history effects - //**************************************** - if (twist_model != TWIST_NONE) { - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist_model == TWIST_MARSHALL) { - k_twist = 0.5*k_tangential*a*a;; //eq 32 of Marshall paper - damp_twist = 0.5*damp_tangential*a*a; - mu_twist = TWOTHIRDS*a*tangential_coeffs[2]; - } - else{ - k_twist = twist_coeffs[0]; - damp_twist = twist_coeffs[1]; - mu_twist = twist_coeffs[2]; - } - if (history_update) { - history[twist_history_index] += magtwist*dt; - } - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = mu_twist*Fncrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - } - // Apply forces & torques + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + if (rollmag != 0.0) { + history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); + history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); + history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); + fr1 *= Frcrit/fr; + fr2 *= Frcrit/fr; + fr3 *= Frcrit/fr; + } else fr1 = fr2 = fr3 = 0.0; + } + } - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; + //**************************************** + // twisting torque, including history effects + //**************************************** - if (peratom_flag) { - contact[1] = fx; - contact[2] = fy; - contact[3] = fz; - } + if (twist_model != TWIST_NONE) { + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + if (twist_model == TWIST_MARSHALL) { + k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper + damp_twist = 0.5*damp_tangential*a*a; + mu_twist = TWOTHIRDS*a*tangential_coeffs[2]; + } + else{ + k_twist = twist_coeffs[0]; + damp_twist = twist_coeffs[1]; + mu_twist = twist_coeffs[2]; + } + if (history_update) { + history[twist_history_index] += magtwist*dt; + } + // M_t torque (eq 30) + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist; + signtwist = (magtwist > 0) - (magtwist < 0); + Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) + if (fabs(magtortwist) > Mtcrit) { + history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - + damp_twist*magtwist); + magtortwist = -Mtcrit * signtwist; // eq 34 + } + } - f[0] += fx; - f[1] += fy; - f[2] += fz; + // apply forces & torques - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; + fx = nx*Fntot + fs1; + fy = ny*Fntot + fs2; + fz = nz*Fntot + fs3; - torque[0] -= radius*tor1; - torque[1] -= radius*tor2; - torque[2] -= radius*tor3; + if (peratom_flag) { + contact[1] = fx; + contact[2] = fy; + contact[3] = fz; + } - if (twist_model != TWIST_NONE) { - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; + f[0] += fx; + f[1] += fy; + f[2] += fz; - torque[0] += tortwist1; - torque[1] += tortwist2; - torque[2] += tortwist3; - } + tor1 = ny*fs3 - nz*fs2; + tor2 = nz*fs1 - nx*fs3; + tor3 = nx*fs2 - ny*fs1; - if (roll_model != ROLL_NONE) { - torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = Reff*(nz*fr1 - nx*fr3); - torroll3 = Reff*(nx*fr2 - ny*fr1); + torque[0] -= radius*tor1; + torque[1] -= radius*tor2; + torque[2] -= radius*tor3; - torque[0] += torroll1; - torque[1] += torroll2; - torque[2] += torroll3; - } + if (twist_model != TWIST_NONE) { + tortwist1 = magtortwist * nx; + tortwist2 = magtortwist * ny; + tortwist3 = magtortwist * nz; + + torque[0] += tortwist1; + torque[1] += tortwist2; + torque[2] += tortwist3; + } + + if (roll_model != ROLL_NONE) { + torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr + torroll2 = Reff*(nz*fr1 - nx*fr3); + torroll3 = Reff*(nx*fr2 - ny*fr1); + + torque[0] += torroll1; + torque[1] += torroll2; + torque[2] += torroll3; + } } - - /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ @@ -1389,9 +1437,10 @@ double FixWallGran::memory_usage() { int nmax = atom->nmax; double bytes = 0.0; - if (use_history) bytes += nmax*size_history * sizeof(double); // shear history - if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid - if (peratom_flag) bytes += nmax*size_peratom_cols*sizeof(double); //store contacts + if (use_history) bytes += nmax*size_history * sizeof(double); // shear history + if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid + // store contacts + if (peratom_flag) bytes += nmax*size_peratom_cols*sizeof(double); return bytes; } @@ -1401,7 +1450,8 @@ double FixWallGran::memory_usage() void FixWallGran::grow_arrays(int nmax) { - if (use_history) memory->grow(history_one,nmax,size_history,"fix_wall_gran:history_one"); + if (use_history) memory->grow(history_one,nmax,size_history, + "fix_wall_gran:history_one"); if (peratom_flag) { memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); } diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index a81cdcb6c8..ee81477ddb 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -66,26 +66,24 @@ class FixWallGran : public Fix { bigint time_origin; double kn,kt,gamman,gammat,xmu; - //For granular - //Model choices + // for granular model choices int normal_model, damping_model; int tangential_model, roll_model, twist_model; - int beyond_contact; - //History flags + // history flags int normal_history, tangential_history, roll_history, twist_history; - //Indices of history entries + // indices of history entries int normal_history_index; int tangential_history_index; int roll_history_index; int twist_history_index; - //Material coefficients + // material coefficients double Emod, poiss, Gmod; - //Contact model coefficients + // contact model coefficients double normal_coeffs[4]; double tangential_coeffs[3]; double roll_coeffs[3]; @@ -97,7 +95,7 @@ class FixWallGran : public Fix { char *idregion; int use_history; // if particle/wall interaction stores history - int history_update; // flag for whether shear history is updated + int history_update; // flag for whether shear history is updated int size_history; // # of shear history values per contact // shear history for single contact per particle @@ -110,7 +108,8 @@ class FixWallGran : public Fix { double *mass_rigid; // rigid mass for owned+ghost atoms int nmax; // allocated size of mass_rigid - // Store particle interactions + // store particle interactions + int store; }; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 95b34e0929..e6f8406be0 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -47,8 +47,9 @@ enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR}; /* ---------------------------------------------------------------------- */ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : - FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), ncontact(NULL), - walls(NULL), history_many(NULL), c2r(NULL) + FixWallGran(lmp, narg, arg), region(NULL), region_style(NULL), + ncontact(NULL), + walls(NULL), history_many(NULL), c2r(NULL) { restart_global = 1; motion_resetflag = 0; @@ -190,7 +191,8 @@ void FixWallGranRegion::post_force(int /*vflag*/) if (!region->match(x[i][0],x[i][1],x[i][2])) continue; if (pairstyle == GRANULAR && normal_model == JKR){ - nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]+pulloff_distance(radius[i])); + nc = region->surface(x[i][0],x[i][1],x[i][2], + radius[i]+pulloff_distance(radius[i])); } else{ nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]); @@ -278,8 +280,9 @@ void FixWallGranRegion::post_force(int /*vflag*/) v[i],f[i],omega[i],torque[i], radius[i],meff,history_many[i][c2r[ic]], contact); else if (pairstyle == GRANULAR) - granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius, v[i],f[i],omega[i],torque[i], radius[i],meff,history_many[i][c2r[ic]], contact); - + granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius, + v[i],f[i],omega[i],torque[i], + radius[i],meff,history_many[i][c2r[ic]],contact); } } } @@ -364,11 +367,11 @@ void FixWallGranRegion::grow_arrays(int nmax) if (use_history) { memory->grow(ncontact,nmax,"fix_wall_gran:ncontact"); memory->grow(walls,nmax,tmax,"fix_wall_gran:walls"); - memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many"); + memory->grow(history_many,nmax,tmax,size_history, + "fix_wall_gran:history_many"); } - if (peratom_flag){ + if (peratom_flag) memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); - } } /* ---------------------------------------------------------------------- diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 344e72f8ef..c86c2b0c90 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -59,7 +59,9 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) comm_forward = 1; - nondefault_history_transfer = 0; //keep default behavior of history[i][j] = -history[j][i] + // keep default behavior of history[i][j] = -history[j][i] + + nondefault_history_transfer = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index 2cb609fd82..b1bb212f89 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -26,6 +26,8 @@ namespace LAMMPS_NS { class PairGranHookeHistory : public Pair { public: + int nondefault_history_transfer; + PairGranHookeHistory(class LAMMPS *); virtual ~PairGranHookeHistory(); virtual void compute(int, int); @@ -42,7 +44,6 @@ class PairGranHookeHistory : public Pair { int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); double memory_usage(); - int nondefault_history_transfer; protected: double kn,kt,gamman,gammat,xmu; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index b3046788e0..c8450c5434 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -11,9 +11,9 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- -Contributing authors: -Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) -Leo Silbert (SNL), Gary Grest (SNL) + Contributing authors: + Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) + Leo Silbert (SNL), Gary Grest (SNL) ----------------------------------------------------------------------- */ #include @@ -52,7 +52,8 @@ using namespace MathConst; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_SDS}; @@ -90,10 +91,10 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) nondefault_history_transfer = 0; tangential_history_index = 0; roll_history_index = twist_history_index = 0; - } /* ---------------------------------------------------------------------- */ + PairGranular::~PairGranular() { delete [] svector; @@ -118,16 +119,17 @@ PairGranular::~PairGranular() memory->destroy(roll_model); memory->destroy(twist_model); - - delete [] onerad_dynamic; delete [] onerad_frozen; delete [] maxrad_dynamic; delete [] maxrad_frozen; } + memory->destroy(mass_rigid); } +/* ---------------------------------------------------------------------- */ + void PairGranular::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; @@ -147,18 +149,18 @@ void PairGranular::compute(int eflag, int vflag) double mi,mj,meff; double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; - // For JKR + // for JKR double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; double t0, t1, t2, t3, t4, t5, t6; double sqrt1, sqrt2, sqrt3; - // Rolling + // rolling double k_roll, damp_roll; double torroll1, torroll2, torroll3; double rollmag, rolldotn, scalefac; double fr, fr1, fr2, fr3; - // Twisting + // twisting double k_twist, damp_twist, mu_twist; double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; @@ -317,7 +319,8 @@ void PairGranular::compute(int eflag, int vflag) t1 = PI27SQ*t0; t2 = 8*dR*dR2*E*E*E; t3 = 4*dR2*E; - sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues + // in case sqrt(0) < 0 due to precision issues + sqrt1 = MAX(0, t0*(t1+2*t2)); t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); t5 = t3/t4 + t4/E; sqrt2 = MAX(0, 2*dR + t5); @@ -328,7 +331,7 @@ void PairGranular::compute(int eflag, int vflag) knfac = normal_coeffs[itype][jtype][0]*a; Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); } else { - knfac = E; //Hooke + knfac = E; // Hooke Fne = knfac*delta; a = sqrt(dR); if (normal_model[itype][jtype] != HOOKE) { @@ -339,7 +342,9 @@ void PairGranular::compute(int eflag, int vflag) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } - //Consider restricting Hooke to only have 'velocity' as an option for damping? + // NOTE: consider restricting Hooke to only have + // 'velocity' as an option for damping? + if (damping_model[itype][jtype] == VELOCITY) { damp_normal = 1; } else if (damping_model[itype][jtype] == VISCOELASTIC) { @@ -354,7 +359,7 @@ void PairGranular::compute(int eflag, int vflag) Fntot = Fne + Fdamp; //**************************************** - //Tangential force, including history effects + // tangential force, including history effects //**************************************** // tangential component @@ -374,7 +379,7 @@ void PairGranular::compute(int eflag, int vflag) vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; vrel = sqrt(vrel); - // If any history is needed: + // if any history is needed if (use_history) { touch[jj] = 1; history = &allhistory[size_history*jj]; @@ -391,45 +396,51 @@ void PairGranular::compute(int eflag, int vflag) } //------------------------------ - //Tangential forces + // tangential forces //------------------------------ k_tangential = tangential_coeffs[itype][jtype][0]; - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; + damp_tangential = tangential_coeffs[itype][jtype][1] * + damp_normal_prefactor; if (tangential_history) { if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { k_tangential *= a; - } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + } else if (tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE) { k_tangential *= a; - if (a < history[3]) { //On unloading, rescale the shear displacements + // on unloading, rescale the shear displacements + if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; history[1] *= factor; history[2] *= factor; } } - // Rotate and update displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 + // rotate and update displacements. + // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (historyupdate) { rsht = history[0]*nx + history[1]*ny + history[2]*nz; if (fabs(rsht) < EPSILON) rsht = 0; if (rsht > 0) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); - scalefac = shrmag/(shrmag - rsht); //if rsht == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! + // if rsht == shrmag, contacting pair has rotated 90 deg + // in one step, in which case you deserve a crash! + scalefac = shrmag/(shrmag - rsht); history[0] -= rsht*nx; history[1] -= rsht*ny; history[2] -= rsht*nz; - //Also rescale to preserve magnitude + // also rescale to preserve magnitude history[0] *= scalefac; history[1] *= scalefac; history[2] *= scalefac; } - //Update history + // update history history[0] += vtr1*dt; history[1] += vtr2*dt; history[2] += vtr3*dt; - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) history[3] = a; + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) + history[3] = a; } // tangential forces = history + tangential velocity damping @@ -444,15 +455,18 @@ void PairGranular::compute(int eflag, int vflag) shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); - history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); - history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + + damp_tangential*vtr1); + history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + + damp_tangential*vtr2); + history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + + damp_tangential*vtr3); fs1 *= Fscrit/fs; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; } else fs1 = fs2 = fs3 = 0.0; } - } else{ //Classic pair gran/hooke (no history) + } else { // classic pair gran/hooke (no history) fs = meff*damp_tangential*vrel; if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; else Ft = 0.0; @@ -462,7 +476,7 @@ void PairGranular::compute(int eflag, int vflag) } //**************************************** - // Rolling resistance + // rolling resistance //**************************************** if (roll_model[itype][jtype] != ROLL_NONE) { @@ -470,12 +484,18 @@ void PairGranular::compute(int eflag, int vflag) relrot2 = omega[i][1] - omega[j][1]; relrot3 = omega[i][2] - omega[j][2]; - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation + // rolling velocity, + // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // this is different from the Marshall papers, + // which use the Bagi/Kuhn formulation // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; + // - 0.5*((radj-radi)/radsum)*vtr1; + // - 0.5*((radj-radi)/radsum)*vtr2; + // - 0.5*((radj-radi)/radsum)*vtr3; + + vrl1 = Reff*(relrot2*nz - relrot3*ny); + vrl2 = Reff*(relrot3*nx - relrot1*nz); + vrl3 = Reff*(relrot1*ny - relrot2*nx); int rhist0 = roll_history_index; int rhist1 = rhist0 + 1; @@ -484,7 +504,7 @@ void PairGranular::compute(int eflag, int vflag) rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; if (historyupdate) { if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { //Rotate into tangential plane + if (rolldotn > 0) { // rotate into tangential plane rollmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + history[rhist2]*history[rhist2]); @@ -492,7 +512,7 @@ void PairGranular::compute(int eflag, int vflag) history[rhist0] -= rolldotn*nx; history[rhist1] -= rolldotn*ny; history[rhist2] -= rolldotn*nz; - //Also rescale to preserve magnitude + // also rescale to preserve magnitude history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; @@ -528,12 +548,14 @@ void PairGranular::compute(int eflag, int vflag) } //**************************************** - // Twisting torque, including history effects + // twisting torque, including history effects //**************************************** + if (twist_model[itype][jtype] != TWIST_NONE) { - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + // omega_T (eq 29 of Marshall) + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; if (twist_model[itype][jtype] == TWIST_MARSHALL) { - k_twist = 0.5*k_tangential*a*a;; //eq 32 of Marshall paper + k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper damp_twist = 0.5*damp_tangential*a*a; mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2]; } else { @@ -544,15 +566,18 @@ void PairGranular::compute(int eflag, int vflag) if (historyupdate) { history[twist_history_index] += magtwist*dt; } - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) + magtortwist = -k_twist*history[twist_history_index] - + damp_twist*magtwist; // M_t torque (eq 30) signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = mu_twist*Fncrit;//critical torque (eq 44) + Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) if (fabs(magtortwist) > Mtcrit) { - history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 + history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - + damp_twist*magtwist); + magtortwist = -Mtcrit * signtwist; // eq 34 } } - // Apply forces & torques + + // apply forces & torques fx = nx*Fntot + fs1; fy = ny*Fntot + fs2; @@ -582,7 +607,7 @@ void PairGranular::compute(int eflag, int vflag) } if (roll_model[itype][jtype] != ROLL_NONE) { - torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr + torroll1 = Reff*(ny*fr3 - nz*fr2); // n cross fr torroll2 = Reff*(nz*fr1 - nx*fr3); torroll3 = Reff*(nx*fr2 - ny*fr1); @@ -619,9 +644,8 @@ void PairGranular::compute(int eflag, int vflag) } } - /* ---------------------------------------------------------------------- -allocate all arrays + allocate all arrays ------------------------------------------------------------------------- */ void PairGranular::allocate() @@ -657,7 +681,7 @@ void PairGranular::allocate() } /* ---------------------------------------------------------------------- - global settings + global settings ------------------------------------------------------------------------- */ void PairGranular::settings(int narg, char **arg) @@ -665,7 +689,7 @@ void PairGranular::settings(int narg, char **arg) if (narg == 1) { cutoff_global = force->numeric(FLERR,arg[0]); } else { - cutoff_global = -1; //Will be set based on particle sizes, model choice + cutoff_global = -1; // will be set based on particle sizes, model choice } normal_history = tangential_history = 0; @@ -705,45 +729,57 @@ void PairGranular::coeff(int narg, char **arg) int iarg = 2; while (iarg < narg) { if (strcmp(arg[iarg], "hooke") == 0) { - if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option"); + if (iarg + 2 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hooke option"); normal_model_one = HOOKE; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // kn + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping iarg += 3; } else if (strcmp(arg[iarg], "hertz") == 0) { int num_coeffs = 2; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); + if (iarg + num_coeffs >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hertz option"); normal_model_one = HERTZ; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // kn + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "hertz/material") == 0) { int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz/material option"); + if (iarg + num_coeffs >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hertz/material option"); normal_model_one = HERTZ_MATERIAL; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "dmt") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hertz option"); normal_model_one = DMT; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio - normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio + normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); // cohesion iarg += 5; } else if (strcmp(arg[iarg], "jkr") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for JKR option"); beyond_contact = 1; normal_model_one = JKR; - normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio - normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); // E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); // damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); // Poisson's ratio + normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); // cohesion iarg += 5; } else if (strcmp(arg[iarg], "damping") == 0) { - if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); + if (iarg+1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, " + "not enough parameters provided for damping model"); if (strcmp(arg[iarg+1], "velocity") == 0) { damping_model_one = VELOCITY; iarg += 1; @@ -753,58 +789,80 @@ void PairGranular::coeff(int narg, char **arg) } else if (strcmp(arg[iarg+1], "tsuji") == 0) { damping_model_one = TSUJI; iarg += 1; - } else error->all(FLERR, "Illegal pair_coeff command, unrecognized damping model"); + } else error->all(FLERR, "Illegal pair_coeff command, " + "unrecognized damping model"); iarg += 1; } else if (strcmp(arg[iarg], "tangential") == 0) { - if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_coeff command, must specify tangential model after 'tangential' keyword"); + if (iarg + 1 >= narg) + error->all(FLERR,"Illegal pair_coeff command, must specify " + "tangential model after tangential keyword"); if (strcmp(arg[iarg+1], "linear_nohistory") == 0) { - if (iarg + 3 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (iarg + 3 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); tangential_model_one = TANGENTIAL_NOHISTORY; tangential_coeffs_one[0] = 0; - tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //gammat - tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //friction coeff. + // gammat and friction coeff + tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); + tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) || (strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strcmp(arg[iarg+1], "linear_history") == 0) tangential_model_one = TANGENTIAL_HISTORY; - else if (strcmp(arg[iarg+1], "mindlin") == 0) tangential_model_one = TANGENTIAL_MINDLIN; - else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); + if (strcmp(arg[iarg+1], "linear_history") == 0) + tangential_model_one = TANGENTIAL_HISTORY; + else if (strcmp(arg[iarg+1], "mindlin") == 0) + tangential_model_one = TANGENTIAL_MINDLIN; + else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; tangential_history = 1; - if ((tangential_model_one == TANGENTIAL_MINDLIN || tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) && + if ((tangential_model_one == TANGENTIAL_MINDLIN || + tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) && (strcmp(arg[iarg+2], "NULL") == 0)) { if (normal_model_one == HERTZ || normal_model_one == HOOKE) { - error->all(FLERR, "NULL setting for Mindlin tangential stiffness requires a normal contact model that specifies material properties"); + error->all(FLERR, "NULL setting for Mindlin tangential " + "stiffness requires a normal contact model that " + "specifies material properties"); } tangential_coeffs_one[0] = -1; } else { - tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt + tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); // kt } - tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + // gammat and friction coeff + tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); + tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else { - error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); + error->all(FLERR, "Illegal pair_coeff command, " + "tangential model not recognized"); } } else if (strcmp(arg[iarg], "rolling") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0) { roll_model_one = ROLL_NONE; iarg += 2; } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for rolling model"); roll_model_one = ROLL_SDS; roll_history = 1; - roll_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kR - roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. + // kR and gammaR and rolling friction coeff + roll_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); + roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); + roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else { - error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized"); + error->all(FLERR, "Illegal pair_coeff command, " + "rolling friction model not recognized"); } } else if (strcmp(arg[iarg], "twisting") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0) { twist_model_one = TWIST_NONE; iarg += 2; @@ -813,24 +871,31 @@ void PairGranular::coeff(int narg, char **arg) twist_history = 1; iarg += 2; } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for twist model"); twist_model_one = TWIST_SDS; twist_history = 1; - twist_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt - twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + // kt and gammat and friction coeff + twist_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); + twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); + twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else { - error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized"); + error->all(FLERR, "Illegal pair_coeff command, " + "twisting friction model not recognized"); } } else if (strcmp(arg[iarg], "cutoff") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); cutoff_one = force->numeric(FLERR,arg[iarg+1]); } else error->all(FLERR, "Illegal pair coeff command"); } - //It is an error not to specify normal or tangential model - if ((normal_model_one < 0) || (tangential_model_one < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); + // error not to specify normal or tangential model + if ((normal_model_one < 0) || (tangential_model_one < 0)) + error->all(FLERR, "Illegal pair coeff command, " + "must specify normal or tangential contact model"); int count = 0; double damp; @@ -849,7 +914,9 @@ void PairGranular::coeff(int narg, char **arg) if (normal_model_one != HERTZ && normal_model_one != HOOKE) { Emod[i][j] = Emod[j][i] = normal_coeffs_one[0]; poiss[i][j] = poiss[j][i] = normal_coeffs_one[2]; - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = FOURTHIRDS*mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]); + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + FOURTHIRDS*mix_stiffnessE(Emod[i][j],Emod[i][j], + poiss[i][j],poiss[i][j]); } else { normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0]; } @@ -860,12 +927,15 @@ void PairGranular::coeff(int narg, char **arg) tangential_model[i][j] = tangential_model[j][i] = tangential_model_one; if (tangential_coeffs_one[0] == -1) { - tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = 8*mix_stiffnessG(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]); + tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = + 8*mix_stiffnessG(Emod[i][j],Emod[i][j],poiss[i][j],poiss[i][j]); } else { - tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = tangential_coeffs_one[0]; + tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = + tangential_coeffs_one[0]; } for (int k = 1; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_one[k]; + tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = + tangential_coeffs_one[k]; roll_model[i][j] = roll_model[j][i] = roll_model_one; if (roll_model_one != ROLL_NONE) @@ -877,13 +947,13 @@ void PairGranular::coeff(int narg, char **arg) for (int k = 0; k < 3; k++) twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k]; - cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; setflag[i][j] = 1; count++; } } + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -902,17 +972,21 @@ void PairGranular::init_style() if (comm->ghost_velocity == 0) error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - // Determine whether we need a granular neigh list, how large it needs to be - use_history = normal_history || tangential_history || roll_history || twist_history; + // determine whether we need a granular neigh list, how large it needs to be + + use_history = normal_history || tangential_history || + roll_history || twist_history; + + // for JKR, will need fix/neigh/history to keep track of touch arrays - //For JKR, will need fix/neigh/history to keep track of touch arrays for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) if (normal_model[i][j] == JKR) use_history = 1; size_history = 3*tangential_history + 3*roll_history + twist_history; - //Determine location of tangential/roll/twist histories in array + // determine location of tangential/roll/twist histories in array + if (roll_history) { if (tangential_history) roll_history_index = 3; else roll_history_index = 0; @@ -1036,7 +1110,7 @@ void PairGranular::init_style() } /* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i + init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairGranular::init_one(int i, int j) @@ -1051,46 +1125,58 @@ double PairGranular::init_one(int i, int j) (twist_model[i][i] != twist_model[j][j])) { char str[512]; - sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j); + sprintf(str,"Granular pair style functional forms are different, " + "cannot mix coefficients for types %d and %d. \n" + "This combination must be set explicitly " + "via pair_coeff command",i,j); error->one(FLERR,str); } if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE) - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); else - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]); + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]); - normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); + normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = + mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); if ((normal_model[i][j] == JKR) || (normal_model[i][j] == DMT)) - normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]); + normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = + mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]); for (int k = 0; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); + tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = + mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); if (roll_model[i][j] != ROLL_NONE) { for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); + roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = + mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); } if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL) { for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); + twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = + mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); } } - // It is possible that cut[i][j] at this point is still 0.0. This can happen when + // It is possible that cut[i][j] at this point is still 0.0. + // This can happen when // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size // To avoid this issue, for cases involving cut[i][j] = 0.0 (possible only // if there is no current information about radius/cutoff of type i and j). // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0. + double pulloff; if (cutoff_type[i][j] < 0 && cutoff_global < 0) { if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - + // radius info about both i and j exist + ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; pulloff = 0.0; if (normal_model[i][j] == JKR) { @@ -1105,7 +1191,13 @@ double PairGranular::init_one(int i, int j) if (normal_model[i][j] == JKR) pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j); cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff); - } else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) + + } else { + + // radius info about either i or j does not exist + // (i.e. not present and not about to get poured; + // set to largest value to not interfere with neighbor list) + double cutmax = 0.0; for (int k = 1; k <= atom->ntypes; k++) { cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); @@ -1123,8 +1215,8 @@ double PairGranular::init_one(int i, int j) } /* ---------------------------------------------------------------------- - proc 0 writes to restart file - ------------------------------------------------------------------------- */ + proc 0 writes to restart file +------------------------------------------------------------------------- */ void PairGranular::write_restart(FILE *fp) { @@ -1149,8 +1241,8 @@ void PairGranular::write_restart(FILE *fp) } /* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts - ------------------------------------------------------------------------- */ + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ void PairGranular::read_restart(FILE *fp) { @@ -1189,7 +1281,6 @@ void PairGranular::read_restart(FILE *fp) } } - /* ---------------------------------------------------------------------- */ void PairGranular::reset_dt() @@ -1216,18 +1307,17 @@ double PairGranular::single(int i, int j, int itype, int jtype, double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; double fs, fs1, fs2, fs3; - //For JKR + // for JKR double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; double delta, t0, t1, t2, t3, t4, t5, t6; double sqrt1, sqrt2, sqrt3; - - //Rolling + // rolling double k_roll, damp_roll; double rollmag; double fr, fr1, fr2, fr3; - //Twisting + // twisting double k_twist, damp_twist, mu_twist; double signtwist, magtwist, magtortwist, Mtcrit; @@ -1338,7 +1428,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, t1 = PI27SQ*t0; t2 = 8*dR*dR2*E*E*E; t3 = 4*dR2*E; - sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues + // in case sqrt(0) < 0 due to precision issues + sqrt1 = MAX(0, t0*(t1+2*t2)); t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); t5 = t3/t4 + t4/E; sqrt2 = MAX(0, 2*dR + t5); @@ -1387,7 +1478,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, } //**************************************** - //Tangential force, including history effects + // tangential force, including history effects //**************************************** // tangential component @@ -1418,7 +1509,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, } //------------------------------ - //Tangential forces + // tangential forces //------------------------------ k_tangential = tangential_coeffs[itype][jtype][0]; damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; @@ -1428,7 +1519,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, k_tangential *= a; } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { k_tangential *= a; - if (a < history[3]) { //On unloading, rescale the shear displacements + // on unloading, rescale the shear displacements + if (a < history[3]) { double factor = a/history[3]; history[0] *= factor; history[1] *= factor; @@ -1457,7 +1549,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs3 *= Fscrit/fs; } else fs1 = fs2 = fs3 = 0.0; } - } else { //Classic pair gran/hooke (no history) + + // classic pair gran/hooke (no history) + } else { fs = meff*damp_tangential*vrel; if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; else Ft = 0.0; @@ -1467,7 +1561,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, } //**************************************** - // Rolling resistance + // rolling resistance //**************************************** if (roll_model[itype][jtype] != ROLL_NONE) { @@ -1476,8 +1570,10 @@ double PairGranular::single(int i, int j, int itype, int jtype, relrot3 = omega[i][2] - omega[j][2]; // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation + // this is different from the Marshall papers, + // which use the Bagi/Kuhn formulation // for rolling velocity (see Wang et al for why the latter is wrong) + vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; @@ -1486,7 +1582,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, int rhist1 = rhist0 + 1; int rhist2 = rhist1 + 1; - // Rolling displacement + // rolling displacement rollmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + history[rhist2]*history[rhist2]); @@ -1512,10 +1608,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, } //**************************************** - // Twisting torque, including history effects + // twisting torque, including history effects //**************************************** + if (twist_model[itype][jtype] != TWIST_NONE) { - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + // omega_T (eq 29 of Marshall) + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; if (twist_model[itype][jtype] == TWIST_MARSHALL) { k_twist = 0.5*k_tangential*a*a;; //eq 32 damp_twist = 0.5*damp_tangential*a*a; @@ -1525,11 +1623,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, damp_twist = twist_coeffs[itype][jtype][1]; mu_twist = twist_coeffs[itype][jtype][2]; } - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) + // M_t torque (eq 30) + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist; signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = mu_twist*Fncrit;//critical torque (eq 44) + Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) if (fabs(magtortwist) > Mtcrit) { - magtortwist = -Mtcrit * signtwist; //eq 34 + magtortwist = -Mtcrit * signtwist; // eq 34 } } @@ -1578,8 +1677,8 @@ void PairGranular::unpack_forward_comm(int n, int first, double *buf) } /* ---------------------------------------------------------------------- - memory usage of local atom-based arrays - ------------------------------------------------------------------------- */ + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ double PairGranular::memory_usage() { @@ -1588,25 +1687,27 @@ double PairGranular::memory_usage() } /* ---------------------------------------------------------------------- - mixing of Young's modulus (E) + mixing of Young's modulus (E) ------------------------------------------------------------------------- */ -double PairGranular::mix_stiffnessE(double Eii, double Ejj, double poisii, double poisjj) +double PairGranular::mix_stiffnessE(double Eii, double Ejj, + double poisii, double poisjj) { return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj); } /* ---------------------------------------------------------------------- - mixing of shear modulus (G) + mixing of shear modulus (G) ------------------------------------------------------------------------ */ -double PairGranular::mix_stiffnessG(double Eii, double Ejj, double poisii, double poisjj) +double PairGranular::mix_stiffnessG(double Eii, double Ejj, + double poisii, double poisjj) { return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj)); } /* ---------------------------------------------------------------------- - mixing of everything else + mixing of everything else ------------------------------------------------------------------------- */ double PairGranular::mix_geom(double valii, double valjj) @@ -1616,10 +1717,11 @@ double PairGranular::mix_geom(double valii, double valjj) /* ---------------------------------------------------------------------- - Compute pull-off distance (beyond contact) for a given radius and atom type + compute pull-off distance (beyond contact) for a given radius and atom type ------------------------------------------------------------------------- */ -double PairGranular::pulloff_distance(double radi, double radj, int itype, int jtype) +double PairGranular::pulloff_distance(double radi, double radj, + int itype, int jtype) { double E, coh, a, Reff; Reff = radi*radj/(radi+radj); @@ -1631,11 +1733,12 @@ double PairGranular::pulloff_distance(double radi, double radj, int itype, int j } /* ---------------------------------------------------------------------- - Transfer history during fix/neigh/history exchange - Only needed if any history entries i-j are not just negative of j-i entries + transfer history during fix/neigh/history exchange + only needed if any history entries i-j are not just negative of j-i entries ------------------------------------------------------------------------- */ -void PairGranular::transfer_history(double* source, double* target) { + +void PairGranular::transfer_history(double* source, double* target) +{ for (int i = 0; i < size_history; i++) target[i] = history_transfer_factors[i]*source[i]; } - diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 4743d271f5..935e676487 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -66,29 +66,29 @@ class PairGranular : public Pair { int size_history; int *history_transfer_factors; - //Model choices + // model choices int **normal_model, **damping_model; int **tangential_model, **roll_model, **twist_model; - //History flags + // history flags int normal_history, tangential_history, roll_history, twist_history; - //Indices of history entries + // indices of history entries int normal_history_index; int tangential_history_index; int roll_history_index; int twist_history_index; - //Per-type material coefficients + // per-type material coefficients double **Emod, **poiss, **Gmod; - //Per-type coefficients, set in pair coeff command + // per-type coefficients, set in pair coeff command double ***normal_coeffs; double ***tangential_coeffs; double ***roll_coeffs; double ***twist_coeffs; - //Optional user-specified global cutoff, per-type user-specified cutoffs + // optional user-specified global cutoff, per-type user-specified cutoffs double **cutoff_type; double cutoff_global; diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index e33ebe57dc..207c409596 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -408,7 +408,8 @@ void FixNeighHistory::pre_exchange_newton() m = npartner[j]++; partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; - if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues); + if (pair->nondefault_history_transfer) + pair->transfer_history(onevalues,jvalues); else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } @@ -521,7 +522,8 @@ void FixNeighHistory::pre_exchange_no_newton() m = npartner[j]++; partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; - if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues); + if (pair->nondefault_history_transfer) + pair->transfer_history(onevalues, jvalues); else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } diff --git a/src/pair.h b/src/pair.h index 0911bff706..035ebc8773 100644 --- a/src/pair.h +++ b/src/pair.h @@ -98,7 +98,7 @@ class Pair : protected Pointers { enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options - int beyond_contact, nondefault_history_transfer; //for granular styles + int beyond_contact, nondefault_history_transfer; // for granular styles // KOKKOS host/device flag and data masks