Plot contact pressure along deformed edge

Pernelle Marone-Hitz
Pernelle Marone-Hitz Member, Moderator, Employee Posts: 734
First Comment First Anniversary Ansys Employee Solution Developer Community of Practice Member

I have a model as follows:

I'd like to plot the pressure variation for nodes in contact with the upper plate, along with the deformation profile.

Is this possible?

Tagged:

Answers

  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 734
    First Comment First Anniversary Ansys Employee Solution Developer Community of Practice Member
    edited January 16

    Credits should go to @Pierre Thieffry who created the example.
    This can be achieved thanks to a Python Code object in Mechanical.

    The logic behind is as follows:

    • Create a named selection to get the nodes of the profile: - Note - when passed to the solver, the name of the named selection gets converted to capital letters.
    • Use DPF to retrieve:
    • *the contact pressure on the contact elements. This is done by selecting contact elements (type 172 in MAPDL) and extracting SMISC 5 to get the contact pressure.
      • the node coordinates in deformed position for the profile
    • The abscissa and pressure values are then exported to a file in the user_files directory of the project

    • To plot the graph in Mechanical, several options are possible. One way is to use matplotlib to do this, but matplotlib is a CPython library whereas Python Code operates in IronPython. The trick is to write a python script in the user_files directory that will be executed by a CPython engine thanks to subprocess.

      This script generates an image in the user_files directory


      and Mechanical can then grab this image and add it to the Mechanical tree.

      The code is a bit long so I will have to post this in several comments because of restrictions in the number of characters.

  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 734
    First Comment First Anniversary Ansys Employee Solution Developer Community of Practice Member

    Code up to helper functions:

    def after_post(this, solution):# Do not edit this line
        """
        Called after post processing.
        Keyword Arguments : 
            this -- the datamodel object instance of the python code object you are currently editing in the tree
            solution -- Solution
        """
    
        import math
        import os
        import mech_dpf
        import Ans.DataProcessing as dpf
    
        def get_user_files_directory():
            user_files = ''
            if ExtAPI.Context == 'Project':
                user_files = GetUserFilesDirectory()
            elif ExtAPI.Context == 'Mechanical':
                import wbjn
                user_files = wbjn.ExecuteCommand(ExtAPI, "returnValue(GetUserFilesDirectory())")
            return user_files
    
        def findClosestNode(ref_coord,node_coords):
            mindist=1e31
            minnode=-1
            for nd,coord in node_coords.items():
                dist=(ref_coord[0]-coord[0])**2+(ref_coord[1]-coord[1])**2+(ref_coord[2]-coord[2])**2
                if dist<mindist:
                    mindist=dist
                    minnode=nd
            return minnode,mindist
    
        def getCurvilinearAbscissa(node_coords):
            local_node_coords={nid:coord for nid,coord in node_coords.items()}
            curv_absc=[]
            distances={}
            # Find maximum node location in Y direction
            ymax=-1e31
            for nd,coord in local_node_coords.items():
                if ymax<coord[1]:
                    ymax=coord[1]
            # Filter nodes to those at Ymax - this will focus the results to the target metal sheet 
            # where maximum pressure should occur
            to_remove=[]
            for ni,coord in local_node_coords.items():
                if abs(coord[1]-ymax)>5e-3:
                    to_remove.append(ni)
            for ni in to_remove:                
                del local_node_coords[ni]
            # Now find chain of nodes from lowest X to largest X and
            # compute distance between nodes
            # find node with lowest X
            xmin=1e31
            for nd,coord in local_node_coords.items():
                if xmin>coord[0]:
                    xmin=coord[0]
                    ndmin=nd
            dist_idx=0
            ref_coord=local_node_coords[ndmin]
            distances[dist_idx]=[ndmin,0.]
            dist_idx+=1
            del local_node_coords[ndmin]
            # Now compute distance between nodes
            while len(local_node_coords.keys())>0:
                next_node,dist=findClosestNode(ref_coord,local_node_coords)
                distances[dist_idx]=[next_node,math.sqrt(dist)]
                ref_coord=node_coords[next_node]
                del local_node_coords[next_node]
                dist_idx+=1
            curv_absc=[[distances[0][0],xmin]]
            # Create curvilinear abscissa by adding distances between nodes
            for i in range(1,max(distances.keys())+1):
                curv_absc.append([distances[i][0],curv_absc[-1][1]+distances[i][1]])
            return curv_absc
    
  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 734
    First Comment First Anniversary Ansys Employee Solution Developer Community of Practice Member

    Main part of the code (up to exporting the values in text format):

        analysis=this.Parent.Parent
        mech_dpf.setExtAPI(ExtAPI)
    
        dataSource = dpf.DataSources(analysis.ResultFileName)
    
        # Retrieve scoping for all contact elements (type 172 in MAPDL)
        scoping_contact_elem=dpf.operators.scoping.on_property(requested_location='Elemental',
            property_name='mapdl_element_type',property_id=172,data_sources=dataSource)
    
        # Nodes from user named selection
        scoping_from_ns=dpf.operators.scoping.on_named_selection(requested_location='Nodal',
            named_selection_name='EDGEFORPLOT',data_sources=dataSource,int_inclusive=0)
    
        # Actual scopings
        elem_scoping=scoping_contact_elem.outputs.mesh_scoping.GetData()
        nodal_scoping=scoping_from_ns.outputs.mesh_scoping.GetData()
    
        # Operator to retrieve contact pressure (SMISC 5 for conta172 elements)
        cont_pres_elem= dpf.operators.result.mapdl.smisc(data_sources=dataSource,mesh_scoping=elem_scoping,item_index=5)
        cont_pres_nodal=dpf.operators.averaging.elemental_nodal_to_nodal_fc(cont_pres_elem)
    
        # Retrieve location of user_files directory for this project
        userFilesDir = get_user_files_directory()
    
        # Get pressure for all contact elements
        pressure_all_contacts=cont_pres_nodal.outputs.getfields_container()[0]
    
        # Rescope to nodes from selections
        pressure_selection_op=dpf.operators.scoping.rescope(fields=pressure_all_contacts,
            mesh_scoping=nodal_scoping)
    
        # Now this is the pressure field we want to plot over curvilinear abscissa
        pl=pressure_selection_op.outputs.getfields_as_field()
    
        # Get node ids and pressure values
        nodes=pl.ScopingIds
        values=pl.Data
    
        # Find node coordinates (in deformed position1)
        value_by_node={}
        node_coords={}
        rst_mesh=dpf.Model(dataSource).Mesh
        # Displacement field to update node coordinates with actual displacements
        u=dpf.operators.result.displacement(data_sources=dataSource,mesh_scoping=nodal_scoping)
        u_field=u.outputs.getfields_container()[0]
    
        for i in range(values.Count):
            val=values[i]
            value_by_node[nodes[i]]=val
            u_node=u_field.GetEntityDataById(nodes[i])
            nd=rst_mesh.NodeById(nodes[i])
            node_coords[nodes[i]]=[nd.X+u_node[0],nd.Y+u_node[1],nd.Z+u_node[2]]
    
        # Compute curvilinear abscissa
        curv_abs=getCurvilinearAbscissa(node_coords)
    
        # Reorder node coordinates on X values
        nids=[]
        xcoords=[]
        ycoords=[]
        for nid, coord in node_coords.items():
            nids.append(nid)
            xcoords.append(coord[0])
            ycoords.append(coord[1])
    
        sorted_coords = zip(xcoords,ycoords)
        sorted_coords.sort()
    
    
        # Export pressure vs abscissa
        myPath = os.path.join(userFilesDir,'pressure_on_path_{0}.txt'.format(dp))
        with open(myPath,'w') as f:
            f.write('Abscissa\tPressure}\n')
            for ca in curv_abs:
                f.write('{0:12.6e}\t{1:12.6e}\n'.format(ca[1],value_by_node[ca[0]]))
            # In same file, export node profile after deformation
            f.write('$\n')
            for row in sorted_coords:
                f.write('{0:12.6e}\t{1:12.6e}\n'.format(row[0],row[1]))
    
  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 734
    First Comment First Anniversary Ansys Employee Solution Developer Community of Practice Member

    Code to create the plot with Matplotlib and add the image to Mechanical:

        # Create plot using CPython code
        pngPath=os.path.join(userFilesDir,'pressure_on_path_{0}.png'.format(dp))
    
        c_python_script='''    
    import os
    import matplotlib.pyplot as plt
    
    with open(r"'''+myPath+'''","r") as f:
        x_for_plot=[]
        y_for_plot=[]
    
        x_profile=[]
        y_profile=[]
    
        line=f.readline()
        line=f.readline()
        while not('$' in line):
            sp_line=line.strip("\\n").split("\\t")
            x_for_plot.append(float(sp_line[0]))
            y_for_plot.append(float(sp_line[1]))    
            line=f.readline()
    
        line=f.readline()
        while line:
            sp_line=line.strip("\\n").split("\\t")
            x_profile.append(float(sp_line[0]))
            y_profile.append(float(sp_line[1]))    
            line=f.readline()
    
    fig, ax = plt.subplots(figsize = (10, 5))    
    
    # using the twinx() for creating another
    # axes object for secondary y-Axis
    ax2 = ax.twinx()
    ax.plot(x_for_plot,y_for_plot,'b')
    ax2.plot(x_profile,y_profile,'r-')
    
    # giving labels to the axises
    ax.set_xlabel('abscissa')
    ax.set_ylabel('pressure',color='b')
    ax2.set_ylabel('y location of node',color='r')
    
    plt.tight_layout()
    
    myPath = r"'''+pngPath+'''"
    fig.savefig(myPath)
    '''    
    
        scriptPath = os.path.join(userFilesDir,'generate_plot.py'.format(dp))
        with open(scriptPath,'w') as f:
            f.write(c_python_script)
    
        import sys,subprocess
        from System.Diagnostics import Process
    
        currentVersion=DataModel.Project.ProductVersion
        versionSplit=currentVersion.split(' ')    
        root_folder='AWP_ROOT'+versionSplit[0][2:]+versionSplit[-1][-1]    
        path1 = os.path.join( os.getenv(root_folder),r'commonfiles\CPython\3_7\winx64\Release\python', 'python.exe') # path to Mechanical's Python
        path2 = scriptPath
        # Argument for Python.exe call
        CMDARGS = '"{0}"'.format(path2)
        # Command for python exe including path to exe
        command = '"{0}"'.format(path1)
    
        ExtAPI.Log.WriteMessage(command)
        ExtAPI.Log.WriteMessage(CMDARGS)
        # Launch command with arguments
        proc = Process()
        proc.StartInfo.FileName = command
        proc.StartInfo.Arguments = CMDARGS
        proc.Start()
        proc.WaitForExit()
    
        # Add plot to Mechanical Tree
        with Transaction():
            images=ExtAPI.DataModel.GetObjectsByType(Ansys.ACT.Automation.Mechanical.Image)
            try:
                for im in images:
                    if im.Name=="Pressure on path":
                        im.Delete()
            except:
                pass
            image = solution.AddImage(pngPath)
            image.Name="Pressure on path"