How do I calculate linearized maximum principal stress (Membrane, Bending) using PyDPF?

Ayush Kumar

• How do I calculate linearized maximum principal stress (Membrane, Bending) using PyDPF?

Ayush Kumar

1. Linearize the component stresses.
2. Calculate maximum principal from the linearized component stresses.

import matplotlib.pyplot as plt
 
from ansys.dpf import core as dpf
from ansys.dpf.core import examples
 
# model = dpf.Model(examples.find_static_rst())
model = dpf.Model(r"\Path\to\file.rst")
nid_1 = 465
nid_2 = 673
 
mesh = model.metadata.meshed_region
n1 = mesh.nodes.node_by_id(nid_1)
n2 = mesh.nodes.node_by_id(nid_2)
 
x_0 = n1.coordinates[0]
y_0 = n1.coordinates[1]
z_0 = n1.coordinates[2]
 
x_1 = n2.coordinates[0]
y_1 = n2.coordinates[1]
z_1 = n2.coordinates[2]
 
path_length = ((x_1 - x_0) ** 2 + (y_1 - y_0) ** 2 + (z_1 - z_0) ** 2) ** 0.5
 
n_points = 49  # A linearized stress has fixed a number of 48 points.
delta = path_length / (n_points - 1)
 
line_unit_vector = [(x_1 - x_0) / path_length, (y_1 - y_0) / path_length, (z_1 - z_0) / path_length]
 
# Line equation
fx = lambda t: x_0 + line_unit_vector[0] * t
fy = lambda t: y_0 + line_unit_vector[1] * t
fz = lambda t: z_0 + line_unit_vector[2] * t
 
coordinates = [[fx(t * delta), fy(t * delta), fz(t * delta)] for t in range(n_points)]
flat_coordinates = [entry for data in coordinates for entry in data]
 
field_coord = dpf.fields_factory.create_3d_vector_field(n_points)
field_coord.data = flat_coordinates
field_coord.scoping.ids = list(range(1, n_points + 1))
 
# Stress Tensor
s = model.operator("S")
s.inputs.requested_location.connect("Nodal")
s_f = s.outputs.fields_container.get_data()
 
mapping_operator = dpf.operators.mapping.on_coordinates(
    fields_container=s_f,
    coordinates=field_coord,
    create_support=True,
    mesh=mesh
)
fields_mapped = mapping_operator.outputs.fields_container.get_data()
 
# Membrane Stress
membrane_stress = (fields_mapped[0].get_entity_data(0) / 2 + fields_mapped[0].get_entity_data(48) / 2 + sum(
    fields_mapped[0].data[1:-1])) / (n_points - 1)
 
# Bending stress
path_1 = -1 * path_length / 2
path_n = path_length / 2
path_range = [path_1 + delta * i for i in range(n_points)]
path_range_field = dpf.fields_factory.create_scalar_field(n_points, location=dpf.locations.nodal)
path_range_field.data = path_range
path_range_field.scoping.ids = range(1, 50)
 
# Function to be integrated
stress_scaled = dpf.operators.math.scale_by_field(fieldA=fields_mapped[0], fieldB=path_range_field)
stress_scaled_data = list(stress_scaled.outputs.field.get_data().data)
 
# Use extended Simpson's rule for Numerical Integration of `stress_scaled_data`
stress_scaled_integral = (17 * stress_scaled_data[0] + 59 * stress_scaled_data[1] + 43 * stress_scaled_data[2] + 49 *
                          stress_scaled_data[3] + 48 * sum(stress_scaled_data[4:-4]) + 49 * stress_scaled_data[
                              n_points - 4] + 43 * stress_scaled_data[n_points - 3] + 59 * stress_scaled_data[
                              n_points - 2] + 17 * stress_scaled_data[n_points - 1]) / 48.0
 
 
# Bending Stress at Node N1
b1 = stress_scaled_integral * (-6.0 / path_length) / 48.0
 
# Bending stress tensor @ N1
b1_f = dpf.fields_factory.create_tensor_field(1, "Nodal")
b1_f.data = b1
 
# Membrane stress tensor
m_f = dpf.fields_factory.create_tensor_field(1, "Nodal")
m_f.data = membrane_stress
 
# Calculate Membrane S1
s1_m = dpf.operators.invariant.principal_invariants()
s1_m.inputs.field.connect(m_f)
s1_m_val = s1_m.outputs.field_eig_1.get_data().data
print("Membrane stress S1 - %s" % s1_m_val[0])
 
# Calculate Bending S1 @ N1
b1_m = dpf.operators.invariant.principal_invariants()
b1_m.inputs.field.connect(b1_f)
b1_m_val = b1_m.outputs.field_eig_1.get_data().data
print("Bending stress S1 - %s" % b1_m_val[0])