From d8e8a0d2d24fcbbf26b5b60fb876ea3d85ba023b Mon Sep 17 00:00:00 2001 From: "Dan S. Bolintineanu" Date: Mon, 18 Feb 2019 09:58:34 -0700 Subject: [PATCH] More changes to pair granular: - tangential damping now set by scaling the normal damping - some fixes to the twisting coefficients for the Marshall twist model - progress (completion?) of doc page --- doc/src/pair_granular.txt | 529 ++++++++++++++++++++++++++++----- src/GRANULAR/pair_granular.cpp | 364 +++++++++++------------ 2 files changed, 619 insertions(+), 274 deletions(-) diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index d6b18217d8..ff2ee94be0 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -17,69 +17,79 @@ pair_style granular command :h3 [Syntax:] -pair_style style cutoff :pre +pair_style granular cutoff :pre -style = {granular} or {granular/multi} :ulb,l cutoff = global cutoff (optional). See discussion below. :l :ule [Examples:] pair_style granular -pair_coeff * * hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.4 :pre +pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 :pre pair_style granular -pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall -pair_coeff 2 2 hertz 200.0 20.0 tangential mindlin 300.0 50.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre +pair_coeff * * hertz 1000.0 50.0 tangential linear_history 800.0 1.0 0.4 :pre pair_style granular -pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall -pair_coeff 2 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall -pair_coeff 1 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre +pair_coeff * * hertz 1000.0 0.3 tangential linear_history 800.0 1.0 0.4 damping tsuji :pre +pair_style granular +pair_coeff 1 1 jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre + +pair_style granular +pair_coeff 1 1 hertz 1000.0 50.0 tangential linear_history 800.0 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall +pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre [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 options. The total computed forces and torques depend on the combination -of models selected for these various modes of motion. +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 options and parameters are entered in the "pair_coeff"_pair_coeff.html command, as described below. +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 -various combinations of particle types, as determined by the "pair_coeff"_pair_coeff.html command. -For {pair_style granular}, coefficients as well as model options can vary between particle types. -As shown in the second example, -type 1- type 1 interactions are based on a Hertzian normal contact model and 2-2 interactions are based on a {dmt} model (see below). +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 command, otherwise an error would result. +{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: -{hooke} : \(k_n\), damping -{hertz} : \(k_n\), damping -{hertz/material} : E, damping, \(\nu\) -{dmt} : E, damping, \(\nu\), cohesion -{jkr} : E, damping, \(\nu\), cohesion :ol +{hooke} : \(k_n\), \(\eta_\{n0\}\) (or \(e\)) +{hertz} : \(k_n\), \(\eta_\{n0\}\) (or \(e\)) +{hertz/material} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\) +{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); -damping is a damping constant or a coefficient of restitution, depending on -the choice of damping model; E is Young's modulus in units of force/length^2; \(\nu\) is Poisson's ratio -and cohesion is a surface energy density, in units of energy/length^2. +\(\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 between two particles {i} and {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\}_j - \mathbf\{r\}_i\) is the vector separating the -two particle centers +\(\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}. @@ -90,7 +100,7 @@ For the {hertz} model, the normal component of force is given by: \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}/{distance}^2, or equivalently +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: @@ -101,17 +111,17 @@ For the {hertz/material} model, the force is given by: 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 and shear moduli of the +if the elastic and shear moduli 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 model, +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 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\} @@ -131,11 +141,25 @@ will not experience force until they come into contact \(\delta \geq 0\); as the and (\(\delta < 0\)), they experience a tensile force up to \(3\pi\gamma R\), at which point they lose contact. -In addition to the above options, the normal force is augmented by a damping term. The optional -{damping} keyword to the {pair_coeff} command followed by the model choice determines the form of the damping. -The damping coefficient that was specified for the normal model -settings is used in computing the damping term, as described below. Note this damping parameter -may be interpreted differently depending on the model choice. +:line + +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\}\). + +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} @@ -144,91 +168,414 @@ The options for the damping model currently supported are: If the {damping} keyword is not specified, the {viscoelastic} model is used by default. -For {damping velocity}, the normal damping is simply proportional to the velocity: +For {damping velocity}, the normal damping is simply equal to the user-specified damping +coefficient in the {normal} model: + \begin\{equation\} -F_\{N,damp\} = -\gamma_N\mathbf\{v\}_\{N,rel\} +\eta_n = \eta_\{n0\}\ \end\{equation\} -Here, \(\gamma_N\) is the damping coefficient, in units of {mass}/{time}, -\(\mathbf\{v\}_\{N,rel\} = (\mathbf\{v\}_i - \mathbf\{v\}_j) \cdot \mathbf\{n\}\) -is the component of relative velocity along the direction of the vector \(\mathbf\{n\}\) that connects the centers of -particles {i} and {j}. +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: \begin\{equation\} -F_\{N,damp\} = -\gamma_N a m_\{eff\} \mathbf\{v\}_\{N,rel\} +\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, \(\gamma_N\) is the damping coefficient, in units of 1/({time}*{distance}). +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 term is given by: +damping coefficient specified as part of the normal model is intepreted +as a restitution coefficient \(e\). The damping constant \(\eta_n\) is given by: + \begin\{equation\} -F_\{N,damp\} = -\alpha (m_\{eff\}k_N)^\{l/2\} \mathbf\{v\}_\{N,rel\} +\eta_n = \alpha (m_\{eff\}k_n)^\{1/2\} \end\{equation\} -For normal contact models based on material parameters, \(k_N = 4/3Ea\). +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} -For further details, see "(Tsuji et al)"_#Tsuji1992. +\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 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\} +\end\{equation\} :line -Following the normal contact model settings, the {pair_coeff} command 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} : \(\gamma_t\), \(\mu_s\) -{linear_history} : \(k_t\), \(\gamma_t\), \(\mu_s\) :ol +{linear_nohistory} : \(x_\{\gamma,t\}\), \(\mu_s\) +{linear_history} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\) :ol -Here, \(\gamma_t\) is the tangential damping coefficient, \(\mu_s\) is the tangential (or sliding) friction -coefficient, and \(k_t\) is the tangential stiffness coefficient. +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 {linear_nohistory}, a simple velocity-dependent Coulomb friction criterion is used, which reproduces the behavior +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_s \|\mathbf\{F\}_n\|, \gamma_t m_\{eff\}\|\mathbf\{v\}_\{t, rel\}\|) \mathbf\{t\} +\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|\mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\} \end\{equation\} -Where \(\|\mathbf\{F\}_n\) is the magnitude of the normal force, -\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf\{n\}\) is the relative tangential -velocity at the point of contact, \(\mathbf\{v\}_\{t\} = \mathbf\{v\}_R - \mathbf\{v\}_R\cdot\mathbf\{n\}\), -\(\mathbf\{v\}_R = \mathbf\{v\}_i - \mathbf\{v\}_j. +The tangential damping force \(\mathbf\{F\}_\mathrm\{t,damp\}\) is given by: -For {linear_history}, the total tangential displacement \(\mathbf\{\xi\}\) is accumulated during the entire +\begin\{equation\} +\mathbf\{F\}_\mathrm\{t,damp\} = -\eta_t \mathbf\{v\}_\{t,rel\} +\end\{equation\} + +The tangetial 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 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: + +\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"_#. +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}. + +For {tangential linear_history}, the tangential force is given by: + +\begin\{equation\} +\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: \begin\{equation\} \mathbf\{\xi\} = \int_\{t0\}^t \mathbf\{v\}_\{t,rel\}(\tau) \mathrm\{d\}\tau \end\{equation\} -The tangential displacement must in the frame of reference of the contacting pair of particles, - -\begin\{equation\} -\mathbf\{\xi\} = \mathbf\{\xi'\} - \mathbf\{n\}(\mathbf\{n\} \cdot \mathbf\{\xi'\}) -\end\{equation\} - -\noindent Since the goal here is a `rotation', the equation above should be accompanied by a rescaling, so that at each step, -the displacement is first rotated into the current frame of reference $\mathbf\{\xi\}$, then updated: +This accumlated 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 folllows 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. -a simple velocity-dependent Coulomb friction criterion is used, which reproduces the behavior -of the {pair gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by: +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\} - :link(Brill1996) +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 to induce a torque on each +particle according to: + +\begin\{equation\} +\mathbf\{\tau\}_i = -R_i \mathbf\{n\} \times \mathbf\{F\}_t +\end\{equation\} + +\begin\{equation\} +\mathbf\{\tau\}_j = -R_j \mathbf\{n\} \times \mathbf\{F\}_t +\end\{equation\} + +:line + +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 referenece 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: + +\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: +\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. + +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\} +\end\{equation\} + +\begin\{equation\} +\tau_\{roll,j\} = -\tau_\{roll,i\} +\end\{equation\} + +: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: + +{none} +{sds} : \(k_\{twist\}\), \(\gamma_\{twist\}\), \(\mu_\{twist\}\) +{marshall} :ol + +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: + +\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: + +\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: + +\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_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 +\end\{equation\} + +\begin\{equation\} +\eta_\{twist\} = 0.5\eta_ta^2 +\end\{equation\} + +\begin\{equation\} +\mu_\{twist\} = \frac\{2\}\{3\}a\mu_t +\end\{equation\} + +Finally, the twisting torque on each particle is given by: + +\begin\{equation\} +\mathbf\{\tau\}_\{twist,i\} = \tau_\{twist\}\mathbf\{n\} +\end\{equation\} + +\begin\{equation\} +\mathbf\{\tau\}_\{twist,j\} = -\mathbf\{\tau\}_\{twist,i\} +\end\{equation\} + +:line + +LAMMPS automatically sets pairwise cutoff values for {pair_style granular} based on particle radii (and in the case +of {jkr} pulloff 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 + +Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are +functionally the same as the corresponding style without the suffix. +They have been optimized to run faster, depending on your available +hardware, as discussed on the "Speed packages"_Speed_packages.html doc +page. The accelerated styles take the same arguments and should +produce the same results, except for round-off and precision issues. + +These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, +USER-OMP and OPT packages, respectively. They are only enabled if +LAMMPS was built with those packages. See the "Build +package"_Build_package.html doc page for more info. + +You can specify the accelerated styles explicitly in your input script +by including their suffix, or you can use the "-suffix command-line +switch"_Run_options.html when you invoke LAMMPS, or you can use the +"suffix"_suffix.html command in your input script. + +See the "Speed packages"_Speed_packages.html doc page for more +instructions on how to use the accelerated styles effectively. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +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 explictly 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 explictly 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\} +\end\{equation\} + +or + +\begin\{equation\} +E_\{eff,ij\} = \frac\{E_\{ij\}\}\{2(1-\nu_\{ij\})\} +\end\{equation\} + +These pair styles write their information to "binary restart +files"_restart.html, so a pair_style command does not need to be +specified in an input script that reads a restart file. + +These pair styles can only be used via the {pair} keyword of the +"run_style respa"_run_style.html command. They do not support the +{inner}, {middle}, {outer} keywords. + +The single() function of these pair styles returns 0.0 for the energy +of a pairwise interaction, since energy is not conserved in these +dissipative potentials. It also returns only the normal component of +the pairwise interaction force. However, the single() function also +calculates 10 extra pairwise quantities. The first 3 are the +components of the tangential force between particles I and J, acting +on particle I. The 4th is the magnitude of this tangential force. +The next 3 (5-7) are the components of the relative velocity in the +normal direction (along the line joining the 2 sphere centers). The +last 3 (8-10) the components of the relative velocity in the +tangential direction. + +These extra quantities can be accessed by the "compute +pair/local"_compute_pair_local.html command, as {p1}, {p2}, ..., +{p10}. + +:line + +[Restrictions:] + +All the granular pair styles are part of the GRANULAR package. It is +only enabled if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +These pair styles require that atoms store torque and angular velocity +(omega) as defined by the "atom_style"_atom_style.html. They also +require a per-particle radius is stored. The {sphere} atom style does +all of this. + +This pair style requires you to use the "comm_modify vel +yes"_comm_modify.html command so that velocities are stored by ghost +atoms. + +These pair styles will not restart exactly when using the +"read_restart"_read_restart.html command, though they should provide +statistically similar results. This is because the forces they +compute depend on atom velocities. See the +"read_restart"_read_restart.html command for more details. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] + +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. @@ -236,5 +583,33 @@ Model for collisions in granular gases. Physical review E, 53(5), 5382. [(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. - \ No newline at end of file + :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. + + :link(Luding2008) + [(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. + + :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. + + :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. + + :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. + + :link(Thornton1991) + [(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the presence of adhesion. + J. Phys. D: Appl. Phys. 24 1942 \ No newline at end of file diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 3713b9251c..5631240fea 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -52,8 +52,8 @@ using namespace MathConst; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, VISCOELASTIC, TSUJI}; enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN}; -enum {TWIST_NONE, TWIST_NOHISTORY, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_NOHISTORY, ROLL_SDS}; +enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; +enum {ROLL_NONE, ROLL_SDS}; /* ---------------------------------------------------------------------- */ @@ -136,9 +136,9 @@ void PairGranular::compute(int eflag, int vflag) double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; - double knfac, damp_normal; + double knfac, damp_normal, damp_normal_prefactor; double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; + double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; double fs, fs1, fs2, fs3; double mi,mj,meff,damp,ccel,tor1,tor2,tor3; @@ -324,10 +324,12 @@ void PairGranular::compute(int eflag, int vflag) } else{ knfac = E; //Hooke - Fne = knfac*delta; a = sqrt(dR); - if (normal_model[itype][jtype] != HOOKE) + if (normal_model[itype][jtype] != HOOKE){ Fne *= a; + knfac *= a; + } + Fne = knfac*delta; if (normal_model[itype][jtype] == DMT) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } @@ -343,7 +345,8 @@ void PairGranular::compute(int eflag, int vflag) damp_normal = sqrt(meff*knfac); } - Fdamp = -normal_coeffs[itype][jtype][1]*damp_normal*vnnr; + damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; @@ -376,17 +379,21 @@ void PairGranular::compute(int eflag, int vflag) if (normal_model[itype][jtype] == JKR){ F_pulloff = 3*M_PI*coh*Reff; - Fcrit = fabs(Fne + 2*F_pulloff); + Fncrit = fabs(Fne + 2*F_pulloff); + } + else if (normal_model[itype][jtype] == DMT){ + F_pulloff = 4*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); } else{ - Fcrit = fabs(Fne); + Fncrit = fabs(Fntot); } //------------------------------ //Tangential forces //------------------------------ k_tangential = tangential_coeffs[itype][jtype][0]; - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; + damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; if (tangential_history){ shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + @@ -419,7 +426,7 @@ void PairGranular::compute(int eflag, int vflag) fs3 = -k_tangential*history[2] - damp_tangential*vtr3; // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { if (shrmag != 0.0) { @@ -460,64 +467,53 @@ void PairGranular::compute(int eflag, int vflag) if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; else vrlmaginv = 0.0; - if (roll_history){ - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; - // Rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); + // Rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + 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 - 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; - } - - - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][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; - - // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; - - 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; + if (historyupdate){ + 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; } - else{ // - fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; - if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; - else fr = 0.0; - fr1 = -fr*vrl1; - fr2 = -fr*vrl2; - fr3 = -fr*vrl3; + + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][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; + + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; + + 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; } } @@ -527,30 +523,24 @@ void PairGranular::compute(int eflag, int vflag) if (twist_model[itype][jtype] != TWIST_NONE){ magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) if (twist_model[itype][jtype] == TWIST_MARSHALL){ - k_twist = 0.5*k_tangential*a*a;; //eq 32 + 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; + mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2]; } else{ k_twist = twist_coeffs[itype][jtype][0]; damp_twist = twist_coeffs[itype][jtype][1]; mu_twist = twist_coeffs[itype][jtype][2]; } - if (twist_model[itype][jtype] > 1){ - if (historyupdate){ - 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 = TWOTHIRDS*a*Fscrit;//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 - } + if (historyupdate){ + history[twist_history_index] += magtwist*dt; } - else{ - if (magtwist > 0) magtortwist = -damp_twist*magtwist; - else magtortwist = 0; + 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 @@ -747,7 +737,7 @@ void PairGranular::coeff(int narg, char **arg) normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } - else if (strcmp(arg[iarg], "damp") == 0){ + 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 (strcmp(arg[iarg+1], "velocity") == 0){ damping_model_one = VELOCITY; @@ -757,10 +747,12 @@ void PairGranular::coeff(int narg, char **arg) damping_model_one = VISCOELASTIC; iarg += 1; } - else if (strcmp(arg[iarg], "tsuji") == 0){ + 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"); + iarg += 1; } else if (strcmp(arg[iarg], "tangential") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); @@ -785,23 +777,18 @@ void PairGranular::coeff(int narg, char **arg) roll_model_one = ROLL_NONE; iarg += 2; } - else{ + 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 (strcmp(arg[iarg+1], "nohistory") == 0){ - roll_model_one = ROLL_NOHISTORY; - } - else if (strcmp(arg[iarg+1], "sds") == 0){ - roll_model_one = ROLL_SDS; - roll_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized"); - } + 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. iarg += 5; } + else{ + 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"); @@ -814,22 +801,17 @@ void PairGranular::coeff(int narg, char **arg) twist_history = 1; iarg += 2; } - else{ + 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 (strcmp(arg[iarg+1], "nohistory") == 0){ - twist_model_one = TWIST_NOHISTORY; - } - else if (strcmp(arg[iarg+1], "sds") == 0){ twist_model_one = TWIST_SDS; twist_history = 1; - } - else{ + 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. + iarg += 5; + } + else{ error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized"); - } - 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. - iarg += 5; } } else if (strcmp(arg[iarg], "cutoff") == 0){ @@ -1063,7 +1045,6 @@ double PairGranular::init_one(int i, int j) 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]); - 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]); @@ -1082,44 +1063,44 @@ double PairGranular::init_one(int i, int j) // 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_global < 0){ - if (cutoff_type[i][j] < 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 - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - if (normal_model[i][j] == JKR){ - pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); - cutoff += pulloff; - } - else{ - pulloff = 0; - } - if (normal_model[i][j] == JKR) - pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); - cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+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 + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + if (normal_model[i][j] == JKR){ + pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); + cutoff += pulloff; + } + else{ + pulloff = 0; + } - 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) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); - } - cutoff = cutmax; - } + if (normal_model[i][j] == JKR) + pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); + cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+pulloff); + + 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{ - cutoff = cutoff_type[i][j]; + 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]); + cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + } + cutoff = cutmax; } } - else{ + else if (cutoff_type[i][j] > 0){ + cutoff = cutoff_type[i][j]; + } + else if (cutoff_global > 0){ cutoff = cutoff_global; } + return cutoff; } @@ -1211,9 +1192,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, double mi,mj,meff,damp,ccel,tor1,tor2,tor3; double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double knfac, damp_normal; + double knfac, damp_normal, damp_normal_prefactor; double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; + double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; double fs, fs1, fs2, fs3; //For JKR @@ -1357,10 +1338,12 @@ double PairGranular::single(int i, int j, int itype, int jtype, } else{ knfac = E; - Fne = knfac*delta; a = sqrt(dR); - if (normal_model[itype][jtype] != HOOKE) + if (normal_model[itype][jtype] != HOOKE){ Fne *= a; + knfac *= a; + } + Fne = knfac*delta; if (normal_model[itype][jtype] == DMT) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } @@ -1376,7 +1359,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); } - Fdamp = -damp_normal*vnnr; + damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; @@ -1416,20 +1400,21 @@ double PairGranular::single(int i, int j, int itype, int jtype, if (normal_model[itype][jtype] == JKR){ F_pulloff = 3*M_PI*coh*Reff; - Fcrit = fabs(Fne + 2*F_pulloff); + Fncrit = fabs(Fne + 2*F_pulloff); + } + else if (normal_model[itype][jtype] == DMT){ + F_pulloff = 4*M_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); } else{ - Fcrit = fabs(Fne); + Fncrit = fabs(Fntot); } //------------------------------ //Tangential forces //------------------------------ k_tangential = tangential_coeffs[itype][jtype][0]; - if (normal_model[itype][jtype] != HOOKE){ - k_tangential *= a; - } - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; + damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; if (tangential_history){ shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + @@ -1441,7 +1426,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs3 = -k_tangential*history[2] - damp_tangential*vtr3; // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { if (shrmag != 0.0) { @@ -1482,47 +1467,38 @@ double PairGranular::single(int i, int j, int itype, int jtype, if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; else vrlmaginv = 0.0; - if (roll_history){ - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; - // Rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); + // Rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][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; + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][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; - // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; - 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; - } - } - else{ // - fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; - if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; - else fr = 0.0; - fr1 = -fr*vrl1; - fr2 = -fr*vrl2; - fr3 = -fr*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; } + } //**************************************** @@ -1533,25 +1509,19 @@ double PairGranular::single(int i, int j, int itype, int jtype, 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; - mu_twist = TWOTHIRDS*a; + mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2];; } else{ k_twist = twist_coeffs[itype][jtype][0]; damp_twist = twist_coeffs[itype][jtype][1]; mu_twist = twist_coeffs[itype][jtype][2]; } - if (twist_history){ - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = TWOTHIRDS*a*Fscrit;//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 - } - } - else{ - if (magtwist > 0) magtortwist = -damp_twist*magtwist; - else magtortwist = 0; + 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 } }