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

Patrick Herberich

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

Patrick Herberich

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

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

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

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

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

Pernelle Marone-Hitz

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

Jimmy He

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()
    ######################################################################