From a103a5d465c0f484b7bc26d250c293d4ca2af8da Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 31 Mar 2018 08:39:54 +0200
Subject: [PATCH] nonlinear support + rework generic groups
---
Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro | 117 +++++++++++++++-----
Magnetodynamics/electromagnet.pro | 10 +-
Magnetodynamics/transfo.pro | 40 ++-----
3 files changed, 109 insertions(+), 58 deletions(-)
diff --git a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
index cf925ed..aedde5b 100644
--- a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
+++ b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
@@ -9,31 +9,40 @@
DefineConstant[
Flag_FrequencyDomain = 1, // frequency-domain or time-domain simulation
- Flag_CircuitCoupling = 0 // consider coupling with external electric circuit
+ Flag_CircuitCoupling = 0, // consider coupling with external electric circuit
+ Flag_NewtonRaphson = 1, // Newton-Raphson or Picard method for nonlinear iterations
CoefPower = 0.5, // coefficient for power calculations
Freq = 50, // frequency (for harmonic simulations)
TimeInit = 0, // intial time (for time-domain simulations)
TimeFinal = 1/50, // final time (for time-domain simulations)
DeltaTime = 1/500, // time step (for time-domain simulations)
- FE_Order = 1 // finite element order
+ FE_Order = 1, // finite element order
Val_Rint = 0, // interior radius of annulus shell transformation region (Vol_Inf_Mag)
Val_Rext = 0 // exterior radius of annulus shell transformation region (Vol_Inf_Mag)
+ NL_tol_abs = 1e-6, // absolute tolerance on residual for noninear iterations
+ NL_tol_rel = 1e-6, // relative tolerance on residual for noninear iterations
+ NL_iter_max = 20 // maximum number of noninear iterations
];
Group {
DefineGroup[
- Vol_CC_Mag, // the non-conducting part
- Vol_C_Mag, // the conducting part
- Vol_V_Mag, // a moving conducting part, with invariant mesh
+ // The full magnetic domain:
+ Vol_Mag,
+
+ // Subsets of Vol_Mag:
+ Vol_C_Mag, // massive conductors
+ Vol_S0_Mag, // stranded conductors with imposed current densities js0
+ Vol_S_Mag, // stranded conductors with imposed current, voltage or circuit coupling
+ Vol_NL_Mag, // nonlinear magnetic materials
+ Vol_V_Mag, // moving massive conducting parts (with invariant mesh)
Vol_M_Mag, // permanent magnets
- Vol_S0_Mag, // current source domain with imposed current densities js0
- Vol_S_Mag, // current source domain with imposed current, imposed voltage or
- // circuit coupling
Vol_Inf_Mag, // annulus where a infinite shell transformation is applied
+
+ // Boundaries:
Sur_FluxTube_Mag, // boundary with Neumann BC
- Sur_Perfect_Mag, // boundary of perfect conductors (i.e. non-meshed)
+ Sur_Perfect_Mag, // boundary of perfect conductors (non-meshed)
Sur_Imped_Mag // boundary of conductors approximated by a surface impedance
- // (i.e. non-meshed)
+ // (non-meshed)
];
If(Flag_CircuitCoupling)
DefineGroup[
@@ -53,12 +62,13 @@ Function {
sigma, // conductivity (in Vol_C_Mag and Vol_S_Mag)
br, // remanent magnetic flux density (in Vol_M_Mag)
js0, // source current density (in Vol_S0_Mag)
+ dhdb, // Jacobian for Newton-Raphson method (in Vol_NL_Mag)
nxh, // n x magnetic field (on Sur_FluxTube_Mag)
- Velocity, // velocity of moving part Vol_V_Mag
+ Velocity, // velocity of moving part (in Vol_V_Mag)
Ns, // number of turns (in Vol_S_Mag)
Sc, // cross-section of windings (in Vol_S_Mag)
- CoefGeos, // geometrical coefficient for 2D or 2D axi model
- Ysur // surface admittance (inverse of surface impedance Zsur) on Sur_Imped_Mag
+ CoefGeos, // geometrical coefficient for 2D or 2D axi model (in Vol_Mag)
+ Ysur // surface admittance (on Sur_Imped_Mag)
];
If(Flag_CircuitCoupling)
DefineFunction[
@@ -72,8 +82,8 @@ Function {
// End of definitions.
Group{
- // all volumes
- Vol_Mag = Region[ {Vol_CC_Mag, Vol_C_Mag, Vol_V_Mag, Vol_M_Mag, Vol_S0_Mag, Vol_S_Mag} ];
+ // all linear materials
+ Vol_L_Mag = Region[ {Vol_Mag, -Vol_NL_Mag} ];
// all volumes + surfaces on which integrals will be computed
Dom_Mag = Region[ {Vol_Mag, Sur_FluxTube_Mag, Sur_Perfect_Mag, Sur_Imped_Mag} ];
If(Flag_CircuitCoupling)
@@ -230,12 +240,24 @@ Formulation {
}
Equation {
Integral { [ nu[] * Dof{d a} , {d a} ];
- In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
+ In Vol_L_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ If(Flag_NewtonRaphson)
+ Integral { [ nu[{d a}] * {d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Integral { [ dhdb[{d a}] * Dof{d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Integral { [ - dhdb[{d a}] * {d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Else
+ Integral { [ nu[{d a}] * Dof{d a}, {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ EndIf
- Integral { [ -nu[] * br[] , {d a} ];
+ Integral { [ - nu[] * br[] , {d a} ];
In Vol_M_Mag; Jacobian Vol; Integration Gauss_v; }
- Integral { [ -js0[] , {a} ];
+ Integral { [ - js0[] , {a} ];
In Vol_S0_Mag; Jacobian Vol; Integration Gauss_v; }
Integral { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
@@ -270,8 +292,21 @@ Formulation {
}
Equation {
Integral { [ nu[] * Dof{d a} , {d a} ];
- In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
- Integral { [ -nu[] * br[] , {d a} ];
+ In Vol_L_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ If(Flag_NewtonRaphson)
+ Integral { [ nu[{d a}] * {d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Integral { [ dhdb[{d a}] * Dof{d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Integral { [ - dhdb[{d a}] * {d a} , {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ Else
+ Integral { [ nu[{d a}] * Dof{d a}, {d a} ];
+ In Vol_NL_Mag; Jacobian Vol; Integration Gauss_v; }
+ EndIf
+
+ Integral { [ - nu[] * br[] , {d a} ];
In Vol_M_Mag; Jacobian Vol; Integration Gauss_v; }
// Electric field e = -Dt[{a}]-{ur},
@@ -284,7 +319,7 @@ Formulation {
Integral { [ - sigma[] * (Velocity[] /\ Dof{d a}) , {a} ];
In Vol_V_Mag; Jacobian Vol; Integration Gauss_v; }
- Integral { [ -js0[] , {a} ];
+ Integral { [ - js0[] , {a} ];
In Vol_S0_Mag; Jacobian Vol; Integration Gauss_v; }
Integral { [ nxh[] , {a} ];
@@ -312,7 +347,6 @@ Formulation {
// js[0] should be of the form: Ns[]/Sc[] * Vector[0,0,1]
Integral { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
In Vol_S_Mag; Jacobian Vol; Integration Gauss_v; }
-
Integral { DtDof [ Ns[]/Sc[] * Dof{a} , {ir} ];
In Vol_S_Mag; Jacobian Vol; Integration Gauss_v; }
Integral { [ Ns[]/Sc[] / sigma[] * (js0[]*Vector[0,0,1]) * Dof{ir} , {ir} ];
@@ -320,7 +354,7 @@ Formulation {
GlobalTerm { [ Dof{Us} / CoefGeos[] , {Is} ]; In Vol_S_Mag; }
// Attention: CoefGeo[.] = 2*Pi for Axi
- GlobalTerm { [ -Dof{I_perfect} , {A_floating} ]; In Sur_Perfect_Mag; }
+ GlobalTerm { [ - Dof{I_perfect} , {A_floating} ]; In Sur_Perfect_Mag; }
If(Flag_CircuitCoupling)
GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ]; In Resistance_Cir; }
@@ -365,7 +399,21 @@ Resolution {
InitSolution[Sys]; // provide initial condition
TimeLoopTheta[TimeInit, TimeFinal, DeltaTime, 1.]{
// Euler implicit (1) -- Crank-Nicolson (0.5)
- Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+ Generate[Sys]; Solve[Sys];
+ If(NbrRegions[Vol_NL_Mag])
+ Generate[Sys]; GetResidual[Sys, $res0];
+ Evaluate[ $res = $res0, $iter = 0 ];
+ Print[{$iter, $res, $res / $res0},
+ Format "Residual %03g: abs %14.12e rel %14.12e"];
+ While[$res > NL_tol_abs && $res / $res0 > NL_tol_rel &&
+ $res / $res0 <= 1 && $iter < NL_iter_max]{
+ Solve[Sys]; Generate[Sys]; GetResidual[Sys, $res];
+ Evaluate[ $iter = $iter + 1 ];
+ Print[{$iter, $res, $res / $res0},
+ Format "Residual %03g: abs %14.12e rel %14.12e"];
+ }
+ EndIf
+ SaveSolution[Sys];
}
EndIf
}
@@ -375,7 +423,22 @@ Resolution {
{ Name Sys; NameOfFormulation MagSta_a_2D; }
}
Operation {
- Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+ InitSolution[Sys];
+ Generate[Sys]; Solve[Sys];
+ If(NbrRegions[Vol_NL_Mag])
+ Generate[Sys]; GetResidual[Sys, $res0];
+ Evaluate[ $res = $res0, $iter = 0 ];
+ Print[{$iter, $res, $res / $res0},
+ Format "Residual %03g: abs %14.12e rel %14.12e"];
+ While[$res > NL_tol_abs && $res / $res0 > NL_tol_rel &&
+ $res / $res0 <= 1 && $iter < NL_iter_max]{
+ Solve[Sys]; Generate[Sys]; GetResidual[Sys, $res];
+ Evaluate[ $iter = $iter + 1 ];
+ Print[{$iter, $res, $res / $res0},
+ Format "Residual %03g: abs %14.12e rel %14.12e"];
+ }
+ EndIf
+ SaveSolution[Sys];
}
}
}
@@ -399,6 +462,10 @@ PostProcessing {
Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
}
}
+ { Name norm_of_b; Value {
+ Term { [ Norm[{d a}] ]; In Vol_Mag; Jacobian Vol; }
+ }
+ }
{ Name h; Value {
Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; }
Term { [ -nu[] * br[] ]; In Vol_M_Mag; Jacobian Vol; }
diff --git a/Magnetodynamics/electromagnet.pro b/Magnetodynamics/electromagnet.pro
index 4c48dca..4e9bf06 100644
--- a/Magnetodynamics/electromagnet.pro
+++ b/Magnetodynamics/electromagnet.pro
@@ -29,11 +29,11 @@ Group {
// Abstract regions used in the "Lib_MagStaDyn_av_2D_Cir.pro" template file
// that is included below:
- Vol_CC_Mag = Region[{Air, AirInf}]; // Non-conducting regions
- Vol_C_Mag = Region[{Core}]; // Massive conducting regions
- Vol_S_Mag = Region[{Ind}]; // Stranded conductors, i.e., coils
- Vol_Inf_Mag = Region[{AirInf}]; // Annulus for infinite shell transformation
- Val_Rint = rInt; Val_Rext = rExt; // Interior and exterior radii of annulus
+ Vol_Mag = Region[{Air, Core, Ind, AirInf}]; // full magnetic domain
+ Vol_C_Mag = Region[Core]; // massive conductors
+ Vol_S_Mag = Region[Ind]; // stranded conductors (coils)
+ Vol_Inf_Mag = Region[AirInf]; // annulus for infinite shell transformation
+ Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of annulus
}
Function {
diff --git a/Magnetodynamics/transfo.pro b/Magnetodynamics/transfo.pro
index 17ae2a6..8880426 100644
--- a/Magnetodynamics/transfo.pro
+++ b/Magnetodynamics/transfo.pro
@@ -28,44 +28,28 @@ DefineConstant[
];
Group {
- /* Abstract regions that will be used in the "Lib_MagStaDyn_av_2D_Cir.pro"
- template file included below; the regions are first intialized as empty,
- before being filled with physical groups */
-
- Vol_CC_Mag = Region[{}]; // Non-conducting regions
- Vol_C_Mag = Region[{}]; // Massive conductors
- Vol_S_Mag = Region[{}]; // Stranded conductors, i.e., coils
-
- // air physical groups
+ // Physical regions:
Air = Region[{AIR_WINDOW, AIR_EXT}];
- Vol_CC_Mag += Region[Air];
-
- // exterior boundary
- Sur_Air_Ext = Region[{SUR_AIR_EXT}];
-
- // magnetic core of the transformer, assumed to be non-conducting
- Core = Region[CORE];
- Vol_CC_Mag += Region[Core];
-
- Coil_1_P = Region[COIL_1_PLUS];
- Coil_1_M = Region[COIL_1_MINUS];
+ Sur_Air_Ext = Region[SUR_AIR_EXT]; // exterior boundary
+ Core = Region[CORE]; // magnetic core of the transformer, assumed non-conducting
+ Coil_1_P = Region[COIL_1_PLUS]; // 1st coil, positive side
+ Coil_1_M = Region[COIL_1_MINUS]; // 1st coil, negative side
Coil_1 = Region[{Coil_1_P, Coil_1_M}];
-
- Coil_2_P = Region[COIL_2_PLUS];
- Coil_2_M = Region[COIL_2_MINUS];
+ Coil_2_P = Region[COIL_2_PLUS]; // 2nd coil, positive side
+ Coil_2_M = Region[COIL_2_MINUS]; // 2nd coil, negative side
Coil_2 = Region[{Coil_2_P, Coil_2_M}];
-
Coils = Region[{Coil_1, Coil_2}];
+ // Abstract regions that will be used in the "Lib_MagStaDyn_av_2D_Cir.pro"
+ // template file included below;
+ Vol_Mag = Region[{Air, Core, Coils}]; // full magnetic domain
If (type_Conds == 1)
- Vol_C_Mag += Region[{Coils}];
+ Vol_C_Mag = Region[{Coils}]; // massive conductors
ElseIf (type_Conds == 2)
- Vol_S_Mag += Region[{Coils}];
- Vol_CC_Mag += Region[{Coils}];
+ Vol_S_Mag = Region[{Coils}]; // stranded conductors (coils)
EndIf
}
-
Function {
mu0 = 4e-7*Pi;
--
GitLab