# Validation of Friction Model Parameters Identified Using the IHB Method Using Finite Element Method

Authors: Abdallah Hadji; Njuki Mureithi
Journal: Shock and Vibration
 (2019)

---

## Abstract

A hybrid friction model was recently developed by Azizian and Mureithi (2013) to simulate the friction behavior of tube-support interaction. However, identification and validation of the model parameters remains unresolved. In previous work, the friction model parameters were identified using the reverse harmonic method, where the following quantities were indirectly obtained by measuring the vibration response of a beam: friction force, sliding speed of the force of impact, and local displacement at the contact point. In the present work, the numerical simulation by the finite element method (FEM) of a beam clamped at one end and simply supported with the consideration of friction effect at the other is conducted. This beam is used to validate the inverse harmonic balance method and the parameters of the friction models identified previously. Two static friction models (the Coulomb model and Stribeck model) are tested. The two models produce friction forces of the correct order of magnitude compared to the friction force calculated using the inverse harmonic balance method. However, the models cannot accurately reproduce the beam response; the Stribeck friction model is shown to give the response closest to experiments. The results demonstrate some of the challenges associated with accurate friction model parameter identification using the inverse harmonic balance method. The present work is an intermediate step toward identification of the hybrid friction model parameters and, longer-term, improved analysis of tube-support dynamic behavior under the influence of friction.


---

## Body

## 1. Introduction

The friction model is an essential element in the detailed analysis of the dynamics of steam generator tubes in the nuclear industry [1]. Most of the friction models currently used to simulate tube-support interaction are cited in [2]. These are special cases of static friction models. The velocity-limited friction model (VLFM) [3] is a continuous Coulomb model without Stribeck effects [4]. The force balance friction model (FBFM) [5] and the spring damper friction model (SDFM) [6] are two models based on springs and dampers. The principle of the Karnopp friction model [7] is used in the FBFM friction model. Recently, a hybrid friction model [1, 8] has been developed to model tube-support interaction. In the hybrid model, all the properties and benefits of other friction models are included, namely, the dynamics of the Dahl model [9] (i.e., hysteresis effect), the dynamics of bristles [10], and the Stribeck effect [4] (transition from the static friction limitation to the kinetic friction limitation). This model also indirectly takes into account the distribution of the stresses in the contact area according to the principle of Cattaneo–Mindlin [11] to model the presliding phenomenon.In finite element modeling (FEM) codes, the Coulomb friction model and the decay friction model are widely used; these are combined with the principle of the velocity limit [3] to model the presliding (sticking) regime [12]. Most of the numerical algorithms for analyzing dynamic friction in FEM codes are presented by Oden and Martins [13]. Diehl [14] has also numerically investigated the friction effect of a circular rigid body in sliding contact with a flexible beam.In the present work, we validate the friction coefficient identified experimentally in previous work [15] using the finite element method (FEM) using the Abaqus software.Four standard friction model parameters have been identified by direct friction force measurement. In the present ongoing work an indirect approach, based on acceleration measurement, has been proposed [15–17]. To identify the friction model parameters experimentally using this indirect method, accurate nonlinear normal modes are needed. Nonlinear normal modes (NNMs) and the principle of inverse harmonic balance (IHB) based on the harmonic balance (HB) method [18] were developed and presented in previous work [15–17].In the previous work [15–17], the parameters of the Dahl [9] and LuGre [19] friction models, which are, respectively, based on Coulomb and Stribeck friction models were reported. In this paper, therefore, we analyze these two friction models to validate their friction coefficients identified also in the work [15], where the basic friction models identified parameters have been used to identify the parameters of the more complex models.


## 2. Principle of the Inverse Harmonic Balance (IHB) Method

Figure1 shows a schematic of the test rig used to extract the contact force (friction and impact) and displacement. The specimen is a beam simply supported at one end (but allowing sliding hence friction effects) and clamped at the other. In Figure 2 all the forces acting at the contact point are represented; for more details see [15–17, 20].Figure 1
Schematic of a nonlinear beam [17].Figure 2
Forces acting at the contact point [17].The beam equation of motion is(1)ρAy¨x,t+EIy″″x,t−Tty″x,t=Ftδx−Xf,with the following boundary conditions:(2)x=0:y0,t=0,y′0,t=0,x=L:EIy″L,t=r+e2Tt,EIy‴L,t=−Nt+my¨L,t.The characteristics of the beam are presented in Table1 below. The beam response is measured by six accelerometers, where their type and position are cited in the previous works [15–17].Table 1
Finite element model characteristics.
Property
Symbol
Quantity
Units
Young’s modulus
E
210
GPa
Density
ρ
7800
kg/m3
Length
L
606
mm
Width
l
38.88
mm
Thick
e
4.78
mm
Radius
r
5
mm
Damping ratio
ζ
1
0.01
%
Proportional mass damping
α
1
0.06
Normal force
N
83
N
Excitation force amplitude
F
ex
6.5
NBased on modal superposition, the response of the system may be written in the following form:(3)yjxj,t=∑i=1nϕixjqit,where qiti=1,2,…,n are the generalized coordinates, ϕixi=1,2,…,n are the normal modes of the system.Using the modal superposition principle, equation (3), the boundary conditions equations, equation (2), and the Galerkin method, equation (1) can be rewritten as follows:(4)q¨it+ωi2qit−FtϕiXf=−NtϕiL+rϕi′L+∑j=1nqjt∫0Lϕj″xϕ˜ixdx−ϕi′Lϕ˜jLTt,i=1−n,where Nt=P−RtModal orthogonality may be expressed as follows:(5)ρA∫0Lϕix,aϕjx,adx+mϕiL,aϕjL,a=δij,EI∫0Lϕì″x,aϕj″x,adx=ωi2δij,and m is the half cylinder mass (contact element).In our previous work, different methods were proposed to identify the nonlinear normal modes (NNMs). Here, we present only a summary of the most accurate method, where the NNMs are computed based on the coupled subharmonic approach. The NNMs are calculated in the following manner:Firstly, the Fourier series of the response signal is written as follows:(6)yjxj,t=∑r=1nαrjxjcosr⋅ωextt+∑r=1nβrjxjsinr⋅ωextt.Next this equation is written as(7)yjxj,t=∑r=1nϕ˜˜ir1xjγr1cosr⋅ωextt+∑r=1nϕ˜˜ir2xjγr2sinr⋅ωextt,(8)q¨˜˜ir1t=γr1cosr⋅ωextt,q¨˜˜ir2t=γr2sinr⋅ωextt,(9)ϕ˜˜ir1xj=αrjxjγr1,ϕ˜˜ir2xj=βrjxjγr2,where ϕ˜˜ir1 is the first subharmonic form (or cosines subharmonic form) and ϕ˜˜ir2 is the second subharmonic form (or sine subharmonic form). The modal equation of motion (equation (4)) becomes(10)∑s=1narskq¨˜˜iskt+arslq¨˜˜islt+brskq˜˜ikt+brslq˜˜ilt−Ftϕ˜˜irkXf=−Ntϕ˜˜irkL+∑s=1nq˜˜is1t∫0Lϕ˜˜is1″xϕ˜irkxdx+r+e/2ϕ˜˜irk′L+∑s=1nq˜˜is2t∫0Lϕ˜is2″xϕ˜irkxdx−∑s=1nϕ˜˜is1′Lq˜˜is1t+∑s=1nϕ˜˜is2′Lq˜˜is2tϕ˜˜irkLTtr=1−n,where k and l are the subharmonic form indices; k represents the subharmonic form index used to normalize the modal equation of motion (equation (10)) (as the test function of Galerkin method [21]), and l is the second subharmonic form. The values of k and l can be chosen as follows: k=1 (first subharmonic form), l=2 (second subharmonic form) or k=2, l=1. In the present work, k=1 and l=2 were chosen. The modal orthogonality equation (5) becomes(11)ρA∫0Lϕ˜˜iskx,aϕ˜˜irkx,adx+mϕ˜˜irkL,aϕ˜˜iskL,a=arsk,EI∫0Lϕ˜˜ìrk″x,aϕ˜˜isk″x,adx=brsk,(12)ρA∫0Lϕ˜˜iskx,aϕ˜˜irlx,adx+mϕ˜˜irkL,aϕ˜˜islL,a=arsl,EI∫0Lϕ˜˜ìrk″x,aϕ˜˜isl″x,adx=brsl,where(13)arrk=δrs,arrl=δrs.Equations (8) and (9) are used to calculate γr1 and γr2 by the normalization of ϕirk from αrj and βrj. Equation (13) is used in the harmonic forms normalization.Finally, the nonlinear normal modesϕ˜˜ir1xandϕ˜˜ir2x can be represented by a natural smoothing spline or a series of trigonometric functions (equation (14)):(14)ϕ˜˜irx=∑k=1naik1sinβikx+aik2sinhβikx+aik3cosβikx+aik4coshβikx,where n is chosen equal to 3 to minimize the fitting error.


## 3. Friction Models

Coulomb and decay friction models are widely used in finite element (FEM) software. In previous work [22], we analyzed these two friction models and validated the friction coefficient identified experimentally in the works [16, 17] using a 1D element. The results showed that the 1D element is not capable of modeling the contact problem with friction. Following the improvement of the inverse harmonic balance method [15], five new parameters of the hybrid friction model were reidentified using the newly calculated experimental results. In the present work, we analyze the different implementations of these models to validate the friction coefficient identified experimentally in [15] using a plane stress element. Firstly, we present the principles and formulation of FEM friction models.


### 3.1. Coulomb Friction Model

In the Coulomb friction model, the friction force is computed in an analytical formulation:(15)T=μcN−sigv,ifv≠0,RT,ifv=0,where μc is the kinetic friction coefficient, N is the resultant of the normal forces at the contact point, v the velocity at the contact point, and RT is the resultant of the tangential force at the contact point in the sliding direction.Generally, in the sticking regime, the resultant of the tangential forceRT is less than the kinetic friction force μcN. The sliding regime begins when the resultant of the tangential force reaches and surpasses this limit (RT≥μcN). This formulation is computed using the Lagrange method in the FEM software Abaqus [12, 23]. The formulation of the Coulomb friction model using the penalty method in Abaqus is similar in form to the Karnopp friction model [7] or the velocity-limited friction model [3]. The latter is used to model the sticking (or presliding) regime (called the “elastic slip” regime in Abaqus [12, 24]). The Coulomb friction model with a FEM penalty formulation is given by(16)T=ksu,ifu≤ucrit,μcN−sigv,ifu>ucrit,(17)ks=Tcritucrit,(18)Tcrit=μcN,where Tcrit is the critical friction force and ucrit the critical elastic slip.In general,ks or ucrit is defined by the user in most FEM software. However, in Abaqus, there are two options to define the limiting values: in the first, ks or ucrit are defined by the user (as above), while in the second, automatic option, ucrit is calculated during the simulations using the formula:(19)ucrit=Ffl¯i,where l¯i is the characteristic contact surface length and Ff is the slip tolerance; the default value, Ff=0.005 [12].The slip toleranceFf is also defined as the maximum allowable “elastic slip” (presliding displacement), expressed as a fraction of a characteristic length. The characteristic length is the average length of the contact surface elements (half cylinder). A value of 0.5% is considered typical and used by default in the code. Alternatively, an absolute value (e.g., based on known experimental values) of the characteristic length may be specified. Note that the “elastic slip” is, strictly, not slip at all but rather refers to a relative presliding displacement between the contact surfaces. The relative surface displacement is possible due the contact area asperities which deform elastically without breaking. Gross slip occurs once all the asperities are broken.


### 3.2. Decay Friction Model (Stribeck Friction Model)

The decay friction model is a special case of the Stribeck friction model [4] (equation (20)), with the exponent δ is equal to one. Equation (20) represents the general Stribeck friction model [4] in the slip regime. However, in the presliding (or elastic slip) regime, the friction force takes the same value in equations (16)–(18) with the replacement of μc by μs in these equations (16) and (18):(20)T=Nμc+μs−μce−v/νsδ−sigv,where μs is the static friction coefficient, vs is the Stribeck velocity, and δ is the Stribeck exponent.To validate the friction model parameters using a FEM code, we must validate the other parameters used as well as the contact algorithms (surface to surface or nodes to surface) and the slip toleranceFf. In the present work, the surface-to-surface contact algorithm with the hard contact model is considered. One of the important parameters investigated is the slip tolerance Ff.


### 3.3. Hybrid Friction Model

This paper presents intermediate results in a longer-term project on hybrid friction model parameter identification and validation using the simple friction models used to create the general hybrid friction model.In the previous work [15], five parameters of the hybrid friction model were identified using Dahl [9] and LuGre friction models [19]. Some of these parameters are the same as those of the Coulomb and decay (Stribeck) friction model. Starting then with these simpler models, some of the challenges associated with FEM-based parameter identification have been highlighted. Importantly, it is clear that friction models, as implemented in commercial FEM codes, should be used with caution. When quantities intimately related to the details of the friction model are considered, the resulting physical outputs may be significantly affected and far from true values. For instance, the slip tolerance is an adjustable parameter that is found in the numerical friction models. This parameter has, however, been found to have little effect on the computed friction force and displacement when physically realistic values are used.The hybrid friction model [1, 8] (Figure 3) is based on the principle of the LuGre model [19]. All properties and benefits of other friction models are included, namely, the dynamics of the Dahl model [9] (i.e., hysteresis effect), the dynamics of bristles [10], and the Stribeck effect [4]. This model also takes into account the distribution of the stresses in the contact area according to the principle of Cattaneo–Mindlin [11] to model the presliding phenomenon.Figure 3
Schematic of the hybrid friction model [8].The formulation of the hybrid friction model [1] is represented in the following equations:(21)T=keze+cpz˙p+csz˙s,keze+kepze−zp=0,kepzp−ze+kpszp−zs+cpz˙p=0,kpszp−zs+czz˙s−x˙+csz˙s=0,where ze, zp, and zs are, respectively, the elastic, plastic, and partial slip relative displacements. The presliding or sticking regime is modeled by stiffness ke, plastic damping cp and presliding damping cs, and two transition stiffness coefficients, elastic-plastic kep and plastic presliding kps, and the slip regime is modeled by the coefficients of the Stribeck damping cz which is given by(22)cz=σ0gx˙x˙,(23)σ0gx˙=Nμc+μs−μce−x˙/νs2,where gx˙ is the Stribeck function (or Stribeck friction model [4]), μc is the kinetic friction coefficient, μs is the static friction coefficient, and vs is the Stribeck velocity.


## 3.1. Coulomb Friction Model

In the Coulomb friction model, the friction force is computed in an analytical formulation:(15)T=μcN−sigv,ifv≠0,RT,ifv=0,where μc is the kinetic friction coefficient, N is the resultant of the normal forces at the contact point, v the velocity at the contact point, and RT is the resultant of the tangential force at the contact point in the sliding direction.Generally, in the sticking regime, the resultant of the tangential forceRT is less than the kinetic friction force μcN. The sliding regime begins when the resultant of the tangential force reaches and surpasses this limit (RT≥μcN). This formulation is computed using the Lagrange method in the FEM software Abaqus [12, 23]. The formulation of the Coulomb friction model using the penalty method in Abaqus is similar in form to the Karnopp friction model [7] or the velocity-limited friction model [3]. The latter is used to model the sticking (or presliding) regime (called the “elastic slip” regime in Abaqus [12, 24]). The Coulomb friction model with a FEM penalty formulation is given by(16)T=ksu,ifu≤ucrit,μcN−sigv,ifu>ucrit,(17)ks=Tcritucrit,(18)Tcrit=μcN,where Tcrit is the critical friction force and ucrit the critical elastic slip.In general,ks or ucrit is defined by the user in most FEM software. However, in Abaqus, there are two options to define the limiting values: in the first, ks or ucrit are defined by the user (as above), while in the second, automatic option, ucrit is calculated during the simulations using the formula:(19)ucrit=Ffl¯i,where l¯i is the characteristic contact surface length and Ff is the slip tolerance; the default value, Ff=0.005 [12].The slip toleranceFf is also defined as the maximum allowable “elastic slip” (presliding displacement), expressed as a fraction of a characteristic length. The characteristic length is the average length of the contact surface elements (half cylinder). A value of 0.5% is considered typical and used by default in the code. Alternatively, an absolute value (e.g., based on known experimental values) of the characteristic length may be specified. Note that the “elastic slip” is, strictly, not slip at all but rather refers to a relative presliding displacement between the contact surfaces. The relative surface displacement is possible due the contact area asperities which deform elastically without breaking. Gross slip occurs once all the asperities are broken.


## 3.2. Decay Friction Model (Stribeck Friction Model)

The decay friction model is a special case of the Stribeck friction model [4] (equation (20)), with the exponent δ is equal to one. Equation (20) represents the general Stribeck friction model [4] in the slip regime. However, in the presliding (or elastic slip) regime, the friction force takes the same value in equations (16)–(18) with the replacement of μc by μs in these equations (16) and (18):(20)T=Nμc+μs−μce−v/νsδ−sigv,where μs is the static friction coefficient, vs is the Stribeck velocity, and δ is the Stribeck exponent.To validate the friction model parameters using a FEM code, we must validate the other parameters used as well as the contact algorithms (surface to surface or nodes to surface) and the slip toleranceFf. In the present work, the surface-to-surface contact algorithm with the hard contact model is considered. One of the important parameters investigated is the slip tolerance Ff.


## 3.3. Hybrid Friction Model

This paper presents intermediate results in a longer-term project on hybrid friction model parameter identification and validation using the simple friction models used to create the general hybrid friction model.In the previous work [15], five parameters of the hybrid friction model were identified using Dahl [9] and LuGre friction models [19]. Some of these parameters are the same as those of the Coulomb and decay (Stribeck) friction model. Starting then with these simpler models, some of the challenges associated with FEM-based parameter identification have been highlighted. Importantly, it is clear that friction models, as implemented in commercial FEM codes, should be used with caution. When quantities intimately related to the details of the friction model are considered, the resulting physical outputs may be significantly affected and far from true values. For instance, the slip tolerance is an adjustable parameter that is found in the numerical friction models. This parameter has, however, been found to have little effect on the computed friction force and displacement when physically realistic values are used.The hybrid friction model [1, 8] (Figure 3) is based on the principle of the LuGre model [19]. All properties and benefits of other friction models are included, namely, the dynamics of the Dahl model [9] (i.e., hysteresis effect), the dynamics of bristles [10], and the Stribeck effect [4]. This model also takes into account the distribution of the stresses in the contact area according to the principle of Cattaneo–Mindlin [11] to model the presliding phenomenon.Figure 3
Schematic of the hybrid friction model [8].The formulation of the hybrid friction model [1] is represented in the following equations:(21)T=keze+cpz˙p+csz˙s,keze+kepze−zp=0,kepzp−ze+kpszp−zs+cpz˙p=0,kpszp−zs+czz˙s−x˙+csz˙s=0,where ze, zp, and zs are, respectively, the elastic, plastic, and partial slip relative displacements. The presliding or sticking regime is modeled by stiffness ke, plastic damping cp and presliding damping cs, and two transition stiffness coefficients, elastic-plastic kep and plastic presliding kps, and the slip regime is modeled by the coefficients of the Stribeck damping cz which is given by(22)cz=σ0gx˙x˙,(23)σ0gx˙=Nμc+μs−μce−x˙/νs2,where gx˙ is the Stribeck function (or Stribeck friction model [4]), μc is the kinetic friction coefficient, μs is the static friction coefficient, and vs is the Stribeck velocity.


## 4. Numerical Simulation

### 4.1. Finite Element Modeling

In previous work [22], a 1D FEM model (Figure 4(a)) was used to validate the friction coefficient identified experimentally in the works [16, 17]. This model was also compared and validated using the 2D and 3D models. However, in this work, the beam of Figure 1 is modeled by 2D elements FEM model (Figure 4(b)). Then, in conjunction with the friction models presented above, it is used to validate the friction coefficient identified experimentally in the previous work [15]. The beam is discretized using 2D plane strain elements. The validation of this FEM model is done using a multistep mesh analysis. This includes (i) cantilever beam end static deflection calculation, (ii) modal analysis, and (iii) model size optimization. However, only the modal analysis and the model size (calculation time) were used in the choice of the meshing method and the element size of the half-cylinder (contact element). The results of this analysis are presented in Figures 5 and 6. The description of the elements presented in legend of those figures is defined in Table 2.Figure 4
Beam meshing: (a) 1D model; (b) 2D model.
(a)
(b)Figure 5
Cantilever beam end deflection (b); zoom-in of the superimposed lines in (a).
(a)
(b)Figure 6
Computed 1st natural frequency (b); is a zoom-in of the superimposed lines in (a).
(a)
(b)Table 2
Abaqus plane strain elements type.
Element symbol
Descriptions
CPS4
4-node bilinear quadrilateral element
CPS4R
4-node bilinear quadrilateral element with reduced integration
CPS4I
4-node bilinear quadrilateral with incompatible modes
CPS8
8-node biquadratic quadrilateral element
CPS8R
8-node quadrilateral biquadratic element with reduced integration
CPS3
3-node triangle bilinear
CPS6
6-node triangle biquadratic
CPS6M
6-node modified biquadratic triangle elementIn general, the triangular elements are the most commonly used to mesh the cylindrical form because of geometrical compatibility. On the other hand, we have solved the problem of the quadrilateral elements’ incompatibility for this geometry type by dividing the circular section of the cylinder into four parts. This means that our contact element (half cylinder) should be divided in two parts (two quarter cylinders) (Figure4(b)). The other contacting surface is a rigid surface for both models.


### 4.2. Finite Element Model Validation

These analyses are necessary in this work to optimize our model mesh to reduce the simulation time. For example, using our optimal mesh (Figure4(b)), the simulation for each iteration (one frequency excitation of 1.2 seconds using direct implicit dynamic integration) takes 90 minutes using a PC having an Intel i7 3 GHz CPU and 24 GB RAM. To obtain one case of parameters analysis results this translates to 30 hours.Figures5 and 6 show that the quadrilateral elements give better results than triangular elements, in both objectives cases analysis (deflection the beam free end, and a modal analysis) using the cantilever beam. The CPS4I (4-node bilinear quadrilateral with incompatible modes) element is the best among all the elements for both cases. For this reason, the element type was chosen to model the beam and the half cylinder (the contact element). The triangular element form is the most commonly used to model the cylindrical structure. In the half cylinder element (contact element) meshing analysis, the solutions adopted to obtain the best results are to divide the half cylinder in two parts making it possible to use the CPS4I element for this geometry, and maintain a space ratio less than two to optimize the number of elements with bias mesh controller concentrated in the contact zone.The characteristics of the beam modeled by FEM are presented in Table1. The friction model parameters are presented in Table 3. The parameters μc, μs, and vs in equation (23) were identified experimentally using the Dahl friction model [9] and LuGre friction model [19]. These models are detailed in the previous works [15–17]. The objective function previously used for parameter identification is that given in equation (24). Optimization based on an objective function is used to validate the simulation results here as well. The objective function is more general than the more commonly used simpler form employed previously in [17]. This objective function is more useful to compare signals with different behaviors such as the appearance of the beating phenomenon in one response signal or of signals having different offsets (difference between the min/max values).(24)Obj=maxTexp.−minTexp.−maxTsuml.−minTsuml.2,where Obj is the objective function, Texp. the experimental friction force, and Tsuml. the simulated friction force.Table 3
Friction models parameters (6.5 N) [15].
Parameter
Coulomb friction model
Dahl friction model
Stribeck friction model
LuGre friction model
μ
c
0.520
0.500
0.161
0.160
μ
s
—
—
0.477
0.492
v
s (m/s)
—
—
6.101 10−4
7.152 10−4
u
crit (m)
—
—
4.146 10−6
4.060 10−6
σ
0 (N/m)
—
1.003 10+7
—
1.140 10+4
σ
1 (N·s/m)
—
—
—
347.030
σ
2 (N·s/m)
—
—
—
0.000
δ
2All the results are obtained at the excitation force amplitude of 6.5 N. A proportional (Rayleigh) damping is used to model the beam structural damping. The mass proportional damping coefficientαn is given by the following equation:(25)αn=2ωnζn,where ωn is the angular natural frequency of mode n and ζn the corresponding damping ratio.The value of the damping ratioζ in Table 1 is identified experimentally using an impact test and the logarithmic decrement method.


## 4.1. Finite Element Modeling

In previous work [22], a 1D FEM model (Figure 4(a)) was used to validate the friction coefficient identified experimentally in the works [16, 17]. This model was also compared and validated using the 2D and 3D models. However, in this work, the beam of Figure 1 is modeled by 2D elements FEM model (Figure 4(b)). Then, in conjunction with the friction models presented above, it is used to validate the friction coefficient identified experimentally in the previous work [15]. The beam is discretized using 2D plane strain elements. The validation of this FEM model is done using a multistep mesh analysis. This includes (i) cantilever beam end static deflection calculation, (ii) modal analysis, and (iii) model size optimization. However, only the modal analysis and the model size (calculation time) were used in the choice of the meshing method and the element size of the half-cylinder (contact element). The results of this analysis are presented in Figures 5 and 6. The description of the elements presented in legend of those figures is defined in Table 2.Figure 4
Beam meshing: (a) 1D model; (b) 2D model.
(a)
(b)Figure 5
Cantilever beam end deflection (b); zoom-in of the superimposed lines in (a).
(a)
(b)Figure 6
Computed 1st natural frequency (b); is a zoom-in of the superimposed lines in (a).
(a)
(b)Table 2
Abaqus plane strain elements type.
Element symbol
Descriptions
CPS4
4-node bilinear quadrilateral element
CPS4R
4-node bilinear quadrilateral element with reduced integration
CPS4I
4-node bilinear quadrilateral with incompatible modes
CPS8
8-node biquadratic quadrilateral element
CPS8R
8-node quadrilateral biquadratic element with reduced integration
CPS3
3-node triangle bilinear
CPS6
6-node triangle biquadratic
CPS6M
6-node modified biquadratic triangle elementIn general, the triangular elements are the most commonly used to mesh the cylindrical form because of geometrical compatibility. On the other hand, we have solved the problem of the quadrilateral elements’ incompatibility for this geometry type by dividing the circular section of the cylinder into four parts. This means that our contact element (half cylinder) should be divided in two parts (two quarter cylinders) (Figure4(b)). The other contacting surface is a rigid surface for both models.


## 4.2. Finite Element Model Validation

These analyses are necessary in this work to optimize our model mesh to reduce the simulation time. For example, using our optimal mesh (Figure4(b)), the simulation for each iteration (one frequency excitation of 1.2 seconds using direct implicit dynamic integration) takes 90 minutes using a PC having an Intel i7 3 GHz CPU and 24 GB RAM. To obtain one case of parameters analysis results this translates to 30 hours.Figures5 and 6 show that the quadrilateral elements give better results than triangular elements, in both objectives cases analysis (deflection the beam free end, and a modal analysis) using the cantilever beam. The CPS4I (4-node bilinear quadrilateral with incompatible modes) element is the best among all the elements for both cases. For this reason, the element type was chosen to model the beam and the half cylinder (the contact element). The triangular element form is the most commonly used to model the cylindrical structure. In the half cylinder element (contact element) meshing analysis, the solutions adopted to obtain the best results are to divide the half cylinder in two parts making it possible to use the CPS4I element for this geometry, and maintain a space ratio less than two to optimize the number of elements with bias mesh controller concentrated in the contact zone.The characteristics of the beam modeled by FEM are presented in Table1. The friction model parameters are presented in Table 3. The parameters μc, μs, and vs in equation (23) were identified experimentally using the Dahl friction model [9] and LuGre friction model [19]. These models are detailed in the previous works [15–17]. The objective function previously used for parameter identification is that given in equation (24). Optimization based on an objective function is used to validate the simulation results here as well. The objective function is more general than the more commonly used simpler form employed previously in [17]. This objective function is more useful to compare signals with different behaviors such as the appearance of the beating phenomenon in one response signal or of signals having different offsets (difference between the min/max values).(24)Obj=maxTexp.−minTexp.−maxTsuml.−minTsuml.2,where Obj is the objective function, Texp. the experimental friction force, and Tsuml. the simulated friction force.Table 3
Friction models parameters (6.5 N) [15].
Parameter
Coulomb friction model
Dahl friction model
Stribeck friction model
LuGre friction model
μ
c
0.520
0.500
0.161
0.160
μ
s
—
—
0.477
0.492
v
s (m/s)
—
—
6.101 10−4
7.152 10−4
u
crit (m)
—
—
4.146 10−6
4.060 10−6
σ
0 (N/m)
—
1.003 10+7
—
1.140 10+4
σ
1 (N·s/m)
—
—
—
347.030
σ
2 (N·s/m)
—
—
—
0.000
δ
2All the results are obtained at the excitation force amplitude of 6.5 N. A proportional (Rayleigh) damping is used to model the beam structural damping. The mass proportional damping coefficientαn is given by the following equation:(25)αn=2ωnζn,where ωn is the angular natural frequency of mode n and ζn the corresponding damping ratio.The value of the damping ratioζ in Table 1 is identified experimentally using an impact test and the logarithmic decrement method.


## 5. Parameter Validation

All the results of the simulation are obtained using the software Abaqus version 6.12-3. Simulations are carried out for 1.2 seconds, with a 10−4 second fixed time step, using the implicit integration method. Data storage was done in 8 × 10−4 second time intervals in order to optimize the total simulation time. Two models were tested, the Coulomb friction model with two different formulations and Lagrange and Penalty formulations, and the decay friction model. We begin by presenting the modal analysis results.Table4 presents the first resonance frequency. The simulations using the two Coulomb (Lagrange and penalty) models give approximately equal first resonance frequencies (with a difference of 0.2 Hz); the predicted frequencies are approximately 3 Hz higher than the experimental frequency of 50 Hz. However, the decay friction model is closest to the experimental results when compared to the Coulomb models. For the resonance frequency, the Lagrange formulation is better than penalty formulation, while the decay friction model frequencies are the closest to the experimental values.Table 4
1st resonance frequency.
Test
1st resonance frequency (Hz)
Experimental
50
Coulomb–Lagrange formulation
52.8
Coulomb-penalty formulation
53
Stribeck (δ=1) friction model
51.6In the following figures, from Figures7–13, the simulations of two cases for each Coulomb friction model formulation are presented. For the Lagrange formulation, simulations using the two friction coefficients identified via Coulomb and LuGre friction models (Table 3) are presented. However, for the penalty formulation, we present the simulation results obtained using the friction coefficient identified via the Coulomb model with two values of the slip tolerance Ff: the default value (0.005) and the optimal value (0.0015) identified during this analysis using the beam response envelope validation criterion.Figure 7
1ST Nonlinear normal mode: (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure 8
Beam response envelope (maximum deflection): (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure 9
FRF at driving point, Coulomb friction model: (a) Lagrange formulation, (b) penalty formulation, and (c) decay friction model.
(a)
(b)
(c)Figure 10
Slip displacement: (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure 11
Slip velocity: (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure 12
Friction force: (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure 13
Contact normal force: (a) Lagrange formulation; (b) penalty formulation.
(a)
(b)Figure7 presents the first nonlinear normal mode (NNM). Simulation results are compared with the experimental normal modes identified using the inverse harmonic balance method (IHB) and the linear normal mode (LNM) of the clamped-simply supported beam. The simulation results using both friction models have approximately the same NNMs; the NNM form is also very close to the first NNM identified experimentally using the IHB. However, the NNM of the simulation results using the penalty formulation are even closer to the experiments NNM compared to the Lagrange formulation results.Note that while the simulated 1st nonlinear modes (NNMs) appear similar to the experimental NNM, the small difference is important. This is because the modal spatial derivatives are particularly critical for accurate determination of the friction dynamics at the contact point; this is evident in equation (10). Hence, even though the NNMs look similar, the slight differences lead to significantly different responses (Figure 8). This is largely due to the importance of the first and the second derivatives of the mode form at the contact point which directly affects the resulting computed friction force Tt in equations (4) and (10).To better compare the results, we can use data shown in Figure8. In this figure, the beam response envelopes (maximal deflections) using both formulations of the Coulomb friction model are presented. From Figure 8(a), Coulomb friction model with the Lagrange formulation, we find that the friction coefficient identified using the LuGre model μc=0.49 yields the closer results to experiments than the value identified using the Coulomb friction model. Contrarywise for the penalty formulation (Figure 8(b)), the value identified using the Coulomb friction model (μc=0.52) gave the best results but with the optimal value of the slip toleranceFf=0.0015. The default value of the slip tolerance is unable to produce the correct results.To provide more details on the beam response, the frequency response (FRF) of the system at the driving point (position 555 mm from the clamped end) is presented in Figure9 for the Coulomb and decay friction models. Similarly to the observations made from Figure 8, the case with μc=0.49 for the Coulomb friction model with Lagrange formulation (Figure 9(a)) and the case with μc=0.52,Ff=0.0015 for the penalty formulation (Figure 9(b)) are the closest results to the experimental response. These two cases can be considered as the optimal cases (having optimal parameters) for each model formulation.From Figure9, we can also observe that the penalty formulation gives smooth system response (without multiple-resonance peaks phenomenon), in the interval 45–56 Hz, compared to the Lagrange formulation, where this phenomenon is not observed experimentally in the interval (from 45 Hz up to the second natural frequency (160 Hz)). However, for both formulations, there are other resonances above 56 Hz, while not presented in the results, we can see the beginning of a second possible resonance after the 1st natural frequency in Figure 9(b). Furthermore, the experimental response occurs over a wider frequency bandwidth. The large response predicted by the Coulomb friction model with both formulations shows that the models do not correctly capture the energy dissipation due to friction. This is partly expected due to the known inadequacies of the Coulomb friction model.Note, however, that there is gradual improvement in the predictive behavior between the more basic Coulomb and the improved decay friction model. The preliminary result of decay friction model peak frequency, in general, is closer to the experimental resonance frequency value than Coulomb model (Table3); however, the beam response amplitude (Figure 9(c)) is higher than the experimental result. The results of Figure 9(c) are obtained using the optimal parameters of the penalty formulation μs=0.52,Ff=0.0015 and the optimal value of μc=0.16 presented in Table 3 from the previous work [15] with the Abaqus decay friction model. In the next step, a complete analysis of the decay and Stribeck models must be done to improve our decay friction model.Figures10 and 11 present the slip displacement and slip velocity at the contact point. For both formulations, the simulated resonance slip displacement (Figure 10) of the optimal parameters defined above: μc=0.49 for Lagrange formulation and μc=0.52,Ff=0.0015 for penalty formulation, leads to the closest values to experimental results. However, all the simulated resonance slip velocities at the contact point (Figure 11), for both formulations, are higher than the experimental values. This may come from the integration and the derivative methods used in Abaqus and/or due the inability of the Coulomb friction model to represent the real behavior of friction especially in the sticking condition. Thus, for both formulations at low frequencies, below the resonance, the system is in the sticking condition (zero slip displacement for the Lagrange formulation and slip displacement equal to the critical displacement for the penalty formulation). For the higher frequencies above the 1st resonance frequency, the slip displacement (Figure 10) values are closer to the experimental values. The slip velocity (as shown in Figure 11) for the Lagrange formulation is higher than the penalty formulation case. Both formulations thus give higher results with the penalty formulation being 1.5 times the experimental values.We can also make the same observations in Figure8, where for the beam frequency response at the driving point, for the low frequencies, the difference between the simulations and the experimental results is quite large. However, for the higher frequencies (above the resonance frequency), there is a good match between the simulations and the experimental results. We conclude that the Coulomb friction model is valid and applicable only for the cases with slipping conditions. Despite the penalty formulation improvement, the Coulomb friction model accuracy remains insufficient to represent the real behavior of the sticking or presliding condition.Figures12 and 13 present the contact force, tangential (friction) force, and normal (impact) force results. The friction force results, Figure 12, support the conclusion above for both formulations. The simulated friction force levels at the resonance frequency have closest values to the experimental results when the optimal parameters defined above are used. However, the simulated normal force levels obtained using the optimal parameters (Figure 13) at the resonance frequency are slightly higher than the experimental values, with an error of 4.39 N (4.93%) for the Lagrange formulation and 3.56 N (4%) for the penalty formulation. Thus, even when the results of the friction force, the system response, and the slip displacement are close to the experimental values, both formulations of the Coulomb model lead to slightly different values for the normal force. The continuous nature of the penalty formulation slightly improves the normal force result compared to the discontinuous Lagrange formulation.Finally, we can use these results (Figures8, 9, and 11) to obtain the optimal Coulomb friction model parameters for the FEM model. The optimal parameters are (Ff between 0.0015 and 0.005 with μc=0.52) for the penalty formulation. For the Lagrange formulation, μc is less than 0.49 for the response system validation criterion and between 0.52 and 0.49 for the others validation criteria. In conclusion, we note that the accurate optimal parameters of the different Coulomb friction model formulations affect the system response, the friction force, and the slip displacement in the same manner and differently to the normal force.


## 6. Parameter Sensitivity Analysis

In this section, we investigate the parameter sensitivity of the two Coulomb friction model formulations. In Figures14–19, we present the simulation results using the parameters of the two cases presented in Figure 7 to Figure 13 plus other cases near the optimal parameters of both formulations. However, we present only the figures (or the validation criteria) where the parameter variation effect is clear. Furthermore, the sensitivity of the two parameters in the penalty formulation is presented separately: subpart “b” in Figures 14–19 for the slip tolerance Ff analysis and subpart “c” in Figures 14–19 for the friction coefficient μc analysis.Figure 14
Beam response envelope (deflation maximal): (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)Figure 15
FRF at driving point: (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)Figure 16
Slip displacement: (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)Figure 17
Slip velocity: (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)Figure 18
Friction force: (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)Figure 19
Contact normal force: (a) Lagrange formulation and penalty formulation; (b)Ff sensitivity; (c) μc sensitivity.
(a)
(b)
(c)In the Lagrange formulation (Figure14(a)), when the optimal parameter (μc=0.49) is increased by 2%, the beam response envelope decreases by 17% relative to the optimal parameter response. However, the beam response envelope increases by 28% of the optimal parameters response when μc is decreased by 2% below the optimal value μc=0.49.For the penalty formulation (Figures14(b) and 14(c)), firstly, the slip tolerance Ff sensitivity was tested by fixing the friction coefficient μc at 0.52 (Figure 14(b)). It was found that the beam response envelope increased with increasing slip tolerance Ff (which means with increasing critical displacement ucrit). Additionally, it was noted that the optimal slip tolerance Ff could be improved within the interval (0.0011 to 0.0015) by using the experimentally identified μc=0.52. The friction coefficient μc sensitivity was also tested (Figure 14(c)) by fixing the slip tolerance Ff to two values: the default value Ff=0.0050 and the optimal value Ff=0.0015. These two cases are analyzed to understand the behavior of the μc variation. It was found that the beam response envelope increased with decreasing μc for both cases.The same approach was used to analyze the other criteria presented in Figures13–17.For the following two criteria, FRF at driving point (Figure15) and slip displacement (Figure 16), the observations of the parameters sensitivity are similar to the beam response envelope (Figure 14) where, the simulation result values increase with decreasing μc for both criteria using both formulations (Figures 15(a), 15(c), 16(a), and 16(c)), and the simulated values increase with increasing Ff for the penalty formulation. Furthermore, for the penalty formulation, using the slip tolerance optimal value Ff=0.0015, the simulation result approaches the experimental value, unlike the case using the default value Ff=0.0050, with the appearance of a multi resonance, confirming the stability response of the slip tolerance Ff optimal value.For the slip velocity criterion (Figure17), the observations of the parameters sensitivity are similar as the FRF at driving point (Figure 15) and slip displacement (Figure 16) for the penalty formulation (Figures 17(b) and 17(c)). However, for the Lagrange formulation, the simulated slip tolerance (Figure 17(a)) behaves differently around the optimal friction coefficient, μc=0.49. The simulated value increases away from the experiment at value for increasing or decreasing friction coefficient μc. This is another confirmation that μc=0.49 is the optimal value for the Lagrange formulation.For the following two criteria, the friction force (Figure18) and the normal force (Figure 19), the simulation values decrease for decreasing μc in the Lagrange formulation (Figures 18(a) and 19(a)) and also by decreasing Ff using the penalty formulation (Figures 18(b) and 19(b)). Furthermore, for the penalty formulation, using the default slip tolerance Ff leads to values far from the experiments, and there is a small difference for the values around the optimal value Ff=0.0015. For the penalty formulation, the friction force (Figure 18(c)) increases with increasing friction coefficient μc. However, the normal force (Figure 19(c)) decreases with increasing μc using the optimal value μc=0.52, which leads to the closest values in experiments for both criteria.Finally, we conclude from this analysis that using the optimal parameter values (μc=0.49 for the Lagrange formulation and Ff=0.0015 with μc=0.52 for the penalty formulation) yields the smoothest (the most stable) simulated results and the closest to experimental values for the six validation criteria presented above. This validates and confirms the accuracy of the friction model parameters identified using the inverse harmonic balance method IHB with LuGre model for the Lagrange formulation and with Coulomb friction model for the penalty formulation. We can conclude also that the accurate optimal parameters of the different Coulomb friction model formulations affect the system response, the slip displacement, and the slip velocity in the same manner. The friction force and normal forces are affected differently.We can use the results of Figures7–18 to obtain the best (optimal) Coulomb friction model parameters for the FEM model. The optimal parameters are for Ff between 0.0015 and 0.0020 with μc=0.52 for the penalty formulation and μc between 0.48 and 0.49 for the Lagrange formulation.The friction coefficient identified experimentally using the Coulomb model is the optimal value for the penalty formulation of the Coulomb model. However, the friction coefficient identified experimentally, using the LuGre friction model, is the optimal value for the Lagrange formulation of the Coulomb model. This can be generalized only after testing different cases in future work.Modeling the Stribeck effect also significantly improves the resonance frequency prediction. The next step involves moving to more sophisticated models—the Dahl, LuGre, and, eventually, the Hybrid friction models. Prior to this, the role played by the numerical implementation of these models in the FEM codes needs to be verified. This is important because the implementation of the friction models in FEM codes is not mathematically trivial and involves adjustments and additional parameters.Finally, a word concerning the slip displacement: this displacement is an important quantity used to calculate the wear work rate (26). The slip displacement value at the contact point is plotted versus frequency in Figures 10 and 16 for the Coulomb model. The values of the slip displacement generated by the Coulomb friction model are close to the experimental ones at resonance and too small (even equal to zero for the Lagrange formulation) for frequencies lower than the resonance frequency. This leads to smaller computed work rates compared to the experimental values in this frequency interval.(26)W˙=∫0tNt⋅u˙tdτ∫0tdτ,where W˙ is the work rate, Nt is the normal force, u˙t is the sliding velocity, and t is the duration of contact.


## 7. Conclusion

This paper presents the first step to validate, using 2D plane strain FEM elements, the parameters identified in previous work [15] using the inverse harmonic balance method results and the basic models of the hybrid friction model.The accuracy of the friction model parameters identified using the inverse harmonic balance method (IHB) with the LuGre model is validated using the Lagrange formulation. The parameters identified with Coulomb friction model are validated using the penalty formulation. The accuracy of the optimal parameters of the different Coulomb friction model formulations affects the system response, the slip displacement, and the slip velocity in the same manner but differently from the friction and normal forces.Both models produce friction forces, normal forces, and slip displacements of the correct order of magnitude compared to the friction force calculated using the inverse harmonic balance method at the resonance frequency. However, their FRF bandwidths are significantly far from the experimental results. Furthermore, Coulomb and decay friction models yield modal parameters (resonance frequency and NNMs) that are close to those of the experiments. The decay friction model yields the results closest to experiments for the NNMs and resonance frequency.The Coulomb friction model is valid and applicable only for cases with slipping conditions. Despite the improvement gained using the penalty formulation in the Coulomb friction model, it remains insufficiently accurate to represent the real physical behavior for sticking or presliding condition.The present work demonstrates that both static friction models produce results closest to the experiments at resonance in the slipping regime, but are incapable of accurately representing all the behaviors of friction, especially in the sticking regime, ultimately leading to incorrect work rate estimates in the complete frequency range. The results presented here represent a part of ongoing work aimed at modeling the detailed physics underlying friction.



---

*Source: 3493052-2019-01-22.xml*