Visualize Hertzian Contact Stress#

The following example demonstrates how to use PyVista to visualize Hertzian contact stress between a cylinder and a flat plate.

This example loads a dataset, constructs a line to represent the point of contact between the cylinder and the block, and samples the stress along that line. Finally, it plots the dataset and the stress distribution.

Background Hertzian contact stress refers to the stress that occurs between two curved surfaces that are in contact with each other. It is named after Heinrich Rudolf Hertz, a German physicist who first described the phenomenon in the late 1800s. Hertzian contact stress is an important concept in materials science, engineering, and other fields where the behavior of materials under stress is a critical consideration.

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

import pyvista as pv
from pyvista import examples

Load the dataset#

Start by loading the dataset using pyvista.examples module. This module provides access to a range of datasets, including FEA (finite element analysis) datasets that are useful for stress analysis.

mesh = examples.download_fea_hertzian_contact_cylinder()
mesh
HeaderData Arrays
UnstructuredGridInformation
N Cells132258
N Points34185
X Bounds-0.000e+00, 2.000e-01
Y Bounds0.000e+00, 2.500e-02
Z Bounds0.000e+00, 3.000e-01
N Arrays16
NameFieldTypeN CompMinMax
DisplacementPointsfloat643-5.000e-042.945e-05
vonMisesPointsfloat6411.316e+021.707e+09
StressPointsfloat646-3.055e+095.274e+08
StrainPointsfloat646-1.937e-021.894e-02
PrincipalStress 1Pointsfloat641-1.787e+092.586e+08
PrincipalStress 2Pointsfloat641-2.594e+093.531e+07
PrincipalStress 3Pointsfloat641-3.078e+091.995e+06
PrincipalStrain 1Pointsfloat641-6.976e-041.448e-02
PrincipalStrain 2Pointsfloat641-1.145e-028.690e-04
PrincipalStrain 3Pointsfloat641-2.787e-02-8.375e-10
StrainEnergyDensityPointsfloat6411.172e-072.307e+07
PlasticStrainPointsfloat6460.000e+000.000e+00
EquivalentPlasticStrainPointsfloat6410.000e+000.000e+00
RankCellsfloat6410.000e+001.500e+01
MaterialCellsfloat6410.000e+001.000e+00
PartIDCellsint3211.000e+002.000e+00


Plot the Dataset#

Plot the dataset by part ID.

mesh.plot(scalars='PartID', cmap=['green', 'blue'], show_scalar_bar=False)
fea hertzian contact pressure

Creating a Line to Denote the Point of Contact#

Construct a line to represent the point of contact between the cylinder and the plate.

ypos = 0.024
a = [0.1, ypos, 0.2 - 1e-4]
b = [0.095, ypos, 0.2 - 1e-4]
line = pv.Line(a, b, resolution=100)
line.clear_data()
line
PolyDataInformation
N Cells1
N Points101
N Strips0
X Bounds9.500e-02, 1.000e-01
Y Bounds2.400e-02, 2.400e-02
Z Bounds1.999e-01, 1.999e-01
N Arrays0


Sampling the Stress along the Line#

We can sample the Z component stress along the contact edge and compare it with expected pressure.

The expected values array is the Hertzian contact pressure and is the analytical solution to the non-adhesive contact problem. Computation of these values is an exercise left up to the reader (the radius of the cylinder is 0.05). See Contact Mechanics

# Sample the stress
sampled = line.sample(mesh, tolerance=1e-3)
x_coord = 0.1 - sampled.points[:, 0]
samp_z_stress = -sampled['Stress'][:, 2]

# Expected Hertzian contact pressure
h_pressure = np.array(
    [
        [0.0000, 1718094092],
        [0.0002, 1715185734],
        [0.0004, 1703502649],
        [0.0006, 1683850714],
        [0.0008, 1655946243],
        [0.001, 1619362676],
        [0.0012, 1573494764],
        [0.0014, 1517500856],
        [0.0016, 1450208504],
        [0.0018, 1369953775],
        [0.002, 1274289906],
        [0.0022, 1159408887],
        [0.0024, 1018830677],
        [0.0026, 839747409.8],
        [0.0028, 587969605.2],
        [0.003, 0],
        [0.005, 0],
    ],
)

plt.plot(x_coord, samp_z_stress, '.', label='Z Component Stress')
plt.plot(h_pressure[:, 0], h_pressure[:, 1], label='Hertzian contact pressure')
plt.legend()
plt.show()
fea hertzian contact pressure

Visualizing the Z Stress Distribution#

You can now visualize the Z stress distribution. Use pyvista.Plotter to create a plot window and add the dataset to it.

pl = pv.Plotter()
z_stress = np.abs(mesh['Stress'][:, 2])
pl.add_mesh(
    mesh,
    scalars=z_stress,
    clim=[0, 1.2e9],
    cmap='gouldian',
    scalar_bar_args={'title': 'Z Component Stress (Pa)', 'color': 'w'},
    lighting=True,
    show_edges=False,
    ambient=0.2,
)
pl.camera_position = 'xz'
pl.set_focus(a)
pl.camera.zoom(2.5)
pl.show()
fea hertzian contact pressure

Total running time of the script: (0 minutes 2.590 seconds)

Gallery generated by Sphinx-Gallery