Comparative behavior of damping terms of viscoelastic contact force models with consideration on relaxation time

Powder Technology 356 (2019) 735–749

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Comparative behavior of damping terms of viscoelastic contact force models with consideration on relaxation time B. Jian a, G.M. Hu a,⁎, Z.Q. Fang a, H.J. Zhou a, R. Xia a,b,⁎ a b

Department of Mechanical Engineering, Faculty of Engineering, Wuhan University, Wuhan 430072, China Hubei Key Laboratory of Province for Waterjet Theory and New Technology (Wuhan University), Wuhan 430072, China

a r t i c l e

i n f o

Article history: Received 22 May 2019 Received in revised form 26 August 2019 Accepted 31 August 2019 Available online 02 September 2019 Keywords: Contact model Damping Viscoelasticity Coefficient of restitution Discrete element method

a b s t r a c t This paper presents a comparative investigation into behaviors of different damping terms of viscoelastic contact models, and introduces a novel damping term with dimensional reasoning analysis and by considering the effect of relaxation time of viscoelastic spheres. Behaviors of different contact models are compared in terms of predictions of non-dimensional variables, where the influence of damping parameters on behaviors of viscoelastic models is analyzed and the relationship between non-dimensional damping parameters and the coefficient of restitution is discussed. The damping parameter related to relaxation time is found to affect behaviors of the current viscoelastic contact model significantly. Viscoelastic models of different damping terms exhibit different behaviors, especially for a heavily damped case. Comparison between simulation and free fall impact test of viscoelastic balls demonstrates that the current contact model involving relaxation time fits experimental results better than other contact models. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The Discrete Element Method (DEM) initially proposed by Cundall and Strack [1] is a valid numerical tool to study the dynamics of multi-particle system. It has been applied extensively to a wide range of areas including the modelling of particle flows [2] and fluidized beds [3,4], powder sintering [5] and compaction [6], and multiphase flows [7] in comminution [8,9], etc. In DEM, a set of particles is regarded as a multi-body assembly, where the deformation of a particle is described as the overlap at the contact point. The instantaneous motion of each particle is concentrated on its mass center and governed by the Newton's second law. The explicit integration of motion equations is performed step by step in time to give rise to all dynamic data of each particle like position, velocity, acceleration and orientation, etc. These dynamic data are completely determined by the contact forces exerted on particles. Therefore, the accurate evaluation of contact forces is of great importance to the evolution of each particle's motion and consequently to the simulation of the whole system. Numerous contact models applied in DEM applications can be found in [10–13] and the references therein. Normal contact ⁎ Corresponding authors at: Department of Mechanical Engineering, Faculty of Engineering, Wuhan University, Wuhan 430072, China. E-mail addresses: [email protected] (G.M. Hu), [email protected] (R. Xia).

https://doi.org/10.1016/j.powtec.2019.08.110 0032-5910/© 2019 Elsevier B.V. All rights reserved.

laws of elastic bodies were pioneered by Hertz [14], and the Hertz theory has been applied to resolving elastic contact problems. However, this elastic theory is inadequate for inelastic contacts which have energy dissipation. In inelastic collisions, the dissipation mechanism is quite complicated, and the two dissipation mechanisms, i.e. viscoelasticity and plasticity, are generally considered in literatures regardless of the vibration of bodies. The energy dissipation of viscoelastic contacts is caused by viscosity or internal friction [15], while the energy loss of plastic collision is mainly caused by the permanently deformations. Walton and Braun [16] considered elasto-perfect plasticity and proposed a normal contact force model with two different linear springs respectively for loading and unloading stages. Thornton [17] and Vu-Quoc et al. [18,19] proposed more accurate and complex elastic–plastic contact models. By comparison, viscoelastic contact models are more prevalent in DEM applications, owing to the plastic deformation being much smaller than the elastic deformation for most particles, such as sand and glass beads. Viscoelastic models are commonly represented with a spring and a dashpot in parallel. The spring is used to model the elasticity and account for the elastic force, while the dashpot is to model the energy dissipation and account for the dissipative force. The simplest viscoelastic model consists of a linear spring and a linear dashpot. Because of its simplicity, the linear model is usually used in early DEM studies. Nevertheless, this linear model has some drawbacks, such as non-zero contact forces at both the

736

B. Jian et al. / Powder Technology 356 (2019) 735–749

beginning and the end of collision [20,21], the artificial attractive behavior [22,23], and the independence of both the contact duration and the coefficient of restitution on the impact velocity [24]. The prevailing nonlinear contact models in DEM commonly have a nonlinear Hertzian spring based on the Hertz theory, and they are called Hertz-type models [12]. Hertz-type contact models mainly differ in their damping terms, and can be classified with the exponent over the inter-particle displacement of the damping term [13,25]. Lee and Herrmann [26] combined the Hertzian spring with a linear dashpot having the displacement exponent of zero to calculate the viscoelastic contact force of spheres, where there are still several unphysical behaviors as the linear contact model [12]. Tsuji et al. [2] heuristically obtained a damping term with the displacement exponent of 1/4, where the coefficient of restitution can be calculated analytically. Lichtensteiger [27] assumed an infinite number of tiny linear dampers distributed over the contact zone to dissipate the energy, and suggested a nonlinear damping model with the displacement exponent of one. Hunt and Grossley [20] focused on the coefficient of restitution interpreted as the damping in vibroimpact and proposed a damping term with the displacement exponent taking 3/2 for spheres, which avoids the physically incorrect attraction [28]. Besides, Alizadeh et al. [29] introduced a non-Newtonian liquid concept to evaluate the dissipative force, and developed a nonlinear damping term having the variable exponent over the displacement rate and the fixed exponent of one over the displacement. It is worth noting that damping parameters of these previous viscoelastic contact models do not always possess clear physical meanings. From continuum mechanics theory, some researchers tried to derive the dissipative force of viscoelastic spheres analytically and to correlate damping parameters with material properties. Kuwabara and Kono [30] extended the Hertz theory to viscoelastic contacts of particles, and obtained a refined nonlinear damping related to the shear and bulk viscosity of particles. With a quasistatic approximation, Brilliantov et al. [31] explored an analogous formulation of the dissipative force from continuum mechanics. Using an identical assumption, Zheng et al. [32] derived an expression of the dissipative force to calculate the contact force of viscoelastic spheres, and extended it to investigate the contact between ellipsoidal particles in Zheng et al. [33]. Through a perturbation approach, Brilliantov et al. [34] discussed the form of the dissipative force in a more rigorous way. However, those expressions of the dissipative force mentioned above neglect the effect of relaxation time, which should be essential to viscoelastic contact problems. A recent study [35] tried to include the effect of relaxation time into the contact force model of viscoelastic spheres. Apart from the model by Hunt and Crossley [20], viscoelastic contact models mentioned above all suffer from an un-physical attractive force towards the end of contact. To overcome this defect, Zhang and Whiten [21] and Schwager and Pöschel [22] suggested that the collision should be terminated when a zero force instead of a zero overlap is reached, though such a modification would increase the actual coefficient of restitution and reduce the contact time [21,22,36]. Kumar et al. [37] investigated the influence of different spring term on predictions from DEM simulations of spherocylindrical particles. As a core component of viscoelastic contact models, the damping term should be of critical importance to DEM simulation. Different damping terms result in different behaviors of models, such as the dissipation properties and mechanical behaviors, consequently leading to the differentiations in computational dynamics of particle systems. Therefore, it is meaningful to better model the dissipation properties of viscoelastic particles to improve the damping term from both theoretical and practical points of view. Considering the difficulty in measuring the physical dissipative properties of materials directly, a classical concept of the coefficient of restitution is extensively used to

characterize the dissipation of objects. The coefficient of restitution is relatively easily obtained from impact experiments and can be used to determine or adjust the damping parameters of models. This paper aims to investigate and compare the behaviors of different damping terms of viscoelastic contact force models, and attempts to improve the damping term by including the effect of relaxation time of viscoelastic spheres. We try to correlate the coefficient of restitution with the non-dimensional damping parameters. Experimental results of viscoelastic balls are used as a benchmark to compare and assess these diverse viscoelastic contact models. The paper is organized as follows. In Section 2, we apply viscoelastic contact models to addressing the collision of spheres with the motion equation transformed into the non-dimensional form, and review common viscoelastic models of different damping terms. In Section 3, the form of the dissipative force of viscoelastic spheres is explored through dimensional reasoning analysis, and then a novel damping term involving the effect of relaxation time is introduced in connection with the previous study [35]. In Section 4, in terms of predictions of nondimensional variables, such as the non-dimensional displacement, displacement rate and contact force, viscoelastic contact models of different damping terms are compared, and the influences of damping parameters on the behavior of contact models are analyzed, then the relationships between the coefficient of restitution and damping parameters are discussed, finally the comparison of the models against the experiment results of viscoelastic balls by Lichtensteiger et al. [38] is also conducted.

2. Viscoelastic contact model 2.1. Normal impact of viscoelastic particles As shown in Fig. 1, we consider a normal impact between two viscoelastic particles of radii R1, R2 and mass m1, m2. The two particles begin to touch each other at the time t = 0 with a relative velocity vi, and produce a normal contact force F(t) and the contact area of radius a(t) at the contact point. During contact, the displacements of two particle centers are δ1(t) and δ2 (t) respectively, and accordingly their relative displacement is δ(t) = δ1(t) _ + δ2(t), and their displacement rate is δðtÞ. Assuming the mass center and the sphere center of each particle remains coincident during the entire period of contact, the motion

Fig. 1. Spring dashpot model for normal impact between two viscoelastic spheres.

B. Jian et al. / Powder Technology 356 (2019) 735–749

2.3. Non-dimensional form of motion equation governed by Hertz-type contact model

equation of impact particles according to Newton's second law is δðt Þ ¼ −F ðt Þ m €

ð1Þ

where m⁎ is the effective mass with m⁎ = m1 m2/(m1 + m2); €δðtÞrepresents the relative acceleration of two particle centers. The initial condi_ tions for Eq. (1) are: δ(0) = 0, δð0Þ ¼ vi . 2.2. Evaluation of contact force In DEM, the normal contact between two particles is generally modelled by a spring and a dashpot connected in parallel, as shown in Fig. 1. The spring is adopted to represent the elasticity of particles, and the dashpot to dissipate energy in collisions. Then F(t) between two particles consists of two parts F ðt Þ ¼ F e ðt Þ þ F d ðt Þ ¼ kδðt Þw þ cδ_ ðt Þq

ð2Þ

where w is the exponent of displacement in the spring term, and q is the exponent of displacement rate in the damping term; k is the spring stiffness and c is the damping coefficient. The damping coefficient c generally reflects dissipative properties of contacting particles, and commonly correlates with a classical definition of coefficient of restitution e in terms of e = vo / vi, with vo denoting the separating speed of particles at the end of contact. The simplest viscoelastic model is composed of a linear spring (w = 1) and a linear dashpot (q = 1), which was originally used by Cundall and Strack [1] with k and c both adjustable. In such a linear model, the damping term is independent of inter-particle displacement δ(t), and there are some non-physical drawbacks: non-zero contact force at both the initial and the end of contact [20,21], and artificial attraction before the end of the contact [22,23]. To overcome physical inconsistency of the linear model, Hunt and Crossley [20] added δ(t) to the damping term, which is the case for most non-linear viscoelastic contact models that have been developed on the basis of the Hertz theory [14,39]. These nonlinear models based on the Hertz theory commonly have a Hertzian spring to calculate the elastic force Fe(t) of spherical particles, with w = 3/2 and k equal to 4  pffiffiffiffiffi E R 3

kHertz ¼

ð3Þ

where R⁎ is the effective radius of particles and R⁎ = R1 R2/(R1 + R2); E⁎ is the effective Young's modulus given by E⁎ = E1 E2 / (E1 (1−υ22) + E2 (1−υ21)). υ1 and υ2 are Poisson's ratio of particles 1 and 2, E1 and E2 are their Young moduli. As for the damping term, the displacement rate exponent q generally remains as one in most models [11]. Then the damping coefficient c of the damping term for such Hertz-type contact models (w = 3/2, q = 1, k = kHertz) can be generalized into one common equation through dimensional analysis as [12,25]. ð2þ2pÞ=5

c ¼ αmð3−2pÞ=5 vi ð1−4pÞ=5 kHertz

δðt Þp ¼ αmð1−2pÞ kHertz T ð4p−1Þ δðt Þp ð4Þ 2p

with the introduced time parameter T defined by

T¼

m2 2

kHertz vi

737

!1=5

Submitting Eq. (2) into Eq. (1) and making w = 3/2, k = kHertz, and q = 1, the motion equation of particles governed by Hertz-type contact model may be obtained as m €δðt Þ þ kHertz δðt Þ3=2 þ cδ_ ðt Þ ¼ 0

ð6Þ

To proceed non-dimensional analysis, we introduce the nondimensional time ^t with ^t ¼ t T

ð7Þ

On the basis of Eqs. (5) and (7), the non-dimensional motion variables, such as the non-dimensional displacement ^δð^tÞ, displacement ^_ ^tÞ and acceleration €^δð^tÞ, may be defined as rate δð ^δ ¼ δ Tvi

ð8Þ

_ d^δ dδ=ðTvi Þ δ_ ¼ ¼ δ^ ¼ dt=T vi d^t

ð9Þ

_ _ i €δ €^ d^δ dδ=v ¼ ¼ δ¼ dt=T vi =T d^t

ð10Þ

From Eq. (1), we additionally introduce a non-dimensional contact force ^Fð^tÞ with the definition of ^F ¼ F=ðm vi =TÞ ¼ −€δ=ðvi =TÞ, which € € can be related to ^δð^tÞ with ^Fð^tÞ ¼ −^δð^tÞ according to Eq. (10). It is inter-

_ esting to find in Eq. (9) that ^δð^tÞ at the end of contact is connected with e through the equation   _ ^δ_ ^t c ¼ δðt c Þ ¼ −vo ¼ −e vi vi

ð11Þ

By substituting the general form Eq. (4) of c and the new nondimensional variables of Eqs. (7)–(10) into Eq. (6), and after some basic mathematical manipulation, the dimensional motion equation governed by Hertz-type models of different p can be transformed into 2   ^€δ ^t þ ^δ ^t 3=2 kHertz vi T 5 m2

!1=2

 p _   k2Hertz vi 5 þ α ^δ ^t ^δ ^t T m2

!p ¼0

ð12Þ

Using the definition of T of Eq. (5), the above equation may be further reduced to the following non-dimensional form  p _   €^^ ^^3=2 þ α ^δ ^t ^δ ^t ¼ 0 δ t þδ t

ð13Þ

_ with the original initial conditions transformed into: ^δð0Þ ¼ 0, ^δð0Þ ¼ 1. Note that Eq. (13) corresponds to non-dimensional form of the motion equation governed by the Hertz-type model of the damping term with the same p as in Eq. (13), such as p equal to 0, 1/4, 1/2, 1 and 3/2 respectively.

ð5Þ

where the non-dimensional damping constant α depends on the coefficient of restitution e, and the superscript p is the displacement exponent. More definitely the time parameter is related to duration of Hertz contact through tHertz/T = (4/5)3/5 √π Γ(2/5)/Γ(9/10) ≈ 3.2181 [39,40]. Note that Hertz-type models of different damping terms can be classified with the displacement exponent p in Eq. (4).

2.4. Damping terms 2.4.1. Model LH (p = 0) Lee and Hermann [26] combined the linear dashpot (p = 0) and the Hertzian spring (w = 3/2) to calculate the normal contact force between spherical particles. They defined c as the product of m⁎ and an extra parameter γ that controls the amount of dissipation. Making p = 0 in Eq. (4), the damping coefficient of the linear damping term

738

B. Jian et al. / Powder Technology 356 (2019) 735–749

of Model LH may be written as c ¼ γm ¼ αm =T

ð14Þ

2.4.2. Model Ts (p = 1/4) Tsuji et al. [2] suggested a nonlinear damping term with p = 1/4, and the damping coefficient c of the damping term in Model Ts was heuristically defined as c¼α

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m kHertz δðt Þ1=4

ð15Þ

If the moment of displacement backing to zero is regarded as the end of contact, i.e. δ(tc) = 0, there is an analytical relation between α and e given as [28,36]. pffiffiffi h pffiffiffii 5 ln ½e α ¼ −qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; when α∈ 0; 5 2 ln ½e þ π2

ð16Þ

2.4.3. Model KK (p = 1/2) Kuwabara and Kono [30] extended the Hertz theory for viscoelastic contacts of particles. With the dissipation pressure assumed to be proportional to the displacement rate and to be in compliance with the charge distribution arising from the uniform electrostatic potential, they developed an expression of the dissipative force with the exponents p = 1/2 and q = 1, and including a dissipative coefficient D′ related to shear viscosity η and bulk viscosity ξ of materials. A theoretical explanation of the dissipative force was given by Brilliantov et al. [31] independently, revealing analytically the relationship between the dissipative force and viscoelastic material properties like viscosity and Poisson's ratio. However, due to the lack of information on viscosity, the damping coefficient is generally treated as an adjustable parameter that can be determined by the coefficient of restitution. Making p = 1/2 in Eq. (4), the damping coefficient c of the Model KK is rewritten as c¼

2 pffiffiffiffiffi R δðt Þ1=2 ¼ αkHertz Tδðt Þ1=2 D0

ð17Þ

2.4.4. Model Li (p = 1) Lichtensteiger [27] assumed the energy dissipation property of spheres may be modelled with an infinite number of tiny linear dampers over the contact zone, and proposed a nonlinear damping model with p = 1 by introducing a parameter B that is a function of vi and controls the amount of dissipation. After putting p = 1 into Eq. (4), the c of the damping term of Model Li is 2

c ¼ Bδðt Þ ¼ α

kHertz T 3 δðt Þ m

ð18Þ

3

c ¼ λδðt Þ3=2

3. Proposal for damping term with consideration on relaxation time of viscoelasticity 3.1. Damping term related to viscoelasticity Those viscoelastic contact models mentioned in Section 2 generally use a damping term of different exponent p in Eq. (4) to model energy dissipation of viscoelastic contacts. However, the damping term in those models has not always been clearly related to viscoelasticity in solids (like viscosity), except Model KK [30]. Some researchers have focused on the relationship between the dissipative force and the viscoelasticity in solids on the basis of continuum mechanics theory. A derivation of the dissipative force was explored by Brilliantov et al. [31] using an assumption of the quasi-static approximation that neglects the effect of relaxation time in dissipative processes, and recently Brilliantov et al. [34] obtained an expression of the dissipative force through a perturbation approach in a more rigorous way. With an identical quasi-static approximation, Zheng et al. [32] analytically derived a similar expression of the dissipative force for a viscoelastic sphere contacting on a rigid plate. In these studies [30–34], the dependency _ of the dissipative force on the δ(t) and δðtÞ of spheres is the same (p = 1/2 and q = 1), though the dependency on material properties is different. However, it is reasonable that the dissipative force of viscoelastic spheres should reflect the influence of relaxation behavior on energy dissipation from an intuitive reasoning. Namely, the damping term related to viscoelasticity should include the effect of relaxation time besides viscosity, while those viscoelastic contact models mentioned above do not take into consideration relaxation time. 3.2. Dissipative force of viscoelastic particles by dimensional reasoning It is instructive to see how the dissipative force grows during the contact of viscoelastic spheres on the basis of elementary dimensional reasoning. For a linear viscoelastic solid, the stress tensor may be combined of an elastic part and a dissipative part [15,31]. The elastic part is linearly related to the strain and the dissipative part is to the strain rate in a solid. As a component (along the z-axis) of the stress tensor, the contact pressure on the contact zone of viscoelastic bodies is also treated as a combination of an elastic pressure and a dissipative one. Therefore, if pðtÞ is the average contact pressure acting mutually on each solid, it consists of the elastic part pe ðtÞ and the dissipative part pd ðtÞ: pðt Þ ¼ pe ðt Þ þ pd ðt Þ

2.4.5. Model Hunt (p = 3/2) In a study of vibro-impact, Hunt and Grossley [20] introduced a damping factor λ associated with energy loss in an impact, and proposed their damping term with the displacement exponent p equal to the exponent w of the spring term, which fixes its value at 3/2 for spheres. So making p = 3/2 in Eq. (4), the damping coefficient of the Model Hunt can be expressed as k T5 ¼ α Hertz δðt Þ3=2 2 m

equation available for other Hertz-type models mentioned above. For these models, researchers tried to build approximate equations between α and e by some compromise method. For example, some studies have obtained such formula connecting α and e for Model Hunt from certain approximate derivations [41,42] or numerical fitting methods [28,43].

ð19Þ

Note that there is an analytical Eq. (16) that can be used to determine α with e for Model Ts, while there is not such an analytical

ð20Þ

With the quasi-static approximation that the displacement field in a deformed viscoelastic solid completely coincides with that in an elastic solid [31–34], the elastic part pe ðtÞ can be determined by the Hertz theory, while pd ðtÞ may be explored by elementary dimensional reasoning. At any time t, the geometric boundary condition for the contact problem in Fig. 1 may be given as [39]. r2 uz1 ðr; t Þ þ uz2 ðr; t Þ ¼ δðt Þ−  2R

ð21Þ

where uz1 and uz2 are the normal displacements of the surface points on the contact area of the bodies 1 and 2; r denotes the radius coordinate of the point with respect to the center of the contact area.

B. Jian et al. / Powder Technology 356 (2019) 735–749

Provided that the ‘deformation’ within the contact zone is small at any time t, as in [39], we characterize the average state of strain εðtÞ in each solid as the ratio d(t)/a(t) with d(t) = uz(0, t) − uz(a(t), t). Then we may write the average strains ε1 ðtÞ and ε2 ðtÞ of the two contacting bodies 1 and 2 with the following non-dimensional equation

uz1 ð0; t Þ þ uz2 ð0; t Þ−½uz1 ðaðt Þ; t Þ þ uz2 ðaðt Þ; t Þ ¼ aðt Þ

ð22Þ

Making r = 0 and r = a(t) in Eq. (21) respectively, and then putting the two resulting equations into Eq. (23) yields

ε1 ðt Þ þ ε2 ðt Þ ¼

aðt Þ

¼

aðt Þ 2R

ð23Þ

Take the derivative of Eq. (23) with respect to t to obtain the sum of average strain rates ε_ 1 ðtÞ and ε_ 2 ðtÞ, a_ ðt Þ ε_ 1 ðt Þ þ ε_ 2 ðt Þ ¼  2R

ð24Þ

It is reasonable for linear viscoelastic solids to assume the average strain εðtÞ in each solid is proportional to the elastic part pe ðtÞ and the _ average strain rate εðtÞ to the dissipative part p ðtÞ. Then Eq. (24) bed

comes a_ ðt Þ pd ðt Þ∝ε_ 1 ðt Þ þ ε_ 2 ðt Þ ¼  2R

length F d ðtÞ ¼ 2aðtÞpd ðtÞ, then from Eq. (25) F d ðt Þ∝

aðt Þ a_ ðt Þ R

ð25Þ

Eq. (25) indicates that the dissipative part of the contact pressure increases in proportion to the time derivative of the dimension of contact area. If Fig. 2 represents two-dimensional contact of cylindrical bodies with a long strip of 2a(t), the dissipative force per unit axial

F d ðt Þ∝

πaðt Þ2 a_ ðt Þ 2R

ð27Þ

We assume that the relation between a(t) and δ(t) for elastic spheres, i.e. a(t)2 = R⁎δ(t), is still valid for viscoelastic cases, as adopted by Olsson and Jelagin [44]. Taking the time derivative of the equation to _ _ as _ obtain aðtÞ ¼ R δðtÞ=2aðtÞ, Eq. (27) may be expressed by δ(t) and δðtÞ F d ðt Þ∝aðt Þδ_ ðt Þ

ð28Þ

If a dashpot is used to dissipate energy of viscoelastic spheres, the displacement rate exponent q should be one in Eq. (2), and the damping coefficient should be proportional to a(t): pffiffiffiffiffi c ¼ Aðη; ξ; υ; τ; ⋯Þaðt Þ ¼ Aðη; ξ; υ; τ; ⋯Þ R δðt Þ1=2

ð29Þ

with the proportional coefficient A(η,ξ,υ,τ, …) reflecting dissipative properties of viscoelastic solids generally. A(η,ξ,υ,τ, …) should be as a function of material properties, such as shear viscosity η and bulk viscosity ξ, Poisson's ratio υ and relaxation time τ, etc. Dimensional reasoning demonstrates that the displacement exponent p and the displacement rate exponent q should be 1/2 and 1 respectively in the damping term, which is consistent with the previous analytical derivations of the dissipative force from continuum mechanics in literatures [30–34]. In spite of dimensional reasoning analysis presented here being less mathematically rigorous than those previous approaches [30–34] based on continuum mechanics theory, this method is relatively simple without complicated derivations and does not require the same material of contacting spheres. The dimensional reasoning method is expected to be another feasible way to explore

Fig. 2. Comparison of non-dimensional variables predicted by models of different damping terms with α = 0.3 ( ^ ^t, (b) δ– ^_ ^t, (c) ^F–^t, (d) ^F–δ. ^ Li with p = 1, Hunt with p=3/2): (a) δ–

model with β = 0.5,

ð26Þ

For spheres under our consideration, Fd(t) equals to the product of pd ðtÞ and the contact area πa(t)2. Hence from Eq. (25),

d1 ðt Þ d2 ðt Þ þ ε1 ðt Þ þ ε2 ðt Þ ¼ aðt Þ aðt Þ

h i δðt Þ− δðt Þ−aðt Þ2 =2R

739

LH with p = 0,

Ts with p=1/4,

KK with p = 1/2,

Current

740

B. Jian et al. / Powder Technology 356 (2019) 735–749

an appropriate damping term of viscoelastic contact models, for particles with spherical or non-spherical geometry like cylindrical and flatsurface bodies. 3.3. Current damping term proposed to include relaxation time (p = 1/2, q = 1) Taking into account the effect of relaxation time τ, a previous study [35] presented a normal contact force approach for viscoelastic spheres of the same material on the basis of continuum mechanics. In that approach, the viscoelastic contact force was decomposed into an elastic force and a dissipative force, where the elastic force is still calculated by the Hertz theory, and the dissipative force has an expression definitely related to viscosity, relaxation time and Poisson's ratio. The dissipative force for spheres of the standard linear solid is [35]. pffiffiffiffiffi 2η  F d ðt Þ ¼ 1− exp−t=τ R δðt Þ1=2 δ_ ðt Þ 1−υ

ð31Þ −t/τ

where the proportional coefficient A(η,ξ,υ,τ, …) = Avisco (1−exp ), and Avisco is a dissipative coefficient concerned with viscosity and Poisson's ratio of viscoelastic materials, and the term (1−exp−t/τ) is included to concern the effect of relaxation time. For contacting spheres of the same material, we have Avisco = 2η/(1−υ) according to Eqs. (30) and (31), while we cannot directly obtain an exact relationship between the damping coefficient and material properties for spheres of different materials, due to Eq. (30) being no longer applicable. In this case, we may extend Avisco and τ to a concept of effective dissipative coefficient and relaxation time which can reflect generally dissipation behaviors of different viscoelastic materials during contact duration, analogous to the effective Young's modulus E⁎ of particles of different materials for the kHertz in Eq. (3). By introducing the concept of effective Avisco and τ, it is reasonable to assume that the equation A(η,ξ,υ,τ, …) = Avisco (1−exp−t/τ) is still valid to reflect generally dissipative properties of viscoelastic spheres of different materials. Then we can generalize the novel damping term with Eq. (31) of c from the case of the same material to the case of different materials with the effective Avisco and τ. Through dimensional analysis similar to Eq. (4), we rewrite Eq. (31) of c for the novel damping term as   c ¼ α 1− exp−t=ðβT Þ kHertz Tδðt Þ1=2

ð32Þ

with the non-dimensional damping constant α defined by pffiffiffiffiffi α ¼ Avisco R =ðkHertz T Þ

ð33Þ

and β defined with the ratio of τ to T, i.e. β ¼ τ=T

   1=2   ^ €^^ ^^3=2 ^δ_ ^t ¼ 0 δ t þδ t þ α 1− exp−t=β ^δ ^t

ð35Þ

ð30Þ

When a dashpot is used to dissipate energy and to calculate the dissipative force of viscoelastic spheres, a novel damping term may be proposed from Eq. (30), considering relaxation time to well reflect dissipative properties of viscoelastic solids. This damping term should have p = 1/2 and q = 1 from Eq. (30), and its damping coefficient c has to be in conformity to Eq. (29) of c based on dimensional reasoning with the equation  pffiffiffiffiffi pffiffiffiffiffi c ¼ Aðη; ξ; υ; τ; ⋯Þ R δðt Þ1=2 ¼ Avisco 1− exp−t=τ R δðt Þ1=2

Eq. (31) and Eq. (32) respectively, and in this case, the current model reduces to the form identical to the Model KK (Eq. (17)). The c by Eq. (31) for the novel damping term links up with viscoelastic material properties, and c may be determined with material parameters via Eq. (31) directly, but the practical use of Eq. (31) is limited by the lack of information on viscosity and relation time. By contrast, the current Eq. (32) of c is expected to be more practical, for the damping constants α and β may be both related to e and the Hertzian spring stiffness kHertz and the time parameter T are both the same as other Hertztype models. By putting the non-dimensional variables of Eqs. (7)–(10) into Eq. (6), the motion equation Eq. (6) with the damping coefficient c of Eq. (32) can be transformed into the following non-dimensional form,

ð34Þ

where the damping constants α and β are both related to the coefficient of restitution e. Compared to the general formula Eq. (4) for damping coefficient, the current Eq. (32) includes the effect of relaxation time on dissipation of viscoelastic spheres with an additional term (1−exp−t/(βT)). It is interesting to note that as relaxation time τ (or β) approaches zero, the term exp−t/τ and exp−t/(βT) may vanish in

with the same initial conditions as the ones of Eq. (13). Eq. (35) corresponds to the non-dimensional form of the motion equation governed by the current contact model with c of Eq. (32). Compared with Eq. (13), the above non-dimensional equation Eq. (35) has an additional ^ term ð1− exp−t=β Þ to involve influence of relaxation time on viscoelastic collision. 4. Results and discussion 4.1. Predictions of non-dimensional variables of contact models With p respectively set as 0, 1/4, 1/2, 1 and 3/2, Eq. (13) respectively corresponds to the non-dimensional motion equations of the Hertztype Models LH, Ts, KK, Li and Hunt. By contrast, Eq. (35) with the parameter β concerning relaxation time is the non-dimensional motion equation governed by the current model. With given values of α (and β), we may numerically solved Eq. (13) (and Eq. (35)) to calculate the _ non-dimensional variables of displacement ^δ, displacement rate ^δ and ^ contact force F of these Hertz-type models. To compare these models of different damping terms and investi_ gate their characteristics of ^δ, ^δ and ^F during the whole process of con_ tact, the graphs of these variables (i.e. ^δ, ^δ and ^F versus ^t, and ^F–^δ loops) were plotted in Fig. 2 with α = 0.3 and β = 0.5. ^ ^t curves predicted by Hertz-type models with Fig. 2(a)shows the δ– α = 0.3 approximate sinusoidal functions of semi-cycle, while their curved amplitudes and periods are different from each other. During the contact, δ^ increases from zero to a maximum (corresponding to the compressing or approaching phase of particles), afterward it reduces (corresponding to the recovering or separating phase of particles) and eventually backs to zero. The model of a larger p predicts a higher ^ ^tÞ (i.e. HuntNLi N KK N Ts N LH) and a slightly shorter amplitude of δð ^ ^tÞ (i.e. HuntbLi b KK b Ts b LH), except for the current period of δð model (p = 1/2 and β = 0.5). The current model involving the effect ^ ^tÞ close to the Model KK (p = of relaxation time has a period of δð ^ ^tÞ is apparently larger than the Model KK's. 1/2), while its amplitude of δð _ In Fig. 2(b), ^δ decreases from one to zero, then its sign reverses to separate the contacting particles, and afterwards its magnitude increases to a maximum which occurs at the moment of ^F returning to _ zero. The final ^δ is bigger than negative one due to dissipation. Figs. 2(c)and (d)indicate Hertz-type models predict ^F starting from zero and ending at zero, apart from Model LH (p = 0) with an initial ^F of α and a final ^F of non-zero value. Kruggel-Emden et al. [12] also observed that Model LH has non-zero forces at the beginning and the end of contact, which is unphysical in reality. This unphysical behavior

B. Jian et al. / Powder Technology 356 (2019) 735–749

741

Fig. 3. ^ δ–^t predicted by models with α of different values: (a) LH with p=0, (b) Ts with p=1/4, (c) KK with p=1/2, (d) Current model with β=1/2, (e) Li with p=1, (f) Hunt with p=3/2.

should result from the linear damping of Model LH as in the linear contact model. Except for the current model, a smaller p predicts a lower peak of ^F (HuntNLi N KK N Ts N LH) and a steeper initial slope of the ^F curve, as shown in Figs. 2(c) and (d). As p b 1, the initial slope of ^F becomes infinite, such as Model KK (p = 1/2) and Model Ts (p = 1/4). Moreover for Model LH (p = 0), there is a step jump of ^F occurring at the start point,

where the ^F curve was assumed to be still changing from zero with an initial slope of infinity and this slope of infinity was such that a step _ jump appears. Therefore, the decline rate of ^δ for a model with a smaller p is faster during the approaching phase, i.e. HuntbLi b KK b Ts b LH in Fig. 2(b). We observed in Figs. 2(c) and (d) that for most models ^F reverses the sign, i.e. becoming attractive, before ^δ backs to zero. The predicted

_ Fig. 4. ^ δ–^t predicted by models with α of different values: (a) LH with p = 0, (b) Ts with p = 1/4, (c) KK with p = 1/2, (d) Current model with β = 1/2, (e) Li with p = 1, (f) Hunt with p = 3/2.

742

B. Jian et al. / Powder Technology 356 (2019) 735–749

attraction in Figs. 2(c) and (d) is alleviated with p increasing until p = 3/2 when the attraction disappears in Model Hunt. The attractive force may be avoided by Model Hunt in Fig. 2(d), for it predicts ^F and ^ δ returning to zero simultaneously [28]. Zhang and Whiten [21] observed the attraction in Model Ts, and Kruggel-Emden et al. [12] observed in the Models LH and KK. They [12,21] pointed that the attractive force has an effect of decelerating contact particles, and cor_ respondingly in Fig. 2(b) the ^δ–^t curve slightly goes ‘upward’ before ^δ backs to zero, which is inconsistent with the physical collision between cohesionless particles. Model Li and the current model both show this unfeasible aspect of attraction in Figs. 2(c) and (d). To avoid the unfeasible attractive force, as suggested by Zhang and Whiten [21], we may define a collision terminates when a zero force reaches, though this treatment may cause a shorter contact time directly [21,22,36] and give rise to a ‘residual’ displacement at the end of contact. The ^F–^ δ hysteretic loop in Fig. 2(d) demonstrates that there is energy dissipation during a contact. The area enclosed by the ^F–^δ loop is related to the coefficient of restitution e. With an increase of p, the hysteretic ^ Apart from loop becomes thinner at lower δ^ and thicker at higher δ. Model Hunt, the treatment of ending contact with zero contact force would omit the part below the δ^ axis, making the area of the hysteretic loop shrink and causing e to rise. 4.2. Influences of damping parameters on behaviors of contact models 4.2.1. Influences of damping constant α on behaviors of contact models To explore the influences of α on behaviors of viscoelastic contact models, we solved Eqs. (13) and (35) (β = 0.5) with several chosen _ values of α, and obtained the variables ^δ, ^δ and ^F versus ^t of each model, which were respectively plotted in Figs. 3–5. As α = 0 with the damping vanishing, the Hertz-type model reduces to an ideal elastic case, and the curves of δ^ and ^F are perfectly

_ symmetrical along the axis of ^t as shown in Figs. 3 and 5, and the final ^δ is equal to negative one in Fig. 4. The ^δ–^t curves in Fig. 3 indicate that the increased α lowers the maximum δ^ and brings forward its occurrence that corresponds to the mo-

_ ment when ^δ decreases to zero in Fig. 4, and also puts off ^δ returning to zero for each model. Correspondingly the symmetry of curves ^δ–^t tends to move to the left. Fig. 4 illustrates that the growing α in each model makes the magni_ tude of the final ^δ (i.e. e via Eq. (11)) smaller, and also makes the decline _ rate of ^δ faster during the approaching phase. As α increases in Fig. 4(a)– _ (e), the ‘upward’ segment of the ^δ–^t curve before the end of contact be-

comes longer and more noticeable, which means the deceleration effect caused by attractive force should be more obvious in Fig. 5, except the Model Hunt in Figs. 4(f) and 5(f). Fig. 5 shows the influence of α on the ^F of contact models. With α increasing from zero, the peak of ^F occurs earlier, and its magnitude is firstly lowered and then raised for each model. This is because the contribution of the damping force on the contact process grows with α. As α takes a small value, the damping effect is insignificant and the spring dominates the contact, then the peak of ^F should decrease with α induc-

ing the decline of the maximum ^δ in Fig. 3. However, as α goes up to a larger value with the damping instead of the spring dominant in the contact, the damping force grows big enough to compensate the reduction of the maximum ^δ, and then makes the peak point of ^F rise up and move to the left. The initial slope of the ^F–^t for each model in Fig. 5 becomes steeper with α, apart from Model LH in Fig. 5(a). For Model LH, the growing α increases the initial ‘jump’ force, and as α is up to a big enough value, such as α = 1.0, the initial ‘jump’ force may become the maximum point during the whole contact. Moreover, we find in Fig. 5(a)–(e) that the growth of α makes the artificial attraction last longer, and consequently the resulting

Fig. 5. ^F–^t predicted by models with α of different values: (a) LH with p=0, (b) Ts with p=1/4, (c) KK with p=1/2, (d) Current model with β=1/2, (e) Li with p=1, (f) Hunt with p=3/2.

B. Jian et al. / Powder Technology 356 (2019) 735–749

Fig. 6. Predictions of the current model with α=0.3 and β of different values (

β=1.0,

deceleration becomes more obvious before the end of contact in Fig. 4 (a)–(e), aside from the Model Hunt in Fig. 5(f). 4.2.2. Influences of damping constant β on behaviors of contact model involving relaxation time To investigate the influence of relaxation time on behaviors of the _ current model, ^ δ, ^ δ and ^F were calculated from Eq. (35) with α = 0.3 and β of several chosen values: β → 0, β = 0.1, β = 0.5 and β = 1.0, _ and the curves of ^ δ–^t, ^ δ–^t, ^F–^t and ^F–^δ were plotted in Fig. 6.

The influence of β on the prediction of the ^δ-^t by the current model is indicated in Fig. 6(a). The growth of β tends to increase the maximum ^δ of the current model, while it seems to have little influence on the moment when the maximum ^δ and zero ^δ occur respectively. The curve

_ δ–^t is also changed with β in Fig. 6(b). As β grows, the decline shape of ^ _ _ rate of ^ δ during approaching phase is slightly slower, and the ^δ–^t before the end of contact is ‘lowered down’, which means that the coefficient of restitution e rises up according to Eq. (11). Fig. 6(c)shows the influence of β on the contact force predicted by the current model with α = 0.3. β increasing from zero makes the initial slope of ^F reduce from the infinite, and raises the peak of ^F with its occurrence delayed. However as β varies, the contact time predicted by the current model is little affected, and the curve shape of ^F in Fig. 6 (c) or Fig. 6(d) before the end of contact is changed insignificantly either. The predicted ^F–^ δ loops for different β coupled with α = 0.3 were plotted in Fig. 6(d). It shows that the curved changes between different β are mainly focused on the loading path and the initial unloading path. The augment of β lowers the initial loading path of ^F, but enlarges the

maximum of both ^ δ and ^F, leading to the ^F–^δ loop becoming thinner and longer in Fig. 6(d).

4.2.3. Relationship of coefficient of restitution and damping parameters As we regard the impact to be terminated when a zero force is reached _ firstly, we can obtain e with δ^ at the end of contact according to Eq. (11). By numerically solving Eq. (13) with a list of α starting from zero, a corresponding list of e can be calculated, and the relationships of e and α for

β=0.5,

β=0.1,

743

_ β→0 identical to the Model KK): (a) ^δ–^t, (b) ^δ–^t, (c) ^F–^t, (d) ^F–^δ.

models (LH, Ts, KK, Li and Hunt) were found, and listed in Table 1 and plotted in Fig. 7(a). As for the current model, the parameter β concerning relaxation time should be additionally specified in Eq. (35) to calculate e. Therefore, with several chosen values of β like 0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 2.5 and 3.2, the relationships between e and α of the current model were also obtained in Table 2 and plotted in Fig. 7(b). Table 1 The coefficient of restitution e calculated with the damping parameter α for the Hertz-type contact models of different p. α

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

e LH (p = 0)

Ts (p = 1/4)

KK (p = 1/2)

Li (p = 1)

Hunt (p = 3/2)

Eq. (14)

Eq. (15)

Eq. (17)

Eq. (18)

Eq. (19)

1.0000 0.9656 0.9332 0.9023 0.8730 0.8449 0.8180 0.7923 0.7677 0.7441 0.7214 0.6207 0.5374 0.4680 0.4097 0.3604 0.3185 0.2827 0.2518 0.2252 0.2022 0.1645 0.1355 0.1128 0.0949 0.0806 0.0690 0.0595 0.0517 0.0452

1.0000 0.9724 0.9459 0.9205 0.8959 0.8723 0.8495 0.8276 0.8065 0.7861 0.7664 0.6778 0.6030 0.5394 0.4849 0.4380 0.3972 0.3617 0.3305 0.3030 0.2787 0.2379 0.2051 0.1784 0.1565 0.1383 0.1230 0.1100 0.0990 0.0894

1.0000 0.9773 0.9552 0.9338 0.9130 0.8929 0.8733 0.8544 0.8360 0.8182 0.8010 0.7221 0.6541 0.5953 0.5440 0.4992 0.4597 0.4249 0.3939 0.3663 0.3416 0.2993 0.2646 0.2358 0.2116 0.1911 0.1736 0.1584 0.1452 0.1337

1.0000 0.9833 0.9671 0.9512 0.9358 0.9208 0.9061 0.8918 0.8779 0.8642 0.8509 0.7890 0.7336 0.6840 0.6394 0.5991 0.5627 0.5297 0.4997 0.4724 0.4474 0.4036 0.3666 0.3349 0.3077 0.2841 0.2635 0.2454 0.2293 0.2150

1.0000 0.9868 0.9740 0.9615 0.9494 0.9375 0.9259 0.9146 0.9035 0.8927 0.8822 0.8329 0.7885 0.7484 0.7120 0.6786 0.6480 0.6197 0.5936 0.5694 0.5469 0.5062 0.4704 0.4388 0.4107 0.3855 0.3629 0.3424 0.3238 0.3069

744

B. Jian et al. / Powder Technology 356 (2019) 735–749

Fig. 7. Relationships between e and α predicted by models: (a) Models without β, (b) Current model with β of different values.

4.3. Comparison of simulation results of contact models against experimental results

from the height of h onto an aluminum plate that was glued to the top of a force transducer. The transducer was then mounted on a steel base placed on a heavy concrete block to measure the impact force over time. The impact force versus time was recorded and used to calculate the acceleration of a test sphere according to the Newton's second law. Ignoring air resistance, the first and the second integration of the acceleration over time yield the displacement rate history and the displacement history of the test sphere respectively. Since the aluminum plate is much more rigid than test spheres and contact duration reported by Lichtensteiger et al. [38] is generally located on the millisecond scale, a test sphere dropping onto the plate may be approximately treated as the sphere in contact with a rigid plate of infinite mass, and the air resistance and the gravity can be neglected with a reasonable tolerance during a contact. Therefore the motion equation Eq. (1) is applicable to the test spheres, and we can apply the Hertztype contact models governing Eq. (6) to the free fall impact.

The comparison against experimental results was implemented to assess the adequacy of these viscoelastic contact models to reality. Lichtensteiger et al. [38] performed a free fall impact test with various spheres of viscoelastic materials, where the test objects were dropped

4.3.1. Comparison of contact models and experiment with uniform impact velocity The results of indoor lacrosse ball (ILB) and outdoor lacrosse ball (OLB) dropped from the reference drop height href = 15 cm

Fig. 7 shows that the coefficient of restitution e is a monotonically decreasing function of α for each model. As α = 0 with the damping vanishing, we have the maximum e = 1 which corresponds to an ideal elastic contact. In Fig. 7(a), the growth of p makes the e–α curve move to the right with a fixed starting point. If α is fixed for Hertz-type models, a larger p predicts a higher e, i.e. Hunt N Li N KK N Ts N LH; reversely if e is fixed, the model of a larger p has a bigger α. For the current model in Fig. 7(b), the parameter β related to relaxation time of materials influences the relationship of e and α. The growing β makes the e–α curve move to the right with a fixed starting point. So as we fix α, a larger β predicts a higher e; in turn if e is fixed, a larger β corresponds to a larger α for the current model.

Table 2 The coefficient of restitution e calculated with the damping parameter α for the current model with Eq. (32) of different β. α

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

e β→0

β = 0.1

β = 0.3

β = 0.5

β = 1.0

β = 1.5

β = 2.0

β = 2.5

β = 3.2

1.0000 0.9773 0.9552 0.9338 0.9130 0.8929 0.8733 0.8544 0.8360 0.8182 0.8010 0.7221 0.6541 0.5953 0.5440 0.4992 0.4597 0.4249 0.3939 0.3663 0.3416 0.2993 0.2646 0.2358 0.2116 0.1911 0.1736 0.1584 0.1452 0.1337

1.0000 0.9778 0.9562 0.9353 0.9150 0.8953 0.8761 0.8576 0.8396 0.8221 0.8052 0.7276 0.6607 0.6026 0.5519 0.5074 0.4683 0.4336 0.4028 0.3753 0.3506 0.3083 0.2735 0.2446 0.2202 0.1995 0.1817 0.1664 0.1530 0.1412

1.0000 0.9796 0.9597 0.9404 0.9215 0.9032 0.8854 0.8681 0.8513 0.8349 0.8190 0.7457 0.6818 0.6259 0.5767 0.5333 0.4949 0.4606 0.4299 0.4024 0.3776 0.3348 0.2993 0.2695 0.2443 0.2227 0.2040 0.1878 0.1736 0.1611

1.0000 0.9812 0.9629 0.9450 0.9276 0.9105 0.8939 0.8778 0.8619 0.8466 0.8316 0.7622 0.7011 0.6472 0.5995 0.5570 0.5191 0.4852 0.4547 0.4271 0.4022 0.3589 0.3228 0.2922 0.2662 0.2438 0.2244 0.2074 0.1925 0.1792

1.0000 0.9845 0.9692 0.9542 0.9395 0.9251 0.9110 0.8972 0.8837 0.8704 0.8574 0.7966 0.7420 0.6930 0.6489 0.6091 0.5731 0.5405 0.5109 0.4838 0.4591 0.4155 0.3785 0.3468 0.3194 0.2956 0.2746 0.2561 0.2397 0.2250

1.0000 0.9869 0.9740 0.9612 0.9487 0.9364 0.9243 0.9123 0.9006 0.8891 0.8779 0.8244 0.7757 0.7313 0.6908 0.6539 0.6201 0.5891 0.5606 0.5344 0.5102 0.4671 0.4299 0.3976 0.3694 0.3444 0.3224 0.3027 0.2850 0.2691

1.0000 0.9887 0.9776 0.9666 0.9558 0.9451 0.9345 0.9241 0.9139 0.9038 0.8939 0.8465 0.8028 0.7625 0.7253 0.6910 0.6594 0.6301 0.6029 0.5777 0.5542 0.5120 0.4751 0.4426 0.4139 0.3883 0.3654 0.3449 0.3263 0.3095

1.0000 0.9902 0.9804 0.9708 0.9612 0.9518 0.9425 0.9333 0.9242 0.9153 0.9064 0.8641 0.8245 0.7878 0.7536 0.7217 0.6920 0.6643 0.6385 0.6143 0.5917 0.5506 0.5142 0.4820 0.4531 0.4272 0.4039 0.3828 0.3636 0.3461

1.0000 0.9917 0.9834 0.9752 0.9670 0.9590 0.9510 0.9431 0.9353 0.9276 0.9200 0.8831 0.8484 0.8157 0.7850 0.7561 0.7290 0.7034 0.6794 0.6567 0.6353 0.5960 0.5607 0.5291 0.5005 0.4745 0.4510 0.4294 0.4097 0.3916

B. Jian et al. / Powder Technology 356 (2019) 735–749

745

Table 3 Test balls extracted from Lichtensteiger et al. [38]. Object

Indoor Lacrosse Ball (ILB) Outdoor Lacrosse Ball (OLB)

Mass

Radius

Incident velocity

Rebound velocity

Coefficient of restitution

Contact time

m [g]

R [mm]

vi [cm/s]

vo [cm/s]

e

tc [ms]

154.7 147.2

31.75 31.75

171.6 (1) 171.6 (1)

79.4 148.2

0.463 0.864

3.83 4.15

Note: (1) Ignoring the air resistance, href = 15 cm corresponds to the incident velocity vi = 171.6 cm/s with the gravity acceleration g = 981 cm/s2.

were extracted from Lichtensteiger etal. [38] and for our comparison. ILB and OLB are of the same size and nearly of the same mass, while they quite differ in damping characteristics. ILB is heavily-damped (e = 0.463), while OLB is relatively lightlydamped (e = 0.864). The test data related to ILB and OLB were summarized in Table 3, and their experimental results of the displacement δ, displacement rate δ_ and contact force F were presented in Figs. 8 and 9 (with the symbol •) respectively for comparison with simulation results. Using the test data of ILB and OLB in Table 3, we estimated the necessary parameters of each contact model for simulation and listed them in Table 4. The effective mass m⁎ and radius R⁎ respectively equaled to the mass and radius of test spheres, and they were given in the second and third columns of Table 4. Because Young's modulus and Poisson's ratio are not available for ILB and OLB in Lichtensteiger [27,38], one cannot obtain kHertz for a Hertz-type contact model via Eq. (3) directly. The kHertz may be determined according to the fact that the contact force was equal to the elastic force when the sphere is at the point of maximum displacement. With the recorded contact force corresponding to the maximum displacement from Lichtensteiger [27,38], we evaluated k Hertz for simulations of ILB and OLB and listed the values of k Hertz in the fourth column of Table 4. With the determined parameters m⁎ and kHertz and the incident velocity vi = 171.6 cm/s corresponding to href = 15 cm, the time parameter T may be evaluated via Eq. (5) for ILB and OLB, as shown in the fifth column of Table 4.

Appropriate values of α for all Hertz-type contact models and β for the current model can be read off from Fig. 7 or from Tables 1 and 2 with a given e. For those Hertz-type models without β, we obtained the values of α for simulations of ILB and OLB from Fig. 7(a) or Table 1 at the condition of e = 0.463 for ILB and e = 0.864 for OLB, and listed them in the seventh column of Table 4. As for the current model, one should first specify β in order to determine α from Fig. 7(b) or Table 2. With e = 0.463 (ILB) and e = 0.864 (OLB), we used β = 0.5 for ILB and β = 2.0 for OLB respectively to obtain the values of α of the current model for simulations and listed them in Table 4. One can specify β on the basis of the fact that the curve shapes of the contact force, displacement rate and displacement versus time of the current model being affected by the value of β besides α, as discussed in Section 4.2. We fitted the proper values of β for simulations, i.e. β = 0.5 (ILB) and β = 2.0 (OLB), from experimental results of the contact force, displacement rate and displacement history. We obtained the relaxation time τ of ILB and OLB by multiplying β and T according to Eq. (34), as shown in the eighth column of Table 4. The damping coefficient c may be determined with its own definition for each Hertz-type contact model, as given in the last column of Table 4. Comparisons of simulation results against the experiment of ILB and OLB were respectively presented in Figs. 8 and 9, where variables of displacement δ, displacement rate δ_ and contact force F versus time t were respectively plotted in the

Fig. 8. Comparison of predictions of different contact models against single-drop test with ILB dropped from the height of 15 cm ( LH with p = 0, _ (c) F–t, (d) F–δ. Current model with β = 0.5, Li with p = 1, Hunt with p = 3/2, • experiment [38]): (a) δ–t, (b) δ–t,

p = 1/2,

Ts with p = 1/4,

KK with

746

B. Jian et al. / Powder Technology 356 (2019) 735–749

Fig. 9. Comparison of predictions of different contact models against single-drop test with OLB dropped from the height of 15 cm ( LH with p=0, _ (c) F–t, (d) F–δ. Current model with β=2.0, Li with p=1, Hunt with p=3/2, • experiment [38]): (a) δ–t, (b) δ–t,

Ts with p=1/4,

KK with

p=1/2,

graphs (a)–(c), and the loops of the contact force F versus displacement δ were in the graph (d). Among these Hertz-type models, the current model can best match the experimental results of ILB in Fig. 8 and OLB in Fig. 9. This is due to that the damping term of the current model involves the effect of relaxation time during a viscoelastic contact. As discussed in Subsection 4.2.2, the involved relaxation time in the current model increases the peak force and the maximum displacement simultaneously, but insignificantly affects the contact time of the current model, compared with Model KK. For those models regardless of relaxation time, there are a few unsuitable predictions of the force history in Fig. 8 (c) for ILB and in Fig. 9 (c) for OLB, and of the displacement history in Figs. 8 (a) and 9 (a) and the displacement rate history in Figs. 8(b) and 9(b). Since the initial slope of the measured force is finite, the model with a smaller p predicts a steeper initial slope of contact force, leading to larger deviations of results of F and δ_ between the simulation and experiment during the initial contact process,

i.e. deviations of Model Hunt b Li b KK b Ts b LH in Fig. 8 (b)– (d) for ILB and Fig. 9 (b)–(d) for OLB. Comparing Fig. 8 with Fig. 9, differences among simulation results by different Hertz-type contact models for heavily-damped ILB are much greater than those for lightly-damped OLB. This is due to that for lightly-damped OLB, the spring dominates the contact and the damping insignificantly affects the predictions, then the differences of damping term result in insignificant differences of simulation results; while for heavily-damped ILB, the damping plays a relatively important role, and the damping term has a notable effect on the predictions, then the differences of damping term result in significant changes of simulation results. Since zero contact force instead of zero displacement is used as the contact ending criteria for the simulation, there should be a ‘residual’ displacement existing at the end of contact, aside from Model Hunt, as shown in Figs. 8 (a), (d) and 9 (a), (d). The ‘residual’ displacement has also been found by Lichtensteiger et al. [38] and Cross [45] independently in impact experiments of viscoelastic objects. Comparing Figs. 8

Table 4 Simulation parameters estimated for each contact model. Object

m⁎ [g]

R⁎ [mm]

kHertz [10 6 N/m3/2]

T [ms]

Model

α

τ [ms]

c [N·s / m3/2]

ILB

154.7

31.75

1.688

1.374

OLB

147.2

31.75

1.724

1.335

LH Ts KK current model Li Hunt LH Ts KK current model Li Hunt

0.508 0.646 0.792 0.973 (β = 0.5) 1.138 1.647 0.087 0.108 0.130 0.263 (β = 2.0) 0.181 0.236

– – – 0.687 – – – – – 2.671 – –

57.256 329.892 δ(t)1/4 1836.91 δ(t)1/2 2255.88 (1 − exp−t/τ) δ(t)1/2 54,348.8 δ(t) 1.62*106 δ(t)3/2 9.535 54.158 δ(t)1/4 299.704 δ(t)1/2 605.257 (1 − exp−t/τ) δ(t)1/2 8701.3 δ(t) 2.374*105 δ(t)3/2

B. Jian et al. / Powder Technology 356 (2019) 735–749

Fig. 10. Experimental and simulated relationships of coefficient of restitution e versus impact velocity vi ( LH with p=0, Ts with p=1/4, KK with p=1/2, Current model, Li with p=1, Hunt with p=3/2, • experiment [38]).

and 9, the heavily-damped ILB tends to have a larger ‘residual’ displacement than the lightly-damped OLB, and this tendency also appears in most Hertz-type models. Besides, growing p of Hertz-type models tends to reduce the amount of the predicted ‘residual’ displacement until p goes up to 3/2 when the ‘residual’ displacement vanishes in Model Hunt. Furthermore, the ‘residual’ displacement of the current model is also affected by relaxation time in Figs. 8 and 9.

4.3.2. Influence of impact velocity As the simulation parameters of viscoelastic contact models, i.e. m⁎, kHertz and c, are fixed, the changing incident velocity vi directly affects T through Eq. (5), and consequently changes the non-dimensional parameters α and β that characterize the motion equations of models.

Fig. 11. Experimental and simulated relationships of contact time tc versus impact velocity vi ( p=1, Hunt with p=3/2, • experiment [38]) for: (a) ILB, (b) OLB.

747

Therefore, the predictions of these contact models should have a dependence on the impact velocity. Using the parameters of m⁎, kHertz and c in Table 4, we simulated the drop impact tests of ILB and OLB by varying the incident velocity vi to investigate its influence on each Hertz-type contact model. The velocity range of the simulation was selected from 30 cm/s to 260 cm/s for the convenience of comparison between simulation and experiment results, as the experimental velocity of Lichtensteiger et al. [38] was in the range of 44.3 cm/s to 242.6 cm/s. The impact velocity dependence of the coefficient of restitution e, contact time tc and peak contact force Fmax were presented in Figs. 10–12 respectively. Fig. 10 shows the effect of vi on e of ILB and OLB. For Hertz-type models with p N 1/4, the coefficient of restitution e tends to monotonically decrease with vi, and the reducing slope of the e-vi curve is steeper with p growing (i.e. Model Hunt N Li N KK and the current model). Their e decreasing with respect to vi coincides with the impact experiment. However, for other models with p = 1/4 or p = 0, the predicted e−vi curve is inconsistent with the experiment in Fig. 10, where Model Ts (p = 1/4) has e being independent of vi, and Model LH (p = 0) predicts e increasing with vi. As shown in Fig. 10, the e-vi curve predicted by the current model is generally in good agreement with the free fall test, though the predicted e at low vi is slightly overestimated for ILB. The slope of the predicted e-vi curve of the current model seems to be influenced by relaxation time. Fig. 11 compares the experimental and simulated results of tc versus vi. With the collision end criterion of zero contact force, each Hertz-type model predicts a monotonically decreasing tendency of tc with respect to vi, which coincides with experimental results of ILB and OLB. The simulated results of tc−vi by most contact models agree well with the experimental results, except Model Hunt and Model Li, of both which the simulated curves tc−vi are overestimated especially for heavilydamped ILB in Fig. 11(a). Among these Hertz-type models without involving relaxation time, the contact time tc predicted by a model with a bigger value of p is generally longer than that by another model with a smaller value of p, though their predicted curves tc−vi becomes closer to each other with vi increasing generally.

LH with p=0,

Fig. 12. Experimental and simulated relationships of peak contact force Fmax versus impact velocity vi ( Li with p=1, Hunt with p=3/2, • experiment [38]) for: (a) ILB, (b) OLB.

Ts with p=1/4,

LH with p=0,

KK with p=1/2,

Ts with p=1/4,

Current model,

KK with p=1/2,

Li with

Current model,

748

B. Jian et al. / Powder Technology 356 (2019) 735–749

The effect of vi on the peak force Fmax for ILB and OLB can be seen from Fig. 12, which indicates that Fmax tends to monotonically increase with vi. For heavily-damped ILB in Fig. 12(a), simulation results of Fmax−vi by Mode Li and the current model agree well with experimental results generally, while Model Hunt overestimates the results of Fmax−vi and other models underestimate the results of Fmax−vi compared with the experiment. For lightly-damped OLB in Fig. 12(b), simulation results of Fmax−vi by the current model match experimental results better than other Hertz-type models, which slightly underestimate the results of Fmax−vi compared with the experiment. Comparing Fig. 12(a) with Fig. 12(b), differences among simulation results of Fmax−vi by different Hertz-type contact models for heavily-damped ILB are much greater than those for lightly-damped OLB. Among those contact models regardless of relaxation time, the peak force Fmax predicted by a model with a bigger value of p is generally higher than that by another model with a smaller value of p, though their predicted curves Fmax−vi generally get closer to each other with vi decreasing. 5. Conclusion In this work, dimensional reasoning demonstrates that the dissipative force of viscoelastic particles increases in proportion to the displacement rate and the square root of the displacement of particles, which is consistent with previous derivations from continuum mechanics. The viscoelastic contact force model involving the effect of relaxation time behaves quite differently from other models. Except for the Model Hunt, other contact models all predict attractive force before displacement backing to zero. To avoid unfeasible attraction, the collision end criterion of a zero force instead of a zero displacement makes the coefficient of restitution higher and the contact time shorter, and causes a ‘residual’ displacement at the end of contact. The resulting ‘residual’ displacement should eventually recover after a viscoelastic contact, which is different from the permanent deformation of plastic impacts. The following conclusions can be drawn: (1) Involving the effect of relaxation time is to the advantage of the current damping term over others. The damping parameter β related to relaxation time significantly influences behaviors of the current model. The growth of β reduces the initial slope of the contact force, makes the peak force augmented and delayed, and increases both the coefficient of restitution e and the maximum displacement. (2) At fixed values of simulation parameters like e, m⁎ and kHertz, the predictions of viscoelastic contact models of different damping terms are quite different between each other, especially for a heavily damped body like ILB. Apart from the current model, a model with a larger displacement exponent p has higher predictions of both the maximum displacement and the peak contact force, and predicts a slower initial slope of contact force, and produces a longer contact time and a smaller residual displacement simultaneously. Besides, predictions of different contact models exhibit different dependence on impact velocity. (3) The simulation results of the current viscoelastic contact model match the impact test of ILB and OLB better than other models. This is due to that the damping term of the current model considers the effect of relaxation time, but others do not. So the current contact model is more appropriate to solving collisions of viscoelastic bodies. Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgements This research is supported by the National Natural Science Foundation of China under grant number 51575404. We are thankful to the

researchers of Jiangxi Naipu Mining Machinery & New Materials Co., Ltd of China for their useful discussions References [1] P.A. Cundall, O.D.L. Strack, Discrete numerical-model for granular assemblies, Geotechnique 29 (1) (1979) 47–65. [2] Y. Tsuji, T. Tanaka, T. Ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239–250. [3] T. Tanaka, T. Kawaguchi, Y. Tsuji, Discrete particle simulation of flow patterns in two-dimensional gas fluidized beds, Int. J. Mod. Phys. B 7 (09n10) (1993) 1889–1898. [4] S. Shrestha, S.B. Kuang, Z.Y. Zhou, Particle scale modelling of bubble properties in central air jet gas-solid fluidized beds, Powder Technol. 339 (2018) 70–80. [5] S. Nosewicz, J. Rojek, K. Pietrzak, M. Chmielewski, Viscoelastic discrete element model of powder sintering, Powder Technol. 246 (2013) 157–168. [6] E. Olsson, P.-L. Larsson, A numerical analysis of cold powder compaction based on micromechanical experiments, Powder Technol. 243 (2013) 71–78. [7] P.W. Cleary, J.E. Hilton, M.D. Sinnott, Modelling of industrial particle and multiphase flows, Powder Technol. 314 (2017) 232–252. [8] N.S. Weerasekara, M.S. Powell, P.W. Cleary, et al., The contribution of DEM to the science of comminution, Powder Technol. 248 (2013) 3–24. [9] P.W. Cleary, R.D. Morrison, Comminution mechanisms, particle shape evolution and collision energy partitioning in tumbling mills, Miner. Eng. 86 (2016) 75–95. [10] A.B. Stevens, C.M. Hrenya, Comparison of soft-sphere models to measurements of collision properties during normal impacts, Powder Technol. 154 (2005) 99–109. [11] H.P. Zhu, Z.Y. Zhou, R.Y. Yang, A.B. Yu, Discrete particle simulation of particulate systems: theoretical developments, Chem. Eng. Sci. 62 (13) (2007) 3378–3396. [12] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, V. Scherer, Review and extension of normal force models for the discrete element method, Powder Technol. 171 (2007) 157–173. [13] R. Kačianauskas, H. Kruggel-Emden, E. Zdancevičius, D. Markauskas, Comparative evaluation of normal viscoelastic contact force models in low velocity impact situations, Adv. Powder Technol. 27 (4) (2016) 1367–1379. [14] H. Hertz, Ueber die Berührung fester elastischer Körper, Journal für die reine und angewandte Mathematik (Crelle's Journal) 92 (1882) 156–171. [15] L.D. Landau, E.M. Lifshitz, Theory of Elasticity, 3rd English ed. Pergamon Press, New York, 1986. [16] O.R. Walton, R.L. Braun, Viscosity, granular temperature and stress calculations for shearing assemblies of inelastic, frictional disks, J. Rheol. 30 (1986) 949–980. [17] C. Thornton, Coefficient of restitution for collinear collisions of elastic-perfectly plastic spheres, J. Appl. Mech. 64 (2) (1997) 383–386. [18] L. Vu-Quoc, X. Zhang, An elastoplastic contact force-displacement model in the normal direction: displacement-driven version, PCOR. R. Soc. Lond. A 455 (1999) 4013–4044. [19] L. Vu-Quoc, X. Zhang, L. Lesburg, A normal force-displacement model for contacting spheres accounting for plastic deformation: force-driven formulation, ASME. J. Appl. Mech. 67 (2) (1999) 363–371. [20] K.H. Hunt, F.R.E. Crossley, Coefficient of restitution interpreted as damping in vibroimpact, J. Appl. Mech. 7 (1975) 440–445. [21] D. Zhang, W.J. Whiten, The calculation of contact forces between particles using spring and damping models, Powder Technol. 88 (1996) 59–64. [22] T. Schwager, T. Pöschel, Coefficient of restitution and linear–dashpot model revisited, Granul. Matter 9 (2007) 465–469. [23] C. Thornton, S.J. Cummins, P.W. Cleary, An investigation of the comparative behaviour of alternative contact force models during inelastic collisions, Powder Technol. 233 (2013) 30–46. [24] J. Shäfer, S. Dippel, D.E. Wolf, Force schemes in simulations of granular materials, J. Phys. I 6 (1) (1996) 5–20. [25] E. Zdancevičius, R. Kačianauskas, D. Zabulionis, Improvement of viscoelastic damping for the hertz contact of particles due to impact velocity, Procedia Eng. 172 (2017) 1286–1290. [26] J. Lee, H.J. Herrmann, Angle of repose and angle of marginal stability: moleculardynamics of granular particles, J. Phys. A Math. Gen. 26 (2) (1993) 373–383. [27] M.J. Lichtensteiger, Impact Analysis of Viscoelastic Spheres, Fruits and Vegetables with Rigid, Plane Surfaces, Unpublished Ph.D. Dissertation, The Ohio State University, Columbus, 1982. [28] G.M. Hu, Z.Y. Hu, B. Jian, L.P. Liu, H. Wan, On the determination of the damping coefficient of non-linear spring-dashpot system to model hertz contact for simulation by discrete element method, J. Comput. 6 (2011) 984–988. [29] E. Alizadeh, F. Bertrand, J. Chaouki, Development of a granular normal contact force model based on a non-Newtonian liquid filled dashpot, Powder Technol. 237 (2013) 202–212. [30] G. Kuwabara, K. Kono, Restitution coefficient in a collision between two spheres, Jpn. J. Appl. Phys. 26 (1987) 1230–1233. [31] N.V. Brilliantov, F. Spahn, J.-M. Hertzsch, T. Pö schel, Model for collisions in granular gases, Phys. Rev. E 53 (1996) 5382–5392. [32] Q.J. Zheng, H.P. Zhu, A.B. Yu, Finite element analysis of the contact forces between a viscoelastic sphere and rigid plane, Powder Technol. 226 (2012) 130–142. [33] Q.J. Zheng, Z.Y. Zhou, A.B. Yu, Contact forces between viscoelastic ellipsoidal particles, Powder Technol. 248 (2013) 25–33. [34] N.V. Brilliantov, A.V. Pimenova, D.S. Goldobin, A dissipative force between colliding viscoelastic bodies: Rrigorous approach, EPL Europhys. Lett. 109 (1) (2015) 14005.

B. Jian et al. / Powder Technology 356 (2019) 735–749 [35] B. Jian, G.M. Hu, Z.Q. Fang, H.J. Zhou, R. Xia, A normal contact force approach for viscoelastic spheres of the same material, Powder Technol. 350 (2019) 51–61. [36] D. Antypov, J. Elliott, On an analytical solution for the damped Hertzian spring, EPL Europhys. Lett. 94 (2011), 50004. [37] R. Kumar, A. Sarkar, W. Ketterhagen, B. Hancock, J. Curtis, C. Wassgren, Influence of normal contact force model on simulations of spherocylindrical particles, AICHE J. 64 (6) (2018) 1986–2001. [38] M.J. Lichtensteiger, R.G. Holmes, M.Y. Hamdy, J.L. Blaisdell, Impact parameters of spherical viscoelastic objects and tomatoes, Am. Soc. Agric. Eng. 31 (2) (1988) 595–602. [39] K.L. Johnson, Contact Mechanics, first ed. Cambridge University Press, Cambridge, 1987.

749

[40] W. Goldsmith, Impact, Edward Arnold Ltd, London, 1960. [41] H.M. Lankarani, P.E. Nikravesh, A contact force model with hysteresis damping for impact analysis of multibody systems, J. Mech. Des. 112 (1990) 369–376. [42] P. Flores, M. MacHado, M.T. Silva, J.M. Martins, On the continuous contact force models for soft materials in multibody dynamics, Multibody Syst. Dyn. 25 (3) (2011) 357–375. [43] A. Jordam Caserta, H.A. Navarro, L. Cabezas-Gómez, Damping coefficient and contact duration relations for continuous nonlinear spring-dashpot contact model in DEM, Powder Technol. 302 (2016) 462–479. [44] E. Olsson, D. Jelagin, A contact model for the normal force between viscoelastic particles in discrete element simulations, Powder Technol. 342 (2019) 985–991. [45] R. Cross, The bounce of a ball, Am. J. Phys. 67 (3) (1999) 222–227.