Hertzian contact models

Hertz contact model is non-linear model which takes in account elastic deformation of spherical solids and changing contact radius. Pure Hertzian model only describes normal deformation of spheres. Later extensions of the model add viscosity, shear response including friction and adhesion. The nomenclature usually used in literature on these models is the following (see [Mod13]):

\(R\)	effective radius
\(E\)	equivalent elastic modulus
\(K\)	effective elastic modulus (\(K=\frac{3}{4}E\))
\(m^{*}\)	effective mass
\(\delta\)	normal geometrical overlap of spheres; \(\delta=r_1+r_2-d\) for two spheres (positive when there is overlap), numerically corresponds to woo.dem.L6Geom.uN with the sign inversed.
\(a\)	contact radius, i.e. radius of the contacting area
\(\epsilon\)	coefficient of restitution
\(\eta\)	viscous coefficient (as in \(F=\eta\dot u\))
\(\dot\delta\)	rate of overlap (normal displacement rate)
\(p\)	function fo contact pressure
\(P_{ne}\)	normal elastic force; integral of \(p\) over the contact area
\(P_{nv}\)	normal viscous force
\(P_n\)	total normal external force acting on the particle

Equivalent radius is computed as harmonic means of particle’s radii

\[\frac{1}{R}=\frac{1}{R_1}+\frac{1}{R_2}\]

which gives correct result also for a sphere-plane contacts (plane having \(R_2=\infty\)).

Equivalent mass

\[\frac{1}{m^{*}}=\frac{1}{m_1}+\frac{1}{m_2}\]

is used for in models of viscosity; \(m_i\) is taken to be zero if the particles is not subject to contact forces (such as when woo.dem.DemData.blocked is set for massless objects like facets or walls).

Moduli (different authors use either \(E\) or \(K=\frac{3}{4}E\)) are computed as

\[ \begin{align}\begin{aligned}\frac{1}{E}=\frac{1-\nu_1^2}{E_1}+\frac{1-\nu_2^2}{E_2},\\\frac{1}{K}=\frac{4}{3}\left(\frac{1-\nu_1^2}{E_1}+\frac{1-\nu_2^2}{E_2}\right)\end{aligned}\end{align} \]

where \(E_i\) are Young’s moduli of contacting particles (woo.dem.ElastMat.young) and \(\nu_i\) their Poisson’s ratios (woo.dem.Cp2_FrictMat_HertzPhys.poisson).

A good summary of various contact models is given at Wikipedia.

Warning

Literature on contact mechanics mostly gives positive sign to compressive stress. This is opposite than the convension used throughout Woo. Therefore the overlap \(\delta\) is positive when spheres are overlapping (whereas woo.dem.L6Geom.uN is negative) and contact force \(P_{n}\) is positive for compression (repulsion) whereas the normal (\(x\)) component of woo.dem.CPhys.force would be negative. Keep this in mind when reading sources.

Hertz contact

Hertz gives contact pressure as

\[p(r)=p_0\sqrt{1-\left(\frac{r}{a}\right)^2}\]

with \(r\in(0,r)\) and \(p_0\) being the maximum pressure in the middle of the contact area; it follows by integration

\[P_{ne}=\int_0^a p(r) 2\pi r \mathrm{d}\,r=\frac{2}{3}p_0\pi a^2.\]

Using the geometrical relationship for non-adhesive contact

\[a^2=R\delta,\]

we can express substitute to obtain

(5)\[P_{ne}=\underbrace{K\sqrt{R}}_{k_{n0}}\delta^{\frac{3}{2}}=k_{n0}\delta^{\frac{3}{2}}\]

where the \(k_{n0}\) term is often separated (woo.dem.HertzPhys.kn0) as it is constant throughout the contact duration. Secant stiffness of the contact is expressed as

\[k_n=\frac{\partial P_{ne}}{\partial \delta}=\frac{3}{2}k_{n0}\delta^{\frac{1}{2}}.\]

By combining the above, we also obtain:

(6)\[a=\sqrt[3]{\frac{R}{K}P_{ne}}.\]

[AE11] also references Landau & Lifschitz’s analytical solution for collision time of purely elastic collision as

\[\tau_{\rm Hertz}=2.214 \left(\frac{\rho}{E}\right)^{\frac{2}{5}}\frac{r_1+r_2}{v_0^{1/5}}\]

where \(\rho\) is density of particles. There is a regression test verifying that simulated collision gives the same result.

Viscosity

Coefficient of restitution

Observable collisions of particles usually result in some energy dissipation of which measure is coefficient of restitution which can be expressed as the ratio of relative velocity before and after collision, \(\epsilon=v_0/v_1\) where \(v_0\), \(v_1\) are relative velocities before and after collision respectively.

The coefficient must be squared when reasoning in terms of energy with kinetic energy \(E_k\) and \(E_d\) the dissipated energy, since kinertic energy \(E_k=\frac{1}{2}mv^2\) contains the velocity squared:

\[\epsilon=\sqrt{\frac{E_{k0}-E_{k1}}{E_{k0}}}=\sqrt{\frac{E_d}{E_{k_0}}}.\]

Similarly, potential field (gravity) can be used to determine coefficient of restitution based on potential energy (\(E_p=mgh\)) for particle falling from intial height \(h_0\) with zero initial velocity, rebounding from horizontal plane and again reaching zero velocity at \(h_1\):

\[\epsilon=\sqrt{\frac{h_1}{h_0}}.\]

Viscous damping

Viscosity adds the \(P_{nv}\) force linearly related to the current rate of overlap \(\dot\delta\) via the linear term \(\eta\),

\[P_{nv}=\eta\dot\delta .\]

The viscous coefficient \(\eta\) is not straight-forwardly realted to the coefficient of restitution \(\epsilon\), which is an integral measure over the whole collision time.

This problem is treated by [AE11] in detail, which derives analytical relationships for viscous coefficient and coefficient of restitution in Hertzian contact. In order to obtain dissipative force resulting in velocity-independent coefficient of restitution, the force must take the form

(7)\[P_{nv}(\delta)=\alpha(\epsilon)\sqrt{m^* k_{n0}}\delta^{\frac{1}{4}}\dot\delta.\]

(this result was previously known, but \(\alpha(\epsilon)\) was only numerically evaluated) and \(\alpha(\epsilon)\) is shown to be

\[\alpha(\epsilon)=\frac{-\sqrt{5}\ln\epsilon}{\sqrt{\ln^2\epsilon+\pi^2}}.\]

The part which is constant throughout the contact life, \(\alpha(\epsilon)\sqrt{m^*k_{n0}}\) is stored in woo.dem.HertzPhys.alpha_sqrtMK.

The total contact force \(P_n\) is superposition of the elastic response \(P\) (5) and normal viscous force \(P_{nv}\) (7), i.e.

(8)\[P_n=P_{ne}+P_{nv}.\]

Nonphysical attraction

The presence of viscous force can lead to attraction between spheres towards the end of collision (\(\dot\delta<0\)) when \(P_n=P_{ne}+P_{nv}<0\). This effect is non-physical in the absence of adhesion. For this reason, such attractive force is ignored, though it can be changed by setting woo.dem.Law2_L6Geom_HertzPhys_DMT.noAttraction.

This effect is discussed in [AE11] and should be taken in account when measuring coefficients of restitution in simulations (see this discussion).

Implementation note:

Hertz contact model can be used via woo.dem.Law2_L6Geom_HertzPhys_DMT, when woo.dem.Cp2_FrictMat_HertzPhys.gamma (surface energy \(\gamma\)) and woo.dem.Cp2_FrictMat_HertzPhys.alpha (Carpick-Ogletree-Salmeron transition parameter \(\alpha\)) are zero. Appropriate engines are returned from woo.dem.DemField.minimalEngines when model is an instance of woo.models.ContactModelSelector with name equal to "Hertz".

Tip

Report issues or inclarities to github.