Powder Technology 276 (2015) 175–182

Contents lists available at ScienceDirect

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

A viscoelastic model for flexible fibers with material damping

Wenguang Nan a, Yueshe Wang a,⁎, Huiping Tang b
a
State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, Xi'an 710049, China
b
State Key Laboratory of Porous Metal Materials, Northwest Institute for Non-ferrous Metal Research, Xi'an 710016, China

a r t i c l e i n f o a b s t r a c t

Article history: In order to quantify the material damping of deformed flexible fibers, a linear viscoelastic sphere-chain model
Received 27 October 2014 based on the discrete element method is presented. In the model the interaction between any two neighboring
Received in revised form 5 January 2015 spheres connected by a visual bond is modeled as springs and dash-pots which can be characterized by the bond
Accepted 23 February 2015
stiffness and bond damping coefficient, respectively. Then the model is applied to simulate the bending of canti-
Available online 1 March 2015
lever beam under static load and the motion of fiber in the simple shear flow, respectively, and the predictions
Keywords:
agree well with those of classic analytic solutions, available experimental data and previous simulation results.
Discrete element method Regarding to the damping behavior of deformed flexible fibers, the results suggest that the damping coefficient
Flexible fiber of an individual fiber mainly depends on the aspect ratio of the fiber and the viscous bond damping coefficient. A
Sphere-chain correlation is formulated to quantify the relationship between the damping coefficient of the local bond and that
Damping of the flexible fiber. Generally, to enhance the global damping of fibers with extremely large aspect ratio, the bond
rolling friction torque is more effective than the viscous bond damping. Those findings will be meaningful for the
development of numerical modeling of flexible fibers.
© 2015 Elsevier B.V. All rights reserved.

1. Introduction schemes have been proposed, such as immersed boundary method

[5–7], slender body theory [8,9] and discrete element method [10–31].
Flexible fibers are critical raw materials in the manufacture of fiber- In the frame of discrete element method, an individual flexible fiber is
reinforced composites, paper and fibrous materials [1–4]. The prop- usually represented by a number of rigid segments which are consecu-
erties of these products such as tensile strength and permeability de- tively connected by elastic bonds, ball-socket joints or hinges. And the
pend strongly on the orientation and spatial distribution of flexible deformation of flexible fiber can be realized by changing the relative
fibers, and the network formation. For example, the entanglement of displacements and rotations of the connected segments. Thus, the dy-
wood pulp fibers with each other could lead to small and concentrated namics of the whole fiber can be easily obtained, provided that the mo-
flocs, and this heterogeneous microstructure impacts greatly the paper- tion of each segment is known. According to the shape of rigid segment,
making processes such as screening and sheet formation. Regarding the the flexible fiber model can be generalized into three major categories:
fabrication of porous metal fiber sintered sheet, the distribution homo- sphere-chain model, spheroid-chain model and rod-chain model. In the
geneity of flexible fibers in the gas conveying duct has a decisive influ- early work of the sphere-chain model [10–19], spheroid-chain model
ence on the pore structure of the packed sheet and then the macro [23,24] and rod-chain model [25–31], a set of constraint equation had
properties of the final product. Thus, it is essential and practical to un- to be solved to acquire the tangential force or constraint force between
derstand the dynamics of flexible fibers under different conditions, each pair of connected segments. And those early work was mainly
such as air flow and collisions. Those information may help provide limited to the simple shear flow with zero Reynolds number where
some guidance for optimizing manufacturing procedures and improv- the inertias of the fluid and fiber could be neglected. Compared to the
ing product quality. spheroid-chain model and rod-chain model, the sphere-chain model
Due to the difficulties resulting from the experiments in particle has various advantages. For example, the detection of the collision
scale, numerical methods have become a proven and effective approach and the calculation of the collision force as well as the theory of the hy-
to explore the dynamic mechanisms of flexible fibers. As known, such draulic force and near-field force (i.e., Van der Waals force and electro-
microscopic information as the contact and entanglement between static force) for spheres have been established thoroughly, while it is
flexible fibers can be obtained, which may be hardly to be depicted not the case for the spheroids and rods. Consequently, the sphere-
by means of experimental analysis even with the aid of advanced tech- chain model is more attractive than other numerical models. Recently,
nology. To describe the motion of flexible fibers, various numerical to minimize the computational effort of sphere-chain model, some re-
searchers [20–22] presented new approaches to predict the constraint
⁎ Corresponding author. force directly instead of the solving of constraint equation. Hence,
E-mail address: wangys@mail.xjtu.edu.cn (Y. Wang). the flexible fiber model can be easier to be combined with the contact

http://dx.doi.org/10.1016/j.powtec.2015.02.037
0032-5910/© 2015 Elsevier B.V. All rights reserved.
176 W. Nan et al. / Powder Technology 276 (2015) 175–182

model of spheres and extended to the fiber flow under various condi- In this work, a linear viscoelastic sphere-chain scheme based on the
tions. In the model adopted by Nguyen et al. [20], the Euler–Bernoulli discrete element method is presented and utilized to explore the inter-
beam theory was used to determine the bond forces acting on the nal damping of flexible fiber. A detailed description and validation of the
spherical elements. Base on the Timoshenko beam theory, Obermayr numerical model are illustrated in Section 2 and Section 3, respectively.
et al. [21] calculated the tangential force directly through the quater- Then the relationship between the simulated damping parameter of the
nion algebra and bond orientation. Guo et al. [22] initially applied bond and the actual damping property of the fiber material is examined
the bonded-particle model, which is originally proposed for the anal- in Section 4. And this is followed by the conclusions in Section 5.
ysis of rock mechanics to the modeling of flexible fiber, and the interac-
tion forces/torques were calculated in an incremental method. Prof. Guo 2. Numerical method
and his co-authors' work [22] is pioneering in the field of modeling of
flexible fiber in the frame of discrete element. They employed the Similar to the previous study presented by Yamamoto and Matsuoka
time step criterion for the first time into the flexible fiber model, the va- [11], the flexible fiber is modeled to be made up of a number of identical
lidity of which was also comprehensively verified by the cases of the spheres with their centers consecutive one-by-one on the symmetry, as
cantilever beam under static and dynamic load. illustrated in Fig. 1. Apparently, the number of sphere elements depends
Similar to the energy dissipation in the collision between spheres, on the aspect ratio of flexible fiber rf which is defined as the ratio of total
the deformed fiber will experience internal damping which can be length of the bonds L to the diameter 2R. Thus, the dynamic characteris-
caused by various combinations of fundamental physical mechanisms. tics of any individual fiber can be obtained by tracking the motion of
For metals, these mechanisms include thermoelasticity on both micro each sphere. It is assumed that the interaction between any two neigh-
and macro scales, grain boundary viscosity, point-defect relaxations, boring rigid spheres is a kind of virtually elastic bond which obeys the
eddy-current effects, stress-induced ordering, and electronic effects linear elastic material law in every time step. And the relative displace-
[32]. Thus, when the flexible fiber is deformed, energy can be absorbed ments and rotations between the connected spheres lead to the defor-
and dissipated by the material itself. Analogous to the importance of the mation of the bonds and also the flexible fiber. In response, the bond
collision damping in the packing and pneumatic conveying of spheres forces/torques are induced and exerted on the spheres to resist the de-
[33,34], the material damping may have noticeable effects on the defor- formation. The translational and rotational motion of any constituent
mation and motion of flexible fibers. However, only few previous liter- sphere in an individual fiber can be described by:
atures have taken into account the internal damping of flexible fibers.
Though Obermayr et al. [21] had implemented the viscous damping dV k b c h
mk ¼ F k þ F k þ F k þ mk g ð1Þ
force/torque into the numerical model, no further results and detailed dt
discussions about the damping of flexible fibers were presented. By
analogy with the viscous collision damping model between spheres, dωk b c
Jk ¼ Mk þ Mk ð2Þ
Guo et al. [35] proposed a bond damping to account the energy dissipa- dt
tion of deformation. The results showed that the bond damping had a
significant effect on the collision between two flexible fibers and the where mk and Jk are the mass and moment of inertia, respectively.
shear stress of the shear flow of flexible fibers. Unlike the damping pa- Vk and ωk are the translational velocity and angular velocity, respec-
rameter in sphere–sphere collision which can be predicted by the resti- tively. Fkb is the bond force which can be divided into the stretch and
tution coefficient, the simulation parameters related to the bond shear components. Mkb is the bond torque, arising from the bond
damping are hypothetical in those work, without any links with the shear forces and the bond bending/twisting torques. Fkc is the contact
real material properties. To date, it is not clear how the local bond force due to the contacts with neighboring spheres from other fibers
damping affects the whole damping of flexible fiber. or the same fiber without mutual bond. The contact force Fkc can be

Fig. 1. Model of the elastic force and torque between two connected spheres.
(Figure reproduced from Guo et al. [22]).
 W. Nan et al. / Powder Technology 276 (2015) 175–182 177

easily determined by using Hertz–Mindlin contact model [36,37]. Mkc is linear spring contact model for sphere collision, the critical value of
the contact torque, originating from the tangential contact force and the bond damping coefficient is not 1.0 anymore and may depend on the
rolling friction torque. Fhk is the hydraulic force if there is fluid around deformation process, as an incremental scheme is used to calculate
the fiber. For simplicity, the subscript k is omitted in the following the bond forces/torques due to deformation. It is also important to
paragraphs. keep in mind that the bonded spheres are always the interior of the
Elastic deformation is modeled as linear springs (as shown in Fig. 1) continuous fiber. Thus, in the view of the integrity of flexible fiber, the
which exhibit restoring forces/torques linear to the displacements/ exact value of ξb makes not much sense, as long as the sphere-chain rep-
rotations. To consider the large deformation conveniently, the calcula- resented fiber has the same damping behavior as the real fiber. If a much
tion is decomposed into many small steps, where each step is geometri- large value of ξb is used, the sphere-chain represented fiber can quickly
cally linear. According to the sphere-chain model developed by reach the state of mechanical equilibrium.
Yamamoto and Matsuoka [11] and the bonded-particle model proposed The time step is determined from the view that the elastic wave
by Potyondy and Cundall [38], the elastic restoring forces/torques are should not travel longer than a single bond length within a time step.
calculated in an incremental scheme as follows: Thus, the critical bond time step can be given as follows [22]:

πR rffiffiffi
b b
dF n ¼ K n dδn
br
¼
br
EV n dt ð3Þ ρ
2 t b ¼ 0:8165d ð11Þ
E
b b br πR br
dF t ¼ K t dδt ¼ GV t dt ð4Þ where d and ρ are the diameter and density of the flexible fiber.
2
It should be noted that critical bond time step tb is much smaller than
the critical Rayleigh time step tc of sphere–sphere collision (tb =
b b br πR3 br
dM n ¼ K twist dθn ¼ Gωn dt ð5Þ 27% ~ 32% tc). Hence, to ensure the numerical stability and maintain cor-
4
rectly dynamic behavior of the flexible fibers under various conditions,
πR
3 the critical bond time step should be used instead of the critical Rayleigh
b b br br
dM t ¼ K bend dθt ¼ Eωt dt ð6Þ time step. Generally, 0.25tb is small enough to acquire reasonable re-
8
sults. Of course, if the bond damping coefficient is much larger, a smaller
where dFnb, dFtb, dMnb and dMtb are the increments of the stretch force Fnb, time step needs to be employed.
shear force Ftb, twisting torque Mnb and bending torque Mtb, respectively.
Knb, Ktb, Ktwist
b b
and Kbend are the stretch stiffness, shear stiffness, twisting 3. Validation of flexible fiber model
stiffness and bending stiffness, respectively, which can be represented
by Young's modulus E or shear modulus G. dδbr br br br
n , dδt , dθn and dθt are 3.1. Static mechanics of cantilever beam
the increments of the stretch deformation, shear deformation, twisting
deformation and bending deformation, respectively. Vnbr and Vtbr are To validate the flexible fiber model, the deformation of a uniform
the normal and tangential components of the relative translational ve- cantilever beam with a concentrated vertical load at the free end is ex-
locity, respectively. ωbr br
n and ωt are the normal and tangential compo- amined, as shown in Fig. 2. In this figure, ymax and Δ are the vertical
nents of the relative angular velocity, respectively. It should be noted and horizontal displacements of the free end respectively, while y is
that the stretch force Fnb can also be calculated directly using the total the vertical deflection of a member at any x. In the small deflection the-
stretch deformation. And our tests by far show no difference between ory, the horizontal movement of the free end is ignored. Nevertheless,
those two schemes, though Ting et al. [39] noted that summing incre- this is not the case for the beam with large deformation, which should
mental normal contact forces of sphere–sphere collisions at each time be described by the large deflection theory.
step could lead to significant accumulated error. According to the small deflection theory, the deflection of the free
In order to take into account the energy dissipation resulting from end can be given as:
the material damping of deformation, the viscous bond damping is in-
troduced and modeled as dashpots (as shown in Fig. 1) which produce PL3
ymax ¼ ð12Þ
damping forces/torques linear to the translational/rotational velocities. 3EI
Based on the vibration theory [40] and the damping models proposed
by Obermayr et al. [21] and Guo et al. [35], the viscous damping where P and I are the vertical load force and moment of inertia,
forces/torques are given as: respectively.
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi For the large deflection of a cantilever beam, the beam is assumed to
b
F dn ¼ ξb 2mK n b V n
br
ð7Þ be inextensible and the vertical load remained vertical during the

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b br
F dt ¼ ξb 2mK t b V t ð8Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b br
M dn ¼ ξb 2JK twist b ωn ð9Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b br
M dt ¼ ξb 2 JK bend b ωt ð10Þ

where ξb is the bond damping coefficient. In the vibration theory, many

other material damping models have been proposed, such as the hyster-
etic damping model and convolution integrals model [32], which are
more superior and accurate than the viscous damping model. However,
they are all calculated in a complex domain, and could not be used
directly in the frame of discrete element method. It should be noted
that the bond damping coefficient is just a model parameter to control
the damping character of sphere-chain represented fiber. Compared to Fig. 2. Bending of the flexible fiber under a point force.
178 W. Nan et al. / Powder Technology 276 (2015) 175–182

deformation of the member. Thus, according to the Elastica solution

[41], the deflection of the free end can be expressed by Eq. (13)–
Eq. (15), which can be solved through numerical method. An extremely
small value of Δ in Eq. (14) is assumed and the integration in Eq. (15) is
carried out numerically by using Simpson's rule with sufficient accura-
cy. If the value of left part (estimated value of beam length) is much
less or more than that of the right part (actual value of beam length),
a new value of Δ will be assumed. And this trial-and-error procedure
will be repeated in a bisection method until the relative error is less
than a critical value. Then the theoretical deflection can be easily ac-
quired through the integration of Eq. (13).

0 GðxÞ
y ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð13Þ
1−½GðxÞ2

P h 2 2
i
GðxÞ ¼ x −ðL−ΔÞ ð14Þ
2EI
Z L−Δ h  0 2 i1=2
1þ y dx ¼ L ð15Þ
0
Fig. 4. The errors of the simulation results compared to the large deflection theory.

The quasi-static loading process is achieved by increasing the load

gradually to its set value. To produce static deformation, a large bond and a similar phenomenon can be also found in Guo et al. [22]. On
damping is applied to dissipate the kinetic energy quickly. The simula- the other hand, considering the applicability of flexible fiber model
tion results of the normalized deflection for fiber with aspect ratio of under various conditions, the error criteria in the calibration process
10 are depicted in Fig. 3, in comparison with the analytic solutions of should be more rigorous. Therefore, the flexible fiber model in this
the deflection theory. It can be seen that the results in this work agree study is recommended for flexible fibers with aspect ratio larger than 6.
well with both small and large deflection theory when the deflection
is small (ymax/L b 0.2), with errors generally less than 0.25% (as shown 3.2. Motion of fibers in shear flow
in Fig. 4). As the applied load increases, large deflection of the free end
arises, and the simulation results can be better predicted by the large The analytical solution proposed by Jeffery [42] showed that a single
deflection theory. The minimum and maximum departure from the ellipsoidal particle in simple shear flow with low Reynolds number will
large deflection theory is 0.24% and 0.66%, respectively. And the discrep- exhibit periodic motion. The dimensionless orbit period of the rotation
ancy is much smaller than that of the tests carried out by Nguyen et al. that solely depends on the aspect ratio can be given by:
[20] and Obermayr et al. [21], suggesting that an incremental approach
!
may be more superior than that of the scheme based on the beam the-  1
ory and bond orientation. It should be noted that the extremely large T γ ¼ 2π r f þ ð16Þ
rf
bending (i.e., ymax/L N 0.6, as shown in Fig. 3) appears rarely in the actual
flow composed of numerous fibers. Thus, the flexible fiber model in this

work is sufficiently precise to describe the mechanical behavior of flex- where T is the period of rotation, γ is the shear rate. Bretherton [43]
ible fiber. However, as shown in Fig. 4, with the decrease of the aspect showed that the rotation of rigid fiber can also be predicted by
ratio, the error increases rapidly. For the case of rf = 6, the error exceeds Eq. (16), provided that an equivalent ellipsoidal aspect ratio calculated
1% even at small deflections. The sudden increase of error for the case of from the measured period of rotation is used, instead of the physical
rf = 6 (PL2/EI = 1–2.5) may be one drawback of the flexible fiber model, aspect ratio.
In the simulation, the free draining approximation [11] is utilized to
evaluate the hydrodynamic force and torque exerted on the sphere in
creeping flow. And the shear modulus is chosen large enough to ensure
there is no significant deformation occurring. Fig. 5 shows the variations
of the dimensionless orbit period of fiber rotation with the aspect ratio. It
is clear that the trend of the prediction results is similar with Jeffery's
theoretic solution and the experimental data of Anczurowski and
Mason [44]. As the sphere chain model has to be used and the shear
flow is considered undisturbed by the fiber motion, the drag and torque
of a fiber may be underestimated [45]. Consequently, the rotational mo-
tion is slowdown and the orbit period is longer than it should be. Thus, a
deviation of the simulated orbit period from the experimental data of
Anczurowski and Mason [44] can be observed, especially for fibers
with a large aspect ratio. Compared with the simulation data published
in the literature, there is also an excellent agreement between our results
and those from Ross and Klingenberg [23] and Kabanemi and Hétu [15].

4. Damping of flexible fiber

As a result of the material damping, the deformation of flexible fiber

would decay intrinsically and gradually. And the damping intensity is
Fig. 3. The normalized deflection of the free end for fiber with aspect ratio of 10. described by the material damping coefficient, which only depends on
 W. Nan et al. / Powder Technology 276 (2015) 175–182 179

Fig. 6. The normalized displacement of the free end under damped vibration.
Fig. 5. The relationship between the orbit period and the aspect ratio.

damping coefficient can be easily obtained based on Eq. (17). Further

the type or constituent of the material. Its value may be measured by details can be found in Meirovitch [40].
various methods such as peak-amplitude method [32], and can be easily Fig. 7 shows the relationship between the bond damping coefficient
accessed through the articles and national standards for material. For ξb and fiber damping coefficient ξf. It is expected that ξf and ξb are not
example, the material damping coefficient is about 1E−4 to 1E−3 for identical. In the numerical model, the damping of the whole fiber is
metals and 0.01 to 0.05 for woods. In the numerical model, the damping weakened by the virtual bond. Thus, to produce a large damping of
of sphere-chain represented fiber is achieved through the damping the whole fiber, the bond with much larger damping capacities is
force and torque on the bond which can be transmitted to its constitu- demanded. In fact, ξf is about two orders of magnitude smaller than ξb,
ent spheres. The damping intensity of sphere-chain represented fiber as shown in Fig. 7. It is also clear that the fiber damping coefficient is ap-
characterized by the fiber damping coefficient can be controlled by proximately proportional to the bond damping coefficient for all cases.
the viscous bond damping coefficient ξb. Overall, the fiber damping This may be caused by the linear viscous bond damping model used in
coefficient is used to quantify the intensity of fiber damping caused by this work.
the bond damping in the simulation while the material damping coeffi- The major difference of the damping characters between the bond
cient is referred to as the intrinsic material property of the flexible fiber and whole fiber comes from the aspect ratio, as the influence of the
in reality. However, those two coefficients are both named as ξf for modulus of fiber material and the external force can be directly reflected
simplicity, as they should have the same value. Only in this way, can by the dynamics of the bond. Hence, it is believed that changes of the
the damping dynamics of the real fiber be properly represented by the fiber damping due to the load force and shear modulus could be ignor-
sphere-chain represented fiber. It may be possible to estimate the able, compared to that of the aspect ratio. It is consistent with the results
simulated parameter ξb according to the actual damping property ξf of shown in Fig. 7. For example, as the load force increases to 0.02 N, the
the fiber material. This can be achieved by examining the damped vibra- variation of the slope of ξf/ξb is less than 3.3%. To further explore the
tion of the cantilever beam. Firstly, a load force is applied to the free end effect of the load force and shear modulus on the fiber damping, the
and a static bending of the beam is acquired, which is illustrated in damped vibration for fiber with rf = 10 and ξb = 1.0 under different
Section 3.1. Then the load force is released and the beam will exhibit values of the load force and shear modulus is also simulated. As
damped vibration under the bond damping. In the vibration process,
the displacement y(t) of the free end normalized by the maximum de-
flection ymax is tracked, as shown in Fig. 6. According to the vibration
theory, the fiber damping coefficient ξf can be given as:

δ
ξ f ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð17Þ
ð2πÞ2 þ δ2

where δ is known as the logarithmic decrement, which can be calculat-

ed based on the peak displacement:

1 A
δ¼ ln 1 ð18Þ
j−1 Aj

where Aj is the displacement of the jth peak of the curve shown in Fig. 6.
According to Meirovitch [40], the calculation of δ should be implement-
ed in an error minimum scheme:

ln A j ¼ ln A1 −δð j−1Þ ð19Þ

Based on Eq. (19), a plot of lnAj versus j can be obtained. Then δ is Fig. 7. The relationship between the bond damping coefficient and fiber damping
acquired through a least square fit. In this way, the simulated fiber coefficient.
180 W. Nan et al. / Powder Technology 276 (2015) 175–182

shown in Fig. 8, with the increasing of the load force, the fiber damping
coefficient ξf decrease slowly. For example, when the force changes
from 0.001 N to 0.008 N, the variation of ξf is 5%. Though the
asymptotical trend is not observed in the whole range of the load
force considered in this work, the overall variation of ξf is small. Com-
pared to the load force, the influence of the shear modulus is much
smaller. And ξf tends to be constant when the shear modulus is larger
than 0.2 MPa. Therefore, it demonstrates once again that the changes
of the fiber damping can be ignorable when the external force and
shear modulus vary. Hence, the bond damping coefficient needs not to
be changed in the case that the shear modulus of fiber has to be reduced
to manage properly computational effort. As discussed above, the
discrepancy of the whole fiber damping from the local bond damping
decreases dramatically as the aspect ratio becomes smaller, which can
be concluded from Fig. 7. With the decreasing of the aspect ratio, the
slope of ξf/ξb increases sharply, and so the difference between ξf and
ξb is reduced. It can also be concluded from Fig. 9, where the relationship
between the aspect ratio and the fiber damping with ξb = 1.0 is illus-
trated. Interestingly, the fiber damping coefficient shows an exponential
drop with the increase of the aspect ratio. Fig. 9. The relationship between the fiber damping coefficient and the aspect ratio.
Therefore, the fiber damping coefficient could be considered as an
exponential function of the aspect ratio while showing linear increase enough to estimate the viscous bond damping coefficient based on the
with the viscous bond damping coefficient. As discussed above, the material damping coefficient.
simulated fiber damping coefficient should be equal to the material When the aspect ratio is large enough (i.e., rf = 30), the simulated
damping coefficient which is the intrinsic property of material. It may fiber damping may be much smaller than the practical value of the
be possible to establish a link between the simulation parameter fiber material even with a very large bond damping coefficient, as
(i.e., ξb) and the actual damping properties of fiber material (i.e., ξf). shown in Fig. 9 and Eq. (20). Thus, to obtain appropriate fiber damping,
Quantitatively, this correlation can be expressed in a form as follows: the viscous bond damping may be not enough. Analogous to the
collision damping between spheres, a damped rolling friction torque
  Trb is complemented to the damping model of the bond, which is
ξf
¼ a exp −br f ð20Þ given as:
ξb

b R  b  ωbr
where the left part represents the linearity between ξf and ξb while the T r ¼ −ξr  F n   br  ð21Þ
2 ω 
right term stands for the influence of the aspect ratio. Based on the curve
fitting of the simulation results, the values of parameters a and b
are 0.03709 and 0.1615, respectively. For fibers with aspect ratio of 8 where ξr is the bond rolling friction coefficient and R/2 is the equivalent
and 10, the values of ξf/ξb calculated from Eq. (20) are 0.01019 and radius of two connected spheres. Fig. 10 shows the relationship be-
0.00738, respectively, consistent with the slope of 0.01045 and tween the fiber damping and bond rolling friction coefficient for fiber
0.00698 shown in Fig. 7. Additionally, the predicted results based on with ξb = 1.0. It can be seen that the damping of the whole fiber is im-
Eq. (20) are also in good agreement with the simulated data depicted proved, as the bond rolling friction torque is implemented. To compare
in Fig. 9, which is case of ξb = 1.0. It is also true for other values of ξb, the effectiveness of the rolling friction damping and viscous bond
which is not shown here for simplicity. Therefore, Eq. (20) is accurate damping, we also plot ξf based on hypothetical linearity that ξr and ξb
have the same effect on the fiber damping. Clearly, the enhancement

Fig. 8. The variation of the fiber damping coefficient with the load force and shear
modulus. Fig. 10. The effect of the bond rolling friction coefficient on the fiber damping.
 W. Nan et al. / Powder Technology 276 (2015) 175–182 181

of fiber damping can be more significant for fibers with larger aspect
ratio. Regarding fiber with aspect ratio of 30, the rolling friction
damping can produce much larger fiber damping than the viscous
bond damping. However, it is not obvious for flexible fiber with rf =
20. Moreover, the bond rolling friction damping will tend to asymptotic
value when ξr is large enough (i.e., ξr N 2.0).
To examine the effect of bond damping on the systems composed
of numerous flexible fibers, the dynamic characteristics of the flexible
fibers under gravity in the packing process are also investigated. To re-
duce the computational effort, 2000 flexible fibers with aspect ratio of
20 are used. And the restitution coefficient and friction coefficient are
0.4 and 0.6, respectively. Fig. 11 shows the packing process of flexible
fibers which are colored by the vertical component of the velocity. As il-
lustrated in Fig. 11(a)–(b), the deformation of flexible fibers occurs as
soon as the fibers reach the wall and collide with it. Then the deformed
fibers begin to pack on the bottom, and their motion is hindered by the
wall or packed fibers, resulting a rising bed. However, as shown in
Fig. 11(c), fibers are bounced back due to the large deformation energy.
Thus, the packed bed has a period of expansion after the densification
process, which is much different to the packing of rigid rods. Finally, a Fig. 12. The variation of the mean bending angle with time.
static packing bed composed of fibers with different deformation is ob-
tained, as shown in Fig. 11(d). The mean value of the bending angle for
every bond is plotted in Fig. 12. The bending angle increases firstly as relationship between the bond damping and fiber damping is also
the deformation occurs, then has a sharp decrease as the fiber is unfold- established. The main results from this study can be summarized as
ed in the expansion period. Finally, the bending angle again has a small follows:
increase due to packing under gravity and a plateau is obtained. The ten-
dency is consistent with the packing process shown in Fig. 11. It is ex- 1) The flexible fiber model is sufficiently precise to characterize the
pected that a larger bond damping coefficient ξb corresponds to flexible fiber. However, its accuracy drops dramatically with the de-
smaller bending angle. However, the effect of the bond rolling friction crease of the aspect ratio, and it is recommended for the flexible fiber
damping is almost equal to the bond damping. It is understandable with aspect ratio larger than 6.
that the difference between the rolling friction torque and viscous 2) The influence of the external force and shear modulus on the fiber
bond damping is not obvious for fibers with aspect ratio of 20, as damping could be ignorable. And the damping coefficient of the
shown in Fig. 10. whole fiber is linear with the viscous bond damping coefficient
while showing an exponential decrease with the aspect ratio,
5. Conclusions which can also be predicted by Eq. (20).
3) The rolling friction damping is more effective than the viscous bond
A linear viscoelastic sphere-chain model based on the discrete ele- damping for fibers with very large aspect ratio.
ment method is proposed for flexible fibers with material damping. In 4) The packing of flexible fibers is much different from that of the rigid
the verification, the bending of a cantilever beam under static load rods. And the bond damping is found to have an important influence
and the motion of fiber in the simple shear flow is examined. And the on the fiber deformation of the packing structure.

Fig. 11. The packing process of flexible fibers.

182 W. Nan et al. / Powder Technology 276 (2015) 175–182

Acknowledgments [20] D.H. Nguyen, N. Kang, J. Park, Validation of partially flexible rod model based on dis-
crete element method using beam deflection and vibration, Powder Technol. 237
(2013) 147–152.
The authors are grateful to The Key Program of National Natural Sci- [21] M. Obermayr, K. Dressler, C. Vrettos, P. Eberhard, A bonded-particle model for
ence Foundation of China (grant no. 51134003) and Science Fund for cemented sand, Comput. Geotech. 49 (2012) 299–313.
[22] Y. Guo, C. Wassgren, B. Hancock, W. Ketterhagen, J. Curtis, Validation and time step
Creative Research Groups of the National Natural Science Foundation determination of discrete element modeling of flexible fibers, Powder Technol. 249
of China (grant no. 51121092) for the financial support of this work. (2013) 386–395.
[23] R.F. Ross, D.J. Klingenberg, Dynamic simulation of flexible fibers composed of linked
rigid bodies, J. Chem. Phys. 106 (1997) 2949–2960.
[24] A. Vakil, S.I. Green, Flexible fiber motion in the flow field of a cylinder, Int. J. Multi-
References phase Flow 37 (2011) 173–186.
[25] C.F. Schmid, D.J. Klingenberg, Mechanical flocculation in flowing fiber suspensions,
[1] S.B. Lindström, T. Uesaka, Particle-level simulation of forming of the fiber network in Phys. Rev. Lett. 84 (2000) 290–293.
papermaking, Int. J. Eng. Sci. 46 (2008) 858–876. [26] C.F. Schmid, D.J. Klingenberg, Properties of fiber flocs with frictional and attractive
[2] A. Londono-Hurtado, T.A. Osswald, J.P. Hernandez-Ortiz, Modeling the behavior of interfiber forces, J. Colloid Interface Sci. 226 (2000) 136–144.
fiber suspensions in the molding of polymer composites, J. Reinf. Plast. Compos. [27] C.F. Schmid, L.H. Switzer, D.J. Klingenberg, Simulations of fiber flocculation: effects of
30 (2011) 781–790. fiber properties and interfiber friction, J. Rheol. 44 (2000) 781–809.
[3] R. Jarabo, E. Fuente, A. Moral, A. Blanco, L. Izquierdo, C. Negro, Effect of sepiolite on [28] L.H. Switzer, D.J. Klingenberg, Rheology of sheared flexible fiber suspensions via
the flocculation of suspensions of fibre-reinforced cement, Cem. Concr. Res. 40 fiber-level simulations, J. Rheol. 47 (2003) 759–778.
(2010) 1524–1530. [29] L.H. Switzer, D.J. Klingenberg, Flocculation in simulations of sheared fiber suspen-
[4] Z.P. Xi, J.L. Zhu, H.P. Tang, Q.B. Ao, H. Zhi, J.Y. Wang, C. Li, Progress of application re- sions, Int. J. Multiphase Flow 30 (2004) 67–87.
searches of porous fiber metals, Materials 4 (2011) 816–824. [30] S.B. Lindström, T. Uesaka, Simulation of the motion of flexible fibers in viscous fluid
[5] L.D. Zhu, C.S. Peskin, Simulation of a flapping flexible filament in a flowing soap film flow, Phys. Fluids 19 (2007) 113307.
by the immersed boundary method, J. Comput. Phys. 179 (2002) 452–468. [31] G. Wang, W. Yu, C.X. Zhou, Optimization of the rod chain model to simulate the mo-
[6] L.D. Zhu, Scaling laws for drag of a compliant body in an incompressible viscous tions of a long flexible fiber in simple shear flows, Eur. J. Mech. B. Fluids 25 (2006)
flow, J. Fluid Mech. 607 (2008) 387–400. 337–347.
[7] J.M. Stockie, S.I. Green, Simulating the motion of flexible pulp fibres using the im- [32] C.W. Bert, Material damping: an introductory review of mathematic measures and
mersed boundary method, J. Comput. Phys. 147 (1998) 147–165. experimental technique, J. Sound Vib. 29 (1973) 129–153.
[8] A.K. Tornberg, M.J. Shelley, Simulating the dynamics and interactions of flexible fi- [33] K. Li, S.B. Kuang, R.H. Pan, A.B. Yu, Numerical study of horizontal pneumatic convey-
bers in Stokes flows, J. Comput. Phys. 196 (2004) 8–40. ing: Effect of material properties, Powder Technol. 251 (2014) 15–24.
[9] E.J. Hinch, The distortion of a flexible inextensible thread in a shearing flow, J. Fluid [34] Z.P. Zhang, L.F. Liu, Y.D. Yuan, A.B. Yu, A simulation study of the effects of dynamic
Mech. 74 (1976) 317–333. variables on the packing of spheres, Powder Technol. 116 (2001) 23–32.
[10] S. Yamamoto, T. Matsuoka, Viscosity of dilute suspensions of rodlike particles: A nu- [35] Y. Guo, J. Curtis, C. Wassgren, W. Ketterhagen, B. Hancock, Granular shear flows of
merical simulation method, J. Chem. Phys. 100 (1993) 3317–3324. flexible rod-like particles, Powders and Grains 2013, Proceedings of the 7th Interna-
[11] S. Yamamoto, T. Matsuoka, A method for dynamic simulation of rigid and flexible fi- tional Conference on Micromechanics of Granular Media, 2013, pp. 491–494.
bers in a flow field, J. Chem. Phys. 98 (1993) 644–650. [36] Y. Tsuji, T. Tanaka, T. Ishida, Lagrangian numerical-simulation of plug flow of cohe-
[12] S. Yamamoto, T. Matsuoka, Dynamic simulation of fiber suspensions in shear-flow, J. sionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239–250.
Chem. Phys. 102 (1995) 2254–2260. [37] A.O. Raji, Discrete Element Modelling of the Deformation of Bulk Agricultural Partic-
[13] M. Yamanoi, J.M. Maia, Stokesian dynamics simulation of the role of hydrodynamic ulates, Newcastle University upon Tyne, UK, 1999.
interactions on the behavior of a single particle suspending in a Newtonian fluid. [38] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock Mech.
Part 1. 1D flexible and rigid fibers, J. Non-Newtonian Fluid Mech. 166 (2011) Min. Sci. 41 (2004) 1329–1364.
457–468. [39] J.M. Ting, M. Khwaja, L.R. Meachum, J.D. Rowell, An ellipse-based discrete element
[14] D. Qi, Direct simulations of flexible cylindrical fiber suspensions in finite Reynolds model for granular materials, Int. J. Numer. Anal. Methods Geomech. 17 (1993)
number flows, J. Chem. Phys. 125 (2006) 114901. 603–623.
[15] K.K. Kabanemi, J.F. Hétu, Effects of bending and torsion rigidity on deformation and [40] L. Meirovitch, Fundamentals of Vibrations, McGraw-Hill, New York, 2001.
breakage of flexible fibers: a direct simulation study, J. Chem. Phys. 136 (2012) [41] D.G. Fertis, Nonlinear Structural Engineering, Springer Berlin Heidelberg, Germany,
074903. 2006.
[16] C.G. Joung, N. Phan-Thien, X.J. Fan, Direct simulation of flexible fibers, J. Non- [42] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proceed-
Newtonian Fluid Mech. 99 (2001) 1–36. ings of the Royal Society of London. Series A, Containing Papers of a Mathematical
[17] A.M. Slowicka, M.L. Ekiel-Jezewska, K. Sadlej, E. Wajnryb, Dynamics of fibers in a and Physical Character 102 (1922) 161–179.
wide microchannel, J. Chem. Phys. 136 (2012) 044904. [43] F.P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds num-
[18] H.F. Guo, C.W. Yu, B.G. Xu, S.Y. Li, Effect of the geometric parameters on a flexible ber, J. Fluid Mech. 14 (1962) 284–304.
fiber motion in a tangentially injected divergent swirling tube flow, Int. J. Eng. Sci. [44] E. Anczurowski, S.G. Mason, Particle motions in sheared suspensions. XXIV. Rotation
49 (2011) 1033–1046. of rigid spheroids and cylinders, Trans. Soc. Rheol. 12 (1968) 209–215.
[19] H.F. Guo, B.G. Xu, C.W. Yu, S.Y. Li, Simulating the motion of a flexible fiber in 3D tan- [45] Z.M. Ning, J.R. Melrose, A numerical model for simulating mechanical behavior of
gentially injected swirling airflow in a straight pipe—effects of some parameters, Int. flexible fibers, J. Chem. Phys. 111 (1999) 10717–10726.
J. Heat Mass Transf. 54 (2011) 4570–4579.