Local stress approach

FE based failure probability calculation

FE Data

we are using VMAP data format and rst file formats. It is also possible to use odb data,

import numpy as np
import pandas as pd
import pickle
import pylife.vmap
import pylife.mesh
import pylife.mesh.meshsignal
import pylife.stress.equistress
import pylife.stress
import pylife.strength.fatigue
import pylife.utils.histogram as psh
import pyvista as pv

VMAP

For plotting of VMAP data we are using pyVista.

pyLife_mesh = (pylife.vmap.VMAPImport("plate_with_hole.vmap").make_mesh('1', 'STATE-2')
               .join_coordinates()
               .join_variable('STRESS_CAUCHY')
               .to_frame())

pyLife_mesh['mises'] = pyLife_mesh.equistress.mises()
grid = pv.UnstructuredGrid(*pyLife_mesh.mesh.vtk_data())
plotter = pv.Plotter()
plotter.add_mesh(grid, scalars=pyLife_mesh.groupby('element_id')['mises'].mean().to_numpy(),
                show_edges=True, cmap='jet')
plotter.add_scalar_bar()
plotter.show()

Now we want to apply the collectives to the mesh

mises = pyLife_mesh.groupby('element_id')[['S11', 'S22', 'S33', 'S12', 'S13', 'S23']].mean().equistress.mises()
mises /= mises.max()  # the nominal load level in the FEM analysis is set, that s_max = 1
collectives = pickle.load(open("collectives.p", "rb"))
collectives = collectives.unstack().T.fillna(0)
collectives_sorted = psh.combine_histogram([collectives[col] for col in collectives],
                                             method="sum")

scaled_collectives = collectives_sorted.load_collective.scale(mises)
display(scaled_collectives.to_pandas().sample(5))

           range                                     element_id
frequency  (78.65448616078119, 88.48629693087884]    3469           136.0
           (22.325070973361967, 26.790085168034363]  2703          5660.0
           (11.075553710196642, 22.151107420393284]  2169             0.0
           (36.96782896221106, 43.12913378924623]    641           1654.0
           (40.848997892571006, 49.0187974710852]    208           1850.0
Name: cycles, dtype: float64

Define the material parameters

mat = pd.Series({
    'k_1': 8.,
    'ND': 1.0e6,
    'SD': 200.0, # amplitude of the stress signal
    'TN': 1./12.,
    'TS': 1./1.1
})

Damage Calculation

damage = mat.fatigue.miner_haibach().damage(scaled_collectives)
print("Max damage : %f" % damage.max())
damage = damage.groupby(['element_id']).sum()

Max damage : 0.001484

grid = pv.UnstructuredGrid(*pyLife_mesh.mesh.vtk_data())
plotter = pv.Plotter()
plotter.add_mesh(grid, scalars=damage.to_numpy(),
                show_edges=True, cmap='jet')
plotter.add_scalar_bar()
plotter.show()

ANSYS

#%% Ansys (license is necessary)
# For Ansys  *.rst files we are using pymapdl
# from ansys.mapdl import reader as pymapdl_reader
# # for more information please go to pymapdl
# # rst_input = post.load_solution("beam_3d.rst")
# # # pymapdl has some nice features
# # rst_input.plot_nodal_displacement(0)
# # rst_input.plot_nodal_stress(0,"X")
# ansys_mesh = pymapdl_reader.read_binary('beam_3d.rst')
# grid_ansys = ansys_mesh.grid
# plotter = pv.Plotter(window_size=[1920, 1080])
# _, volume, _  = ansys_mesh.element_solution_data(0,"ENG")
# volume = pd.DataFrame(volume)[1]

# nodes, ansys_mesh_mises = ansys_mesh.nodal_stress(0)
# ansys_mesh_mises = pd.DataFrame(data = ansys_mesh_mises,
#                                 columns=['S11', 'S22', 'S33', 'S12', 'S13', 'S23']).equistress.mises()


# test = pd.DataFrame(ansys_mesh.mesh.elem).iloc[:, 8:]
# #%%
# plotter.add_mesh(grid_ansys, scalars=ansys_mesh_mises,
#                 show_edges=True, cmap='jet')
# plotter.add_scalar_bar()
# plotter.show()