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

sliding surface 2D for a linear motor

parent 69fe0859
No related branches found
No related tags found
No related merge requests found
Pipeline #5194 canceled
Include "mrvl_common.pro";
Mesh.SurfaceEdges=0;
lc1 = lc; // fine
lc2 = 3*lc;
lc3 = 8*lc; // coarse
// 1. R A I L - S T A T O R
bnd_slft() = {};
bnd_srgh() = {};
bnd_sbot() = {};
// 1.1 Stator core
pA = newp; Point(pA) = {0,0,0,lc3};
pB = newp; Point(pB) = {L,0,0,lc3};
lAB = newl; Line(lAB) = {pA,pB};
bnd_sbot() += { lAB };
pD = newp; Point(pD) = {L,cs,0,lc2};
lBD = newl; Line(lBD) = {pB,pD};
bnd_srgh() += { lBD };
p5 = pD;
For i In {0:NbrStaPer-1}
p1 = newp; Point(p1) = {L-i*(ws+vs)-vs/2,cs,0,lc2};
Line(newl) = {p5,p1};
p2 = newp; Point(p2) = {L-i*(ws+vs)-vs/2,cs+hs,0,lc1};
Line(newl) = {p1,p2};
p3 = newp; Point(p3) = {L-i*(ws+vs)-vs/2-ws,cs+hs,0,lc1};
Line(newl) = {p2,p3};
p4 = newp; Point(p4) = {L-i*(ws+vs)-vs/2-ws,cs,0,lc2};
Line(newl) = {p3,p4};
p5 = newp; Point(p5) = {L-i*(ws+vs)-vs-ws,cs,0,lc2};
Line(newl) = {p4,p5};
EndFor
pC = p5;
lCA = newl; Line(lCA)={pC,pA};
bnd_slft() += { lCA };
ll_sc = newll; Line Loop(ll_sc) = { lAB, lBD ... lCA };
sStatorCore = news; Plane Surface(sStatorCore) = {ll_sc};
// 1.2 Air stator
pE = newp; Point(pE) = {0,cs+hs+ag/3,0,lc};
pF = newp; Point(pF) = {L,cs+hs+ag/3,0,lc};
lDF = newl; Line(lDF)={pD,pF};
lFE = newl; Line(lFE)={pF,pE};
lEC = newl; Line(lEC)={pE,pC};
bnd_srgh() += {lDF};
bnd_slft() += {lEC};
lags = lFE;
ll_as = newll; Line Loop(ll_as) = { (-(lCA-1)) ... (-(lBD+1)), lDF, lFE, lEC };
sAirStator = news; Plane Surface(sAirStator) = {ll_as};
pGz = newp; Point(pGz) = {0,cs+hs+2*ag/3,SlidingSurfaceShift,lc1};
pHz = newp; Point(pHz) = {L,cs+hs+2*ag/3,SlidingSurfaceShift,lc1};
l1 = -lFE;
l2 = newl; Line(l2) = {pF,pHz};
l3 = newl; Line(l3) = {pHz,pGz};
l4 = newl; Line(l4) = {pGz,pE};
bnd_srgh() += {l2};
bnd_slft() += {l4};
lags = l1;
lagsz = l3;
ll_al = newll; Line Loop(ll_al) = { l1, l2 ... l4 };
sAirLayer = news; Plane Surface(sAirLayer) = {ll_al};
// 2. T R A N S L A T E U R ( M O V E R )
bnd_mrlft() = {};
bnd_mrgh() = {};
bnd_mtop() = {};
// 2.1 Air gap layer
pG = newp; Point(pG) = {0,cs+hs+2*ag/3,0,lc1};
pH = newp; Point(pH) = {L,cs+hs+2*ag/3,0,lc1};
pI = newp; Point(pI) = {0,cs+hs+hb ,0,lc3};
pJ = newp; Point(pJ) = {L,cs+hs+hb,0,lc3};
l1 = newl; Line(l1) = {pG,pH};
l2 = newl; Line(l2) = {pH,pJ};
l3 = newl; Line(l3) = {pJ,pI};
l4 = newl; Line(l4) = {pI,pG};
bnd_mrgh() += {l2};
bnd_mtop() += {l3};
bnd_mlft() += {l4};
lagm = l1;
ll_am = newll; Line Loop(ll_am) = { l1, l2 ... l4 };
// 2.2 Mover core and coils
ym = cs+hs+ag;
Lm = 6*wt+5*vt;
pI = newp; Point(pI) = {(L-Lm)/2,ym+ht+ct,0,lc2};
pJ = newp; Point(pJ) = {(L+Lm)/2,ym+ht+vt,0,lc2};
llCoils() = {};
llCore() = {};
pA = newp+1;
For i In {0:5}
p1 = newp; Point(p1) = {(L-Lm)/2-vw +i*(wt+vt),ym+ht,0,lc};
If(i>0)
l0 = newl; Line(l0) = {p4,p1}; // connect with previous tooth
EndIf
p2 = newp; Point(p2) = {(L-Lm)/2 +i*(wt+vt),ym+ht,0,lc2};
p3 = newp; Point(p3) = {(L-Lm)/2+wt +i*(wt+vt),ym+ht,0,lc2};
p4 = newp; Point(p4) = {(L-Lm)/2+wt+vw+i*(wt+vt),ym+ht,0,lc2};
p5 = newp; Point(p5) = {(L-Lm)/2-vw +i*(wt+vt),ym+ht-hw,0,lc2};
p6 = newp; Point(p6) = {(L-Lm)/2 +i*(wt+vt),ym+ht-hw,0,lc2};
p7 = newp; Point(p7) = {(L-Lm)/2+wt +i*(wt+vt),ym+ht-hw,0,lc2};
p8 = newp; Point(p8) = {(L-Lm)/2+wt+vw+i*(wt+vt),ym+ht-hw,0,lc2};
p9 = newp; Point(p9) = {(L-Lm)/2 +i*(wt+vt),ym,0,lc1};
p10 = newp; Point(p10) = {(L-Lm)/2+wt+i*(wt+vt),ym,0,lc1};
l1 = newl; Line(l1) = {p5,p6};
l2 = newl; Line(l2) = {p6,p2};
l3 = newl; Line(l3) = {p2,p1};
l4 = newl; Line(l4) = {p1,p5};
ll1 = newll; Line Loop(ll1) = { l1 ... l4 };
l5 = newl; Line(l5) = {p7,p8};
l6 = newl; Line(l6) = {p8,p4};
l7 = newl; Line(l7) = {p4,p3};
l8 = newl; Line(l8) = {p3,p7};
ll2 = newll; Line Loop(ll2) = { l5 ... l8 };
l9 = newl; Line(l9) = {p6,p9};
l10 = newl; Line(l10) = {p9,p10};
l11 = newl; Line(l11) = {p10,p7};
llCoils() += {ll1, ll2};
llCore() += { -l2, l9 ... l11, -l8};
If(i>0)
llCore() += {l0,-l3};
EndIf
If(i<5)
llCore() += {-l7};
EndIf
EndFor
pB = p3;
l1 = newl; Line(l1) = {pB,pJ};
l2 = newl; Line(l2) = {pJ,pI};
l3 = newl; Line(l3) = {pI,pA};
sCoils() = {};
For i In {0:#llCoils()-1}
sc = news; Plane Surface(sc) = {llCoils(i)};
sCoils() += {sc};
EndFor
llc = newll; Line Loop(llc) = { llCore(), l1 ... l3 };
sMoverCore = news; Plane Surface(sMoverCore) = {llc()};
sAirMover = news; Plane Surface(sAirMover) = {ll_am,-llc,-llCoils()};
Transfinite Line{ lags, lagsz, lagm } = NbrDiv+1 Using Power 1.;
// P H Y S I C A L R E G I O N S
Physical Surface("StatorCore", 1) = {sStatorCore};
Physical Surface("MoverCore" , 2) = {sMoverCore};
Physical Surface("AirStator" , 3) = {sAirStator};
Physical Surface("AirLayer" , 4) = {sAirLayer};
Physical Surface("AirMover" , 5) = {sAirMover};
Physical Surface("COILAP", 6) = {sCoils(0),sCoils(7)};
Physical Surface("COILAN", 7) = {sCoils(1),sCoils(6)};
Physical Surface("COILBP", 8) = {sCoils(2),sCoils(9)};
Physical Surface("COILBN", 9) = {sCoils(3),sCoils(8)};
Physical Surface("COILCP", 10) = {sCoils(4),sCoils(11)};
Physical Surface("COILCN", 11) = {sCoils(5),sCoils(10)};
Physical Line("Noflux" , 12) = {bnd_sbot(),bnd_mtop()};
Physical Line("MoverPerMaster" , 13) = {bnd_mlft()};
Physical Line("MoverPerSlave" , 14) = {bnd_mrgh()};
Physical Line("StatorPerMaster", 15) = {bnd_slft()};
Physical Line("StatorPerSlave" , 16) = {bnd_srgh()};
Physical Line("SlidingMaster" , 17) = {lagm()};
Physical Line("SlidingSlave" , 18) = {lagsz()};
Physical Point("Sliding_Submaster", 20) = {pHz};
Physical Point("Sliding_Subslave", 21) = {pH};
// Ensure identical meshes on the periodic faces
For num In {0:#bnd_mlft()-1}
Periodic Line { bnd_mrgh(num) }
= { bnd_mlft(num) } Translate {L,0,0} ;
EndFor
For num In {0:#bnd_slft()-1}
Periodic Line { bnd_srgh(num) }
= { bnd_slft(num) } Translate {L,0,0} ;
EndFor
nicepos[] = Boundary{ Surface { Surface '*'};};
Physical Line("NICEPOS", 19) = {nicepos[]};
// ALSTOM linear reluctance motor
// Fr. Henrotte
Include "mrvl_common.pro";
DefineConstant[
R_ = {"Static", Choices {"Static", "QuasiStatic"}, Visible 1,
Name "GetDP/1ResolutionChoices"},
C_ = {"-solve -v 4 -pos -2", Name "GetDP/9ComputeCommand", Visible 0}
P_ = { "", Choices {"Fields", "Check_Periodicity"}, Visible 0,
Name "GetDP/2PostOperationChoices"}
];
Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
MoverPosition = mm*
DefineNumber[0, Name "Options/Mover position [mm]",
Visible !Flag_QuasiStatic];
TranslationStep = mm*
DefineNumber[1, Name "Options/Translation step [mm]",
Visible Flag_QuasiStatic];
NumStep =
DefineNumber[20, Name "Options/Number of steps",
Visible Flag_QuasiStatic];
velocity = 1 ; // m/s
Mag_DTime = TranslationStep/velocity ;
Mag_TimeMax = Mag_DTime*NumStep;
Group {
// Geometrical regions
StatorCore = Region[ 1 ] ;
MoverCore = Region[ 2 ] ;
AirStator = Region[ 3 ] ;
AirLayer = Region[ 4 ] ;
AirMover = Region[ 5 ] ;
CoilAP = Region[ 6 ] ;
CoilAN = Region[ 7 ] ;
CoilBP = Region[ 8 ] ;
CoilBN = Region[ 9 ] ;
CoilCP = Region[ 10 ] ;
CoilCN = Region[ 11 ] ;
NoFlux = Region[ 12 ] ;
Sur_MoverPerMaster = Region[ 13 ] ;
Sur_MoverPerSlave = Region[ 14 ] ;
Sur_StatorPerMaster = Region[ 15 ] ;
Sur_StatorPerSlave = Region[ 16 ] ;
Sur_SlidingSlave = Region[ 17 ] ;
Sur_SlidingMaster = Region[ 18 ] ;
NICEPOS = Region[ 19 ] ;
Lin_SlidingSubslave = Region[ 20 ] ;
Lin_SlidingSubmaster = Region[ 21 ] ;
// Abstract regions
Inductors = Region[ {CoilAP, CoilBP, CoilCP, CoilAN, CoilBN, CoilCN} ] ;
Iron = Region[ {StatorCore, MoverCore} ] ;
Air = Region[ {AirStator, AirMover, AirLayer} ] ;
Vol_Conductors = Region[ { Inductors } ] ;
Vol_Mag = Region[ { Air, Iron, Inductors } ] ;
Sur_Dirichlet = Region[ { NoFlux } ] ;
Vol_Moving = Region[ { MoverCore, AirMover, Inductors } ] ;
Sur_Link = Region[ {Sur_SlidingMaster, Sur_SlidingSlave,
Sur_StatorPerMaster, Sur_StatorPerSlave,
Sur_MoverPerMaster, Sur_MoverPerSlave} ] ;
Dom_Hcurl_a = Region[ {Vol_Mag, Sur_Dirichlet, Sur_Link} ] ;
}
Function {
mu0 = 4.e-7*Pi ;
murFer = 1000 ;
nu [ Air ] = 1. / mu0 ;
nu [ Inductors ] = 1. / mu0 ;
nu [ Iron ] = 1. / (murFer * mu0) ;
NbTurns = 120 ;
Iphase = 1;
phase = -60*deg;
js[ CoilAP ] = Vector[ 0,0, Iphase * Cos[phase + 0*deg] * NbTurns / SurfaceArea[]{6} ];
js[ CoilBP ] = Vector[ 0,0, Iphase * Cos[phase + 60*deg] * NbTurns / SurfaceArea[]{7} ];
js[ CoilCP ] = Vector[ 0,0, Iphase * Cos[phase + 120*deg] * NbTurns / SurfaceArea[]{8} ];
js[ CoilAN ] = Vector[ 0,0,-Iphase * Cos[phase + 0*deg] * NbTurns / SurfaceArea[]{6} ];
js[ CoilBN ] = Vector[ 0,0,-Iphase * Cos[phase + 60*deg] * NbTurns / SurfaceArea[]{7} ];
js[ CoilCN ] = Vector[ 0,0,-Iphase * Cos[phase + 120*deg] * NbTurns / SurfaceArea[]{8} ];
}
Function {
// Sliding surface:
// a1 is the span of Region and RegionRef.
// a0 is the mesh step of the meshes of the sliding surfaces
// a2[] is the position of the mover
// 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
TranslatePZ[] = Vector[ X[]+$1,Y[], Z[]] ;
Tol = 1e-8 ; fFloor[] = Floor[ $1 + Tol ] ; fRem[] = $1 - fFloor[ $1 / $2 ] * $2;
a1 = L ;
a0 = a1/NbrDiv ; // mesh stepsize on the sliding surface
a2[] = $MoverPosition ;
a3[] = ( ( fRem[ a2[], a0 ]#2 <= 0.5*a0 ) ? #2 : #2-a0 ) ;
AlignWithMaster[] = TranslatePZ[ - fFloor[ ( X[] - a3[] )/a1 ] * a1 ] ;
RestoreRef[] = XYZ[] - Vector[ 0, 0, SlidingSurfaceShift ] ;
Coef[] = ( fFloor[ ((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 ; } }
}
}
}
}
Constraint {
{ Name MagneticVectorPotential_2D ; Type Assign ;
Case {
{ Region Sur_Dirichlet ; Value 0. ; }
// Periodicity boundary condition on lateral faces
{ Type Link ; Region Sur_StatorPerSlave ; RegionRef Sur_StatorPerMaster ;
ToleranceFactor 1e-8;
Coefficient Periodicity ; Function TranslatePZ[ -a1 ] ; }
{ Type Link ; Region Sur_MoverPerSlave ; RegionRef Sur_MoverPerMaster ;
ToleranceFactor 1e-8;
Coefficient Periodicity ; Function TranslatePZ[ -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[] ; }
}
}
}
FunctionSpace {
{ Name Hcurl_a ; Type Form1P ;
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ;
Support Dom_Hcurl_a ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef ae ; EntityType NodesOf ; //EntitySubType Not ;
NameOfConstraint MagneticVectorPotential_2D ; }
}
}
}
Formulation {
{ Name Magnetostatics ; 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 { [ -js[] , {a} ] ; In Vol_Conductors ;
Jacobian Vol ; Integration Int ; }
}
}
}
Resolution {
{ Name Static ;
System {
{ Name Sys_Mag ; NameOfFormulation Magnetostatics ;}
}
Operation {
Evaluate[$MoverPosition = MoverPosition];
ChangeOfCoordinates[ NodesOf[ Vol_Moving ], TranslatePZ[ a2[] ] ] ;
ChangeOfCoordinates[ NodesOf[ Sur_SlidingMaster ], TranslatePZ[ a3[] ] ] ;
Evaluate[ $a2 = a2[] ];
Evaluate[ $a3 = a3[] ];
Evaluate[ $aa = Fmod[ a2[], a0 ] ];
Evaluate[ $bb = fRem[ a2[], a0 ] ];
Print[ {$MoverPosition, a0, $a3, $aa, $bb},
Format "wt=a2=%e a0=%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 Magnetostatics ; }
}
Operation {
Evaluate[$MoverPosition = MoverPosition];
Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
PostOperation [Fields];
TimeLoopTheta{
Time0 0 ; TimeMax Mag_TimeMax ; DTime Mag_DTime ; Theta 1. ;
Operation {
Evaluate[$MoverPosition = MoverPosition + velocity*$Time];
ChangeOfCoordinates[ NodesOf[ Vol_Moving ], TranslatePZ[ velocity*$DTime ] ] ;
ChangeOfCoordinates[ NodesOf[ Sur_SlidingMaster ], TranslatePZ[ a3[] ] ] ;
UpdateConstraint[Sys_Mag];
Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
PostOperation [Fields];
ChangeOfCoordinates[ NodesOf[ Sur_SlidingMaster ], TranslatePZ[ -a3[] ] ] ;
}
}
}
}
}
PostProcessing {
{ Name MagSta_a ; NameOfFormulation Magnetostatics ;
PostQuantity {
{ Name a ; Value { Term { [ {a} ] ;
In Vol_Mag ; Jacobian Vol ; } } }
{ Name az ; Value { Term { [ CompZ[{a}] ] ;
In Vol_Mag ; Jacobian Vol ; } } }
{ Name b ; Value { Term { [ {d a} ] ;
In Vol_Mag ; Jacobian Vol ; } } }
{ Name h ; Value { Term { [ nu[] * {d a} ] ;
In Vol_Mag ; } } }
{ Name bsurf ; Value {
Local { [ {d a} ]; In Sur_Link ; Jacobian Sur; } } }
{ Name un ; Value { Term { [ 1 ] ;
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].LineWidth = 2;",
"View[l].VectorType = 4;" ] ,
File "tmp.geo", LastTimeStepOnly] ;
*/
Print[ un, OnElementsOf Region[ NICEPOS ],
LastTimeStepOnly, File "un.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].LineWidth = 2;",
"View[l].ColormapNumber = 0;",
"View[l].ShowScale = 0;",
"Geometry.Lines = 0;",
"Geometry.Points = 0;"] ,
File "tmp.geo", LastTimeStepOnly] ;
Print[ az, OnElementsOf Region[ Vol_Mag ],
LastTimeStepOnly, File "a.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 1;",
"View[l].NbIso = 40;"] ,
File "tmp.geo", LastTimeStepOnly] ;
}
PostOperation Check_Periodicity UsingPost MagSta_a {
Print[ bsurf, OnElementsOf Region[ {Sur_SlidingMaster, Sur_SlidingSlave} ],
LastTimeStepOnly, File "bs.pos"] ;
Echo[ Str["l=PostProcessing.NbViews-1;",
"View[l].ArrowSizeMax = 100;",
"View[l].CenterGlyphs = 0;",
"View[l].GlyphLocation = 2;",
"View[l].VectorType = 2;" ] ,
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
// Geometrical parameters
wd = 40*mm;
ag = 2*mm; // 1*mm
// Rail-stator
ws = 15*mm;
cs = 15*mm;
vs = 21*mm;
hs = 10*mm;
Ps = 36*mm;
// Translateur - rotor
wt = 12*mm;
ct = 12*mm;
vt = 12*mm;
ht = 33*mm;
Pt = 24*mm;
// Bobines - windings
hw = 20*mm; // max 30*mm
vw = 4*mm; // max 5.5*mm
// Air box
hb = ag+ht+ct+20*mm;
NbrStaPer = 7;
L = NbrStaPer*(ws+vs);
Lm = 6*wt+5*vt;
ym = cs+hs+ag;
// Number of elements on the sliding surface in rotation direction
NbrDiv = L/lc;
SlidingSurfaceShift = 1.e-4;
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment