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

Options
Ayush Kumar
Ayush Kumar Member, Moderator, Employee Posts: 416
First Anniversary Ansys Employee Solution Developer Community of Practice Member First Answer
edited June 2023 in Structures
  • How do I calculate linearized maximum principal stress (Membrane, Bending) using PyDPF?
Tagged:

Comments

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 416
    First Anniversary Ansys Employee Solution Developer Community of Practice Member First Answer
    edited May 2023
    Options
    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])