Validation of normal and frictional contact models of spherical bodies by FEM analysis

Contact forces between two spheres are computed, including the contact pressure (normal) and the frictional stress (tangential) using a finite element method (FEM). A CAD model of a part of a sphere was developed. A mesh was created using ANSYS® Solid 186, 20-Noded hexahedral element and analyzed for its sensitivity. ANSYS® Contact 174 and Target 170, 8-Noded surface elements were used. Contact pressure and frictional stress contours were calculated by varying the displacements. Normal and Tangential contact forces were computed by integrating contact pressure and frictional stress over the contact surface. The values obtained for the normal force were compared with the non-linear spring model as given by Hertz [1]. Similarly values of the tangential force were compared with the model of Mindlin and Deresiewicz (MD) [2]. The FEM results were found to be in agreement with the models. INTRODUCTION There are two components of the contact force: (1) the normal contact force and (2) the tangential contact force. The normal contact force arises from the normal overlap of particles as given in Figure 1. Hertz [1] model has certain limitations as highlighted by Johnson [3]. The assumptions made in the Hertz [1] model are as follows: • The surfaces are continuous and non conforming • The strains are small • Each solid can be considered as an elastic half-space • The surfaces are frictionless The normal contact force between spherical particles is given by Hertz [1] based on theory of classical mechanics. For two spheres in contact with radii, R1 and R2, and relative normal displacement, 2a, the radius of contact area, c, is given by


INTRODUCTION
There are two components of the contact force: (1) the normal contact force and (2) the tangential contact force.The normal contact force arises from the normal overlap of particles as given in Figure 1.Hertz [1] model has certain limitations as highlighted by Johnson [3].The assumptions made in the Hertz [1] model are as follows:

•
The surfaces are continuous and non conforming • The strains are small • Each solid can be considered as an elastic half-space • The surfaces are frictionless The normal contact force between spherical particles is given by Hertz [1] based on theory of classical mechanics.For two spheres in contact with radii, R1 and R2, and relative normal displacement, 2a, the radius of contact area, c, is given by (1) And the normal force, P, by

Int. Jnl. of Multiphysics Volume 4 • Number 2 • 2010
Where the Equivalent Young's modulus 'E * ' is (3) And R * is (4) The normal traction p, i.e. the distribution of pressure on the contact area, is given by (5) Where r is the distance from the centre of the contact area, and p o is the maximum normal pressure, which can be calculated using Eqn.6.

(6)
MD [2] studied the frictional forces using classical mechanics during contact and developed relations which are not only dependent on the current loading state but also on the history of Figure 1 Normal contact of spheres (reproduced from Zhang and Quoc [4]).
the loading.The main assumptions highlighted by MD [2] are given as follows: • The normal force and contact area are computed using the model given by Hertz [1] • Whenever there is change in normal traction it is accompanied by change in tangential traction, and if the change is more than the product of the coefficient of friction and normal traction, slip will occur

•
The annulus of slip progresses concentrically inwards • When slip occurs, then the product of normal traction and the coefficient of friction will be equal to tangential traction.• Contact parameters are computable if every previous step of loading is known from the equilibrium state Figure 2 shows the tangential contact between spherical particles.Denoting the tangential displacement of particle 1 and 2 as d 1 and d 2 respectively, tangential force felt by each sphere due to contact is Q.An expression for Q is given in Eqn. 7.
There are number of cases that have been considered by MD [2] to calculate the tangential forces during contact.All these cases vary based on the loading history.Zhang and Quoc [4] have simplified the MD [2] model based on a simple loading history and considered the following cases.

•
Normal force is constant • Tangential force is increasing • Tangential force is decreasing In all cases given above, it is possible to calculate the change in the tangential force using a spring constant, which differs between cases, as Where Q′ is the tangential force after and Q is the tangential force before tangential displacment ∆d = |d 1 − d 2 |, K T is the tangential spring constant (varies based on loading cases).
Under different scenarios of loading and unloading as discussed above; the tangential force non-linear spring constants can be computed.Figure 3 shows the relationships of non-linear Figure 3 Variation of tangential force with tangential displacement; normal force constant, tangential force increasing and decreasing (reproduced from Zhang and Quoc [4]).
spring constants when normal force is constant.O is the point of zero tangential loading.A is the maximum tangential loading in unit direction.B shows the point of zero tangential force but non-zero tangential displacement.C is the point of zero tangential displacement and nonzero tangential force.D shows the point of negative tangential force (opposite direction in 1-D) equivalent to point A. E is tangential force above A but in negative direction.F shows the point of positive recovery from maximum negative loading.The values of spring constants between each location of loading are shown in Figure 3. K T,0 is constant based on material properties as given in Eqn. 8. (8) Where v is the Poisson ratio and G is the Shear Modulus of the material and m is the coefficient of friction.
Similarly, Walten and Braun [5] presented a simplified model for normal and tangential contact of spheres.The advantage of using this model is the computational efficiency where there is no need to store the contact history which is required otherwise, hence accuracy is compromised.
Thornton [6] presented a simplified relation of measuring tangential force based on MD [2].This model is true for loading, unloading and reloading conditions.Moreover, Thornton [6] model also suggests how to compute the contact forces in 2-D contact.They have presented a displacement vector model to keep track of tangential direction, which may change as the particle moves.
In order to study the tangential force displacements model in detail, FEM analysis was conducted.For this purpose FEM model was built and solved using ANSYS® FEM commercial package [7,8], where the Augmented Lagrangian method [3,10,11] is used to compute the contact pressure and Coulomb law of friction [3,10,11] is used to compute frictional stress.Within the material, a linear static analysis has been performed.
Vu-Quoc, Zhang and Lesburg [9] have also presented the FEM analysis of spheres in contact using the ABAQUS® FEM commercial package for FEM analysis.They built a 2-D model for analysis and applied forces as initial conditions.In this study, a 3-D model of a sphere has been built and analyzed with the displacements as initial conditions.

FINITE ELEMENT MODELING
Figure 4 and figure 5 shows the CAD and FEM model built using ANSYS® Multiphysics.Solid 186, 20-noded hexahedral elements were used for the analysis based on the recommendation for non-linear analysis in ANSYS® Technical Manuals [9].Targe 174 8-noded quadrilateral target and Contac 170 8-noded quadrilateral contact elements were used for capturing contact variables.Problem specifications are shown in Table 1.
FEM mesh sensitivity tests were carried out to insure the quality of results.

FE LOADS AND BOUNDARY CONDITIONS
One part of spheres was constrained to zero displacements and other part was displaced to initiate the contact as shown in Figure 6.The values of normal and tangential displacements were controlled as per the loadings conditions as given in Figure 7.As shown in above Figure 5 loads (initial conditions) are applied.Each load is initiated from quasi-static state of the FEM model.

RESULTS AND DISCUSSION
Contact pressure and frictional stress were integrated over the areas to compute the net normal and tangential forces.Normal forces were calculated using the model given by Hertz [1] as described in Eqns. 1 to 6. Similarly, tangential contact force is computed using the MD [2] model using Eqns.7 and 8 and loading spring constant given by Zhang and Quoc [4].Comparison within normal and tangential force results is shown in Figure 8 and 9 respectively.As shown in Figure 8, FEM results for the normal force are in agreement with the expression of Hertz [1] as given Eqns. 1 to 6.
As shown in Figure 9, FEM results are in agreement with MD [2] tangential force contact model.
The normal pressure contour plot, tangential stress contour plot and contact status, at loading position 5 as given in Figure 7 are shown in Figure 10, 11 and 12 respectively.
The contact pressure distribution can be computed using equations 1 to 6.At 17 µm of relative normal displacement of spheres; according to Hertz [1] model maximum contact pressure p o is 581 MPa and FEM results gave the maximum stress of 585 MPa as shown in figure 10.Similar results were obtained for various loading conditions.This proves that FEM results are in good agreement with the Hertz [1] normal force model.
The MD [2] model assumes that the tangential stress distribution is axis-symmetric, in contrast to that found from the FEM analysis shown in Figure 11.In order to investigate it further contact status has been plotted as given in Figure 12.The contact status differentiates between stick and sliding regions.It is well-understood that frictional force is equal to coefficient of friction times the normal force in the sliding region whereas in the stick region its value has to be computed by tangential force models as discussed in section 2. It has been found that contact status is axisymmetric.In addition, it has been observed that frictional stress contours are axisymmetric in the sliding region (figure 11) and equal to the value of the pressure distribution (figure 10) multiplied with the coefficient of friction.In the non-stick region frictional stress contours are found to be non-axisymmetric.
According to MD [2] model maximum frictional stress is at the boundary of sticking and sliding region, which has also been found to be untrue as frictional stress distribution is non-axis-symmetric.

CONCLUSION
The Hertz [1] model for normal contact force has been validated using finite element methods using displacement as loading condition.Normal contact pressure distribution via finite element methods has been found match the parabolic relation given by Hertz [1].
MD [2] tangential force model and its derivatives as given in literature by Thornton [6] and Zhang and Quoc [4] have been validated using finite element methods.Despite the fact that the values for the tangential force are found to be in agreement with the theory, there is a mismatch in tangential stress distribution which requires further investigation.

Figure 5 FEMFigure 4
Figure 5 FEM model of 2 parts of spheres in contact; figure 6 shows the point of contact and boundary conditions applied.

Figure 9
Figure 9 Tangential force due to normal and tangential contact of spheres.

Figure 8
Figure 8 Normal force due to normal contact of spheres.

Figure 10 Figure 11
Figure 10 Contact pressure contour plot (MPa) for FEM load position 5.The area of contact has pressure > 0.
Authors acknowledge the support of Dr. Stuart Scott and Christopher Bohn from Department of Engineering, University of Cambridge for proof reading and helpful suggestions.

Figure 12
Figure 12 Contact status for the FEM load position 5.

Table 1
Specification of the Contact Problem