Frictional Contact Finite Element Analysis of Pneumatic Cylinders

Friction is present in most mechanical systems. The aim of this paper is to analyse friction between a pneumatic cylinder and its piston seal using the finite element method (FEM). The high-fidelity numerical model used real geometry and material characteristics of a double acting pneumatic cylinder. The simulations required fine tuning: the piston seal was created with a hyperelastic, nearly incompressible material model utilizing mixed u-p elements. Three friction characteristics were analysed with slow velocity motion of the piston: Coulomb, an extended Coulomb and the Benson-Hallquist model. The result of the Coulomb model has given a constant friction force after an oscillation at first. The extended Coulomb friction model is used to show the effect of cohesion between two surfaces. After that, the Benson-Hallquist model was simulated by FE and then experimental and FE results were compared. As a result, an exponentially decreasing friction force was obtained. It was shown that the friction behaviors can be simulated with FE and this method will allow studying the effects of different types of piston seal and cylinder barrel on the friction force.


Introduction
For industrial automation speed and reliability with simple construction and effortless maintenance are the most wanted properties.These demands can be ensured by pneumatic actuators.These devices transform pneumatic energy into linear motion.This paper examines a double acting pneumatic cylinder (Fig. 1).The compressed air enters into the tube and moves the piston in both extend and retract strokes.The pressure p 2 on the piston imparted by the air creates the driving force F. The double acting cylinders are used when the driving force is important in both directions.The seals in cylinders are o-rings or membranes.Studying the factors affecting the piston motion are important for the control.Friction has a high priority, because this nonlinear factor causes the difficulty of the exact location / velocity control.The first steps on pneumatic system control and modelling are dated to the 1950's, to Shearer's work [1].Since then, there are numerous dynamic models for pneumatic actuators.For example, Araki et al., [2] produced linear models of servopneumatic systems.Choi et al., [3] modelled the nonlinear friction property of an electro pneumatic servo valve.The friction model had two parts: a static and a dynamic one.The static friction model represented the stick-slip motion while the dynamic model incorporated other friction properties such as the "presliding" displacement.Persson [4] sets up a mathematical model for the friction hysteresis.
In this paper the possibility and the limits of the FEM are investigated with respect to modelling the friction with a servopneumatic cylinder's piston, tube and seal geometry.An experimental setup on a Festo pneumatic cylinder provided the basis to observe friction at low velocities (see Széll et al., [5]).Different friction models were examined in FEM to generate the modality of the experimental results and to prove the adaptability of using the model in other simulations.For an accurate simulation the geometry data of the pneumatic piston was needed.In this paper a Hoerbiger-Origa rodless double acting pneumatic cylinder's geometry was used.
Besides the piston, the viscoelastic piston seal is another research area: Váradi and Pálfi [6] investigated the viscoelastic material properties of an EDPM 75 rubber by experiment and FE simulations.Annavarapu et al. [7] examined different finite element methods for frictional sliding, such as the Nitsche and the Lagrange method.Available FE friction models were compared with each other.Gaul et al. [8] developed a special simulation framework which allows efficiently modelling interfaces with friction and normal contact by appropriate constitutive laws which are implemented by contact elements in a finite element model.
The structure of the paper is as follows.In the second section the pneumatic cylinder will be introduced.The third section presents the theoretical background of the examined friction models and the experimental results.In the fourth section the chosen material model, the FE geometry, the mesh, the simulation's stability problems and their solutions will be discussed.The fifth section reports the results.The sixth section concludes the work and describes the direction for further simulations.

The modelled cylinder 2.1 Geometry and material properties
The Festo pneumatic cylinder's catalogue only contains the fastening, stroke length, the possible loads and the operating properties.Thus a properly similar Hoerbiger-Origa cylinder was disassembled and its geometrical properties were recorded.The cylinder is shown in Fig. 2 (the piston's velocity is v), the stroke length is L=700 mm, and d=32 mm is the diameter of the cylinder.The cylinder's thickness is h=3 mm (shown in Fig. 3).Since the influence of friction on the motion was examined, we had to ensure the right amount of normal force with a certain overlap between the cylinder and the piston seal during the assembly.In the FE model for the 6 bar maximum pressure s*=0.11mm extrusion gap had to be generated.The geometry and material properties are summarized in Table 1.The double acting pneumatic cylinders have two pistons with one piston seal on each side (installation is shown in Fig. 3).According to the piston seals' installation catalogue the material is NBR (acrylonitrile butadiene rubber).The density ρ of the seal is 1.31 g/cm 3 (see Table 2).The friction coefficient µ was chosen to be 0.05.For a hyperelastic material such as the rubber, the Poisson-ratio is ν@0.5.For the NBR seal the Young's modulus E is 5 MPa, while for the bulk modulus K=2980 MPa was assumed.3 Friction This synopsis is only to briefly review some of the most widely applied friction models.Detailed overview can be found for example in Dupont et al. [9].To understand the physics of friction in pneumatic cylinders, see Fig. 4 (and details in Subsection 3.1).On the free body diagram of the seal only the result forces are shown.Here F n is the normal force, F f is the friction force, Δp is the pressure difference between the sides of the piston and A is the piston's area.
Gu, Ji-Cheng, et al. [10] analysed quasi-static slip motion and instability for single degree of freedom elastic systems.Benson and Hallquist [11] presented a contact algorithm in FEM and showed how to calculate the contact force included the friction.Fiedler et al. [12] used a user written subroutine for contact recognition.Zhao et al. [13] gave finite element formations and iteration processes for static and dynamic contact problems.Kennedy et al. [14] analysed the surface temperature of sliding systems with FEM.

Friction models
The most known friction model is the Coulomb friction between two solid surfaces in contact.Between non-moving surfaces (v=0) the Coulomb friction force (F c ) is where F n is the normal force exerted by each surface on the other and µ s is the static coefficient of friction.
If the force -which prevents motion between the surfacesis overcome (so the body moves v≠0) then the friction force is equal to the product of the dynamic coefficient of friction and the normal force.The dynamic friction coefficient is lower than the static coefficient of friction.
where µ d is the dynamic coefficient of friction.Hence the friction force as the function of velocity in the Coulomb friction model: where F f is the friction force between the two bodies and v is the relative velocity between the surfaces (see Fig. 5a).This is the simplest model of the friction, since the dry friction is not the function of the contact surfaces' size.
In the extended version of the Coulomb model, the cohesion (provides sliding resistance even with zero normal pressure) between the two bodies was taken into account.On the occasion that the cohesion is zero, we can re-create the Coulomb model (for details see Section 5).
Another friction model is an exponential friction model, according to Benson and Hallquist [11]:

Experimental data
The experiments were carried out on a standard Festo pneumatic cylinder.As excitation the left chamber was loaded with a slowly rising pressure while the right chamber was exhausted.The measurements were taken right in the area of the friction hysteresis.From the measured acceleration the velocity was calculated by integration.The difference of the force needed for the motion and the force acting on the piston was dissipated by the friction.Thus, the friction force as the function of velocity was measured, which is shown in Fig. 6: Figure 6 shows that the friction force is decreasing at low velocity.The parameters of the Benson-Hallquist model were determined by a nonlinear fit to be µ s =0.0551, µ d =0.0395 and c=0.6297.The fit is also shown in Fig. 7.
The normal force could not be defined from the experiments, thus the value 68 N was used based on the results of the Finite Element simulations using the Coulomb friction model (see details in Section 5).The FE analysis presented in this paper is capable of recreating the experimental curve at low velocity and positive acceleration.

The numerical simulation
Table 1 and Table 2 contain the geometry data of the FE model.

FE materials
For the materials of the piston and the cylinder see Tables 1-2.The piston seal was made of NBR, but the engineering database of the FE software does not contain this type of material.
To define it, a material model had to be determined.For the simulation of rubber-like substances in FEM the Neo-Hookean hyperelastic model can be used because of the advantageous mathematical manageability.The Neo-Hookean model is similar to Hooke's law, and can predict the nonlinear stress-strain behaviour of materials undergoing large deformations.
To determine this type of material model's stress-strain curve two data are needed: initial shear modulus (G 0 ) and incompressibility parameter (D 1 ).Using Table 2 data:

Meshing, elements and properties
When evaluating the results of the Coulomb friction model, grid independence tests were made with three different grid densities (Table 3).A structured mesh was used, setting the cell size 2, 0.8 and 0.4 mm.Comparing the friction force calculated from the coarse mesh showed 27 % difference to the one from the medium mesh.The friction force difference between the medium and fine mesh is negligibly small (3 %).Thus we conclude that the medium mesh is acceptable for our simulation.This mesh is shown in Fig. 8.

Fig. 8 The utter mesh
Contact friction is a material property used with contact elements.In 2D the CONTA175 element (see Fig. 9) can be used to represent the contact and sliding between two surfaces.This element located on the surfaces of solid, beam and shell elements.Contact occurs when the element surface penetrates one of the target segment elements so at least one contact element covers another.Contact friction is allowed with this type of element.On the target side the TARGE169 element was chosen.The contact and target elements shared common real set parameters.Since the pneumatic cylinder has two piston seals, thus the surface of the cylinder takes part in defining two contacts at the same time.(6) (7) Fig. 9 The CONTA175 element In the 2D model to treat the penetration after making the planar geometry the cylinder was offset.The cylinder is not connected to the piston seal but the contact pairs can be handled by the program.Contact friction can be specified either through the coefficients of friction models or as user defined friction properties.

Initialization
The simulation was done in two phases (load steps).During the first load step the cylinder was pressed onto the piston seal.This means a displacement load for the cylinder, at 1 s time with 0.01 s time step the cylinder moved 0.5 mm towards the -x direction (Fig. 10).At the second load step the velocity load was added and the selected friction models were solved.

Convergence problems
There were 3 major issues for obtaining convergence: nonlinear material property (instability), contact stiffness was too high and the choice of the "symmetric" contact.
We can determine the unbalanced force locations by plotting Newton-Raphson residuals (Fig. 11).

Fig. 11 The Newton-Raphson residuals
We can optimize the structural analysis convergence by using different solutions: The convergence problem of the Neo-Hookean material model -used for the piston seal material (NBR) -can be solved with mixed u-p elements.These elements use both displacement and pressure as primary unknown variables.
The contact stiffness can be reduced to correct the force unbalance.The stiffness allows penetration between the associated elements.The contact problem in this research was solved with a contact stiffness factor of 0.1.
When defining the contacts not just the denotation of the curves but choosing the contact type is also significant.There are two possible options: asymmetric contact is defined as having all contact elements on one surface and all target elements on the other surface.Symmetric means when each surface can be both a target and a contact surface.The asymmetric contact option was chosen: the cylinder was the target surface and the piston seal was the contact surface.
Even so, transient problems have appeared: non-physical solutions can be detected for a short period of time (see 0.03-0.05s of the solution in Fig. 13).

Simulations and results
During the simulations three increasingly complex model were examined.The loads were also defined accordingly: while in the first model the set-up had a linear velocity load, in the next case, the Benson-Hallquist model had the velocity curve load (with positive acceleration) taken from the experimental data, see Fig. 12.

Coulomb model -Model A
Coulomb friction is the simplest model, using this in simulations also serves to check whether a known solution is reproduced.In this case the friction coefficient was constant 0.05 and a=0.5 m/s 2 corresponding to the linear velocity profile v=a•t.The final result of the simulation using the Coulomb model can be seen in Fig. 13.The crescent velocity of the piston affected the friction force (Fig. 14).After the compression from 0 to 0.03 s when v=0 m/s the awakening residual force was the static friction.When the piston moved at very low velocity there was an oscillation in the friction force, which was smoothed down within 0.02 s.The cause of the oscillation is the determination of the contact regions by the software.After 0.05 s a friction force with permanent magnitude (3.68 N) is visible.This is expected; since the normal force was about F n ≈68 N. Based on these results, in all other simulations the transients (0.02 s) were discarded.The result of the friction force with Model B can be seen in Fig. 16.For a slight change in the value of cohesion, the results are fairly insensitive.When b was 20000 Pa, the friction force curve had the same nature but its value showed a significant difference.Simulations were run with different τ lim values but the results were completely the same as the cohesion runs.
Changing the cohesion value to b=0, the results are nearly the same as the Coulomb friction model (see Fig. 17

Benson-Hallquist model -Model C
During FE simulations there is an option to define the Benson-Hallquist model.The friction force can be calculated as the follows: , where F f is the force of friction, the superscripts n and n+1 denote the n and n+1 time steps, and the superscript t denotes a trial value of F f .The maximum possible friction force is F y .
In the exponential friction model (according to Benson-Hallquist see Subsection 3.1 and 3.2) c=0.6297,R f =1.3949.As we see the friction force curve from the finite element simulation is slightly different from the curve character obtained during the measurement (see Fig. 18).At low velocity (0<v<0.005m/s) the finite element simulation result is a much steeper decline compared to the measured data set.At higher velocities, the friction force in FEM is smaller than expected.The oscillation detailed in Subsection 5.1 is present in this case too, but the curve shows an exponential decay character after all, so the FE program applied the friction model correctly.Solving the problem of oscillation requires further investigations but with this model the described experimental data results -showing a decreased friction force -can be simulated.

Summary and Conclusion
In this paper the piston motion of a Hoerbiger-Origa double acting rodless pneumatic cylinder was studied with FEM.The simulations were carried out with a 2D model after mesh-independence examinations.Among the factors affecting the piston movement the friction has been highlighted: different friction models were introduced.Solutions are shown for convergence problems due to the contact model.During the simulations three velocity loads were used.
In Model A (Coulomb friction) the friction force showed a velocity-dependence at first but after 0.03 s it reached a constant 3.68 N value regardless of the velocity.
In the extended Coulomb model the friction force was velocity dependent.When the cohesion value was set to zero, the results of the Coulomb friction were obtained.
For the Benson-Hallquist model the missing coefficients were obtained from the measurement data with curve fitting.The results show that the friction force varies exponentially with velocity but it is steeper than the experimental data.The model is capable of giving back the character of measurement data.
The set of the user defined friction models had outlined.Further simulations may be needed such as FSI simulations.Based on these simulations, the piston movement can be examined and the effects of thermal processes in motion can be observed.

Acknowledgement
The author would like to express her thanks to András Czmerk, Zsolt Szabó and Tamás Kalmár-Nagy from Budapest University of Technology and Economics for the experimental results and all the support.The project presented in this article is also supported by János Szegletes from eCon Engineering Kft.The work relates to the scientific programs "Development of quality-oriented and harmonized R+D+I strategy and the functional model at BME" (Project ID: TÁMOP-4.2.1/B-09/1/KMR-2010-0002) and "Talent care and cultivation in the scientific workshops of BME" (Project ID: TÁMOP-4.2.2/B-10/1-2010-0009). (9)

Fig. 4
Fig. 4 The schematic figure of the cylinder and the seal's free body diagram d is the dynamic coefficient of friction, μ s =R f •μ d is the static friction coefficient, R f is the ratio of static and dynamic friction and c is the decay coefficient (see Fig.5c).Further, complex friction models are the viscous, Stribeck, Trustin and the LuGre model.Couillard et al.[15] proposed an identification procedure to evaluate the ability of the LuGre friction model to predict small amplitude frictionally damped vibrations.The damping effect of friction was introduced in[16].Mingfu et al. designed a dry friction damper which significantly reduces vibration in a rotor system.Calculating the friction force in FEM is detailed in Section 5.

Fig. 5
Fig. 5 Different friction curves (Friction force as the function of velocity) a) Coulomb friction, b) viscous friction c) Benson-Hallquist model

Fig. 12 Fig. 13
Fig. 12 Velocity load as the function of time

Fig. 14
Fig. 14 Friction force as the function of velocity, in Model A (Coulomb friction)

Fig. 15
Fig. 15 Contact friction stress as the function of contact pressure

Fig. 16
Fig.16 The friction force as the function of velocity -Model B (cohesion model) ).

( 8 )Fig. 17
Fig. 17 The friction force as the function of velocity in Model B with 0 Pa

Fig. 18
Fig. 18 Friction force as the function of velocity in Model C (Benson-Hallquist)

Table 2
Piston seal material properties

Table 3
Details of the used meshes with different mesh density