How do I calculate linearized stress (Membrane, Bending and Peak) using DPF?

Ayush Kumar
Ayush Kumar Member, Moderator, Employee Posts: 450
100 Answers 250 Likes 100 Comments First Anniversary
✭✭✭✭
edited June 2023 in Structures

How do I calculate linearized stress (Membrane, Bending and Peak) using DPF?

Tagged:

Comments

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 450
    100 Answers 250 Likes 100 Comments First Anniversary
    ✭✭✭✭
    edited May 2023
    import mech_dpf
    import Ans.DataProcessing as dpf
    
    
    path = ExtAPI.DataModel.GetObjectsByName("Path")[0]
    x_0 = path.StartXCoordinate.Value
    y_0 = path.StartYCoordinate.Value
    z_0 = path.StartZCoordinate.Value
    
    
    x_1 = path.EndXCoordinate.Value
    y_1 = path.EndYCoordinate.Value
    z_1 = path.EndZCoordinate.Value
    
    
    path_length = ((x_1 - x_0) ** 2 + (y_1 - y_0) ** 2 + (z_1 - z_0) ** 2) ** 0.5
    
    
    n_points = path.DiscretizationPoints + 2
    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.FieldsFactory.Create3DVectorField(n_points)
    field_coord.Data = flat_coordinates
    field_coord.Scoping.Ids = list(range(1, n_points + 1))
    
    
    model = dpf.Model(Model.Analyses[0].ResultFileName)
    
    
    ts = dpf.Scoping()
    ts.Ids = [2]
    
    
    # SX
    s = model.CreateOperator("SX")
    s.inputs.requested_location.Connect("Nodal")
    s.inputs.time_scoping.Connect(ts)
    s2_f = s.outputs.fields_container.GetData()
    
    
    mapping_operator = dpf.operators.mapping.on_coordinates(
        fields_container=s2_f,
        coordinates=field_coord,
        create_support=True,
        mesh=model.Mesh
    )
    fields_mapped = mapping_operator.outputs.fields_container.GetData()
    
    
    ls = list(fields_mapped[0].Data)
    # Membrane Stress
    membrane_stress = (ls[0] / 2 + ls[-1] / 2 + sum(ls[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.FieldsFactory.CreateScalarField(numEntities=n_points, location=dpf.locations.nodal)
    path_range_field.Data = path_range
    path_range_field.ScopingIds = 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.GetData().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
    
    
    b1 = stress_scaled_integral * (-6.0 / path_length) / 48.0
    b2 = (-1.0) * b1
    <?WORD FLAGGED?>_bending = (b2 - b1) / (path_length)
    fb = lambda t: b1 + <?WORD FLAGGED?>_bending * t
    bending_stress = [fb(t * delta) for t in range(n_points)]
    
    
    bend_f = dpf.FieldsFactory.CreateScalarField(numEntities=49, location=dpf.locations.nodal)
    bend_f.Data = bending_stress
    bend_f.Scoping.Ids = list(range(1, n_points + 1))
    
    
    # Add Membrane and Bending Stress
    bend_mem_f = dpf.operators.math.add_constant(field=bend_f, ponderation=membrane_stress)
    bend_mem_f_neg = dpf.operators.math.scale(field=bend_mem_f, ponderation=-1.0)
    
    
    # Peak stress
    ps = dpf.operators.math.add(fieldA=fields_mapped[0], fieldB=bend_mem_f_neg)
    
    
    
    
    
  • dafedin
    dafedin Member Posts: 21
    10 Comments First Anniversary Name Dropper
    **

    Hello Ayush,

    Thanks for this comprehensive example. I wonder if there is a way to solve the reversed task: find elements of the model whose results were used for mapping to path points. Maybe there is another option rather than screening over the whole mesh and finding elements whose centroids are situated in the path direction.

    Regards

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 450
    100 Answers 250 Likes 100 Comments First Anniversary
    ✭✭✭✭

    Hi @dafedin,

    That information is not provided as an output of the operator at the moment.

  • Jim Kosloski
    Jim Kosloski Member, Employee Posts: 24
    10 Comments Name Dropper First Anniversary Ansys Employee
    ✭✭✭✭

    Is there a more efficient way than looping where we could get these results (bending and membrane) for ALL timepoints? Most of our DPF examples focus on a single time point, but one of the big benefits of DPF is the ability to process multiple timepoints simultaneously. All of these calculations use fields_mapped[0].Data how can we use all the data in the fields_mapped container?

    Also how can I use this in a Python Result and plot the path results on the path. If I send bend_f to the dpf_workflow.Add I get an error.

  • Jim Kosloski
    Jim Kosloski Member, Employee Posts: 24
    10 Comments Name Dropper First Anniversary Ansys Employee
    ✭✭✭✭

    Is there a more efficient way than looping where we could get these results (bending and membrane) for ALL timepoints? Most of our DPF examples focus on a single time point, but one of the big benefits of DPF is the ability to process multiple timepoints simultaneously. All of these calculations use fields_mapped[0].Data how can we use all the data in the fields_mapped container?

    Also how can I use this in a Python Result and plot the path results on the path. If I send bend_f to the dpf_workflow.Add I get an error.

  • Jim Kosloski
    Jim Kosloski Member, Employee Posts: 24
    10 Comments Name Dropper First Anniversary Ansys Employee
    ✭✭✭✭

    @Ayush Kumar a few more questions on this. 1. What does the math.scale_by_field operator do? I do not see it documented anywhere.
    2. Do the field coordinates input into the mapping operator need to be in SI?

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 450
    100 Answers 250 Likes 100 Comments First Anniversary
    ✭✭✭✭
    edited August 31

    @Jim Kosloski
    1. As most of the calculation is done on the Client side, I don't see a DPF way to calculate for all time-steps at once (in Mechanical DPF). For PyDPF, NumPy can be used for vectorized operations.
    2. You can try forwarding the field to an operator (dpf.operators.utility.forward()) and then pass that operator to workflow. Also make sure you have the new mesh defined using SetOutputMesh(...).
    3. scale_by_field: https://dpf.docs.pyansys.com/version/stable/api/ansys.dpf.core.operators.math.scale_by_field.html
    4. Coordinate field input in the mapping operator is not unit dependent.