Advances in Tribology

Volume 2012 (2012), Article ID 832859, 11 pages

http://dx.doi.org/10.1155/2012/832859

## Improved Pressure Distribution in Elliptic Elastic Contacts between High-Order Surfaces

Department of Mechanics and Technologies, Stefan cel Mare University of Suceava, 13th University Street, 720229 Suceava, Romania

Received 30 December 2011; Revised 5 April 2012; Accepted 19 April 2012

Academic Editor: Shyam Bahadur

Copyright © 2012 Marilena Glovnea et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The improvement of mechanical contacts or microcontacts seeks a nearly uniform current density over most of contact area. When microtopography is homogeneous, this aim is achieved if nominal shape of contacting surfaces yields a nearly uniform central pressure which decreases monotonously to zero in contour points. These authors derived recently this shape for circular contacts by employing high-order surfaces. This paper extends this result to elliptical contacts. Some results are used to this end, derived for elliptical elastic contacts between high-order surfaces. As homogeneous high order surfaces lead to a highly nonuniform pressure distribution, central pressure is flattened by making the first derivatives of pressure vanish in contact center. Then, the contacts between fourth, sixth, and eighth, order surfaces are analyzed and recurrence relations for pressure distribution and contact parameters are proposed.

#### 1. Introduction

Both high-load-carrying capacity of the contact and the avoidance of pressure and stress risers require a contact pressure distribution as even as possible over most of contact area. Shape improvement of circular contacts uses a pressure distribution made of a central plateau surrounded by a narrow annulus of monotonous decrease to zero [1]. Such a pressure distribution yields if the equivalent rigid punch is bounded by polynomial surface of a higher order than two [1]. Hertz theory fails when dealing with such surfaces. Therefore, a new Hertz kind of theory is needed for high-order bounding surface.

For revolution surfaces, Shtaerman [2] proofed that a rigid punch of equation , pressed against an elastic half-space, generates a pressure expressed by the product between an even order polynomial in radius , of degree , and typical Hertz square root , being the outer contact radius. Klubin, and later on Popov, quoted in [3], found that a pressure given by the ratio of a Legendre polynomial of order , , to the Hertz square root generates a normal surface displacement which, within the circle, is proportional to .

For elliptical contact domains, Shtaerman [2] showed that a pressure written as the ratio between an even order polynomial of degree in and typical Hertz square root gives rise to a surface potential expressed by an even polynomial of degree in and when applied to an elastic half-space. This potential leads to an explicit integral expression of normal displacement in the points of bounding plane. In 1953, Galin [3] proved the following theorem: if a punch of front surface described by a polynomial of degree is pressed against an elastic half-space over a domain bounded by an ellipse of semiaxes and , the contact pressure can be written as the ratio between another polynomial, of the same degree , and Hertz square root. Gladwell [4] supplied an alternative proof of Galin’s theorem for a general anisotropic half-space. For transversely isotropic half-spaces, Gladwell found a polynomial for the displacement, expressed by definite integrals, when the pressure is the ratio between associated Legendre function and typical Hertz square root. The problem of normal indentation of an elastic half-space by a rigid frictionless axisymmetric punch described by a fractional power series of radial coordinate was analyzed by Borodich [5].

Later on, these authors [6] found that circular contacts bounded by even monomial surfaces of order lead to a nonuniform pressure distribution, having a local minimum in contact center, a circumferential maximum near contact boundary, and zero value in contour points. This pressure distribution was improved by imposing that the difference between current and central pressures possesses a multiple root in the origin of multiplicity order . Because all odd order derivatives vanish in the origin, this requirement is satisfied if all even order derivatives of pressure, up to order , are zero for . These conditions yield the required central pressure plateau.

This paper derives a similar procedure for elliptic elastic contacts between high-order surfaces aiming to get a flat central pressure. A generalized Hertz pressure of order is chosen to this end. This is defined as the product between typical Hertz square root and an even polynomial of degree in , , and of order in and . First, a recurrence relation is derived for the coefficients of polynomial which yield pressure distributions possessing a flat central plateau. Then the paper derives the surfaces of equivalent punches which generate the proposed pressure distribution when the punch is pressed against an elastic half-space. Finally, equations are derived for contact parameters, namely, contact ellipse half-axes, maximum central pressure, and normal approach or indentation of half-space.

#### 2. Basic Equations

The solution for elliptic elastic contacts is based on results establishing the class of high-order surfaces leading to an elliptical contact area and the correlation between pressure distribution and equations of bounding surfaces [6]. A generalized Hertz pressure of order has following form: The sum is an even polynomial of degree in , and of degree in and . This pressure is expressed in terms of elliptic parameter as follows: When applied over an elastic half-space boundary, generalized Hertz pressure generates following polynomial normal displacement: in which: being simple integrals depending on ellipse eccentricity: Equation (5) leads to combinations of complete elliptical integrals of first and second kind. Generalized Hertz pressure reduces to Hertz pressure when .

If the pressure is generated by a rigid punch pressed against the elastic half-space by a normal force , the deformed half-space boundary coincides with front surface of the rigid punch in the points of contact area. According to interference equation, front surface of the punch generating the pressure given by one of (1), (2) is expressed by following even polynomial of order : being normal approach or half-space penetration, maximum number of terms in an even polynomial of order , and the coefficients of punch surface.

Equation (6), in which is given by (4), defines an explicit one to one correspondence between generalized Hertz pressure and punch surface. A generalized Hertz pressure of order yields an even polynomial punch surface of degree with respect to co-ordinates and , with no free term.

#### 3. Contact of Homogeneous Surfaces

Applied to the contact of homogeneous even order surfaces, the above equations lead following pressure distributions expressed in terms of elliptic parameter : Axial pressure profiles predicted by (7) exhibit a non-uniform pressure central plateau, which is unsatisfactory for both electric current density distribution and contact strength. The degree of nonuniformity increases with surface order.

#### 4. Flat Central Pressure

All pressure distributions given by (7) possess an average negative concavity in central region. This requires important pressure maxima in peripheral zone to get a vanishing pressure in contour points. Intuitively, one feels that central pressure dimple decreases when adding lower order terms in surface equation, especially second order terms, which generate a Hertz-like pressure component compensating central crater. Of course, this means that contacting surfaces are not any longer homogeneous. Mathematically, the avoidance of peripheral pressure maxima requires central concavity of pressure distribution be zero or the difference between current generalized Hertz pressure and central pressure possess a multiple root in the origin of multiplicity order . This means that the first derivatives of pressure must be zero for . All odd order pressure derivatives contain the as factor and consequently they vanish in the origin. Therefore, outer pressure maxima vanish if all even order derivatives of pressure up to order take zero value when . These conditions yield the coefficients in (1) or (2) of pressure distribution which assure a flat central plateau. It is more convenient to use (2) because this implies fewer calculations. For instance, the second derivative of (2) with respect to elliptic parameter takes in the origin following expression: The derivative in (8) vanishes if the coefficient is: In a similar way, the forth derivative of pressure vanishes in the origin if coefficient has the following value: whereas make the sixth, eighth and tenth derivatives of pressure to vanish in the origin, respectively.

The above values of coefficients yield simply following general recurrence relation:

Equation (12) allows writing the equation of pressure distribution of any desired order, which possesses a flat central plateau surrounded by a peripheral zone of monotonously decreasing pressure. Several such pressure distributions are:

These pressure profiles exhibit a well-defined central pressure plateau surrounded by a monotonous decrease to zero in contour points. Moreover, as increases, the central flat region of pressure distribution extends, whereas maximum pressure decreases. This means that either maximum pressure decreases for the same load or load carrying capacity increases for the same maximum pressure. It is thus obvious that newly proposed pressure distributions are superior with respect to classical Hertz pressure.

#### 5. Punch Surface and Contact Parameters

Once pressure distribution is known, the coefficients of polynomial surface of equivalent punch result by coefficient identification in (6). Several surfaces leading to flat central pressure are derived below for , and by using the expressions of integrals given in [6].

##### 5.1. Fourth-Order Surfaces

In the case , which means fourth order surfaces, (6) takes following form: As , (19) becomes According to [6], the integrals and are: The integrals are expressed in terms of complete elliptical integrals of first and second kind in [6]. By using these expressions of , coefficient identification in (20) leads to following equations of normal approach and surface coefficients: Coefficients depend all on the eccentricity of contact ellipse. This eccentricity is found by aid of two of these coefficients, arbitrarily chosen. Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity by involving and . The division of (23) and (25), member by member, yields the following transcendental equation having the eccentricity as unknown: The eccentricity is the solution of (28) which is solved numerically. Once knowing the eccentricity, remaining coefficients , , and result from (24), (26), and (27). Therefore, surface polynomial possesses only two independent coefficients, in this case and , which yield ellipse eccentricity. The remaining coefficients must have specified values given by (24), (26), (27).

Force balance equation on direction yields following relation between central pressure and normal load :

By substituting (30) for central pressure in (23) and (25), one finds following expressions of ellipse half-axes:

Maximum contact pressure, which occurs now in contact center, yields by substituting (31) for ellipse half-axes in (30):

Finally, (22) leads to following expression of normal approach: In all above equations SI system is used.

As guessed initially, the surfaces of equivalent punch leading to a flat central pressure are non-homogeneous. They contain all even powers of coordinates, up to .

Punch surface is described by following equation: where the coefficients and are arbitrary imposed and , , given result from (24), (26), (27).

##### 5.2. Sixth-Order Surfaces

As in previous case, the coefficients of polynomial surface of equivalent punch leading to flat central pressure result by coefficient identification in (6). Now , which means sixth order surface, and (6) takes following form: As and , (35) becomes As established in [6], the integral is given by following expression:

The substitution of expressions (21), (37) for , , and , respectively, in (36) and identification of coefficients yields following results:

All coefficients depend on the eccentricity of contact ellipse. As in previous case, this eccentricity is found by aid of two of these coefficients, arbitrary chosen. Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity by involving and . Dividing member by member equations (39) and (42), following transcendental equation having the eccentricity as unknown results:

The solution of (48) is found numerically. Once knowing the eccentricity, all remaining coefficients , , and to result from (40), (41) and (43)–(47). Therefore, surface polynomial possesses only two independent coefficients, in this case and , which define ellipse eccentricity. The remaining coefficients must have specific values given by above-mentioned equations.

Force balance equation on direction yields now following relation:

By substituting (49) for central pressure in (39), (42), one finds following expressions of contact ellipse half-axes:

Central or maximum contact pressure yields by substituting (50), for ellipse half-axes in (49):

Finally, (38) leads to following normal approach:

Surface equation is now straightforward: where, excepting the coefficients and which can be chosen arbitrarily, all coefficients result from (40), (41) and (43)–(47).

##### 5.3. Eighth-Order Surfaces

As in previous cases, the coefficients of polynomial surface of equivalent punch leading to flat central pressure result by coefficient identification in (6). This time , meaning an eighth order surface, and (6) becomes where is given by following relation [6]:

Because coefficients are now , and , (54) takes a simpler form, as follows:

The substitution of expressions (21), (37), (55) for , , , and respectively, in (56) and then identification of resulting coefficients yields following relations between surface coefficients and contact parameters , , , , and :

All coefficients depend on the eccentricity of contact ellipse. As in previous cases, this eccentricity is found by aid of two of these coefficients, arbitrary chosen. Following the procedure derived in [6] for homogeneous surfaces, it is convenient to find eccentricity by involving and . Dividing member by member (59)–(61) and (42), following transcendental equation having the eccentricity as unknown results:The solution of (72) is found numerically. Once knowing the eccentricity, the remaining coefficients result from (40), (41) and (63)–(71). Therefore, surface polynomial possesses only two independent coefficients, in this case and , which define ellipse eccentricity. The remaining coefficients must have specified values given by above-mentioned equations.

Force balance equation on direction yields now following correlation between central pressure and applied normal load:

Following equations for half-axes and result from (58), (62), and (73):

Central or maximum contact pressure yields by substituting (75) and (76) for ellipse half-axes in (73):

Finally, (57) leads to following normal approach: Surface equation is now straightforward: where, excepting the coefficients and which can be chosen arbitrarily, all coefficients result from (59) to (61) and (63) to (71).

#### 6. Discussion

An improved pressure distribution between elastic bodies bounded by high-order symmetrical surfaces is a generalized Hertz pressure in which coefficients result from recurrence (12).

Maximum pressure , which occurs now in contact center, relates to average pressure by following equation:

Naturally, if , (80) yields the well-known maximum Hertz pressure. Equations (22), (38), (57) lead to following recurrence relation for :

Contact parameters , , and are expressed as roots of order as follows: where is a coefficient depending on : and , , , and are coefficients which depend both on degree and contact ellipse eccentricity. Once is found, the remaining coefficients result from following equations:

As shown in [7], the newly advanced framework can also be used to get a nearly flat maximum pressure along the line contact. The severe end effects in the contact between parallel cylinders of different lengths can be attenuated by various methods such as partial crowning, hollow ended roller or generatrix profiling, the latter usually based on the logarithmic profile proposed by Lundberg [8]. Reusner [9] evidenced the advantages of a new special logarithmic profile used by SKF, while Teutsch and Sauer [10] advanced a fast method for roller-race contact analysis in roller bearings, based on a theoretical and implicit load-deflection relationship.

In the framework reviewed and extended herein, maximum extension of flat central pressure results if pressure coefficients for a given polynomial order are adequately chosen. Following this idea, a different approach to end effect attenuation in elastic finite length line contact between revolution bodies was advanced in [7], by allowing the roller generatrix to be a polynomial yielding a nearly flat maximum pressure along most of the contact length. A special choice of polynomial coefficients results in a contact width constant along most of the contact length, which can be interpreted as a contact area in a modified line contact. The central region of contact width increases with the degree of the surface polynomial. Eventually, the generatrix of the equivalent rigid roller, found numerically by aid of interference equation, proves to be a high-order even polynomial. For imposed load, contact area extents and degree of the polynomial, the procedure advanced in [7] establish a convenient pressure distribution in line contacts, which in its turn yields to a symmetrical polynomial roller generatrix and a central roller cross radius. While coefficients yielding the equivalent punch surface need to be recomputed every time a perturbation is introduced in the input (i.e., in load, contact half-axes, or surface degree), as they depend explicitly on central pressure and on contact ellipse eccentricity, formulas (12) for hold in all cases.

#### 7. Numerical Validation

The frictionless contact of elastic bodies assimilated to elastic half-spaces can be simulated numerically for an irregular initial clearance using the well-known algorithm based on the conjugate gradient (CG) method advanced by Polonsky and Keer, [11]. The method is fast because the rate of convergence of CG is superlinear and robust because there is mathematical proof of convergence for the CG when the system matrix is symmetrical and positive definite.

The system arising from digitization of geometrical condition of deformation with respect to direction is essentially nonlinear due to presence of rigid-body approach. Early attempts to linearize it resulted in an additional outer loop, in which normal approach was iterated with respect to static force equilibrium. Another difficulty stems from the fact that the contact area, which determines the size of the system having the nodal pressures as unknowns, is also a priori unknown, and a trial-and-error technique is required. Different techniques to overcome these difficulties are overviewed and benchmarked by Allwood, [12].

The algorithm developed by the authors is based on the work reported in [11]. Linearization is achieved by assessing, in every iteration, estimates for the rigid-body translations and rotations, the latter being related to bending moments transmitted through conformal contacts [13]. These estimates are numerically derived through the least square method, as the best fit of an overdetermined system of equations, assembled from equations corresponding to grid cells included in the current contact area. In order to force the solution to verify the static equilibrium equations, a correction of pressure is imposed after every iteration in the CG algorithm. The size of the system, that is, the contact area, is also adjusted at every CG cycle, according to boundary conditions expressed in terms of pressure and gap. Although the algorithm is essentially based on one level of iteration only, a reset of descent directions in the CG minimization process is required every time the contact area changes.

On every iteration of the CG, two convolution type products must be calculated. The order of computation of when using direct multiplication—summation can be reduced significantly to if convolution is performed in the frequency domain, as element-wise multiplication. The errors introduced by problem periodization, implicitly assumed when the spatial series of pressure and influence coefficients are transferred to frequency domain via fast Fourier transform are avoided with a minimal additional computational cost using the Discrete Convolution Fast Fourier Transform (DCFFT) technique advanced by Liu et al. [14]. The computer code obtained using the CG combined with the DCFFT technique is fast enough to solve grids of up to points in the contact area in a reasonable amount of time, allowing for simulation of deterministic rough contact scenarios.

An implementation of this algorithm is used herein to validate the framework for the smooth contact between high-order surfaces. The initial contact clearance fully authorizes the use of half-space assumption, allowing for the use of appropriate Green functions when computing free-surface deflections.

All simulations are performed with a normal load . The data used for validation was generated by imposing a major half-axis of contact ellipse of 1 mm and an aspect ratio of contact ellipse , using appropriate formulas advanced in this work (see Table 1). With these parameters fixed (but otherwise arbitrarily chosen), contact geometry was readily available as input for the numerical program, as all terms can be computed explicitly, using corresponding closed-form relations advanced in this work.

To assure that the discretization error is reduced to an acceptable level, a uniformly spaced rectangular grid, usually assumed for smooth contact scenarios, was imposed in a surface domain exceeding with 20% symmetrically the contact axes predicted by the analytical framework. The imposed precision for pressure convergence was fixed at , according to numerical experimentations reported in [11].

The relative error of pressure distribution was calculated using the following formula: where subscripts are used to index the grid cells in the numerical model, is the elementary patch area, denote the nodal pressures predicted by the numerical program, and discrete pressures computed using the newly advanced closed-form relations, at coordinates matching grid cells control points.

The predictions of the numerical program are found to agree well with the closed-form relations advanced in the theoretical framework, as shown in Tables 2, 3, and 4. The most important errors are found when assessing contact half-axis and can be attributed to grid resolution, as the numerical prediction for this parameter can only vary with integers of grid steps. The analytical and numerical solutions for pressure profiles along contact major half-axis are depicted in Figure 1. Dimensionless pressure is defined as ratio to Hertzian pressure corresponding to the imposed load and contact area.

As shown in [15], a central plateau of uniform pressure can also be found in the elastic-plastic nonconforming contact, when the residual term of displacement and of subsurface stresses becomes significant. The residual print due to permanent deformation of the surface, acting together with the hardening of the elastic-plastic material, decrease the central pressure computed according to the elastic model. The elastic-plastic material responds to loading by developing residual stresses, which decrease stresses induced by pressure, as to oppose further plastic yielding. In the framework proposed herein, the flattened pressure plateau results in purely elastic conditions, and is related only to the fine-tuning of the surface profile. It is clear that the newly advanced indenter will accommodate larger loadings prior to yield inception compared to the quadratic one, when the same contact area is established.

#### 8. Conclusions

Pressure distribution in elastic elliptic contacts between high-order homogeneous surfaces is nonuniform, possessing a local minimum in contact center, maxima in peripheral region, and zero value in contour points, thus reducing the load carrying capacity of the mechanical contact.

A more advantageous pressure distribution must possess a flat central region surrounded by a peripheral monotonous decrease to zero. Such pressure distribution yields from a generalized Hertz pressure if coefficients result from a recurrence equation. As the degree of generalized Hertz pressure increases, the central plateau of pressure extends and maximum pressure decreases, giving overall a better pressure distribution.

The surfaces of equivalent rigid punch generating these favorable pressures when pressed against an elastic half-space are nonhomogeneous and contain all even power terms up to . Analytical expressions of surface coefficients involve , , , , and algebraic combinations of complete elliptical integrals of first and second kind. Only two of these coefficients can be independent to get an elliptical contact area. Although these can be arbitrarily chosen, it is convenient to take and as the independent coefficients yielding contact eccentricity. All remaining parameters result from equations proposed in the paper.

Analytical expressions of contact half-axes and , maximum pressure , and normal approach are derived by similarity to Hertz equations for fourth, sixth, and eighth order surfaces. These are expressed as roots of order of involved parameters and they can be applied directly to calculate contact elements.

General or recurrence equations are established for main contact parameters, which reduce to the results derived in here when , , and . The newly advanced formulas are verified well by a numerical program for the elastic contact with known, but otherwise arbitrarily distributed initial clearance, giving confidence in the proposed framework.

#### Nomenclature

: | Contact domain or contact area; |

: | Typical Hertz square root; |

: | Complete elliptical integrals of first and second order, respectively; |

: | Symbol for the number of combinations; |

: | Major and minor half-axis of contact ellipse, respectively; |

: | Eccentricity of contact ellipse; |

: | Coordinates in the plane bounding the half-space; |

: | Coordinate along normal in initial contact point; |

: | Aspect ratio of contact ellipse; |

: | Elastic constant of the contact; |

: | Elliptic parameter: ; |

, : | Poisson’s ratios and Young moduli, respectively, of contacting materials. |

#### Acknowledgments

This paper was supported by the project “Progress and development through post-doctoral research and innovation in engineering and applied sciences—PRiDE—Contract no. POSDRU/89/1.5/S/57083”, project co-funded from European Social Fund through Sectorial Operational Program Human Resources 2007–2013.

#### References

- M. Glovnea and E. Diaconescu, “Improvement of punch profiles for elastic circular contacts,”
*Journal of Tribology*, vol. 128, no. 3, pp. 486–492, 2006. View at Publisher · View at Google Scholar · View at Scopus - I. Shtaerman,
*Contact Problems in the Theory of Elasticity*, Gostehizdat, Moscow, Russia, 1949, English Translation in 1970 at British Library, FTD-MT-24-61-70. - L. A. Galin,
*Contact Problems in the Theory of Elasticity*, Gostehizdat, Moscow, Russia, 1953. - G. M. L. Gladwell, “Polynomial solutions for an ellipse on an anisotropic elastic half-space,”
*Quarterly Journal of Mechanics and Applied Mathematics*, vol. 31, no. 2, pp. 251–260, 1978. View at Publisher · View at Google Scholar · View at Scopus - F. M. Borodich, “Hertz contact problems for elastic anisotropic halfspace with initial stress,”
*Soviet Applied Mechanics*, vol. 26, no. 2, pp. 126–133, 1990. View at Publisher · View at Google Scholar · View at Scopus - E. N. Diaconescu, “Elliptic elastic contact between high order symmetrical surfaces,”
*Journal of Tribology*, vol. 128, no. 4, pp. 908–914, 2006. View at Publisher · View at Google Scholar · View at Scopus - E. Diaconescu, “A new approach to end effect attenuation in elastic finite length line contacts between revolution bodies,” in
*Proceedings of the STLE/ASME International Joint Tribology Conference (IJTC '08)*, pp. 633–635, Miami, Fla, USA, October 2008, Paper no. IJTC2008-71240. View at Scopus - G. Lundberg, “Elastische Berührung zweier Halbräume,”
*Forschung auf dem Gebiete des Ingenieurwesens*, vol. 10, no. 5, pp. 201–211, 1939. View at Publisher · View at Google Scholar · View at Scopus - H. Reusner, “The logarithmic roller profile - the key to superior performance of cylindrical and taper roller bearings,”
*Ball Bearing Journal*, no. 230, pp. 2–10, 1987. View at Google Scholar · View at Scopus - R. Teutsch and B. Sauer, “An alternative slicing technique to consider pressure concentrations in non-Hertzian line contacts,”
*Journal of Tribology*, vol. 126, no. 3, pp. 436–442, 2004. View at Publisher · View at Google Scholar · View at Scopus - I. A. Polonsky and L. M. Keer, “A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques,”
*Wear*, vol. 231, no. 2, pp. 206–219, 1999. View at Publisher · View at Google Scholar · View at Scopus - J. Allwood, “Survey and performance assessment of solution methods for elastic rough contact problems,”
*Journal of Tribology*, vol. 127, no. 1, pp. 10–23, 2005. View at Publisher · View at Google Scholar · View at Scopus - S. Spînu and E. Diaconescu, “Numerical simulation of elastic conforming contacts under eccentric loading,” in
*Proceedings of the STLE/ASME International Joint Tribology Conference (IJTC '08)*, pp. 649–651, Miami, Fla, USA, October 2008. View at Scopus - S. Liu, Q. Wang, and G. Liu, “A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses,”
*Wear*, vol. 243, no. 1-2, pp. 101–111, 2000. View at Publisher · View at Google Scholar · View at Scopus - S. Spinu, G. Frunza, and E. Diaconescu, “Numerical simulation of elastic-plastic non-conforming contact,” in
*Numerical Simulations—Applications, Examples and Theory*, L. Angermann, Ed., pp. 253–280, InTech, 2011. View at Google Scholar