Note
Go to the end to download the full example code.
Stress sampling example#
This example shows how to sample the stresses in the rotor bridges.
In addition to calculating the average rotor bridge stress, it also allows post-processing to consider non-linear plastic effects. Optionally, the results from the non-linear calculation are used to replace the standard result in the outputs tab.
This script should be run from the scripting tab after the stress calculation has been run in Motor-CAD.
Standard imports#
import math
import os
import sys
import numpy as np
import ansys.motorcad.core as pymotorcad
Utility functions to check non-linear material data#
def check_youngs_modulus(non_linear_strain, non_linear_stress, youngs_modulus):
"""Check the initial slope of the non-linear stress strain curve matches the Young's modulus.
Parameters
----------
non_linear_strain : list or ndarray
Strain values
non_linear_stress : list or ndarray
Corresponding stress values in MPa
youngs_modulus : float
Young's modulus in MPa
Returns
-------
None
Raises
------
ValueError
If the initial slope is not consistent with the Young's modulus.
"""
initial_youngs_modulus = (non_linear_stress[1] - non_linear_stress[0]) / (
non_linear_strain[1] - non_linear_strain[0]
)
# Error if notably different
if not math.isclose(initial_youngs_modulus, youngs_modulus, rel_tol=0.01, abs_tol=0.1):
raise ValueError(
"Youngs Modulus and initial slope of non-linear data are different, "
"please check the inputs. Initial slope is "
+ str(initial_youngs_modulus)
+ " MPa "
+ "Young's modulus is "
+ str(youngs_modulus)
+ " MPa"
)
def find_divergence_point(non_linear_strain, non_linear_stress, youngs_modulus):
"""Find last point where the non-linear stress strain curve is on the linear region.
Parameters
----------
non_linear_strain : list or ndarray
Strain values
non_linear_stress : list or ndarray
Corresponding stress values in MPa
youngs_modulus : float
Young's modulus in MPa
Returns
-------
float
The last stress point in the linear range
"""
for i in range(1, len(non_linear_stress)):
if not math.isclose(
non_linear_stress[i] / non_linear_strain[i], youngs_modulus, rel_tol=0.0001
):
return non_linear_stress[i - 1]
# Return maximum stress in list if no divergence found
return non_linear_stress[-1]
Functions to estimate plastic results#
def apply_neuber_correction(elastic_stress, youngs_modulus, non_linear_strain, non_linear_stress):
"""Update the Neuber correction estimates from the Von Mises stress.
Parameters
----------
elastic_stress : float
The linear stress input
youngs_modulus : float
The Young's modulus used in the linear stress calculation.
non_linear_strain, non_linear_stress : np.ndarray
Corresponding arrays of strain (ratio) and stress (MPa) for the non-linear correction.
"""
elastic_strain = elastic_stress / youngs_modulus
elastic_stress_strain_product = elastic_stress * elastic_strain
# Check input data is sensible
check_youngs_modulus(non_linear_strain, non_linear_stress, youngs_modulus)
# If on the elastic portion, just return inputs, so we don't get interpolation errors
if elastic_stress < find_divergence_point(non_linear_strain, non_linear_stress, youngs_modulus):
return {
"nonlinear_strain": elastic_strain,
"nonlinear_stress": elastic_stress,
"plastic_strain": 0,
}
# Find a matching stress-strain product in the non-linear response
# (This is the Neuber correction)
non_linear_stress_strain_product = non_linear_stress * non_linear_strain
# Error if out of range:
if elastic_stress_strain_product > np.max(non_linear_stress_strain_product):
raise ValueError(
"Input too large (elastic stress strain product > maximum stress strain product in "
"non-linear data). "
"Elastic stress is "
+ str(elastic_stress)
+ ", elastic stress strain product is "
+ str(elastic_stress_strain_product)
+ ", maximum plastic stress strain product is "
+ str(np.max(non_linear_stress_strain_product))
)
# Lookup to find non-linear strain at matching stress strain product:
equivalent_non_linear_strain = np.interp(
elastic_stress_strain_product, non_linear_stress_strain_product, non_linear_strain
)
# Lookup to find non-linear stress at this strain:
equivalent_non_linear_stress = np.interp(
equivalent_non_linear_strain, non_linear_strain, non_linear_stress
)
# Find plastic strain
plastic_strain = equivalent_non_linear_strain - elastic_strain
return {
"nonlinear_strain": float(equivalent_non_linear_strain),
"nonlinear_stress": float(equivalent_non_linear_stress),
"plastic_strain": float(plastic_strain),
}
def apply_glinka_correction(elastic_stress, youngs_modulus, non_linear_strain, non_linear_stress):
"""Update the Glinka correction estimates from the Von Mises stress.
Parameters
----------
elastic_stress : float
The linear stress input
youngs_modulus : float
The Young's modulus used in the linear stress calculation.
non_linear_strain, non_linear_stress : np.ndarray
Corresponding arrays of strain (ratio) and stress (MPa) for the non-linear correction.
"""
elastic_strain = elastic_stress / youngs_modulus
elastic_stress_strain_integral = 0.5 * elastic_strain * elastic_stress
# Check input data is sensible
check_youngs_modulus(non_linear_strain, non_linear_stress, youngs_modulus)
# If on the elastic portion, just return inputs, so we don't get interpolation errors
if elastic_stress < find_divergence_point(non_linear_strain, non_linear_stress, youngs_modulus):
return {
"nonlinear_strain": elastic_strain,
"nonlinear_stress": elastic_stress,
"plastic_strain": 0,
}
# Find a matching stress-strain integral in the non-linear response
# (This is the Glinka correction)
non_linear_stress_strain_integral = np.zeros(len(non_linear_stress))
for i in range(1, len(non_linear_stress)):
# This assumes that our stress strain curve starts at zero,
# and uses simple trapezium integration
non_linear_stress_strain_integral[i] = (
non_linear_stress_strain_integral[i - 1]
+ (non_linear_strain[i] - non_linear_strain[i - 1])
* (non_linear_stress[i] + non_linear_stress[i - 1])
/ 2
)
# Error if out of range:
if elastic_stress_strain_integral > np.max(non_linear_stress_strain_integral):
raise ValueError(
"Input too large (elastic stress strain integral > "
"maximum stress strain integral in non-linear data). "
"Elastic stress is "
+ str(elastic_stress)
+ ", elastic stress strain integral is "
+ str(elastic_stress_strain_integral)
+ ", maximum plastic stress strain integral is "
+ str(np.max(non_linear_stress_strain_integral))
)
# Lookup to find non-linear strain at matching stress strain integral:
equivalent_non_linear_strain = np.interp(
elastic_stress_strain_integral, non_linear_stress_strain_integral, non_linear_strain
)
# Lookup to find non-linear stress at this strain:
equivalent_non_linear_stress = np.interp(
equivalent_non_linear_strain, non_linear_strain, non_linear_stress
)
# Find plastic strain
plastic_strain = equivalent_non_linear_strain - elastic_strain
return {
"nonlinear_strain": float(equivalent_non_linear_strain),
"nonlinear_stress": float(equivalent_non_linear_stress),
"plastic_strain": float(plastic_strain),
}
Functions to find stress sample points for average stresses#
def find_bridge_stress_sample_points(mc):
all_points = []
# Sample points is hardcoded to 15 in Motor-CAD for stress averaging
sample_points = 15
# Check the rotor type
rotor_type = mc.get_variable("BPMRotor")
# U shape is 13, V web is 11
if rotor_type == 11:
layers = mc.get_variable("VMagnet_Layers")
elif rotor_type == 13:
layers = mc.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 = mc.get_variable("AvStressRadialLocation_Bridge")
rotor_diameter = mc.get_variable("RotorDiameter")
poles = mc.get_variable("Pole_Number")
pole_pairs = poles / 2
for layer in range(layers):
if rotor_type == 11:
# V Web template
bridge_thickness = mc.get_array_variable("BridgeThickness_Array", layer)
web_thickness = mc.get_array_variable("WebThickness_Array", layer)
pole_arc = math.radians(mc.get_array_variable("PoleArc_Array", layer) / pole_pairs)
theta_4 = math.asin(web_thickness / (2 * (rotor_diameter / 2 - bridge_thickness)))
theta_0 = math.radians(180 / poles) + pole_arc / 2
theta_1 = math.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 = mc.get_array_variable("UShape_BridgeThickness_Array", layer)
web_thickness = mc.get_array_variable("UShape_WebThickness_Array", layer)
outer_thickness = mc.get_array_variable("UShape_Thickness_Outer_Array", layer)
theta_offset = math.radians(
mc.get_array_variable("UShape_OuterAngleOffset_Array", layer)
)
inner_rad = rotor_diameter / 2 - bridge_thickness
# The start angle of the FEA model
theta_0 = math.radians(360 / poles)
# The angle to the end of the web
theta_1 = math.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 = math.radians(0.01)
iteration = 0
last_err = math.cos(theta_offset - theta_1 - test_phi / 2) * math.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 = math.cos(theta_offset - theta_1 - test_phi / 2) * math.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)
layer_points = []
for point in range(sample_points):
this_theta = theta + point * delta_theta
x = r0 * math.cos(this_theta)
y = r0 * math.sin(this_theta)
layer_points.append([x, y])
all_points.append(layer_points)
return all_points
Non-linear stress strain data for rotor material#
Our non-linear stress-strain data, in this case for a region with Young’s modulus of 185 GPa, and with plastic deformation above 480 MPa. Strain is a ratio (dimensionless), Stress is in MPa.
non_linear_strain = np.array(
[
0,
0.0002,
0.0004,
0.0006,
0.0008,
0.001,
0.0012,
0.0014,
0.0016,
0.0018,
0.002,
0.0022,
0.0024,
0.0026,
0.0028,
0.003,
0.0032,
0.0034,
0.0036,
0.0038,
0.004,
0.0042,
0.0044,
0.0046,
0.0048,
0.005,
0.0052,
0.0054,
0.0056,
0.0058,
]
)
non_linear_stress = np.array(
[
0,
37,
74,
111,
148,
185,
222,
259,
296,
333,
370,
407,
444,
481,
500,
520,
540,
560,
570,
580,
590,
600,
605,
610,
615,
620,
625,
630,
635,
640,
]
)
Main script using the functions defined previously#
Note
To run the calculation automatically when a stress calculation is completed, the code below can be placed in method called final(self) within a mechanical_stress() class. For example:
class mechanical_stress():
def final(self):
# Code below goes here
# Option whether to overwrite output values
overwrite_stress_outputs = True
# Connect to Motor-CAD
mc = 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:
mc.set_variable("MessageDisplayState", 2)
mc.load_template("e10")
# Set extreme overspeed condition to show yield
mc.set_variable("ShaftSpeed", 25000)
mc.do_mechanical_calculation()
# The material properties of the rotor
rotor_youngs_modulus = mc.get_variable("YoungsCoefficient_RotorLam")
# Find points to sample for bridge
bridge_samples = find_bridge_stress_sample_points(mc)
# Iterate over layers and then points to get results
for layer_index, layer_samples in enumerate(bridge_samples):
stresses = []
non_linear_stresses = []
for xy_point in layer_samples:
stress_von_mises = mc.get_point_value("SVM", xy_point[0], xy_point[1])[0]
non_linear_result = apply_neuber_correction(
stress_von_mises, rotor_youngs_modulus, non_linear_strain, non_linear_stress
)
stresses.append(stress_von_mises)
non_linear_stresses.append(non_linear_result["nonlinear_stress"])
# Display results for this layer
print("Point positions: " + str(layer_samples))
print("Stress at points (original) : " + str(stresses))
print("Stress at points (corrected): " + str(non_linear_stresses))
print("Mean stress (original) : " + str(sum(stresses) / len(stresses)))
corrected_mean_stress = sum(non_linear_stresses) / len(non_linear_stresses)
print("Mean stress (corrected): " + str(corrected_mean_stress))
if overwrite_stress_outputs:
mc.set_array_variable("AvStress_MagnetBridge", layer_index, corrected_mean_stress)
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 (original) : [439, 512.5, 512.5, 437.3, 578, 473.1, 571.9, 523.4, 523.4, 632.3, 650, 650, 505.2, 505.2, 429.5]
Stress at points (corrected): [439, 502.47043918918916, 502.47043918918916, 437.3, 548.8476658476659, 473.1, 544.5390724815725, 510.0997027027027, 510.0997027027027, 577.1775568990042, 585.1143451143452, 585.1143451143452, 497.4063374217591, 497.4063374217591, 429.5]
Mean stress (original) : 529.5533333333334
Mean stress (corrected): 509.30972960561564
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 (original) : [432.6, 432.6, 380.3, 503.5, 503.5, 424.6, 538.2, 538.2, 433.1, 433.1, 598.3, 598.3, 325.6, 328.3, 328.3]
Stress at points (corrected): [432.6, 432.6, 380.3, 496.2275317486161, 496.2275317486161, 424.6, 520.6816216216216, 520.6816216216216, 433.1, 433.1, 562.0901716581446, 562.0901716581446, 325.6, 328.3, 328.3]
Mean stress (original) : 453.23333333333335
Mean stress (corrected): 445.0999100037843
Total running time of the script: (0 minutes 40.401 seconds)