How do I calculate stress along a path of desired length normal to a surface at a node?

Options
Ayush Kumar
Ayush Kumar Member, Moderator, Employee Posts: 414
First Anniversary Ansys Employee Solution Developer Community of Practice Member First Answer
edited June 2023 in Structures

How do I calculate stress along a path of desired length normal to a surface at a node?

Tagged:

Best Answer

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 414
    First Anniversary Ansys Employee Solution Developer Community of Practice Member First Answer
    Answer ✓
    Options

    With ACT you can only define path through the complete length of the body of the chosen node and interpolate to get value a certain point. With PyDPF can define path of desired length.

    ACT - define a named selection face and node. Node where you want the path to start and the face to which the node belongs.

    # Create face and node named selections
    geom_id = ExtAPI.DataModel.GetObjectsByName("face")[0].Ids[0]
    node_id = ExtAPI.DataModel.GetObjectsByName("node")[0].Ids[0]
    
    mesh_data = DataModel.MeshDataByName(DataModel.MeshDataNames[0])
    my_face = DataModel.GeoData.GeoEntityById(geom_id)
    mesh_node = mesh_data.NodeById(node_id)
    base_point = mesh_node.X, mesh_node.Y, mesh_node.Z
    
    # Get Normal vector at the Node
    u, v = my_face.ParamAtPoint(base_point)
    vec_normal_out = my_face.NormalAtParam(u, v)  # Normal Vector to the face at given node
    vec_normal_in = [-vec_normal_out[0], -vec_normal_out[1], -vec_normal_out[2]]
    
    vec_global = [1, 0, 0]  # Global X axis
    
    
    # Cross Product function definition
    def cross_product(vec_1, vec_2):
        vec_cross_product = [vec_1[1] * vec_2[2] - vec_1[2] * vec_2[1],
                             vec_1[2] * vec_2[0] - vec_1[0] * vec_2[2],
                             vec_1[0] * vec_2[1] - vec_1[1] * vec_2[0]]
        return vec_cross_product
    
    
    # Dot product of two vectors
    dot_product = lambda vec_1, vec_2: sum(map(lambda x, y: x * y, vec_1, vec_2))
    
    
    # Vector magnitude
    def magnitude(vec):
        return (vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2) ** 0.5
    
    
    # Euler axis - angle
    rot_rad = acos(dot_product(vec_global, vec_normal_in) / magnitude(vec_global) * magnitude(vec_normal_in))
    # Euler axis - Axis of rotation
    axis = cross_product(vec_global, vec_normal_in)
    # Unit vector normal to Surface and Global X axis vector
    axis_unit_vec = [axis[0] / magnitude(axis), axis[1] / magnitude(axis), axis[2] / magnitude(axis)]
    
    
    def get_euler_angles_from_axis(euler_axis, rot_angle):
        # Euler axis–angle ↔ quaternions
        q1 = euler_axis[0] * sin(rot_angle / 2)
        q2 = euler_axis[1] * sin(rot_angle / 2)
        q3 = euler_axis[2] * sin(rot_angle / 2)
        q4 = cos(rot_angle / 2)
    
        # Quaternion → Euler angles (z-y′-x″ intrinsic)
        roll = atan2((2 * q4 * q1 + 2 * q2 * q3), (1 - 2 * (q1 ** 2 + q2 ** 2)))
        pitch = asin(2 * (q4 * q2 - q1 * q3))
        yaw = atan2((2 * q4 * q3 + 2 * q1 * q2), (1 - 2 * (q2 ** 2 + q3 ** 2)))
    
        # Euler angles (degrees)
        roll_deg = roll * 180 / pi
        pitch_deg = pitch * 180 / pi
        yaw_deg = yaw * 180 / pi
    
        return yaw_deg, pitch_deg, roll_deg
    
    
    yaw_deg, pitch_deg, roll_deg = get_euler_angles_from_axis(axis_unit_vec, rot_rad)
    
    # New coordinate system creation and orientation change to match surface normal vector
    new_cs = ExtAPI.DataModel.Project.Model.CoordinateSystems.AddCoordinateSystem()
    selection = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshNodes)
    selection.Ids = [mesh_node.Id]
    new_cs.OriginLocation = selection
    new_cs.RotateZ()
    new_cs.RotateY()
    new_cs.RotateX()
    new_cs.SetTransformationValue(1, yaw_deg)
    new_cs.SetTransformationValue(2, pitch_deg)
    new_cs.SetTransformationValue(3, roll_deg)
    
    # Creation of new path based on oriented coordinate system & X-axis intersection method
    new_path = ExtAPI.DataModel.Project.Model.AddConstructionGeometry().AddPath()
    new_path.PathType = PathScopingType.Ray
    new_path.PathCoordinateSystem = new_cs
    new_path.NumberOfSamplingPoints = 100
    
    # Add Stress object
    vms = Model.Analyses[0].Solution.AddEquivalentStress()
    vms.ScopingMethod = GeometryDefineByType.Path
    vms.Location = new_path
    vms.EvaluateAllResults()
    
    vms.Activate()
    
    paneTabular=ExtAPI.UserInterface.GetPane(MechanicalPanelEnum.TabularData)
    control = paneTabular.ControlUnknown
    col = control.ColumnsCount
    row = control.RowsCount
    abcissa = [float(control.cell(row_nr, 2).Text) for row_nr in range(2, row + 1)]
    ordinate = [float(control.cell(row_nr, 3).Text) for row_nr in range(2, row + 1)]
    
    def linear_interpolation(p1, p2, x):
        slope = (x - p1[0]) / (p2[0] - p1[0])
        p_int = p1[1] + (p2[1] - p1[1]) * slope
        return p_int
        
    at_x = 0.6
    temp_val = filter(lambda x: (x > at_x), abcissa)[0]
    index = abcissa.IndexOf(temp_val)
    
    p1 = abcissa[index - 1], ordinate[index - 1]
    p2 = abcissa[index], ordinate[index]
    stress_at_x = linear_interpolation(p1, p2, at_x)
    
    

    PyDPF

    """
    .. _stress_gradient_path:
    
    Stress gradient normal to a defined node.
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    This example shows how to plot a stress gradient normal to a selected node.
    A path is created of a defined length.
    
    """
    
    ###############################################################################
    # First, import the DPF-Core module as ``dpf_core`` and import the
    # included examples file and ``DpfPlotter``
    #
    import matplotlib.pyplot as plt
    from ansys.dpf import core as dpf
    from ansys.dpf.core import operators as ops
    from ansys.dpf.core.plotter import DpfPlotter
    from ansys.dpf.core import examples
    
    ###############################################################################
    # Next, open an example and print out the ``model`` object.  The
    # :class:`Model <ansys.dpf.core.model.Model> class helps to organize access
    # methods for the result by keeping track of the operators and data sources
    # used by the result
    # file.
    #
    # Printing the model displays:
    #
    # - Analysis type
    # - Available results
    # - Size of the mesh
    # - Number of results
    # - Unit
    #
    model = dpf.Model(examples.hemisphere)
    print(model)
    ###############################################################################
    # Define the `node_id` normal to which a stress gradient should be plotted.
    #
    node_id = 1928
    ###############################################################################
    # The following command prints the mesh unit
    #
    unit = model.metadata.meshed_region.unit
    print("Unit: %s" % unit)
    ###############################################################################
    # `depth` defines the path length / depth to which the path will penetrate.
    # While defining `depth` make sure you use the correct mesh unit.
    # `delta` defines distance between consecutive points on the path.
    depth = 10  # in mm
    delta = 0.1  # in mm
    ###############################################################################
    # Get the meshed region
    #
    mesh = model.metadata.meshed_region
    ###############################################################################
    # Get Equivalent stress fields container.
    #
    stress_fc = model.results.stress().eqv().outputs.fields_container()
    ###############################################################################
    # Define Nodal scoping.
    # Make sure to define ``"Nodal"`` as the requested location, important for the
    # `normals` operator.
    #
    nodal_scoping = dpf.Scoping(location="Nodal")
    nodal_scoping.ids = [node_id]
    ###############################################################################
    # Get Skin Mesh because `normals` operator requires Shells as input.
    #
    skin_mesh = ops.mesh.skin(mesh=mesh)
    skin_meshed_region = skin_mesh.outputs.mesh.get_data()
    ###############################################################################
    # Get normal at a node using `normals` operator.
    #
    normal = ops.geo.normals()
    normal.inputs.mesh.connect(skin_meshed_region)
    normal.inputs.mesh_scoping.connect(nodal_scoping)
    normal_vec_out_field = normal.outputs.field.get_data()
    ###############################################################################
    # Normal vector is along the surface normal. We need to invert the vector
    # using `math.scale` operator inwards in the geometry, to get the path
    # direction.
    #
    normal_vec_in_field = ops.math.scale(field=normal_vec_out_field,
                                         ponderation=-1.0)
    normal_vec_in = normal_vec_in_field.outputs.field.get_data().data[0]
    ###############################################################################
    # Get Nodal coordinates, they serve as the first point on the line.
    #
    node = mesh.nodes.node_by_id(node_id)
    line_fp = node.coordinates
    ###############################################################################
    # Create 3D line equation.
    #
    fx = lambda t: line_fp[0] + normal_vec_in[0] * t
    fy = lambda t: line_fp[1] + normal_vec_in[1] * t
    fz = lambda t: line_fp[2] + normal_vec_in[2] * t
    ###############################################################################
    # Create coordinates using 3D line equation-
    #
    coordinates = [[fx(t * delta), fy(t * delta), fz(t * delta)] for t in
                   range(int(depth / delta))]
    flat_coordinates = [entry for data in coordinates for entry in data]
    ###############################################################################
    # Create Field for coordinates of the path.
    #
    field_coord = dpf.fields_factory.create_3d_vector_field(len(coordinates))
    field_coord.data = flat_coordinates
    field_coord.scoping.ids = list(range(1, len(coordinates) + 1))
    ###############################################################################
    # Let's now map results on the path.
    mapping_operator = ops.mapping.on_coordinates(
        fields_container=stress_fc,
        coordinates=field_coord,
        create_support=True,
        mesh=mesh)
    fields_mapped = mapping_operator.outputs.fields_container()
    ###############################################################################
    # Here, we request the mapped field data and its mesh
    field_m = fields_mapped[0]
    mesh_m = field_m.meshed_region
    ###############################################################################
    # Create stress vs length chart.
    #
    x_initial = 0.0
    length = [x_initial + delta * index for index in range(len(field_m.data))]
    plt.plot(length, field_m.data, "r")
    plt.xlabel("Length (%s)" % mesh.unit)
    plt.ylabel("Stress (%s)" % field_m.unit)
    plt.show()
    ###############################################################################
    # To create a plot we need to add both the meshes
    # `mesh_m` - mapped mesh
    # `mesh` - original mesh
    pl = DpfPlotter()
    pl.add_field(field_m, mesh_m)
    pl.add_mesh(mesh, style="surface", show_edges=True,
                color="w", opacity=0.3)
    pl.show_figure(show_axes=True, cpos=[
        (62.687, 50.119, 67.247),
        (5.135, 6.458, -0.355),
        (-0.286, 0.897, -0.336)])
    

Answers