A viscoelastic model for flexible fibers with material damping

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 b

State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, Xi'an 710049, China 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

Article history: Received 27 October 2014 Received in revised form 5 January 2015 Accepted 23 February 2015 Available online 1 March 2015 Keywords: Discrete element method Flexible fiber Sphere-chain Damping

a b s t r a c t In order to quantify the material damping of deformed flexible fibers, a linear viscoelastic sphere-chain model based on the discrete element method is presented. In the model the interaction between any two neighboring spheres connected by a visual bond is modeled as springs and dash-pots which can be characterized by the bond stiffness and bond damping coefficient, respectively. Then the model is applied to simulate the bending of cantilever beam under static load and the motion of fiber in the simple shear flow, respectively, and the predictions agree well with those of classic analytic solutions, available experimental data and previous simulation results. Regarding to the damping behavior of deformed flexible fibers, the results suggest that the damping coefficient of an individual fiber mainly depends on the aspect ratio of the fiber and the viscous bond damping coefficient. A correlation is formulated to quantify the relationship between the damping coefficient of the local bond and that 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 Flexible fibers are critical raw materials in the manufacture of fiberreinforced composites, paper and fibrous materials [1–4]. The properties of these products such as tensile strength and permeability depend strongly on the orientation and spatial distribution of flexible fibers, and the network formation. For example, the entanglement of wood pulp fibers with each other could lead to small and concentrated flocs, and this heterogeneous microstructure impacts greatly the papermaking processes such as screening and sheet formation. Regarding the fabrication of porous metal fiber sintered sheet, the distribution homogeneity of flexible fibers in the gas conveying duct has a decisive influence on the pore structure of the packed sheet and then the macro properties of the final product. Thus, it is essential and practical to understand the dynamics of flexible fibers under different conditions, such as air flow and collisions. Those information may help provide some guidance for optimizing manufacturing procedures and improving product quality. Due to the difficulties resulting from the experiments in particle scale, numerical methods have become a proven and effective approach to explore the dynamic mechanisms of flexible fibers. As known, such microscopic information as the contact and entanglement between flexible fibers can be obtained, which may be hardly to be depicted by means of experimental analysis even with the aid of advanced technology. To describe the motion of flexible fibers, various numerical ⁎ Corresponding author. E-mail address: [email protected] (Y. Wang).

http://dx.doi.org/10.1016/j.powtec.2015.02.037 0032-5910/© 2015 Elsevier B.V. All rights reserved.

schemes have been proposed, such as immersed boundary method [5–7], slender body theory [8,9] and discrete element method [10–31]. In the frame of discrete element method, an individual flexible fiber is usually represented by a number of rigid segments which are consecutively connected by elastic bonds, ball-socket joints or hinges. And the deformation of flexible fiber can be realized by changing the relative displacements and rotations of the connected segments. Thus, the dynamics of the whole fiber can be easily obtained, provided that the motion of each segment is known. According to the shape of rigid segment, the flexible fiber model can be generalized into three major categories: sphere-chain model, spheroid-chain model and rod-chain model. In the early work of the sphere-chain model [10–19], spheroid-chain model [23,24] and rod-chain model [25–31], a set of constraint equation had to be solved to acquire the tangential force or constraint force between each pair of connected segments. And those early work was mainly limited to the simple shear flow with zero Reynolds number where the inertias of the fluid and fiber could be neglected. Compared to the spheroid-chain model and rod-chain model, the sphere-chain model has various advantages. For example, the detection of the collision and the calculation of the collision force as well as the theory of the hydraulic force and near-field force (i.e., Van der Waals force and electrostatic force) for spheres have been established thoroughly, while it is not the case for the spheroids and rods. Consequently, the spherechain model is more attractive than other numerical models. Recently, to minimize the computational effort of sphere-chain model, some researchers [20–22] presented new approaches to predict the constraint force directly instead of the solving of constraint equation. Hence, the flexible fiber model can be easier to be combined with the contact

176

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

model of spheres and extended to the fiber flow under various conditions. In the model adopted by Nguyen et al. [20], the Euler–Bernoulli beam theory was used to determine the bond forces acting on the spherical elements. Base on the Timoshenko beam theory, Obermayr et al. [21] calculated the tangential force directly through the quaternion algebra and bond orientation. Guo et al. [22] initially applied the bonded-particle model, which is originally proposed for the analysis of rock mechanics to the modeling of flexible fiber, and the interaction forces/torques were calculated in an incremental method. Prof. Guo 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 time step criterion for the first time into the flexible fiber model, the validity of which was also comprehensively verified by the cases of the cantilever beam under static and dynamic load. Similar to the energy dissipation in the collision between spheres, the deformed fiber will experience internal damping which can be caused by various combinations of fundamental physical mechanisms. For metals, these mechanisms include thermoelasticity on both micro and macro scales, grain boundary viscosity, point-defect relaxations, eddy-current effects, stress-induced ordering, and electronic effects [32]. Thus, when the flexible fiber is deformed, energy can be absorbed and dissipated by the material itself. Analogous to the importance of the collision damping in the packing and pneumatic conveying of spheres [33,34], the material damping may have noticeable effects on the deformation and motion of flexible fibers. However, only few previous literatures have taken into account the internal damping of flexible fibers. Though Obermayr et al. [21] had implemented the viscous damping force/torque into the numerical model, no further results and detailed discussions about the damping of flexible fibers were presented. By analogy with the viscous collision damping model between spheres, Guo et al. [35] proposed a bond damping to account the energy dissipation of deformation. The results showed that the bond damping had a significant effect on the collision between two flexible fibers and the shear stress of the shear flow of flexible fibers. Unlike the damping parameter in sphere–sphere collision which can be predicted by the restitution coefficient, the simulation parameters related to the bond damping are hypothetical in those work, without any links with the real material properties. To date, it is not clear how the local bond damping affects the whole damping of flexible fiber.

In this work, a linear viscoelastic sphere-chain scheme based on the discrete element method is presented and utilized to explore the internal damping of flexible fiber. A detailed description and validation of the numerical model are illustrated in Section 2 and Section 3, respectively. Then the relationship between the simulated damping parameter of the bond and the actual damping property of the fiber material is examined in Section 4. And this is followed by the conclusions in Section 5. 2. Numerical method Similar to the previous study presented by Yamamoto and Matsuoka [11], the flexible fiber is modeled to be made up of a number of identical spheres with their centers consecutive one-by-one on the symmetry, as illustrated in Fig. 1. Apparently, the number of sphere elements depends on the aspect ratio of flexible fiber rf which is defined as the ratio of total length of the bonds L to the diameter 2R. Thus, the dynamic characteristics of any individual fiber can be obtained by tracking the motion of each sphere. It is assumed that the interaction between any two neighboring rigid spheres is a kind of virtually elastic bond which obeys the linear elastic material law in every time step. And the relative displacements and rotations between the connected spheres lead to the deformation of the bonds and also the flexible fiber. In response, the bond forces/torques are induced and exerted on the spheres to resist the deformation. The translational and rotational motion of any constituent sphere in an individual fiber can be described by: mk

Jk

dV k b c h ¼ F k þ F k þ F k þ mk g dt

dωk b c ¼ Mk þ Mk dt

ð1Þ

ð2Þ

where mk and Jk are the mass and moment of inertia, respectively. Vk and ωk are the translational velocity and angular velocity, respectively. Fkb is the bond force which can be divided into the stretch and shear components. Mkb is the bond torque, arising from the bond shear forces and the bond bending/twisting torques. Fkc is the contact force due to the contacts with neighboring spheres from other fibers 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

easily determined by using Hertz–Mindlin contact model [36,37]. Mkc is the contact torque, originating from the tangential contact force and the rolling friction torque. Fhk is the hydraulic force if there is fluid around the fiber. For simplicity, the subscript k is omitted in the following paragraphs. Elastic deformation is modeled as linear springs (as shown in Fig. 1) which exhibit restoring forces/torques linear to the displacements/ rotations. To consider the large deformation conveniently, the calculation is decomposed into many small steps, where each step is geometrically linear. According to the sphere-chain model developed by Yamamoto and Matsuoka [11] and the bonded-particle model proposed by Potyondy and Cundall [38], the elastic restoring forces/torques are calculated in an incremental scheme as follows: b

b

dF n ¼ K n dδn b

b

dF t ¼ K t dδt

br

br

b

b

b

b

¼

¼

dM n ¼ K twist dθn

dM t ¼ K bend dθt

πR br EV n dt 2

πR br GV t dt 2

br

br

¼

πR3 br Gωn dt 4

πR br Eωt dt 8

ð4Þ

where d and ρ are the diameter and density of the flexible fiber. 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 = 27% ~ 32% tc). Hence, to ensure the numerical stability and maintain correctly dynamic behavior of the flexible fibers under various conditions, the critical bond time step should be used instead of the critical Rayleigh time step. Generally, 0.25tb is small enough to acquire reasonable results. Of course, if the bond damping coefficient is much larger, a smaller time step needs to be employed.

ð5Þ

ð6Þ

b

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi br 2mK n b V n

ð7Þ

b

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi br 2mK t b V t

ð8Þ

F dt ¼ ξb b

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi br 2JK twist b ωn

ð9Þ

b

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi br 2 JK bend b ωt

ð10Þ

M dn ¼ ξb M dt ¼ ξb

rffiffiffi ρ E

t b ¼ 0:8165d

where dFnb, dFtb, dMnb and dMtb are the increments of the stretch force Fnb, shear force Ftb, twisting torque Mnb and bending torque Mtb, respectively. b b Knb, Ktb, Ktwist and Kbend are the stretch stiffness, shear stiffness, twisting stiffness and bending stiffness, respectively, which can be represented br br br by Young's modulus E or shear modulus G. dδbr n , dδt , dθn and dθt are the increments of the stretch deformation, shear deformation, twisting deformation and bending deformation, respectively. Vnbr and Vtbr are the normal and tangential components of the relative translational vebr locity, respectively. ωbr n and ωt are the normal and tangential components of the relative angular velocity, respectively. It should be noted that the stretch force Fnb can also be calculated directly using the total stretch deformation. And our tests by far show no difference between those two schemes, though Ting et al. [39] noted that summing incremental normal contact forces of sphere–sphere collisions at each time step could lead to significant accumulated error. In order to take into account the energy dissipation resulting from the material damping of deformation, the viscous bond damping is introduced and modeled as dashpots (as shown in Fig. 1) which produce damping forces/torques linear to the translational/rotational velocities. Based on the vibration theory [40] and the damping models proposed by Obermayr et al. [21] and Guo et al. [35], the viscous damping forces/torques are given as: F dn ¼ ξb

linear spring contact model for sphere collision, the critical value of bond damping coefficient is not 1.0 anymore and may depend on the deformation process, as an incremental scheme is used to calculate the bond forces/torques due to deformation. It is also important to keep in mind that the bonded spheres are always the interior of the continuous fiber. Thus, in the view of the integrity of flexible fiber, the exact value of ξb makes not much sense, as long as the sphere-chain represented fiber has the same damping behavior as the real fiber. If a much large value of ξb is used, the sphere-chain represented fiber can quickly reach the state of mechanical equilibrium. The time step is determined from the view that the elastic wave should not travel longer than a single bond length within a time step. Thus, the critical bond time step can be given as follows [22]:

ð3Þ

3

¼

177

where ξb is the bond damping coefficient. In the vibration theory, many other material damping models have been proposed, such as the hysteretic 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

ð11Þ

3. Validation of flexible fiber model 3.1. Static mechanics of cantilever beam To validate the flexible fiber model, the deformation of a uniform cantilever beam with a concentrated vertical load at the free end is examined, as shown in Fig. 2. In this figure, ymax and Δ are the vertical and horizontal displacements of the free end respectively, while y is the vertical deflection of a member at any x. In the small deflection theory, the horizontal movement of the free end is ignored. Nevertheless, this is not the case for the beam with large deformation, which should be described by the large deflection theory. According to the small deflection theory, the deflection of the free end can be given as: ymax ¼

PL3 3EI

ð12Þ

where P and I are the vertical load force and moment of inertia, respectively. For the large deflection of a cantilever beam, the beam is assumed to be inextensible and the vertical load remained vertical during the

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 accuracy. 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 acquired through the integration of Eq. (13). GðxÞ 0 y ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1−½GðxÞ2 GðxÞ ¼ Z

L−Δ h 0

i P h 2 2 x −ðL−ΔÞ 2EI  0 2 i1=2 1þ y dx ¼ L

ð13Þ

ð14Þ ð15Þ 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 damping is applied to dissipate the kinetic energy quickly. The simulation results of the normalized deflection for fiber with aspect ratio of 10 are depicted in Fig. 3, in comparison with the analytic solutions of the deflection theory. It can be seen that the results in this work agree 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 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 deflection theory. The minimum and maximum departure from the large deflection theory is 0.24% and 0.66%, respectively. And the discrepancy is much smaller than that of the tests carried out by Nguyen et al. [20] and Obermayr et al. [21], suggesting that an incremental approach may be more superior than that of the scheme based on the beam theory and bond orientation. It should be noted that the extremely large 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 flexible fiber. However, as shown in Fig. 4, with the decrease of the aspect ratio, the error increases rapidly. For the case of rf = 6, the error exceeds 1% even at small deflections. The sudden increase of error for the case of rf = 6 (PL2/EI = 1–2.5) may be one drawback of the flexible fiber model,

and a similar phenomenon can be also found in Guo et al. [22]. On the other hand, considering the applicability of flexible fiber model under various conditions, the error criteria in the calibration process should be more rigorous. Therefore, the flexible fiber model in this study is recommended for flexible fibers with aspect ratio larger than 6. 3.2. Motion of fibers in shear flow The analytical solution proposed by Jeffery [42] showed that a single ellipsoidal particle in simple shear flow with low Reynolds number will exhibit periodic motion. The dimensionless orbit period of the rotation that solely depends on the aspect ratio can be given by:



T γ ¼ 2π r f þ

1 rf

! ð16Þ



where T is the period of rotation, γ is the shear rate. Bretherton [43] showed that the rotation of rigid fiber can also be predicted by Eq. (16), provided that an equivalent ellipsoidal aspect ratio calculated from the measured period of rotation is used, instead of the physical 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 motion 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

Fig. 3. The normalized deflection of the free end for fiber with aspect ratio of 10.

As a result of the material damping, the deformation of flexible fiber would decay intrinsically and gradually. And the damping intensity is described by the material damping coefficient, which only depends on

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

Fig. 6. The normalized displacement of the free end under damped vibration.

Fig. 5. The relationship between the orbit period and the aspect ratio.

the type or constituent of the material. Its value may be measured by various methods such as peak-amplitude method [32], and can be easily accessed through the articles and national standards for material. For example, the material damping coefficient is about 1E−4 to 1E−3 for metals and 0.01 to 0.05 for woods. In the numerical model, the damping of sphere-chain represented fiber is achieved through the damping force and torque on the bond which can be transmitted to its constituent spheres. The damping intensity of sphere-chain represented fiber characterized by the fiber damping coefficient can be controlled by the viscous bond damping coefficient ξb. Overall, the fiber damping coefficient is used to quantify the intensity of fiber damping caused by the bond damping in the simulation while the material damping coefficient is referred to as the intrinsic material property of the flexible fiber in reality. However, those two coefficients are both named as ξf for simplicity, as they should have the same value. Only in this way, can the damping dynamics of the real fiber be properly represented by the sphere-chain represented fiber. It may be possible to estimate the simulated parameter ξb according to the actual damping property ξf of the fiber material. This can be achieved by examining the damped vibration of the cantilever beam. Firstly, a load force is applied to the free end and a static bending of the beam is acquired, which is illustrated in Section 3.1. Then the load force is released and the beam will exhibit damped vibration under the bond damping. In the vibration process, the displacement y(t) of the free end normalized by the maximum deflection ymax is tracked, as shown in Fig. 6. According to the vibration theory, the fiber damping coefficient ξf can be given as: δ ξ f ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2πÞ2 þ δ2

179

damping coefficient can be easily obtained based on Eq. (17). Further details can be found in Meirovitch [40]. Fig. 7 shows the relationship between the bond damping coefficient ξb and fiber damping coefficient ξf. It is expected that ξf and ξb are not identical. In the numerical model, the damping of the whole fiber is weakened by the virtual bond. Thus, to produce a large damping of the whole fiber, the bond with much larger damping capacities is demanded. In fact, ξf is about two orders of magnitude smaller than ξb, as shown in Fig. 7. It is also clear that the fiber damping coefficient is approximately proportional to the bond damping coefficient for all cases. This may be caused by the linear viscous bond damping model used in this work. The major difference of the damping characters between the bond and whole fiber comes from the aspect ratio, as the influence of the modulus of fiber material and the external force can be directly reflected by the dynamics of the bond. Hence, it is believed that changes of the fiber damping due to the load force and shear modulus could be ignorable, compared to that of the aspect ratio. It is consistent with the results shown in Fig. 7. For example, as the load force increases to 0.02 N, the variation of the slope of ξf/ξb is less than 3.3%. To further explore the effect of the load force and shear modulus on the fiber damping, the damped vibration for fiber with rf = 10 and ξb = 1.0 under different values of the load force and shear modulus is also simulated. As

ð17Þ

where δ is known as the logarithmic decrement, which can be calculated based on the peak displacement: δ¼

1 A ln 1 j−1 Aj

ð18Þ

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 implemented 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 acquired through a least square fit. In this way, the simulated fiber

Fig. 7. The relationship between the bond damping coefficient and fiber damping 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. Compared 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 illustrated. Interestingly, the fiber damping coefficient shows an exponential drop with the increase of the aspect ratio. Therefore, the fiber damping coefficient could be considered as an exponential function of the aspect ratio while showing linear increase with the viscous bond damping coefficient. As discussed above, the simulated fiber damping coefficient should be equal to the material damping coefficient which is the intrinsic property of material. It may be possible to establish a link between the simulation parameter (i.e., ξb) and the actual damping properties of fiber material (i.e., ξf). Quantitatively, this correlation can be expressed in a form as follows:   ξf ¼ a exp −br f ξb

ð20Þ

where the left part represents the linearity between ξf and ξb while the 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 and 10, the values of ξf/ξb calculated from Eq. (20) are 0.01019 and 0.00738, respectively, consistent with the slope of 0.01045 and 0.00698 shown in Fig. 7. Additionally, the predicted results based on Eq. (20) are also in good agreement with the simulated data depicted in Fig. 9, which is case of ξb = 1.0. It is also true for other values of ξb, which is not shown here for simplicity. Therefore, Eq. (20) is accurate

Fig. 8. The variation of the fiber damping coefficient with the load force and shear modulus.

Fig. 9. The relationship between the fiber damping coefficient and the aspect ratio.

enough to estimate the viscous bond damping coefficient based on the material damping coefficient. When the aspect ratio is large enough (i.e., rf = 30), the simulated fiber damping may be much smaller than the practical value of the fiber material even with a very large bond damping coefficient, as shown in Fig. 9 and Eq. (20). Thus, to obtain appropriate fiber damping, 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 given as: R  b  ωbr b T r ¼ −ξr  F n   br  ω  2

ð21Þ

where ξr is the bond rolling friction coefficient and R/2 is the equivalent radius of two connected spheres. Fig. 10 shows the relationship between the fiber damping and bond rolling friction coefficient for fiber with ξb = 1.0. It can be seen that the damping of the whole fiber is improved, as the bond rolling friction torque is implemented. To compare the effectiveness of the rolling friction damping and viscous bond 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. 10. The effect of the bond rolling friction coefficient on the fiber damping.

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

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 reduce 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 illustrated 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 static packing bed composed of fibers with different deformation is obtained, 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 the deformation occurs, then has a sharp decrease as the fiber is unfolded in the expansion period. Finally, the bending angle again has a small increase due to packing under gravity and a plateau is obtained. The tendency is consistent with the packing process shown in Fig. 11. It is expected that a larger bond damping coefficient ξb corresponds to smaller bending angle. However, the effect of the bond rolling friction damping is almost equal to the bond damping. It is understandable that the difference between the rolling friction torque and viscous bond damping is not obvious for fibers with aspect ratio of 20, as shown in Fig. 10. 5. Conclusions A linear viscoelastic sphere-chain model based on the discrete element method is proposed for flexible fibers with material damping. In the verification, the bending of a cantilever beam under static load and the motion of fiber in the simple shear flow is examined. And the

181

Fig. 12. The variation of the mean bending angle with time.

relationship between the bond damping and fiber damping is also established. The main results from this study can be summarized as follows: 1) The flexible fiber model is sufficiently precise to characterize the flexible fiber. However, its accuracy drops dramatically with the decrease of the aspect ratio, and it is recommended for the flexible fiber with aspect ratio larger than 6. 2) The influence of the external force and shear modulus on the fiber damping could be ignorable. And the damping coefficient of the whole fiber is linear with the viscous bond damping coefficient while showing an exponential decrease with the aspect ratio, which can also be predicted by Eq. (20). 3) The rolling friction damping is more effective than the viscous bond damping for fibers with very large aspect ratio. 4) The packing of flexible fibers is much different from that of the rigid rods. And the bond damping is found to have an important influence 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 The authors are grateful to The Key Program of National Natural Science Foundation of China (grant no. 51134003) and Science Fund for Creative Research Groups of the National Natural Science Foundation of China (grant no. 51121092) for the financial support of this work.

References [1] S.B. Lindström, T. Uesaka, Particle-level simulation of forming of the fiber network in papermaking, Int. J. Eng. Sci. 46 (2008) 858–876. [2] A. Londono-Hurtado, T.A. Osswald, J.P. Hernandez-Ortiz, Modeling the behavior of fiber suspensions in the molding of polymer composites, J. Reinf. Plast. Compos. 30 (2011) 781–790. [3] R. Jarabo, E. Fuente, A. Moral, A. Blanco, L. Izquierdo, C. Negro, Effect of sepiolite on the flocculation of suspensions of fibre-reinforced cement, Cem. Concr. Res. 40 (2010) 1524–1530. [4] Z.P. Xi, J.L. Zhu, H.P. Tang, Q.B. Ao, H. Zhi, J.Y. Wang, C. Li, Progress of application researches of porous fiber metals, Materials 4 (2011) 816–824. [5] L.D. Zhu, C.S. Peskin, Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method, J. Comput. Phys. 179 (2002) 452–468. [6] L.D. Zhu, Scaling laws for drag of a compliant body in an incompressible viscous flow, J. Fluid Mech. 607 (2008) 387–400. [7] J.M. Stockie, S.I. Green, Simulating the motion of flexible pulp fibres using the immersed boundary method, J. Comput. Phys. 147 (1998) 147–165. [8] A.K. Tornberg, M.J. Shelley, Simulating the dynamics and interactions of flexible fibers in Stokes flows, J. Comput. Phys. 196 (2004) 8–40. [9] E.J. Hinch, The distortion of a flexible inextensible thread in a shearing flow, J. Fluid Mech. 74 (1976) 317–333. [10] S. Yamamoto, T. Matsuoka, Viscosity of dilute suspensions of rodlike particles: A numerical simulation method, J. Chem. Phys. 100 (1993) 3317–3324. [11] S. Yamamoto, T. Matsuoka, A method for dynamic simulation of rigid and flexible fibers in a flow field, J. Chem. Phys. 98 (1993) 644–650. [12] S. Yamamoto, T. Matsuoka, Dynamic simulation of fiber suspensions in shear-flow, J. Chem. Phys. 102 (1995) 2254–2260. [13] M. Yamanoi, J.M. Maia, Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers, J. Non-Newtonian Fluid Mech. 166 (2011) 457–468. [14] D. Qi, Direct simulations of flexible cylindrical fiber suspensions in finite Reynolds number flows, J. Chem. Phys. 125 (2006) 114901. [15] K.K. Kabanemi, J.F. Hétu, Effects of bending and torsion rigidity on deformation and breakage of flexible fibers: a direct simulation study, J. Chem. Phys. 136 (2012) 074903. [16] C.G. Joung, N. Phan-Thien, X.J. Fan, Direct simulation of flexible fibers, J. NonNewtonian Fluid Mech. 99 (2001) 1–36. [17] A.M. Slowicka, M.L. Ekiel-Jezewska, K. Sadlej, E. Wajnryb, Dynamics of fibers in a wide microchannel, J. Chem. Phys. 136 (2012) 044904. [18] H.F. Guo, C.W. Yu, B.G. Xu, S.Y. Li, Effect of the geometric parameters on a flexible fiber motion in a tangentially injected divergent swirling tube flow, Int. J. Eng. Sci. 49 (2011) 1033–1046. [19] H.F. Guo, B.G. Xu, C.W. Yu, S.Y. Li, Simulating the motion of a flexible fiber in 3D tangentially injected swirling airflow in a straight pipe—effects of some parameters, Int. J. Heat Mass Transf. 54 (2011) 4570–4579.

[20] D.H. Nguyen, N. Kang, J. Park, Validation of partially flexible rod model based on discrete element method using beam deflection and vibration, Powder Technol. 237 (2013) 147–152. [21] M. Obermayr, K. Dressler, C. Vrettos, P. Eberhard, A bonded-particle model for cemented sand, Comput. Geotech. 49 (2012) 299–313. [22] Y. Guo, C. Wassgren, B. Hancock, W. Ketterhagen, J. Curtis, Validation and time step determination of discrete element modeling of flexible fibers, Powder Technol. 249 (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. Multiphase Flow 37 (2011) 173–186. [25] C.F. Schmid, D.J. Klingenberg, Mechanical flocculation in flowing fiber suspensions, Phys. Rev. Lett. 84 (2000) 290–293. [26] C.F. Schmid, D.J. Klingenberg, Properties of fiber flocs with frictional and attractive interfiber forces, J. Colloid Interface Sci. 226 (2000) 136–144. [27] C.F. Schmid, L.H. Switzer, D.J. Klingenberg, Simulations of fiber flocculation: effects of fiber properties and interfiber friction, J. Rheol. 44 (2000) 781–809. [28] L.H. Switzer, D.J. Klingenberg, Rheology of sheared flexible fiber suspensions via fiber-level simulations, J. Rheol. 47 (2003) 759–778. [29] L.H. Switzer, D.J. Klingenberg, Flocculation in simulations of sheared fiber suspensions, Int. J. Multiphase Flow 30 (2004) 67–87. [30] S.B. Lindström, T. Uesaka, Simulation of the motion of flexible fibers in viscous fluid flow, Phys. Fluids 19 (2007) 113307. [31] G. Wang, W. Yu, C.X. Zhou, Optimization of the rod chain model to simulate the motions of a long flexible fiber in simple shear flows, Eur. J. Mech. B. Fluids 25 (2006) 337–347. [32] C.W. Bert, Material damping: an introductory review of mathematic measures and experimental technique, J. Sound Vib. 29 (1973) 129–153. [33] K. Li, S.B. Kuang, R.H. Pan, A.B. Yu, Numerical study of horizontal pneumatic conveying: Effect of material properties, Powder Technol. 251 (2014) 15–24. [34] Z.P. Zhang, L.F. Liu, Y.D. Yuan, A.B. Yu, A simulation study of the effects of dynamic variables on the packing of spheres, Powder Technol. 116 (2001) 23–32. [35] Y. Guo, J. Curtis, C. Wassgren, W. Ketterhagen, B. Hancock, Granular shear flows of flexible rod-like particles, Powders and Grains 2013, Proceedings of the 7th International Conference on Micromechanics of Granular Media, 2013, pp. 491–494. [36] 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. [37] A.O. Raji, Discrete Element Modelling of the Deformation of Bulk Agricultural Particulates, Newcastle University upon Tyne, UK, 1999. [38] D.O. Potyondy, P.A. Cundall, A bonded-particle model for rock, Int. J. Rock Mech. Min. Sci. 41 (2004) 1329–1364. [39] J.M. Ting, M. Khwaja, L.R. Meachum, J.D. Rowell, An ellipse-based discrete element model for granular materials, Int. J. Numer. Anal. Methods Geomech. 17 (1993) 603–623. [40] L. Meirovitch, Fundamentals of Vibrations, McGraw-Hill, New York, 2001. [41] D.G. Fertis, Nonlinear Structural Engineering, Springer Berlin Heidelberg, Germany, 2006. [42] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 102 (1922) 161–179. [43] F.P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, J. Fluid Mech. 14 (1962) 284–304. [44] E. Anczurowski, S.G. Mason, Particle motions in sheared suspensions. XXIV. Rotation of rigid spheroids and cylinders, Trans. Soc. Rheol. 12 (1968) 209–215. [45] Z.M. Ning, J.R. Melrose, A numerical model for simulating mechanical behavior of flexible fibers, J. Chem. Phys. 111 (1999) 10717–10726.