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

Ayush Kumar
Member, Moderator, Employee Posts: 472
✭✭✭✭
How do I calculate stress along a path of desired length normal to a surface at a node?
Tagged:
2
Best Answer
-
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
andnode
. 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): = (x - p1[0]) / (p2[0] - p1[0]) p_int = p1[1] + (p2[1] - p1[1]) * 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)])
1
Answers
-
Ayush,
Very useful script, thanks for posting.
Could you show where this method is documented:
stress_fc = model.results.stress().eqv().outputs.fields_container()
If I look at the API I don't see how to do this.
0 -
Jim,
Not sure if this method is documented separately, I found it the generic example:
0