From 430eee1e23c5fffc05c9889b0558c1759b91e745 Mon Sep 17 00:00:00 2001
From: "Ruth Sabariego (U0089683)" <Ruth.Sabariego@esat.kuleuven.be>
Date: Tue, 4 Feb 2020 10:58:31 +0100
Subject: [PATCH] Homogenisation of windings: FD & TD
---
HomogenisationWindings/.DS_Store | Bin 0 -> 6148 bytes
HomogenisationWindings/BH.pro | 91 ++
HomogenisationWindings/cell.geo | 8 +
HomogenisationWindings/cell.pro | 262 ++++
HomogenisationWindings/cell_dat.pro | 67 +
HomogenisationWindings/cell_hexa_circ.geo | 99 ++
HomogenisationWindings/cell_square_circ.geo | 47 +
HomogenisationWindings/coeff/.DS_Store | Bin 0 -> 6148 bytes
.../coeff/CoeffsTD_RS_la0.65.pro | 91 ++
HomogenisationWindings/coeff/pB_RS_la0.65.dat | 33 +
HomogenisationWindings/coeff/pI_RS_la0.65.dat | 33 +
HomogenisationWindings/coeff/qB_RS_la0.65.dat | 33 +
HomogenisationWindings/coeff/qI_RS_la0.65.dat | 33 +
HomogenisationWindings/ind_axi.geo | 215 ++++
HomogenisationWindings/ind_axi.pro | 1096 +++++++++++++++++
HomogenisationWindings/ind_axi_dat.pro | 105 ++
16 files changed, 2213 insertions(+)
create mode 100644 HomogenisationWindings/.DS_Store
create mode 100644 HomogenisationWindings/BH.pro
create mode 100644 HomogenisationWindings/cell.geo
create mode 100644 HomogenisationWindings/cell.pro
create mode 100644 HomogenisationWindings/cell_dat.pro
create mode 100644 HomogenisationWindings/cell_hexa_circ.geo
create mode 100644 HomogenisationWindings/cell_square_circ.geo
create mode 100644 HomogenisationWindings/coeff/.DS_Store
create mode 100644 HomogenisationWindings/coeff/CoeffsTD_RS_la0.65.pro
create mode 100644 HomogenisationWindings/coeff/pB_RS_la0.65.dat
create mode 100644 HomogenisationWindings/coeff/pI_RS_la0.65.dat
create mode 100644 HomogenisationWindings/coeff/qB_RS_la0.65.dat
create mode 100644 HomogenisationWindings/coeff/qI_RS_la0.65.dat
create mode 100644 HomogenisationWindings/ind_axi.geo
create mode 100644 HomogenisationWindings/ind_axi.pro
create mode 100644 HomogenisationWindings/ind_axi_dat.pro
diff --git a/HomogenisationWindings/.DS_Store b/HomogenisationWindings/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..ff8f623880e3299a33c80a2d02ba807591f30a8e
GIT binary patch
literal 6148
zcmeHKOHRWu5PhZv6(K>FEOUk4AgaO%asWO_Q;^cQY7*3Cw%mq$u;w-_cw-yVK-eOL
zW+eOBex9-Yq}Uk%X4&QoU=Co$B50&BVlq5R9fiq48W*F$14^t==k-XSzc?g&KgBIt
z>~W2Q{b%K!wWVx%gnmqSRcWg=Sson6$Y6Qw0>u|*hKfwX&@F1t$n&}nJi$)1D%XoN
zvzx7~b>6(Z5dyQZNs_ls2E6W1o`mn2{0C-x-$%PxJv%3@Q~^~$6;K6Kfdd8P>7vcf
zJZh~9r~<0Mrvh?6L@a`Fz|Nz)I#{_BfH>u_Gqxp{5T7Vu9I*39hGx7g(Yu<EVi@nv
z@hHrT19l$0J4{GEOj!9Nc`<Hv#*f@|nAoG%s(>mmRAAz}b9w)-_{oeW`9l}&sRF9N
zKU2U=lQdcKP~mO;c074&BbFN$F`3s%gU0sYCx9KfM=l-d`AFKF7YFP-N)#D)?nM0&
MaDr&13jBfs?{0oTR{#J2
literal 0
HcmV?d00001
diff --git a/HomogenisationWindings/BH.pro b/HomogenisationWindings/BH.pro
new file mode 100644
index 0000000..f269dea
--- /dev/null
+++ b/HomogenisationWindings/BH.pro
@@ -0,0 +1,91 @@
+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.] ;
+
+}
diff --git a/HomogenisationWindings/cell.geo b/HomogenisationWindings/cell.geo
new file mode 100644
index 0000000..57755ad
--- /dev/null
+++ b/HomogenisationWindings/cell.geo
@@ -0,0 +1,8 @@
+Include "cell_dat.pro";
+
+If(CellType==0)
+ Include "cell_square_circ.geo";
+EndIf
+If(CellType==1)
+ Include "cell_hexa_circ.geo";
+EndIf
diff --git a/HomogenisationWindings/cell.pro b/HomogenisationWindings/cell.pro
new file mode 100644
index 0000000..c2f7706
--- /dev/null
+++ b/HomogenisationWindings/cell.pro
@@ -0,0 +1,262 @@
+//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
diff --git a/HomogenisationWindings/cell_dat.pro b/HomogenisationWindings/cell_dat.pro
new file mode 100644
index 0000000..26fa527
--- /dev/null
+++ b/HomogenisationWindings/cell_dat.pro
@@ -0,0 +1,67 @@
+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;
diff --git a/HomogenisationWindings/cell_hexa_circ.geo b/HomogenisationWindings/cell_hexa_circ.geo
new file mode 100644
index 0000000..794510e
--- /dev/null
+++ b/HomogenisationWindings/cell_hexa_circ.geo
@@ -0,0 +1,99 @@
+//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[];
diff --git a/HomogenisationWindings/cell_square_circ.geo b/HomogenisationWindings/cell_square_circ.geo
new file mode 100644
index 0000000..c5453c3
--- /dev/null
+++ b/HomogenisationWindings/cell_square_circ.geo
@@ -0,0 +1,47 @@
+//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[]};
diff --git a/HomogenisationWindings/coeff/.DS_Store b/HomogenisationWindings/coeff/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
GIT binary patch
literal 6148
zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3
zem<@ulZcFPQ@L2!n>{z**<q8>++&mCkOWA81W14cNZ<zv;LbK1Poaz?KmsK2CSc!(
z0ynLxE!0092;Krf2c+FF_Fe*7ECH>lEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ
zLs35+`xjp>T0<F0fCPF1$Cyrb|F7^5{eNG?83~ZUUlGt@xh*qZDeu<Z%US-OSsOPv
j)R!Z4KLME7ReXlK;d!wEw5GODWMKRea10D2@KpjYNUI8I
literal 0
HcmV?d00001
diff --git a/HomogenisationWindings/coeff/CoeffsTD_RS_la0.65.pro b/HomogenisationWindings/coeff/CoeffsTD_RS_la0.65.pro
new file mode 100644
index 0000000..fff4119
--- /dev/null
+++ b/HomogenisationWindings/coeff/CoeffsTD_RS_la0.65.pro
@@ -0,0 +1,91 @@
+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;
+*/
+}
diff --git a/HomogenisationWindings/coeff/pB_RS_la0.65.dat b/HomogenisationWindings/coeff/pB_RS_la0.65.dat
new file mode 100644
index 0000000..25e7adc
--- /dev/null
+++ b/HomogenisationWindings/coeff/pB_RS_la0.65.dat
@@ -0,0 +1,33 @@
+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
diff --git a/HomogenisationWindings/coeff/pI_RS_la0.65.dat b/HomogenisationWindings/coeff/pI_RS_la0.65.dat
new file mode 100644
index 0000000..8e094a4
--- /dev/null
+++ b/HomogenisationWindings/coeff/pI_RS_la0.65.dat
@@ -0,0 +1,33 @@
+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
diff --git a/HomogenisationWindings/coeff/qB_RS_la0.65.dat b/HomogenisationWindings/coeff/qB_RS_la0.65.dat
new file mode 100644
index 0000000..078e2a3
--- /dev/null
+++ b/HomogenisationWindings/coeff/qB_RS_la0.65.dat
@@ -0,0 +1,33 @@
+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
diff --git a/HomogenisationWindings/coeff/qI_RS_la0.65.dat b/HomogenisationWindings/coeff/qI_RS_la0.65.dat
new file mode 100644
index 0000000..0574466
--- /dev/null
+++ b/HomogenisationWindings/coeff/qI_RS_la0.65.dat
@@ -0,0 +1,33 @@
+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
diff --git a/HomogenisationWindings/ind_axi.geo b/HomogenisationWindings/ind_axi.geo
new file mode 100644
index 0000000..5bdea6b
--- /dev/null
+++ b/HomogenisationWindings/ind_axi.geo
@@ -0,0 +1,215 @@
+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[]};
+}
diff --git a/HomogenisationWindings/ind_axi.pro b/HomogenisationWindings/ind_axi.pro
new file mode 100644
index 0000000..1a7d5c4
--- /dev/null
+++ b/HomogenisationWindings/ind_axi.pro
@@ -0,0 +1,1096 @@
+// 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}
+];
diff --git a/HomogenisationWindings/ind_axi_dat.pro b/HomogenisationWindings/ind_axi_dat.pro
new file mode 100644
index 0000000..e2889f2
--- /dev/null
+++ b/HomogenisationWindings/ind_axi_dat.pro
@@ -0,0 +1,105 @@
+// 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;
--
GitLab