Dear all;
I am making a rotor for an induction eletric machine so there is the outer circle and the iner hole where the shaft goes through. Also there is the induction rotor bars for that I made a code that generates a circular padron to automaticaly generate them but I am new to Gmsh and never handle mutiple surfaces and subdomains. They are being bunch up togheter afther the first one and if I mesh it the rotor bar geometries are being ignored the error is : Error Error Error
: Unknown OpenCASCADE curve with tag 5
: ‘C:\Users\Usuario\Downloads\padrao circular.geo’, line 41: Could not add curve loop : Unknown OpenCASCADE curve loop with tag 16
Here is the output geo file (it continues the padron afther the second rotor):
SetFactory("OpenCASCADE");
// Define Constants
r_rlm = 0.31; // Radius of the rotor
r_shaft = 0.08; // Radius of the shaft
rbr_width1 = 0.0143; // First width of the rotor bar
rbr_width2 = 0.03; // Second width of the rotor bar
rbr_height = 0.0935; // Height of the rotor bar
rbr_y = 0.1786; // Y coordinate of the center of rotor bars
rbr_y2 = 0.0834; // Y coordinate of the displaced center of rotor bars
Circle(1) = {0, 0, 0, 0.31, 0, 2*Pi}; // Rotor circle, centered at (0, 0)
Circle(2) = {0, 0, 0, 0.08, 0, 2*Pi}; // Shaft circle, centered at (0, 0)
// Create line loops
Line Loop(3) = 1; // Outer circle
Line Loop(4) = -2; // Inner circle, note the negative sign
// Create a plane surface for the main domain
Plane Surface(5) = {3, 4};
// Assign a physical surface for the main domain
Physical Surface("MainDomain") = {5};
// Rotor 1 angle = pi/8.0*0
Point(6) = {0.00715, 0.1786, 0, 1.0}; // Point on the right side
Point(7) = {-0.00715, 0.1786, 0, 1.0}; // Point on the left side
Point(8) = {0.015, 0.262, 0, 1.0}; // Point on the right side
Point(9) = {-0.015, 0.262, 0, 1.0}; // Point on the left side
Point(10) = {0.0, 0.2721, 0, 1.0}; // Top point
//
Line(11) = {7, 6};
Line(12) = {6, 8};
Line(13) = {8, 10};
Line(14) = {10, 9};
Line(15) = {9, 7};
Line Loop(16) = {11, 12, 13, -5, -6};
Plane Surface(17) = {16};
Physical Surface("Rotor1") = {17};
Physical Surface("Rotor1_Subdomain") = 17;
// End Rotor 1
// Rotor 2 angle = pi/8.0*1
Point(18) = {-0.00715, 0.1786, 0, 1.0}; // Point on the right side
Point(19) = {0.015, 0.262, 0, 1.0}; // Point on the left side
Point(20) = {-0.015, 0.262, 0, 1.0}; // Point on the right side
Point(21) = {0.0, 0.2721, 0, 1.0}; // Point on the left side
Point(22) = {-0.06174152236294933, 0.1677410710479262, 0, 1.0}; // Top point
//
Line(23) = {19, 18};
Line(24) = {18, 20};
Line(25) = {20, 22};
Line(26) = {22, 21};
Line(27) = {21, 19};
Line Loop(28) = {23, 24, 25, -5, -6};
Plane Surface(29) = {28};
Physical Surface("Rotor2") = {29};
Physical Surface("Rotor2_Subdomain") = 29;
// End Rotor 2
And here is the python code :
import math
import matplotlib.pyplot as plt
import numpy as np
# Define the initial and final angles in terms of pi notation
n=16/2 #numero de barras de rotor
initial_angle_pi = 1/n # Initial angle in pi notation
final_angle_pi = 2 # Final angle in pi notation
step_angle_pi = 1/n # Angle increment per step in pi notation
# Convert pi notation angles to radians
initial_angle = initial_angle_pi * math.pi
final_angle = final_angle_pi * math.pi
step_angle = step_angle_pi * math.pi
def cartesian_to_polar(x, y):
"""
Transforms Cartesian coordinates to polar coordinates.
Args:
x: The x-coordinate.
y: The y-coordinate.
Returns:
A tuple containing the polar coordinates (radius, angle in radians).
"""
radius = math.sqrt(x**2 + y**2)
angle = math.atan2(y, x)
return radius, angle
def polar_to_cartesian(radius, angle):
"""
Converts polar coordinates to Cartesian coordinates.
Args:
radius: The radius in polar coordinates.
angle: The angle in radians in polar coordinates.
Returns:
A tuple containing the Cartesian coordinates (x, y).
"""
x = radius * math.cos(angle)
y = radius * math.sin(angle)
return x, y
r_rlm = 0.31 # Radius of the rotor
r_shaft = 0.08 # Radius of the shaft
rbr_width1 = 0.0143 # First width of the rotor bar
rbr_width2 = 0.03 # Second width of the rotor bar
rbr_height = 0.0935 # Height of the rotor bar
rbr_y = 0.1786 # Y coordinate of the center of rotor bars
rbr_y2 = 0.0834 # Y coordinate of the displaced center of rotor bars
# Create the points
Point_3 = (rbr_width1/2, rbr_y) # Point on the right side
Point_4 = (-rbr_width1/2, rbr_y) # Point on the left side
Point_5 = (rbr_width2/2, rbr_y+rbr_y2) # Point on the right side
Point_6 = (-rbr_width2/2, rbr_y+rbr_y2) # Point on the left side
Point_7 = (0, rbr_y + rbr_height) # Top point
# Create the array of points
points = [Point_3, Point_4, Point_5, Point_6, Point_7]
# Convert the points to polar coordinates
polar_points = []
for point in points:
radius, angle = cartesian_to_polar(point[0], point[1])
polar_points.append((radius, angle))
rotated_points = []
# Append original points to rotated_points
for point in points:
rotated_points.append(point)
# Rotate the points
current_angle = initial_angle
while current_angle <= final_angle:
for i in range(len(points)):
radius, angle = polar_points[i]
rotated_x, rotated_y = polar_to_cartesian(radius, angle + current_angle)
rotated_points.append((rotated_x, rotated_y))
current_angle += step_angle
points.clear()
# Convert the points back to cartesian coordinates
for point in rotated_points:
x, y = point
points.append((x, y))
points = np.array(points)
x_coords = points[:, 0]
y_coords = points[:, 1]
# Plot the rotated points
plt.plot(points[:, 0], points[:, 1], 'o')
plt.xlim(-r_rlm, r_rlm)
plt.ylim(-r_rlm, r_rlm)
plt.show()
output_text = ""
# Initial block of text
initial_block = f"""
SetFactory("OpenCASCADE");
// Define Constants
r_rlm = {r_rlm}; // Radius of the rotor
r_shaft = {r_shaft}; // Radius of the shaft
rbr_width1 = {rbr_width1}; // First width of the rotor bar
rbr_width2 = {rbr_width2}; // Second width of the rotor bar
rbr_height = {rbr_height}; // Height of the rotor bar
rbr_y = {rbr_y}; // Y coordinate of the center of rotor bars
rbr_y2 = {rbr_y2}; // Y coordinate of the displaced center of rotor bars
Circle(1) = {{0, 0, 0, {r_rlm}, 0, 2*Pi}}; // Rotor circle, centered at (0, 0)
Circle(2) = {{0, 0, 0, {r_shaft}, 0, 2*Pi}}; // Shaft circle, centered at (0, 0)
// Create line loops
Line Loop(3) = 1; // Outer circle
Line Loop(4) = -2; // Inner circle, note the negative sign
// Create a plane surface for the main domain
Plane Surface(5) = {{3, 4}};
// Assign a physical surface for the main domain
Physical Surface("MainDomain") = {{5}};
"""
output_text += initial_block
# Loop through angles
index = 0
current_angle = initial_angle
angle_index = 1 # Starting index
point_index = 6 # Starting point index
while current_angle <= final_angle:
cos_val = math.cos(current_angle)
sin_val = math.sin(current_angle)
# point block
point_block = f"""
// Rotor {angle_index} angle = pi/{n}*{angle_index-1}
Point({point_index}) = {{{points[index, 0]}, {points[index, 1]}, 0, 1.0}}; // Point on the right side
Point({point_index+1}) = {{{points[index+1, 0]}, {points[index+1, 1]}, 0, 1.0}}; // Point on the left side
Point({point_index+2}) = {{{points[index+2, 0]}, {points[index+2, 1]}, 0, 1.0}}; // Point on the right side
Point({point_index+3}) = {{{points[index+3, 0]}, {points[index+3, 1]}, 0, 1.0}}; // Point on the left side
Point({point_index+4}) = {{{points[index+4, 0]}, {points[index+4, 1]}, 0, 1.0}}; // Top point
"""
output_text += point_block
# line block and surfaces
line_block = f"""
//
Line({point_index+5}) = {{{point_index+1}, {point_index}}};
Line({point_index+6}) = {{{point_index}, {point_index+2}}};
Line({point_index+7}) = {{{point_index+2}, {point_index+4}}};
Line({point_index+8}) = {{{point_index+4}, {point_index+3}}};
Line({point_index+9}) = {{{point_index+3}, {point_index+1}}};
Line Loop({point_index+10}) = {{{point_index+5}, {point_index+6}, {point_index+7}, -5, -6}};
Plane Surface({point_index+11}) = {{{point_index+10}}};
Physical Surface("Rotor{angle_index}") = {{{point_index+11}}};
Physical Surface("Rotor{angle_index}_Subdomain") = {point_index+11};
"""
output_text += line_block
output_text += f"\n// End Rotor {angle_index}\n"
index += 1
angle_index += 1
point_index += 12
current_angle += step_angle
# Print the output
print(output_text)
# Write the output to a file
with open("padrao circular.geo", "w") as file:
file.write(output_text)
I tried mutiple ways do separate into surfaces and check the indices and the order for anti clock wise rotation. I don’t have experience with gmsh errors. But I am glad to learn, it must have better ways to do what I have done so far. Also any sugestion into the 2D mesh genaration would be more them welcome.