Lua automation script for simulation of radial magnetic bearings

 -- Using manual4.2 pdf document to write this script

-- Run command:

-- Home Gaming laptop

--runCommand       = C:\femm42\bin\femm.exe -lua-script=C:\06_NAL\2024_RC_AMB\02_elecDesign\2024oct30_radialActuator\2024oct30_FEMMsimulation\2024nov12a_radialAMB.lua

mydir              = "C:\\06_NAL\\2024_RC_AMB\\02_elecDesign\\2024oct30_radialActuator\\2024oct30_FEMMsimulation\\"

outputFemmFilePath = "2024nov12a_ramb.fem";

--mydir = "C:\\Users\\vigne\\OneDrive\\Desktop\\2024nov03_femm\\"

--mydir = "F:\\2024_RC_amb\\02_elecDesign\\02_radialActuator\\2024oct30_FEMMsimulation\\"

-- Vicky computer

-- C:\femm42\bin\femm.exe -lua-script=C:\Users\vigne\OneDrive\Desktop\2024nov03_femm\2024nov07b_radialAMB.lua

-- Office workstation

-- C:\femm42\bin\femm.exe -lua-script=F:\2024_RC_amb\02_elecDesign\02_radialActuator\2024oct30_FEMMsimulation\2024nov07b_radialAMB.lua

-- Reading from manual 4.2.pdf in Vikki house, arcot;

-- Hints from manual

-- block labels denote what material properties and  mesh size are associated with each solution region

-- group mode, is meant to glue different objects together into parts so that entire parts can  be manipulated more easily

-- Prescribed A  boundary condition, the vector potential, A, is prescribed  along a given boundary, can be made zero on the outer od of the ramb stator to make flux not leave the bearing stator.

-- Lua variables are global, unless you specify explicitly that they are local!

-- https://www.codecademy.com/resources/docs/lua/variables

function mirrorAbt_67p5(inp_x1_x2,inp_y1_y2)

-- Get the mirrored point, by rotating to align the mirror line to x axis, getting reflection, and rotating back.

-- First rotate to align the mirror line to x axis

local temp_x1_x2 = inp_x1_x2 * cos(rad(-67.5 )) - inp_y1_y2  * sin(rad(-67.5 ));

local temp_y1_y2 = inp_x1_x2 * sin(rad(-67.5 )) + inp_y1_y2  * cos(rad(-67.5 ));

-- Flip the y axis for mirroring about x axis

temp_y1_y2 = -temp_y1_y2;

-- Rotate back to the original mirror line

local return_x1_x2= temp_x1_x2 * cos(rad(67.5 )) - temp_y1_y2  * sin(rad(67.5 ));

local return_y1_y2= temp_x1_x2 * sin(rad(67.5 )) + temp_y1_y2  * cos(rad(67.5 ));

-- return

return return_x1_x2,return_y1_y2;

end

function createGeometry()

--Create the shaft circle

local angle_increment = 360 / numSegs_noUni;

for i=0,numSegs_noUni-1 do

-- Calculate angles in radians for the start and end points of each segment

local angle1 = i * angle_increment * (pi / 180)

local angle2 = (i + 1) * angle_increment * (pi / 180)

-- Calculate the (x, y) coordinates for the start and end points

local x1 = 0.5*shaftOD_m * cos(angle1)

local y1 = 0.5*shaftOD_m  * sin(angle1)

local x2 = 0.5*shaftOD_m  * cos(angle2)

local y2 = 0.5*shaftOD_m  * sin(angle2)

-- Add nodes at the start and end points

mi_addnode(x1, y1)

mi_addnode(x2, y2)

-- Draw an arc segment between the two nodes to form the circle

mi_addarc(x1, y1, x2, y2, angle_increment, numSegs_noUni)

end

mi_zoomnatural();

--Create the  outer radius of the stator OD

angle_increment = 360 / 2;

for i=0,numSegs_noUni-1 do

-- Calculate angles in radians for the start and end points of each segment

local angle1 = i * angle_increment * (pi / 180)

local angle2 = (i + 1) * angle_increment * (pi / 180)

-- Calculate the (x, y) coordinates for the start and end points

local x1 = 0.5*odOuter_m * cos(angle1)

local y1 = 0.5*odOuter_m * sin(angle1)

local x2 = 0.5*odOuter_m * cos(angle2)

local y2 = 0.5*odOuter_m * sin(angle2)

-- Add nodes at the start and end points

mi_addnode(x1, y1)

mi_addnode(x2, y2)

-- Draw an arc segment between the two nodes to form the circle

mi_addarc(x1, y1, x2, y2, angle_increment, numSegs_noUni)

end

mi_zoomnatural();

--Create the inner radius of the stator  OD

local odInnAngArr_deg={   -22.5+thetaLegSubtOdInner_deg/2, 22.5-thetaLegSubtOdInner_deg/2, 22.5+thetaLegSubtOdInner_deg/2, 67.5-thetaLegSubtOdInner_deg/2, 67.5+thetaLegSubtOdInner_deg/2, 112.5-thetaLegSubtOdInner_deg/2, 112.5+thetaLegSubtOdInner_deg/2, 157.5-thetaLegSubtOdInner_deg/2, 157.5+thetaLegSubtOdInner_deg/2, 202.5-thetaLegSubtOdInner_deg/2, 202.5+thetaLegSubtOdInner_deg/2, 247.5-thetaLegSubtOdInner_deg/2, 247.5+thetaLegSubtOdInner_deg/2, 292.5-thetaLegSubtOdInner_deg/2, 292.5+thetaLegSubtOdInner_deg/2, 337.5-thetaLegSubtOdInner_deg/2, 337.5+thetaLegSubtOdInner_deg/2}

-- Draw the arc segments between legs

for i=1,16,2 do

local x1 = 0.5*odInner_m * cos((pi/180)*odInnAngArr_deg[i])

local y1 = 0.5*odInner_m * sin((pi/180)*odInnAngArr_deg[i])

local x2 = 0.5*odInner_m * cos((pi/180)*odInnAngArr_deg[i+1])

local y2 = 0.5*odInner_m * sin((pi/180)*odInnAngArr_deg[i+1])

-- Add nodes at the start and end points

mi_addnode(x1, y1)

mi_addnode(x2, y2)

-- Draw an arc segment between the two nodes to form the circle

mi_addarc(x1, y1, x2, y2, (odInnAngArr_deg[i+1]-odInnAngArr_deg[i]), numSegs_noUni)

end

--Create the inner radius of the OD

local idInnAngArr_deg={

   -22.5+thetaLegSubtIdInner_deg/2,

22.5-thetaLegSubtIdInner_deg/2,

22.5+thetaLegSubtIdInner_deg/2,

67.5-thetaLegSubtIdInner_deg/2,

67.5+thetaLegSubtIdInner_deg/2,

   112.5-thetaLegSubtIdInner_deg/2,

   112.5+thetaLegSubtIdInner_deg/2,

   157.5-thetaLegSubtIdInner_deg/2,

   157.5+thetaLegSubtIdInner_deg/2,

   202.5-thetaLegSubtIdInner_deg/2,

   202.5+thetaLegSubtIdInner_deg/2,

   247.5-thetaLegSubtIdInner_deg/2,

   247.5+thetaLegSubtIdInner_deg/2,

   292.5-thetaLegSubtIdInner_deg/2,

   292.5+thetaLegSubtIdInner_deg/2,

   337.5-thetaLegSubtIdInner_deg/2,

   337.5+thetaLegSubtIdInner_deg/2

   }

-- Draw the arc segments at the tips of the legs

for i=2,17,2 do

local x1 = 0.5*idInner_m * cos((pi/180)*idInnAngArr_deg[i])

local y1 = 0.5*idInner_m * sin((pi/180)*idInnAngArr_deg[i])

local x2 = 0.5*idInner_m * cos((pi/180)*idInnAngArr_deg[i+1])

local y2 = 0.5*idInner_m * sin((pi/180)*idInnAngArr_deg[i+1])

-- Add nodes at the start and end points

mi_addnode(x1, y1)

mi_addnode(x2, y2)

-- Draw an arc segment between the two nodes to form the circle

mi_addarc(x1, y1, x2, y2, (idInnAngArr_deg[i+1]-idInnAngArr_deg[i]), numSegs_noUni)

end

-- Draw the parallel lines of  the legs

for i=1,16,1 do

local x1 = 0.5*idInner_m * cos((pi/180)*idInnAngArr_deg[i])

local y1 = 0.5*idInner_m * sin((pi/180)*idInnAngArr_deg[i])

local x2 = 0.5*odInner_m * cos((pi/180)*odInnAngArr_deg[i])

local y2 = 0.5*odInner_m * sin((pi/180)*odInnAngArr_deg[i])

-- Add nodes at the start and end points

mi_addnode(x1, y1)

mi_addnode(x2, y2)

-- Draw an arc segment between the two nodes to form the circle

mi_addsegment(x1, y1, x2, y2)

end

-- Create the top quadrant, right leg, left side coil manually first.

local x1coil_TopQuad_RightLeg_LeftSide = 0.01531; local y1coil_TopQuad_RightLeg_LeftSide = 0.07484;

local x2coil_TopQuad_RightLeg_LeftSide = 0.02039; local y2coil_TopQuad_RightLeg_LeftSide = 0.07274;

local x3coil_TopQuad_RightLeg_LeftSide = 0.00642; local y3coil_TopQuad_RightLeg_LeftSide = 0.03901;

local x4coil_TopQuad_RightLeg_LeftSide = 0.00134; local y4coil_TopQuad_RightLeg_LeftSide = 0.04112;

mi_addnode(x1coil_TopQuad_RightLeg_LeftSide,y1coil_TopQuad_RightLeg_LeftSide);

mi_addnode(x2coil_TopQuad_RightLeg_LeftSide,y2coil_TopQuad_RightLeg_LeftSide);

mi_addnode(x3coil_TopQuad_RightLeg_LeftSide,y3coil_TopQuad_RightLeg_LeftSide);

mi_addnode(x4coil_TopQuad_RightLeg_LeftSide,y4coil_TopQuad_RightLeg_LeftSide);

mi_addsegment(x1coil_TopQuad_RightLeg_LeftSide,y1coil_TopQuad_RightLeg_LeftSide,x2coil_TopQuad_RightLeg_LeftSide,y2coil_TopQuad_RightLeg_LeftSide);

mi_addsegment(x2coil_TopQuad_RightLeg_LeftSide,y2coil_TopQuad_RightLeg_LeftSide,x3coil_TopQuad_RightLeg_LeftSide,y3coil_TopQuad_RightLeg_LeftSide);

mi_addsegment(x3coil_TopQuad_RightLeg_LeftSide,y3coil_TopQuad_RightLeg_LeftSide,x4coil_TopQuad_RightLeg_LeftSide,y4coil_TopQuad_RightLeg_LeftSide);

mi_addsegment(x4coil_TopQuad_RightLeg_LeftSide,y4coil_TopQuad_RightLeg_LeftSide,x1coil_TopQuad_RightLeg_LeftSide,y1coil_TopQuad_RightLeg_LeftSide);

-- Get mid points of the coil to select the segments later for rotating.

midPnt_x1_x2 = 0.5*(x1coil_TopQuad_RightLeg_LeftSide+x2coil_TopQuad_RightLeg_LeftSide);

midPnt_x2_x3 = 0.5*(x2coil_TopQuad_RightLeg_LeftSide+x3coil_TopQuad_RightLeg_LeftSide);

midPnt_x3_x4 = 0.5*(x3coil_TopQuad_RightLeg_LeftSide+x4coil_TopQuad_RightLeg_LeftSide);

midPnt_x4_x1 = 0.5*(x4coil_TopQuad_RightLeg_LeftSide+x1coil_TopQuad_RightLeg_LeftSide);

midPnt_y1_y2 = 0.5*(y1coil_TopQuad_RightLeg_LeftSide+y2coil_TopQuad_RightLeg_LeftSide);

midPnt_y2_y3 = 0.5*(y2coil_TopQuad_RightLeg_LeftSide+y3coil_TopQuad_RightLeg_LeftSide);

midPnt_y3_y4 = 0.5*(y3coil_TopQuad_RightLeg_LeftSide+y4coil_TopQuad_RightLeg_LeftSide);

midPnt_y4_y1 = 0.5*(y4coil_TopQuad_RightLeg_LeftSide+y1coil_TopQuad_RightLeg_LeftSide);

-- Mirror left side to get the right side of the top quadrant right leg to make the other side of the top right coil

    mi_clearselected();

mi_selectsegment(midPnt_x1_x2,midPnt_y1_y2);

mi_selectsegment(midPnt_x2_x3,midPnt_y2_y3);

mi_selectsegment(midPnt_x3_x4,midPnt_y3_y4);

mi_selectsegment(midPnt_x4_x1,midPnt_y4_y1);

local mirrorLineSegments = 1;

    mi_mirror(0,0,cos(rad(67.5)),sin(rad(67.5)),mirrorLineSegments )

    

-- Rotate the left side coil 8 times

    mi_clearselected();

mi_selectsegment(midPnt_x1_x2,midPnt_y1_y2);

mi_selectsegment(midPnt_x2_x3,midPnt_y2_y3);

mi_selectsegment(midPnt_x3_x4,midPnt_y3_y4);

mi_selectsegment(midPnt_x4_x1,midPnt_y4_y1);

mi_copyrotate(0,0,45,8,1)

mi_clearselected();

    -- Create the mid point of top quad, right leg, left coil, required for assigning material.

    avg_x_TopQuad_RightLeg_LeftSide  = 0.25*(x1coil_TopQuad_RightLeg_LeftSide+x2coil_TopQuad_RightLeg_LeftSide+x3coil_TopQuad_RightLeg_LeftSide+x4coil_TopQuad_RightLeg_LeftSide);

    avg_y_TopQuad_RightLeg_LeftSide  = 0.25*(y1coil_TopQuad_RightLeg_LeftSide+y2coil_TopQuad_RightLeg_LeftSide+y3coil_TopQuad_RightLeg_LeftSide+y4coil_TopQuad_RightLeg_LeftSide);

-- Rotate to get the mid point of the top quad, right leg, right coil.

    local rotAngleInsideSameCoil_rad = -2*(atan(avg_y_TopQuad_RightLeg_LeftSide/avg_x_TopQuad_RightLeg_LeftSide)-rad(67.5));

    avg_x_TopQuad_RightLeg_RightSide = avg_x_TopQuad_RightLeg_LeftSide * cos((rotAngleInsideSameCoil_rad )) - avg_y_TopQuad_RightLeg_LeftSide * sin((rotAngleInsideSameCoil_rad ));

avg_y_TopQuad_RightLeg_RightSide = avg_x_TopQuad_RightLeg_LeftSide * sin((rotAngleInsideSameCoil_rad )) + avg_y_TopQuad_RightLeg_LeftSide * cos((rotAngleInsideSameCoil_rad ));

print("Top quadrant right leg left coil avg  (" .. avg_x_TopQuad_RightLeg_LeftSide   .. ", " ..  avg_y_TopQuad_RightLeg_LeftSide  .. ")")

print("Rotation angle = " .. deg(rotAngleInsideSameCoil_rad));

print("Top quadrant right leg right coil avg (" .. avg_x_TopQuad_RightLeg_RightSide  .. ", " ..  avg_y_TopQuad_RightLeg_RightSide .. ")")

-- Mirror the mid points to get their coordinates to select them later

mirrored_midPnt_x1_x2,mirrored_midPnt_y1_y2 = mirrorAbt_67p5(midPnt_x1_x2,midPnt_y1_y2);

mirrored_midPnt_x2_x3,mirrored_midPnt_y2_y3 = mirrorAbt_67p5(midPnt_x2_x3,midPnt_y2_y3);

mirrored_midPnt_x3_x4,mirrored_midPnt_y3_y4 = mirrorAbt_67p5(midPnt_x3_x4,midPnt_y3_y4);

mirrored_midPnt_x4_x1,mirrored_midPnt_y4_y1 = mirrorAbt_67p5(midPnt_x4_x1,midPnt_y4_y1);

    -- Rotate the mirrored coil 8 times

    mi_clearselected();

mi_selectsegment(mirrored_midPnt_x1_x2,mirrored_midPnt_y1_y2);

mi_selectsegment(mirrored_midPnt_x2_x3,mirrored_midPnt_y2_y3);

mi_selectsegment(mirrored_midPnt_x3_x4,mirrored_midPnt_y3_y4);

mi_selectsegment(mirrored_midPnt_x4_x1,mirrored_midPnt_y4_y1);

mi_copyrotate(0,0,45,8,1)

mi_clearselected();

-- Get the mid points of the top quad, left leg, left and right coils

    rotAngle_deg = 45;

avg_x_TopQuad_LeftLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg));

avg_y_TopQuad_LeftLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg));

avg_x_TopQuad_LeftLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_TopQuad_LeftLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg));

print("Top quadrant left leg left side       (" .. avg_x_TopQuad_LeftLeg_LeftSide    .. ", " ..  avg_y_TopQuad_LeftLeg_LeftSide   .. ")")

print("Top quadrant left leg right side      (" .. avg_x_TopQuad_LeftLeg_RightSide   .. ", " ..  avg_y_TopQuad_LeftLeg_RightSide  .. ")")

rotAngle_deg = 90;

avg_x_LeftQuad_TopLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg))  - avg_y_TopQuad_RightLeg_LeftSide * sin(rad(rotAngle_deg));

avg_y_LeftQuad_TopLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg))  + avg_y_TopQuad_RightLeg_LeftSide * cos(rad(rotAngle_deg));

avg_x_LeftQuad_TopLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg))  - avg_y_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_LeftQuad_TopLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg))  + avg_y_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg));

avg_x_LeftQuad_BotLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * cos(rad(rotAngle_deg))  - avg_y_TopQuad_LeftLeg_LeftSide * sin(rad(rotAngle_deg));

avg_y_LeftQuad_BotLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * sin(rad(rotAngle_deg))  + avg_y_TopQuad_LeftLeg_LeftSide * cos(rad(rotAngle_deg));

avg_x_LeftQuad_BotLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * cos(rad(rotAngle_deg))  - avg_y_TopQuad_LeftLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_LeftQuad_BotLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * sin(rad(rotAngle_deg))  + avg_y_TopQuad_LeftLeg_RightSide * cos(rad(rotAngle_deg));

print("Left quadrant top    leg  left side (" .. avg_x_LeftQuad_TopLeg_LeftSide    .. ", " .. avg_y_LeftQuad_TopLeg_LeftSide  .. ")")

print("Left quadrant top    leg right side (" .. avg_x_LeftQuad_TopLeg_RightSide   .. ", " .. avg_y_LeftQuad_TopLeg_RightSide .. ")")

print("Left quadrant bottom leg  left side (" .. avg_x_LeftQuad_BotLeg_LeftSide    .. ", " .. avg_y_LeftQuad_BotLeg_LeftSide  .. ")")

print("Left quadrant bottom leg right side (" .. avg_x_LeftQuad_BotLeg_RightSide   .. ", " .. avg_y_LeftQuad_BotLeg_RightSide .. ")")

rotAngle_deg = 180;

avg_x_BotQuad_LeftLeg_LeftSide   = avg_x_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg));

avg_y_BotQuad_LeftLeg_LeftSide   = avg_x_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg));

avg_x_BotQuad_LeftLeg_RightSide  = avg_x_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_BotQuad_LeftLeg_RightSide  = avg_x_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg));

avg_x_BotQuad_RightLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * cos(rad(rotAngle_deg)) - avg_y_TopQuad_LeftLeg_LeftSide   * sin(rad(rotAngle_deg));

avg_y_BotQuad_RightLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * sin(rad(rotAngle_deg)) + avg_y_TopQuad_LeftLeg_LeftSide   * cos(rad(rotAngle_deg));

avg_x_BotQuad_RightLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * cos(rad(rotAngle_deg)) - avg_y_TopQuad_LeftLeg_RightSide  * sin(rad(rotAngle_deg));

avg_y_BotQuad_RightLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * sin(rad(rotAngle_deg)) + avg_y_TopQuad_LeftLeg_RightSide  * cos(rad(rotAngle_deg));

print("Bottom quadrant  left leg  left side coil avg (" .. avg_x_BotQuad_LeftLeg_LeftSide  .. ", " .. avg_y_BotQuad_LeftLeg_LeftSide   .. ")")

print("Bottom quadrant  left leg right Side coil avg (" .. avg_x_BotQuad_LeftLeg_RightSide .. ", " .. avg_y_BotQuad_LeftLeg_RightSide  .. ")")

print("Bottom quadrant right leg  left side coil avg (" .. avg_x_BotQuad_RightLeg_LeftSide .. ", " .. avg_y_BotQuad_RightLeg_LeftSide  .. ")")

print("Bottom quadrant right leg right Side coil avg (" .. avg_x_BotQuad_RightLeg_LeftSide .. ", " .. avg_y_BotQuad_RightLeg_RightSide .. ")")

rotAngle_deg = 270;

avg_x_RightQuad_TopLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * cos(rad(rotAngle_deg)) - avg_y_TopQuad_LeftLeg_LeftSide  * sin(rad(rotAngle_deg));

avg_y_RightQuad_TopLeg_LeftSide  = avg_x_TopQuad_LeftLeg_LeftSide   * sin(rad(rotAngle_deg)) + avg_y_TopQuad_LeftLeg_LeftSide  * cos(rad(rotAngle_deg));

avg_x_RightQuad_TopLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * cos(rad(rotAngle_deg)) - avg_y_TopQuad_LeftLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_RightQuad_TopLeg_RightSide = avg_x_TopQuad_LeftLeg_RightSide  * sin(rad(rotAngle_deg)) + avg_y_TopQuad_LeftLeg_RightSide * cos(rad(rotAngle_deg));

avg_x_RightQuad_BotLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg));

avg_y_RightQuad_BotLeg_LeftSide  = avg_x_TopQuad_RightLeg_LeftSide  * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_LeftSide  * cos(rad(rotAngle_deg));

avg_x_RightQuad_BotLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg)) - avg_y_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg));

avg_y_RightQuad_BotLeg_RightSide = avg_x_TopQuad_RightLeg_RightSide * sin(rad(rotAngle_deg)) + avg_y_TopQuad_RightLeg_RightSide * cos(rad(rotAngle_deg));

print("Right quadrant top    leg  left side (" .. avg_x_RightQuad_TopLeg_LeftSide      .. ", " .. avg_y_RightQuad_TopLeg_LeftSide  .. ")")

print("Right quadrant top    leg right side (" .. avg_x_RightQuad_TopLeg_RightSide     .. ", " .. avg_y_RightQuad_TopLeg_RightSide .. ")")

print("Right quadrant bottom leg  left side (" .. avg_x_RightQuad_BotLeg_LeftSide      .. ", " .. avg_y_RightQuad_BotLeg_LeftSide  .. ")")

print("Right quadrant bottom leg right side (" .. avg_x_RightQuad_BotLeg_RightSide    .. ", " .. avg_y_RightQuad_BotLeg_RightSide .. ")")

end

-- -----------------------------------------------------------------------------------------------------------

showconsole()

--clearconsole()

--hideconsole()

--quit()

newdocument(0) --magnetics problem

setcompatibilitymode(0) -- set compatability mode to FEMM 4.2, the latest.

-- Page 86/161 femm manual

mi_probdef(0,"meters","planar",1E-8,40e-3)

-- 0 Hz DC

-- 40e-3 is the depth into the plane of the paper

-- Problem type: "planar" or "axi"

messagebox("Balaji The GREATESTEST\nRadial active magnetic bearing simulation\n Nov10a\n Flux and rotor forces for a given current")

-- Geometry parameters

legWidth_m   = 12e-3;

odOuter_m    = 200e-3;

odInner_m    = 160e-3;

idInner_m    = 60e-3;

shaftOD_m    = 59e-3;

asslyHoles_m = 6e-3;

Nperleg_turns= 50;

numSegs_noUni= 8;-- number of segments in the circle, for getting smooth circles.

-- Angle of subtended by one leg 

-- legWidth_m = 2*R*sin(theta/2), From right angle triangle with hypotenuse as radius

thetaLegSubtOdInner_rad = 2*asin(legWidth_m/(odInner_m));

thetaLegSubtOdOuter_rad = 2*asin(legWidth_m/(odOuter_m));

thetaLegSubtIdInner_rad = 2*asin(legWidth_m/(idInner_m));

thetaLegSubtOdInner_deg = thetaLegSubtOdInner_rad*180/pi;

thetaLegSubtOdOuter_deg = thetaLegSubtOdOuter_rad*180/pi;

thetaLegSubtIdInner_deg = thetaLegSubtIdInner_rad*180/pi;

-- Create the Geometry

createGeometry()

-- Use the inbuilt materials

mi_getmaterial('Air')

mi_getmaterial('18 AWG')

mi_getmaterial('Cold rolled low carbon strip steel')

mi_getmaterial('M-43 Steel')-- EN24 is 4340 USA/ASTM in equivalent american grading chart

-- Create block labels

mi_addblocklabel(0.035,0)-- Air

mi_addblocklabel(avg_x_TopQuad_RightLeg_LeftSide   ,  avg_y_TopQuad_RightLeg_LeftSide );-- copper windings

mi_addblocklabel(avg_x_TopQuad_RightLeg_RightSide  ,  avg_y_TopQuad_RightLeg_RightSide);-- copper windings

mi_addblocklabel(avg_x_TopQuad_LeftLeg_LeftSide    ,  avg_y_TopQuad_LeftLeg_LeftSide  );-- copper windings

mi_addblocklabel(avg_x_TopQuad_LeftLeg_RightSide   ,  avg_y_TopQuad_LeftLeg_RightSide );-- copper windings

mi_addblocklabel(avg_x_LeftQuad_TopLeg_LeftSide    ,  avg_y_LeftQuad_TopLeg_LeftSide  );-- copper windings

mi_addblocklabel(avg_x_LeftQuad_TopLeg_RightSide   ,  avg_y_LeftQuad_TopLeg_RightSide );-- copper windings

mi_addblocklabel(avg_x_LeftQuad_BotLeg_LeftSide    ,  avg_y_LeftQuad_BotLeg_LeftSide  );-- copper windings

mi_addblocklabel(avg_x_LeftQuad_BotLeg_RightSide   ,  avg_y_LeftQuad_BotLeg_RightSide );-- copper windings

mi_addblocklabel(avg_x_BotQuad_LeftLeg_LeftSide    ,  avg_y_BotQuad_LeftLeg_LeftSide  );-- copper windings

mi_addblocklabel(avg_x_BotQuad_LeftLeg_RightSide   ,  avg_y_BotQuad_LeftLeg_RightSide );-- copper windings

mi_addblocklabel(avg_x_BotQuad_RightLeg_LeftSide   ,  avg_y_BotQuad_RightLeg_LeftSide );-- copper windings

mi_addblocklabel(avg_x_BotQuad_RightLeg_RightSide  ,  avg_y_BotQuad_RightLeg_RightSide);-- copper windings

mi_addblocklabel(avg_x_RightQuad_TopLeg_LeftSide   ,  avg_y_RightQuad_TopLeg_LeftSide );-- copper windings

mi_addblocklabel(avg_x_RightQuad_TopLeg_RightSide  ,  avg_y_RightQuad_TopLeg_RightSide);-- copper windings

mi_addblocklabel(avg_x_RightQuad_BotLeg_LeftSide   ,  avg_y_RightQuad_BotLeg_LeftSide );-- copper windings

mi_addblocklabel(avg_x_RightQuad_BotLeg_RightSide  ,  avg_y_RightQuad_BotLeg_RightSide);-- copper windings

mi_addblocklabel(0.09,0)-- Laminations stack

mi_addblocklabel(0,0)-- Steel EN24 shaft

--Add a "circuit property" so that we can calculate the properties of the

--coil as seen from the terminals.

seriesConnection   = 1;

parallelConnection = 0;

iTopAct_Ampere     = 1;

iLeftAct_Ampere    = -1;

iBottomAct_Ampere  = 1;

iRightAct_Ampere   = 1;

--iTopAct_Ampere  = tonumber(prompt("Enter current in top actuator:"));

nTurns             = 50;

mi_addcircprop('topCoil',     iTopAct_Ampere   , seriesConnection);

mi_addcircprop('leftCoil',    iLeftAct_Ampere  , seriesConnection);

mi_addcircprop('bottomCoil',  iBottomAct_Ampere, seriesConnection);

mi_addcircprop('rightCoil',   iRightAct_Ampere , seriesConnection);

-- Select block labels and apply materials

meshSize = 1000e-6;

mi_selectlabel(0.035,0)-- Air

mi_setblockprop('Air',0,meshSize,'<None>',0,0,0)

mi_clearselected();

-- Assign coil and winding direction

mi_selectlabel(avg_x_TopQuad_RightLeg_LeftSide   ,  avg_y_TopQuad_RightLeg_LeftSide ); mi_setblockprop('18 AWG',0,meshSize,'topCoil'    ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_TopQuad_RightLeg_RightSide  ,  avg_y_TopQuad_RightLeg_RightSide); mi_setblockprop('18 AWG',0,meshSize,'topCoil'    ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_TopQuad_LeftLeg_LeftSide    ,  avg_y_TopQuad_LeftLeg_LeftSide  ); mi_setblockprop('18 AWG',0,meshSize,'topCoil'    ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_TopQuad_LeftLeg_RightSide   ,  avg_y_TopQuad_LeftLeg_RightSide ); mi_setblockprop('18 AWG',0,meshSize,'topCoil'    ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_LeftQuad_TopLeg_LeftSide    ,  avg_y_LeftQuad_TopLeg_LeftSide  ); mi_setblockprop('18 AWG',0,meshSize,'leftCoil'   ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_LeftQuad_TopLeg_RightSide   ,  avg_y_LeftQuad_TopLeg_RightSide ); mi_setblockprop('18 AWG',0,meshSize,'leftCoil'   ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_LeftQuad_BotLeg_LeftSide    ,  avg_y_LeftQuad_BotLeg_LeftSide  ); mi_setblockprop('18 AWG',0,meshSize,'leftCoil'   ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_LeftQuad_BotLeg_RightSide   ,  avg_y_LeftQuad_BotLeg_RightSide ); mi_setblockprop('18 AWG',0,meshSize,'leftCoil'   ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_BotQuad_LeftLeg_LeftSide    ,  avg_y_BotQuad_LeftLeg_LeftSide  ); mi_setblockprop('18 AWG',0,meshSize,'bottomCoil' ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_BotQuad_LeftLeg_RightSide   ,  avg_y_BotQuad_LeftLeg_RightSide ); mi_setblockprop('18 AWG',0,meshSize,'bottomCoil' ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_BotQuad_RightLeg_LeftSide   ,  avg_y_BotQuad_RightLeg_LeftSide ); mi_setblockprop('18 AWG',0,meshSize,'bottomCoil' ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_BotQuad_RightLeg_RightSide  ,  avg_y_BotQuad_RightLeg_RightSide); mi_setblockprop('18 AWG',0,meshSize,'bottomCoil' ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_RightQuad_TopLeg_LeftSide   ,  avg_y_RightQuad_TopLeg_LeftSide ); mi_setblockprop('18 AWG',0,meshSize,'rightCoil'  ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(avg_x_RightQuad_TopLeg_RightSide  ,  avg_y_RightQuad_TopLeg_RightSide); mi_setblockprop('18 AWG',0,meshSize,'rightCoil'  ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_RightQuad_BotLeg_LeftSide   ,  avg_y_RightQuad_BotLeg_LeftSide ); mi_setblockprop('18 AWG',0,meshSize,'rightCoil'  ,0,0,-nTurns);  mi_clearselected();

mi_selectlabel(avg_x_RightQuad_BotLeg_RightSide  ,  avg_y_RightQuad_BotLeg_RightSide); mi_setblockprop('18 AWG',0,meshSize,'rightCoil'  ,0,0, nTurns);  mi_clearselected();

mi_selectlabel(0.09,0)-- Laminations stack

mi_setblockprop('Cold rolled low carbon strip steel',0,meshSize,'<None>',0,0,0)

mi_clearselected();

mi_selectlabel(0,0)-- Steel EN24 shaft

mi_setblockprop('M-43 Steel',0,meshSize,'<None>',0,0,0)

mi_clearselected();

-- Add boundary condition, dirichlet at the outer rim of the amb

mi_addboundprop("A=0",0,0,0,0,0,0,0,0,0);

mi_selectarcsegment(0,odOuter_m*0.5)

notHiddenInPostProcessor = 0;

boundaryGroupNum = 1;

mi_setarcsegmentprop(1,'A=0',notHiddenInPostProcessor,boundaryGroupNum);

mi_clearselected();

mi_selectarcsegment(0,-odOuter_m*0.5)

mi_setarcsegmentprop(1,'A=0',notHiddenInPostProcessor,boundaryGroupNum);

mi_clearselected();

-- Save before solving

mi_saveas(mydir .. outputFemmFilePath)

-- --Now,analyze the problem and load the solution when the analysis is finished

mi_analyze()

mi_loadsolution()

-- Outputplots :

----------------

-- Show the flux density plot

-- Parameters: legend (1 to show, 0 to hide), gscale (0 for color, 1 for grayscale), upper_B, lower_B, type (0 for B, 1 for H)

colour = 0;

showPlot = 0;

upperB_Tesla = 1.6;

lowerB_Tesla = 0;

fluxDensityPlot = 0;

fluxIntensityPlot = 1;

--mo_showvector(1) 

mo_showdensityplot(showPlot, colour, upperB_Tesla, lowerB_Tesla, fluxDensityPlot)

--mo_savebitmap("fluxDensity_plot.png")

pause()

--mo_showdensityplot(showPlot, colour, upperB_Tesla, lowerB_Tesla, fluxIntensityPlot)

--mo_savebitmap("fluxIntensity_plot.png")

pause()

-- Contour along the top 90 degree sector of the shaft

-- Left side of y axis

mo_seteditmode("contour")

mo_addcontour(-0.5*shaftOD_m*sin(45/57.3),0.5*shaftOD_m*cos(45/57.3))

mo_addcontour(-0.5*shaftOD_m*sin(35/57.3),0.5*shaftOD_m*cos(35/57.3))

mo_addcontour(-0.5*shaftOD_m*sin(25/57.3),0.5*shaftOD_m*cos(25/57.3))

mo_addcontour(-0.5*shaftOD_m*sin(15/57.3),0.5*shaftOD_m*cos(15/57.3))

mo_addcontour(-0.5*shaftOD_m*sin(05/57.3),0.5*shaftOD_m*cos(05/57.3))

-- Right side of y axis

mo_addcontour( 0.5*shaftOD_m*sin(45/57.3),0.5*shaftOD_m*cos(45/57.3))

mo_addcontour( 0.5*shaftOD_m*sin(35/57.3),0.5*shaftOD_m*cos(35/57.3))

mo_addcontour( 0.5*shaftOD_m*sin(25/57.3),0.5*shaftOD_m*cos(25/57.3))

mo_addcontour( 0.5*shaftOD_m*sin(15/57.3),0.5*shaftOD_m*cos(15/57.3))

mo_addcontour( 0.5*shaftOD_m*sin(05/57.3),0.5*shaftOD_m*cos(05/57.3))

mo_makeplot(2,200,mydir .. "top90degSect_shaft_Potential.emf")

-- --pause()

-- --quit()

 

----If wetop90degSect_shaft were interested in the flux density at specific positions, 

----we could inquire at specific points directly:

--b0=mo_getb(0.0114,0.0266);

--messagebox('Flux density at the tip of the top left  leg is ' .. b0[1] .. 'Tesla');

--b1=mo_getb(-0.0114,0.0266);

--messagebox('Flux density at the tip of the top right leg is ' .. b1[1] .. 'Tesla');

--

----The program will report the terminal properties of the circuit:

----current, voltage, and flux linkage 

--vals = mo_getcircuitproperties('topCoil');

----If we were interested in inductance, it could be obtained by

----dividing flux linkage by current

--L = 1000*vals[2]/vals[0];

--messagebox('The self-inductance of the coil is  ' .. L .. 'mHenry');