Skip to content
Snippets Groups Projects
Commit f9bfbe28 authored by François Henrotte's avatar François Henrotte
Browse files

new model to illustrate Sliding surface 3D

parent b5561896
No related branches found
No related tags found
No related merge requests found
Pipeline #5195 passed
Include "rfpm_common.pro";
// Tolerance for the geometric identification.
tol = 1.1*SlidingSurfaceShift;
// Auxiliary Gmsh functions for OpenCascade solid modelling
Function CheckNumbering
// Check whether BooleanFragments has preserved numbering as expected
For num In {0:#Volumes()-1}
If( Volumes[num] != vv(num) )
Error("BooleanFragment changed volume tag %g -> %g",
Volumes[num], vv(num));
Abort;
EndIf
EndFor
/* EXPLAIN
For num In {0:#Surfaces()-1}
If( Surfaces(num) != vv(num+#Volumes()) )
Error("BooleanFragment changed surface tag %g -> %g",
Surfaces(num), vv(num+#Volumes()));
Abort;
EndIf
EndFor
*/
Return
Function ChangeCoordinates
If( CoordinateSystem == 1) // cylindrical coordinates
coord() = {Hypot(point(0),point(1)),Atan2(point(1),point(0)),point(2)};
Else // Cartesian coordinates
coord() = point();
EndIf
Return
/*
Due to internal relabelling that cannot be controlled by the user,
the consistent way to identify the surface labels (tags)
of a solid modelled with the OpenCascade factory is to proceed
geometrically.
The surface tags that are needed are in general
pieces of the outer boundary of the solid
where specific boundary conditions are to be imposed.
Very often, the volume occupied by a FE model
is a parallelepipedic block in some coordinate system
(Cartesian, cylindrical, ...)
The boundary conditions are thus associated
with the faces of that block and can therefore be identified
by means of thin bounding boxes in the coordinate space.
The function 'Volume2Cube' returns
the parallelepipedic coordinate block of a list of volumes
'Volumes'.
The function 'Cube2Face' returns the list of faces of the model
that have all their points located inside a thin box encosing
the 'FaceId'th face of that coordinate block.
*/
Function Volume2Cube
//
points() = PointsOf{ Volume{ Volumes() }; };
point() = Point { points[0] };
Call ChangeCoordinates;
cube() = {coord(0),coord(1),coord(2),coord(0),coord(1),coord(2)};
For p In { 1:#points()-1 }
point() = Point { points(p) };
Call ChangeCoordinates;
cube(0) = (coord(0) < cube(0) ? coord(0) : cube(0) );
cube(1) = (coord(1) < cube(1) ? coord(1) : cube(1) );
cube(2) = (coord(2) < cube(2) ? coord(2) : cube(2) );
cube(3) = (coord(0) > cube(3) ? coord(0) : cube(3) );
cube(4) = (coord(1) > cube(4) ? coord(1) : cube(4) );
cube(5) = (coord(2) > cube(5) ? coord(2) : cube(5) );
EndFor
Return
Function Cube2Face
bbox() = cube();
bbox(Modulo(FaceId+3,6)) = bbox(FaceId);
//Printf("bbox=", bbox());
face()={};
surfaces() = Boundary{ Volume{ Volumes() }; };
For s In { 0:#surfaces()-1 }
NbIn = 0;
points() = PointsOf { Surface { surfaces(s) } ; };
For p In { 0:#points()-1 }
point() = Point { points(p) };
Call ChangeCoordinates;
If( (CoordinateSystem == 1) && (coord(0)==0) )
coord(1) = bbox(1);
EndIf
If( (coord(0)>bbox(0)-tol) &&
(coord(1)>bbox(1)-tol) &&
(coord(2)>bbox(2)-tol) &&
(coord(0)<bbox(3)+tol) &&
(coord(1)<bbox(4)+tol) &&
(coord(2)<bbox(5)+tol) )
NbIn++;
EndIf
EndFor
If( NbIn == #points() )
face() += surfaces(s);
EndIf
EndFor
Return
SetFactory("OpenCASCADE");
Solver.AutoMesh = 2;
Mesh.CharacteristicLengthMin = lc;
Mesh.CharacteristicLengthMax = lc;
// Hide mesh in post-processing views
Mesh.VolumeEdges = 0;
Mesh.SurfaceEdges = 0;
L_r = 44*mm;
L_m = 40*mm;
L_s = 44*mm;
pp = 3;
anglePole = Pi/pp;
e_in = 0;
e_r = 3.5*mm;
b_depth = 0.5*mm;
e_m = 5.3*mm;
e_sl = 0;
e_ar = 1*mm;
e_ss = 0.1*mm; // air ring width
e_a = 1*mm; // width stator air gap
e_s = 5*mm; // width stator core
sRM = 0.75;
A_infinity = 66.75*mm;
R_infinity = 21.75*mm;
// Radii
R_in = e_in;
R_r = R_in + e_r;
R_m_in = R_r - b_depth;
R_m = R_m_in + e_m;
R_sl = R_m + e_sl;
R_ar = R_sl + e_ar;
R_ss = R_ar + e_ss;
R_s = R_ar + e_a;
R_o = R_s + e_s;
// R O T O R
p1r = newp; Point(p1r) = {R_in,0,0};
p2r = newp; Point(p2r) = {R_in,0,L_r/2};
p3r = newp; Point(p3r) = {R_r,0,L_r/2};
p4r = newp; Point(p4r) = {R_r,0,0};
l1r = newl; Line(l1r) = {p1r,p2r};
l2r = newl; Line(l2r) = {p2r,p3r};
l3r = newl; Line(l3r) = {p3r,p4r};
l4r = newl; Line(l4r) = {p4r,p1r};
ll_r = newll; Line Loop(ll_r) = {l1r, l2r, l3r, l4r};
sRotor = news; Plane Surface(sRotor) = {ll_r};
vRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
Surface{sRotor};
};
// Magnets
p1m = newp; Point(p1m) = {R_m_in,0,0};
p2m = newp; Point(p2m) = {R_m_in,0,L_m/2};
p3m = newp; Point(p3m) = {R_m,0,L_m/2};
p4m = newp; Point(p4m) = {R_m,0,0};
l1m = newl; Line(l1m) = {p1m,p2m};
l2m = newl; Line(l2m) = {p2m,p3m};
l3m = newl; Line(l3m) = {p3m,p4m};
l4m = newl; Line(l4m) = {p4m,p1m};
ll_m = newll; Line Loop(ll_m) = {l1m, l2m, l3m, l4m};
sMagnet = news; Plane Surface(sMagnet) = {ll_m};
Rotate{ {0,0,1}, {0,0,0}, anglePole*(1-sRM)/2 } {
Surface{sMagnet};
}
vMagnetExt() = Extrude {{0,0,1}, {0,0,0}, anglePole*sRM} {
Surface{sMagnet};
};
// Air rotor
p1ar = newp; Point(p1ar) = {R_in,0,0};
p2ar = newp; Point(p2ar) = {R_in,0,A_infinity/2};
p3ar = newp; Point(p3ar) = {R_ar,0,A_infinity/2};
p4ar = newp; Point(p4ar) = {R_ar,0,0};
l1ar = newl; Line(l1ar) = {p1ar,p2ar};
l2ar = newl; Line(l2ar) = {p2ar,p3ar};
l3ar = newl; Line(l3ar) = {p3ar,p4ar};
l4ar = newl; Line(l4ar) = {p4ar,p1ar};
ll_ar = newll; Line Loop(ll_ar) = {l1ar, l2ar, l3ar, l4ar};
sAirRotor = news; Plane Surface(sAirRotor) = {ll_ar};
vAirRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
Surface{sAirRotor};
};
// Create volumes with boolean operations
vMagnet = vMagnetExt[1];
vRotor = BooleanDifference{
Volume{vRotorExt[1]}; Delete;
}{
Volume{vMagnet};
};
vAirRotor = BooleanDifference{
Volume{vAirRotorExt[1]}; Delete;
}{
Volume{vRotor}; Volume{vMagnet};
};
// Glue up rotor pieces with 'BooleanFragment'.
// This operation is necessary, as it connects the different
// parts of the model together,
// but it entails an internal renumbering of surfaces, lines, etc...
// Hence the need for a geometrical identification afterwards (see below).
Volumes() = { vMagnet, vRotor, vAirRotor };
VolumesRotor() = BooleanFragments{
Volume { Volumes() }; Delete; }{};
vv() = VolumesRotor(); Call CheckNumbering;
// S T A T O R
// Air ring
// This region is specific to the sliding surface mechanism
// to represent rotor motion.
// It is a thin volumic region of width 'e_ss' in contact
// with the 'SlidingMaster' surface.
p1ss = newp; Point(p1ss) = {R_ar,0,SlidingSurfaceShift};
p2ss = newp; Point(p2ss) = {R_ar,0,SlidingSurfaceShift+A_infinity/2};
p3ss = newp; Point(p3ss) = {R_ss,0,A_infinity/2};
p4ss = newp; Point(p4ss) = {R_ss,0,0};
l1ss = newl; Line(l1ss) = {p1ss,p2ss};
l2ss = newl; Line(l2ss) = {p2ss,p3ss};
l3ss = newl; Line(l3ss) = {p3ss,p4ss};
l4ss = newl; Line(l4ss) = {p4ss,p1ss};
ll_ss = newll; Line Loop(ll_ss) = {l1ss, l2ss, l3ss, l4ss};
sAirRing = news; Plane Surface(sAirRing) = {ll_ss};
vAirRingExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
Surface{sAirRing};
};
vAirRing = vAirRingExt[1] ;
// Air stator
p1as = newp; Point(p1as) = {R_ss,0,0};
p2as = newp; Point(p2as) = {R_ss,0,A_infinity/2};
p3as = newp; Point(p3as) = {R_infinity,0,A_infinity/2};
p4as = newp; Point(p4as) = {R_infinity,0,0};
l1as = newl; Line(l1as) = {p1as,p2as};
l2as = newl; Line(l2as) = {p2as,p3as};
l3as = newl; Line(l3as) = {p3as,p4as};
l4as = newl; Line(l4as) = {p4as,p1as};
ll_as = newll; Line Loop(ll_as) = {l1as, l2as, l3as, l4as};
sAirStator = news; Plane Surface(sAirStator) = {ll_as};
vAirStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
Surface{sAirStator};
};
// Stator core
p1s = newp; Point(p1s) = {R_s,0,0};
p2s = newp; Point(p2s) = {R_s,0,L_s/2};
p3s = newp; Point(p3s) = {R_o,0,L_s/2};
p4s = newp; Point(p4s) = {R_o,0,0};
l1s = newl; Line(l1s) = {p1s,p2s};
l2s = newl; Line(l2s) = {p2s,p3s};
l3s = newl; Line(l3s) = {p3s,p4s};
l4s = newl; Line(l4s) = {p4s,p1s};
ll_s = newll; Line Loop(ll_s) = {l1s, l2s, l3s, l4s};
sStator = news; Plane Surface(sStator) = {ll_s};
vStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
Surface{sStator};
};
// Create volumes by boolean operations
vStator = vStatorExt[1];
vAirStator = BooleanDifference{
Volume{vAirStatorExt[1]}; Delete;
}{
Volume{vStator};
};
// Glue up stator pieces
Volumes() = { vAirRing, vStator, vAirStator };
VolumesStator() = BooleanFragments{
Volume { Volumes() }; Delete; }{};
Printf("VolumesStator=", VolumesStator());
Printf("VolumesRotor=", VolumesRotor());
// Recover Physical Surfaces geometrically
CoordinateSystem = 1; // cylindrical coordinates
Volumes() = VolumesRotor(); Call Volume2Cube;
Printf("cube rotor", cube());
FaceId = 1; Call Cube2Face; // no statement here FIXME
RotorPerMaster() = face();
FaceId = 2; Call Cube2Face;
RotorBottom() = face();
FaceId = 3; Call Cube2Face;
SlidingSlave() = face();
FaceId = 4; Call Cube2Face;
RotorPerSlave() = face();
FaceId = 5; Call Cube2Face;
RotorTop() = face();
/* control if the geometrical identification succeeded
Physical Surface("Test1", 101) = { RotorPerMaster() };
Physical Surface("Test2", 102) = { RotorTop() };
Physical Surface("Test3", 103) = { SlidingSlave() };
Physical Surface("Test4", 104) = { RotorBottom() };
Physical Surface("Test5", 105) = { RotorPerSlave() };
*/
Volumes() = VolumesStator(); Call Volume2Cube;
Printf("cube stator", cube());
FaceId = 0; Call Cube2Face;
SlidingMaster() = face();
FaceId = 1; Call Cube2Face;
StatorPerMaster() = face();
FaceId = 2; Call Cube2Face;
StatorBottom() = face();
FaceId = 3; Call Cube2Face;
StatorOuter() = face();
FaceId = 4; Call Cube2Face;
StatorPerSlave() = face();
FaceId = 5; Call Cube2Face;
StatorTop() = face();
/* control if the geometrical identification succeeded
Physical Surface("Test11", 111) = { SlidingMaster() };
Physical Surface("Test12", 112) = { StatorPerMaster() };
Physical Surface("Test13", 113) = { StatorBottom() };
Physical Surface("Test14", 114) = { StatorOuter() };
Physical Surface("Test15", 115) = { StatorPerSlave() };
Physical Surface("Test16", 116) = { StatorTop() };
*/
// Make AirRing a transfinite region
// The width of the air ring is arbitrary, but it is in general thin,
// especially in machine models where it has to fit into the air gap.
// To improve the quality of the mesh,
// it is thus appropriate to have it transfinite in practice.
surfaces() = Boundary { Volume{ vAirRing }; };
//Printf("%g %g", vAirRing, #surfaces());
lines() = Boundary { Surface{ SlidingSlave() }; };
For num In { 0:#surfaces()-1 }
lines() += Boundary { Surface{ surfaces[num] }; };
Printf("%g %g", surfaces[num], #lines());
EndFor
For num In { 0:#lines()-1 }
points() = PointsOf { Line { lines(num) }; };
If( #points() != 2 )
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
// Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
If( Fabs( Atan2( ptA(1), ptA(0) ] - Atan2[ ptB(1), ptB(0) ) ) > tol )
// line along theta-axis
Transfinite Line{ lines[num] } = NbrDiv+1;
ElseIf( Fabs[ ptA(2) - ptB(2) ] > tol )
// line along z-axis
Transfinite Line{ lines(num) } = Fabs( ptA(2) - ptB(2) )/lc-1;
Else
// line along r-axis
Transfinite Line{ lines(num) } = 2;
EndIf
EndFor
Transfinite Surface{ surfaces() } ;
Transfinite Surface{ SlidingSlave() } Right;
Transfinite Volume{ vAirRing };
// Identify 'SlidingSubmaster' Curve
AngleSubSurfaces = (ModelAngleMax-ModelAngleMin)*deg;
lines() = Boundary { Surface{ SlidingMaster() }; };
For l In { 0:#lines()-1 }
points() = PointsOf { Line { lines(l) }; };
If( #points() != 2 )
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
// Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces) < tol &&
Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces) < tol )
SlidingSubmaster = lines(l);
EndIf
EndFor
// Identify 'SlidingSubslave' curve
lines() = Boundary { Surface{ SlidingSlave() }; };
For l In { 0:#lines()-1 }
points() = PointsOf { Line { lines(l) }; };
If( #points() != 2 )
Error("This is unexpected"); Abort;
EndIf
ptA() = Point { points(0) };
ptB() = Point { points(1) };
If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces ) < tol &&
Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces ) < tol )
SlidingSubslave = lines(l);
EndIf
EndFor
// P H Y S I C A L R E G I O N S
Physical Volume("Rotor" , 1) = {vRotor};
Physical Volume("Magnet" , 2) = {vMagnet};
Physical Volume("Stator" , 3) = {vStator};
Physical Volume("AirRotor" , 4) = {vAirRotor};
Physical Volume("AirStator", 5) = {vAirStator};
Physical Volume("AirRing" , 6) = {vAirRing};
Physical Surface("Outer" , 10) = {StatorOuter()};
Physical Surface("Top" , 11) = {RotorTop(),StatorTop()};
Physical Surface("Bottom" , 12) = {RotorBottom(),StatorBottom()};
Physical Surface("RotorPerMaster" , 13) = {RotorPerMaster()};
Physical Surface("RotorPerSlave" , 14) = {RotorPerSlave()};
Physical Surface("StatorPerMaster", 15) = {StatorPerMaster()};
Physical Surface("StatorPerSlave" , 16) = {StatorPerSlave()};
Physical Surface("SlidingMaster" , 17) = {SlidingMaster()};
Physical Surface("SlidingSlave" , 18) = {SlidingSlave()};
Physical Line("Sliding_Submaster", 20) = {SlidingSubmaster};
Physical Line("Sliding_Subslave", 21) = {SlidingSubslave};
// Ensure identical meshes on the periodic faces
// The 'Air ring' being meshed transfinite with prisms,
// the periodic surface routine issues a harmless warning
// Info : No periodic vertices in surface xxx
For num In {0:#RotorPerMaster()-1}
Periodic Surface { RotorPerSlave(num) }
= { RotorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
EndFor
For num In {0:#StatorPerMaster()-1}
Periodic Surface { StatorPerSlave(num) }
= { StatorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
EndFor
treeLines() = CombinedBoundary { Physical Surface { 10 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 11 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 12 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 13 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 14 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 15 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 16 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 17 }; } ;
treeLines() += CombinedBoundary { Physical Surface { 18 }; } ;
Physical Line("treeLines", 22) = { treeLines() };
/*
3D Magnetostatics with
- OpenCascade solid modelling
- geometrical identfication of faces
- antiperiodic boundary conditions
- sliding surface
*/
Include "rfpm_common.pro";
DefineConstant[
R_ = {"Static", Choices {"Static", "QuasiStatic"}, Visible 1,
Name "GetDP/1ResolutionChoices"},
C_ = {"-solve -v 4 -pos", Name "GetDP/9ComputeCommand", Visible 0}
P_ = { "Check_Periodicity", Choices {"Fields", "Check_Periodicity"}, Visible 0,
Name "GetDP/2PostOperationChoices"}
];
Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
RotorPosition = deg*
DefineNumber[0, Name "Options/Rotor position [deg]",
Visible !Flag_QuasiStatic];
AngularStep = deg*
DefineNumber[1, Name "Options/1Angular step [deg]",
Visible Flag_QuasiStatic];
NumStep =
DefineNumber[20, Name "Options/1Number of steps",
Visible Flag_QuasiStatic];
rpm = 2.*Pi/60. ;
w1 = 1000*rpm ;
Mag_DTime = AngularStep/Fabs[w1] ;
Mag_TimeMax = Mag_DTime*NumStep;
Group {
// Geometrical regions
Rotor = Region[1];
Magnet = Region[2];
Stator = Region[3];
AirRotor = Region[4];
AirStator = Region[5];
AirRing = Region[6];
Outer = Region[10];
Top = Region[11];
Bottom = Region[12];
Sur_RotorPerMaster = Region[13];
Sur_RotorPerSlave = Region[14];
Sur_StatorPerMaster = Region[15];
Sur_StatorPerSlave = Region[16];
Sur_SlidingMaster = Region[17];
Sur_SlidingSlave = Region[18];
Lin_SlidingSubmaster = Region[20];
Lin_SlidingSubslave = Region[21];
Lin_TreeLines = Region[22];
// Abstract regions
Iron = Region[ { Rotor, Stator } ];
Air = Region[ { AirRotor, AirStator, AirRing } ];
Vol_Conductors = Region[ {} ] ;
Vol_Magnets = Region[ Magnet ];
Vol_Mag = Region[ { Air, Magnet, Iron } ];
Vol_Moving = Region[ { Rotor, AirRotor, Magnet } ] ;
Sur_Dirichlet = Region[ { Outer, Top, Bottom } ];
Sur_Link = Region[ {Sur_SlidingMaster, Sur_SlidingSlave,
Sur_StatorPerMaster, Sur_StatorPerSlave,
Sur_RotorPerMaster, Sur_RotorPerSlave} ] ;
Dom_Hcurl_a = Region[ { Vol_Mag, Sur_Link } ]; // No Neumann surface
Sur_Tree_Sliding = Region[ { Sur_SlidingMaster, Sur_SlidingSlave } ];
Lin_Tree = Region[ Lin_TreeLines ];
Sur_Tree = Region[ { Sur_Link, Sur_Dirichlet } ];
Vol_Tree = Region[ Vol_Mag ];
}
Function{
mu0 = Pi*4e-7;
br = 1.2;
mur = 3;
nu[ Iron ] = 1/(mur*mu0);
nu[ Air ] = 1/mu0;
nu[ Magnet ] = 1/mu0;
br[] = br*Unit[ Vector[X[], Y[], 0] ];
js[] = 0;
}
Function {
// Sliding surface:
// a1 is the angular span of Region and RegionRef.
// a0 is the mesh step of the meshes of
// a2[] is the angular position (in radian) of the rotor
// relative to its reference position (as in the msh file)
// assumed to be aligned with the stator.
// a3[] is the shear angle of region AIRBM
// (smaller than half the discretization step of the sliding surface
// to adapt to a2[] values that are not multiple of this step).
// AlignWithMaster[] maps a point of coordinates {X[], Y[], Z[]} onto
// its image in the open set RegionRef-SubRegionRef by the symmetry mapping.
// Coef[] is evaluated on Master nodes
// fFloor[] is a safe version of Floor[] for real represented int variables
// fRem[a,b] substracts from a multiple of b so that the result is in [0,1[
Periodicity = -1. ; // -1 for antiperiodicity, 1 for periodicity
RotatePZ[] = Rotate[ XYZ[], 0, 0, $1 ] ;
Tol = 1e-8 ; fFloor[] = Floor[ $1 + Tol ] ; fRem[] = $1 - fFloor[ $1 / $2 ] * $2;
deg = Pi/180;
a1 = (ModelAngleMax-ModelAngleMin)*deg ;
a0 = a1/NbrDiv ; // angular span of one moving band element
a2[] = $RotorPosition ;
a3[] = ( ( fRem[ a2[], a0 ]#2 <= 0.5*a0 ) ? #2 : #2-a0 ) ;
AlignWithMaster[] = RotatePZ[ - fFloor[ (Atan2[ Y[], X[] ] - a3[] )/a1 ] * a1 ] ;
RestoreRef[] = XYZ[] - Vector[ 0, 0, SlidingSurfaceShift ] ;
Coef[] = ( fFloor[ ((Atan2[ Y[], X[] ] - a2[] )/a1) ] % 2 ) ? Periodicity : 1. ;
}
Group {
DefineGroup[ DomainInf ];
DefineVariable[ Val_Rint, Val_Rext ];
}
Jacobian {
{ Name Vol ;
Case { { Region DomainInf ; Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
{ Region All ; Jacobian Vol ; }
}
}
{ Name Sur ;
Case { { Region All ; Jacobian Sur ; }
}
}
}
Integration {
{ Name Int ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 21 ; }
}
}
}
}
}
Constraint {
{ Name a ;
Case {
{ Region Sur_Dirichlet ; Value 0. ; }
// Periodicity condition on lateral faces
{ Type Link ; Region Sur_StatorPerSlave ; RegionRef Sur_StatorPerMaster ;
ToleranceFactor 1e-8;
Coefficient Periodicity ; Function RotatePZ[ -a1 ] ; }
{ Type Link ; Region Sur_RotorPerSlave ; RegionRef Sur_RotorPerMaster ;
ToleranceFactor 1e-8;
Coefficient Periodicity ; Function RotatePZ[ -a1 ] ; }
// Sliding surface
{ Type Link ; Region Sur_SlidingSlave ; SubRegion Lin_SlidingSubslave ;
RegionRef Sur_SlidingMaster ; SubRegionRef Lin_SlidingSubmaster ;
ToleranceFactor 1e-8;
Coefficient Coef[] ; Function AlignWithMaster[] ;
FunctionRef RestoreRef[] ;
}
}
}
{ Name GaugeCondition_a ; Type Assign ;
Case {
{ Region Vol_Tree ; Value 0. ;
SubRegion Region[ { Sur_Tree, Lin_Tree} ] ;
SubRegion2 Region[ Sur_Tree_Sliding, AlignedWith Z ] ;
}
}
}
}
FunctionSpace {
{ Name Hcurl_a; Type Form1;
BasisFunction {
{ Name se; NameOfCoef ae; Function BF_Edge;
Support Dom_Hcurl_a ; Entity EdgesOf[ All ]; }
}
Constraint {
{ NameOfCoef ae; EntityType EdgesOf ;
NameOfConstraint a; }
If( Gauge==1 ) // spanning tree gauging
{ NameOfCoef ae; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
NameOfConstraint GaugeCondition_a ; }
EndIf
}
}
}
Formulation {
{ Name MagSta_a; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Vol_Mag ; Jacobian Vol ; Integration Int ; }
Galerkin { [ -nu[] * br[] , {d a} ] ;
In Vol_Magnets ; Jacobian Vol ; Integration Int ; }
Galerkin { [ -js[] , {a} ] ; In Vol_Conductors ;
Jacobian Vol ; Integration Int ; }
}
}
}
Resolution {
{ Name Static ;
System {
{ Name Sys_Mag ; NameOfFormulation MagSta_a ;}
}
Operation {
Evaluate[$RotorPosition = RotorPosition];
ChangeOfCoordinates [ NodesOf[ Vol_Moving ], RotatePZ[ a2[] ] ] ;
ChangeOfCoordinates [ NodesOf[ Sur_SlidingMaster ], RotatePZ[ a3[] ] ] ;
Evaluate[ $a2 = a2[] ];
Evaluate[ $a3 = a3[] ];
Evaluate[ $aa = Fmod[ a2[], a0 ] ];
Evaluate[ $bb = fRem[ a2[], a0 ] ];
Print[ {$RotorPosition/deg, $a2/deg, $a3/deg, $aa/deg, $bb/deg},
Format "wt=%e a2=%e a3=%e %e %e"] ;
UpdateConstraint[Sys_Mag] ;
Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
PostOperation [Fields];
}
}
{ Name QuasiStatic ;
System {
{ Name Sys_Mag ; NameOfFormulation MagSta_a ; }
}
Operation {
Evaluate[$RotorPosition = RotorPosition];
Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
PostOperation [Fields];
TimeLoopTheta{
Time0 0 ; TimeMax Mag_TimeMax ; DTime Mag_DTime ; Theta 1. ;
Operation {
Evaluate[$RotorPosition = RotorPosition + w1*$Time];
ChangeOfCoordinates[NodesOf[ Vol_Moving ], RotatePZ[ w1*$DTime ] ] ;
ChangeOfCoordinates[NodesOf[ Sur_SlidingMaster ], RotatePZ[ a3[] ] ] ;
UpdateConstraint[Sys_Mag] ;
Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
PostOperation [Fields];
ChangeOfCoordinates [ NodesOf[ Sur_SlidingMaster ], RotatePZ[ -a3[] ] ] ;
}
}
}
}
}
PostProcessing {
{ Name MagSta_a ; NameOfFormulation MagSta_a ;
PostQuantity {
{ Name b ; Value { Local { [ {d a} ];
In Vol_Mag ; Jacobian Vol; } } }
{ Name bsurf ; Value { Local { [ {d a} ];
In Sur_Link ; Jacobian Sur; } } }
{ Name asurf ; Value { Local { [ {a} ];
In Sur_Link ; Jacobian Sur; } } }
{ Name br ; Value { Local { [ br[] ];
In Vol_Magnets ; Jacobian Vol; } } }
{ Name a ; Value { Local { [ {a} ];
In Vol_Mag ; Jacobian Vol; } } }
}
}
}
PostOperation Fields UsingPost MagSta_a {
Print[ b, OnElementsOf Region[ {Vol_Mag} ],
LastTimeStepOnly, File "b.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ArrowSizeMax = 100;",
"View[l].CenterGlyphs = 1;",
"View[l].GlyphLocation = 1;",
"View[l].RangeType = 1;",
"View[l].SaturateValues = 1;",
"View[l].CustomMax = 2;",
"View[l].CustomMin = 1e-06;",
"View[l].LineWidth = 2;",
"View[l].VectorType = 4;" ] ,
File "tmp.geo", LastTimeStepOnly] ;
}
PostOperation {
{ Name pos; NameOfPostProcessing MagSta_a;
Operation {
Print[ b, OnElementsOf Dom_Hcurl_a, File "b.pos" ];
// Print[W_mag_Airgap[AIRGAP], OnGlobal, Format TimeTable,
// File StrCat[resPath, "Energy_airgap.pos"]];
// Print[Flux[PM], OnGlobal, Format TimeTable,
// File >> StrCat[resPath, "Flux.pos"]];
}
}
}
PostOperation Check_Periodicity UsingPost MagSta_a {
PrintGroup[ EdgesOfTreeIn[ { Vol_Tree }, StartingOn { Sur_Tree } ],
In Vol_Tree, File "Tree.pos"];
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ColorTable = { DarkRed };",
"View[l].LineWidth = 5;"] ,
File "tmp.geo", LastTimeStepOnly] ;
Print[ bsurf, OnElementsOf Region[ {Sur_SlidingMaster} ],
LastTimeStepOnly, File "bsm.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ArrowSizeMax = 100;",
"View[l].CenterGlyphs = 0;",
"View[l].GlyphLocation = 1;",
"View[l].RangeType = 1;",
"View[l].SaturateValues = 1;",
"View[l].CustomMax = 2;",
"View[l].CustomMin = 0;",
"View[l].ScaleType = 1;",
"View[l].LineWidth = 3;",
"View[l].VectorType = 4;" ] ,
File "tmp.geo", LastTimeStepOnly] ;
Print[ bsurf, OnElementsOf Region[ {Sur_SlidingSlave} ],
LastTimeStepOnly, File "bss.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ArrowSizeMax = 100;",
"View[l].CenterGlyphs = 0;",
"View[l].GlyphLocation = 1;",
"View[l].RangeType = 1;",
"View[l].SaturateValues = 1;",
"View[l].CustomMax = 2;",
"View[l].CustomMin = 0;",
"View[l].ScaleType = 1;",
"View[l].LineWidth = 3;",
"View[l].VectorType = 4;" ] ,
File "tmp.geo", LastTimeStepOnly] ;
Print[ bsurf, OnElementsOf Region[ Sur_Link ],
LastTimeStepOnly, File "bsl.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ArrowSizeMax = 100;",
"View[l].CenterGlyphs = 0;",
"View[l].GlyphLocation = 1;",
"View[l].ScaleType = 1;",
"View[l].LineWidth = 3;",
"View[l].VectorType = 4;" ] ,
File "tmp.geo", LastTimeStepOnly] ;
}
mm = 1.e-3 ;
deg = Pi/180.;
lc = mm * DefineNumber[2, Name "Options/Global mesh size (mm)"];
If(lc == 0)
Error("Glabal mesh size cannot be zero. Reset to default."); lc=2*mm;
EndIf
// Number of elements on the sliding surface in rotation direction
NbrDiv = 30*mm/lc;
Gauge = DefineNumber[1, Name "Options/Gauging",
Choices{ 0="MUMPS", 1="Tree"}];
ModelAngleMin = 0. ;
ModelAngleMax = 60. ;
SlidingSurfaceShift = 1.e-5;
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment