Determination of the Normal Contact Stiffness and Integration Time Step for the Finite Element Modeling of Bristle-Surface Interaction

Libardo V. Vanegas-Useche, Magd M. Abdel-Wahab and Graham A. Parker

Abstract: In finite element modeling of impact, it is necessary to define appropriate values of the normal contact stiffness, Kn, and the Integration Time Step (ITS). Because impacts are usually of very short duration, very small ITSs are required. Moreover, the selection of a suitable value of Kn is a critical issue, as the impact behavior depends dramatically on this parameter. In this work, a number of experimental tests and finite element analyses have been performed in order to obtain an appropriate value of Kn for the interaction between a bristle of a gutter brush for road sweeping and a concrete surface. Furthermore, a suitable ITS is determined. The experiments consist of releasing a steel bristle that is placed vertically at a certain distance from a concrete surface and tracking the impact. Similarly, in the finite element analyses, a beam is modeled in free fall and impacting a surface; contact and target elements are attached to the beam and the surface, respectively. The results of the experiments and the modeling are integrated through the principle of conservation of energy, the principle of linear impulse and momentum, and Newton’s second law. The results demonstrate that, for the case studied,Kn and the impact time tend to be independent of the velocity just before impact and that Kn has a very large variation, as concrete is a composite material with a rough surface.Also, the ratio between the largest height of the bristle after impact and the initial height tends to be constant.

Keywords: Brush, street sweeping, finite element modeling, contact mechanics.

1 Introduction

The street sweeping activity, which is sometimes performed by lorry-type sweepers, is an important service both for aesthetic purposes and due to public hygiene [Vanegas-Useche,Abdel-Wahab and Parker (2010)]. These sweepers normally have a suction unit, a wide broom, and a gutter brush; this brush has to sweep the rubbish that is found in the gutter.As about 80% of the debris on roadways is located in the gutter [Michielen and Parker(2000); Peel, Michielen and Parker (2001)], the effective operation of the gutter brush is important.

Because of this, it is of interest to carry out research on the dynamics and performance of gutter brushes. However, the amount of research on this is very limited, and it seems to have been carried out only at the University of Surrey (UK). Research on this area has been performed by means of mathematical models [Peel (2002); Vanegas-Useche, Abdel-Wahab and Parker (2007); Vanegas-Useche, Abdel-Wahab and Parker (2008); Vanegas-Useche, Abdel-Wahab and Parker (2011a)], finite element analyses [Wang (2005);Abdel-Wahab, Parker and Wang (2007); Vanegas-Useche, Abdel-Wahab and Parker(2011b); Vanegas-Useche, Abdel-Wahab and Parker (2011c); Abdel-Wahab, Vanegas-Useche and Parker (2015)], and experimental tests [Vanegas-Useche, Abdel-Wahab and Parker (2010); Peel (2002); Abdel-Wahab, Parker and Wang (2007); Abdel-Wahab,Wang, Vanegas-Useche et al. (2011); Vanegas-Useche, Abdel-Wahab and Parker (2015a);Vanegas-Useche, Abdel-Wahab and Parker (2015b)].

Regarding finite element modeling of gutter brushes, it might involve dynamic analyses with a certain number of bristles that make and lose contact with a road surface and among them. Thus, in order to study bristle tip-road and bristle-bristle interaction, impact and, consequently, contact modeling has to be performed. When modeling impact,several parameters have to be defined; two critical parameters are the Integration Time Step, ITS, of the dynamic analysis and the normal contact stiffness, Kn, of the contact pair.The latter parameter is required by the contact algorithms pure penalty method and augmented Lagrangian method.

A suitable value of Kncannot be defined a priori; it depends on the shape, size, and type of the contact elements, as well as the material properties of the bodies. Large values of the contact stiffness may produce an ill-conditioned numerical problem, and small values may produce higher residual penetrations [Bernakiewicz and Viceconti (2002)]. The value of the normal contact stiffness for finite element analyses is usually obtained by experimentation [Hattori and Serpa (2015)].

Thus, in order to determine an appropriate range or value of Knfor the modeling of the contact between a bristle of a gutter brush and a concrete surface, this article integrates the results of experimental tests, analytical procedures, and finite element analyses involving contact and carried out in ANSYS?. Through this process, the integration time step is also determined. To the knowledge of the authors, the present work is the first article in which these two parameters are determined for the case of the finite element modeling of gutter brushes for street sweeping interacting with a surface. It is considered that this work is necessary, because such modeling is important for studying the behavior of this type of brushes.

2 Contact modeling

Impact is a complex phenomenon that takes place when two or more bodies make contact at a significant speed with each other [Riley and Sturges (1996)]. It is characterized by large forces, very short durations, and, consequently, large accelerations and decelerations [Barkan (1974)]. The contact between two surfaces is conventionally modeled in finite element analyses by means of a contact element and a target element,which constitute a contact pair. The contact pair embodies two contact “springs” that provide a tangential and a normal force when the elements are in contact; this is illustrated in Fig. 1. The tangential force, associated with the tangential contact stiffness,Kt, arises due to friction. ANSYS? defines automatically the value of Kt, which is proportional to the coefficient of friction, μ, and Kn. Regarding the normal contact force,it appears when there is a certain contact penetration, Δc, of the contact element into the target element. Loosely speaking, the normal force, N, is proportional to Δc, and the constant of proportionality is Kn. In order to achieve contact compatibility and convergence,the contact pair tends to reduce the penetration to an acceptable numerical level.

Figure 1: Modeling the contact between two surfaces in the finite element method

According to the ANSYS? documentation [SAS IP (2018)], the value of Knmay be based on convergence patterns. If the convergence difficulty is caused by too much penetration, then Knmight be underestimated, and if it is due to many equilibrium iterations for attaining convergence of displacements and residual forces, Knmay be overestimated. Whereas this recommendation may be appropriate for static problems,dynamic systems such as the one dealt with in this work are very sensitive to Kn. In fact,it is a critical parameter when modeling the impacts between bristle and road or a pair of bristles. A very high value of Knwould produce very high impact forces, decelerations,and accelerations. Conversely, a very low value of Knwould produce very small forces and large penetrations. These may or may not reflect the real contact behavior.

Indeed, the contact stiffness between two bodies depends on a number of aspects such as contact geometry, the characteristics of the surfaces (e.g. the size of asperities), the Young’s moduli of the contacting bodies, and the dimensions. At a macro level, the contact between a bristle tip and a surface could be, in theory, point-to-surface, line-tosurface, or surface-to-surface, as shown in Fig. 2. It could be argued that the contact stiffness is different in every case, as the resistance to compression varies from a minimum in the case of point-to-surface contact (Fig. 2(a)) to a maximum in the case of surface-to-surface contact (Fig. 2(c)). At a micro level, the asperities of the surfaces hinder full contact between them (i.e. the real contact area is smaller than the apparent contact area). The contact stiffness may depend on the real contact area, which in turn depends on the normal force. When the normal force is small, few asperities are in contact, and the contact stiffness tends to be small. As the force increases, more asperities become into contact; this increases the real contact area and the contact stiffness, as the resistance to deformation increases with this area. Moreover, during the impact the tip of the bristle may contact the cement (Fig. 3(a)), the aggregate (Fig. 3(b)), or both; also, it may tend to slide down an irregularity (Fig. 3(b)) if it makes contact with an inclined part of the surface. In the light of all these issues, bristle-road contact may not be characterized by a single value of Kn.

Figure 2: Schematic representation of the types of contact geometry and their contact deformations

Figure 3: Schematic representation of the contact between a bristle and either cement or aggregate of a concrete surface

Contact modeling has been performed in different manners; for example, some authors have used isogeometric formulations. Kruse et al. [Kruse, Nguyen-Thanh, De Lorenzis et al.(2015)] apply the isogeometric collocation method for frictional contact between bodies subjected to large deformation; the advantage of this formulation is that it reduces the computational cost. They validate the developed formulation through two variants of Hertz problems with friction, the classical ironing contact problem with friction, and a 3-D skewed contact patch test. It is shown that implementing the contact formulation in the collocation framework to machine precision in a 3-D problem with inclined non-matching discretization passes the contact patch test. Similarly, Kruse et al. [Kruse, Nguyen-Thanh,Wriggers et al. (2018)] present an isogeometric formulation for frictionless contact between deformable bodies. For this formulation, they apply the concept of the third medium, in which not only the contacting bodies, but also the contact medium relies on continuum formulations in which the bodies can interact. As occurs with the normal contact stiffness,Kn, the parameters of this intermediate medium are important for contact behavior, and the authors examine their role, based on numerical tests. They conclude that, as expected, the stiffness of the third body must be chosen so that it is considerably smaller than its counterparts of the interacting bodies. However, a very small value will lead to convergence problems. An advantage of the method is that the non-smooth contact problem is transformed into a smooth continuum problem. Finally, it is noted that isogeometric analyses have also been dealt with, for example, for the non-linear deformation of thin shells: for multi-patches based on RHT-splines [Nguyen-Thanh, Zhou, Zhuang et al. (2017)]and by coupling the isogeometric approach with the meshfree method [Li, Nguyen-Thanh and Zhou (2018)]; for crack propagation problems, by coupling isogeometric analysis with the meshfree method [Nguyen-Thanh, Huang and Zhou (2018)] and through an adaptive extended isogeometric formulation based on polynomial splines over hierarchical T-meshes[Nguyen-Thanh and Zhou (2017)].

3 Methodology

Experimental tests and finite element analyses were performed to estimate an appropriate value of the normal contact stiffness, Kn, for bristle tip-concrete surface interaction. The experimental tests consisted of releasing a hardened tempered mild steel bristle of rectangular cross section (0.50 mm×2.08 mm), and with a length of 336.5 mm, from a set of heights, y0, over a motorway grade concrete test surface of 1 m×1 m (Fig. 4). The height y0is the vertical distance between the lower tip and the horizontal surface. As shown in Fig. 4(a), the bristle is placed vertically, and it is released from rest (velocity v0=0, at time t0=0); its motion before and after the first impact was recorded by means of a digital camera. It has to be noted that the conditions withstand by the bristles of a gutter brush differ greatly from those of the experiments. The bristles tend to impact the surface with a certain inclination, are clamped at the upper end, and withstand bending. Therefore,the impact forces and times will vary at different fashions. Nonetheless, taking into account that the experiments and the model recreate the same case, it is expected that the results obtained in this work are suitable as a first approximation. Future work may be performed in order to recreate conditions that are similar to those in the actual brushing process.

Figure 4: Main stages in the collision between the steel bristle and the concrete surface

For the experiments, it is necessary to determine suitable heights for releasing the bristle.To obtain adequate heights, a practical range for the component normal to the road of the velocity of the bristle tips (vn) when they are about to contact the surface is estimated. In a gutter brush, this velocity will vary from an unknown maximum value, in the case of a bristle contacting the surface at high speed, to a minimum of zero. When vn=0, the exact value of Knbecomes essentially unimportant. However, when the velocity is high, Knis critical, as this will determine how fast the bristle tip will be decelerated and accelerated.Therefore, the value of Knmay be determined based on a velocity close to the maximum expected velocity. From practical values of the geometric and operating parameters of gutter brushes (see, for example Vanegas-Useche et al. [Vanegas-Useche, Abdel-Wahab and Parker (2010)]), it is estimated that an appropriate value of vncould be 0.9 m/s.Applying the principle of conservation of energy between t0(when the bristle is released)and ti(when it makes contact with the surface):

where m is the mass of the bristle, viis its velocity when contact begins, and g is the acceleration due to gravity. Thus

Therefore, releasing the bristle from a height of 41 mm will produce a velocity vi=0.9 m/s2. Based on this, the following five heights y0are selected for the experiments: 50 mm,40 mm, 30 mm, 20 mm, and 10 mm. The bristle was released 26 times from these different heights. The position of the bristle tip was tracked through a digital video camera that captures 25 frames per second. The images were analyzed frame by frame, to obtain the experimental points (y, t) (tip position, time). Then, the points were fitted with a pair of curves based on the kinematics equations for constant acceleration. Fig. 5 shows four frames of one of the experiments, and Fig. 6 presents its results. The maximum height after the impact is estimated from the curve. For the case in Fig. 6, the initial position is y0=50.5 mm and the highest position after the impact is y1=34 mm. Values for y0and y1were recorded for each experiment.

Figure 5: Estimated positions of the bristle tip at different times

Figure 6: Experimental values and theoretical trend for one of the experiments; the error bars represent the uncertainty from the fuzziness of the images

Regarding the model, it consists of a beam, in a vertical position and free of constraint,that is subjected to gravity and that collides with a rigid surface. The beam is assumed straight, homogeneous, and isotropic, with the dimensions given previously, and with the following material properties: Elastic modulus E=207 GPa, density ρ=7800 kg/m3,Poisson ration ν=0.28, static and dynamic coefficient of friction μ=0.5. For the interaction between the beam and the surface (node-to-surface contact), contact and target elements are used (Fig. 7). The CONTA175 element (contact point or surface) is attached to the lowest node of the beam. The TARGE170 element (3-D rigid or flexible target surface) is attached to a flat horizontal area that corresponds to the surface and is modeled as rigid. In this transient analysis, gravity is modeled by applying an acceleration (ACEL command) of g=9.8 m/s2. The impact is assumed to be perfectly elastic and is solved through the augmented Lagrangian method.

Figure 7: Finite element model of the collision between a beam and a surface

The model was applied with the following characteristics: stiffness proportional damping coefficient βD=0.0107 ms, obtained experimentally [Vanegas-Useche, Abdel-Wahab and Parker (2015a)], mass proportional damping coefficient αD=0 (air resistance is neglected),16 BEAM189 elements. The element type BEAM189 is a 3-D quadratic (3-node) finite strain beam. It is limited by two end nodes and has a midside node. The element has an optional node, which is used to indicate the orientation of the cross section of the beam. It is based on first-order shear-deformation theory (Timoshenko beam theory) and is appropriate for analyzing slender to moderately stubby beams. It has 6 or 7 degrees of freedom at each node (three translations, three rotations, and an optional warping magnitude). This element is suitable for nonlinear, large rotation analysis and includes stress stiffness terms. The values of Knstudied include the set {0.01, 0.1, 0.5, 0.75, 1, 2, 3,4, 5, 10, 100} MN/m. The total simulation time was varied from 0.13 s to 0.3 s,depending on the initial height of the beam, in order to ensure that the position of the lower tip is tracked before and after the first impact. The required sensitivity analyses were carried out. Analyses were made to determine the effect of the number of beam elements and the number of target elements. It was found that the number of beam elements required to produce errors less than 1% is much less than 16; however, it was decided to use 16 elements, because this number was the one used in other works (e.g.[Vanegas-Useche, Abdel-Wahab and Parker (2011c)]) and did not increase the computing time considerably. Regarding the number of target elements, this mesh does not affect the results, as these elements are assumed rigid.

Finally, it is necessary to determine an appropriate substep time, or Integration Time Step(ITS), for modeling the short duration impact; the ITS has to be small enough to appropriately model this contact. It is not suitable to let ANSYS? automatically choose the ITS, because unreliable results are obtained. Therefore, it is necessary to specify an upper limit for it. To estimate an appropriate value for this limit, an analysis is carried out based on the case shown in Fig. 7. In ANSYS?, contact is established when there is a certain amount of contact penetration, Δc(Fig. 8). First, the maximum penetration, Δcmax,is estimated based on the principle of conservation of energy. At time t0(Fig. 4(a)), the bristle has a potential energy that may be approximated to m g y0, where mg is the weight of the bristle and y0is the initial height; the reference for this energy corresponds to position at time tp(Fig. 4(c)), which is the time at which the maximum contact penetration is achieved. At this time, there is a potential energy (spring with constant Kn) given by 0.5 KnΔcmax2. Equating these two energies produces:

The penetration as a function of time may be obtained from Newton’s second law. Fig. 8 shows the free body diagram of the bristle during the impact. The bristle is subjected to the weight mg and the contact force KnΔc; therefore:

where ∑F corresponds to the forces withstood by the bristle and a is its acceleration.

Taking into account that the impact force (KnΔc) tends to be much greater than the weight of the bristle:

Given that Δc(ti)=0 and Δc’(ti)=vi=(2 g y0)1/2(from Eq. (1)), the solution of this homogeneous equation is:

where t is time and tiand tfare the time at the start and end of contact, respectively (Figs.4(b) and 4(d)).

Figure 8: Free body diagram of the bristle during impact (time range [ti, tf])

Finally, the time required for the bristle to stop from the time in which contact starts is obtained from the principle of linear impulse and momentum between times tiand tp:

Given that vp=0, neglecting again mg, and using Eq. (6), Eq. (7) produces:

As an example, Fig. 9 shows curves of contact penetration vs. contact time for three initial heights, y0, and Kn=2 MN/m. From Eq. (8) or Fig. (9), it is concluded that tp– ti(for the case illustrated in Fig. 7) is independent of y0, i.e., it is independent of the velocity at the start of contact. Although the case of the bristle of a gutter brush differs from the problem in Fig. 7, the contact time may be of the order of the value given by Eq.8. Therefore, the maximum ITS can be a constant value based on this equation, regardless of the dynamic characteristics of the bristles.

Figure 9: Contact penetration against contact time for three values of y0; m=2.73 g, g=9.8 m/s2, Kn=2 MN/m

4 Results and discussion

Fig. 10 presents the results of the experiments: initial position of the bristle, y0, versus its highest position after impact, y1. The high dispersion exhibited may be due to effects related to (a) the roughness and irregularities of the concrete surface, as well as the bristle tip surface, (b) the fact that concrete is a composite material whose constituents have different properties, and (c) the different contact geometries. Firstly, experimental observations suggest that when bristle-surface contact occurs in the inclined surface of an irregularity, the maximum height after the impact tends to be smaller. This may occur because the tip tends to slide down the irregularity, and the work of the friction force reduces the kinetic energy of the bristle. Secondly, the concrete surface may be considered as a three-phase composite material that consists of cement, aggregates, and interfacial transition zone [Zheng and Zhou (2006)]. The dynamics of the impact differ,depending on whether the tip contacts an aggregate, the cement paste, or a transition zone.This is because they have different characteristics, in particular, different Young’s moduli(Poisson ratio, hardness, and wear resistance are other characteristics that affect impact dynamics). E.g., it is reported that a particular cement paste, a fine aggregate, and a coarse aggregate have elastic moduli of 12 GPa, 80 GPa, and 69 GPa, respectively[Zheng and Zhou (2006)]. Lastly, impact behavior tends to vary because each impact is different from each other at a macro (Fig. 2) or micro (asperity) level (Fig. 3).

Figure 10: Initial position against highest position after impact of a steel bristle in free fall impacting a concrete surface

The experimental data in Fig. 10 were fitted through different regression functions: linear,quadratic, logarithmic, exponential, and power regressions. The coefficients of determination, R2, for these functions vary between 0.6502 and 0.7394. The highest value of R2(0.7394) is obtained for the quadratic regression; however, the linear regression yields R2=0.7326, which is very close to the highest value. Additionally, taking into account that the trend should intercept the origin of the diagram, both trends (linear and quadratic) were enforced to intersect the point (0, 0). These regressions are provided in Figs. 11 and 12 as solid lines. In these figures, it may be observed that R2for the quadratic function is slightly smaller than the initial value and that R2for the linear regression does not change (the initial trend intercepts the y1axis at y1=0.006 mm≈0). As both trends are very similar and yield almost the same R2value, it is decided to model the experimental data through the linear equation. Therefore, y1is approximately proportional to y0, i.e. the ratio y1/ y0is constant and approximately equal to 65%:

Nevertheless, the results indicate that the impact behavior vary widely. Therefore, the experimental ratios y1/ y0were fitted to a normal distribution; this is characterized by a mean value of 65% and a standard deviation of 21%. The range of values y1/ y0for a confidence of 90% corresponds to [30%, 100%], which in turn corresponds to 65%±1.645 times the standard deviation. Thus, it may be stated that y1/ y0is mostly between 30% and 100% (dashed lines in Fig. 12).

Figure 11: Quadratic regression for the experimental points (y0, y1)

Figure 12: Linear regression for the experimental points (y0, y1)

Regarding the integration time step, the times required to achieve the maximum penetration are, from Eq. (8), tp–ti≈0.82 ms and 0.01 ms, for Kn=0.01 MN/m and 100 MN/m, respectively; these two values are the minimum and maximum values of Knused.Therefore, the ITS should be about equal or smaller than the corresponding value. A number of sensitivity analyses were carried out, and it was found that a maximum ITS of 0.01 ms is appropriate, as smaller upper limits provide practically the same results. It is decided that the upper limit of the ITS is 0.01 ms.

Finally, the results of the simulations are presented in Fig. 13; this shows curves y vs. t for the different initial heights and diverse values of Kn. For all the values of y0, the value of Knthat produces y1=0.65·y0is 2.0 MN/m. Therefore, it could be assumed that an appropriate mean value of Knis 2 MN/m, regardless of the initial velocity at the beginning of the impact, which depends on the initial height. However, there is not a unique appropriate value of Kn, due to the high variability of the impact situations. The normal contact stiffness could mostly vary between about 100 MN/m (y1≈0.3·y0) and 0.01 MN/m (y1≈y0). These two extreme values are different by a factor of 1×104. This high variability of Knis partly due to the reasons already discussed, when analyzing the high dispersion of the experimental results. The implications of this high variability are that the modeling of the bristle as a beam (1-D) element together with the concrete as a flat surface produces approximate results. More accurate models would require that the bristle is modeled as a 3-D element, and the concrete bed as a rough composite body.However, as the distribution, size, and orientation of aggregates, as well as the surface roughness of concrete, vary from one road to another, it may be difficult to obtain accurate solutions through such a complex model. Therefore, it is considered that modeling the road as a flat surface is a suitable approximation; the value of Knmay be varied in the range 0.01 MN/m a 100 MN/m to capture different impact scenarios.

Figure 13: Vertical position of the bristle tip against time for different values of Kn

5 Conclusions

This article presented the methodology and results of experimental tests, as well as the accompanying finite element analyses, to determine the normal contact stiffness for the interaction between a steel bristle of a gutter brush for road sweeping and a concrete surface. Also, an appropriate maximum value of the integration time step was determined,based on the principle of conservation of energy, Newton’s second law, the principle of linear impulse and momentum, and various sensitivity analyses.

The normal contact stiffness, which is a parameter of a contact pair in a finite element model, is a critical parameter in problems involving impacts. A large value of the normal contact stiffness may produce large impact forces, and a small value may produce small forces and large penetrations, which may not reproduce the actual impact behavior. The comparison between the results of experimental tests and finite element analyses suggests that the normal contact stiffness to model the interaction between the steel bristles and a concrete surface is independent of the initial velocity before impact (in the practical velocity range studied) and may be of about 2 MN/m. However, it could vary notably, e.g.between 100 MN/m and 0.01 MN/m, for a 90% confidence, due to the inherent variability of the characteristics of the impact between a bristle tip and a concrete surface.Also, the results indicate that, for the case studied, a constant value of Knmay be related to a constant ratio y1/ y0; besides, the impact time tends to be constant, as expected from the principle of linear impulse and momentum. As for the integration time step, it was found that 0.01 ms is an appropriate maximum value. As a further work, more accurate models will require: (a) Experimental tests that simulates the real conditions of the interaction between a gutter brush and a concrete surface, in order to recreate more accurately the impact conditions and (b) Modeling of the bristle as a solid body and the concrete as a body with a rough surface composed of the different aggregates.

Acknowledgements:The authors acknowledge the support of the Universidad Tecnológica de Pereira (Colombia), the University of Surrey (UK), Ghent University (Belgium), and the Programme Alβan, European Union Programme of High Level Scholarships for Latin America, identification number (E03D04976CO).