# Pellet contact model¶

Note

Development of this model was commercially sponsored.

Pellet contact model is custom-developed model with high plastic dissipation in the normal sense, accumulating plastic dissipation immediately when loaded. Adhesion force is initially zero, though it grows with normal plastic deformation, i.e. is dependent on the history of the contact loading.

Normal plastic dissipation is controlled by a single dimensionless parameter \(\alpha\) (`normPlastCoeff`

, averaged from `woo.dem.PelletMat.normPlastCoeff`

) which determines the measure of deviation from the linear response. Adhesion is controlled by the \(k_a\) parameter (`ka`

, computed from the ratio `woo.dem.PelletMat.kaDivKn`

).

Note

As an extension to the model, there is history-independent adhesion (confinement) and persistent plastic deformation due to rolling (thinning), which are documented below.

The normal yield force takes the form

where \(k_n\) is the `normal stiffness`

at origin (\(\frac{\partial F_n^y}{\partial u_n}(u_n=0)=k_n\)) computed as in Linear (Cundall) contact model, \(d_0\) is the `initial contact length`

and \(u_n\) is the `normal displacement`

. The contact accumulates plastic displacement \(u_n^{\mathrm{pl}}\) (`uNPl`

) which is initially zero.

The normal trial force is evaluated as

and it is used to compute the final normal force \(F_n\) in the following way:

If \(F_n^t<0\) (compression), then

- if \(F_n^t\geq F_n^y\) (elastic regime), \(F_n=F_n^t\).
- otherwise (plastic slip), update \(u_n^{\mathrm{pl}}=u_n-\frac{F_n^y}{k_n}\) and set \(F_n=F_n^y\).

If \(F_n^t>0\) (tension), evaluate the adhesion function

\[h(u_n)=-k_a u_n^{\mathrm{pl}}-4\frac{-k_a}{u_n^{\mathrm{pl}}}\left(u_n-\frac{u_n^{\mathrm{pl}}}{2}\right)^2\]and set \(F_n=\min(F_n^t,h(u_n))\).

The \(h(u_n)\) is constructed as polynomial such that it is zero for \(u_n\in\{u_n^{\mathrm{pl}},0\}\) and reaches its maximum \(-k_a u_n^{\mathrm{pl}}\) for

\[h\left(u_n=\frac{u_n^{\mathrm{pl}}}{2}\right)=-k_a u_n^{\mathrm{pl}}.\]where \(k_a\) is the adhesion modulus introduced above.

This plot shows the deformation-force diagram with two loading-unloading cycles, showing how plastic deformation and adhesion grow.

(Source code, png, hires.png, pdf)

## Energy balance¶

As this model unloads linearly, elastic potential energy is computed as with the linear model, i.e.

Normal plastic dissipation, when normal sliding takes place, is computed using backwards trapezoid integration; the previous yield force is

and the normal plastic dissipation

Tangential plastic dissipation \(\Delta E_{pt}\) is computed the same as for the Linear (Cundall) contact model.

Note

Besides the global `energy tracker`

, the `contact law`

stores the dissipated energy also in `PelletMatState`

, if contacting particles define it. Each particle receives the increment of \(\Delta E_{pn}/2\) (`normPlast`

) and \(\Delta E_{pt}/2\) (`shearPlast`

) . This allows to localize and visualize energy dissipation with the granularity of a single particle – by choosing the appropriate scalar with `woo.gl.Gl1_DemField.matStateIx`

.

## Confinement¶

History-independent adhesion can be added by fake confinement \(\sigma_c\) (`confSigma`

) of which negative values lead to constant adhesion.

If \(\sigma_c\) (`confSigma`

) is non-zero, \(F_n\) is updated as

where \(A\) is the `contact area`

.

To account for the effect of \(\sigma_c\) depending on particle size, \(r_{\mathrm{ref}}\) and \(\beta_c\) (`confRefRad`

and `confExp`

) can be specified, in which case \(F_n\) is updated as

this leads to the confining stress being scaled by \(\dagger\) – increased/decreased for bigger/smaller particles with \(\beta_c>0\), or vice versa, decreased/increased for bigger/smaller particles with \(\beta_c<0\).

## Thinning¶

Larger-scale deformation of spherical pellets is modeled using an ad-hoc algorithm for reducing particle radius during plastic deformation along with rolling (`mass`

and `inertia`

are not changed, as if the `density`

were accordingly increased). When the contact is undergoing plastic deformation (i.e. \(F_n^t<F_n^y\); notice that \(F_n^y<0\) and \(F_n^t<0\) in compression), then the particle radius is updated.

This process is controlled by three parameters:

- \(\theta_t\) (
`thinRate`

) controlling the amount of radius decrease per unit normal plastic deformation (\(\Delta u_N^{\mathrm{pl}}\)) and unit rolling angle (\(\omega_r\Delta t\), where \(\omega_r\) is the tangential part (\(y\) and \(z\)) of`relative angular velocity`

\(\vec{\omega}_r\)); the unit of \(\theta_t\) is therefore \(\mathrm{1/rad}\) (i.e. dimensionless). - \(r_{\min}^{\mathrm{rel}}\) (
`thinRelRMin`

), dimensionless, relative minimum particle radius (the original paticle radius is noted \(r_0\)). - \(\gamma_t\) (
`thinExp`

), exponent for decreasing the thinning rate as the minimum radius is being approached. The valid range is \(\gamma_t\in\langle 0,\infty)\), otherwise, the reduction is not done (equivalent to \(\gamma_t=0\)). Larger values of \(\gamma_t\) make approaching the radius floor \(r_{\min}^{\mathrm{rel}}\) harder.

(Source code, png, hires.png, pdf)

When thinning is active, the radius is updated as follows:

Note

Unlike \(u_n^{\mathrm{pl}}\) which is stored per-contact (`uNPl`

) and is zero-initialized for every new contact, the change of `radius`

is *permanent*. It is possible to recover the original radius in `woo.dem.BoxDeleter`

by setting the `recoverRadius`

flag, which re-computes the radius from mass and density.

### Radius-dependence¶

The \(r_{\min}^{\mathrm{rel}}\) and \(\theta_t\) in the previous can be made effectively size-dependent by setting the \(r_{\rm thinRefRad}\) (`thinRefRad`

) to a positive value. In that case, two additional exponents `thinMinExp`

and `thinRateExp`

are used to multiply the equations \((*)\) and \((**)\) by exponential scaling factors. Note that the dependency is on the original radius \(r_0\), **not** the current radius \(r\):

Tip

Ask questions or report issues at github.