How to Use Python Results Object for the Mohr-Coulomb Failure Criterion?

Faysal
Faysal Member Posts: 1 **
in General Language Questions

Hello everyone,

I have been attempting to check for shear failure in my static structural model by applying the Mohr-Coulomb failure criterion using the Python Result object. To achieve this, I followed these steps:

Extracted the minimum and maximum principal stress results.
Adjusted the sign of these stresses, as I am conducting a geomechanical analysis where compressive stress is considered positive and tensile stress negative.
Evaluated shear failure for each element using the Mohr-Coulomb failure criterion.

I wrote a code based on a discussion I found here:
https://discuss.ansys.com/discussion/3068/how-to-get-stress-failure-criteria-with-temperature-dependent-orthotropic-strength-limits

Unfortunately, my code is not producing the expected results. As I am new to this, I am struggling to identify the issue. I would greatly appreciate any guidance or suggestions.

Thank you in advance for your help.

Best regards,
Faysal

Here's my code:

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

def define_dpf_workflow(analysis):
import mech_dpf
import Ans.DataProcessing as dpf
import math

mech_dpf.setExtAPI(ExtAPI)

# User Input: Mohr-Coulomb Material Properties
cohesion = 8.2  # MPa
friction_angle = 48  # degrees
phi_rad = math.radians(friction_angle)  # Convert to radians

# Load Model
analysis = ExtAPI.DataModel.Project.Model.Analyses[0]
my_data_sources = dpf.DataSources(analysis.ResultFileName)
model = dpf.Model(my_data_sources)

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

# Get Principal Stresses
op1 = dpf.operators.result.stress_principal_1()
op1.inputs.time_scoping.Connect(timeScop)
op1.inputs.data_sources.Connect(my_data_sources)
op1.inputs.bool_rotate_to_global.Connect(False)
op1.inputs.requested_location.Connect("Elemental")
s1_elem = op1.outputs.fields_container.GetData()

op3 = dpf.operators.result.stress_principal_3()
op3.inputs.time_scoping.Connect(timeScop)
op3.inputs.data_sources.Connect(my_data_sources)
op3.inputs.bool_rotate_to_global.Connect(False)
op3.inputs.requested_location.Connect("Elemental")
s3_elem = op3.outputs.fields_container.GetData()

# Compute Failure Check
for nset in range(number_sets):
    FE_nset = dpf.FieldsFactory.CreateScalarField(numEntities=len(s1_elem[0].Scoping.Ids), location=dpf.locations.elemental)
    time = model.TimeFreqSupport.TimeFreqs.Data[nset]

    for Id in s1_elem[nset].Scoping.Ids:
        # Extract and change sign to match Geomechanics convention
        sigma_1_ansys = s1_elem[nset].GetEntityDataById(Id)[0]  # Max principal stress (ANSYS)
        sigma_3_ansys = s3_elem[nset].GetEntityDataById(Id)[0]  # Min principal stress (ANSYS)

        # Change sign and swap to match Geomechanics
        sigma_1 = -sigma_3_ansys  # After sign change, S3 becomes S1
        sigma_3 = -sigma_1_ansys  # After sign change, S1 becomes S3

        # Compute Mohr-Coulomb failure line equation for sigma_1 vs sigma_3
        failure_line_sigma_1 = 2 * cohesion * math.cos(phi_rad) / (1 - math.sin(phi_rad)) + sigma_3 * (1 + math.sin(phi_rad)) / (1 - math.sin(phi_rad))

        # Check for failure: If (sigma_1, sigma_3) is above or on the failure line
        failure = float(1 if sigma_1 >= failure_line_sigma_1 else 0)  # Ensure float type

        # Add failure flag to field
        FE_nset.Add(id=Id, data=[failure])

    # Ensure correct initialization of FE container
    if nset == 0:
        FE = dpf.FieldsContainer()
    FE.AddFieldByTimeId(FE_nset, nset + 1)

FE.TimeFreqSupport = s1_elem.TimeFreqSupport

# Create and connect workflow
output = dpf.operators.utility.forward()
output.inputs.any.Connect(FE)

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

Categories