nonlinear

Magnetodynamics/electromagnet.pro

@@ -40,16 +40,25 @@ Group {
   Vol_S_Mag = Region[Ind]; // stranded conductors (coils)
   Vol_Inf_Mag = Region[AirInf]; // ring-shaped shell for infinite transformation
   Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of ring
+
+  DefineConstant[
+    TimeDomain = {0, Choices{0 = "Frequency-domain", 1 = "Time-domain"},
+      Name "Model parameters/03Analysis type"}
+    NonLinearCore = {0, Choices{0, 1}, Visible TimeDomain,
+      Name "Model parameters/Non linear B-H curve in core"}
+  ];
+  If(NonLinearCore)
+    Vol_NL_Mag = Region[Core];
+  EndIf
 }
 
 Function {
   DefineConstant[
-    murCore = {100, Name "Model parameters/Mur core"},
     Current = {0.01, Name "Model parameters/Current"},
-    TimeDomain = {0, Choices{0 = "Frequency-domain", 1 = "Time-domain"},
-      Name "Model parameters/03Analysis type"}
     frequency = {50, Visible !TimeDomain,
       Name "Model parameters/Frequency" }
+    murCore = {100, Visible !NonLinearCore,
+      Name "Model parameters/Mur core"}
   ];
 
   If(TimeDomain)
@@ -69,7 +78,36 @@ Function {
 
   mu0 = 4.e-7 * Pi;
   nu[ Region[{Air, Ind, AirInf}] ]  = 1. / mu0;
-  nu[ Core ]  = 1. / (murCore * mu0);
+
+  If(NonLinearCore)
+    Mat_core_b = {
+      0, 0.0902935, 0.194131, 0.293454, 0.397291,
+      0.496614, 0.600451, 0.690745, 0.794582, 0.893905,
+      0.997743, 1.09707, 1.19639, 1.2912, 1.39955,
+      1.49436, 1.54853, 1.59368, 1.63883, 1.693,
+      1.74266, 1.79684, 1.8465, 1.89616};
+    Mat_core_h = {
+      0, 4.336661847, 7.34801661, 9.72476491, 11.5026052,
+      13.1550314, 14.38338925, 15.72777105, 17.00435285, 17.9766304,
+      19.00397556, 20.54698273, 22.46623219, 25.98451136, 31.43230549,
+      42.54537896, 53.84673274, 66.63960707, 89.21748645, 138.2173146,
+      245.0317673, 454.3420034, 911.3848661, 1869.712483};
+    Mat_core_b2 = Mat_core_b()^2;
+    Mat_core_h2 = Mat_core_h()^2;
+    Mat_core_nu = Mat_core_h() / Mat_core_b();
+    Mat_core_nu(0) = Mat_core_nu(1);
+    Mat_core_nu_b2  = ListAlt[Mat_core_b2(), Mat_core_nu()] ;
+    nu_core[] = InterpolationAkima[ SquNorm[$1] ]{ Mat_core_nu_b2() } ;
+    dnudb2_core[] = dInterpolationAkima[SquNorm[$1]]{ Mat_core_nu_b2() } ;
+    h_core[] = nu_core[$1] * $1 ;
+    dhdb_core[] = TensorDiag[1,1,1] * nu_core[$1#1] +
+      2*dnudb2_core[#1] * SquDyadicProduct[#1]  ;
+
+    nu[ Core ]  = nu_core[$1];
+    dhdb[ Core ] = dhdb_core[$1];
+  Else
+    nu[ Core ]  = 1. / (murCore * mu0);
+  EndIf
 
   sigma[ Core ] = 1e6 / 10;
   sigma[ Ind ] = 5e7;