Stress sampling example

This example shows how to sample the stresses in the rotor bridges.

This script should be run from the scripting tab after the stress calculation has been run in Motor-CAD.

Point positions: [[53.38678975481784, 44.47098476450547], [53.30031932182026, 44.57458655379761], [53.21364784070531, 44.6780202083239], [53.12677563839672, 44.78128533793393], [53.03970304257532, 44.88438155311296], [52.952430381677786, 44.987308464983414], [52.86495798489549, 45.09006568530628], [52.77728618217318, 45.192652826482664], [52.689415304207785, 45.29506950155518], [52.6013456824471, 45.39731532420947], [52.513077649088615, 45.49938990877559], [52.42461153707821, 45.60129287022953], [52.335947680108916, 45.70302382419464], [52.24708641261966, 45.80458238694305], [52.158028069794, 45.905968175397206]]
Stress at points: [23.63, 27.58, 27.58, 23.54, 31.11, 25.46, 30.78, 28.17, 28.17, 34.03, 34.99, 34.99, 27.19, 27.19, 23.12]
Mean stress: 28.502000000000002
Point positions: [[57.956238644795526, 39.068982929524914], [57.89539584745769, 39.15908789369561], [57.834412963433635, 39.24909810651729], [57.773290140280686, 39.33901335019659], [57.71202752589479, 39.428833407169954], [57.65062526851012, 39.51855806010415], [57.58908351669874, 39.60818709189676], [57.52740241937025, 39.69772028567677], [57.46558212577142, 39.78715742480505], [57.40362278548582, 39.876498292874864], [57.34152454843347, 39.965742673712455], [57.279287564870465, 40.054890351377516], [57.21691198538861, 40.143941110163716], [57.1543979609151, 40.23289473459925], [57.091745642712056, 40.32175100944735]]
Stress at points: [23.28, 23.28, 20.47, 27.1, 27.1, 22.85, 28.97, 28.97, 23.31, 23.31, 32.2, 32.2, 17.53, 17.67, 17.67]
Mean stress: 24.394000000000002

from math import asin, cos, radians, sin
import os
import sys

# Standard imports
import ansys.motorcad.core as pymotorcad

# Connect to Motor-CAD
mcApp = pymotorcad.MotorCAD()

# Users should run this script from the scripting tab after the stress calculation
# Trigger this automatically for the automated documentation build
if "PYMOTORCAD_DOCS_BUILD" in os.environ:
    mcApp.set_variable("MessageDisplayState", 2)
    mcApp.load_template("e10")
    mcApp.do_mechanical_calculation()

# Sample points is hardcoded to 15 in Motor-CAD for stress averaging
sample_points = 15

# Check the rotor type
rotor_type = mcApp.get_variable("BPMRotor")
# U shape is 13, V web is 11
if rotor_type == 11:
    layers = mcApp.get_variable("VMagnet_Layers")
elif rotor_type == 13:
    layers = mcApp.get_variable("Magnet_Layers")
else:
    sys.exit("Stress sampling only available for V web and U templates")

# Get variables independent of the rotor type
average_stress_location_bridge = mcApp.get_variable("AvStressRadialLocation_Bridge")
rotor_diameter = mcApp.get_variable("RotorDiameter")
poles = mcApp.get_variable("Pole_Number")
pole_pairs = poles / 2

for layer in range(layers):
    if rotor_type == 11:
        # V Web template
        bridge_thickness = mcApp.get_array_variable("BridgeThickness_Array", layer)
        web_thickness = mcApp.get_array_variable("WebThickness_Array", layer)
        pole_arc = radians(mcApp.get_array_variable("PoleArc_Array", layer) / pole_pairs)
        theta_4 = asin(web_thickness / (2 * (rotor_diameter / 2 - bridge_thickness)))
        theta_0 = radians(180 / poles) + pole_arc / 2
        theta_1 = radians(360 / poles) - theta_4
        theta_bridge_span = theta_1 - theta_0
        # Arc covers half the bridge
        delta_theta = theta_bridge_span / 2 / sample_points
        theta = theta_0 + theta_bridge_span / 2
    elif rotor_type == 13:
        # U template
        bridge_thickness = mcApp.get_array_variable("UShape_BridgeThickness_Array", layer)
        web_thickness = mcApp.get_array_variable("UShape_WebThickness_Array", layer)
        outer_thickness = mcApp.get_array_variable("UShape_Thickness_Outer_Array", layer)
        theta_offset = radians(mcApp.get_array_variable("UShape_OuterAngleOffset_Array", layer))
        inner_rad = rotor_diameter / 2 - bridge_thickness

        # The start angle of the FEA model
        theta_0 = radians(360 / poles)

        # The angle to the end of the web
        theta_1 = asin(web_thickness / (2 * inner_rad))

        # We now need to solve for phi (the arc angle from the centre) with
        # cos(theta_offset - theta_1 - phi/2) * sin(phi/2) = outer_thickness / (2 * inner_rad)
        # This is non-trivial, so do this numerically
        test_phi = 0
        found_phi = False
        phi_step = radians(0.01)
        iteration = 0
        last_err = cos(theta_offset - theta_1 - test_phi / 2) * sin(
            test_phi / 2
        ) - outer_thickness / (2 * inner_rad)
        while found_phi == False and iteration < 36000:
            test_phi = test_phi + phi_step
            iteration = iteration + 1
            err = cos(theta_offset - theta_1 - test_phi / 2) * sin(
                test_phi / 2
            ) - outer_thickness / (2 * inner_rad)
            # Check if error has changed sign, if so we are close to the correct solution
            if err * last_err < 0:
                found_phi = True
            else:
                last_err = err
        theta_bridge_span = test_phi

        # Arc covers half the bridge
        delta_theta = theta_bridge_span / 2 / sample_points
        theta = theta_0 - (theta_1 + theta_bridge_span / 2)

    # Common logic for both rotor types
    r0 = rotor_diameter / 2 - bridge_thickness * (1 - average_stress_location_bridge)
    stresses = []
    points = []
    for point in range(sample_points):
        this_theta = theta + point * delta_theta
        x = r0 * cos(this_theta)
        y = r0 * sin(this_theta)
        points.append([x, y])
        stresses.append(mcApp.get_point_value("SVM", x, y)[0])

    # Display results for this layer
    print("Point positions: " + str(points))
    print("Stress at points: " + str(stresses))
    print("Mean stress: " + str(sum(stresses) / len(stresses)))