Skip to content
Snippets Groups Projects
Commit a800005a authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

No commit message

No commit message
parent b8911444
No related branches found
No related tags found
No related merge requests found
Showing
with 0 additions and 2488 deletions
Group {
// Input groups:
DefineGroup[
Domain_M = {{},
Name "Regions/0Sources/Permanent magnets"},
Domain_S = {{},
Name "Regions/0Sources/Inductor (imposed j_s)"},
Domain_Inf = {{},
Name "Regions/0Special regions/Infinite domain (spherical shell)",
Closed "1"},
Domain_Mag = {{},
Name "Regions/Other regions/Passive magnetic regions"},
Dirichlet_phi_0 = {{},
Name "Regions/0Boundary conditions/h_t = 0", Closed "1"},
Dirichlet_a_0 = {{},
Name "Regions/0Boundary conditions/b_n = 0"} ];
DefineGroup[ Domain = {{Domain_Mag, Domain_M, Domain_S, Domain_Inf},
Name "Regions/Computational domain", Visible 0} ];
}
Function{
// Input constants:
DefineConstant[ Val_Rint, Val_Rext // interior/exterior radius of Domain_Inf
];
// Input functions:
DefineFunction[ mu, // magnetic permeability
nu, // magnetic reluctivity
hc, // coercive magnetic field
js // source current density
];
// remove this: only for demo
//DefineConstant[ hcx = {0, Label "Coercive field h_x", Path "Sources"}];
//DefineConstant[ hcy = {1000, Label "Coercive field h_y", Path "Sources"}];
//hc[] = Vector[hcx,hcy,0];
//mu[] = 4*Pi*10^-7;
//nu[] = 1/mu[];
}
Jacobian {
{ Name JVol ;
Case {
{ Region Domain_Inf ; Jacobian VolSphShell{Val_Rint, Val_Rext} ; }
{ Region All ; Jacobian Vol ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
}
}
}
}
}
/* --------------------------------------------------------------------------
MagSta_phi : Magnetic scalar potential phi formulation
-------------------------------------------------------------------------- */
Constraint {
{ Name phi ;
Case {
{ Region Dirichlet_phi_0 ; Value 0. ; }
}
}
}
FunctionSpace {
{ Name Hgrad_phi ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef phin ; Function BF_Node ;
Support Domain ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef phin ; EntityType NodesOf ; NameOfConstraint phi ; }
}
}
}
Formulation {
{ Name MagSta_phi ; Type FemEquation ;
Quantity {
{ Name phi ; Type Local ; NameOfSpace Hgrad_phi ; }
}
Equation {
Galerkin { [ - mu[] * Dof{d phi} , {d phi} ] ;
In Domain ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ - mu[] * hc[] , {d phi} ] ;
In Domain_M ; Jacobian JVol ; Integration I1 ; }
}
}
}
Resolution {
{ Name MagSta_phi ;
System {
{ Name A ; NameOfFormulation MagSta_phi ; }
}
Operation {
Generate[A] ; Solve[A] ; SaveSolution[A] ;
}
}
}
PostProcessing {
{ Name MagSta_phi ; NameOfFormulation MagSta_phi ;
Quantity {
{ Name b ; Value { Local { [ - mu[] * {d phi} ] ; In Domain ; Jacobian JVol ; }
Local { [ - mu[] * hc[] ] ; In Domain_M ; Jacobian JVol ; } } }
{ Name h ; Value { Local { [ - {d phi} ] ; In Domain ; Jacobian JVol ; } } }
{ Name hc ; Value { Local { [ hc[] ] ; In Domain_M ; Jacobian JVol ; } } }
{ Name phi ; Value { Local { [ {phi} ] ; In Domain ; Jacobian JVol ; } } }
}
}
}
PostOperation {
{ Name MagSta_phi ; NameOfPostProcessing MagSta_phi;
Operation {
Print[ b, OnElementsOf Domain, File "MagSta_phi_b.pos" ] ;
Print[ h, OnElementsOf Domain, File "MagSta_phi_h.pos" ] ;
Print[ hc, OnElementsOf Domain, File "MagSta_a_hc.pos" ] ;
Print[ phi, OnElementsOf Domain, File "MagSta_phi_phi.pos" ] ;
}
}
}
/* --------------------------------------------------------------------------
MagSta_a : Magnetic vector potential a formulation (2D)
-------------------------------------------------------------------------- */
Constraint {
{ Name a ;
Case {
{ Region Dirichlet_a_0 ; Value 0. ; }
}
}
}
FunctionSpace {
{ Name Hcurl_a ; Type Form1P ;
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ;
Support Domain ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef ae ; EntityType NodesOf ; NameOfConstraint a ; }
}
}
}
Formulation {
{ Name MagSta_a ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ hc[] , {d a} ] ;
In Domain_M ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ -js[] , {a} ] ;
In Domain_S ; Jacobian JVol ; Integration I1 ; }
}
}
}
Resolution {
{ Name MagSta_a ;
System {
{ Name A ; NameOfFormulation MagSta_a ; }
}
Operation {
Generate[A] ; Solve[A] ; SaveSolution[A];
}
}
}
PostProcessing {
{ Name MagSta_a ; NameOfFormulation MagSta_a ;
Quantity {
{ Name az ; Value { Local { [ CompZ[{a}] ] ; In Domain ; Jacobian JVol ; } } }
{ Name b ; Value { Local { [ {d a} ] ; In Domain ; Jacobian JVol ; } } }
{ Name a ; Value { Local { [ {a} ] ; In Domain ; Jacobian JVol ; } } }
{ Name h ; Value { Local { [ nu[] * {d a} ] ; In Domain ; Jacobian JVol ; }
Local { [ hc[] ] ; In Domain_M ; Jacobian JVol ; } } }
{ Name hc ; Value { Local { [ hc[] ] ; In Domain_M ; Jacobian JVol ; } } }
}
}
}
PostOperation {
{ Name MagSta_a ; NameOfPostProcessing MagSta_a;
Operation {
Print[ b, OnElementsOf Domain, File "MagSta_a_b.pos" ] ;
Print[ h, OnElementsOf Domain, File "MagSta_a_h.pos" ] ;
Print[ hc, OnElementsOf Domain, File "MagSta_a_hc.pos" ] ;
Print[ a, OnElementsOf Domain, File "MagSta_a_a.pos" ] ;
}
}
}
<?xml version="1.0" encoding="utf-8"?>
<models>
<model>
<title>Magnet</title>
<summary>Simple magnet example.</summary>
<file type="pro">magnet.pro</file>
</model>
</models>
Include "magnet_data.pro";
DefineConstant[ h = {0.14, Min 0.1, Max 0.2, Step 0.01,
Name "Parameters/Geometry/Core height (m)"} ] ;
DefineConstant[ l = {0.14, Min 0.05, Max 0.2, Step 0.01,
Name "Parameters/Geometry/Core width (m)"} ] ;
DefineConstant[ d = {0.03, Min 0.01, Max 0.05, Step 0.002,
Name "Parameters/Geometry/Core thickness (m)"} ] ;
DefineConstant[ e = {5e-3, Min 5e-4, Max d, Step 1e-3,
Name "Parameters/Geometry/Air gap (m)", Highlight "LightYellow"} ] ;
DefineConstant[ ha = {0.03, Min 0.01, Max 0.1, Step 0.01,
Name "Parameters/Geometry/Magnet height (m)"} ] ;
lc0 = d / 5 ;
lc1 = e / 2 ;
lc2 = (Val_Rext - Val_Rint) / 8. ;
Point(1) = {0, 0, 0, lc0};
Point(2) = {-l/2, 0, 0, lc0};
Point(3) = {-l/2, h/2, 0, lc0};
Point(4) = {l/2, 0, 0, lc1};
Point(5) = {l/2, h/2, 0, lc0};
Point(6) = {-l/2, ha/2, 0, lc0};
Point(7) = {-l/2+d, ha/2, 0, lc0};
Point(8) = {-l/2+d, 0, 0, lc0};
Point(9) = {l/2-d, 0, 0, lc1};
Point(10) = {l/2-d, h/2-d, 0, lc0};
Point(11) = {-l/2+d, h/2-d, 0, lc0};
Point(12) = {l/2, e/2, 0, lc1};
Point(13) = {l/2-d, e/2, 0, lc1};
Point(30) = {Val_Rint, 0, 0, lc2};
Point(31) = {Val_Rext, 0, 0, lc2};
Point(32) = {0, Val_Rint, 0, lc2};
Point(33) = {0, Val_Rext, 0, lc2};
Point(34) = {-Val_Rext, 0, 0, lc2};
Point(35) = {-Val_Rint, 0, 0, lc2};
Line(1) = {34, 35};
Line(2) = {35, 2};
Line(3) = {2, 8};
Line(4) = {8, 1};
Line(5) = {1, 9};
Line(6) = {9, 4};
Line(7) = {4, 30};
Line(8) = {30, 31};
Line(9) = {2, 6};
Line(10) = {6, 3};
Line(11) = {3, 5};
Line(12) = {5, 12};
Line(13) = {12, 4};
Line(14) = {9, 13};
Line(15) = {13, 10};
Line(16) = {10, 11};
Line(17) = {11, 7};
Line(18) = {7, 8};
Line(19) = {7, 6};
Line(20) = {13, 12};
Circle(21) = {35, 1, 32};
Circle(22) = {32, 1, 30};
Circle(23) = {34, 1, 33};
Circle(24) = {33, 1, 31};
Line Loop(25) = {21, 22, 8, -24, -23, 1};
Plane Surface(26) = {25};
Line Loop(27) = - {22, -7, -13, -12, -11, -10, -9, -2, 21};
Plane Surface(28) = {27};
Line Loop(29) = - {11, 12, -20, 15, 16, 17, 19, 10};
Plane Surface(30) = {29};
Line Loop(31) = {19, -9, 3, -18};
Plane Surface(32) = {31};
Line Loop(33) = - {20, 13, -6, 14};
Plane Surface(34) = {33};
Line Loop(35) = {15, 16, 17, 18, 4, 5, 14};
Plane Surface(36) = {35};
// physical entities (for which elements will be saved)
Physical Surface(AIR) = {28, 36};
Physical Surface(AIR_INF) = {26};
Physical Surface(AIR_GAP) = {34};
Physical Surface(MAGNET) = {32};
Physical Surface(CORE) = {30};
Physical Line(LINE_INF) = {23, 24};
Physical Line(LINE_X) = {1:8};
/*
To solve the problem
with scalar potential, type 'getdp test -solve MagSta_phi -pos phi'
with vector potential, type 'getdp test -solve MagSta_a -pos a'
*/
Include "magnet_data.pro";
Group {
// AIR, AIR_INF, etc. are variables defined in core.txt, and correspond to the
// tags of physical regions in the mesh
Air = Region[ AIR ];
AirInf = Region[ AIR_INF ];
Core = Region[ CORE ];
AirGap = Region[ AIR_GAP ];
Magnet = Region[ MAGNET ];
// These are the generic group names that are used in "Magnetostatics.pro"
Domain_S = Region[ {} ] ;
Domain_Inf = Region[ AirInf ] ;
Domain_M = Region[ Magnet ] ;
Domain_Mag = Region[ {Air, Core, AirGap} ] ;
Dirichlet_a_0 = Region[ LINE_INF ] ;
Dirichlet_phi_0 = Region[ {LINE_X, LINE_INF} ] ;
}
Function {
mu0 = 4.e-7 * Pi ;
// DefineConstant is used to define a default value for murCore; this value
// can be changed interactively by the ONELAB server
DefineConstant[ murCore = {200., Min 1, Max 1000, Step 10,
Name "Parameters/Materials/Core relative permeability"} ];
nu [ Region[{Air, AirInf, AirGap, Magnet}] ] = 1. / mu0 ;
nu [ Core ] = 1. / (murCore * mu0) ;
mu [ Region[{Air, AirInf, AirGap, Magnet}] ] = mu0 ;
mu [ Core ] = murCore * mu0;
DefineConstant[ Hc = {920000,
Name "Parameters/Materials/hc", Label "Magnet coercive field (A/m)"} ];
hc [ Magnet ] = Rotate[ Vector[Hc, 0, 0.], 0, 0, Pi/2] ;
}
Include "Magnetostatics.pro"
eps = 1.e-5 ;
PostOperation {
{ Name phi ; NameOfPostProcessing MagSta_phi;
Operation {
Print[ phi, OnElementsOf Domain, File "phi.pos" ] ;
Print[ hc, OnElementsOf Domain, File "hc.pos" ] ;
Print[ b, OnElementsOf Domain, File "b_phi.pos" ] ;
Print[ b, OnLine {{-0.07,eps,0}{0.09,eps,0}} {500}, File "b_phi.txt", Format Table ] ;
}
}
{ Name a ; NameOfPostProcessing MagSta_a;
Operation {
Print[ a, OnElementsOf Domain, File "a.pos"] ;
Print[ b, OnElementsOf Domain, File "b_a.pos" ] ;
Print[ h, OnElementsOf Domain, File "h_a.pos" ] ;
Print[ b, OnLine {{-0.07,eps,0}{0.09,eps,0}} {500}, File "b_a.txt" , Format Table ] ;
}
}
}
DefineConstant[ Val_Rint = {0.15, Min 0.2, Max 1, Step 0.1,
Name "Parameters/Geometry/1Internal shell radius (m)"} ];
DefineConstant[ Val_Rext = {0.25, Min Val_Rint, Max 0.5, Step 0.1,
Name "Parameters/Geometry/2External shell radius (m)"}];
AIR = 100;
AIR_INF = 101;
AIR_GAP = 102;
MAGNET = 103;
CORE = 104;
LINE_INF = 105;
LINE_X = 106;
Function{
// nu = 100. + 10. * exp ( 1.8 * b * b )
// analytical
nu_1a[] = 100. + 10. * Exp[1.8*SquNorm[$1]] ;
dnudb2_1a[] = 18. * Exp[1.8*SquNorm[$1]] ;
h_1a[] = nu_1a[$1]*$1 ;
dhdb_1a[] = TensorDiag[1,1,1] * nu_1a[$1#1] + 2*dnudb2_1a[#1] * SquDyadicProduct[#1] ;
dhdb_1a_NL[] = 2*dnudb2_1a[$1#1] * SquDyadicProduct[#1] ;
// interpolated
Mat1_h = {
0.0000e+00, 5.5023e+00, 1.1018e+01, 1.6562e+01, 2.2149e+01, 2.7798e+01, 3.3528e+01,
3.9363e+01, 4.5335e+01, 5.1479e+01, 5.7842e+01, 6.4481e+01, 7.1470e+01, 7.8906e+01,
8.6910e+01, 9.5644e+01, 1.0532e+02, 1.1620e+02, 1.2868e+02, 1.4322e+02, 1.6050e+02,
1.8139e+02, 2.0711e+02, 2.3932e+02, 2.8028e+02, 3.3314e+02, 4.0231e+02, 4.9395e+02,
6.1678e+02, 7.8320e+02, 1.0110e+03, 1.3257e+03, 1.7645e+03, 2.3819e+03, 3.2578e+03,
4.5110e+03, 6.3187e+03, 8.9478e+03, 1.2802e+04, 1.8500e+04, 2.6989e+04, 3.9739e+04,
5.9047e+04, 8.8520e+04, 1.3388e+05, 2.0425e+05, 3.1434e+05, 4.8796e+05, 7.6403e+05
} ;
Mat1_b = {
0.0000e+00, 5.0000e-02, 1.0000e-01, 1.5000e-01, 2.0000e-01, 2.5000e-01, 3.0000e-01,
3.5000e-01, 4.0000e-01, 4.5000e-01, 5.0000e-01, 5.5000e-01, 6.0000e-01, 6.5000e-01,
7.0000e-01, 7.5000e-01, 8.0000e-01, 8.5000e-01, 9.0000e-01, 9.5000e-01, 1.0000e+00,
1.0500e+00, 1.1000e+00, 1.1500e+00, 1.2000e+00, 1.2500e+00, 1.3000e+00, 1.3500e+00,
1.4000e+00, 1.4500e+00, 1.5000e+00, 1.5500e+00, 1.6000e+00, 1.6500e+00, 1.7000e+00,
1.7500e+00, 1.8000e+00, 1.8500e+00, 1.9000e+00, 1.9500e+00, 2.0000e+00, 2.0500e+00,
2.1000e+00, 2.1500e+00, 2.2000e+00, 2.2500e+00, 2.3000e+00, 2.3500e+00, 2.4000e+00
} ;
Mat1_b2 = List[Mat1_b]^2;
Mat1_nu = List[Mat1_h]/List[Mat1_b];
Mat1_nu(0) = Mat1_nu(1);
Mat1_nu_b2 = ListAlt[Mat1_b2, Mat1_nu] ;
nu_1[] = InterpolationLinear[ SquNorm[$1] ]{List[Mat1_nu_b2]} ;
dnudb2_1[] = dInterpolationLinear[SquNorm[$1]]{List[Mat1_nu_b2]} ;
h_1[] = nu_1[$1] * $1 ;
dhdb_1[] = TensorDiag[1,1,1] * nu_1[$1#1] + 2*dnudb2_1[#1] * SquDyadicProduct[#1] ;
dhdb_1_NL[] = 2*dnudb2_1[$1#1] * SquDyadicProduct[#1] ;
// nu = 123. + 0.0596 * exp ( 3.504 * b * b )
// analytical 3kW machine
nu_3kWa[] = 123. + 0.0596 * Exp[3.504*SquNorm[$1]] ;
dnudb2_3kWa[] = 0.0596*3.504 * Exp[3.504*SquNorm[$1]] ;
h_3kWa[] = nu_3kWa[$1]*$1 ;
dhdb_3kWa[] = TensorDiag[1,1,1] * nu_3kWa[$1#1] + 2*dnudb2_3kWa[#1] * SquDyadicProduct[$1] ;
dhdb_3kWa_NL[] = 2*dnudb2_3kWa[$1#1] * SquDyadicProduct[#1] ;
// interpolated
Mat3kW_h = {
0.0000e+00, 6.1465e+00, 1.2293e+01, 1.8440e+01, 2.4588e+01, 3.0736e+01, 3.6886e+01,
4.3037e+01, 4.9190e+01, 5.5346e+01, 6.1507e+01, 6.7673e+01, 7.3848e+01, 8.0036e+01,
8.6241e+01, 9.2473e+01, 9.8745e+01, 1.0508e+02, 1.1150e+02, 1.1806e+02, 1.2485e+02,
1.3199e+02, 1.3971e+02, 1.4836e+02, 1.5856e+02, 1.7137e+02, 1.8864e+02, 2.1363e+02,
2.5219e+02, 3.1498e+02, 4.2161e+02, 6.0888e+02, 9.4665e+02, 1.5697e+03, 2.7417e+03,
4.9870e+03, 9.3633e+03, 1.8037e+04, 3.5518e+04, 7.1329e+04, 1.4591e+05, 3.0380e+05,
6.4363e+05, 1.3872e+06, 3.0413e+06, 6.7826e+06, 1.5386e+07, 3.5504e+07, 8.3338e+07
} ;
Mat3kW_b = {
0.0000e+00, 5.0000e-02, 1.0000e-01, 1.5000e-01, 2.0000e-01, 2.5000e-01, 3.0000e-01,
3.5000e-01, 4.0000e-01, 4.5000e-01, 5.0000e-01, 5.5000e-01, 6.0000e-01, 6.5000e-01,
7.0000e-01, 7.5000e-01, 8.0000e-01, 8.5000e-01, 9.0000e-01, 9.5000e-01, 1.0000e+00,
1.0500e+00, 1.1000e+00, 1.1500e+00, 1.2000e+00, 1.2500e+00, 1.3000e+00, 1.3500e+00,
1.4000e+00, 1.4500e+00, 1.5000e+00, 1.5500e+00, 1.6000e+00, 1.6500e+00, 1.7000e+00,
1.7500e+00, 1.8000e+00, 1.8500e+00, 1.9000e+00, 1.9500e+00, 2.0000e+00, 2.0500e+00,
2.1000e+00, 2.1500e+00, 2.2000e+00, 2.2500e+00, 2.3000e+00, 2.3500e+00, 2.4000e+00
} ;
Mat3kW_b2 = List[Mat3kW_b]^2;
Mat3kW_nu = List[Mat3kW_h]/List[Mat3kW_b];
Mat3kW_nu(0) = Mat3kW_nu(1);
Mat3kW_nu_b2 = ListAlt[Mat3kW_b2, Mat3kW_nu] ;
nu_3kW[] = InterpolationLinear[SquNorm[$1]]{List[Mat3kW_nu_b2]} ;
dnudb2_3kW[] = dInterpolationLinear[SquNorm[$1]]{List[Mat3kW_nu_b2]} ;
h_3kW[] = nu_3kW[$1] * $1 ;
dhdb_3kW[] = TensorDiag[1,1,1]*nu_3kW[$1#1] + 2*dnudb2_3kW[#1] * SquDyadicProduct[#1] ;
dhdb_3kW_NL[] = 2*dnudb2_3kW[$1] * SquDyadicProduct[$1] ;
}
<?xml version="1.0" encoding="utf-8"?>
<models>
<model>
<title>Eight-pole permanent magnet synchronous machine</title>
<summary></summary>
<file type="geo">pmsm.geo</file>
<preview type="png">screenshot1_512.png</preview>
<url>http://onelab.info/wiki/Electric_Machines</url>
</model>
</models>
This diff is collapsed.
Include "pmsm_data.geo";
Mesh.Algorithm = 6; // 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Geometry.CopyMeshingMethod = 1;
//Mesh.CharacteristicLengthFactor = 0.5 ;
fact_trans = Mesh.CharacteristicLengthFactor ;
// Mesh characteristic lengths
s = 0.4 ;
pR1=(rR2-rR1)/6.*s;
pR2=(rR2-rR1)/6.*s;
pS1=(rS7-rS1)/7.*s;
pS2=(rS7-rS1)/12.*s;
pS3=(rS6-rS3)/10.*s;
NbrDivMB = 2*Ceil[2*Pi*rRext/8/pR1/fact_trans] ; //1/8 Moving band
//--------------------------------------------------------------------------------
//--------------------------------------------------------------------------------
cen = newp ; Point(cen)={0,0,0,pR1};
nicepos_rotor[] = {};
nicepos_stator[] = {};
Include "pmsm_rotor.geo";
Include "pmsm_stator.geo";
// For nice visualisation...
//Mesh.Light = 0 ;
//Mesh.SurfaceFaces = 1; Mesh.SurfaceEdges=0;
Hide { Point{ Point '*' }; }
Hide { Line{ Line '*' }; }
Show { Line{ nicepos_rotor[], nicepos_stator[] }; }
Physical Line(NICEPOS) = { nicepos_rotor[], nicepos_stator[] };
//For post-processing...
//View[0].Light = 0;
View[0].NbIso = 25; // Number of intervals
View[0].IntervalsType = 1;
DefineConstant[ Flag_AddInfo = {0, Choices{0,1},
Name "Input/02Add info about phases and axis"} ];
For i In {PostProcessing.NbViews-1 : 0 : -1}
If(StrFind(View[i].Attributes, "tmp"))
Delete View[i];
EndIf
EndFor
If(Flag_AddInfo)
rr = 1.25 * rS3 ;
For k In {0:NbrPolesInModel-1}
xa[] += rr*Cos(1*Pi/24+k*Pi/4) ; ya[] += rr*Sin(1*Pi/24+k*Pi/4) ;
xb[] += rr*Cos(3*Pi/24+k*Pi/4) ; yb[] += rr*Sin(3*Pi/24+k*Pi/4) ;
xc[] += rr*Cos(5*Pi/24+k*Pi/4) ; yc[] += rr*Sin(5*Pi/24+k*Pi/4) ;
EndFor
// Adding some axes
rr0 = 0.3 * rS7 ;
rr1 = 1.3 * rS7 ;
th_d = InitialRotorAngle ;
th_q = th_d + 22.5 * deg2rad ;
th_a = 30 * deg2rad ;
th_b = (30 + 120/4) * deg2rad ;
th_c = (30 + 240/4) * deg2rad ;
ff = 0.9;
xd[0] = rr0*Cos(th_d) ; yd[0] = rr0*Sin(th_d) ;
xd[1] = ff*rr1*Cos(th_d) ; yd[1] = ff*rr1*Sin(th_d) ;
xq[0] = rr0*Cos(th_q) ; yq[0] = rr0*Sin(th_q) ;
xq[1] = ff*rr1*Cos(th_q) ; yq[1] = ff*rr1*Sin(th_q) ;
xaa[0] = rr0*Cos(th_a) ; yaa[0] = rr0*Sin(th_a) ;
xaa[1] = rr1*Cos(th_a) ; yaa[1] = rr1*Sin(th_a) ;
xbb[0] = rr0*Cos(th_b) ; ybb[0] = rr0*Sin(th_b) ;
xbb[1] = rr1*Cos(th_b) ; ybb[1] = rr1*Sin(th_b) ;
xcc[0] = rr0*Cos(th_c) ; ycc[0] = rr0*Sin(th_c) ;
xcc[1] = rr1*Cos(th_c) ; ycc[1] = rr1*Sin(th_c) ;
Include "info_view.geo";
EndIf
//
// Permanent Magnet Synchronous Generator
//
Include "pmsm_data.geo";
DefineConstant[
Flag_AnalysisType = {1, Choices{0="Static", 1="Time domain"},
Name "Input/19Type of analysis", Highlight "Blue",
Help Str["- Use 'Static' to compute static fields created in the machine",
"- Use 'Time domain' to compute the dynamic response of the machine"]} ,
Flag_SrcType_Stator = { 0, Choices{0="None",1="Current"},
Name "Input/41Source type in stator", Highlight "Blue"},
Flag_NL = { 1, Choices{0,1},
Name "Input/60Nonlinear BH-curve"},
Flag_NL_law_Type = { 0, Choices{
0="Analytical", 1="Interpolated",
2="Analytical VH800-65D", 3="Interpolated VH800-65D"},
Name "Input/61BH-curve", Highlight "Blue", Visible Flag_NL}
];
Flag_Cir = !Flag_SrcType_Stator ;
Group {
Stator_Fe = #STATOR_FE ;
Stator_Al = #{};
Stator_Cu = #{};
Stator_Air = #STATOR_AIR ;
Stator_Airgap = #STATOR_AIRGAP ;
Stator_Bnd_A0 = #STATOR_BND_A0 ;
Stator_Bnd_A1 = #STATOR_BND_A1 ;
Rotor_Fe = #ROTOR_FE ;
Rotor_Al = #{};
Rotor_Cu = #{};
Rotor_Fe = #ROTOR_FE ;
Rotor_Air = #ROTOR_AIR ;
Rotor_Airgap = #ROTOR_AIRGAP ;
Rotor_Bnd_A0 = #ROTOR_BND_A0 ;
Rotor_Bnd_A1 = #ROTOR_BND_A1 ;
MovingBand_PhysicalNb = #MOVING_BAND ; // Fictitious number for moving band, not in the geo file
Surf_Inf = #SURF_EXT ;
Surf_bn0 = #SURF_INT ;
Surf_cutA0 = #{STATOR_BND_A0, ROTOR_BND_A0};
Surf_cutA1 = #{STATOR_BND_A1, ROTOR_BND_A1};
Dummy = #NICEPOS;
nbMagnets = NbrPolesTot/SymmetryFactor ;
For k In {1:nbMagnets}
Rotor_Magnet~{k} = Region[ (ROTOR_MAGNET+k-1) ];
Rotor_Magnets += Region[ Rotor_Magnet~{k} ];
EndFor
nbInds = (Flag_Symmetry) ? NbrPolesInModel*NbrSectTotStator/NbrPolesTot : NbrSectTotStator ;
Printf("NbrPolesInModel=%g, nbInds=%g SymmetryFactor=%g",
NbrPolesInModel, nbInds, SymmetryFactor);
Stator_Ind_Ap = #{}; Stator_Ind_Am = #{STATOR_IND_AM};
Stator_Ind_Bp = #{}; Stator_Ind_Bm = #{STATOR_IND_BM};
Stator_Ind_Cp = #{STATOR_IND_CP}; Stator_Ind_Cm = #{};
If(NbrPolesInModel > 1)
Stator_Ind_Ap += #STATOR_IND_AP;
Stator_Ind_Bp += #STATOR_IND_BP;
Stator_Ind_Cm += #STATOR_IND_CM;
EndIf
PhaseA = Region[{ Stator_Ind_Ap, Stator_Ind_Am }];
PhaseB = Region[{ Stator_Ind_Bp, Stator_Ind_Bm }];
PhaseC = Region[{ Stator_Ind_Cp, Stator_Ind_Cm }];
// FIXME: Just one physical region for nice graph in Onelab
PhaseA_pos = Region[{ Stator_Ind_Am }];
PhaseB_pos = Region[{ Stator_Ind_Bm }];
PhaseC_pos = Region[{ Stator_Ind_Cp }];
Stator_IndsP = Region[{ Stator_Ind_Ap, Stator_Ind_Bp, Stator_Ind_Cp }];
Stator_IndsN = Region[{ Stator_Ind_Am, Stator_Ind_Bm, Stator_Ind_Cm }];
Stator_Inds = Region[ {PhaseA, PhaseB, PhaseC} ] ;
Rotor_Inds = Region[ {} ] ;
StatorC = Region[{ }] ;
StatorCC = Region[{ Stator_Fe }] ;
RotorC = Region[{ }] ;
RotorCC = Region[{ Rotor_Fe, Rotor_Magnets }] ;
// Moving band: with or without symmetry, these BND lines must be complete
Stator_Bnd_MB = #STATOR_BND_MOVING_BAND;
For k In {1:SymmetryFactor}
Rotor_Bnd_MB~{k} = Region[ (ROTOR_BND_MOVING_BAND+k-1) ];
Rotor_Bnd_MB += Region[ Rotor_Bnd_MB~{k} ];
EndFor
Rotor_Bnd_MBaux = Region[ {Rotor_Bnd_MB, -Rotor_Bnd_MB~{1}}];
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
Function {
NbrPhases = 3 ;
// For a radial remanent b
For k In {1:nbMagnets}
br[ Rotor_Magnet~{k} ] = (-1)^(k-1) * b_remanent * Vector[ Cos[Atan2[Y[],X[]]], Sin[Atan2[Y[],X[]]], 0 ];
EndFor
//Data for modeling a stranded inductor
NbWires[] = 104 ; // Number of wires per slot
// STATOR_IND_AM comprises all the slots in that phase, we need thus to divide by the number of slots
nbSlots[] = Ceil[nbInds/NbrPhases/2] ;
SurfCoil[] = SurfaceArea[]{STATOR_IND_AM}/nbSlots[] ;//All inductors have the same surface
FillFactor_Winding = 0.5 ; // percentage of Cu in the surface coil side, smaller than 1
Factor_R_3DEffects = 1.5 ; // bigger than Adding 50% of resistance
DefineConstant[ rpm = { rpm_nominal,
Name "Input/7speed in rpm",
Highlight "AliceBlue", Visible (Flag_AnalysisType==1)} ]; // speed in rpm
wr = rpm/60*2*Pi ; // speed in rad_mec/s
// supply at fixed position
DefineConstant[ Freq = { wr*NbrPolePairs/(2*Pi), ReadOnly 1,
Name "Output/1Freq", Highlight "LightYellow" } ];
Omega = 2*Pi*Freq ;
T = 1/Freq ;
DefineConstant[
thetaMax_deg = { 180, Name "Input/21End rotor angle (loop)",
Highlight "AliceBlue", Visible (Flag_AnalysisType==1) }
];
theta0 = InitialRotorAngle + 0. ;
thetaMax = thetaMax_deg * deg2rad ; // end rotor angle (used in doing a loop)
DefineConstant[
NbTurns = { (thetaMax-theta0)/(2*Pi), Name "Input/24Number of revolutions",
Highlight "LightGrey", ReadOnly 1, Visible (Flag_AnalysisType==1)},
delta_theta_deg = { 1., Name "Input/22Step [deg]",
Highlight "AliceBlue", Visible (Flag_AnalysisType==1)}
];
delta_theta[] = delta_theta_deg * deg2rad ;
time0 = 0 ; // at initial rotor position
delta_time = delta_theta_deg * deg2rad/wr;
timemax = thetaMax/wr;
DefineConstant[
NbSteps = { Ceil[(timemax-time0)/delta_time], Name "Input/23Number of steps",
Highlight "LightGrey", ReadOnly 1, Visible (Flag_AnalysisType==1)}
];
RotorPosition[] = InitialRotorAngle + $Time * wr ;
RotorPosition_deg[] = RotorPosition[]*180/Pi;
Flag_ParkTransformation = 1 ;
Theta_Park[] = ((RotorPosition[] + Pi/8) - Pi/6) * NbrPolePairs; // electrical degrees
Theta_Park_deg[] = Theta_Park[]*180/Pi;
DefineConstant[
ID = { 0, Name "Input/50Id stator current",
Highlight "AliceBlue", Visible (Flag_SrcType_Stator==1)},
IQ = { Inominal, Name "Input/51Iq stator current",
Highlight "AliceBlue", Visible (Flag_SrcType_Stator==1)}
] ;
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
If(Flag_SrcType_Stator)
UndefineConstant["Input/8Load resistance"];
EndIf
If(Flag_Cir)
Include "pmsm_8p_circuit.pro" ;
EndIf
Include "machine_magstadyn_a.pro" ;
//
// Circuit for Permanent Magnet Synchronous Generator
//
Group{
R1 = #55551 ;
R2 = #55552 ;
R3 = #55553 ;
Input1 = #10001 ;
Input2 = #10002 ;
Input3 = #10003 ;
Input4 = #10004 ;
Resistance_Cir = Region[{R1, R2, R3}];
DomainZ_Cir = Region[ {Resistance_Cir} ];
DomainSource_Cir = Region[ {} ] ;
If(Flag_SrcType_Stator>1)
DomainSource_Cir += Region[ {Input1, Input2, Input3} ] ;
EndIf
DomainZt_Cir = Region[ {DomainZ_Cir, DomainSource_Cir} ];
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
Function {
// Open circuit - load - short circuit
DefineConstant[ ZR = {200, Choices{1e-8, 200, 1e8},
Name "Input/8Load resistance", Highlight "AliceBlue"} ];
Resistance[#{R1, R2, R3}] = ZR ;
}
// --------------------------------------------------------------------------
Constraint {
If (SymmetryFactor<8)
If(Flag_SrcType_Stator==0)
{ Name ElectricalCircuit ; Type Network ;
Case Circuit1 {
{ Region Stator_Ind_Ap ; Branch {100,102} ; }
{ Region Stator_Ind_Am ; Branch {103,102} ; }
{ Region R1 ; Branch {103,100} ; }
}
Case Circuit2 {
{ Region Stator_Ind_Bp ; Branch {200,202} ; }
{ Region Stator_Ind_Bm ; Branch {203,202} ; }
{ Region R2 ; Branch {203,200} ; }
}
Case Circuit3 {
{ Region Stator_Ind_Cp ; Branch {300,302} ; }
{ Region Stator_Ind_Cm ; Branch {303,302} ; }
{ Region R3 ; Branch {303,300} ; }
}
}
EndIf
If (Flag_SrcType_Stator==2)
{ Name ElectricalCircuit ; Type Network ;
Case Circuit1 {
{ Region Input1 ; Branch {100,101} ; }
{ Region Stator_Ind_Ap ; Branch {101,102} ; }
{ Region Stator_Ind_Am ; Branch {103,102} ; }
{ Region R1 ; Branch {103,100} ; }
}
Case Circuit2 {
{ Region Input2 ; Branch {200,201} ; }
{ Region Stator_Ind_Bp ; Branch {201,202} ; }
{ Region Stator_Ind_Bm ; Branch {203,202} ; }
{ Region R2 ; Branch {203,200} ; }
}
Case Circuit3 {
{ Region Input3 ; Branch {300,301} ; }
{ Region Stator_Ind_Cp ; Branch {301,302} ; }
{ Region Stator_Ind_Cm ; Branch {303,302} ; }
{ Region R3 ; Branch {303,300} ; }
}
}
EndIf
EndIf
If(SymmetryFactor==8)
If(Flag_SrcType_Stator==0) // Only one physical region in geo allow per branch
{ Name ElectricalCircuit ; Type Network ;
Case Circuit1 {
{ Region PhaseA ; Branch {100,102} ; }
{ Region R1 ; Branch {102,100} ; }
}
Case Circuit2 {
{ Region PhaseB ; Branch {200,202} ; }
{ Region R2 ; Branch {202,200} ; }
}
Case Circuit3 {
{ Region PhaseC ; Branch {300,302} ; }
{ Region R3 ; Branch {302,300} ; }
}
}
EndIf
If(Flag_SrcType_Stator==2) // Only one physical region in geo allow per branch
{ Name ElectricalCircuit ; Type Network ;
Case Circuit1 {
{ Region Input1 ; Branch {100,101} ; }
{ Region PhaseA ; Branch {101,102} ; }
{ Region R1 ; Branch {102,100} ; }
}
Case Circuit2 {
{ Region Input2 ; Branch {200,201} ; }
{ Region PhaseB ; Branch {201,202} ; }
{ Region R2 ; Branch {202,200} ; }
}
Case Circuit3 {
{ Region Input3 ; Branch {300,301} ; }
{ Region PhaseC ; Branch {302,301} ; }
{ Region R3 ; Branch {302,300} ; }
}
}
EndIf
EndIf
}
// Permanent magnet synchronous machine
// Example of Prof. Dr. Mauricio Valencia Ferreira da Luz (Florianopolis, August 23, 2010)
// Modified and customised for Onelab by Ruth V. Sabariego (February, 2013)
mm = 1e-3 ;
deg2rad = Pi/180 ;
pp = "Input/Constructive parameters/";
DefineConstant[
NbrPolesInModel = { 1, Choices {1="1", 2="2", 4="4", 8="8"},
Name "Input/20Number of poles in FE model",
Highlight "Blue", Visible 1},
InitialRotorAngle_deg = {7.5,
Name "Input/21Start rotor angle [deg]",
Highlight "AliceBlue"}
] ;
//--------------------------------------------------------------------------------
InitialRotorAngle = InitialRotorAngle_deg*deg2rad ; // initial rotor angle, 0 if aligned
//------------------------------------------------
//------------------------------------------------
NbrPolesTot = 8 ; // number of poles in complete cross-section
NbrPolePairs = NbrPolesTot/2 ;
SymmetryFactor = NbrPolesTot/NbrPolesInModel ;
Flag_Symmetry = (SymmetryFactor==1)?0:1 ;
NbrSectTot = NbrPolesTot ; // number of "rotor teeth"
NbrSect = NbrSectTot*NbrPolesInModel/NbrPolesTot ; // number of "rotor teeth" in FE model
//--------------------------------------------------------------------------------
//------------------------------------------------
// Stator
//------------------------------------------------
NbrSectTotStator = 24; // number of stator teeth
NbrSectStator = NbrSectTotStator*NbrPolesInModel/NbrPolesTot; // number of stator teeth in FE model
//--------------------------------------------------------------------------------
DefineConstant[
lm = {2.352*mm , Name StrCat[pp, "Magnet height [m]"], Closed 0}
];
Th_magnet = 32.67 *deg2rad ; // angle in degrees 0 < Th_magnet < 45
//--------------------------------------------------------------------------------
rRext = 25.6*mm;
rR1 = 10.5*mm;
rR2 = (rRext-lm); //23.243e-03;
rR3 = (rRext-0.7389*lm); //23.862e-03;
rR4 = (rRext-0.72278*lm); //23.9e-03;
rR5 = rRext; //25.6e-03;
//Gap = rS1-rR5;
DefineConstant[
AxialLength = {35*mm, Name StrCat[pp, "Axial length [m]"], Closed 1},
Gap = {(26.02-25.6)*mm, Name StrCat[pp, "Airgap width [m]"], Closed 0}
];
rS1 = rR5 + Gap; //rS1 = 26.02*mm;
rS2 = rS1 + 0.6*mm; //rS2 = 26.62*mm;
rS3 = rS2 + 0.34*mm; //rS3 = 26.96*mm;
rS4 = rS3 + 11.2*mm; //rS4 = 38.16*mm;
rS5 = rS4 + 0.11*mm; //rS5 = 38.27*mm;
rS6 = rS5 + 1.75*mm; //rS6 = 40.02*mm;
rS7 = rS6 + 5.98*mm; //rS7 = 46.00*mm;
rB1 = rR5+Gap/3;
rB1b = rB1;
rB2 = rR5+Gap*2/3;
A0 = 45 * deg2rad ; // with this choice, axis A of stator is at 30 degrees with regard to horizontal axis
A1 = 0 * deg2rad ; // Rotor initial aligned position, current position in angRot
sigma_fe = 0. ; // laminated steel
DefineConstant[
mur_fe = {1000, Name StrCat[pp, "Relative permeability for linear case"]},
b_remanent = {1.2, Name StrCat[pp, "Remanent induction [T]"] }
];
rpm_nominal = 500 ;
Inominal = 3.9 ; // Nominal current
Tnominal = 2.5 ; // Nominal torque
// ----------------------------------------------------
// Numbers for physical regions in .geo and .pro files
// ----------------------------------------------------
// Rotor
ROTOR_FE = 1000 ;
ROTOR_AIR = 1001 ;
ROTOR_AIRGAP = 1002 ;
ROTOR_MAGNET = 1010 ; // Index for first Magnet (1/8 model->1; full model->8)
ROTOR_BND_MOVING_BAND = 1100 ; // Index for first line (1/8 model->1; full model->8)
ROTOR_BND_A0 = 1200 ;
ROTOR_BND_A1 = 1201 ;
SURF_INT = 1202 ;
// Stator
STATOR_FE = 2000 ;
STATOR_AIR = 2001 ;
STATOR_AIRGAP = 2002 ;
STATOR_BND_MOVING_BAND = 2100 ;// Index for first line (1/8 model->1; full model->8)
STATOR_BND_A0 = 2200 ;
STATOR_BND_A1 = 2201 ;
STATOR_IND = 2300 ; //Index for first Ind (1/8 model->3; full model->24)
STATOR_IND_AP = STATOR_IND + 1 ; STATOR_IND_BM = STATOR_IND + 2 ;STATOR_IND_CP = STATOR_IND + 3 ;
STATOR_IND_AM = STATOR_IND + 4 ; STATOR_IND_BP = STATOR_IND + 5 ;STATOR_IND_CM = STATOR_IND + 6 ;
SURF_EXT = 3000 ; // outer boundary
MOVING_BAND = 9999 ;
NICEPOS = 111111 ;
//--------------------------------------------------------------------------------
// Rotor PMSM
//--------------------------------------------------------------------------------
A = InitialRotorAngle-45/2*deg2rad + A1; // with Theta_Park
sinA = Sin(A); cosA = Cos(A);
pntR[]+=newp; Point(newp)={rR1*cosA, rR1*sinA, 0, pR1};
pntR[]+=newp; Point(newp)={rR2*cosA, rR2*sinA, 0, pR1};
pntR[]+=newp; Point(newp)={rR4*cosA, rR4*sinA, 0, pR1};
pntR[]+=newp; Point(newp)={rR5*cosA, rR5*sinA, 0, pR1};
pntR[]+=newp; Point(newp)={rB1*cosA, rB1*sinA, 0, pR2};
For k In {0:#pntR[]-2}
linR0[]+=newl; Line(newl) = {pntR[k], pntR[k+1]};
EndFor
Transfinite Line{linR0[0]} = Ceil[(rR2-rR1)/pR1/fact_trans] ;
Transfinite Line{linR0[1]} = Ceil[(rR4-rR2)/pR1/fact_trans] ;
Transfinite Line{linR0[2]} = Ceil[(rR5-rR4)/pR1/fact_trans] ;
Transfinite Line{linR0[3]} = Ceil[(rB1-rR5)/pR1/fact_trans] ;
For k In {0:#linR0[]-1}
linR1[] += Rotate {{0, 0, 1}, {0, 0, 0}, A0+A1} { Duplicata{Line{linR0[k]};} };
EndFor
AA[] = {(A0-Th_magnet)/2+A1, Th_magnet, (A0-Th_magnet)/2+A1} ;
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[0]} { Point{pntR[0]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[1]} { Point{lin[0]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[2]} { Point{lin[0]}; };
cirR[]+=lin[1];
surfint[]=cirR[{0,1,2}] ; // boundary conditions
pMagnet[] = Rotate {{0, 0, 1}, {0, 0, 0}, AA[0]} { Duplicata{Point{pntR[1]};} };
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[1]} { Point{pMagnet[0]}; };
pMagnet[] += lin[0];
cirR[] += lin[1] ;
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[0]} { Point{pntR[2]}; };
cirR[]+=lin[1]; pMagnet[] += lin[0];
pMagnet[] += Rotate {{0, 0, 1}, {0, 0, 0}, AA[1]} { Duplicata{Point{lin[0]};} };
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[2]} { Point{pMagnet[3]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[0]} { Point{pntR[3]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[1]} { Point{lin[0]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, AA[2]} { Point{lin[0]}; };
cirR[]+=lin[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, A0+A1} { Point{pntR[4]}; };
cirR[]+=lin[1];
linR2[] = Rotate {{0, 0, 1}, {0, 0, 0}, (A0-Th_magnet)/2+A1} { Duplicata{Line{linR0[{1,2}]};} };
linR3[] = Rotate {{0, 0, 1}, {0, 0, 0},-(A0-Th_magnet)/2+A1} { Duplicata{Line{linR1[{1,2}]};} };
// surfaces rotor
Line Loop(newll) = {linR0[{0,1}], cirR[4], -linR2[0], cirR[3], linR3[0], cirR[5], -linR1[{1,0}], -cirR[{2,1,0}]};
srotor[0]=news; Plane Surface(srotor[0]) = {newll-1};
Line Loop(newl) = {linR2[1], cirR[7], -linR3[{1,0}], -cirR[3], linR2[0]};
smagnet[0]=news; Plane Surface(smagnet[0]) = {newll-1};
nn = #cirR[]-1 ;
Line Loop(newll) = {cirR[{nn-5}], linR2[1], -cirR[{nn-3}], -linR0[2]};
sairrotor[]+=news; Plane Surface(news) = -{newll-1};
Line Loop(newll) = {cirR[{nn-4}], linR1[2], -cirR[{nn-1}], -linR3[1]};
sairrotor[]+=news; Plane Surface(news) = -{newll-1};
Line Loop(newll) = {linR0[3], cirR[nn], -linR1[3], -cirR[{nn-1:nn-3:-1}]};
sairrotormb[]+=news; Plane Surface(news) = {newll-1};
// -------------------------------------------------------------------------------
// Moving band == AirGap rotor side
// -------------------------------------------------------------------------------
Transfinite Line{cirR[nn]} = NbrDivMB+1 ;
//Filling the gap for the whole 2*Pi
lineMBrotor[]=cirR[{nn}];
For k In {1:NbrPolesTot-1}
lineMBrotoraux[]+=Rotate {{0, 0, 1}, {0, 0, 0}, k*A0} { Duplicata{Line{lineMBrotor[]};} };
EndFor
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------------------
If(SymmetryFactor<8)
// FULL MODEL ==> Rotation of NbrPolesTot*Pi/4
// For simplicity: rotating first the interior and exterior boundaries
If (SymmetryFactor>1)
For k In {0:#linR1[]-1}
linR1_[] += Rotate {{0, 0, 1}, {0, 0, 0}, 2*Pi/SymmetryFactor-Pi/4} { Duplicata{Line{linR1[k]};} };
EndFor
linR1[] = linR1_[];
EndIf
For k In {1:NbrPolesInModel-1}
surfint[] += Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Line{surfint[{0:2}]};} };
EndFor
For k In {1:NbrPolesInModel-1}
srotor[] += Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{srotor[0]};} };
smagnet[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{smagnet[0]};} };
sairrotor[] += Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sairrotor[{0,1}]};} };
sairrotormb[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sairrotormb[0]};} };
EndFor
EndIf
// -------------------------------------------------------------------------------
// Physical regions
// -------------------------------------------------------------------------------
Physical Surface(ROTOR_FE) = {srotor[]}; // Rotor
Physical Surface(ROTOR_AIR) = {sairrotor[]}; // AirRotor
Physical Surface(ROTOR_AIRGAP) = {sairrotormb[]};// AirRotor for possible torque computation with Maxwell stress tensor
NN = (Flag_Symmetry)?NbrPolesInModel:NbrPolesTot;
For k In {0:NN-1}
Physical Surface(ROTOR_MAGNET+k) = {smagnet[k]}; // Magnets
EndFor
Physical Line(SURF_INT) = {surfint[]}; // SurfInt
If(Flag_Symmetry) //Lines for symmetry link
Physical Line(ROTOR_BND_A0) = linR0[];
Physical Line(ROTOR_BND_A1) = linR1[];
EndIf
lineMBrotor[] += lineMBrotoraux[] ;
If(!Flag_Symmetry)
Physical Line(ROTOR_BND_MOVING_BAND) = {lineMBrotor[]};
EndIf
If(Flag_Symmetry)
nr = #lineMBrotor[];
nnp = nr/(NbrPolesTot/NbrSect) ;
For k In {1:Floor[NbrPolesTot/NbrSect]}
kk= ((k*nnp-1) > nr) ? nr-1 : k*nnp-1 ;
Physical Line(ROTOR_BND_MOVING_BAND+k-1) = lineMBrotor[{(k-1)*nnp:kk}] ;
EndFor
k1 = Floor[NbrPolesTot/NbrSect];
k2 = Ceil[NbrPolesTot/NbrSect];
If (k2 > k1)
Physical Line(ROTOR_BND_MOVING_BAND+k2-1) = lineMBrotor[{(k2-1)*nnp:#lineMBrotor[]-1}] ;
EndIf
EndIf
// For nice visualisation...
linRotor[] = CombinedBoundary{Surface{srotor[]};};
linMagnet[] = Boundary{Surface{smagnet[]};};
nicepos_rotor[] += { linRotor[], linMagnet[] };
Color SteelBlue {Surface{srotor[]};}
Color SkyBlue {Surface{sairrotor[], sairrotormb[]};}
Color Orchid {Surface{smagnet[{0:#smagnet[]-1:2}]};}
If(#smagnet[]>1)
Color Purple {Surface{smagnet[{1:#smagnet[]-1:2}]};}
EndIf
// -------------------------------------------------------------------------------
// Moving band == AirGap stator side
// -------------------------------------------------------------------------------
pntG[]+=newp; Point(newp) = {rB2, 0., 0., pS1}; // aligned with the stator
circ[] = Extrude {{0, 0, 1}, {0, 0, 0}, A0} { Point{pntG[0]}; };
pntG[]+=circ[0];
lineMBstator[]=circ[1];
Transfinite Line{lineMBstator[0]} = NbrDivMB+1 ;
//Filling the gap for the whole 2*Pi
For k In {1:NbrPolesTot-1}
lineMBstatoraux[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*A0} { Duplicata{Line{lineMBstator[0]};} };
EndFor
// -------------------------------------------------------------------------------
// Stator
// -------------------------------------------------------------------------------
pntS[] = newp; Point(newp)={rS1, 0, 0, pS1};
linS[] = newl; Line(newl) = {pntG[0], pntS[0]};
linS[]+= Rotate {{0, 0, 1}, {0, 0, 0}, A0} { Duplicata{Line{linS[0]};} };
pntS[]+=newp; Point(newp)={rS7,0,0,pS2};
points[]=Boundary{Line{linS[1]};};
pntS[]+=points[1];
lin[] = Extrude {{0, 0, 1}, {0, 0, 0}, A0} { Point{pntS[1]}; };
cirS[]= lin[1]; pntS[]+=lin[0];
linS[]+=newl; Line(newl) = {pntS[0], pntS[1]};
linS[]+=newl; Line(newl) = {pntS[2], pntS[3]};
// -------------------------------------------------------------------------------
// Slots
// -------------------------------------------------------------------------------
A2 = 0.0;
AA[]=deg2rad*{2.77+A2, 4.0+A2, 5.52+A2, 5.56+A2, 5.65+A2, 9.35+A2, 9.44+A2, 9.48+A2, 11+A2, 12.23+A2} ;
For k In {0:#AA[]-1}
cosAA[]+=Cos(AA[k]); sinAA[]+=Sin(AA[k]);
EndFor
pntSlot[]+=newp; Point(newp)={rS5*cosAA[0], rS5*sinAA[0], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS3*cosAA[1], rS3*sinAA[1], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS1*cosAA[2], rS1*sinAA[2], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS2*cosAA[3], rS2*sinAA[3], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS4*cosAA[3], rS4*sinAA[3], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS6*cosAA[4], rS6*sinAA[4], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS6*cosAA[5], rS6*sinAA[5], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS4*cosAA[6], rS4*sinAA[6], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS2*cosAA[6], rS2*sinAA[6], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS1*cosAA[7], rS1*sinAA[7], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS3*cosAA[8], rS3*sinAA[8], 0., pS3};
pntSlot[]+=newp; Point(newp)={rS5*cosAA[9], rS5*sinAA[9], 0., pS3};
// air slot 1
linASlot[]+=newl ; Line(newl)={pntSlot[2], pntSlot[3]};
linASlot[]+=newl ; Line(newl)={pntSlot[3], pntSlot[1]};
linASlot[]+=newl ; Circle(newl)={pntSlot[1], cen, pntSlot[10]};
linASlot[]+=newl ; Line(newl)={pntSlot[10], pntSlot[8]};
linASlot[]+=newl ; Line(newl)={pntSlot[8], pntSlot[9]};
linASlot[]+=newl ; Circle(newl)={pntSlot[9], cen, pntSlot[2]};
Line Loop(newll) = {linASlot[]};
sairslot[] += news ; Plane Surface(sairslot[0]) = {newll-1};
// coil slot 1
linSlot[]+=newl ; Line(newl)={pntSlot[1], pntSlot[0]};
linSlot[]+=newl ; Circle(newl)= {pntSlot[0], pntSlot[4], pntSlot[5]};
linSlot[]+=newl ; Line(newl)={pntSlot[5], pntSlot[6]};
linSlot[]+=newl ; Circle(newl)={pntSlot[6], pntSlot[7],pntSlot[11]};
linSlot[]+=newl ; Line(newl)={pntSlot[11], pntSlot[10]};
Line Loop(newll) = {-linASlot[2],linSlot[]};
sslot[] += news ; Plane Surface(sslot[0]) = {newll-1};
// slots 2 and 3
A2 = 15*deg2rad;
pntSlot0[0] = pntSlot[2];
pntSlot1[0] = pntSlot[9];
For k In{1:2}
pntSlot0[] += Rotate {{0, 0, 1}, {0, 0, 0}, A2} { Duplicata{Point{pntSlot0[k-1]};} };
pntSlot1[] += Rotate {{0, 0, 1}, {0, 0, 0}, A2} { Duplicata{Point{pntSlot1[k-1]};} };
EndFor
For k In{1:2}
sslot[] += Rotate {{0, 0, 1}, {0, 0, 0}, A2} { Duplicata{Surface{sslot[k-1]};} };
sairslot[] += Rotate {{0, 0, 1}, {0, 0, 0}, A2} { Duplicata{Surface{sairslot[k-1]};} };
EndFor
cSlot[]+=newl; Circle(newl) = {pntS[0], cen, pntSlot[2]};
cSlot[]+=newl; Circle(newl) = {pntSlot1[0], cen, pntSlot0[1]};
cSlot[]+=newl; Circle(newl) = {pntSlot1[1], cen, pntSlot0[2]};
cSlot[]+=newl; Circle(newl) = {pntSlot1[2], cen, pntS[2]};
linesslot0[] = CombinedBoundary{ Surface{ sslot[0], sairslot[0] } ;};
linesslot1[] = CombinedBoundary{ Surface{ sslot[1], sairslot[1] } ;};
linesslot2[] = CombinedBoundary{ Surface{ sslot[2], sairslot[2] } ;};
Line Loop(newll) = {-lineMBstator[0],linS[0], cSlot[0],-linesslot0[{4}],
cSlot[1],-linesslot1[{9}],
cSlot[2],-linesslot2[{9}], cSlot[3], -linS[1]};
sairgapS[0]=news; Plane Surface(sairgapS[0]) = {newll-1};
linesslot0[] -= linesslot0[{4}];
linesslot1[] -= linesslot1[{9}];
linesslot2[] -= linesslot2[{9}];
Line Loop(newll) = { cSlot[0], linesslot0[],
cSlot[1], linesslot1[],
cSlot[2], linesslot2[],
cSlot[3], linS[3], -cirS[0], -linS[2]};
sstator[0]=news; Plane Surface(sstator[0]) = -{newll-1};
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------------------
auxlink[]=linS[{1,3}]; // A1
If(SymmetryFactor<8)
// FULL MODEL ==> Rotation of NbrPolesTot*Pi/4
// For simplicity: rotating the interior and exterior boundaries
If (SymmetryFactor>1)
For k In {0:#auxlink[]-1}
auxlink_[] += Rotate {{0, 0, 1}, {0, 0, 0}, 2*Pi/SymmetryFactor-Pi/4} { Duplicata{Line{auxlink[k]};} };
EndFor
auxlink[] = auxlink_[];
EndIf
For k In {1:NbrPolesInModel-1}
cirS[] += Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Line{cirS[{0}]};} };
EndFor
For k In {1:NbrPolesInModel-1}
sstator[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sstator[0]};} };
sairgapS[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sairgapS[0]};} };
sairslot[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sairslot[{0:2}]};} };
sslot[]+= Rotate {{0, 0, 1}, {0, 0, 0}, k*Pi/4} { Duplicata{ Surface{sslot[{0:2}]};} };
EndFor
EndIf
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------------------
// Physical regions
// -------------------------------------------------------------------------------
// -------------------------------------------------------------------------------
Physical Surface(STATOR_FE) = {sstator[]}; // Stator
Physical Surface(STATOR_AIR) = {sairslot[]}; // AirStator
Physical Surface(STATOR_AIRGAP) = {sairgapS[]}; // AirStator for possible torque computation with Maxwell stress tensor
NN = (Flag_Symmetry)?NbrSectStator:NbrSectTotStator;
//For k In {0:NN-1}
// Physical Surface(STATOR_IND+k) = {sslot[k]}; //Inds
//EndFor
Physical Surface(STATOR_IND_AM) = {sslot[{0:NN-1:6}]};
Physical Surface(STATOR_IND_CP) = {sslot[{1:NN-1:6}]};
Physical Surface(STATOR_IND_BM) = {sslot[{2:NN-1:6}]};
If(NbrSectStator>2)
Physical Surface(STATOR_IND_AP) = {sslot[{3:NN-1:6}]};
Physical Surface(STATOR_IND_CM) = {sslot[{4:NN-1:6}]};
Physical Surface(STATOR_IND_BP) = {sslot[{5:NN-1:6}]};
EndIf
Color Pink {Surface{ sslot[{0:NN-1:6}] };} // A-
Color ForestGreen {Surface{ sslot[{1:NN-1:6}] };} // C+
Color PaleGoldenrod{Surface{ sslot[{2:NN-1:6}] };} // B-
If (#sslot[]>=6)
Color Red {Surface{ sslot[{3:NN-1:6}] };} // A+
Color SpringGreen{Surface{ sslot[{4:NN-1:6}] };} // C-
Color Gold {Surface{ sslot[{5:NN-1:6}] };} // B+
EndIf
Physical Line(SURF_EXT) = {cirS[]}; // SurfExt
If(Flag_Symmetry) //Lines for symmetry link
Physical Line(STATOR_BND_A0) = linS[{0,2}];
Physical Line(STATOR_BND_A1) = auxlink[] ;
EndIf
lineMBstator[] += lineMBstatoraux[] ;
If(!Flag_Symmetry)
Physical Line(STATOR_BND_MOVING_BAND) = {lineMBstator[]};
EndIf
If(Flag_Symmetry)
ns = #lineMBstator[];
nns = ns/SymmetryFactor ;
For k In {1:SymmetryFactor}
kk= ((k*nns-1) > ns) ? ns-1 : k*nns-1 ;
Physical Line(STATOR_BND_MOVING_BAND+k-1) = {lineMBstator[{(k-1)*nns:kk}]};
EndFor
k1 = Floor[NbrPolesTot/NbrSect];
k2 = Ceil[NbrPolesTot/NbrSect];
If (k2 > k1)
Physical Line(STATOR_BND_MOVING_BAND+k2-1) = lineMBstator[{(k2-1)*nns:#lineMBstator[]-1}] ;
EndIf
EndIf
// For nice visualisation...
linStator[] = CombinedBoundary{Surface{sstator[]};};
linSlot[] = CombinedBoundary{Surface{sslot[]};};
nicepos_stator[] += {linStator[],linSlot[] };
Color SteelBlue {Surface{sstator[]};}
Color SkyBlue {Surface{sairslot[],sairgapS[]};}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment