Revision as of 13:02, 9 May 2017

Click on Solid Mechanics to go back.

Contact of Cylinders - the Hertz problem

Detailed discussions of this problem can be found in Hills and Nowells (1994) as well as Williams and Dwyer-Joyce (2001). ^{[1] [2]}

If two circular cylinders with radii $R_1$ and $R_2$ are pressed together by a force per unit length of magnitude $P$ with their axes parallel, then the contact patch will be of half-width $b$ such that

The resulting pressure distribution $p(x)$ is semielliptical, i.e., of the form

The coordinate $x$ is measured perpendicular to that of the cylinder axes. The cylinder is pressing on the lower halfplane from above. Analytical solutions for stresses in the plane are presented below. The coordinate ranges are $x \in (-\infty, \infty), z\in (-\infty, 0]$.

The surface stresses are given in the following equations. At the contact interface

Note that these stresses are all independent of Poisson's ratio although, for plane strain, the third principal stress

At a general point $(x,z)$ the stresses may be expressed in terms of $m$ and $n$, defined by

Plots of this analytical solution are presented below, using $p_0 = 8.36\cdot10^{7}, b = 4.13\cdot10^{-4}$

The solutions were verified and graphs were generated using this mathematica notebook: File:hertzian_analytical.nb Copy-pasteable Matlab code:

 1 function [sxx, syy, sxy] = analytical(x, z, b, p0)
 2 % Returns the analytical solution for stress for hertzian contact.
 3 % x in (-inf, inf), z in (-inf, 0].
 4 bxz = b^2 - x.^2 + z.^2;
 5 xz = 4*x.^2.*z.^2;
 6 koren = sqrt(bxz.^2 + xz);
 7 m2 = 1/2 * (koren + bxz);
 8 n2 = 1/2 * (koren - bxz);
 9 m = sqrt(m2);
10 n = sign(x) .* sqrt(n2);
11 mpn = m2 + n2;
12 zmn = (z.^2 + n2) ./ mpn;
13 sxx = -p0 / b * (m .* (1 + zmn) + 2*z);
14 syy = -p0 / b * m .* (1 - zmn);
15 sxy = p0 / b * n .* ((m2 - z.^2) ./ mpn);
16 end

The extreme of $\tau_1$ is achieved at $x = 0$, $z = -\sqrt{\frac{1}{2} \left(\sqrt{5}-1\right)} b = -0.786b$.

When observing for $z \to -\infty$ it holds that $s_{xx} = \frac{\frac{b^3 p_0}{4}+b p_0 x^2}{z^3}+O(z^{-5})$, $s_{zz} = \frac{b p_0}{z}-\frac{b p_0 \left(b^2+4 x^2\right)}{2 z^3}+O(z^{-5})$, $s_{xz} = \frac{b p_0 x}{z^2}+\frac{-\frac{3}{2} b^3 p_0 x-2 b p_0 x^3}{z^4}+O(z^{-5})$.

When observing for $x \to \infty$ (and similarly for $-\infty$) it holds that

$s_{xx} = \frac{b p_0 z}{x^2}+\frac{\frac{3}{4} b^3 p_0 z-2 b p_0 z^3}{x^4}+O(x^{-5})$, $s_{zz} = \frac{b p_0 z^3}{x^4}+O(x^{-5})$, $s_{xz} = \frac{b p_0 z^2}{x^3}+O(x^{-5})$.

Contact of cylinders under partial slip

The second case we study is the application of a tangential force $Q$ to the previous problem. When the tangential force is less than the limiting force of friction, i.e.,

Besides the normal traction $p(x)$ we know have an additional shear traction given by

The distributions $q(x)$ and $p(x)$ as well as the widths of the stick and slip zones can be seen in the image below.

The effect of bulk stress

Additionally we might be interested in the addition of bulk stress. This type of stress occurs in fretting fatique experiments like the one shown below.

500px

The previous solution for contact of cylinders under partial slip can be adjusted for the presence of bulk stresses $\sigma_\mathrm{axial}$. These cause an eccentricity $e$ to the solution given above. The shear traction $q(x)$ can be written as:

FreeFem++ numerical solution

For the numerical solution in FreeFem++ we choose parameters similar to those in Pereira et al. (2016): ^[3]

• Modulus of elasticity: $E = 72.1$ GPa
• Poisson's ratio: $\nu = 0.33$
• Normal load: $P = 543$ N
• Coefficient of friction: $\mu = 0.85$
• Cylinder radius: $R = 50$ mm
• Specimen length: $L = 40$ mm
• Specimen height: $H = 10$ mm

We assume that both the pad and specimen are from the same material, therefore the combined modulus is

(this might not be correct, because $p_0$ should be equal to $2P/(\pi b)$? Does $R/H$ perhaps represent a dimensionless radius?)

Out of simplicity we only use traction boundaries at the top surface. On the left, right and bottom sides of the specimen Dirichlet boundaries are used. The center of the coordinate system is placed at $(x,y) = (L/2,0)$. The axises therefore go from $-L/2$ to $L/2$ in $x$-direction and from $0$ to $H$ in $y$-direction.

                top
      ________________________
     |                        |
     |                        |
left |                        | right
     |                        |
     |________________________|
                
               bottom

Boundary conditions:

The numerical solution is obtained for plane strain conditions using quadratic finite elements. The final mesh after adaptation had ~$40000$ dof. The images below show the raw values of the surface stresses, sub-surface stresses along the axis of symmetry and a close-up of the maximum shear stress contours close to the contact. The contact width $b \approx 0.29$ mm for the parameters described at the beginning of this section.

MLSM numerical solution for halfplane contact

Problem parameters used:

• Modulus of elasticity: $E = 72.1$ GPa
• Poisson's ratio: $\nu = 0.33$
• Normal load: $P = 543$ N
• Coefficient of friction: $\mu = 0.85$
• Cylinder radius: $R = 50$ mm
• Specimen length: $L = 40$ mm
• Specimen height: $H = 20$ mm

The boundary conditions are zero except for the top layer which has zero traction in $x$ direction and traction $t_z$ in $z$ direction. The solution matches the analytical solution well. 800px

800px

800px

Execution time is shown below. Program ran on a single core. Most of the time is spent on solving the system of the equations.

The error decreases as the density of nodes increases, up to some limit when the system becomes too unstable.

800px

MLSM numerical solution for bulk load

Case definition

We are solving the equation

• top: $\vec{t}(x) = (q(x), -p(x))$ or componentwise $(2 \mu+\lambda) \frac{\partial u}{\partial x}(x, 0) + \lambda \frac{\partial v}{\partial y}(x, 0) = q(x)$ and $\mu \frac{\partial u}{\partial y}(x, 0) + \mu \frac{\partial v}{\partial x}(x, 0) = -p(x)$.
• left: $\b{u} = 0$ or componentwise $u(-L/2, y) = v(-L/2, y) = 0$.
• bottom: up-down symmetry conditions: componentwise $\frac{\partial u}{\partial y}(x, -H/2) = 0, v(x, -H/2) = 0$.
• right: this part is traction free, ie. $\vec{t}(y) = 0$ or componentwise $\mu \frac{\partial u}{\partial y}(L/2, y) + \mu \frac{\partial v}{\partial x}(L/2, y) = 0$ and $(2 \mu+\lambda) \frac{\partial v}{\partial y}(L/2, y) + \lambda \frac{\partial u}{\partial x}(L/2, y) = 0$.

Phsyical parameters

Basic parameters are:

• $E = \unit{72.1}{GPa}$, Youngs modulus
• $\nu = 0.33$, Poissons ratio
• $L = \unit{40}{mm}$, length, width of the pad
• $H = \unit{4}{mm}$, height, thickness of the pad
• $\sigma_{ax} = \unit{100}{MPa}$, axial pressure
• $F = \unit{543}{N}$, normal force
• $Q = \unit{155}{N}$, tangential force
• $R = \unit{10}{mm}$ or $\unit{50}{mm}$, raduis of curvature of the cylindrical pad
• $COF = \mu = 0.3$ or $0.85$ or $2$, coefficient of friction

Derived parameters are, for choice of $R = \unit{0.1}{mm}$ and $\mu = 0.3$:

• $E^\ast = \frac{E}{2(1-\nu^2)} = \unit{40.4}{GPa}$, combined Young's modulus
• $p_{max} = p_0 = \sqrt{\frac{FE^\ast}{HR \pi}} = \unit{418.10407}{MPa}$, maximal normal pressure
• $a = 2 \sqrt{\frac{FR}{HE^\ast \pi}} = \unit{0.2067}{mm}$, semi contact width, contact region is $[-a, a] \times \{0\}$
• $c = a \sqrt{1 - \frac{Q}{\mu F}} = \unit{0.04504}{mm}$, stick zone semi width
• $e = \frac{a \sigma_{ax}}{4 \mu p_{max}} = \unit{0.041197}{mm}$, eccentricity due to axial load

References

1. Jump up ↑ Hills, D. A. and Nowell, D. (1994). Mechanics of Fretting Fatique, p. 20-25. Springer Science+Business Media, Dordrecht.
2. Jump up ↑ Williams, John A. and Dwyer-Joyce, Rob S. (2001). Contact Between Solid Surfaces, p. 121 in Modern Tribology Handbook: Volume 1, Principles of Tribology, editor: Bushan, Bharat. CRC Press LLC, Boca Raton.
3. Jump up ↑ Pereira et al. (2016). On the convergence of stresses in fretting fatique. Materials, 9, 639.