Skip to content
Snippets Groups Projects
Commit 430eee1e authored by Ruth Sabariego (U0089683)'s avatar Ruth Sabariego (U0089683)
Browse files

Homogenisation of windings: FD & TD

parent 09b04ede
No related branches found
No related tags found
No related merge requests found
Pipeline #5510 passed
Showing
with 2213 additions and 0 deletions
File added
Function{
// analytical
// analytical Brauer law for nonlinear isotropic material:
// nu(b^2) = k1 + k2 * exp ( k3 * b^2 )
// nu = 100. + 10. * exp ( 1.8 * b * b )
k1 = 100.; k2 = 10.; k3 = 1.8;
nu_1a[] = k1 + k2 * Exp[k3*SquNorm[$1]] ;
dnudb2_1a[] = k2 * k3 * Exp[k3*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 = Mat1_b()^2;
Mat1_nu = Mat1_h()/Mat1_b();
Mat1_nu(0) = Mat1_nu(1);
Mat1_nu_b2 = ListAlt[Mat1_b2(), Mat1_nu()] ;
nu_1[] = InterpolationLinear[ SquNorm[$1] ]{Mat1_nu_b2()} ;
dnudb2_1[] = dInterpolationLinear[SquNorm[$1]]{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 = Mat3kW_b()^2;
Mat3kW_nu = Mat3kW_h()/Mat3kW_b();
Mat3kW_nu(0) = Mat3kW_nu(1);
Mat3kW_nu_b2 = ListAlt[Mat3kW_b2(), Mat3kW_nu()] ;
nu_3kW[] = InterpolationLinear[SquNorm[$1]]{Mat3kW_nu_b2()} ;
dnudb2_3kW[] = dInterpolationLinear[SquNorm[$1]]{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] ;
DefineFunction[nu_lin, nu_nonlin, dhdb_NL] ;
// linear case -- testing purposes
h_lin[] = nu_lin[] * $1 ;
dhdb_lin[] = nu_lin[] * TensorDiag[1.,1.,1.] ;
dhdb_lin_NL[] = TensorDiag[0.,0.,0.] ;
}
Include "cell_dat.pro";
If(CellType==0)
Include "cell_square_circ.geo";
EndIf
If(CellType==1)
Include "cell_hexa_circ.geo";
EndIf
//getdp cell -setnumber Mode 1 -setnumber CellType 0 -sol MagDyn_a
Include "cell_dat.pro"
ResDir = "coeff/";
DefineConstant[
// Some plots.... Far from handy !!! :-(
// 1,2,3 -> Interval type
xx_ = "2000", // Active X (second pair corresponds to X')
yy_ = "0200", // Active Y (second pair corresponds to Y')
null_="0000", // Not active for that graph
// Rr = relative equivalent radius of conductor == reduced frequency
Rr = {2, Min 0.25, Max 4, Step 0.25, Loop 0, Name "Input/1Reduced frequency", Graph StrCat[xx_, xx_, xx_, xx_]}
Mode = {2, Choices{
1 = "skin",
2 = "proximity"}, //,3 = "proximity Y"}
Name "Input/Eddy current mode"}
plot_coefs = {0, Choices{0,1}, Name "Output/0Plot skin and proximity coefficients ?"}
];
ExtGmsh = ".pos" ;
ExtGnuplot = ".dat" ;
StringOut = Sprintf("_RS_la%.2g_1layer.dat", Fill) ; // Just two digits
//StringOut = ".dat" ;
Function {
//Sigma = 6e7;
Sigma = 5.67e7;
mu0 = 4.e-7 * Pi ; nu0 = 1./mu0;
delta = Rc/Rr;
Omega = 2/(delta*delta*mu0*Sigma);
Freq = Omega/(2*Pi); Period = 1./Freq;
Printf("delta %g mm ", delta*1000);
Printf("Frequency %g Hz Pulsation %g rad/s", Freq, Omega);
}
Group {
TheCond = Region[{COND}];
TheIsol = Region[{ISOL}];
TheCell = Region[{TheCond, TheIsol}];
Isol = Region[{@ ISOL:ISOL+NbrCond-1 @}];
Cond = Region[{@ COND:COND+NbrCond-1 @}];
Bound = Region[{BOUND}];
DomainC = Region[{Cond}];
DomainCC = Region[{Isol}];
Domain = Region[{DomainC, DomainCC}] ;
DomainDummy = Region[{1234}];
}
Function {
nu[] = nu0;
sigma[Cond] = Sigma;
sigma[TheIsol] = 0;
AreaTheCell[] = SurfaceArea[]{COND}+SurfaceArea[]{ISOL} ;
AreaTheCond[] = SurfaceArea[]{COND} ;
Rdc[] = 1. /AreaTheCond[] /Sigma ;
Req_[] = Sqrt[AreaTheCond[]/Pi];
SkinRef_i[] = mu0/(8*Pi) ;
ProxRef_[] = AreaTheCell[]*Sigma*Fill*(Req_[]*Omega)^2/4 ;
Constraint_a[] = (Mode==1) ? 0. :
((Mode==2) ? Vector[1,0,0] * Vector[$Y,-$X,0] : Vector[0,1,0] * Vector[$Y,-$X,0]) ;
Constraint_I[] = (Mode==1) ? 1. : 0. ;
}
Constraint {
{ Name MagneticVectorPotential_2D ;
Case {
{ Region Bound ; Value Constraint_a[] ; }
}
}
{ Name Current ;
Case {
{ Region DomainC ; Value Constraint_I[] ; }
}
}
{ Name Voltage ;
Case {
}
}
}
Jacobian {
{ Name Vol ; Case { { Region All ; Jacobian Vol ; } } }
}
Integration {
{ Name CurlCurl ; Case { { Type Gauss ; Case { { GeoElement Triangle ; NumberOfPoints 1 ; } } } } }
}
FunctionSpace {
{ Name Hcurl_a_2D ; Type Form1P ;
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ;
Support Domain ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef ae ; EntityType NodesOf ; NameOfConstraint MagneticVectorPotential_2D ; }
}
}
{ Name Hregion_u_2D ; Type Form1P ;
BasisFunction {
{ Name sr ; NameOfCoef ur ; Function BF_RegionZ ;
Support DomainC; Entity DomainC; }
}
GlobalQuantity {
{ Name U ; Type AliasOf ; NameOfCoef ur ; }
{ Name I ; Type AssociatedWith ; NameOfCoef ur ; }
}
Constraint {
{ NameOfCoef U ; EntityType Region ; NameOfConstraint Voltage ; }
{ NameOfCoef I ; EntityType Region ; NameOfConstraint Current ; }
}
}
}
Formulation {
{ Name MagDyn_a ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
{ Name ur ; Type Local ; NameOfSpace Hregion_u_2D ; }
{ Name I ; Type Global ; NameOfSpace Hregion_u_2D [I] ; }
{ Name U ; Type Global ; NameOfSpace Hregion_u_2D [U] ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain ; Jacobian Vol ; Integration CurlCurl ; }
Galerkin { DtDof [ sigma[] * Dof{a} , {a} ] ;
In DomainC ; Jacobian Vol ; Integration CurlCurl ; }
Galerkin { [ sigma[] * Dof{ur} , {a} ] ;
In DomainC ; Jacobian Vol ; Integration CurlCurl ; }
Galerkin { DtDof [ sigma[] * Dof{a} , {ur} ] ;
In DomainC ; Jacobian Vol ; Integration CurlCurl ; }
Galerkin { [ sigma[] * Dof{ur} , {ur} ] ;
In DomainC ; Jacobian Vol ; Integration CurlCurl ; }
GlobalTerm { [ Dof{I} , {U} ] ; In DomainC ; }
}
}
}
Resolution {
{ Name MagDyn_a ;
System {
{ Name A ; NameOfFormulation MagDyn_a ; Frequency Freq ; }
}
Operation {
CreateDir[ResDir];
SetTime[Rr];
Generate[A]; Solve[A]; SaveSolution[A];
PostOperation[Map_local] ;
PostOperation[Get_coeffs] ;
}
}
}
PostProcessing {
{ Name MagDyn_a ; NameOfFormulation MagDyn_a ;
PostQuantity {
{ Name a ; Value { Term { [ CompZ[{a}] ] ; In Domain ; Jacobian Vol ;} } }
{ Name b ; Value { Term { [ {d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name bav ; Value { Integral { [{d a}/AreaTheCell[]] ;
In Domain ; Jacobian Vol ; Integration CurlCurl ; } } }
// Skin-effect coefficients
{ Name pI ; Value { Integral { [ SquNorm[(Sigma*(Dt[{a}]+{ur}))]*AreaTheCond[] ] ;
In TheCell ; Jacobian Vol ; Integration CurlCurl ; } } }
{ Name qI ; Value { Integral { [ nu0 * SquNorm[{d a}]/SkinRef_i[] ] ;
In TheCell ; Jacobian Vol ; Integration CurlCurl ; } } }
// Proximity-effect coefficients
{ Name qB ; Value { Integral { [ SquNorm[{d a}]/AreaTheCell[] ] ;
In TheCell ; Jacobian Vol ; Integration CurlCurl ; } } }
{ Name pB ; Value { Integral { [ SquNorm[sigma[]*(Dt[{a}]+{ur})]/Sigma/ProxRef_[] ] ;
In TheCond ; Jacobian Vol ; Integration CurlCurl ; } } }
{ Name nuRe ; Value { Term { Type Global; [ nu0*$qB ] ; In DomainDummy ; } } }
{ Name nuIm ; Value { Term { Type Global; [ nu0*$pB*Fill*Rr^2/2 ] ; In DomainDummy ; } } }
// Normalized by nu0 or mu0
{ Name muRe ; Value { Term { Type Global; [ Re[1/Complex[$qB, $pB*Fill*Rr^2/2] ]] ; In DomainDummy ; } } }
{ Name muIm ; Value { Term { Type Global; [ Im[1/Complex[$qB, $pB*Fill*Rr^2/2] ]] ; In DomainDummy ; } } }
}
}
}
PostOperation Map_local UsingPost MagDyn_a {
Print[ b, OnElementsOf Domain, File StrCat[ResDir, "b",ExtGmsh]] ;
Print[ a, OnElementsOf Domain, File StrCat[ResDir, "a",ExtGmsh]] ;
Echo[ Str["i=PostProcessing.NbViews-1;
View[i].Light=0;
View[i].LineWidth = 2;
View[i].RangeType=3;
View[i].IntervalsType=3;
View[i].NbIso = 25;"], File StrCat[ResDir,"option_cell.pos"] ];
}
PostOperation Get_coeffs UsingPost MagDyn_a {
If(Mode==1)
Print[ pI[TheCond], OnGlobal, Format TimeTable, File > StrCat[ResDir, "pI", StringOut],
SendToServer "Output/0Skin/pI", StoreInVariable $skin_pI ] ;
Print[ qI[TheCell], OnGlobal, Format TimeTable, File > StrCat[ResDir, "qI", StringOut],
SendToServer "Output/0Skin/qI", StoreInVariable $skin_qI ] ;
EndIf
If(Mode==2)
Print[ qB[TheCell], OnGlobal, Format TimeTable, File > StrCat[ResDir, "qB", StringOut],
SendToServer "Output/1Prox/qB", StoreInVariable $qB] ;
Print[ pB[TheCond], OnGlobal, Format TimeTable, File > StrCat[ResDir, "pB", StringOut],
SendToServer "Output/1Prox/pB", StoreInVariable $pB ] ;
Print[ muRe, OnRegion DomainDummy, Format Table, LastTimeStepOnly,
SendToServer "Output/2Prox/Re(mu)", File > StrCat[ResDir, "muRe", StringOut]];
Print[ muIm, OnRegion DomainDummy, Format Table, LastTimeStepOnly,
SendToServer "Output/2Prox/Im(mu)", File > StrCat[ResDir, "muIm", StringOut]];
EndIf
If(Mode==3)
Print[ qB[TheCell], OnGlobal, Format TimeTable, File > StrCat[ResDir, "qBy", StringOut],
SendToServer "Output/Prox/qBy", StoreInVariable $qBy ] ;
Print[ pB[TheCond], OnGlobal, Format TimeTable, File > StrCat[ResDir, "pBy", StringOut],
SendToServer "Output/Prox/pBy", StoreInVariable $pBy ] ;
EndIf
}
DefineConstant[
R_ = {"MagDyn_a", Name "GetDP/1ResolutionChoices", Visible 1, Closed 1},
C_ = {"-solve -v 4 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 1, Closed 1}
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 1,Closed 1}
];
If(plot_coefs)
DefineConstant[
pI_= {1, Name "Output/0Skin/pI", Graph StrCat[yy_,null_,null_,null_], Visible (Mode==1)}
qI_= {1, Name "Output/0Skin/qI", Graph StrCat[null_,yy_,null_,null_], Visible (Mode==1)}
qB_= {1, Name "Output/1Prox/qB", Graph StrCat[null_,null_,yy_,null_], Visible (Mode==2)}
pB_= {1, Name "Output/1Prox/pB", Graph StrCat[null_,null_,null_,yy_], Visible (Mode==2)}
];
EndIf
mm=1e-3;
r_sh = 9.5*1e-2;
d_sh = 0.2;
DefineConstant[
CellType = {0, Choices{0="Square",1="Hexa"}, Name "Cell/00Type of Cell"}
NbrLayers = {0, Min 0, Max 1, Name "Cell/01Number of layers"}
//Rc = { (CellType==0) ? Sqrt[1/Pi] : Sqrt[1/Pi]/Sqrt[3],
Rc = { (CellType==0) ? r_sh : Sqrt[1/Pi]/Sqrt[3],
Name "Cell/10Conductor radius (mm)"} //Cefc06
// Dex = { 2.2*Rc, Min 2.01*Rc, Max 4*Rc, Name "Cell/11Packing side (mm)", Visible (CellType==0)}
Dex = { d_sh, Min 2.01*Rc, Max 4*Rc, Name "Cell/11Packing side (mm)", Visible (CellType==0)}
];
If(CellType==0)
UndefineConstant["Cell/20Fill factor"];
UndefineConstant["Cell/11Outer radius hexagonal cell (mm)"];
DefineConstant[
Fill = {Pi*Rc^2/(Dex*Dex), Name "Cell/20Fill factor", ReadOnly 1, Highlight "LightGrey"}
];
Dex = Dex*mm;
Dey = Dex;
EndIf
If(CellType==1)
UndefineConstant["Cell/20Fill factor"];
DefineConstant[
Fill = {0.9, Name "Cell/20Fill factor", ReadOnly 0}
Rx = {Rc/Sin[Pi/3]/Fill, Name "Cell/11Outer radius hexagonal cell (mm)",
ReadOnly 1, Visible CellType==1, Highlight "LightGrey"}
];
Rx = Rx*mm;
Ry = Rx;
EndIf
Rc = Rc*mm;
AreaCond = Pi*Rc^2;
AreaCell = AreaCond/Fill;
NbrCond = (NbrLayers==0) ? 1 : ((CellType==0)?9:7) ;
// Printf("Round conductor, radius = %g mm", Rc*1e3);
// Printf("Square packing Fill factor = %g mm^2 / %g mm^2 = %g ",AreaCond*1e6,AreaCell*1e6,Fill);
If (Rc > Dex/2)
Printf("Square packing: Two big fill factor!!! Rc %g > Dex/2 %g",Rc*1e6,Dex/2*1e6);
Printf("limit Fill factor %g",AreaCond/(4*Rc*Rc));
Dex = 2*(Rc+Rc*1e-2) ; Dey = Dex ;
EndIf
// characteristic lengths
// ======================
DefineConstant[
MD = {1/1.5, Name "Cell/99Mesh density factor"}
];
p = Rc/10*MD ;
pc = Rc/10*MD ;
ISOL = 200;
COND = 1000;
BOUND = 2000;
//Include "cell_dat.pro";
//p = Rc/3;
//pc = Rc/3;
Dx = 0.; Dy = 0.;
Geometry.AutoCoherence = 0;
Mesh.Algorithm = 1;
//-----------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------
Function HexCirc_
dP=newp-1;
dR=news-1;
Point(dP+1) = {Dx+Rx*Cos[1*Pi/6], Dy+Ry*Sin[1*Pi/6], 0 , p};
Point(dP+2) = {Dx+Rx*Cos[3*Pi/6], Dy+Ry*Sin[3*Pi/6], 0 , p};
Point(dP+3) = {Dx+Rx*Cos[5*Pi/6], Dy+Ry*Sin[5*Pi/6], 0 , p};
Point(dP+4) = {Dx+Rx*Cos[7*Pi/6], Dy+Ry*Sin[7*Pi/6], 0 , p};
Point(dP+5) = {Dx+Rx*Cos[9*Pi/6], Dy+Ry*Sin[9*Pi/6], 0 , p};
Point(dP+6) = {Dx+Rx*Cos[11*Pi/6], Dy+Ry*Sin[11*Pi/6], 0 , p};
Point(dP+7) = {Dx, Dy, 0 , p};
Line(dR+1) = {dP+1,dP+2};
Line(dR+2) = {dP+2,dP+3};
Line(dR+3) = {dP+3,dP+4};
Line(dR+4) = {dP+4,dP+5};
Line(dR+5) = {dP+5,dP+6};
Line(dR+6) = {dP+6,dP+1};
Line Loop(dR+1) = {dR+6,dR+1,dR+2,dR+3,dR+4,dR+5};
Point(dP+8) = {Dx, Dy, 0 , pc};
Point(dP+9) = {Dx+Rc, Dy, 0 , pc};
Point(dP+10) = {Dx-Rc, Dy, 0 , pc};
Point(dP+11) = {Dx, Dy+Rc, 0 , pc};
Point(dP+12) = {Dx, Dy-Rc, 0 , pc};
Circle(dR+7) = {dP+9,dP+8,dP+11};
Circle(dR+8) = {dP+11,dP+8,dP+10};
Circle(dR+9) = {dP+10,dP+8,dP+12};
Circle(dR+10) = {dP+12,dP+8,dP+9};
Line Loop(dR+2) = {dR+7,dR+8,dR+9,dR+10};
Plane Surface(dR+3) = {dR+1,dR+2};
Plane Surface(dR+4) = {dR+2};
surfCond[] += {dR+4};
surfIsol[] += {dR+3};
Return
//=========================================
surfCond[] = {};
surfIsol[] = {};
Call HexCirc_ ;
If (NbrLayers > 0)
dx1 = 2*Rx*Cos[Pi/6] ;
dy1 = 2*Ry*Cos[Pi/6] ;
iL = 0 ; kk = 0 ; l = 0 ;
For iL In {1:NbrLayers}
Dx = -iL * dx1 ;
Dy = 0;
For k In {1:6}
For l In {1:iL}
kk = 0 ;
If (iL == NbrLayers)
kk = k ;
EndIf
// COND += 1 ;
// ISOL += 1 ;
Dx += dx1 * Cos[(2-k)*Pi/3] ;
Dy += dy1 * Sin[(2-k)*Pi/3] ;
Call HexCirc_ ;
EndFor
EndFor
EndFor
EndIf
Geometry.AutoCoherence = 1;
Coherence;
For k In {0:#surfCond[]-1}
Physical Surface(COND+k) = {surfCond[k]};
Physical Surface(ISOL+k) = {surfIsol[k]};
EndFor
allSurfaces[] = Surface '*' ;
lines_boundary[] = CombinedBoundary{Surface{allSurfaces[]}; };
Physical Line(BOUND) = lines_boundary[];
//Include "cell_square_circ_dat.pro";
cen0 = newp ; Point(cen0) = { 0, 0, 0 , pc};
pw[]+=newp ; Point(pw[0]) = { Rc, 0, 0 , pc};
pw[]+=newp ; Point(pw[1]) = { 0, Rc, 0 , pc};
pw[]+=newp ; Point(pw[2]) = {-Rc, 0, 0 , pc};
pw[]+=newp ; Point(pw[3]) = { 0,-Rc, 0 , pc};
cw[]+=newl ; Circle(cw[0]) = {pw[0], cen0, pw[1]};
cw[]+=newl ; Circle(cw[1]) = {pw[1], cen0, pw[2]};
cw[]+=newl ; Circle(cw[2]) = {pw[2], cen0, pw[3]};
cw[]+=newl ; Circle(cw[3]) = {pw[3], cen0, pw[0]};
llcw[] += newll ; Line Loop(llcw[0]) = {cw[]};
sw[] += news ; Plane Surface(sw[0]) = {llcw[0]};
pb[]+=newp ; Point(pb[0]) = { Dex/2, Dey/2, 0 , p};
pb[]+=newp ; Point(pb[1]) = {-Dex/2, Dey/2, 0 , p};
pb[]+=newp ; Point(pb[2]) = {-Dex/2, -Dey/2, 0 , p};
pb[]+=newp ; Point(pb[3]) = { Dex/2, -Dey/2, 0 , p};
lb[]+=newl ; Line(lb[0]) = {pb[0], pb[1]};
lb[]+=newl ; Line(lb[1]) = {pb[1], pb[2]};
lb[]+=newl ; Line(lb[2]) = {pb[2], pb[3]};
lb[]+=newl ; Line(lb[3]) = {pb[3], pb[0]};
llb[] += newll ; Line Loop(llb[0]) = {lb[]};
is[] += news ; Plane Surface(is[0]) = {llb[0],llcw[0]};
If(NbrLayers==1)
xaux[] = {0, 0, Dex, Dex, Dex, -Dex, -Dex, -Dex};
yaux[] = {Dey, -Dey, Dey, 0, -Dey, Dey, 0, -Dey};
For k In {0:7}
surf[]=Translate {xaux[k], yaux[k], 0} {Duplicata{ Surface{sw[0]}; }}; sw[] += surf[0] ;
surf[]=Translate {xaux[k], yaux[k], 0} {Duplicata{ Surface{is[0]}; }}; is[] += surf[0] ;
EndFor
EndIf
For k In {0:#sw[]-1}
Physical Surface(COND+k) = {sw[k]};
Physical Surface(ISOL+k) = {is[k]};
EndFor
bndi[] =CombinedBoundary{Surface{is[],sw[]};};
Physical Line(BOUND) = {bndi[]};
File added
Function {
// Coefficients for homogenization in time domain
// with initial guess Y = zeros(1,2*Ord-1); Y(1,1) = 1; % ++++ original 27.01.2017
// Proximity effect (matrix P)
Ca_Sq_1_1 = 1. ;
Cb_Sq_1_1 = 0.9966639477003278;
Ca_Sq_2_1 = 1. ;
Ca_Sq_2_2 = 1. ;
Cb_Sq_2_1 = 0.991889311702569;
Cb_Sq_2_2 = 0.5582220430389034;
Cc_Sq_2_1 = 0.7036557691449754;
Ca_Sq_3_1 = 1. ;
Ca_Sq_3_2 = 1. ;
Ca_Sq_3_3 = 1. ;
Cb_Sq_3_1 = 0.9970530831671187;
Cb_Sq_3_2 = 0.5910488876830849;
Cb_Sq_3_3 = 0.1388368497131864;
Cc_Sq_3_1 = 0.7212112634042309;
Cc_Sq_3_2 = 0.09021468658196302;
Ca_Sq_4_1 = 1. ;
Ca_Sq_4_2 = 1. ;
Ca_Sq_4_3 = 1. ;
Ca_Sq_4_4 = 1. ;
Cb_Sq_4_1 = 0.9977778378722477;
Cb_Sq_4_2 = 0.5964651655582043;
Cb_Sq_4_3 = 0.1565496755914985;
Cb_Sq_4_4 = 0.0598158454879603;
Cc_Sq_4_1 = 0.7235590940728726;
Cc_Sq_4_2 = 0.09610768489856511;
Cc_Sq_4_3 = 0.03640901775369249;
// Skin effect (matrix S)
Sa_Sq_1_1 = 1. ;
Sb_Sq_1_1 = 2.758272999411088;
Sa_Sq_2_1 = 1. ;
Sa_Sq_2_2 = 1. ;
Sb_Sq_2_1 = 2.525106665179376;
Sb_Sq_2_2 = 0.6873975017935905;
Sc_Sq_2_1 = 0.8564741317336022;
Sa_Sq_3_1 = 1. ;
Sa_Sq_3_2 = 1. ;
Sa_Sq_3_3 = 1. ;
Sb_Sq_3_1 = 2.523956861106585;
Sb_Sq_3_2 = 0.8230227201355227;
Sb_Sq_3_3 = 0.2709304248426088;
Sc_Sq_3_1 = 0.9182551475941608;
Sc_Sq_3_2 = 0.2264149782781768;
Sa_Sq_4_1 = 1. ;
Sa_Sq_4_2 = 1. ;
Sa_Sq_4_3 = 1. ;
Sa_Sq_4_4 = 1. ;
Sb_Sq_4_1 = 2.50157222488763;
Sb_Sq_4_2 = 0.9527222452515992;
Sb_Sq_4_3 = 0.5445711604362855;
Sb_Sq_4_4 = 0.1673917470007547;
Sc_Sq_4_1 = 0.9407809031805078;
Sc_Sq_4_2 = 0.3647042320179443;
Sc_Sq_4_3 = 0.1504371848184149;
/*
// old values
Sa_Sq_4_1 = 1. ;
Sa_Sq_4_2 = 1. ;
Sa_Sq_4_3 = 1. ;
Sa_Sq_4_4 = 1. ;
Sb_Sq_4_1 = 1.892116674831446;
Sb_Sq_4_2 = 1.388437802179596e-15;
Sb_Sq_4_3 = 4.254427483363831e-16;
Sb_Sq_4_4 = 0.2244497401308106;
Sc_Sq_4_1 = 1.439947715782265;
Sc_Sq_4_2 = 2.088756751297725;
Sc_Sq_4_3 = 0.5862977592702109;
*/
}
0 1
0.25 0.9984797102123941
0.5 0.9965504583511245
0.75 0.9882931098705841
1 0.9668596338010401
1.25 0.9251466324415509
1.5 0.8594785554811457
1.75 0.7729776511320164
2 0.6752841987717687
2.25 0.5780646223656311
2.5 0.4900812379090719
2.75 0.4153193856063418
3 0.3539783124642665
3.25 0.3043201994523901
3.5 0.2640729253772856
3.75 0.2311438203598003
4 0.2038563490018324
4.25 0.1809516029837255
4.5 0.1615069261034517
4.75 0.1448447352371066
5 0.1304581813885478
5.25 0.1179582785716005
5.5 0.1070388417605652
5.75 0.09745385145001674
6 0.08900269179423495
6.25 0.08152013011678975
6.5 0.0748691305811709
6.75 0.0689354289872187
7 0.06362328727611276
7.25 0.05885210633748061
7.5 0.05455370261205555
7.75 0.05067011406056374
8 0.04715183174770684
0 1
0.25 1.000292198465451
0.5 1.004632454799866
0.75 1.022576730096728
1 1.06518934015796
1.25 1.137034415586967
1.5 1.234811391838964
1.75 1.354085103502352
2 1.49189675132252
2.25 1.644618526120431
2.5 1.807143578637712
2.75 1.974241214763512
3 2.142111983254704
3.25 2.308913911935778
3.5 2.474345701972097
3.75 2.638906450909167
4 2.803290943979866
4.25 2.968061231919891
4.5 3.133545545125549
4.75 3.299865123204656
5 3.467007163270016
5.25 3.634895946649011
5.5 3.803442482092406
5.75 3.9725705179669
6 4.142224829905706
6.25 4.312369441796784
6.5 4.482981978531287
6.75 4.654047928036489
7 4.825556459611321
7.25 4.997498058931526
7.5 5.169863573115075
7.75 5.342644091497067
8 5.515831180746527
0 1.0
0.25 1.000215650342192
0.5 1.003442630297233
0.75 1.017259852487556
1 1.053167879435913
1.25 1.123242149948414
1.5 1.234156164920968
1.75 1.381695699352034
2 1.551156898427067
2.25 1.724524664948355
2.5 1.888326196664408
2.75 2.036534635416051
3 2.168876707555387
3.25 2.287815907668703
3.5 2.39632879694952
3.75 2.496855372456543
4 2.591052376804712
4.25 2.67992289322597
4.5 2.764046966615837
4.75 2.843782093019773
5 2.919393816085748
5.25 2.991122211801619
5.5 3.059205528411236
5.75 3.123881769969344
6 3.185382555451071
6.25 3.243926810658909
6.5 3.299717029209635
6.75 3.352938227562404
7 3.403758728985339
7.25 3.452331823645406
7.5 3.498797623273519
7.75 3.543284748277366
8 3.585911723092616
0 1
0.25 2.766777056507074
0.5 2.757741498450176
0.75 2.721105328878783
1 2.638935350011856
1.25 2.516393262268829
1.5 2.380491290104168
1.75 2.253031030379621
2 2.139558796705369
2.25 2.037924521862575
2.5 1.945881796582725
2.75 1.862692915358307
3 1.788259116208388
3.25 1.72230282621608
3.5 1.664152885451796
3.75 1.61288056884093
4 1.567498230260288
4.25 1.527094270354651
4.5 1.49089116506077
4.75 1.458252339865828
5 1.428664696829399
5.25 1.401714151139559
5.5 1.377062788349206
5.75 1.354430674525283
6 1.333582490470254
6.25 1.314318043063593
6.5 1.2964655167714
6.75 1.27987652846567
7 1.26442233974431
7.25 1.249990829857953
7.5 1.236484002507318
7.75 1.223815899448496
8 1.211910845038602
Include "ind_axi_dat.pro";
Mesh.Algorithm = 1; // 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=BAMG, 8=DelQuad)
p=.07e-3 *MD;
pc = Rc/8 *MD;
pr = 0.3e-3*1.5 *MD;
pr2 = 1e-3 *MD;
pr3 = 4e-3 *MD;
pax = 0.5e-3 *MD;
If(Fine)
prW = pr*0.7 ;
Else
prW = MDH*pr*0.7 ;
EndIf
dPo[]+=newp ; Point(dPo[0]) = {X1, Y2, 0 , pax};
dPo[]+=newp ; Point(dPo[1]) = {X2, Y2, 0 , pr2};
dPo[]+=newp ; Point(dPo[2]) = {X3, Y2, 0 , pr2};
dPo[]+=newp ; Point(dPo[3]) = {X4, Y2, 0 , pr2};
dPo[]+=newp ; Point(dPo[4]) = {X4, Y3, 0 , pr2};
dPo[]+=newp ; Point(dPo[5]) = {X4, Y4, 0 , pr};
dPo[]+=newp ; Point(dPo[6]) = {X4, 0, 0 , pr};
dPo[]+=newp ; Point(dPo[7]) = {X3, 0, 0 , pr};
dPo[]+=newp ; Point(dPo[8]) = {X3, Y4, 0 , pr};
dPo[]+=newp ; Point(dPo[9]) = {X3, Y3, 0 , pr};
dPo[]+=newp ; Point(dPo[10]) = {X2, Y3, 0 , pr};
dPo[]+=newp ; Point(dPo[11]) = {X2, Y4, 0 , pr*0.7};
dPo[]+=newp ; Point(dPo[12]) = {X1, Y4, 0 , pax};
dPo[]+=newp ; Point(dPo[13]) = {X1, Y3, 0 , pax};
dPg[]+=newp ; Point(dPg[0]) = {X2, 0, 0 , pr*0.7};
dPg[]+=newp ; Point(dPg[1]) = {X1, 0, 0 , pax};
For k1 In {0:#dPo[]-1}
k2 = (k1 != #dPo[]-1) ? k1+1 : 0 ;
dLo[]+=newl ; Line(dLo[k1]) = {dPo[k1], dPo[k2]};
EndFor
lliron=newll ; Line Loop(lliron) = {dLo[]};
surfiron[] += news ; Plane Surface(surfiron[0]) = {lliron};
dLg[]+=newl ; Line(dLg[0]) = {dPo[11],dPg[0]};
dLg[]+=newl ; Line(dLg[1]) = {dPg[0],dPg[1]};
dLg[]+=newl ; Line(dLg[2]) = {dPg[1],dPo[12]};
llgap1 = newll ;
Line Loop(llgap1) = {dLg[], -dLo[11]};
surfgap1[] += news ; Plane Surface(surfgap1[0]) = {llgap1};
dPw[]+=newp ; Point(dPw[0]) = {Xw1, 0*Yw1, 0., prW};
dPw[]+=newp ; Point(dPw[1]) = {Xw1, Yw2, 0., prW};
dPw[]+=newp ; Point(dPw[2]) = {Xw2, Yw2, 0., prW};
dPw[]+=newp ; Point(dPw[3]) = {Xw2, 0*Yw1, 0., prW};
dLw[]+=newl ; Line(dLw[0]) = {dPg[0], dPw[0]};
For k1 In {0:#dPw[]-2}
k2 = (k1 != #dPw[]-1) ? k1+1 : 0 ;
dLw[]+=newl ; Line(dLw[1+k1]) = {dPw[k1], dPw[k2]};
EndFor
dLw[]+=newl ; Line(dLw[4]) = {dPw[3],dPo[7]};
llgap2 = newll ; Line Loop(llgap2) = {dLg[0],dLw[],dLo[{7:10}]};
surfgap2[] += news ; Plane Surface(surfgap2[0]) = {llgap2};
surfhomo[]={}; sw[]={}; is[]={};
If (Fine)
// Auxiliar elementary cell
cen0 = newp ; Point(cen0) = { 0, 0, 0 , pc};
pw[]+=newp ; Point(pw[0]) = { Rc, 0, 0 , pc};
pw[]+=newp ; Point(pw[1]) = { 0, Rc, 0 , pc};
pw[]+=newp ; Point(pw[2]) = {-Rc, 0, 0 , pc};
pw[]+=newp ; Point(pw[3]) = { 0,-Rc, 0 , pc};
cw[]+=newl ; Circle(cw[0]) = {pw[0], cen0, pw[1]};
cw[]+=newl ; Circle(cw[1]) = {pw[1], cen0, pw[2]};
cw[]+=newl ; Circle(cw[2]) = {pw[2], cen0, pw[3]};
cw[]+=newl ; Circle(cw[3]) = {pw[3], cen0, pw[0]};
llcw0 = newll ; Line Loop(llcw0) = {cw[]};
sw0 = news ; Plane Surface(sw0) = {llcw0};
pb[]+=newp ; Point(pb[0]) = { Dex/2, Dey/2, 0 , p};
pb[]+=newp ; Point(pb[1]) = {-Dex/2, Dey/2, 0 , p};
pb[]+=newp ; Point(pb[2]) = {-Dex/2, -Dey/2, 0 , p};
pb[]+=newp ; Point(pb[3]) = { Dex/2, -Dey/2, 0 , p};
lb[]+=newl ; Line(lb[0]) = {pb[0], pb[1]};
lb[]+=newl ; Line(lb[1]) = {pb[1], pb[2]};
lb[]+=newl ; Line(lb[2]) = {pb[2], pb[3]};
lb[]+=newl ; Line(lb[3]) = {pb[3], pb[0]};
llb0 = newll ; Line Loop(llb0) = {lb[]};
is0 = news ; Plane Surface(is0) = {llb0,llcw0};
gx = (Xw2-Xw1)-NbrLayersX * Dex ;
gy = 0.5*((Yw2-Yw1)-NbrLayersY * Dey) ;
NbrLayersY2 = 0.5*NbrLayersY;
For ix In {0:NbrLayersX-1}
For iy In {0:NbrLayersY2-1}
xaux = Xw1+gx/2+ix*Dex+Dex/2 ;
yaux = 0*(Yw1+gy/2)+iy*Dey+Dey/2 ;
surf[]=Translate {xaux, yaux, 0} {Duplicata{ Surface{sw0};}}; sw[] += surf[0] ;
surf[]=Translate {xaux, yaux, 0} {Duplicata{ Surface{is0};}}; is[] += surf[0] ;
If(ix==0)
bnd[] = Boundary{Surface{is[ix*NbrLayersY2+iy]};}; bndi[] += bnd[{1}];
EndIf
If(ix==NbrLayersX-1)
bnd[] = Boundary{Surface{is[ix*NbrLayersY2+iy]};}; bndi[] += bnd[{3}];
EndIf
If(iy==0)
bnd[] = Boundary{Surface{is[ix*NbrLayersY2+iy]};}; axx[] += bnd[{2}];
EndIf
If(iy==NbrLayersY2-1)
bnd[] = Boundary{Surface{is[ix*NbrLayersY2+iy]};}; bndi[] += bnd[{0}];
EndIf
EndFor
EndFor
//Removing auxiliar elementary cell
Delete{Surface{sw0,is0};Line{cw[],lb[]};Point{pw[],pb[]};}
bnd[] = Boundary{Line{axx[0]};}; dPw[]+= bnd[0];
bnd[] = Boundary{Line{axx[{#axx[]-1}]};}; dPw[]+= bnd[1];
dLw2[]+=newl ; Line(dLw2[0]) = {dPw[4], dPw[0]};
dLw2[]+=newl ; Line(dLw2[1]) = {dPw[3], dPw[5]};
llisol = newll ; Line Loop(llisol) = {bndi[],dLw2[0], dLw[{1:3}],dLw2[1]};
surfgap3[] += news ; Plane Surface(surfgap3[0]) = {llisol};
EndIf
If(!Fine)
dPis[]+=newp ; Point(dPis[0]) = {Xw1_, Yw1_, 0 , prW};
dPis[]+=newp ; Point(dPis[1]) = {Xw1_, Yw2_, 0 , prW};
dPis[]+=newp ; Point(dPis[2]) = {Xw2_, Yw2_, 0 , prW};
dPis[]+=newp ; Point(dPis[3]) = {Xw2_, Yw1_, 0 , prW};
For k1 In {0:#dPw[]-1}
k2 = (k1 != #dPw[]-1) ? k1+1 : 0 ;
bndi[]+=newl ; Line(bndi[k1]) = {dPis[k1], dPis[k2]};
EndFor
llhomo = newll ; Line Loop(llhomo) = {bndi[]};
surfhomo[] += news ; Plane Surface(surfhomo[0]) = {llhomo};
dLw2[]+=newl ; Line(dLw2[0]) = {dPis[0], dPw[0]};
dLw2[]+=newl ; Line(dLw2[1]) = {dPw[3], dPis[3]};
llisol = newll ; Line Loop(llisol) = {bndi[{0:2}], -dLw2[1], -dLw[{3:1}],-dLw2[0]};
surfgap3[] += news ; Plane Surface(surfgap3[0]) = {llisol};
EndIf
line_out[] = {dLo[{0:5,12:13}],dLg[2]};
If(!Flag_HalfModel)
line_out[]+=Symmetry {0,1,0,0} { Duplicata{Line{line_out[]};} }; // For convenience
surfiron[]+=Symmetry {0,1,0,0} { Duplicata{Surface{surfiron[0]};} };
surfgap1[]+=Symmetry {0,1,0,0} { Duplicata{Surface{surfgap1[0]};} };
surfgap2[]+=Symmetry {0,1,0,0} { Duplicata{Surface{surfgap2[0]};} };
surfgap3[]+=Symmetry {0,1,0,0} { Duplicata{Surface{surfgap3[0]};} };
If(Fine)
nn = #sw[];
For k In {0:nn-1}
sw[] += Symmetry {0,1,0,0} { Duplicata{Surface{sw[k]};} };
EndFor
is[] += Symmetry {0,1,0,0} { Duplicata{Surface{is[]};} };
Else
surfhomo[]+=Symmetry {0,1,0,0} { Duplicata{Surface{surfhomo[0]};} };
EndIf
EndIf
// Physical regions
//===================================
Physical Line("Outer boundary", OUTBND) = {line_out[]};
Physical Surface("Iron core", IRON) = {surfiron[]};
Physical Surface("Airgap", AIRGAP) = {surfgap1[]};
Physical Surface("Air",AIR) = {surfgap2[],surfgap3[]};
Reverse Surface {surfiron[],surfgap1[],surfgap3[]}; // Inverting normals
If(Fine)
nn = #sw[];
For k In {0:nn-1}
Physical Surface(Sprintf("Cond %g", k),iCOND+k) = {sw[k]};
EndFor
Physical Surface("Insulation", INSULATION) = {is[]};
Physical Surface("All conductors", ALLCOND) = {sw[]};
Else
Physical Surface("Winding window", iCOND) = {surfhomo[]};
EndIf
Mesh.Light = 0;
Color SteelBlue{
Surface{surfiron[]};
}
Color Red{
Surface{sw[], surfhomo[]};
}
Color Gold{
Surface{is[]};
}
Color SkyBlue{
Surface{surfgap1[],surfgap2[],surfgap3[]};
}
// Time domain + fine model:
// gmsh ind_axi.geo -setnumber Flag_HomogenisedModel 0 -2 -o fine.msh -v 3
// getdp ind_axi -msh fine.msh -setnumber Flag_HomogenisedModel 0 -sn Flag_FD 0 -sn Flag_SrcType 2 -sn NbT 1 -sol Analysis
// Resistance (Ohms): p 'SF_s0.dat' u 1:2 w l lw 2, 'SH.dat' u 1:2 w p ps 2 pt 6 lw 2
// Inductance (mH): p 'SF_s0.dat' u 1:($3*1e3) w l lw 2, 'SH.dat' u 1:($3*1e3) w p ps 2 pt 6 lw 2
Include "ind_axi_dat.pro"
Include "BH.pro"
la = Fill; // fill factor
// new files, finer mesh, max X=8
file_ZSkinRe = Sprintf("coeff/pI_RS_la%.2g.dat", la);
file_ZSkinIm = Sprintf("coeff/qI_RS_la%.2g.dat", la);
file_NuProxRe = Sprintf("coeff/qB_RS_la%.2g.dat", la);
file_NuProxIm = Sprintf("coeff/pB_RS_la%.2g.dat", la);
// time domain ()
file_TDcoeffs = Sprintf("coeff/CoeffsTD_RS_la%.2g.pro", la);
Include file_TDcoeffs;
DirRes = "res/";
po = "{Output/";
DefineConstant[
visu = {0, Choices{0, 1}, AutoCheck 0,
Name StrCat[mfem,"Visu/Real-time visualization"], Highlight "LightPink"}
Flag_FD = { 1, Choices{0,1}, Name StrCat[mfem,"2Frequency domain analysis?"], Highlight "AliceBlue" }
Flag_TD = !Flag_FD
Flag_IronCore = {1, Choices{0,1}, Name StrCat[mfem, "3Core/Iron core?"], Highlight Str[col1]}
Flag_NL = {0, Choices{0,1}, Name StrCat[mfem, "3Core/Nonlinear bh-curve?"],
Highlight Str[col1], Visible Flag_IronCore, ReadOnly Flag_FD}
Nb_max_iter = {100, Name StrCat[mfem, "3Core/Nonlinear solver/Max. num. iterations"],
Visible Flag_NL, Highlight Str[col3]}
iter_max = Nb_max_iter
relaxation_factor = {1., Name StrCat[mfem, "3Core/Nonlinear solver/Relaxation factor"],
Visible Flag_NL, Highlight Str[col3]}
stop_criterion = {1e-8, Name StrCat[mfem, "3Core/Nonlinear solver/Stopping criterion"],
Visible Flag_NL, Highlight Str[col3]} //1e-4
Flag_ImposedVoltage = {1, Choices{0,1}, Name StrCat[mfem,"Source/001Imposed Voltage?"]}
Flag_Circuit = {Flag_ImposedVoltage, Choices{0,1}, Name StrCat[mfem,"Source/002Use circuit"],
ReadOnly (Flag_ImposedVoltage==1)}
Flag_imposedRr = {0, Choices{0,1},
Name StrCat[mfem, "Source/1Imposed reduced frequency X"], Highlight Str[col1]}
Flag_imposedFreq = !Flag_imposedRr
ORDER = {1, Choices{
1="order 1",
2="order 2",
3="order 3",
4="order 4"},
Name StrCat[mfem,"1Order of homog. approx. (time)"], Visible (Flag_TD==1 && Flag_HomogenisedModel==1),
ReadOnly !(Flag_TD==1 && Flag_HomogenisedModel==1), Highlight "Red"}
ExtGmsh = ".pos"
ExtGnuplot = ".dat"
];
If(Flag_imposedRr) // Reduced frequency
DefineConstant[
Rr = {4, Min 0.1, Max 5, Step 0.1, Name StrCat[mfem,"Source/1Reduced frequency"], ReadOnly 0, Highlight "Ivory"}
delta = {Rc/Rr, Name StrCat[mfem,"Source/2Skin depth [m]"], ReadOnly 1, Highlight "LightGrey"}
Freq = {1/(delta*delta*mu0*SigmaCu*Pi), Name StrCat[mfem,"Source/3Frequency [Hz]"], ReadOnly 1, Highlight "LightGrey"}
];
Else
DefineConstant[
Freq = { 50000, Min 0.1, Max 500e3, Name StrCat[mfem,"Source/3Frequency [Hz]"], ReadOnly 0, Highlight "Ivory"}
delta = {1/Sqrt[mu0*SigmaCu*Freq*Pi], Name StrCat[mfem,"Source/2Skin depth [m]"], ReadOnly 1, Highlight "LightGrey"}
Rr = {Rc/delta, Name StrCat[mfem,"Source/1Reduced frequency"], ReadOnly 1, Highlight "LightGrey"}
];
EndIf
DefineConstant[
Omega = 2*Pi*Freq,
Period = 1./Freq,
Vdc = 50 // amplitude of PWM voltage
Val_EE = { Vdc, Name StrCat[mfem,"Source/003source amplitude"], Highlight "Ivory", Visible 1}
//------------------------------------------------------
// Time stepping variables...
NbT = {1, Name StrCat[mfem,"Source/5Nb periods"], Highlight "Ivory", Visible Flag_TD}//5
NbSteps = { 120, Name StrCat[mfem,"Source/6Nb steps per period"],
Highlight "LightGrey", Visible Flag_TD, ReadOnly 1}
tinit = {0., Name StrCat[mfem,"Source/7Initial time"], Highlight "Ivory", Visible Flag_TD}
tend = {tinit+NbT*Period , Name StrCat[mfem,"Source/8Final time"], ReadOnly 1, Highlight "LightGrey", Visible Flag_TD}
deltat = {Period/NbSteps, Name StrCat[mfem,"Source/9Step size"], ReadOnly 1, Highlight "LightGrey", Visible Flag_TD}
thetav = 1.
];
Group{
Air = Region[{AIR, AIRGAP}];
Insulation = Region[{INSULATION}];
If(Flag_IronCore)
Iron = Region[{IRON}];
Else
Iron = Region[{}];
Air += Region[{IRON}];
EndIf
OuterBoundary = Region[{OUTBND}]; // including symmetry
Winding = Region[{}] ;
DomainCC = Region[{Air, Insulation, Iron}] ;
SymFactor = Flag_HalfModel ? 2.:1. ; //half inductor with axisymmetry
nbturns = (Flag_HomogenisedModel==0) ? NbrCond/SymFactor : 1 ; // number of turns
If (Flag_HomogenisedModel==0) // Fine case
For iF In {1:nbturns}
Turn~{iF} = Region[{(iCOND+iF-1)}] ;
Winding += Region[{(iCOND+iF-1)}] ;
EndFor
DomainC = Region[{Winding}] ;
DomainS = Region[{}] ;
EndIf
If (Flag_HomogenisedModel==1) //Homogenised case
Turn~{1} = Region[{iCOND}] ;
Winding = Region[{iCOND}];
DomainC = Region[{}] ;
DomainS = Region[{Winding}] ;
EndIf
DomainCC += Region[{DomainS}] ;
If(Flag_NL)
Domain_Lin = Region[{Air, Insulation, Winding}];
Domain_Lin_NoJs = Region[{Air, Insulation}];
Domain_NonLin = Region[{Iron}];
Else
Domain_Lin = Region[{Air, Insulation, Winding, Iron}];
Domain_Lin_NoJs = Region[{Air, Insulation, Iron}];
Domain_NonLin = Region[{}];
EndIf
Domain = Region[{DomainC, DomainCC}] ;
//--------------------------------------------------------
//--------------------------------------------------------
// Groups related to source circuit
Input = # 12345 ;
// Groups related to the circuit
Input = # 12345 ;
iZH = 10000 ;
iLH = 20000 ; // Zskin in homog. coil
iLHp = 30000 ;
For k In {1:ORDER}
Zh~{k} = Region[{(iZH+k)}];
Lh~{k} = Region[{(iLH+k)}];
EndFor
For k In {2:ORDER}
Lhp~{k-1} = Region[{(iLHp+k-1)}];
EndFor
Resistance_Cir = Region[{}];
If(Flag_FD && Flag_HomogenisedModel==1) // Frequency domain
Resistance_Cir += Region[{Zh~{1}}];
EndIf
If(Flag_TD && Flag_HomogenisedModel==1) // Time domain
For k In {1:ORDER}
Resistance_Cir += Region[{Zh~{k}}];
EndFor
EndIf
Inductance_Cir = Region[{}];
If(Flag_TD && Flag_HomogenisedModel==1) // Time domain
For k In {1:ORDER}
Inductance_Cir += Region[{Lh~{k}}];
EndFor
For k In {2:ORDER}
Inductance_Cir += Region[{Lhp~{k-1}}];
EndFor
EndIf
Capacitance1_Cir = Region[ {} ] ;
Capacitance2_Cir = Region[ {} ] ;
Capacitance_Cir = Region[ {Capacitance1_Cir, Capacitance2_Cir} ] ;
SourceV_Cir = Region[ {Input} ] ;
SourceI_Cir = Region[ {} ] ;
DomainZ_Cir = Region[ {Resistance_Cir, Inductance_Cir, Capacitance_Cir} ] ;
DomainSource_Cir = Region[ {SourceV_Cir, SourceI_Cir} ] ;
DomainZt_Cir = Region[ {DomainZ_Cir, DomainSource_Cir} ] ;
}
Function {
CoefGeo = 2*Pi*SymFactor ; // axisymmetry + symmetry factor
sigma[#{Winding}] = SigmaCu ;
sigma[#{Air, Insulation, Iron}] = 0.;
rho[] = 1/sigma[];
nu[#{Air, Insulation}] = nu0;
If (!Flag_NL)
nu[#{Iron}] = nu0/1000;
Else
nu[ #{Iron} ] = nu_3kW[$1] ;
h[ #{Iron} ] = h_3kW[$1];
dhdb_NL[ #{Iron} ]= dhdb_3kW_NL[$1] ;
dhdb[ #{Iron} ] = dhdb_3kW[$1] ;
EndIf
//==================================================================
If(Flag_FD) // Frequency domain
FSinusoidal[] = Complex_MH[1,0]{Freq} ; //Cos F_Cos_wt_p[]{2*Pi*Freq, 0};
Else // Time domain
FSinusoidal[] = Complex_MH[0,-1]{Freq} ; //Sin F_Sin_wt_p[]{2*Pi*Freq, 0};
EndIf
//------------------------------------------------------
Fct_Src[] = FSinusoidal[];
//==================================================================
// Homogenization coefficients: round conductor & square packing
// Frequency domain
skin_rhor_list() = ListFromFile[ file_ZSkinRe ];
skin_rhoi_list() = ListFromFile[ file_ZSkinIm ];
prox_nur_list() = ListFromFile[ file_NuProxRe ];
prox_nui_list() = ListFromFile[ file_NuProxIm ];
skin_rhor[] = InterpolationLinear[$1]{ skin_rhor_list() };
skin_rhoi[] = InterpolationLinear[$1]{ skin_rhoi_list() };
prox_nur[] = InterpolationLinear[$1]{ prox_nur_list() } ;
prox_nui[] = InterpolationLinear[$1]{ prox_nui_list() } ;
If(Flag_HomogenisedModel==0)
nu[Winding] = nu0 ;
Else
//Proximity effect
nu[Winding] = nu0*Complex[prox_nur[Rr], prox_nui[Rr]*Fill*Rr^2/2];
EndIf
If(Flag_FD) // used if Flag_HomogenisedModel==1
// Skin effect => complex impedance
Zskin[] = 1/SymFactor*Complex[ skin_rhor[Rr]*Rdc, 2*Pi*Freq*skin_rhoi[Rr]*mu0*Len/(8*Pi*Fill)] ;
EndIf
//==================================================================
// Auxiliary functions for post-processing
nuOm[#{Air, Insulation}] = -nu[]*Complex[0.,1.];
nuOm[#{Iron}] = -nu[$1]*Complex[0.,1.];
nuOm[#{Winding}] = Complex[ Omega * Im[nu[]], -Re[nu[]] ];
kkk[] = SymFactor * skin_rhor[Rr] / Fill /SigmaCu ;
//==================================================================
// For time-domain homogenization
// Prox effect
For k In {1:ORDER}
Ca~{k}[] = nu0 * Ca_Sq~{ORDER}~{k} ;
Cb~{k}[] = SigmaCu * Fill * Rc^2/4 * Cb_Sq~{ORDER}~{k} ;
EndFor
For k In {2:ORDER}
Cc~{k-1}[] = SigmaCu * Fill * Rc^2/4 * Cc_Sq~{ORDER}~{k-1} ;
EndFor
// Skin effect
For k In {1:ORDER}
Sa~{k}[] = Rdc * Sa_Sq~{ORDER}~{k} ;
Sb~{k}[] = Rdc * SigmaCu * mu0* Rc^2/(8*Fill) * Sb_Sq~{ORDER}~{k} ;
EndFor
For k In {2:ORDER}
Sc~{k-1}[] = Rdc * SigmaCu * mu0 * Rc^2/(8*Fill) * Sc_Sq~{ORDER}~{k-1} ;
EndFor
//==================================================================
If(Flag_TD) // used if Flag_HomogenisedModel==1
For k In {1:ORDER}
Zskin~{k}[] = 1/SymFactor * Sa~{k}[];
Lskin~{k}[] = 1/SymFactor * Sb~{k}[];
EndFor
For k In {2:ORDER}
Lskin_p~{k-1}[] = 2*1/SymFactor * Sc~{k-1}[];
EndFor
Zskin[] = Zskin~{1}[];
Lskin[] = Lskin~{1}[];
EndIf
DefineFunction[
Resistance, Inductance, Capacitance
];
Ns[] = (Flag_HomogenisedModel==0) ? 1 : NbrCond/SymFactor ;
Sc[] = SurfaceArea[]{iCOND};
If (Flag_HomogenisedModel==1)
// Accounting for eddy currents (homogenization)
If(Flag_FD)
Resistance[Zh~{1}] = Zskin[] ;
EndIf
If(Flag_TD)
For k In {1:ORDER}
Resistance[Zh~{k}] = Zskin~{k}[] ;
Inductance[Lh~{k}] = Lskin~{k}[] ;
EndFor
For k In {2:ORDER}
Inductance[Lhp~{k-1}] = Lskin_p~{k-1}[] ;
EndFor
EndIf
EndIf
// List of nodes related to circuit
N1() = {1:nbturns}; // Node 1 for each turn
N2() = {2:nbturns+1}; // Node 2 for each turn
}
Constraint {
{ Name MVP_2D ;
Case {
{ Region OuterBoundary ; Type Assign ; Value 0. ; }
}
}
// Massive/stranded conductor constraints
{ Name Current_2D ;
Case {
If(Flag_Circuit==0)
{ Region Winding ; Value Val_EE; TimeFunction Fct_Src[] ; }
EndIf
}
}
{ Name Voltage_2D ;
Case{
}
}
{ Name Voltage_Cir ;
Case {
If(Flag_Circuit && Flag_ImposedVoltage)
{ Region Input ; Value Val_EE; TimeFunction Fct_Src[] ; }
EndIf
}
}
{ Name Current_Cir ;
Case {
If(Flag_Circuit && !Flag_ImposedVoltage)
{ Region Input ; Value -Val_EE; TimeFunction Fct_Src[] ; }
EndIf
}
}
{ Name ElectricalCircuit ; Type Network ;
// Common to fine and homogenised models
Case Circuit1 {
If(Flag_HomogenisedModel==0)
{ Region Input ; Branch {N1(0), N2(nbturns-1)} ; }
EndIf
If(Flag_HomogenisedModel==1 && Flag_FD)
{ Region Input ; Branch {777, N2(nbturns-1)} ; }
{ Region Zh~{1} ; Branch {777, N1(0)}; } // Complex impedance: Zskin
EndIf
If(Flag_HomogenisedModel==1 && Flag_TD)
{ Region Input ; Branch {777, N2(nbturns-1)} ; }
If(ORDER==1)
{ Region Zh~{1}; Branch {777, 800}; }
{ Region Lh~{1}; Branch {800, N1(0)}; }
EndIf
If(ORDER==2)
{ Region Zh~{1} ; Branch {777, 800}; }
{ Region Lh~{1} ; Branch {800, 801}; }
{ Region Lhp~{1}; Branch {801, N1(0)}; }
{ Region Zh~{2} ; Branch {801, 802}; }
{ Region Lh~{2} ; Branch {802, N1(0)}; }
EndIf
If(ORDER==3)
{ Region Zh~{1} ; Branch {777, 800}; }
{ Region Lh~{1} ; Branch {800, 801}; }
{ Region Lhp~{1}; Branch {801, N1(0)}; }
{ Region Zh~{2} ; Branch {801, 802}; }
{ Region Lh~{2} ; Branch {802, 803}; }
{ Region Lhp~{2}; Branch {803, N1(0)}; }
{ Region Zh~{3} ; Branch {803, 804}; }
{ Region Lh~{3} ; Branch {804, N1(0)}; }
EndIf
If(ORDER==4)
{ Region Zh~{1} ; Branch {777, 800}; }
{ Region Lh~{1} ; Branch {800, 801}; }
{ Region Lhp~{1}; Branch {801, N1(0)}; }
{ Region Zh~{2} ; Branch {801, 802}; }
{ Region Lh~{2} ; Branch {802, 803}; }
{ Region Lhp~{2}; Branch {803, N1(0)}; }
{ Region Zh~{3} ; Branch {803, 804}; }
{ Region Lh~{3} ; Branch {804, 805}; }
{ Region Lhp~{3}; Branch {805, N1(0)}; }
{ Region Zh~{4} ; Branch {805, 806}; }
{ Region Lh~{4} ; Branch {806, N1(0)}; }
EndIf
EndIf
For k In {0:nbturns-1} // list indexes start at 0
{ Region Turn~{k+1} ; Branch {N1(k), N2(k)} ; }
EndFor
}
}
}
//-----------------------------------------------------------------------------
Jacobian {
{ Name Vol ; Case { { Region All ; Jacobian VolAxiSqu ; } } }
{ Name Sur ; Case { { Region All ; Jacobian SurAxi ; } } }
}
Integration {
{ Name II ; Case {
{ Type Gauss ; Case {
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
} }
} }
}
//-----------------------------------------------------------------------------
FunctionSpace {
{ Name Hcurl_a_2D ; Type Form1P ; // Split for TD + homog: subspace isolating unknowns
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ;
Support Domain ; Entity NodesOf[ All, Not Winding ] ; }
{ Name seh ; NameOfCoef aeh ; Function BF_PerpendicularEdge ;
Support Domain ; Entity NodesOf[ Winding ] ; }
}
SubSpace {
{ Name aH ; NameOfBasisFunction {seh}; } // Subspace only used in TD homog
}
Constraint {
{ NameOfCoef ae ; EntityType NodesOf ; NameOfConstraint MVP_2D ; }
}
}
{ Name Hregion_u_2D ; Type Form1P ; // Gradient of Electric scalar potential (2D)
BasisFunction {
{ Name sr ; NameOfCoef ur ; Function BF_RegionZ ;
Support DomainC ; Entity DomainC ; }
}
GlobalQuantity {
{ Name U ; Type AliasOf ; NameOfCoef ur ; }
{ Name I ; Type AssociatedWith ; NameOfCoef ur ; }
}
Constraint {
{ NameOfCoef U ; EntityType Region ; NameOfConstraint Voltage_2D ; }
{ NameOfCoef I ; EntityType Region ; NameOfConstraint Current_2D ; }
}
}
{ Name Hregion_i_2D ; Type Vector ;
BasisFunction {
{ Name sr ; NameOfCoef ir ; Function BF_RegionZ ;
Support DomainS ; Entity DomainS ; }
}
GlobalQuantity {
{ Name Is ; Type AliasOf ; NameOfCoef ir ; }
{ Name Us ; Type AssociatedWith ; NameOfCoef ir ; }
}
Constraint {
{ NameOfCoef Us ; EntityType Region ; NameOfConstraint Voltage_2D ; }
{ NameOfCoef Is ; EntityType Region ; NameOfConstraint Current_2D ; }
}
}
// For circuit equations
{ Name Hregion_Z ; Type Scalar ;
BasisFunction {
{ Name sr ; NameOfCoef ir ; Function BF_Region ;
Support DomainZt_Cir ; Entity DomainZt_Cir ; }
}
GlobalQuantity {
{ Name Iz ; Type AliasOf ; NameOfCoef ir ; }
{ Name Uz ; Type AssociatedWith ; NameOfCoef ir ; }
}
Constraint {
{ NameOfCoef Uz ; EntityType Region ; NameOfConstraint Voltage_Cir ; }
{ NameOfCoef Iz ; EntityType Region ; NameOfConstraint Current_Cir ; }
}
}
// Time domain basis functions with homogenisation
For k In {2:ORDER}
{ Name Hcurl_a~{k} ; Type Form1P ;
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_PerpendicularEdge ;
Support Winding ; Entity NodesOf[ All ] ; }
}
}
EndFor
}
Formulation {
{ Name MagDyn_a ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
{ Name ur ; Type Local ; NameOfSpace Hregion_u_2D ; }
{ Name I ; Type Global ; NameOfSpace Hregion_u_2D[I] ; }
{ Name U ; Type Global ; NameOfSpace Hregion_u_2D[U] ; }
{ Name ir ; Type Local ; NameOfSpace Hregion_i_2D ; }
{ Name Us ; Type Global ; NameOfSpace Hregion_i_2D[Us] ; }
{ Name Is ; Type Global ; NameOfSpace Hregion_i_2D[Is] ; }
{ Name Uz ; Type Global ; NameOfSpace Hregion_Z [Uz] ; }
{ Name Iz ; Type Global ; NameOfSpace Hregion_Z [Iz] ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain_Lin ; Jacobian Vol ; Integration II ; }
If(Flag_NL)
Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
EndIf
Galerkin { DtDof [ sigma[] * Dof{a} , {a} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { [ sigma[] * Dof{ur}/CoefGeo , {a} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof [ sigma[] * Dof{a} , {ur} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { [ sigma[] * Dof{ur}/CoefGeo , {ur} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
GlobalTerm { [ Dof{I}, {U} ] ; In DomainC ; }
Galerkin { [ -Ns[]/Sc[] * Dof{ir}, {a} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof [ CoefGeo*Ns[]/Sc[] * Dof{a}, {ir} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
Galerkin { [ Ns[]/Sc[] / sigma[] * Ns[]/Sc[]* Dof{ir} , {ir} ] ; // resistance term
In DomainS ; Jacobian Vol ; Integration II ; }
GlobalTerm { [ Dof{Us}/CoefGeo , {Is} ] ; In DomainS ; }
If(Flag_Circuit)
GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { NeverDt[ Resistance[] * Dof{Iz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { [ Dof{Uz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { [ Dof{Iz} , {Iz} ] ; In Capacitance1_Cir ; }
GlobalTerm { NeverDt[ Dof{Iz} , {Iz} ] ; In Capacitance2_Cir ; }
GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ] ; In Capacitance_Cir ; }
GlobalEquation {
Type Network ; NameOfConstraint ElectricalCircuit ;
{ Node {I}; Loop {U}; Equation {I}; In DomainC ; }
{ Node {Is}; Loop {Us}; Equation {Us}; In DomainS ; }
{ Node {Iz}; Loop {Uz}; Equation {Uz}; In DomainZt_Cir ; }
}
EndIf
}
}
{ Name MagDyn_a_Homog ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
{ Name ir ; Type Local ; NameOfSpace Hregion_i_2D ; }
{ Name Us ; Type Global ; NameOfSpace Hregion_i_2D[Us] ; }
{ Name Is ; Type Global ; NameOfSpace Hregion_i_2D[Is] ; }
{ Name Uz ; Type Global ; NameOfSpace Hregion_Z [Uz] ; }
{ Name Iz ; Type Global ; NameOfSpace Hregion_Z [Iz] ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain_Lin ; Jacobian Vol ; Integration II ; }
If(Flag_NL)
Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
EndIf
Galerkin { [ -1/AreaCell * Dof{ir}, {a} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof [ 1/AreaCell * Dof{a}, {ir} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
GlobalTerm { [ Dof{Us}/CoefGeo, {Is} ] ; In DomainS ; }
// Circuit equations
If(Flag_Circuit)
GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { NeverDt[ Resistance[] * Dof{Iz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { [ Dof{Uz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { [ Dof{Iz} , {Iz} ] ; In Capacitance1_Cir ; }
GlobalTerm { NeverDt[ Dof{Iz} , {Iz} ] ; In Capacitance2_Cir ; }
GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ] ; In Capacitance_Cir ; }
GlobalEquation {
Type Network ; NameOfConstraint ElectricalCircuit ;
{ Node {Is}; Loop {Us}; Equation {Us}; In DomainS ; }
{ Node {Iz}; Loop {Uz}; Equation {Uz}; In DomainZt_Cir ; }
}
EndIf
}
}
{ Name MagDyn_a_Homog_Time ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
{ Name a_1 ; Type LocalQuantity ; NameOfSpace Hcurl_a_2D[aH] ; }
For i In {2:ORDER}
{ Name a~{i} ; Type Local ; NameOfSpace Hcurl_a~{i} ; }
EndFor
{ Name ir ; Type Local ; NameOfSpace Hregion_i_2D ; }
{ Name Us ; Type Global ; NameOfSpace Hregion_i_2D[Us] ; }
{ Name Is ; Type Global ; NameOfSpace Hregion_i_2D[Is] ; }
{ Name Uz ; Type Global ; NameOfSpace Hregion_Z [Uz] ; }
{ Name Iz ; Type Global ; NameOfSpace Hregion_Z [Iz] ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain_Lin_NoJs ; Jacobian Vol ; Integration II ; }
// Proximity effect
For i In {1:ORDER}
Galerkin { [ Ca~{i}[] * Dof{d a~{i} } , {d a~{i}} ] ;
In DomainS; Jacobian Vol; Integration II ; }
Galerkin { DtDof [ Cb~{i}[] * Dof{d a~{i}} , {d a~{i} } ] ;
In DomainS; Jacobian Vol; Integration II ; }
EndFor
For i In {2:ORDER}
Galerkin { DtDof [ Cc~{i-1}[] * Dof{d a~{i-1}} , {d a~{i}} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof [ Cc~{i-1}[] * Dof{d a~{i}} , {d a~{i-1}} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
EndFor
If(Flag_NL)
Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NonLin ; Jacobian Vol ; Integration II ; }
EndIf
Galerkin { [ -1/AreaCell * Dof{ir}, {a} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof [ 1/AreaCell * Dof{a}, {ir} ] ;
In DomainS ; Jacobian Vol ; Integration II ; }
GlobalTerm { [ Dof{Us}/CoefGeo, {Is} ] ; In DomainS ; }
// Circuit equations
If(Flag_Circuit)
GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { NeverDt[ Resistance[] * Dof{Iz} , {Iz} ] ; In Resistance_Cir ; }
GlobalTerm { [ Dof{Uz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ] ; In Inductance_Cir ; }
GlobalTerm { [ Dof{Iz} , {Iz} ] ; In Capacitance1_Cir ; }
GlobalTerm { NeverDt[ Dof{Iz} , {Iz} ] ; In Capacitance2_Cir ; }
GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ] ; In Capacitance_Cir ; }
GlobalEquation {
Type Network ; NameOfConstraint ElectricalCircuit ;
{ Node {Is}; Loop {Us}; Equation {Us}; In DomainS ; }
{ Node {Iz}; Loop {Uz}; Equation {Uz}; In DomainZt_Cir ; }
}
EndIf
}
}
}
Resolution {
{ Name Analysis ;
System {
If(Flag_HomogenisedModel==0)
If(Flag_FD) // Frequency domain
{ Name A ; NameOfFormulation MagDyn_a ; Type ComplexValue ; Frequency Freq ; }
Else // Time domain
{ Name A ; NameOfFormulation MagDyn_a ; }
EndIf
EndIf
If(Flag_HomogenisedModel==1)
If(Flag_FD) // Frequency domain
{ Name A ; NameOfFormulation MagDyn_a_Homog ; Type ComplexValue ; Frequency Freq ; }
Else // Time domain
{ Name A ; NameOfFormulation MagDyn_a_Homog_Time ; }
EndIf
EndIf
}
Operation {
CreateDir[DirRes];
InitSolution[A]; SaveSolution[A];
If(Flag_FD) // Frequency domain
Generate[A] ; Solve[A] ; SaveSolution[A] ;
If(Flag_HomogenisedModel==0)
PostOperation [Map_local];
PostOperation [Get_global];
Else
PostOperation [Map_local_Homog];
PostOperation [Get_global_Homog];
EndIf
Else //Fine reference in time domain
TimeLoopTheta[tinit, tend, deltat, thetav]{
If(!Flag_NL)
Generate[A] ; Solve[A] ;
Else
IterativeLoop[Nb_max_iter, stop_criterion, relaxation_factor] {
GenerateJac[A]; SolveJac[A]; }
EndIf
SaveSolution[A] ;
Test[ GetNumberRunTime[visu]{StrCat[mfem,"Visu/Real-time visualization"]} ]{
If(Flag_HomogenisedModel==0)
PostOperation[Map_local];
Else
PostOperation[Map_local_Homog_Time];
EndIf
}
}
EndIf
}
}
}// Resolution
PostProcessing {
{ Name MagDyn_a ; NameOfFormulation MagDyn_a ; NameOfSystem A;
PostQuantity {
{ Name a ; Value { Term { [ {a} ] ; In Domain ; Jacobian Vol ;} } }
{ Name az ; Value { Term { [ CompZ[{a}] ] ; In Domain ; Jacobian Vol ;} } }
{ Name raz ; Value { Term { [ CompZ[{a}]*X[] ] ; In Domain ; Jacobian Vol ;} } }
{ Name b ; Value { Term { [ {d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name h ; Value { Term { [ nu[{d a}]*{d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name j ; Value { Term { [ -sigma[]*(Dt[{a}]+{ur}/CoefGeo) ] ; In DomainC ; Jacobian Vol ; } } }
{ Name jz ; Value {
Term { [ CompZ[ -sigma[]*(Dt[{a}]+{ur}/CoefGeo) ] ] ; In DomainC ; Jacobian Vol ; }
Term { [ CompZ[ {ir}/AreaCond ] ] ; In DomainS ; Jacobian Vol ; }
} }
{ Name b2av ; Value { Integral { [ CoefGeo*{d a}/AreaCell ] ;
In Domain ; Jacobian Vol ; Integration II ; } } }
{ Name j2F ; Value { Integral { // Joue losses
[ CoefGeo*sigma[]*SquNorm[(Dt[{a}]+{ur}/CoefGeo)] ] ;
In DomainC ; Jacobian Vol ; Integration II ; } } }
{ Name b2F ; Value { Integral { // Magnetic Energy
[ CoefGeo*nu[{d a}]*SquNorm[{d a}] ] ;
In Domain ; Jacobian Vol ; Integration II ; } } }
{ Name SoF ; Value {
Integral {
[ CoefGeo * Complex[ sigma[] * SquNorm[(Dt[{a}]+{ur}/CoefGeo)], nu[]*SquNorm[{d a}] ] ] ;
In Domain ; Jacobian Vol ; Integration II ; }
} }//Complex power
{ Name U ; Value {
Term { [ {U} ] ; In DomainC ; }
Term { [ {Us} ] ; In DomainS ; }
Term { [ {Uz} ] ; In DomainZt_Cir ; }
} }
{ Name I ; Value {
Term { [ {I} ] ; In DomainC ; }
Term { [ {Is} ] ; In DomainS ; }
Term { [ {Iz} ] ; In DomainZt_Cir ; }
} }
{ Name MagEnergy ; Value {
Integral { [ CoefGeo*nu[{d a}]*({d a}*{d a})/2 ] ;
In Domain ; Jacobian Vol ; Integration II ; } } }
}
}
{ Name MagDyn_a_Homog ; NameOfFormulation MagDyn_a_Homog ;
PostQuantity {
{ Name a ; Value { Term { [ {a} ] ; In Domain ; Jacobian Vol ;} } }
{ Name az ; Value { Term { [ CompZ[{a}] ] ; In Domain ; Jacobian Vol ;} } }
{ Name raz ; Value { Term { [ CompZ[{a}]*X[] ] ; In Domain ; Jacobian Vol ;} } }
{ Name b ; Value { Term { [ {d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name j ; Value { Term { [ -1/AreaCell*{ir} ] ; In DomainS ; Jacobian Vol ; } } }
{ Name jz ; Value { Term { [ CompZ[ -1/AreaCell*{ir} ] ] ; In DomainS ; Jacobian Vol ; } } }
{ Name j2H ; Value { Integral { // Joule losses
[ CoefGeo*(Re[{d a}*Conj[nuOm[]*{d a}]]+kkk[]*SquNorm[-1/AreaCell*{ir}]) ] ;
In DomainS ; Jacobian Vol ; Integration II ; } } }
{ Name SoH ; Value { Integral { // Complex power = Active power +j * Reactive power => S = P+j*Q
[ CoefGeo * ({d a}*Conj[nuOm[{d a}]*{d a}] + kkk[]*SquNorm[-1/AreaCell*{ir}]) ] ;
In Domain ; Jacobian Vol ; Integration II ; } } } //Complex power
{ Name U ; Value {
Term { [ {Us} ] ; In DomainS ; }
Term { [ {Uz} ] ; In DomainZt_Cir ; }
} }
{ Name I ; Value {
Term { [ {Is} ] ; In DomainS ; }
Term { [ {Iz} ] ; In DomainZt_Cir ; }
} }
{ Name conjI ; Value {
Term { [ Conj[{Is}] ] ; In DomainS ; }
Term { [ Conj[{Iz}] ] ; In DomainZt_Cir ; }
} }
{ Name squnormI ; Value {
Term { [ SquNorm[{Is}] ] ; In DomainS ; }
Term { [ SquNorm[{Iz}] ] ; In DomainZt_Cir ; }
} }
{ Name reI ; Value {
Term { [ Re[{Is}] ] ; In DomainS ; }
Term { [ Re[{Iz}] ] ; In DomainZt_Cir ; }
} }
{ Name imI ; Value {
Term { [ Im[{Is}] ] ; In DomainS ; }
Term { [ Im[{Iz}] ] ; In DomainZt_Cir ; }
} }
}
}
{ Name MagDyn_a_Homog_Time ; NameOfFormulation MagDyn_a_Homog_Time ;
PostQuantity {
{ Name a ; Value { Term { [ {a} ] ; In Domain ; Jacobian Vol; } } }
{ Name az ; Value { Term { [ CompZ[{a}] ] ; In Domain ; Jacobian Vol; } } }
{ Name raz; Value { Term { [ CompZ[{a}]*X[]] ; In Domain ; Jacobian Vol; } } }
For i In {1:ORDER}
{ Name a~{i} ; Value { Term { [ {a~{i}} ] ; In DomainS ; Jacobian Vol; } } }
{ Name az~{i} ; Value { Term { [ CompZ[{a~{i}}] ] ; In DomainS ; Jacobian Vol; } } }
{ Name raz~{i}; Value { Term { [ CompZ[{a~{i}}]*X[] ] ; In DomainS ; Jacobian Vol; } } }
EndFor
{ Name b ; Value { Term { [ {d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name j ; Value { Term { [ -1/AreaCell*{ir} ] ; In DomainS ; Jacobian Vol ; } } }
{ Name jz ; Value { Term { [ CompZ[ -1/AreaCell*{ir} ] ] ; In DomainS ; Jacobian Vol ; } } }
{ Name jI ; Value { Term { [ Rdc*{Is}^2 ] ; In DomainS ; } } }
{ Name j2HH ; Value { // Joule Losses -- total proximity losses
For i In {1:ORDER}
Integral { [ CoefGeo * Cb~{i}[] * Dt[{d a~{i}}] * Dt[{d a~{i}}] ];
In Winding ; Jacobian Vol ; Integration II ; }
EndFor
For i In {1:ORDER-1}
Integral { [ CoefGeo * 2 * Cc~{i}[] * Dt[{d a~{i}}] * Dt[{d a~{i+1}}] ];
In Winding ; Jacobian Vol ; Integration II ; }
EndFor
}
}
{ Name j2HH_skin ; Value { // Joule Losses -- total skin losses
If(Flag_FD)
Term { [ SymFactor*Zskin[] * SquNorm[{Iz}] ]; In Zh~{1}; }
EndIf
If(Flag_TD)
For k In {1:ORDER}
Term { [ SymFactor*Zskin~{k}[] * SquNorm[{Iz}] ]; In Zh~{k}; }
EndFor
EndIf
}
}
{ Name j2F ; Value { // Joule losses
Integral { [ CoefGeo *SquNorm[1/AreaCell*{ir}]/(Fill*sigma[]) ] ;
In Winding ; Jacobian Vol ; Integration II ; }
} }
{ Name MagEnergy ; Value {
For i In {1:ORDER}
Integral { [ CoefGeo * nu0/2 * {d a~{i}}*{d a~{i}} ] ;
In Winding ; Jacobian Vol ; Integration II ; }
EndFor
}
}
{ Name U ; Value {
Term { [ {Us} ] ; In DomainS ; }
Term { [ {Uz} ] ; In DomainZt_Cir ; }
} }
{ Name I ; Value {
Term { [ {Is} ] ; In DomainS ; }
Term { [ {Iz} ] ; In DomainZt_Cir ; }
} }
}
}
}
//---------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------
PostOperation Map_local UsingPost MagDyn_a {
Print[ jz, OnElementsOf Region[{DomainC,DomainS}], File StrCat[DirRes, "jz", ExtGmsh], LastTimeStepOnly ] ;
Print[ b, OnElementsOf Domain, File StrCat[DirRes, "b", ExtGmsh], LastTimeStepOnly ] ;
Print[ raz,OnElementsOf Domain, File StrCat[DirRes, "a", ExtGmsh], LastTimeStepOnly ] ;
Echo[ Str["View[PostProcessing.NbViews-1].Light=0;
View[PostProcessing.NbViews-1].LineWidth = 2;
View[PostProcessing.NbViews-1].RangeType=3;
View[PostProcessing.NbViews-1].IntervalsType=1;
View[PostProcessing.NbViews-1].NbIso = 25;"],
File "res/option.pos"];
// RangeType = 1; // Value scale range type (1=default, 2=custom, 3=per time step)
// IntervalsType = 2; // Type of interval display (1=iso, 2=continuous, 3=discrete, 4=numeric)
If(!Flag_Circuit)
Print[ I, OnRegion Winding, Format TimeTable, >File Sprintf("res/I_f%g.dat", Freq),
SendToServer StrCat[po,"I [A]"], Color "Pink", LastTimeStepOnly];
Print[ U, OnRegion Winding, Format TimeTable, >File Sprintf("res/U_f%g.dat", Freq),
SendToServer StrCat[po,"V [V]"], Color "LightGreen", LastTimeStepOnly];
Else
Print[ I, OnRegion Input, Format TimeTable, File >Sprintf("res/I_f%g.dat", Freq),
SendToServer StrCat[po,"I [A]"], Color "Pink", LastTimeStepOnly];
Print[ U, OnRegion Input, Format TimeTable, File >Sprintf("res/U_f%g.dat", Freq),
SendToServer StrCat[po,"V [V]"], Color "LightGreen", LastTimeStepOnly];
EndIf
}
PostOperation Get_global UsingPost MagDyn_a {
Print[ j2F[ Winding ], OnGlobal, Format TimeTable, File > Sprintf("res/j2F_iron%g.dat", Flag_IronCore)] ;// Joule losses
Print[ SoF[ Domain ], OnGlobal, Format TimeTable, File > Sprintf("res/SF_iron%g.dat", Flag_IronCore)] ; // Complex power
}
PostOperation Get_allTS UsingPost MagDyn_a {
If(!Flag_Circuit)
Print[ I, OnRegion Winding, Format TimeTable, File Sprintf("res/I_f%g.dat", Freq) ];
Print[ U, OnRegion Winding, Format TimeTable, File Sprintf("res/U_f%g.dat", Freq) ];
Else
Print[ I, OnRegion Input, Format TimeTable, File Sprintf("res/I_f%g.dat", Freq),
SendToServer StrCat[po,"I [A]"]{0}, Color "Pink"];
Print[ U, OnRegion Input, Format TimeTable, File Sprintf("res/U_f%g.dat", Freq),
SendToServer StrCat[po,"V [V]"]{0}, Color "LightGreen"];
Print[ I, OnRegion Turn~{nbturns} , Format TimeTable, File Sprintf("res/Iw_f%g.dat", Freq) ];
EndIf
Print[ j2F[Winding], OnGlobal, Format TimeTable, File Sprintf("res/jl_f%g.dat", Freq) ] ; // Joule losses
}
//---------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------
PostOperation Map_local_Homog UsingPost MagDyn_a_Homog {
// Print[ jz, OnElementsOf DomainS, File StrCat[DirRes,"jH",ExtGmsh] ] ;
Print[ b, OnElementsOf Domain, File StrCat[DirRes,"bH",ExtGmsh] ] ;
Print[ raz, OnElementsOf Domain, File StrCat[DirRes,"aH",ExtGmsh] ] ;
Echo[ Str["View[PostProcessing.NbViews-1].Light=0;
View[PostProcessing.NbViews-1].LineWidth = 2;
View[PostProcessing.NbViews-1].RangeType=3;
View[PostProcessing.NbViews-1].IntervalsType=1;
View[PostProcessing.NbViews-1].NbIso = 25;"],
File "res/option.pos"];
}
PostOperation Get_global_Homog UsingPost MagDyn_a_Homog {
// Complex power: S = P + i*Q, P=active power, Q=reactive power
Print[ SoH[Domain], OnGlobal, Format TimeTable,
File Sprintf("res/SH_f%g.dat", Freq) ] ;
Print[ j2H[Winding], OnGlobal, Format Table,
File Sprintf("res/j2H_f%g.dat", Freq) ] ;
If(!Flag_Circuit)
Print[ U, OnRegion Winding, Format Table, File Sprintf("res/U_f%g.dat", Freq) ] ;
Print[ I, OnRegion Winding, Format Table, File Sprintf("res/I_f%g.dat", Freq) ] ;
Else
Print[ U, OnRegion Input, Format Table, File Sprintf("res/U_f%g.dat", Freq) ];
Print[ I, OnRegion Input, Format Table, File Sprintf("res/I_f%g.dat", Freq) ];
EndIf
}
PostOperation Map_local_Homog_Time UsingPost MagDyn_a_Homog_Time {
Print[ j, OnElementsOf DomainS, File StrCat[DirRes,"jH_time",ExtGmsh], LastTimeStepOnly ] ;
Print[ b, OnElementsOf Domain, File StrCat[DirRes,"bH_time",ExtGmsh], LastTimeStepOnly ] ;
Print[ raz, OnElementsOf Domain, File StrCat[DirRes,"aH_time",ExtGmsh], LastTimeStepOnly ] ;
Echo[ Str["View[PostProcessing.NbViews-1].Light=0;
View[PostProcessing.NbViews-1].LineWidth = 2;
View[PostProcessing.NbViews-1].RangeType=3;
View[PostProcessing.NbViews-1].IntervalsType=1;
View[PostProcessing.NbViews-1].NbIso = 25;"],
File "res/option.pos"];
}
PostOperation Get_global_Homog_Time UsingPost MagDyn_a_Homog_Time {
Print[ j2HH[Winding], OnGlobal, Format Table, LastTimeStepOnly,
File Sprintf("res/j2H_time_f%g.dat", Freq) ] ;
Print[ MagEnergy[Winding], OnGlobal, Format Table, LastTimeStepOnly,
File Sprintf("res/ME_time_f%g.dat", Freq) ] ;
If(!Flag_Circuit)
Print[ U, OnRegion Winding, Format Table, File Sprintf("res/Uh_time_f%g.dat", Freq), LastTimeStepOnly ] ;
Print[ I, OnRegion Winding, Format Table, File Sprintf("res/Ih_time_f%g.dat", Freq), LastTimeStepOnly ] ;
Else
Print[ U, OnRegion Input, Format Table, File Sprintf("res/Uh_time_f%g.dat", Freq), LastTimeStepOnly ];
Print[ I, OnRegion Input, Format Table, File Sprintf("res/Ih_time_f%g.dat", Freq), LastTimeStepOnly ];
EndIf
}
PostOperation Get_allTS_Homog_Time UsingPost MagDyn_a_Homog_Time {
If(!Flag_Circuit)
Print[ I, OnRegion Winding, Format TimeTable, File Sprintf("res/Ih_time_f%g_o%g.dat", Freq, ORDER) ];
Print[ U, OnRegion Winding, Format TimeTable, File Sprintf("res/Uh_time_f%g_o%g.dat", Freq, ORDER) ];
Else
Print[ I, OnRegion Input, Format TimeTable, File Sprintf("res/Ih_time_f%g_o%g.dat", Freq, ORDER),
SendToServer StrCat[po,"I [A]"]{0}, Color "Pink"];
Print[ U, OnRegion Input, Format TimeTable, File Sprintf("res/Uh_time_f%g_o%g.dat", Freq, ORDER),
SendToServer StrCat[po,"V [V]"]{0}, Color "LightGreen"];
EndIf
Print[ j2HH[Winding], OnGlobal, Format TimeTable, File Sprintf("res/jlh_time_f%g_o%g.dat", Freq, ORDER) ] ; // Joule losses (proximity)
}
// This is only for Onelab
DefineConstant[
R_ = {"Analysis", Name "GetDP/1ResolutionChoices", Visible 1, Closed 1}
C_ = {"-solve -v 4 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 1}
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 1}
];
// GUI
mgeo = "{0Geometrical parameters/";
mfem = "{1Analysis parameters/";
col1 = "AliceBlue";
col2 = "Blue";
col3 = "Ivory";
DefineConstant[
Flag_HomogenisedModel = {0, Choices{
0="Fine reference", 1="Homogenized"}, Name StrCat[mfem,"001Model"], Highlight Str[col2] }
Fine = !Flag_HomogenisedModel
Homo = Flag_HomogenisedModel // homogenized windings
];
DefineConstant[
MD = {1, Name StrCat[ mgeo,"01Mesh density"], Highlight Str[col1]}
MDH = {2, Name StrCat[ mgeo,"01Mesh density in winding window"], Highlight Str[col1]}
Flag_HalfModel = {1, Choices{0,1}, Name StrCat[ mgeo,"00Symmetry"], Highlight Str[col1], Closed 1}
NbrLayersX = {6, Min 1, Max 6, Name StrCat[ mgeo,"10Number of layers X"], Highlight Str[col3]}
NbrLayersY = {20, Min 2, Max 20, Step 2, Name StrCat[ mgeo,"11Number of layers Y"],
Highlight Str[col3]} // It has to be even
NbrCond = {NbrLayersX * NbrLayersY, ReadOnly 1, Highlight "LightGrey",
Name StrCat[ mgeo,"12Number of turns (total)"]}
Rc = { Sqrt[1/Pi]*1e-3, Name StrCat[ mgeo,"20Conductor radius"], Highlight Str[col3]}
Dex = { 2.2*Rc, Name StrCat[ mgeo,"21Packing side"], Highlight Str[col3]}
Dey = Dex
AreaCell = Dex*Dey
AreaCond = Pi*Rc^2
Fill = {AreaCond/AreaCell, Name StrCat[ mgeo,"30Fill factor"], ReadOnly 1, Highlight "LightGrey"}
sgap = {1., Min 1, Max 8, Name StrCat[ mgeo,"40Gap size factor"], Highlight Str[col3]}
];
// Some dimensions
X1 = 0. ;
X2 = 10e-3;
X3 = 20e-3;
X4 = 25e-3;
Y2 = 19e-3;
Y3 = 13.5e-3;
Y4 = 1.5e-3*sgap ; //Half height of the gap
// coordinates of rectangular winding window
Xw1 = 11e-3 ;
Xw2 = Xw1 + 8.e-3 ;
Yw1 = -25.5e-3/2 ;
Yw2 = 25.5e-3/2 ;
If(Fmod[NbrLayersY,2])
Printf("Warning: Number of layers along Y has to be even %g => %g", NbrLayersY, NbrLayersY+1);
NbrLayersY = NbrLayersY+1;
EndIf
// Printf("Round conductor of radius %g mm",Rc*1e3);
// Printf("Fill factor = %g mm^2 / %g mm^2 = %g ", AreaCond*1e6, AreaCell*1e6, Fill);
SigmaCu = 6e7;
mu0 = 4.e-7 * Pi ;
nu0 = 1./mu0;
//DC Resistance R_dc = L*rho/A ; L = length (m); A = area (m^2)
//Inductance of a solenoid L = mu0*mur*NbrTurns*NbrTurns*Section/Length
NbrTurns = NbrCond ;
Len = (2*Pi*(Xw1+Xw2)/2)*NbrTurns ; // Length = reserved word?
gx = (Xw2-Xw1)-NbrLayersX * Dex ;
gy = (Yw2-Yw1)-NbrLayersY * Dey ;
// correcting of rectangular winding window coordinates (fine model=> a bit smaller)
Xw1_ = Xw1+gx/2 ; // radius of gap
Xw2_ = Xw2-gx/2 ;
Yw1_ = 0 ;
Yw2_ = Yw2-gy/2 ;
DefineConstant[
// Some values computed directly from the geometry
Rdc = { Len/SigmaCu/AreaCond, Name StrCat[ mgeo,"50Total DC resistance [Ω]"], ReadOnly 1, Highlight "LightGrey"}
L = { mu0*Pi*Xw1_^2*NbrTurns^2/(2*Y4), Name StrCat[ mgeo,"51Inductance [H m⁻¹]"], ReadOnly 1, Highlight "LightGrey"}
tau = { L/Rdc, Name StrCat[ mgeo,"Time constant [s]"], ReadOnly 1, Highlight "LightGrey"}
b_gap = { mu0*NbrTurns*1./(2*Y4), Name StrCat[ mgeo,"53Induction in airgap [T]"], ReadOnly 1, Highlight "LightGrey"}
// Printf("DC resistance %g [Ohm]", Rdc);
// Printf("Inductance %g [H/m]", L);
// Printf("Induction in airgap %g [T]", b_gap);
];
//----------------------------------
// Physical numbers
//----------------------------------
OUTBND = 1111;
AIR = 1000;
AIRGAP = 1100;
IRON = 2000;
INSULATION = 3000;
iCOND = 4000;
ALLCOND = 5000;
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment