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

Homogenisation of laminations in GetDP. Files used for generating results of... 

Homogenisation of laminations in GetDP. Files used for generating results of ICS newsletter July 2020
parent 41f66e93
No related branches found
No related tags found
No related merge requests found
Pipeline #6768 passed
Stage: test
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 4075 additions and 0 deletions
HomogenisationLaminations/lam_hyst/.DS_Store 0 → 100644
+ 0
− 0
View file @ 4cf738ee
File added
HomogenisationLaminations/lam_hyst/BH.pro 0 → 100644
+ 158
− 0
View file @ 4cf738ee
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.] ;
// Parameters for Jiles-Atherton hysteresis model
//Gyselink's paper: Incorporation of a Jiles-Atherton vector hysteresis model in 2D FE magnetic field computations
Msat = 1145220; aaa = 59.5; kkk = 99.2; ccc = 0.54; alpha = 1.3e-4 ;
// FeSi -- data for Jiles-Atherton model from Bastos paper
// Oxz is the lamination plane, Oy is perpendicular to the lamination
// Ox == transverse direction (Oy taken equal to Ox, as the field is negligible)
Msat_x = 1.31e6;
a_x = 233.78;
k_x = 374.975;
c_x = 0.736;
alpha_x = 562e-6 ;
// Oz == rolling direction
Msat_z = 1.33e6;
a_z = 172.856;
k_z = 232.652;
c_z = 0.652;
alpha_z = 417e-6;
// hyst_FeSi = { Msat_x, a_x, k_x, c_x, alpha_x}; // transverse direction
// hyst_FeSi = { Msat_z, a_z, k_z, c_z, alpha_z}; // rolling direction
hyst_FeSi = { Msat, aaa, kkk, ccc, alpha};
//Using anhysteretic curve from Jiles-Atherton model
anhys_b = {
0.00000000, 0.32182278, 0.61073510, 0.79096285, 0.90997687,
0.99429048, 1.05693296, 1.10504060, 1.14290665, 1.17329715,
1.19808489, 1.21858522, 1.23574806, 1.25027452, 1.26269133,
1.27340002, 1.28271066, 1.29086537, 1.29805515, 1.30443212,
1.31011850, 1.31521331, 1.31979742, 1.32393739, 1.32768837,
1.33109640, 1.33420015, 1.33703230, 1.33962062, 1.34198886,
1.34415738, 1.34614376, 1.34796321, 1.34962896, 1.35115252,
1.35254394, 1.35381202, 1.35496447, 1.35600803, 1.35694862,
1.35779139, 1.35854083, 1.35920085, 1.35977480, 1.36026555,
1.36067549, 1.36100662, 1.36126050, 1.36143836, 1.36154102,
1.36156898 } ;
anhys_h = {
0.00000000, 31.48945679, 62.94768133, 94.34347236, 125.64569053,
156.82328930, 187.84534575, 218.68109121, 249.29994180, 279.67152877,
309.76572863, 339.55269297, 369.00287816, 398.08707454, 426.77643550,
455.04250600, 482.85725087, 510.19308256, 537.02288850, 563.32005806,
589.05850886, 614.21271269, 638.75772080, 662.66918868, 685.92340017,
708.49729101, 730.36847168, 751.51524967, 771.91665092, 791.55244067,
810.40314351, 828.45006274, 845.67529883, 862.06176725, 877.59321537,
892.25423862, 906.03029572, 918.90772314, 930.87374864, 941.91650394,
952.02503648, 961.18932028, 969.40026594, 976.64972956, 982.93052091,
988.23641048, 992.56213574, 995.90340628, 998.25690814, 999.62030702,
999.99225068 } ;
anhys_b2 = anhys_b()^2 ;
anhys_nu = anhys_h()/anhys_b() ;
anhys_nu(0) = anhys_nu(1) ;
anhys_nu_b2 = ListAlt[anhys_b2(), anhys_nu()] ;
nu_anhys[] = InterpolationLinear[SquNorm[$1]]{anhys_nu_b2()} ;
dnudb2_anhys[] = dInterpolationLinear[SquNorm[$1]]{anhys_nu_b2()} ;
h_anhys[] = nu_anhys[$1#1] * #1 ;
dhdb_anhys[] = TensorDiag[1,1,1] * nu_anhys[$1#1] + 2*dnudb2_anhys[#1] * SquDyadicProduct[#1] ;
dhdb_anhys_NL[] = 2*dnudb2_anhys[$1#1] * SquDyadicProduct[#1] ;
}
HomogenisationLaminations/lam_hyst/Macro_homogenisation_laminations.pro 0 → 100644
+ 656
− 0
View file @ 4cf738ee
This diff is collapsed.
HomogenisationLaminations/lam_hyst/homog_laminations.pro 0 → 100644
+ 447
− 0
View file @ 4cf738ee
Include "BH.pro"; // nonlinear materials
DefineConstant[
// Symmetry
AxialLength = 1,
SymmetryFactor = 1,
// Simulation parameters
visu = {1, Choices{0, 1}, AutoCheck 0,
Name StrCat[mfem,"Visu/Real-time maps"], Highlight "LawnGreen"}
Clean_Results = {0, Choices {0,1},
Name StrCat[mfem,"000Remove previous result files"], Highlight "Red"},
// Nonlinear iterations
Nb_max_iter = {50, Name StrCat[mfem,"8Nonlinear solver/Max. num. iterations"], Visible Flag_NL, Highlight "AliceBlue", Closed}
relaxation_factor = {1., Name StrCat[mfem,"8Nonlinear solver/Relaxation factor"], Visible Flag_NL, Highlight "AliceBlue"}
stop_criterion = {1e-4, Name StrCat[mfem,"8Nonlinear solver/Stopping criterion"], Visible Flag_NL, Highlight "AliceBlue"}
// time loop, defaut values
time0 = {0, Name StrCat[mfem,"7Time solver/00initial instant"], Visible Flag_AnalysisType==1, Closed}
NbT = {1, Name StrCat[mfem,"7Time solver/02nbr of periods"], Visible Flag_AnalysisType==1, Highlight "AliceBlue"}
timemax = {NbT*T, Name StrCat[mfem,"7Time solver/01final instant"], ReadOnly, Visible Flag_AnalysisType==1, Highlight "LightGrey"}
NbSteps = {2*120, Name StrCat[mfem,"7Time solver/03steps per period"], Visible Flag_AnalysisType==1, Highlight "AliceBlue"}
delta_time = {T/NbSteps, Name StrCat[mfem,"7Time solver/04step or deltatime"], ReadOnly, Visible Flag_AnalysisType==1, Highlight "LightGrey"}
// ResDir = "resFloop/",
ResDir = "res_cost/",
ExtGmsh = ".pos",
ExtGnuplot = ".dat"
R_ = {"Analysis", Name "GetDP/1ResolutionChoices", Visible 1}
C_ = {"-solve -v 3 -v2", Name "GetDP/9ComputeCommand", Visible 1}
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 1}
];
Group{
DefineGroup[DomainLam, DomainLamCC, DomainLamC];
}
Function {
mu0 = 4.e-7 * Pi ;
nu0 = 1./mu0;
mu_fe = mu0*mur_fe ;
nu_fe = nu0/mur_fe ;
nu[ #{Air,Ind} ] = 1./mu0 ;
fillfactor = (Flag_HomogType==0) ? 1.:fillfactor;
Printf("===> Checking fillfactor %g", fillfactor);
sigma [ #{Ind} ] = 5.9e7 ;
sigma [ #{Iron} ] = sigma_fe ;
rho[] = 1/sigma[] ;
nu_lin[] = nu_fe/fillfactor;
If(!NbrRegions[Domain_NL]) // Linear case
If(!Flag_ComplexReluctivity)
Printf("===> using LINEAR law");
nu[#{Iron}] = nu_lin[] ;
dhdb_NL[ #{Iron} ] = dhdb_lin_NL[$1] ;
EndIf
Else
Printf("===> using NON LINEAR law");
If(!Flag_Hysteresis)
If(Flag_NonlinearLawType==0)
Printf("===> testing NON LINEAR case with LINEAR law");
nu[#{Iron}] = nu_lin[] ;
dhdb_NL[ #{Iron} ] = dhdb_lin_NL[$1] ;
EndIf
If(Flag_NonlinearLawType==1)
Printf("===> 3kW machine NON LINEAR law");
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_NonlinearLawType==2)
Printf("===> anhysteretic NON LINEAR law corresponding to Jiles-Atherton params");
nu[ #{Iron} ] = nu_anhys[$1] ;
h[ #{Iron} ] = h_anhys[$1] ;
dhdb_NL[ #{Iron} ]= dhdb_anhys_NL[$1] ;
dhdb[ #{Iron} ] = dhdb_anhys[$1] ;
EndIf
Else
Printf("===> Hysteretic NON LINEAR law => Jiles-Atherton model");
EndIf
EndIf
relaxation_function[] = ($Iteration<Nb_max_iter/2) ? relaxation_factor : 0.2 ;
relax_max = 1;
relax_min = 1/10;
relax_numtest=10;
test_all_factors=0;
relaxation_factors_lin() = LinSpace[relax_max,relax_min,relax_numtest];
}
Jacobian {
{ Name Vol ; Case { { Region All ; Jacobian Vol ; } } }
{ Name Sur ; Case { { Region All ; Jacobian Sur ; } } }
{ Name Lin ; Case { { Region All ; Jacobian Lin ; } } }
}
Integration {
{ Name II ; Case {
{ Type Gauss ;
Case {
{ GeoElement Line ; NumberOfPoints 4 ; } // 1:20
{ GeoElement Triangle ; NumberOfPoints 6 ; } // 1, 3, 4, 6, 7, 12, 13, 16
{ GeoElement Quadrangle ; NumberOfPoints 7 ; } // 1, 3, 4, 7
}
}
}
}
{ Name I1p ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Triangle ; NumberOfPoints 1 ; }
{ GeoElement Quadrangle ; NumberOfPoints 1 ; } // 1, 3, 4, 7
}
}
}
}
}
FunctionSpace {
{ Name Hcurl_a; Type Form1;
BasisFunction {
{ Name se ; NameOfCoef ae ; Function BF_Edge ;
Support Domain ; Entity EdgesOf[ All ] ; }
}
Constraint {
{ NameOfCoef ae; EntityType EdgesOf ; NameOfConstraint MagneticVectorPotential; }
{ NameOfCoef ae; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
NameOfConstraint Gauge_av ; }
}
}
{ Name Hgrad_u; Type Form0;
BasisFunction {
{ Name su ; NameOfCoef un ; Function BF_Node ;
Support DomainC ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint ElectricScalarPotential; }
}
}
// Current in stranded coil (3D)
{ Name Hregion_i ; Type Scalar ;
BasisFunction {
{ Name sr ; NameOfCoef ir ; Function BF_Region ;
Support DomainS0 ; Entity DomainS0 ; }
}
GlobalQuantity {
{ Name Is ; Type AliasOf ; NameOfCoef ir ; }
{ Name Us ; Type AssociatedWith ; NameOfCoef ir ; }
}
Constraint {
{ NameOfCoef Us ; EntityType Region ; NameOfConstraint Voltage ; }
{ NameOfCoef Is ; EntityType Region ; NameOfConstraint Current ; }
}
}
// BFs for case with hysteresis
{ Name H_hysteresis ; Type Vector;
BasisFunction {
{ Name sex ; NameOfCoef aex ; Function BF_VolumeX ; Support Domain ; Entity VolumesOf[ All ] ; }
{ Name sey ; NameOfCoef aey ; Function BF_VolumeY ; Support Domain ; Entity VolumesOf[ All ] ; }
{ Name sez ; NameOfCoef aez ; Function BF_VolumeZ ; Support Domain ; Entity VolumesOf[ All ] ; }
}
}
}
If(Flag_HomogType)
Include "Macro_homogenisation_laminations.pro";
EndIf
Formulation {
{ Name MagStaDyn_av ; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a ; }
// Massive conductor (DomainC)
{ Name v ; Type Local ; NameOfSpace Hgrad_u ; }
// Source: Stranded conductor (DomainB)
{ Name ir ; Type Local ; NameOfSpace Hregion_i ; }
{ Name Us ; Type Global ; NameOfSpace Hregion_i[Us] ; }
{ Name Is ; Type Global ; NameOfSpace Hregion_i[Is] ; }
If(Flag_Hysteresis)
{ Name h ; Type Local ; NameOfSpace H_hysteresis ; }
EndIf
If(Flag_HomogType>0)
// Homogenized laminations
For k In {2:ORDERMAX} // b_2, b_4, b_6
{ Name b~{2*k-2} ; Type Local ; NameOfSpace Hb~{2*k-2} ; }
EndFor
If(Flag_Hysteresis)
For i In {1:nbrQuadraturePnts}
{ Name h~{i} ; Type Local ; NameOfSpace H_hysteresis~{i} ; }
EndFor
EndIf
EndIf
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In Domain_L ; Jacobian Vol ; Integration II ; }
If(NbrRegions[Domain_NL])
If( (Flag_HomogType==0) || (Flag_HomogType && ORDER==1))
If(!Flag_Hysteresis)
Galerkin { [ nu[{d a}/fillfactor] * Dof{d a} , {d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration II ; }
Galerkin { JacNL [ 1/fillfactor*dhdb_NL[{d a}] * Dof{d a} , {d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration II ; }
Else
Galerkin { [ SetVariable[
h_Jiles[ {h}[1],
{d a}[1]/fillfactor,
{d a}/fillfactor ]{List[hyst_FeSi]},
QuadraturePointIndex[] ]{$hjiles}, {d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration I1p ; }
Galerkin { JacNL[ 1/fillfactor*dhdb_Jiles[
{h},
{d a}/fillfactor,
{h}-{h}[1] ]{List[hyst_FeSi]} * Dof{d a} , {d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration I1p ; }
// h_Jiles saved in local quantity {h}
// BF is constant per element => 1 integration point is enough
Galerkin { [ Dof{h} , {h} ] ; // saving h_Jiles in local quantity h
In Domain_NL ; Jacobian Vol ; Integration I1p ; }
Galerkin { [ -GetVariable[QuadraturePointIndex[]]{$hjiles} , {h} ] ;
In Domain_NL ; Jacobian Vol ; Integration I1p ; }
EndIf
EndIf
EndIf
Galerkin { DtDof[ sigma[] * Dof{a} , {a} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { [ sigma[] * Dof{d v} , {a} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { DtDof[ sigma[] * Dof{a} , {d v} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
Galerkin { [ sigma[] * Dof{d v} , {d v} ] ;
In DomainC ; Jacobian Vol ; Integration II ; }
// Galerkin { [ -js[], {a} ] ;
// In DomainS0 ; Jacobian Vol ; Integration II ; }
Galerkin { [ -Idir[] * Dof{ir}, {a} ] ;
In DomainS0 ; Jacobian Vol ; Integration II ; }
If(Flag_HomogType>0)
If(!Flag_ComplexReluctivity && Flag_AnalysisType>0)// Either additional terms or complex nu in DomainLam
Call Macro_Terms_DomainLamC_sigma;
If(!Flag_Hysteresis)
Call Macro_Terms_DomainLam_nu;
Else
Call Macro_Terms_DomainLam_hysteresis;
EndIf
EndIf
EndIf
}
}
}
//======================================================================
Resolution {
{ Name Analysis ;
System {
If( (Flag_AnalysisType<2) || NbrRegions[Domain_NL])
{ Name Sys_Mag ; NameOfFormulation MagStaDyn_av ;}
Else
{ Name Sys_Mag ; NameOfFormulation MagStaDyn_av ; Frequency Freq; }
EndIf
}
Operation {
If(Clean_Results)
// SystemCommand[StrCat["rm -rf ", ResDir]];
SystemCommand[StrCat["rm -rf ", ResDir, "*.pos"]];
SystemCommand[StrCat["rm -rf ", ResDir, "*.dat"]];
EndIf
CreateDir[ResDir];
If(Flag_AnalysisType!=1)
If(!NbrRegions[Domain_NL])
Generate[Sys_Mag] ; Solve[Sys_Mag] ;
Else
IterativeLoop[Nb_max_iter, stop_criterion, relaxation_factor]{
GenerateJac[Sys_Mag] ; SolveJac[Sys_Mag] ; }
EndIf
SaveSolution [Sys_Mag] ;
Test[ GetNumberRunTime[visu]{StrCat[mfem,"Visu/Real-time maps"]} ]{
PostOperation[Get_LocalFields];
}
PostOperation[Get_GlobalQuantities];
Else // time domain
InitSolution[Sys_Mag];
TimeLoopTheta[time0, timemax, delta_time , 1.]{
If(!NbrRegions[Domain_NL])
Generate[Sys_Mag] ; Solve[Sys_Mag] ;
Else
IterativeLoop[Nb_max_iter, stop_criterion, relaxation_function[]]{
GenerateJac[Sys_Mag] ;
//If(!Flag_Hysteresis)
SolveJac[Sys_Mag] ;
//Else
//SolveJac_AdaptRelax[Sys_Mag, relaxation_factors_lin(), test_all_factors ] ;
//EndIf
}
EndIf
SaveSolution [Sys_Mag] ;
Test[ GetNumberRunTime[visu]{StrCat[mfem,"Visu/Real-time maps"]} ]{
PostOperation[Get_LocalFields];
}
PostOperation[Get_GlobalQuantities];
}
EndIf
}
}
}
PostProcessing {
{ Name MagStaDyn_av ; NameOfFormulation MagStaDyn_av ;
PostQuantity {
{ Name a ; Value { Local { [ {a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name b ; Value { Local { [ {d a} ] ; In Domain ; Jacobian Vol ; } } }
{ Name bz ; Value { Local { [ CompZ[{d a}] ] ; In Domain ; Jacobian Vol ; } } }
{ Name h ; Value {
Local { [ nu[]*{d a} ] ; In Domain_L ; Jacobian Vol ; }
If(!Flag_Hysteresis)
Local { [ nu[{d a}]*{d a} ] ; In Domain_NL ; Jacobian Vol ; }
Else
Local { [ {h} ] ; In Domain_NL ; Jacobian Vol ; }
EndIf
}
}
If(Flag_Hysteresis)
{ Name hdb ; Value { // Useful for computing losses in hysteretic case, by integrating in one period (steady state)
Local { [ {h}*({d a}-{d a}[1]) ] ; In Domain_NL ; Jacobian Vol ; }
}
}
EndIf
{ Name v ; Value { Term { [ {v} ] ; In DomainC ; Jacobian Vol ;} } }
{ Name e ; Value { Term { [ -(Dt[{a}]+{d v}) ] ; In DomainC ; Jacobian Vol ;} } }
{ Name j ; Value { Term { [ -sigma[]*(Dt[{a}]+{d v}) ] ; In DomainC ; Jacobian Vol ;} } }
{ Name normj ; Value { Term { [ Norm[-sigma[]*(Dt[{a}]+{d v})] ] ; In DomainC ; Jacobian Vol ;} } }
{ Name Jtot ; Value {
Integral { [ SymmetryFactor*AxialLength*(-sigma[]*(Dt[{a}]+{d v})) ] ;
In DomainC ; Jacobian Vol ; Integration II; } } }
{ Name ir ; Value { Local { [ {ir} ] ; In DomainS0 ; Jacobian Vol ; } } } // Source
{ Name Isrc ; Value { Integral { [ SymmetryFactor*AxialLength*{ir} ] ;
In DomainS0 ; Jacobian Vol ; Integration II ; } } } // Source
// { Name js ; Value { Local { [ js[] ] ; In DomainS0 ; Jacobian Vol ; } } } // Source
{ Name js ; Value { Local { [ Idir[] * {ir} ] ; In DomainS0 ; Jacobian Vol ; } } } // Source
{ Name Flux ; Value { Integral { [ SymmetryFactor*AxialLength * Idir[] * {a} ] ; // /Sc[] => to add => now in matlab files
In DomainS0 ; Jacobian Vol ; Integration II ; } } }
{ Name ComplexPower ; Value { // real part -> eddy current losses; imag part-> magnetic energy
Integral { [ SymmetryFactor*AxialLength*Complex[ SquNorm[sigma[]*(Dt[{a}]+{d v}/SymmetryFactor)]/sigma[], nu[]*SquNorm[{d a}] ] ] ;
In Iron; Jacobian Vol ; Integration II ; }
}
}
{ Name EddyCurrentLosses ;
Value {
Integral { [ SymmetryFactor*AxialLength*sigma[]*SquNorm[Dt[{a}]+{d v}/SymmetryFactor] ] ;
In Iron ; Jacobian Vol ; Integration II ; }
}
}
{ Name MagneticEnergy ;
Value {
If(!Flag_Hysteresis)
Integral {
[ SymmetryFactor*AxialLength* nu[{d a}]*SquNorm[{d a}] ] ;
In Domain ; Jacobian Vol ; Integration II ; }
Else
Integral {
[ SymmetryFactor*AxialLength* nu[{d a}]*SquNorm[{d a}] ] ;
In Domain_L ; Jacobian Vol ; Integration II ; }
Integral {
[ SymmetryFactor*AxialLength* {h}*{d a} ] ;
In Domain_NL ; Jacobian Vol ; Integration II ; }
EndIf
}
}
If(Flag_HomogType>0)
Call PostQuantities_Homog;
EndIf
}
}
}
PostOperation ViewTree UsingPost MagStaDyn_av {
// Print the tree for debugging purposes.
PrintGroup[ EdgesOfTreeIn[ { DomainGauge }, StartingOn { Surface_NoGauge_av } ],
In DomainGauge, File StrCat[ResDir,"Tree.pos"] ];
Echo[ Str["l=PostProcessing.NbViews-1;","View[l].ColorTable = { DarkRed };","View[l].LineWidth = 5;"] ,
File StrCat[ResDir,"tmp.geo"], LastTimeStepOnly] ;
}
HomogenisationLaminations/lam_hyst/lam2d.geo 0 → 100644
+ 443
− 0
View file @ 4cf738ee
Include "lam2d_data.pro" ;
ff = 1 ;
Mesh.CharacteristicLengthFactor = 0.85/ff ;
Mesh.Algorithm = 1; // 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Geometry.CopyMeshingMethod = 1; // Copy meshing method when duplicating geometrical entities?
//Mesh.SurfaceFaces = 1; // Display faces of surface mesh?
// Some common characteristic lengths
//====================================
pind = h/10 ;
lca = pind; // lc for air
Lay1 = ff*14; Pro1 = 0.8; // 14 horizontal divisions of the laminations
Lay3 = (Flag_Fine>1) ? 2/2+1 : ff*11;
Bwidth=1; // vertical divisions of one lamination, no bump if Bwidth == 1
Lay3_tot = 2/2*(nlam/2)+1 ; // 22 vertical divisions of the stack ==> homogenised model
Lay_isol_between_iron =(!Flag_HomogType) ? 4/2 : 2; //+++ only one division for isol
plam = dlam/Lay3 ;
Lay1_isol1 = ff*4 ; Bisolal1 = 1; // Isol with slight conductivity
//=====================================
// First lamination -- up to down
//=====================================
pl0[] += newp; Point(newp) = {xlam, ylam, 0, plam};
pl0[] += newp; Point(newp) = {xlam, (Flag_Fine) ? ylam-dlam : 0, 0, plam};
pl0[] += newp; Point(newp) = {0, (Flag_Fine) ? ylam-dlam : 0, 0, plam};
pl0[] += newp; Point(newp) = {0, ylam, 0, plam};
For k In {0:3}
ll0[] += newl; Line(newl) = {pl0[k], pl0[(k==3?0:k+1)]};
EndFor
Line Loop(newll) = {ll0[]};
surfIron[] +=news; Plane Surface(news) = {newll-1};
laxis_iron[] += ll0[2] ;
lborder_iron[] +=ll0[0];
// Iron
Transfinite Line{-ll0[1],ll0[3]} = Lay1 Using Progression Pro1 ;
Transfinite Line{ll0[{0,2}]} = (Flag_Fine?Lay3:Lay3_tot) Using Bump Bwidth;
// pli[] += newp; Point(newp) = {xlam+w_sc, ylam, 0, plam};
// pli[] += newp; Point(newp) = {xlam+w_sc, (Flag_Fine) ? ylam-dlam : 0, 0, plam};
// lli[] += newl; Line(newl) = {pl0[0], pli[0]};
// lli[] += newl; Line(newl) = {pli[0], pli[1]};
// lli[] += newl; Line(newl) = {pli[1], pl0[1]};
// Line Loop(newll) = {-ll0[0],lli[]};
// surfIsol_side1[] +=news; Plane Surface(news) = {newll-1};
//short circuit bar
// Transfinite Line{lli[{0,2}]} = Lay1_isol1 Using Progression Bisolal1;
// Transfinite Line{lli[{1}]} = (Flag_Fine?Lay3:Lay3_tot) Using Bump Bwidth;
//=====================================
// Inductor
//=====================================
phi = Pi/4;
dla_ind = rla_ind * Cos[phi] ;
xind1 = xlam + dla_ind ;
xind2 = xlam + dla_ind + w_ind ;
yind1 = ylam ;
yind2 = ylam + dla_ind;
pi[]+=newp ; Point(newp) = { xind1, 0, 0, pind};
pi[]+=newp ; Point(newp) = { xind2, 0, 0, pind};
pi[]+=newp ; Point(newp) = { xind2, yind1, 0, pind};
pi[]+=newp ; Point(newp) = { xind1, yind1, 0, pind};
li[]+=newl ; Line(newl) = {pi[0],pi[1]};
li[]+=newl ; Line(newl) = {pi[1],pi[2]};
li[]+=newl ; Line(newl) = {pi[2],pi[3]};
li[]+=newl ; Line(newl) = {pi[3],pi[0]};
Line Loop(newll) = {li[{0:3}]};
surfind[] += news ; Plane Surface(news) = {newll-1};
pi[]+=newp ; Point(newp) = { xlam, yind2, 0, pind};
pi[]+=newp ; Point(newp) = { xlam, yind2+w_ind, 0, pind};
pi[]+=newp ; Point(newp) = { 0, yind2+w_ind, 0, pind};
pi[]+=newp ; Point(newp) = { 0, yind2, 0, pind};
li[]+=newl ; Line(newl) = {pi[0+4],pi[1+4]};
li[]+=newl ; Line(newl) = {pi[1+4],pi[2+4]};
li[]+=newl ; Line(newl) = {pi[2+4],pi[3+4]};
li[]+=newl ; Line(newl) = {pi[3+4],pi[0+4]};
Line Loop(newll) = {li[{4:7}]};
surfind[] += news ; Plane Surface(news) = {newll-1};
ci[]+=newl ; Circle(newl) = {pi[3],pl0[0],pi[4]};
ci[]+=newl ; Circle(newl) = {pi[2],pl0[0],pi[5]};
Line Loop(newll) = {ci[0],li[4],-ci[1],li[2]};
surfind[] += news ; Plane Surface(news) = {newll-1};
// Inductor
Transfinite Line{li[{0,2,4,6}]} = 2; // 3width inductor //Ceil[w_ind/pind] ;
Transfinite Line{ci[]} = 7; // rounded part of inductor //Ceil[w_ind/pind]+1 ;
Transfinite Line{li[{1,3}]} = nlam ;
Transfinite Line{li[{5,7}]} = Lay1+3*2 ;
Transfinite Surface "*" ;
If(Flag_Recombine)
Recombine Surface "*" ;
EndIf
bndind[] += Abs(CombinedBoundary{Surface{surfind[]};});
bndind_out[]+= {-bndind[3],bndind[7],-bndind[1]};
bndind_in[]+= {bndind[5], -bndind[6], bndind[2]};
cen0 = newp ; Point(cen0) = {0, 0, 0, pind};
If(Flag_Fine)
//=====================================
// First isolation -- up to down
//=====================================
pe0[] += newp; Point(newp) = {xlam, ylam-e, 0, plam};
pe0[] += newp; Point(newp) = {0, ylam-e, 0, plam};
le0[] += newl; Line(newl) = {pl0[1], pe0[0]};
le0[] += newl; Line(newl) = {pe0[0], pe0[1]};
le0[] += newl; Line(newl) = {pe0[1], pl0[2]};
laxis[] += le0[2];
Line Loop(newll) = {le0[], -ll0[1]};
surfIsol[] +=news; Plane Surface(news) = {newll-1};
// Isolation between iron
Transfinite Line{le0[{0,2}]} = Lay_isol_between_iron ;
Transfinite Line{-le0[1]} = Lay1 Using Progression Pro1 ;
If(nlam==2)
ldown[] = le0[1];
ledge[] = le0[0];
EndIf
For ii In {1:nlam/2-1}
surfIron[] += Translate {0, -ii*e,0} { Duplicata { Surface{surfIron[0]}; } };
aux[] = Boundary{Surface{surfIron[ii]};};
skinIron[] += aux[] ;
laxis_iron[] += aux[2] ;
lborder_iron[] += aux[0];
EndFor
For ii In {1:nlam/2-1}
surfIsol[] += Translate {0, -ii*e,0} { Duplicata { Surface{surfIsol[0]}; } };
aux[] = Boundary{Surface{surfIsol[ii]};};
laxis[] += aux[2] ;
ldown[] += aux[1] ;
ledge[] += aux[0];
EndFor
nn=#ldown[];
pstomv[] = Boundary{Line{ldown[{nn-1}]};};
Translate {0, de/2,0} { Point{pstomv[]}; }
pntedge[] = Boundary{Line{ledge[#ledge[]-1]};};
lai[]+=newl; Line(newl) = {pntedge[1], pi[0]};
lai[]+=newl; Line(newl) = {pi[7], pl0[3]};
bndironisol[] = CombinedBoundary{
Surface{surfIron[],surfIsol[]};};
Transfinite Surface "*";
If(Flag_Recombine)
Recombine Surface "*" ;
EndIf
bndironisol[] -={laxis[], laxis_iron[],ldown[{#ldown[]-3:#ldown[]-1}]};
llairout = newll ; Line Loop(newll) = {bndironisol[], lai[], -bndind_in[]} ;
surfair[] += news ; Plane Surface(news) = {llairout};
lines_sym[] += {ldown[{#ldown[]-1}],lai[0]};
EndIf
If(!Flag_Fine)
surfIsol[]={};
volisol[]={};
laxis[]={};
lai[]+= newl; Line(newl) = {pl0[1], pi[0]};
lai[]+= newl; Line(newl) = {pi[7], pl0[3]};
llairout = newll ; Line Loop(newll) = {bndind_in[], -lai[{0}],-ll0[{3,0}],-lai[1]} ;
surfair[] += news ; Plane Surface(news) = {llairout};
lines_sym[] += {ll0[1],lai[0]};
EndIf
lines_symy0[] = {lines_sym[],li[0]}; // middle
lines_symx0[] = {laxis[], laxis_iron[], lai[1], li[6]}; // axis
lines_bnd[] = {li[{1,5}],ci[1]}; // without air outside the coil
// Symmetry with regard to y-axisO
//====================================================
surfIronD[] = {}; lborder_ironD[]={};
If(!Flag_Symmetry || Flag_Symmetry==1)
// For convenience, symmetry of the lines where bnd conditions are imposed
lborder_iron[] += Symmetry {1,0,0,0} { Duplicata{Line{lborder_iron[]};} };
lines_symy0[] += Symmetry {1,0,0,0} { Duplicata{Line{lines_symy0[]};} };
lines_bnd[] += Symmetry {1,0,0,0} { Duplicata{Line{lines_bnd[]};} };
nnk = (Flag_Fine)?nlam/2:1;
For k In {0:nnk-1}
surfIron[] += Symmetry {1,0,0,0} { Duplicata{Surface{surfIron[k]};} };
If(Flag_Fine)
surfIsol[] += Symmetry {1,0,0,0} { Duplicata{Surface{surfIsol[k]};} };
EndIf
EndFor
surfind[] += Symmetry {1,0,0,0} { Duplicata{Surface{surfind[]};} };
surfair[] += Symmetry {1,0,0,0} { Duplicata{Surface{surfair[]};} };
// Complete model => Symmetry of the half model
// ==============
If(Flag_Symmetry==0)
lborder_ironD[] = Symmetry {0,1,0,0} { Duplicata{Line{lborder_iron[]};} };
lines_bnd[] += Symmetry {0,1,0,0} { Duplicata{Line{lines_bnd[]};} };
kk = #surfIron[];
For k In {0:kk-1}
surfIronD[] += Symmetry {0,1,0,0} { Duplicata{Surface{surfIron[k]};} };
If(Flag_Fine)
surfIsol[] += Symmetry {0,1,0,0} { Duplicata{Surface{surfIsol[k]};} };
EndIf
EndFor
surfind[] += Symmetry {0,1,0,0} { Duplicata{Surface{surfind[]};} };
surfair[] += Symmetry {0,1,0,0} { Duplicata{Surface{surfair[]};} };
EndIf
EndIf
// Physical regions
//====================================================
Physical Surface("air",AIR) = {surfair[]} ;
Physical Surface("inductor",IND) = {surfind[]} ;
Physical Surface("insulation",ISOL) = {surfIsol[]} ;
If(Flag_Fine)
nni = nlam/2;
For ii In {0:nni-1}
// 1/4 of geometry => common to all cases (Flag_Symmetry==2)
Physical Surface(Sprintf("iron %g", ii),IRON+ii) = surfIron[{ii}];
Physical Line(Sprintf("skin iron (right side) %g",ii), SKINIRON_R+ii) = lborder_iron[{ii}];
If(Flag_Symmetry==2)
Physical Line(Sprintf("elec iron (middle) %g",ii), ELECIRON+ii) = laxis_iron[{ii}] ;
EndIf
If(Flag_Symmetry<2) // Complete ==0 or half geometry ==1
Physical Surface(Sprintf("iron %g", ii), IRON+ii) += surfIron[{ii+nni}]; // left side
Physical Line(Sprintf("skin iron (left side) %g",ii), SKINIRON_L+ii) = lborder_iron[{ii+nni}];
If(Flag_Symmetry==0)
Physical Surface(Sprintf("iron %g",nlam-ii-1),IRON+nlam-ii-1) = surfIronD[{ii, ii+nni}];
Physical Line(Sprintf("skin iron (right side) %g",nlam-ii-1), SKINIRON_R+nlam-ii-1) = lborder_ironD[{ii}];
Physical Line(Sprintf("skin iron (left side) %g",nlam-ii-1), SKINIRON_L+nlam-ii-1) = lborder_ironD[{ii+nni}];
EndIf
EndIf
allSkinIron_i[] = Abs( CombinedBoundary{ Physical Surface{IRON+ii}; } );
allSkinIron_i[] -= {laxis_iron[],lborder_iron[]};
Physical Line(Sprintf("skin iron %g", ii), SKINIRON+ii) = allSkinIron_i[] ;
If(Flag_Symmetry==0)
allSkinIron_iD[] = Abs( CombinedBoundary{ Physical Surface{IRON+nlam-ii-1}; } );
allSkinIron_iD[] -= {laxis_iron[],lborder_iron[], lborder_ironD[]};
Physical Line(Sprintf("skin iron %g", nlam-ii-1), SKINIRON+nlam-ii-1) = allSkinIron_iD[] ;
EndIf
EndFor
allIronIsol[] = Abs( CombinedBoundary{Surface{surfIron[],surfIronD[],surfIsol[]};} );
allIronIsol[] -= {lines_symy0[],lines_symx0[]};
Physical Line(Sprintf("skin iron + isol"), SKINMETAL) = allIronIsol[] ;
EndIf
If(!Flag_Fine)
Physical Surface("iron (homog)", IRON) = {surfIron[],surfIronD[]};
allSkinIron[] = Abs(CombinedBoundary{Surface{surfIron[],surfIronD[]};});
skinmetal[] = allSkinIron[];
skinmetal[] -= {lines_symy0[], lines_symx0[]};
Physical Line(Sprintf("skin metal"), SKINMETAL) = skinmetal[] ;
allSkinIron[] -= {laxis_iron[],lborder_iron[],lborder_ironD[]};
Physical Line("elec iron", ELECIRON) = laxis_iron[] ;
If(Flag_Symmetry==0) // Complete model model
Physical Line("skin iron (top, bottom)", SKINIRON) = allSkinIron[];
Physical Line("skin iron (right)", SKINIRON_R) = {lborder_iron[0], lborder_ironD[0]} ;
Physical Line("skin iron (left)" , SKINIRON_L) = {lborder_iron[1], lborder_ironD[1]} ;
EndIf
If(Flag_Symmetry==1) //1/2 of the model
Physical Line("skin iron (top)", SKINIRON) = allSkinIron[{1,3}];
Physical Line("skin iron (right)", SKINIRON_R) = lborder_iron[0] ;
Physical Line("skin iron (left)", SKINIRON_L) = lborder_iron[1] ;
EndIf
If(Flag_Symmetry==2) //1/4 of the model
Physical Line("skin iron (top)", SKINIRON) = allSkinIron[{1}];
Physical Line("skin iron (right)", SKINIRON_R) = lborder_iron[0] ;
EndIf
EndIf
//======================================================
Physical Line("outer boundary", OUTERBND) = lines_bnd[];
If(Flag_Symmetry!=0) // Flag_Symmetry==0 => Complete model: middle does not have to be define in case of complete geo
Physical Line("sym y=0", SYMMETRY_Y0) = lines_symy0[]; // All but inductor and conductors
EndIf
If(Flag_Symmetry==2)
lines_symx0[] -= laxis_iron[]; // Avoiding repetition of elements in physicals...
Physical Line("symmetry at x=0, not iron", SYMMETRY_X0) = {lines_symx0[]};
EndIf
// Not used... just testing purposes (gauging)
Physical Point(CORNER) = pi[1+1]; // 1 is at x-axis
Color Red { Surface{ surfind[] };}
Color Cyan { Surface{ surfair[] };}
Color Cyan { Surface{ surfIsol[] };}
Color SteelBlue { Surface{ surfIron[], surfIronD[] };}
// Orchid, LightBlue, Gold, LightGreen
// Inverting normals
Reverse Surface{surfIron[],surfIsol[],surfind[2]};
// Gauss integration points along the thickness of the lamination
//===============================================================
If(0)
Ngp = 5; // nbrQuadraturePoints (half lamination)
y_1 = 0.07443716949081559;
y_2 = 0.2166976970646236;
y_3 = 0.3397047841495122;
y_4 = 0.4325316833444923;
y_5 = 0.4869532642585858;
aaa = 0.5;
For ii In {0:nlam/2-1}
For kk In {1:Ngp}
pg[]+=newp ; Point(newp) = { aaa*xlam, h/2-d/2+y~{kk}*d-ii*e, 0};
pg[]+=newp ; Point(newp) = { aaa*xlam, h/2-d/2-y~{kk}*d-ii*e, 0};
EndFor
EndFor
EndIf
// Checking line for post-processing
//==================================
If(0)
dist_cen = 1e-3;
x0 = dist_cen; y0 = 0.;
x1 = w_lam/2; y1 = h/2;
npts = 30 ;
For ll In {1:nlam-1:2}
For ii In {1:npts-1}
If(ll==3)
Printf("nlam %g x1=%g",nlam, e*ll/2-dlam/2+ii*dlam/npts);
EndIf
pl[]+=newp ; Point(newp) = { x0, e*ll/2-dlam/2+ii*dlam/npts ,0};
EndFor
EndFor
EndIf
// Improving the mesh of the air around the laminations
//=============================================================
ddl = 0.3*dla_ind;
pmesh[]+=newp ; Point(newp) = {xlam/8, ylam + ddl, 0., plam};
pmesh[]+=newp ; Point(newp) = {xlam/4, ylam + ddl, 0., plam};
pmesh[]+=newp ; Point(newp) = {xlam/2-xlam/8, ylam + ddl, 0., plam};
pmesh[]+=newp ; Point(newp) = {xlam/2, ylam + ddl, 0., plam};
pmesh[]+=newp ; Point(newp) = {xlam/2+xlam/8, ylam + ddl, 0., plam};
Point{pmesh[]} In Surface {surfair[0]};
Characteristic Length pl0[3] = (Flag_Fine==0 ? 4:1) * plam;
If(!Flag_Symmetry || Flag_Symmetry==1)
pmesh_[] += Symmetry {1,0,0,0} { Duplicata{Point{pmesh[]};} };
Point{pmesh_[]} In Surface {surfair[1] };
If(!Flag_Symmetry) // Just testing
pmeshD[] += Symmetry {0,1,0,0} { Duplicata{Point{pmesh[]};} };
pmeshD_[] += Symmetry {0,1,0,0} { Duplicata{Point{pmesh_[]};} };
Point{pmeshD[]} In Surface {surfair[2]};
Point{pmeshD_[]} In Surface {surfair[3]};
EndIf
EndIf
//-------------------------------------------------------------------------------
// For nice visualization
//-------------------------------------------------------------------------------
// Hide { Point{ Point {:} }; Line{ Line {:} }; }
// Hide { Point{ Point {:} }; Line{9,13};}
// Show { Line{ linStator[], linRotor[] }; }
HomogenisationLaminations/lam_hyst/lam2d.pro 0 → 100644
+ 315
− 0
View file @ 4cf738ee
Include "lam2d_data.pro" ;
//--------------------------------------------------------------------------
// gmsh lam2d.geo -2 -setnumber Flag_HomogType 0 -o fine.msh
// gmsh lam2d.geo -2 -setnumber Flag_HomogType 1 -o block.msh (homogenized block)
// gmsh lam2d.geo -2 -setnumber Flag_HomogType 2 -o hom2.msh (homogenized lamination)
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
DefineConstant[
Flag_AnalysisType = {1,
Choices{0="Static", 1="Time domain", 2="Frequency domain"},
ServerAction Str["Reset", StrCat[mfem, "00Approx. order"]],
Name StrCat[mfem, "100Type of analysis"], Highlight "Blue",
Help Str["- Use 'Static' to compute static fields",
"- Use 'Time domain' to compute the dynamic response (transient)",
"- Use 'Frequency domain' to compute the dynamic response (steady state)"]}
Freq = {100, Min 0.1, Max 500e3, Name StrCat[mfem,"101Frequency [Hz]"], Visible Flag_AnalysisType>0, Highlight "LightGreen"}
Omega = 2*Pi*Freq
T = 1/Freq
Flag_NL = {!(Flag_AnalysisType==2), Choices{0,1},
Name StrCat[mfem,"20Nonlinear BH-curve?"], Highlight "LightPink", ReadOnly (Flag_AnalysisType==2)}
Flag_NonlinearLawType = {1,
Choices{
0="linear law (testing)",
1="3 kW machine law",
2="anhysteretic law",
3="Jiles-Atherton hysteretic law"},
Name StrCat[mfem,"21 Choose BH-curve?"], Highlight "Magenta", Visible Flag_NL,
ReadOnly (Flag_AnalysisType==2) }
Flag_Hysteresis = {(Flag_NonlinearLawType==3)?1:0, Choices{0,1},
Name StrCat[mfem,"22Jiles-Atherton hysteretic law?"],
Highlight "LightPink", ReadOnly, Visible Flag_NL}
ORDER = {(Flag_AnalysisType==0 || (Flag_Hysteresis==1))?1:(Flag_AnalysisType==2)?-2:2,
Choices {
0="exact",
-1="-1 (homog FD zero order)",
-2="-2 (homog FD 2nd order)",
-3="-3 (homog FD 4th order)",
1="1 (homog TD zero order)",
2="2 (homog TD 2nd order)",
3="3 (homog TD 4th order)",
4="4 (homog TD 6th order)" },
Name StrCat[mfem, "00Approx. order"],
Help Str["- With 'Time domain', use order>0",
"- With 'Frequency domain' all possibities hold"],
Highlight "Red", Visible Flag_HomogType, ReadOnly Flag_AnalysisType==0}
Flag_ComplexReluctivity = Flag_HomogType && !(ORDER>0) // 0 3D fine reference model OR Homog real; 1 Homog with nu complex
];
Group{
NN = (Flag_HomogType==1)? 1 : nlam/(Flag_Symmetry==0 ? 1:2) ; // Testing complete geometry
For k In {0:NN-1}
Iron~{k} = Region[{(IRON+k)}] ;
Iron += Region[{(IRON+k)}] ;
ElecIron~{k} = Region[{(ELECIRON+k)}] ;
ElecIron += Region[{ElecIron~{k}}] ;
EndFor
Ind = Region[{IND}];
AirOut = Region[{AIR}];
AirIn = Region[{ISOL}]; // between laminations
SymX0 = Region[{}];
SurfaceElec = Region[{}]; // av-formulation
SymMiddle = Region[{SYMMETRY_Y0}];
If(1)
OuterBnd = Region[{}];
Else
// Equivalent to no BC at outer boundary (with or without air outside ind)
// Using this condition or the previous gives exactly the same result
OuterBnd0 = #OUTERBND ;
PntOuterBnd = #CORNER;
OutLayerElements = ElementsOf[OuterBnd0, OnOneSideOf PntOuterBnd];
OuterBnd = Region[{OuterBnd0, -OutLayerElements}];
EndIf
If(Flag_Symmetry==2) // 1/4 of the model
SymX0 += Region[{SYMMETRY_X0}];
SymX0 += Region[{ElecIron}];
SurfaceElec += Region[{ElecIron}]; // av-formulation
EndIf
SurfaceGe0 = Region[{SymMiddle, SymX0}];
Surface_FixedMVP = Region[{SurfaceGe0, OuterBnd}];
If (Flag_HomogType<2)
Air = Region[{AirOut, AirIn}];
Else
// laminated structure is there but min mesh
// for testing purposes
Air = Region[{AirOut}];
Iron += Region[{AirIn}];
EndIf
DomainS0 = Region[{Ind}];
DomainCC = Region[{Air, Ind}];
If(Flag_HomogType==0) // Fine reference model
If(Flag_AnalysisType)
DomainC = Region[{Iron}];
Else
DomainCC += Region[{Iron}];
DomainC = Region[{}];
EndIf
Else
DomainLamC = Region[{Iron}];
DomainLamCC = Region[{}];
DomainLam = Region[{DomainLamC, DomainLamCC}];
DomainNoLamC = Region[ {} ] ;
DomainC = Region[ {DomainNoLamC} ] ;
// SkinDomainC = Region[{SkinMetal}]; //+++ for gauging => not needed when using av formulation
EndIf
If(!Flag_NL)
Domain_NL = Region[{}];
Domain_L = Region[{Air, Ind, Iron}];
Else
Domain_NL = Region[{Iron}];
Domain_L = Region[{Air, Ind}];
EndIf
Domain = Region[{Domain_L, Domain_NL}];
DomainGauge = Region[{Domain}];
Surface_NoGauge_av = Region[{Surface_FixedMVP}];
}
Function {
DefineConstant[
mur_fe = {2000, Name StrCat[mfem,"31relative mu (iron)"], Highlight "Ivory", Visible !Flag_NL}
sigma_fe = {2e6, Name StrCat[mfem,"32conductivity (iron)"], Highlight "Ivory"}
// sigma_sc = {sigma_fe/1e3, Name StrCat[mfem,"{23conductivity (sc)"], Highlight "Cyan"}
];
//Fct_Src[] = F_Cos_wt_p[]{2*Pi*Freq, (Flag_AnalysisType==1) ? Pi/2 : 0};
// With hysteresis: Damped start necessary
Trelax = 1/Freq/8;
Frelax[] = 1;//($Time < Trelax) ? 0.5 * (1. - Cos [Pi*$Time/Trelax] ) : 1. ;
Fct_Src[] = (Flag_Hysteresis==1 ? Frelax[]:1) * ((Flag_AnalysisType==1) ? F_Sin_wt_p[]{2*Pi*Freq, 0} : F_Cos_wt_p[]{2*Pi*Freq, 0}) ;
xcen = xlam;
ycen = ylam;
phi0[] = Atan2[Y[]-(Y[]>0.?1.:-1.)*ycen, X[]-(X[]>0?1:-1)*xcen];
Idir[] = (Fabs[X[]]>=xcen && Fabs[Y[]]>=ycen) ?
Vector[ -Sin[phi0[]#1], Cos[#1], 0. ] :
((Fabs[X[]]<xcen) ? Vector[(Y[]>0?-1:1), 0., 0.] : Vector[0.,(X[]>0?1:-1)*1., 0.]);
Sc[] = SurfaceArea[]{IND};
val_current = 8e5; //anhysteretic law is saturated but converges well
//js[] = val_current * Idir[] * Fct_Src[] ; // --- 21/12/2018
}
// Constraints
//============================================================
Constraint {
{ Name MagneticVectorPotential ;
Case {
{ Region SurfaceGe0 ; Value 0. ; }
{ Region OuterBnd ; Value 0. ; }
}
}
{ Name Gauge_av; Type Assign;
Case {
{ Region DomainGauge ; SubRegion Surface_NoGauge_av ; Value 0. ; }
}
}
{ Name ElectricScalarPotential ;
Case {
For k In {0:NN-1}
{ Region ElecIron~{k} ; Value 0. ; }
EndFor
}
}
{ Name Current ; Type Assign ;
Case {
{ Region Ind ; Value val_current; TimeFunction Fct_Src[] ; }
}
}
{ Name Voltage ; Type Assign ;
Case {
}
}
}
//======================================================================
dist_cen = 1e-4;
//dist_cen = w_lam/2-w_lam/2/8;
x0 = dist_cen; y0 = 0.;
x1 = w_lam/2; y1 = h/2;
NptsLam = 10*3 ;
//Npts = NptsLam * nlam/2;
newappend = Sprintf("_TD%g_nl%g_m%g", (Flag_AnalysisType==1), Flag_NonlinearLawType, Flag_HomogType);
If(Flag_HomogType==1)
newappend = Sprintf("_TD%g_nl%g_m%g_o%g", (Flag_AnalysisType==1), Flag_NonlinearLawType, Flag_HomogType, ORDER);
EndIf
ExtGmsh = StrCat[ newappend, Sprintf("_de%g_f%g.pos", de/mm, Freq) ];
ExtGplt = StrCat[ newappend, Sprintf("_de%g_f%g.dat", de/mm, Freq) ];
ExtGpltGlobal = StrCat[ newappend, Sprintf("_de%g_f%g.dat", de/mm, Freq)];
Include "homog_laminations.pro";
ResId="";
po_mag = StrCat["{9Output - Magnetics/", ResId];
//===========================================================================================
PostOperation {
{ Name Get_LocalFields ; NameOfPostProcessing MagStaDyn_av ; LastTimeStepOnly ;
Operation {
Print[ js, OnElementsOf Ind, File StrCat[ ResDir, "js", ExtGmsh] ];
// Print[ a, OnElementsOf Domain, File StrCat[ ResDir, "a", ExtGmsh] ];
If(NbrRegions[DomainC])
Print[ v, OnElementsOf DomainC, File StrCat[ ResDir, "v", ExtGmsh] ];
Print[ j, OnElementsOf DomainC, File StrCat[ ResDir, "j", ExtGmsh] ];
Echo[Str["k=PostProcessing.NbViews-1;","View[k].RangeType = 3;",
"View[k].CenterGlyphs = 0;","View[k].GlyphLocation = 1;"], File "res/maps.opt"];
EndIf
// If(Flag_Hysteresis)
// Print[ h, OnElementsOf Iron, File StrCat[ ResDir, "h", ExtGmsh] ];
// EndIf
Print[ bz, OnElementsOf Iron, File StrCat[ ResDir, "bz", ExtGmsh] ];
Echo[Str["k=PostProcessing.NbViews-1;","View[k].RangeType = 3;"], File "res/maps.opt"];
/*
If(Flag_HomogType>0)
If (ORDER>0) // Only time domain (or FD with b_k)
For k In {1:ORDER}
Print[ bz~{2*k-2}, OnElementsOf DomainLam, File StrCat[ResDir,Sprintf("bz%g",2*k-2),ExtGmsh] ] ;
Echo[ Str["k=PostProcessing.NbViews-1;","View[k].RangeType = 3;"], File "res/maps.opt"];
Print[ j~{2*k-2}, OnElementsOf DomainLam, File StrCat[ResDir,Sprintf("j%g",2*k-2),ExtGmsh] ] ;
Echo[ Str["k=PostProcessing.NbViews-1;","View[k].RangeType = 3;"], File "res/maps.opt"];
EndFor
EndIf
EndIf
*/
If(Flag_AnalysisType==2) // a choice
If(Flag_HomogType==0)
For ll In {1:nlam-1:2}
y0~{ll} = e*ll/2-dlam/2;
y1~{ll} = y0~{ll} + dlam;
Print[ b, OnLine {{x0, y0~{ll}, 0.}{x0, y1~{ll}, 0.}}{NptsLam}, Format TimeTable, LastTimeStepOnly,
File >> StrCat[ ResDir, "bl", ExtGplt] ];
Print[ j, OnLine {{x0, y0~{ll}, 0.}{x0, y1~{ll}, 0.}}{NptsLam}, Format TimeTable, LastTimeStepOnly,
File >> StrCat[ ResDir, "jl", ExtGplt] ];
EndFor
EndIf
If(Flag_HomogType>0)
If (ORDER>0) // Only time domain (or FD with b_k)
For k In {1:ORDER}
For ll In {1:nlam-1:2}
y0~{ll} = e*ll/2-dlam/2;
y1~{ll} = y0~{ll} + dlam;
Print[ bz~{2*k-2}, OnLine {{x0, y0~{ll}, 0.}{x0, y1~{ll}, 0.}}{NptsLam}, Format TimeTable, LastTimeStepOnly,
File >> StrCat[ ResDir, Sprintf("bl%g_o%g",2*k-2,ORDER), ExtGplt] ];
Print[ j~{2*k-2}, OnLine {{x0, y0~{ll}, 0.}{x0, y1~{ll}, 0.}}{NptsLam}, Format TimeTable, LastTimeStepOnly,
File >> StrCat[ ResDir, Sprintf("jl%g_o%g",2*k-2,ORDER), ExtGplt] ];
EndFor
EndFor
EndIf
EndIf
EndIf
}
}
{ Name Get_GlobalQuantities ; NameOfPostProcessing MagStaDyn_av ; Format TimeTable; LastTimeStepOnly ;
Operation {
Print[ Flux[DomainS0], OnGlobal, SendToServer StrCat[po_mag, "Flux linkage [Wb]"], Color "Ivory",
File > StrCat[ ResDir, "fl", ExtGpltGlobal] ];
Print[ EddyCurrentLosses[Iron], OnGlobal, SendToServer StrCat[po_mag, "Joule losses [W]"], Color "Ivory",
File > StrCat[ ResDir, "ecl", ExtGpltGlobal] ];
// Print[ MagneticEnergy[Iron], OnGlobal, SendToServer StrCat[po_mag, "Magnetic energy [J]"], Color "Ivory",
// File > StrCat[ ResDir, "me", ExtGpltGlobal] ];
}
}
}
HomogenisationLaminations/lam_hyst/lam2d_data.pro 0 → 100644
+ 101
− 0
View file @ 4cf738ee
mgeo = "{0Geometrical parameters/";
mgeo2 = "{0Geometrical parameters/dimensions/";
mfem = "{1Analysis parameters/";
DefineConstant[
Flag_HomogType = {0, Choices{
0="Fine reference model",
1="Homogenized block",
2="Homogenized lamination"}, Name StrCat[mgeo,"Model"] }
nlam = {2*10, Name StrCat[mgeo,"Number of laminations"]}
Flag_Symmetry = {2, Choices{
0="No symmetry",
1="1/2 of the model (at y=0)",
2="1/4 of the model (at y=0 & x=0)"},
Name StrCat[mgeo,"Symmetry at x=0"]}
Flag_Recombine = {1, Choices{0,1}, Name StrCat[mgeo, "Recombine"]}
];
// in geo file => Flag_Fine =0 homog. block, =1 ref. fine, =2 homog. lamination,
Flag_Fine = (Flag_HomogType==2) ? Flag_HomogType : !Flag_HomogType ;
SymmetryFactor = (!Flag_Symmetry) ? 1:Flag_Symmetry*2 ;
mm = 1e-3 ;
DefineConstant[
dlam = { 0.50, Name StrCat[mgeo2,"0lamination thickness"], Units 'mm', Closed 1}
de = { 0.05, Name StrCat[mgeo2,"1isolation thickness"], Units 'mm'} // 0.07
fillfactor = { dlam/(dlam+de), Name StrCat[mgeo2,"2fill factor"], Units '', ReadOnly 1}
w_lam = { 10, Name StrCat[mgeo2,"3lamination width"], Units 'mm'}
];
dlam = dlam*mm;
de = de *mm;
e = dlam+de;
h = nlam*dlam + (nlam-1)*(e-dlam); // height lamination stack
w_lam = w_lam*mm;
// inductor dimensions
r1 = 8.4*mm;
r2 = 9*mm;
w_ind = r2-r1; //1*mm ;
rla_ind = 2*e; // minimum distance between lamination and inductor (radius)
//---------------------------------------------------------------------------
xlam = w_lam/2 ; //h/2 ;
ylam = h/2 ;
x_air = 1.5*(h/2+rla_ind+w_ind);
y_air = 1.5*(h/2+rla_ind+w_ind) ;
area_ind = w_ind*w_lam/2 + w_ind*h/2 + Pi/4*((rla_ind+w_ind)^2-rla_ind^2);
Printf("Area of 1/4 of the inductor %g", area_ind);
// Some material characteristics
// ==================================
//sigmaLam = 5e6 ;
//sigmaCu = 5.9e7 ; // conductivity of copper [S/m]
//sigmaIsol = sigmaLam/1e6 ; // slightly conducting, as if we had a shortcircuit
//murIron = 2500 ; // To use in linear case
//============================
// Physical numbers
//============================
IRON = 1000 ; // all laminations or complete block (if homog); first lamination 10001 and so on.
ELECIRON = 1500 ;
SKINIRON =1100 ; // skin of all lamination or of complete block (if homog); skin of first lamination 20001 and so on
SKINIRON_R =1300; // right side
SKINIRON_L =1400; // left side
IND = 2000;
SKININD = 2100;
SKININDIN = 2101;
AIR = 6000;
ISOL = 3000; // Between laminations
SYMMETRY_Y0 = 11111; // middle
SYMMETRY_X0 = 22222; // axis
OUTERBND = 33333;
CORNER = 4444;
SKINMETAL = 5555; // Homog cases
HomogenisationLaminations/lam_hyst/matlab/.DS_Store 0 → 100644
+ 0
− 0
View file @ 4cf738ee
File added
HomogenisationLaminations/lam_hyst/matlab/PQ.pro 0 → 100644
+ 25
− 0
View file @ 4cf738ee
Function {
p_0 = 1 ;
p_2 = 1/5 ;
p_4 = 1/9 ;
p_6 = 1/13 ;
p_8 = 1/17 ;
q_0_0 = 1/12 ;
q_0_2 = -1/60 ;
q_0_4 = 0 ;
q_0_6 = 0 ;
q_0_8 = 0 ;
q_2_2 = 1/210 ;
q_2_4 = -1/1260 ;
q_2_6 = 0 ;
q_2_8 = 0 ;
q_4_4 = 1/1386 ;
q_4_6 = -1/5148 ;
q_4_8 = 0 ;
q_6_6 = 1/4290 ;
q_6_8 = -1/13260 ;
q_8_8 = 1/9690 ;
}
HomogenisationLaminations/lam_hyst/matlab/gausspoints_Ngp5.pro 0 → 100644
+ 40
− 0
View file @ 4cf738ee
Function {
nbrQuadraturePnts = 5 ;
y_1 = 0.4869532642585859 ;
y_2 = 0.4325316833444923 ;
y_3 = 0.3397047841495122 ;
y_4 = 0.2166976970646236 ;
y_5 = 0.07443716949081558 ;
w_1 = 0.06667134430868803 ;
w_2 = 0.1494513491505806 ;
w_3 = 0.2190863625159821 ;
w_4 = 0.2692667193099962 ;
w_5 = 0.2955242247147529 ;
c_2_1 = 0.9227408894325529 ;
c_2_2 = 0.6225019425809211 ;
c_2_3 = 0.1923960422444001 ;
c_2_4 = -0.2182526485213317 ;
c_2_5 = -0.4667546467891736 ;
c_4_1 = 0.7540759623195408 ;
c_4_2 = 0.01876577623813898 ;
c_4_3 = -0.4237995624971215 ;
c_4_4 = -0.1750153257920291 ;
c_4_5 = 0.2940357210203751 ;
c_6_1 = 0.5198876842062001 ;
c_6_2 = -0.37631042528707 ;
c_6_3 = -0.05814566531984765 ;
c_6_4 = 0.3212307651150516 ;
c_6_5 = -0.1765653629252029 ;
c_8_1 = 0.2555659454712284 ;
c_8_2 = -0.3351473744834802 ;
c_8_3 = 0.31798282660715 ;
c_8_4 = -0.224720190317701 ;
c_8_5 = 0.08085046072097912 ;
}
HomogenisationLaminations/lam_hyst/matlab/get_b.m 0 → 100644
+ 45
− 0
View file @ 4cf738ee
function b = get_b(b1, h1, h2, n)
% b1 - induction at previous step
% h1 - magnetic field at previous step
% h2 - magnetic field at current step
% n - number of integration points for computing b2
global mu0 Msat aa kk cc alpha
b = b1 ;
dh = (h2-h1)/n ;
for k = 0:n-1
h = (n-k)/n*h1 + k/n*h2 ;
m = b/mu0 - h ; % magnetisation
he = h + alpha * m ; % effective field
% anhysteretic magnetisation
if(abs(he)<0.01*aa)
man = Msat*he/(3*aa) ;
else
man = Msat*(cosh(he/aa)/sinh(he/aa)-aa/he) ;
end
%Irreversible magnetisation
mi = (m-cc*man) / (1-cc) ;
if (abs(he) < 0.01*aa)
dmandhe=Msat/(3.*aa) ;
else
dmandhe=Msat/aa*(1-(cosh(he/aa)/sinh(he/aa)).^2+(aa/he).^2) ;
end
if (~abs(man-mi) || (man-mi)*dh < 0.)
dmidhe = 0 ;
else
dmidhe = abs(man-mi)/(kk) ;
end
% differential susceptibility
dmdh = (cc*dmandhe + (1-cc)*dmidhe) / ...
(1 - alpha*cc*dmandhe - alpha*(1-cc)*dmidhe) ;
dbdh = mu0*(1+dmdh) ;
b = b + dbdh * dh ;
end
\ No newline at end of file
HomogenisationLaminations/lam_hyst/matlab/jiles_atherton.m 0 → 100644
+ 85
− 0
View file @ 4cf738ee
clear all; close all
% Material parameters: e.g. hyst_FeSi = { Msat, a, k, c, alpha}
global mu0 Msat aa kk cc alpha ;
Msat = 1145220 ;
aa = 59.5 ;
kk = 99.2 ;
cc = 0.54 ;
alpha = 1.3e-4 ;
mu0 = 4*pi*1e-7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
ha = 1000 ;
freq = 1 ;
T = 1/freq ;
t0 = 0 ;
nbrT = 2 ;
tmax = T*nbrT ;
nbrsteps = 200 ;
N = nbrsteps*nbrT ;
time = linspace(t0,tmax,N) ;
h = ha*sin(2*pi*freq*time) ;
b = zeros(size(h)) ;
nint = 10 ;
for k=1:N-1
b(k+1) = get_b(b(k), h(k), h(k+1), nint) ;
end
figure(1), hold on, plot(time, h, 'b','LineWidth',2)
xlabel('time(s)'), ylabel('h')
figure(2), plot(h, b, 'r','LineWidth',2), grid on
xlabel('h'), ylabel('b')
% Save anhysteretic curve for first nonlinear computation
%===============================================================================
if(1) % one for creating file
N = nbrsteps*1/4 ; % anhysteretic curve ==> increasing h, till
% max sin
figure(3), plot(h(1:N), b(1:N), 'm','LineWidth',2), grid on
xlabel('h'), ylabel('b')
filename = sprintf('BH_anhysteretic.pro');
fid = fopen(filename, 'wt');
fprintf(fid, 'Function { \n\n');
fprintf(fid, 'anhys_b = { \n' );
fprintf(fid, '%.8f, %.8f, %.8f, %.8f, %.8f,\n', b(1:N));
fprintf(fid, '%.8f } ; \n\n', b(N+1));
fprintf(fid, 'anhys_h = { \n' );
fprintf(fid, '%.8f, %.8f, %.8f, %.8f, %.8f, \n', h(1:N));
fprintf(fid, '%.8f } ; \n\n', h(N+1));
fprintf(fid, 'anhys_b2 = List[anhys_b]^2 ;\n');
fprintf(fid, 'anhys_nu = List[anhys_h]/List[anhys_b] ;\n');
fprintf(fid, 'anhys_nu(0) = anhys_nu(1) ;\n');
fprintf(fid, 'anhys_nu_b = ListAlt[anhys_b, anhys_nu] ;\n');
fprintf(fid, 'anhys_nu_b2 = ListAlt[anhys_b2, anhys_nu] ;\n');
fprintf(fid, '\n');
fprintf(fid, 'nu_anhys[] = InterpolationLinear[SquNorm[$1]]{List[anhys_nu_b2]} ;\n');
fprintf(fid, 'dnudb2_anhys[] = dInterpolationLinear[SquNorm[$1]]{List[anhys_nu_b2]} ;\n');
fprintf(fid, 'h_anhys[] = nu_anhys[$1#1] * #1 ;\n');
fprintf(fid, ['dhdb_anhys[] = TensorDiag[1,1,1] * nu_anhys[$1#1] ' ...
'+ 2*dnudb2_anhys[#1] * SquDyadicProduct[#1] ;\n']);
fprintf(fid, 'dhdb_anhys_NL[] = 2*dnudb2_anhys[$1#1] * SquDyadicProduct[#1] ;\n');
fprintf(fid, '} \n');
%================================================================
fclose(fid);
end
\ No newline at end of file
HomogenisationLaminations/lam_hyst/runHyst.sh 0 → 100755
+ 68
− 0
View file @ 4cf738ee
#!/bin/sh
homog=1
getmesh=0
if [ $homog = 0 ]; then
mymesh="fine.msh"; nn=1; else
mymesh="homog.msh"; nn=3;
fi
if [ $getmesh = 1 ]; then
gmsh lam2d.geo -nt 4 -v 4 -3 -setnumber Flag_HomogType $homog -setnumber Flag_Symmetry 2 -o $mymesh
fi
ff=5000
ordermin=1 ordermax=$nn;
law=3 # 0 linear; 1 3kw machine; 2 anhysteretic law; 3 hysteretic law
nbt=20 # number of periods
nbs=240 # number of steps per period
runthis=1
if [ $runthis = 1 ]; then
for (( i=$ordermin ; i<=$ordermax ; i++ )) ; do
echo "====================== ORDER $i ==========================" ;
getdp -v 3 -cpu lam2d.pro -msh $mymesh \
-setnumber visu 0 \
-setnumber Flag_HomogType $homog \
-setnumber ORDER $i \
-setnumber Flag_Symmetry 2 \
-setnumber Flag_AnalysisType 1 \
-setnumber Flag_NL 1 \
-setnumber Flag_NonlinearLawType $law \
-setnumber Freq $ff \
-setnumber NbT $nbt \
-setnumber NbSteps $nbs \
-sol Analysis
done
fi
# frequency loop for losses figure as function of frequency
fmin=5
fmax=5000
nbf=20
runthis=0
if [ $runthis = 1 ]; then
for (( k=1 ; k<=$nbf ; k++ )) ; do
echo "e(l($fmin)+($k-1)*l($fmax/$fmin)/($nbf-1))" | bc -l
fk=`echo "e(l($fmin)+($k-1)*l($fmax/$fmin)/($nbf-1))" | bc -l`
echo "======================> $k freq $fk <==========================" ;
for (( i=$ordermin ; i<=$ordermax ; i++ )) ; do
echo "====================== ORDER $i ==========================" ;
getdp -v 3 -cpu lam2d.pro -msh $mymesh \
-setnumber visu 0 \
-setnumber Flag_HomogType $homog \
-setnumber ORDER $i \
-setnumber Flag_Symmetry 2 \
-setnumber Flag_AnalysisType 1 \
-setnumber Flag_NL 1 \
-setnumber Flag_NonlinearLawType $law \
-setnumber Freq $fk \
-setnumber NbT $nbt \
-setnumber NbSteps $nbs \
-sol Analysis
done
done
fi
HomogenisationLaminations/srm/.DS_Store 0 → 100644
+ 0
− 0
View file @ 4cf738ee
File added
HomogenisationLaminations/srm/BH.pro 0 → 100644
+ 135
− 0
View file @ 4cf738ee
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.] ;
// Parameters for Jiles-Atherton hysteresis model
//Gyselink's paper: Incorporation of a Jiles-Atherton vector hysteresis model in 2D FE magnetic field computations
Msat = 1145220; aaa = 59.5; kkk = 99.2; ccc = 0.54; alpha = 1.3e-4 ;
hyst_FeSi = { Msat, aaa, kkk, ccc, alpha};
//Using anhysteretic curve from Jiles-Atherton model
anhys_b = {
0.00000000, 0.32182278, 0.61073510, 0.79096285, 0.90997687,
0.99429048, 1.05693296, 1.10504060, 1.14290665, 1.17329715,
1.19808489, 1.21858522, 1.23574806, 1.25027452, 1.26269133,
1.27340002, 1.28271066, 1.29086537, 1.29805515, 1.30443212,
1.31011850, 1.31521331, 1.31979742, 1.32393739, 1.32768837,
1.33109640, 1.33420015, 1.33703230, 1.33962062, 1.34198886,
1.34415738, 1.34614376, 1.34796321, 1.34962896, 1.35115252,
1.35254394, 1.35381202, 1.35496447, 1.35600803, 1.35694862,
1.35779139, 1.35854083, 1.35920085, 1.35977480, 1.36026555,
1.36067549, 1.36100662, 1.36126050, 1.36143836, 1.36154102,
1.36156898 } ;
anhys_h = {
0.00000000, 31.48945679, 62.94768133, 94.34347236, 125.64569053,
156.82328930, 187.84534575, 218.68109121, 249.29994180, 279.67152877,
309.76572863, 339.55269297, 369.00287816, 398.08707454, 426.77643550,
455.04250600, 482.85725087, 510.19308256, 537.02288850, 563.32005806,
589.05850886, 614.21271269, 638.75772080, 662.66918868, 685.92340017,
708.49729101, 730.36847168, 751.51524967, 771.91665092, 791.55244067,
810.40314351, 828.45006274, 845.67529883, 862.06176725, 877.59321537,
892.25423862, 906.03029572, 918.90772314, 930.87374864, 941.91650394,
952.02503648, 961.18932028, 969.40026594, 976.64972956, 982.93052091,
988.23641048, 992.56213574, 995.90340628, 998.25690814, 999.62030702,
999.99225068 } ;
anhys_b2 = anhys_b()^2 ;
anhys_nu = anhys_h()/anhys_b() ;
anhys_nu(0) = anhys_nu(1) ;
anhys_nu_b2 = ListAlt[anhys_b2(), anhys_nu()] ;
nu_anhys[] = InterpolationLinear[SquNorm[$1]]{anhys_nu_b2()} ;
dnudb2_anhys[] = dInterpolationLinear[SquNorm[$1]]{anhys_nu_b2()} ;
h_anhys[] = nu_anhys[$1#1] * #1 ;
dhdb_anhys[] = TensorDiag[1,1,1] * nu_anhys[$1#1] + 2*dnudb2_anhys[#1] * SquDyadicProduct[#1] ;
dhdb_anhys_NL[] = 2*dnudb2_anhys[$1#1] * SquDyadicProduct[#1] ;
}
HomogenisationLaminations/srm/Macro_homogenisation_laminations.pro 0 → 100644
+ 660
− 0
View file @ 4cf738ee
This diff is collapsed.
HomogenisationLaminations/srm/MagStaDyn_SlidingSurface.pro 0 → 100644
+ 789
− 0
View file @ 4cf738ee
This diff is collapsed.
HomogenisationLaminations/srm/generate_3d.geo 0 → 100644
+ 83
− 0
View file @ 4cf738ee
DefineConstant[
NumSlices = (Flag_AirLayer)?2:1,
SlicePhysOffset = _Offset,
SlicePhysOffset_aux = _Offset_top
];
AxialLength~{1} = dlam/(Flag_HalfLam?2:1); // width of 1 lamination (or axial length)
one~{1}[] = one[];
layer_top~{1}[] = layer_top[];
If(Flag_AirLayer)
AxialLength~{2} = zair; // adding air on top
one~{2}[] = 1;
layer_top~{2}[] = 1;
EndIf
Geometry.AutoCoherence = 0;
ps~{0}[] = Physical Surface {:};
pl~{0}[] = Physical Line {:};
pp~{0}[] = Physical Point {:};
For slice In{1:NumSlices}
For s In {0:#ps~{0}[]-1}
pv~{slice}[] += ps~{slice-1}[s] + SlicePhysOffset * slice;
ps~{slice}[] += ps~{slice-1}[s] + SlicePhysOffset_aux * slice;
ele[] = Physical Surface{ps~{slice-1}[s]};
name = StrCat( Physical Surface{ps~{0}[s]}, Sprintf(" (slice %g)", slice) );
Physical Volume(pv~{slice}[s]) = {};
Physical Surface(ps~{slice}[s]) = {};
For e In {0:#ele[]-1}
p[] = Extrude{0,0, AxialLength~{slice}}{
Surface{ ele[e] }; Recombine; Layers{one~{slice}[],layer_top~{slice}[]};
};
Physical Volume(Str(name), pv~{slice}[s]) += p[1];
Physical Surface(Str(name), ps~{slice}[s]) += p[0];
EndFor
EndFor
For l In {0:#pl~{0}[]-1}
psl~{slice}[] += pl~{slice-1}[l] + SlicePhysOffset * slice;
pl~{slice}[] += pl~{slice-1}[l] + SlicePhysOffset_aux * slice;
ele[] = Physical Line{pl~{slice-1}[l]};
name = StrCat( Physical Line{pl~{0}[l]}, Sprintf(" (slice %g)", slice) );
Physical Surface(psl~{slice}[l]) = {};
Physical Line(pl~{slice}[l]) = {};
For e In {0:#ele[]-1}
p[] = Extrude{0,0,AxialLength~{slice}}{
Line{ ele[e] }; Recombine; Layers{one~{slice}[], layer_top~{slice}[]};
};
Physical Surface(Str(name), psl~{slice}[l]) += p[1];
Physical Line(Str(name), pl~{slice}[l]) += p[0];
EndFor
EndFor
// Extruding points
For k In {0:#pp~{0}[]-1}
plp~{slice}[] += pp~{slice-1}[k] + SlicePhysOffset * slice; // lines from physical points
pp~{slice}[] += pp~{slice-1}[k] + SlicePhysOffset_aux * slice; // pnts -> aux for next layer
ele[] = Physical Point{pp~{slice-1}[k]};
name = StrCat( Physical Point{pp~{0}[k]}, Sprintf(" (slice %g)", slice) );
Physical Line(plp~{slice}[k]) = {};
Physical Point(pp~{slice}[k]) = {};
For e In {0:#ele[]-1}
p[] = Extrude{0,0,AxialLength~{slice}}{
Point{ ele[e] }; Recombine; Layers{one~{slice}[], layer_top~{slice}[]};
};
Physical Line(Str(name), plp~{slice}[k]) += p[1];
Physical Point(Str(name), pp~{slice}[k]) += p[0];
EndFor
EndFor
EndFor
Geometry.AutoCoherence = 1;
Coherence;
HomogenisationLaminations/srm/matlab/.DS_Store 0 → 100644
+ 0
− 0
View file @ 4cf738ee
File added
HomogenisationLaminations/srm/matlab/PQ.pro 0 → 100644
+ 25
− 0
View file @ 4cf738ee
Function {
p_0 = 1 ;
p_2 = 1/5 ;
p_4 = 1/9 ;
p_6 = 1/13 ;
p_8 = 1/17 ;
q_0_0 = 1/12 ;
q_0_2 = -1/60 ;
q_0_4 = 0 ;
q_0_6 = 0 ;
q_0_8 = 0 ;
q_2_2 = 1/210 ;
q_2_4 = -1/1260 ;
q_2_6 = 0 ;
q_2_8 = 0 ;
q_4_4 = 1/1386 ;
q_4_6 = -1/5148 ;
q_4_8 = 0 ;
q_6_6 = 1/4290 ;
q_6_8 = -1/13260 ;
q_8_8 = 1/9690 ;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment