Is it possible to consider temperature dependent strength data for safety factor calculation?

Patrick Herberich
Patrick Herberich Member, Employee Posts: 2
First Comment First Anniversary Ansys Employee Photogenic
✭✭✭

Is it possible to consider temperature dependent strength data for safety factor calculation within the stress tool?

Comments

  • Patrick Herberich
    Patrick Herberich Member, Employee Posts: 2
    First Comment First Anniversary Ansys Employee Photogenic
    ✭✭✭

    It is possible via a Python Result object inside Mechanical. The following script calculates the safety factor considering temperature dependent limit stresses:

    def post_started(sender, analysis):# Do not edit this line
        define_dpf_workflow(analysis)
    
    
    def LinearInterpolation(x,x1,y1,x2,y2):
        # Simple linear interpolation
        # Returns -1 if x is outside of [x1,x2]
        if x<x1 or x>x2:
            return -1 # out of range
        if y1==y2:
            return y1
        else:
            a=(y2-y1)/(x2-x1)
            b=y1-a*x1
            return a*x+b
    
    def GetLimitStress(temperature,stress_data):
        # Find limit stress value for a given temperature
        # This is done by simple linear interpolation between two 
        # points of a limit stress/temperature curve
    
        # Find closest temperature in stress_data
        num_points=len(stress_data)/2
        for i in range(num_points):
            temp=stress_data[2*i]
            if temp>temperature:
                index=i
                break
        # Get temperature/limit stress points between which to interpolat          
        if index!=num_points:
            x2=stress_data[index*2]
            x1=stress_data[(index-1)*2]
            y2=stress_data[index*2+1]
            y1=stress_data[(index-1)*2+1]
        else:
            return -2 # temperature value is out of bounds
    
        return LinearInterpolation(temperature,x1,y1,x2,y2) # return limit stress
    
    def GetLimitStressField(stress,temperatures,stress_data):
        # From a field of stresses and temperatures, find limit stress
        # at each node. Stress_data defines the limit stress/temperature curve
        import mech_dpf
        import Ans.DataProcessing as dpf
    
        stress_nodes=stress.outputs.getfields_container()[0].ScopingIds
        limit_field=dpf.FieldsFactory.CreateScalarField(1,'Nodal')
    
        local_temp=temperatures.outputs.getfields_container()[0]
    
        for nid in stress_nodes:
            temp=local_temp.GetEntityDataById(nid)[0]
            limit_field.Add(nid,[float(GetLimitStress(temp,stress_data))])
    
        return limit_field
    
    def define_dpf_workflow(analysis):
        import mech_dpf
        import Ans.DataProcessing as dpf
    
        # Define limit stress data as list in form of [temp1,limit_stress1,temp2,limit_stress2,...])
        limit_stress_data=[-100.,200.,0.,220.,50.,210.,100,200.,1000.,150.]
    
        my_data_sources = dpf.DataSources(analysis.ResultFileName)
        my_time_scoping = dpf.Scoping()
        my_time_scoping.Ids = [1] # the first set
    
        # Retrieve structural temperature at nodes
        temperature=dpf.operators.result.structural_temperature(requested_location='Nodal')
        temperature.inputs.data_sources.Connect(my_data_sources)
        temperature.inputs.time_scoping.Connect(my_time_scoping)
    
        # Retrieve Mises stress at nodes
        s_eqv_op = dpf.operators.result.stress_von_mises()
        s_eqv_op.inputs.requested_location.Connect('Nodal')
        s_eqv_op.inputs.data_sources.Connect(my_data_sources)
        s_eqv_op.inputs.time_scoping.Connect(my_time_scoping)
    
        # Retrieve limit_stress for each node based on Mises and temperature
        limit_stress=GetLimitStressField(s_eqv_op,temperature,limit_stress_data)
    
        # Divide limit stress by Mises stress
        limit_factor=dpf.operators.math.component_wise_divide() 
        limit_factor.inputs.fieldA.Connect(limit_stress)
        limit_factor.inputs.fieldB.Connect(s_eqv_op)
    
    
        dpf_workflow = dpf.Workflow()
        dpf_workflow.Add(limit_factor)
        dpf_workflow.SetOutputContour(limit_factor)
        dpf_workflow.Record('wf_id', True)
        this.WorkflowId = dpf_workflow.GetRecordedId()
    
  • Mike.Thompson
    Mike.Thompson Member, Employee Posts: 339
    25 Answers 100 Comments 25 Likes First Anniversary
    ✭✭✭✭

    Only for clarification and full context, I wanted to point out that mechanical has this general capability via a stress tool post processing object that will plot factors of safety that are temperature dependent based on material property.

    This forum is for scripted solutions, so obviously the solution is on point, and it is a good starting point for implementing more complex or specific processing routines, then what is offered in a native solution.

  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 871
    100 Answers 500 Comments 250 Likes First Anniversary
    ✭✭✭✭

    My understanding is that Mechanical's stress tool will indeed take into account temperature-dependent stress (as we can define E = f(T)),

    but natively it is not possible to define that the Tensile Yield Stress is temperature-dependent:

    which is what the code from @Patrick Herberich does.

  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 871
    100 Answers 500 Comments 250 Likes First Anniversary
    ✭✭✭✭

    Piggybacking on this solution. I've slightly modified the code so that instead of having the Yield Strenght vs Temp values hard-coded, they are read from a comma delimited csv file placed in the user_files directory of the project.


    Code is:

    def post_started(sender, analysis):# Do not edit this line
        define_dpf_workflow(analysis)
    
    def 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 LinearInterpolation(x,x1,y1,x2,y2):
        # Simple linear interpolation
        # Returns -1 if x is outside of [x1,x2]
        if x<x1 or x>x2:
            return -1 # out of range
        if y1==y2:
            return y1
        else:
            a=(y2-y1)/(x2-x1)
            b=y1-a*x1
            return a*x+b
    
    def GetLimitStress(temperature,stress_data):
        # Find limit stress value for a given temperature
        # This is done by simple linear interpolation between two 
        # points of a limit stress/temperature curve
    
        # Find closest temperature in stress_data
        num_points=len(stress_data)/2
        for i in range(num_points):
            temp=stress_data[2*i]
            if temp>temperature:
                index=i
                break
        # Get temperature/limit stress points between which to interpolat          
        if index!=num_points:
            x2=stress_data[index*2]
            x1=stress_data[(index-1)*2]
            y2=stress_data[index*2+1]
            y1=stress_data[(index-1)*2+1]
        else:
            return -2 # temperature value is out of bounds
    
        return LinearInterpolation(temperature,x1,y1,x2,y2) # return limit stress
    
    def GetLimitStressField(stress,temperatures,stress_data):
        # From a field of stresses and temperatures, find limit stress
        # at each node. Stress_data defines the limit stress/temperature curve
        import mech_dpf
        import Ans.DataProcessing as dpf
    
        stress_nodes=stress.outputs.getfields_container()[0].ScopingIds
        limit_field=dpf.FieldsFactory.CreateScalarField(1,'Nodal')
    
        local_temp=temperatures.outputs.getfields_container()[0]
    
        for nid in stress_nodes:
            temp=local_temp.GetEntityDataById(nid)[0]
            limit_field.Add(nid,[float(GetLimitStress(temp,stress_data))])
    
        return limit_field
    
    def get_stress_data():
        import csv
        import os
        limit_stress_data = []
        user_dir = user_files_directory()
        if os.path.exists(os.path.join(user_dir, 'TemperatureDependentYieldStrength.csv')):
            with open(os.path.join(user_dir, 'TemperatureDependentYieldStrength.csv')) as csv_file:
                csv_reader = csv.reader(csv_file, delimiter=',')
                line_count = 0
                for row in csv_reader:
                    if line_count == 0:
                        line_count += 1
                    else:
                        limit_stress_data.append(float(row[0]))
                        limit_stress_data.append(float(row[1]))
                        line_count += 1
        else:
            msg = Ansys.Mechanical.Application.Message('Temperature Dependent Yield Strenght file does not exist', MessageSeverityType.Error)
            ExtAPI.Application.Messages.Add(msg)
        return limit_stress_data
    
    def define_dpf_workflow(analysis):
        import mech_dpf
        import Ans.DataProcessing as dpf
    
        # Define limit stress data as list in form of [temp1,limit_stress1,temp2,limit_stress2,...])
        limit_stress_data = get_stress_data()
    
        my_data_sources = dpf.DataSources(analysis.ResultFileName)
        my_time_scoping = dpf.Scoping()
        my_time_scoping.Ids = [1] # the first set
    
        # Retrieve structural temperature at nodes
        temperature=dpf.operators.result.structural_temperature(requested_location='Nodal')
        temperature.inputs.data_sources.Connect(my_data_sources)
        temperature.inputs.time_scoping.Connect(my_time_scoping)
    
        # Retrieve Mises stress at nodes
        s_eqv_op = dpf.operators.result.stress_von_mises()
        s_eqv_op.inputs.requested_location.Connect('Nodal')
        s_eqv_op.inputs.data_sources.Connect(my_data_sources)
        s_eqv_op.inputs.time_scoping.Connect(my_time_scoping)
    
        # Retrieve limit_stress for each node based on Mises and temperature
        limit_stress=GetLimitStressField(s_eqv_op,temperature,limit_stress_data)
    
        # Divide limit stress by Mises stress
        limit_factor=dpf.operators.math.component_wise_divide() 
        limit_factor.inputs.fieldA.Connect(limit_stress)
        limit_factor.inputs.fieldB.Connect(s_eqv_op)
    
    
        dpf_workflow = dpf.Workflow()
        dpf_workflow.Add(limit_factor)
        dpf_workflow.SetOutputContour(limit_factor)
        dpf_workflow.Record('wf_id', True)
        this.WorkflowId = dpf_workflow.GetRecordedId()
    
  • Abath
    Abath Member Posts: 1
    First Comment
    **

    I am planning to run this code for multiple materials. How can I proceed?

  • Pernelle Marone-Hitz
    Pernelle Marone-Hitz Member, Moderator, Employee Posts: 871
    100 Answers 500 Comments 250 Likes First Anniversary
    ✭✭✭✭

    You'll have to identify which body uses which material and restrict the scoping to the node Ids of each body.

  • Jimmy He
    Jimmy He Member, Employee Posts: 19
    10 Comments First Anniversary Name Dropper Photogenic
    ✭✭✭

    Hi @Abath , FYI this is an example I created for multiple bodies and multiple materials:

    def post_started(sender, analysis):# Do not edit this line
        define_dpf_workflow(analysis)
    
    
    def define_dpf_workflow(analysis):
        ######################################################################
        # Usually no need to change this part
        # Import libraries
        import mech_dpf
        import Ans.DataProcessing as dpf
        import math
        mech_dpf.setExtAPI(ExtAPI)
        my_data_sources = dpf.DataSources(analysis.ResultFileName)
    
        # Read mesh in results file
        mesh_op = dpf.operators.mesh.mesh_provider() # operator instanciation
        mesh_op.inputs.data_sources.Connect(my_data_sources)
        mesh = mesh_op.outputs.mesh.GetData()
        ######################################################################
    
        # import materials module
        import materials 
    
        # Credit: https://stackoverflow.com/questions/7343697/how-to-implement-linear-interpolation
        from bisect import bisect_right
        class Interpolate:
            def __init__(self, x_list, y_list):
                if any(y - x <= 0 for x, y in zip(x_list, x_list[1:])):
                    raise ValueError("x_list must be in strictly ascending order!")
                self.x_list = x_list
                self.y_list = y_list
                intervals = zip(x_list, x_list[1:], y_list, y_list[1:])
                self.slopes = [(y2 - y1) / (x2 - x1) for x1, x2, y1, y2 in intervals]
    
            def __call__(self, x):
                if self.x_list[0] >= x:
                    return self.y_list[0]
                elif x >= self.x_list[-1]:
                    return self.y_list[-1]
                i = bisect_right(self.x_list, x) - 1
                return self.y_list[i] + self.slopes[i] * (x - self.x_list[i])
    
        # Prepare all interpolation objects
        inpt_objs = {}
        for mat in ExtAPI.DataModel.Project.Model.Materials.Children:
            # get engineering data material properties for this material
            matED = mat.GetEngineeringDataMaterial()
            props = materials.GetMaterialPropertyByName(matED,"Isotropic Hardening")
            Temp = props['Temperature'][1:]
            YS = props['Yield Strength'][1:] # Unit at index 0
    
            inpt_objs[ mat.Name ] = Interpolate( Temp , YS )
    
    
        # Getting mesh info
        meshData = DataModel.AnalysisList[0].MeshData # reference meshData
        #print( len(meshData.Elements) )
        #print( len(meshData.Nodes) )
    
        mats , nodes = [] , []
        for p in DataModel.GetObjectsByType(Ansys.ACT.Automation.Mechanical.Body):
            bid = p.GetGeoBody().Id
            mat = p.Material
            bodyMesh = meshData.MeshRegionById(bid) # get mesh data
    
            nodes += list( bodyMesh.NodeIds )
            mats += [mat] * bodyMesh.NodeCount
    
        node2mat = dict( zip( nodes , mats ))
    
    
        ############################################################################
        # Get fields from DPF
        # vm stress nodal
        vm = dpf.operators.result.stress_von_mises()
        vm.inputs.requested_location.Connect('Nodal')
        vm.inputs.data_sources.Connect(my_data_sources)
        nIds = vm.outputs.fields_container.GetData()[0].ScopingIds
        vm = vm.outputs.fields_container.GetData()[0].Data
    
        T = dpf.operators.result.structural_temperature()
        T.inputs.requested_location.Connect('Nodal')
        T.inputs.data_sources.Connect(my_data_sources)
        temp = T.outputs.fields_container.GetData()[0].Data
    
    
        # Create a blank field and assign values
        my_field = dpf.FieldsFactory.CreateScalarField(numEntities=0, location='Nodal')
        my_field.MeshedRegionSupport = mesh
    
        # Calculate field value
        for i , s , t in zip( nIds , vm , temp ):
            my_mat = node2mat[i]
            my_ys = inpt_objs[ my_mat ]( t ) / 1e6 # Convert to MPa, hard coded for demo
    
            fos = my_ys / ( s + 1e-6 )
    
            # Insert field value
            my_field.Add(i,[ fos ])
        my_field.ScopingIds = nIds
    
    
        ######################################################################
        # Usually no need to change this part
        # Forward field to an operator to plot it
        forward = dpf.operators.utility.forward_field()
        forward.inputs.field.Connect(my_field)
    
        dpf_workflow = dpf.Workflow()
        dpf_workflow.Add(forward)
        dpf_workflow.SetOutputContour(forward,dpf.enums.GFXContourType.FENodalScoping)
        dpf_workflow.Record('wf_id', False)
        this.WorkflowId = dpf_workflow.GetRecordedId()
        ######################################################################