How to get stress failure criteria with temperature-dependent orthotropic strength limits?

Javier Vique

As we know, MAPDL allows us to get the stress failure criteria using FC command in POST1. However, if we want to stay at Mechanical GUI and have an interactive plot, how can it be done?

Javier Vique

mech-dpf allows us to plot the stress failure criteria. The script below can be used in a Python Result object to get the stress failure criteria, following a custom criteria:
-Only stresses on x and y (local 1 & 2) directions are considered.
-Stress strength on x is different for tensile and compression. For y direction, it is considered the same.
-User input is 3 dictionaries with temperatures and stress strength values. The values given below are dummy.
-It is assumed that the model only contains one material. Otherwise, the user should introduce a mesh_scoping on the specific material.
-A basic interpolation function is used, considering that the temperature is always within the values given in the dictionary previously introduced.
-Safety factor is got at centroid (elemental results)

def post_started(sender, analysis):# Do not edit this line
    define_dpf_workflow(analysis)

# Uncomment this function to enable retrieving results from the table/chart
def table_retrieve_result(value):# Do not edit this line
    import mech_dpf
    import Ans.DataProcessing as dpf
    wf = dpf.Workflow(this.WorkflowId)
    wf.Connect('contour_selector', value)
    this.Evaluate()

def define_dpf_workflow(analysis):

    def interpolate(dictionary, x):
        # Find the two keys closest to x
        keys = sorted(dictionary.keys())
        low_key = max(k for k in keys if k <= x)
        high_key = min(k for k in keys if k >= x)

        # If x is equal to one of the keys, return the corresponding value
        if low_key == high_key:
            return dictionary[low_key]

        # Perform linear interpolation
        low_value = dictionary[low_key]
        high_value = dictionary[high_key]

        # Interpolate
        interpolated_value = low_value + (high_value - low_value) * (x - low_key) / (high_key - low_key)
        return interpolated_value

    import mech_dpf
    import Ans.DataProcessing as dpf
    mech_dpf.setExtAPI(ExtAPI)

    #User Input
    data_0_t = {23.0: 100., 100.0: 150.}
    data_0_c = {23.0: 100., 100.0: 150.}
    data_90 = {23.0: 100., 100.0: 150.}

    analysis = ExtAPI.DataModel.Project.Model.Analyses[0]

    my_data_sources = dpf.DataSources(analysis.ResultFileName)
    model=dpf.Model(my_data_sources)

    # Set time scoping
    timeScop = dpf.Scoping()
    number_sets = model.TimeFreqSupport.NumberSets
    timeScop.Ids = range(1,number_sets+1)

    op = dpf.operators.result.structural_temperature()
    op.inputs.time_scoping.Connect(timeScop)
    op.inputs.data_sources.Connect(my_data_sources)
    op.inputs.bool_rotate_to_global.Connect(False)
    op.inputs.requested_location.Connect('Elemental')
    temp_elem = op.outputs.fields_container.GetData()

    op = dpf.operators.result.stress_X()
    op.inputs.time_scoping.Connect(timeScop)
    op.inputs.data_sources.Connect(my_data_sources)
    op.inputs.bool_rotate_to_global.Connect(False)
    op.inputs.requested_location.Connect('Elemental')
    sX_elem = op.outputs.fields_container.GetData()

    op = dpf.operators.result.stress_Y()
    op.inputs.time_scoping.Connect(timeScop)
    op.inputs.data_sources.Connect(my_data_sources)
    op.inputs.bool_rotate_to_global.Connect(False)
    op.inputs.requested_location.Connect('Elemental')
    sY_elem = op.outputs.fields_container.GetData()

    for nset in range(number_sets):
        SF_nset = dpf.FieldsFactory.CreateScalarField(numEntities=len(temp_elem[0].Scoping.Ids),location=dpf.locations.elemental)
        time = model.TimeFreqSupport.TimeFreqs.Data[nset]
        for Id in temp_elem[0].Scoping.Ids:
            temp_i = temp_elem[nset].GetEntityDataById(Id)[0]
            data_0_t_i = interpolate(data_0_t, temp_i)
            data_0_c_i = interpolate(data_0_c, temp_i)
            data_90_i = interpolate(data_90, temp_i)
            sXi = sX_elem[nset].GetEntityDataById(Id)[0]
            sYi = sY_elem[nset].GetEntityDataById(Id)[0]
            if sXi < 0.0:
                SF_x = abs(sXi)/data_0_c_i
            else:
                SF_x = sXi/data_0_t_i
            SF_y = abs(sYi)/data_90_i
            SF_i = max(SF_x,SF_y)
            SF_nset.Add(id=Id,data=[SF_i])
        if nset == 0:
            SF = dpf.FieldsContainerFactory.OverTimeFreqFieldsContainer([SF_nset])
        else:
            SF.AddFieldByTimeId(SF_nset,nset+1)

    SF.TimeFreqSupport = sX_elem.TimeFreqSupport

    output = dpf.operators.utility.forward()
    output.inputs.any.Connect(SF)

    dpf_workflow = dpf.Workflow()
    dpf_workflow.Add(output)
    dpf_workflow.SetOutputContour(output)
    dpf_workflow.Record('wf_id', False)
    this.WorkflowId = dpf_workflow.GetRecordedId()