Derived stress quantities

CalculiX writes a six-component stress tensor. It does not write von Mises stress or the principal values, because those are functions of what is already there. This reader computes them on the way through, under the names PyVista’s own reader uses, so a script written against one works against the other.

For any array whose name contains STRESS or STRAIN and which has six components, five arrays are appended.

from pathlib import Path

import numpy as np
import pyvista as pv

import pyvista_frd

# Sphinx-Gallery runs each example from its own directory; the decks and
# their solved output live one level up, in ``doc/_data``.
DATA = Path('../_data').resolve()
mesh = pyvista_frd.read(DATA / 'cantilever.frd')

The tensor and what comes from it.

print('stored :', mesh.point_data['STRESS'].shape, '(xx, yy, zz, xy, yz, zx)')
for name in ('STRESS_Mises', 'STRESS_sgMises', 'STRESS_PS1', 'STRESS_PS2', 'STRESS_PS3'):
    print(f'derived: {name:16s} {mesh.point_data[name].shape}')

stored : (900, 6) (xx, yy, zz, xy, yz, zx)
derived: STRESS_Mises     (900,)
derived: STRESS_sgMises   (900,)
derived: STRESS_PS1       (900,)
derived: STRESS_PS2       (900,)
derived: STRESS_PS3       (900,)

_Mises is the equivalent stress and _sgMises the same magnitude carrying the sign of the trace, which separates tension from compression in a single scalar. The cantilever is in bending, so the two faces disagree.

warped = mesh.warp_by_vector('DISP', factor=200)

pl = pv.Plotter(shape=(2, 1), window_size=(1000, 640))
for row, (name, cmap) in enumerate((('STRESS_Mises', 'inferno'), ('STRESS_sgMises', 'coolwarm'))):
    pl.subplot(row, 0)
    pl.add_text(name, font_size=10)
    pl.add_mesh(warped.copy(), scalars=name, cmap=cmap)
    pl.camera.tight(padding=0.12, adjust_render_window=False, view='xz')
pl.show()

The principal values are the eigenvalues of the tensor, ascending, so PS1 is the largest and PS3 the smallest. In a bent beam the largest principal stress picks out the tension side.

pl = pv.Plotter(shape=(3, 1), window_size=(1000, 840))
for row, name in enumerate(('STRESS_PS1', 'STRESS_PS2', 'STRESS_PS3')):
    pl.subplot(row, 0)
    pl.add_text(name, font_size=9)
    pl.add_mesh(warped.copy(), scalars=name, cmap='coolwarm')
    pl.camera.tight(padding=0.12, adjust_render_window=False, view='xz')
pl.show()

These are not approximations of the solver’s own numbers – CalculiX never wrote them. They are computed here, and the conformance suite requires the von Mises arrays to be bit-identical to PyVista’s, which is why the expression order in the C++ core is fixed rather than merely equivalent.

The principal values come from an eigensolver, where bit-identity is not an achievable goal, so those are held to a stated band instead. doc/ divergences.md records the band and why it is measured against the magnitude of the tensor rather than of the eigenvalue.

mises = mesh.point_data['STRESS_Mises']
tensor = mesh.point_data['STRESS']
trace = tensor[:, :3].sum(axis=1)
print(f'peak von Mises      {mises.max():10.3f}')
print(f'nodes in tension    {int(np.count_nonzero(trace > 0)):10d} of {len(trace)}')

peak von Mises           5.980
nodes in tension           450 of 900