From 282c890abb318f78d17f38c1151f5a5fc76d08fd Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 14 Mar 2018 12:39:02 +0100
Subject: [PATCH] first pass at tutorial 7
---
Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro | 468 ++++++++++++++++++++
Magnetodynamics/electromagnet.geo | 54 +++
Magnetodynamics/electromagnet.pro | 77 ++++
Magnetodynamics/electromagnet_common.pro | 4 +
Magnetodynamics/transfo.geo | 173 ++++++++
Magnetodynamics/transfo.pro | 260 +++++++++++
Magnetodynamics/transfo_common.pro | 49 ++
Magnetostatics/electromagnet.pro | 8 +-
8 files changed, 1089 insertions(+), 4 deletions(-)
create mode 100644 Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
create mode 100644 Magnetodynamics/electromagnet.geo
create mode 100644 Magnetodynamics/electromagnet.pro
create mode 100644 Magnetodynamics/electromagnet_common.pro
create mode 100644 Magnetodynamics/transfo.geo
create mode 100644 Magnetodynamics/transfo.pro
create mode 100644 Magnetodynamics/transfo_common.pro
diff --git a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
new file mode 100644
index 0000000..feb65ed
--- /dev/null
+++ b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
@@ -0,0 +1,468 @@
+// This is a template .pro file containing a general formulation for 2D
+// magnetostatic and magnetodynamic problems in terms of the magnetic vector
+// potential a (potentially coupled with the electric scalar potential v), with
+// optional circuit coupling.
+
+// Below are definitions of the constants (inside "DefineConstant"), groups
+// (inside "DefineGroup") and functions (inside "DefineFunction") that can be
+// redefined from outside this template.
+
+DefineConstant[
+ Flag_FrequencyDomain = 1, // frequency-domain or time-domain simulation
+ Flag_CircuitCoupling = 0 // consider coupling with external electric circuit
+ CoefPower = 0.5, // coefficient for power calculations
+ Freq = 50, // frequency (for harmonic simulations)
+ TimeInit = 0, // intial time (for time-domain simulations)
+ TimeFinal = 1, // final time (for time-domain simulations)
+ DeltaTime = 0.01, // time step (for time-domain simulations)
+ InterpolationOrder = 1 // finite element order
+ Val_Rint = 0, // interior radius of annulus shell transformation region (VolInf_Mag)
+ Val_Rext = 0 // exterior radius of annulus shell transformation region (VolInf_Mag)
+];
+
+Group {
+ DefineGroup[
+ VolCC_Mag, // the non-conducting part
+ VolC_Mag, // the conducting part
+ VolV_Mag, // a moving conducting part, with invariant mesh
+ VolM_Mag, // permanent magnets
+ VolS0_Mag, // current source domain with imposed current densities js0
+ VolS_Mag, // current source domain with imposed current, imposed voltage or
+ // circuit coupling
+ VolInf_Mag, // annulus where a infinite shell transformation is applied
+ SurFluxTube_Mag, // boundary with Neumann BC
+ SurPerfect_Mag, // boundary of perfect conductors (i.e. non-meshed)
+ SurImped_Mag // boundary of conductors approximated by a surface impedance
+ // (i.e. non-meshed)
+ ];
+ If(Flag_CircuitCoupling)
+ DefineGroup[
+ SourceV_Cir, // voltage sources
+ SourceI_Cir, // current sources
+ Resistance_Cir, // resistances (linear)
+ Inductance_Cir, // inductances
+ Capacitance_Cir, // capacitances
+ Diode_Cir // diodes (nonlinear resistances)
+ ];
+ EndIf
+}
+
+Function {
+ DefineFunction[
+ nu, // reluctivity (in Vol_Mag)
+ sigma, // conductivity (in VolC_Mag and VolS_Mag)
+ br, // remanent magnetic flux density (in VolM_Mag)
+ js0, // source current density (in VolS0_Mag)
+ nxh, // n x magnetic field (on SurFluxTube_Mag)
+ Velocity, // velocity of moving part VolV_Moving
+ Ns, // number of turns (in VolS_Mag)
+ Sc, // cross-section of windings (in VolS_Mag)
+ CoefGeos, // geometrical coefficient for 2D or 2D axi model
+ Ysur // surface admittance (inverse of surface impedance Zsur) on SurImped_Mag
+ ];
+ If(Flag_CircuitCoupling)
+ DefineFunction[
+ Resistance, // resistance values
+ Inductance, // inductance values
+ Capacitance // capacitance values
+ ];
+ EndIf
+}
+
+// End of definitions.
+
+Group{
+ // all volumes
+ Vol_Mag = Region[ {VolCC_Mag, VolC_Mag, VolV_Mag, VolM_Mag, VolS0_Mag, VolS_Mag} ];
+ // all volumes + surfaces on which integrals will be computed
+ VolWithSur_Mag = Region[ {Vol_Mag, SurFluxTube_Mag, SurPerfect_Mag, SurImped_Mag} ];
+ If(Flag_CircuitCoupling)
+ // all circuit impedances
+ DomainZ_Cir = Region[ {Resistance_Cir, Inductance_Cir, Capacitance_Cir} ];
+ // all circuit sources
+ DomainSource_Cir = Region[ {SourceV_Cir, SourceI_Cir} ];
+ // all circuit elements
+ Domain_Cir = Region[ {DomainZ_Cir, DomainSource_Cir} ];
+ EndIf
+}
+
+Jacobian {
+ { Name Vol;
+ Case {
+ { Region VolInf_Mag ;
+ Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
+ { Region All; Jacobian Vol; }
+ }
+ }
+ { Name Sur;
+ Case {
+ { Region All; Jacobian Sur; }
+ }
+ }
+}
+
+Integration {
+ { Name Gauss_v;
+ Case {
+ { Type Gauss;
+ Case {
+ { GeoElement Point; NumberOfPoints 1; }
+ { GeoElement Line; NumberOfPoints 5; }
+ { GeoElement Triangle; NumberOfPoints 7; }
+ { GeoElement Quadrangle; NumberOfPoints 4; }
+ { GeoElement Tetrahedron; NumberOfPoints 15; }
+ { GeoElement Hexahedron; NumberOfPoints 14; }
+ { GeoElement Prism; NumberOfPoints 21; }
+ }
+ }
+ }
+ }
+}
+
+// Same FunctionSpace for both static and dynamic formulations
+FunctionSpace {
+ { Name Hcurl_a_2D; Type Form1P; // 1-form (circulations) on edges
+ // perpendicular to the plane of study
+ BasisFunction {
+ // \vec{a}(x) = \sum_{n \in N(Domain)} a_n \vec{s}_n(x)
+ // without nodes on perfect conductors (where a is constant)
+ { Name s_n; NameOfCoef a_n; Function BF_PerpendicularEdge;
+ Support VolWithSur_Mag; Entity NodesOf[All, Not SurPerfect_Mag]; }
+
+ // global basis function on boundary of perfect conductors
+ { Name s_skin; NameOfCoef a_skin; Function BF_GroupOfPerpendicularEdges;
+ Support VolWithSur_Mag; Entity GroupsOfNodesOf[SurPerfect_Mag]; }
+
+ // additional basis functions for 2nd order interpolation
+ If(InterpolationOrder == 2)
+ { Name s_e; NameOfCoef a_e; Function BF_PerpendicularEdge_2E;
+ Support Vol_Mag; Entity EdgesOf[All]; }
+ EndIf
+ }
+ GlobalQuantity {
+ { Name A; Type AliasOf; NameOfCoef a_skin; }
+ { Name I; Type AssociatedWith; NameOfCoef a_skin; }
+ }
+ Constraint {
+ { NameOfCoef a_n;
+ EntityType NodesOf; NameOfConstraint MagneticVectorPotential_2D; }
+
+ { NameOfCoef I;
+ EntityType GroupsOfNodesOf; NameOfConstraint Current_2D; }
+
+ If(InterpolationOrder == 2)
+ { NameOfCoef a_e;
+ EntityType EdgesOf; NameOfConstraint MagneticVectorPotential_2D_0; }
+ EndIf
+ }
+ }
+}
+
+FunctionSpace {
+ // Gradient of Electric scalar potential (2D)
+ { Name Hregion_u_2D; Type Form1P; // same as for \vec{a}
+ BasisFunction {
+ { Name sr; NameOfCoef ur; Function BF_RegionZ;
+ // constant vector (over the region) with nonzero z-component only
+ Support Region[{VolC_Mag, SurImped_Mag}];
+ Entity Region[{VolC_Mag, SurImped_Mag}]; }
+ }
+ 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; }
+ }
+ }
+
+ // Current in stranded coil (2D)
+ { Name Hregion_i_2D; Type Vector;
+ BasisFunction {
+ { Name sr; NameOfCoef ir; Function BF_RegionZ;
+ Support VolS_Mag; Entity VolS_Mag; }
+ }
+ 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; }
+ }
+ }
+}
+
+If(Flag_CircuitCoupling)
+ // UZ and IZ for impedances
+ FunctionSpace {
+ { Name Hregion_Z; Type Scalar;
+ BasisFunction {
+ { Name sr; NameOfCoef ir; Function BF_Region;
+ Support Domain_Cir; Entity Domain_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; }
+ }
+ }
+ }
+EndIf
+
+
+// Static Formulation
+Formulation {
+ { Name MagSta_a_2D; Type FemEquation;
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+ { Name ir; Type Local; NameOfSpace Hregion_i_2D; }
+ }
+ Equation {
+ Galerkin { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ -nu[] * br[] , {d a} ];
+ In VolM_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ -js0[] , {a} ];
+ In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
+ In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ nxh[] , {a} ];
+ In SurFluxTube_Mag; Jacobian Sur; Integration Gauss_v; }
+ }
+ }
+}
+
+// Dynamic Formulation (eddy currents)
+Formulation {
+ { Name MagDyn_a_2D; Type FemEquation;
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+ { Name A_floating; Type Global; NameOfSpace Hcurl_a_2D [A]; }
+ { Name I_perfect; Type Global; NameOfSpace Hcurl_a_2D [I]; }
+
+ { 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]; }
+
+ If(Flag_CircuitCoupling)
+ { Name Uz; Type Global; NameOfSpace Hregion_Z [Uz]; }
+ { Name Iz; Type Global; NameOfSpace Hregion_Z [Iz]; }
+ EndIf
+ }
+ Equation {
+ Galerkin { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
+ Galerkin { [ -nu[] * br[] , {d a} ];
+ In VolM_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ // Electric field e = -Dt[{a}]-{ur},
+ // with {ur} = Grad v constant in each region of VolC
+ Galerkin { DtDof [ sigma[] * Dof{a} , {a} ];
+ In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+ Galerkin { [ sigma[] * Dof{ur} / CoefGeos[] , {a} ];
+ In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ - sigma[] * (Velocity[] /\ Dof{d a}) , {a} ];
+ In VolV_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ -js0[] , {a} ];
+ In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { [ nxh[] , {a} ];
+ In SurFluxTube_Mag; Jacobian Sur; Integration Gauss_v; }
+
+ Galerkin { DtDof [ Ysur[] * Dof{a} , {a} ];
+ In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+ Galerkin { [ Ysur[] * Dof{ur} / CoefGeos[] , {a} ];
+ In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+
+ // When {ur} act as a test function, one obtains the circuits relations,
+ // relating the voltage and the current of each region in VolC
+ Galerkin { DtDof [ sigma[] * Dof{a} , {ur} ];
+ In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+ Galerkin { [ sigma[] * Dof{ur} / CoefGeos[] , {ur} ];
+ In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+ GlobalTerm { [ Dof{I} *(CoefGeos[]/Fabs[CoefGeos[]]) , {U} ]; In VolC_Mag; }
+
+ Galerkin { DtDof [ Ysur[] * Dof{a} , {ur} ];
+ In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+ Galerkin { [ Ysur[] * Dof{ur} / CoefGeos[] , {ur} ];
+ In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+ GlobalTerm { [ Dof{I} *(CoefGeos[]/Fabs[CoefGeos[]]) , {U} ]; In SurImped_Mag; }
+
+ // js[0] should be of the form: Ns[]/Sc[] * Vector[0,0,1]
+ Galerkin { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
+ In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Galerkin { DtDof [ Ns[]/Sc[] * Dof{a} , {ir} ];
+ In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+ Galerkin { [ Ns[]/Sc[] / sigma[] * (js0[]*Vector[0,0,1]) * Dof{ir} , {ir} ];
+ In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+ GlobalTerm { [ Dof{Us} / CoefGeos[] , {Is} ]; In VolS_Mag; }
+ // Attention: CoefGeo[.] = 2*Pi for Axi
+
+ GlobalTerm { [ -Dof{I_perfect} , {A_floating} ]; In SurPerfect_Mag; }
+
+ If(Flag_CircuitCoupling)
+ 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 { NeverDt[ Dof{Iz} , {Iz} ]; In Capacitance_Cir; }
+ GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ]; In Capacitance_Cir; }
+
+ GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ]; In Diode_Cir; }
+ GlobalTerm { NeverDt[ Resistance[{Uz}] * Dof{Iz} , {Iz} ]; In Diode_Cir; }
+
+ GlobalTerm { [ 0. * Dof{Iz} , {Iz} ]; In DomainSource_Cir; }
+
+ GlobalEquation {
+ Type Network; NameOfConstraint ElectricalCircuit;
+ { Node {I}; Loop {U}; Equation {I}; In VolC_Mag; }
+ { Node {Is}; Loop {Us}; Equation {Us}; In VolS_Mag; }
+ { Node {Iz}; Loop {Uz}; Equation {Uz}; In Domain_Cir; }
+ }
+ EndIf
+
+ }
+ }
+}
+
+Resolution {
+ { Name MagDyn_a_2D;
+ System {
+ { Name Sys; NameOfFormulation MagDyn_a_2D;
+ If(Flag_FrequencyDomain)
+ Type ComplexValue; Frequency Freq;
+ EndIf
+ }
+ }
+ Operation {
+ If(Flag_FrequencyDomain)
+ Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+ Else
+ InitSolution[Sys]; // provide initial condition
+ TimeLoopTheta[Time0, TimeMax, DeltaTime, 1.]{
+ // Euler implicit (1) -- Crank-Nicolson (0.5)
+ Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+ }
+ EndIf
+ }
+ }
+ { Name MagSta_a_2D;
+ System {
+ { Name Sys; NameOfFormulation MagSta_a_2D; }
+ }
+ Operation {
+ Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+ }
+ }
+}
+
+// Same PostProcessing for both static and dynamic formulations (both refer to
+// the same FunctionSpace from which the solution is obtained)
+PostProcessing {
+ { Name MagDyn_a_2D; NameOfFormulation MagDyn_a_2D;
+ PostQuantity {
+ // In 2D, a is a vector with only a z-component: (0,0,az)
+ { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+ // The equilines of az are field lines (giving the magnetic field direction)
+ { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name h; Value {
+ Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; }
+ Term { [ -nu[] * br[] ]; In VolM_Mag; Jacobian Vol; }
+ }
+ }
+ { Name js; Value {
+ Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+ Term { [ (js0[]*Vector[0,0,1])*{ir} ]; In VolS_Mag; Jacobian Vol; }
+ Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; } // to force a vector result out of sources
+ }
+ }
+ { Name j;
+ // Only correct in MagDyn
+ Value {
+ Term { [ -sigma[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In VolC_Mag; Jacobian Vol; }
+ Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+ Term { [ (js0[]*Vector[0,0,1])*{ir} ]; In VolS_Mag; Jacobian Vol; }
+ Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }
+ // Current density in A/m
+ Term { [ -Ysur[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In SurImped_Mag; Jacobian Sur; }
+ }
+ }
+
+ { Name JouleLosses;
+ Value {
+ Integral { [ CoefPower * sigma[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
+ In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+ Integral { [ CoefPower * 1./sigma[]*SquNorm[js0[]] ];
+ In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Integral { [ CoefPower * 1./sigma[]*SquNorm[(js0[]*Vector[0,0,1])*{ir}] ];
+ In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+ Integral { [ CoefPower * Ysur[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
+ In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+ }
+ }
+
+ { Name U; Value {
+ Term { [ {U} ]; In VolC_Mag; }
+ Term { [ {Us} ]; In VolS_Mag; }
+ If(Flag_CircuitCoupling)
+ Term { [ {Uz} ]; In Domain_Cir; }
+ EndIf
+ }
+ }
+
+ { Name I; Value {
+ Term { [ {I} ]; In VolC_Mag; }
+ Term { [ {Is} ]; In VolS_Mag; }
+ If(Flag_CircuitCoupling)
+ Term { [ {Iz} ]; In Domain_Cir; }
+ EndIf
+ }
+ }
+
+ }
+ }
+
+ { Name MagSta_a_2D; NameOfFormulation MagSta_a_2D;
+ PostQuantity {
+ // In 2D, a is a vector with only a z-component: (0,0,az)
+ { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+ // The equilines of az are field lines (giving the magnetic field direction)
+ { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name h; Value { Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name j;
+ Value {
+ Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+ Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }
+ }
+ }
+ }
+ }
+}
diff --git a/Magnetodynamics/electromagnet.geo b/Magnetodynamics/electromagnet.geo
new file mode 100644
index 0000000..7bd410c
--- /dev/null
+++ b/Magnetodynamics/electromagnet.geo
@@ -0,0 +1,54 @@
+/* -------------------------------------------------------------------
+ File "electromagnet.geo"
+
+ This file is the geometrical description used by GMSH to produce
+ the file "electromagnet.msh".
+ ------------------------------------------------------------------- */
+
+dxCore = 50.e-3; dyCore = 100.e-3;
+xInd = 75.e-3; dxInd = 25.e-3; dyInd = 50.e-3;
+
+Include "electromagnet_common.pro";
+
+s = DefineNumber[1, Name "Model parameters/Global mesh size",
+ Help "Reduce for finer mesh"];
+p0 = 12.e-3 *s;
+pCorex = 4.e-3 *s; pCorey0 = 8.e-3 *s; pCorey = 4.e-3 *s;
+pIndx = 5.e-3 *s; pIndy = 5.e-3 *s;
+pInt = 12.5e-3*s; pExt = 12.5e-3*s;
+
+Point(1) = {0,0,0,p0};
+Point(2) = {dxCore,0,0,pCorex};
+Point(3) = {dxCore,dyCore,0,pCorey};
+Point(4) = {0,dyCore,0,pCorey0};
+Point(5) = {xInd,0,0,pIndx};
+Point(6) = {xInd+dxInd,0,0,pIndx};
+Point(7) = {xInd+dxInd,dyInd,0,pIndy};
+Point(8) = {xInd,dyInd,0,pIndy};
+Point(9) = {rInt,0,0,pInt};
+Point(10) = {rExt,0,0,pExt};
+Point(11) = {0,rInt,0,pInt};
+Point(12) = {0,rExt,0,pExt};
+
+Line(1) = {1,2}; Line(2) = {2,5}; Line(3) = {5,6};
+Line(4) = {6,9}; Line(5) = {9,10}; Line(6) = {1,4};
+Line(7) = {4,11}; Line(8) = {11,12}; Line(9) = {2,3};
+Line(10) = {3,4}; Line(11) = {6,7}; Line(12) = {7,8};
+Line(13) = {8,5};
+
+Circle(14) = {9,1,11}; Circle(15) = {10,1,12};
+
+Line Loop(16) = {-6,1,9,10}; Plane Surface(17) = {16};
+Line Loop(18) = {11,12,13,3}; Plane Surface(19) = {18};
+Line Loop(20) = {7,-14,-4,11,12,13,-2,9,10}; Plane Surface(21) = {-20}; // "-" for orientation
+Line Loop(22) = {8,-15,-5,14}; Plane Surface(23) = {-22};
+
+Physical Surface(101) = {21}; /* Air */
+Physical Surface(102) = {17}; /* Core */
+Physical Surface(103) = {19}; /* Ind */
+Physical Surface(111) = {23}; /* AirInf */
+
+Physical Line(1100) = {1,2,3,4,5}; /* Surface ht = 0 */
+Physical Line(1101) = {6,7,8}; /* Surface bn = 0 */
+Physical Line(1102) = {15}; /* Surface Inf */
+
diff --git a/Magnetodynamics/electromagnet.pro b/Magnetodynamics/electromagnet.pro
new file mode 100644
index 0000000..eb3d8a9
--- /dev/null
+++ b/Magnetodynamics/electromagnet.pro
@@ -0,0 +1,77 @@
+/* -------------------------------------------------------------------
+ Tutorial 7a : magnetostatic field of an electromagnet, bis
+
+ Features:
+ - Same as Tutorial 2, but using a generic template formulation library
+
+ To compute the solution in a terminal:
+ getdp electromagnet -solve MagSta_a_2D
+ getdp electromagnet -pos Map_a
+
+ To compute the solution interactively from the Gmsh GUI:
+ File > Open > electromagnet.pro
+ Run (button at the bottom of the left panel)
+ ------------------------------------------------------------------- */
+
+Include "electromagnet_common.pro";
+
+Group {
+ // Physical regions
+ Air = Region[ 101 ]; Core = Region[ 102 ];
+ Ind = Region[ 103 ]; AirInf = Region[ 111 ];
+ Surface_ht0 = Region[ 1100 ];
+ Surface_bn0 = Region[ 1101 ];
+ Surface_Inf = Region[ 1102 ];
+
+ // Abstract regions used in the "Lib_MagStaDyn_av_2D_Cir.pro" template file
+ // that is included below:
+ VolCC_Mag = Region[{Air, Core, AirInf}]; // Non-conducting regions
+ VolS_Mag = Region[{Ind}]; // Stranded conductors, i.e., coils
+ VolInf_Mag = Region[{AirInf}]; // annulus for infinite shell transformation
+ Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of annulus
+}
+
+Function {
+ DefineConstant[
+ murCore = {100, Name "Model parameters/Mur core"},
+ Current = {0.01, Name "Model parameters/Current"}
+ ];
+
+ mu0 = 4.e-7 * Pi;
+ nu[ Region[{Air, Ind, AirInf}] ] = 1. / mu0;
+ nu[ Core ] = 1. / (murCore * mu0);
+
+ Ns[ Ind ] = 1000 ; // number of turns in coil
+ Sc[ Ind ] = SurfaceArea[] ; // surface (cross section) of coil
+ // Current density in each coil portion for a unit current (will be multiplied
+ // by the actual total current in the coil)
+ js0[ Ind ] = Ns[]/Sc[] * Vector[0,0,-1];
+}
+
+Constraint {
+ { Name MagneticVectorPotential_2D;
+ Case {
+ { Region Surface_bn0; Value 0; }
+ { Region Surface_Inf; Value 0; }
+ }
+ }
+ { Name Current_2D;
+ Case {
+ { Region Ind; Value Current; }
+ }
+ }
+ { Name Voltage_2D;
+ Case {
+ }
+ }
+}
+
+Include "Lib_MagStaDyn_av_2D_Cir.pro";
+
+PostOperation {
+ { Name Map_a; NameOfPostProcessing MagSta_a_2D;
+ Operation {
+ Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
+ }
+ }
+}
diff --git a/Magnetodynamics/electromagnet_common.pro b/Magnetodynamics/electromagnet_common.pro
new file mode 100644
index 0000000..56879cd
--- /dev/null
+++ b/Magnetodynamics/electromagnet_common.pro
@@ -0,0 +1,4 @@
+// Parameters shared by Gmsh and GetDP
+
+rInt = 200.e-3;
+rExt = 250.e-3;
diff --git a/Magnetodynamics/transfo.geo b/Magnetodynamics/transfo.geo
new file mode 100644
index 0000000..ae9875b
--- /dev/null
+++ b/Magnetodynamics/transfo.geo
@@ -0,0 +1,173 @@
+
+Include "transfo_common.dat";
+
+Solver.AutoShowLastStep = 0; // don't automatically show the last time step
+
+// CORE
+
+p_Leg_1_L_0=newp; Point(newp) = {-width_Core/2., 0, 0, c_Core};
+p_Leg_1_R_0=newp; Point(newp) = {-width_Core/2.+width_Core_Leg, 0, 0, c_Core};
+
+p_Leg_1_L_1=newp; Point(newp) = {-width_Core/2., height_Core/2., 0, c_Core};
+p_Leg_1_R_1=newp; Point(newp) = {-width_Core/2.+width_Core_Leg, height_Core/2.-width_Core_Leg, 0, c_Core};
+
+
+p_Leg_2_L_0=newp; Point(newp) = {width_Core/2.-width_Core_Leg, 0, 0, c_Core};
+p_Leg_2_R_0=newp; Point(newp) = {width_Core/2., 0, 0, c_Core};
+
+p_Leg_2_L_1=newp; Point(newp) = {width_Core/2.-width_Core_Leg, height_Core/2.-width_Core_Leg, 0, c_Core};
+p_Leg_2_R_1=newp; Point(newp) = {width_Core/2., height_Core/2., 0, c_Core};
+
+
+l_Core_In[]={};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_1_R_0, p_Leg_1_R_1};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_1_R_1, p_Leg_2_L_1};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_2_L_1, p_Leg_2_L_0};
+
+l_Core_Out[]={};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_2_R_0, p_Leg_2_R_1};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_2_R_1, p_Leg_1_L_1};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_1_L_1, p_Leg_1_L_0};
+
+l_Core_Leg_1_Y0[]={};
+l_Core_Leg_1_Y0[]+=newl; Line(newl) = {p_Leg_1_L_0, p_Leg_1_R_0};
+
+l_Core_Leg_2_Y0[]={};
+l_Core_Leg_2_Y0[]+=newl; Line(newl) = {p_Leg_2_L_0, p_Leg_2_R_0};
+
+
+ll_Core=newll; Line Loop(newll) = {l_Core_In[], l_Core_Leg_2_Y0[], l_Core_Out[], l_Core_Leg_1_Y0[]};
+s_Core=news; Plane Surface(news) = {ll_Core};
+
+Physical Surface("CORE", CORE) = {s_Core};
+
+// COIL_1_PLUS
+
+x_[]=Point{p_Leg_1_R_0};
+p_Coil_1_p_int_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_p_ext_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_1+width_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_p_int_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+p_Coil_1_p_ext_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_1+width_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+
+l_Coil_1_p[]={};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_ext_0, p_Coil_1_p_ext_1};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_ext_1, p_Coil_1_p_int_1};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_int_1, p_Coil_1_p_int_0};
+
+l_Coil_1_p_Y0[]={};
+l_Coil_1_p_Y0[]+=newl; Line(newl) = {p_Coil_1_p_int_0, p_Coil_1_p_ext_0};
+
+ll_Coil_1_p=newll; Line Loop(newll) = {l_Coil_1_p[], l_Coil_1_p_Y0[]};
+s_Coil_1_p=news; Plane Surface(news) = {ll_Coil_1_p};
+
+Physical Surface("COIL_1_PLUS", COIL_1_PLUS) = {s_Coil_1_p};
+
+
+// COIL_1_MINUS
+
+x_[]=Point{p_Leg_1_L_0};
+p_Coil_1_m_int_0=newp; Point(newp) = {x_[0]-gap_Core_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_m_ext_0=newp; Point(newp) = {x_[0]-(gap_Core_Coil_1+width_Coil_1), 0, 0, c_Coil_1};
+p_Coil_1_m_int_1=newp; Point(newp) = {x_[0]-gap_Core_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+p_Coil_1_m_ext_1=newp; Point(newp) = {x_[0]-(gap_Core_Coil_1+width_Coil_1), height_Coil_1/2, 0, c_Coil_1};
+
+l_Coil_1_m[]={};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_ext_0, p_Coil_1_m_ext_1};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_ext_1, p_Coil_1_m_int_1};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_int_1, p_Coil_1_m_int_0};
+
+l_Coil_1_m_Y0[]={};
+l_Coil_1_m_Y0[]+=newl; Line(newl) = {p_Coil_1_m_int_0, p_Coil_1_m_ext_0};
+
+ll_Coil_1_m=newll; Line Loop(newll) = {l_Coil_1_m[], l_Coil_1_m_Y0[]};
+s_Coil_1_m=news; Plane Surface(news) = {ll_Coil_1_m};
+
+Physical Surface("COIL_1_MINUS", COIL_1_MINUS) = {s_Coil_1_m};
+
+
+// COIL_2_PLUS
+
+x_[]=Point{p_Leg_2_R_0};
+p_Coil_2_p_int_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_p_ext_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_2+width_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_p_int_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+p_Coil_2_p_ext_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_2+width_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+
+l_Coil_2_p[]={};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_ext_0, p_Coil_2_p_ext_1};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_ext_1, p_Coil_2_p_int_1};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_int_1, p_Coil_2_p_int_0};
+
+l_Coil_2_p_Y0[]={};
+l_Coil_2_p_Y0[]+=newl; Line(newl) = {p_Coil_2_p_int_0, p_Coil_2_p_ext_0};
+
+ll_Coil_2_p=newll; Line Loop(newll) = {l_Coil_2_p[], l_Coil_2_p_Y0[]};
+s_Coil_2_p=news; Plane Surface(news) = {ll_Coil_2_p};
+
+Physical Surface("COIL_2_PLUS", COIL_2_PLUS) = {s_Coil_2_p};
+
+
+// COIL_2_MINUS
+
+x_[]=Point{p_Leg_2_L_0};
+p_Coil_2_m_int_0=newp; Point(newp) = {x_[0]-gap_Core_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_m_ext_0=newp; Point(newp) = {x_[0]-(gap_Core_Coil_2+width_Coil_2), 0, 0, c_Coil_2};
+p_Coil_2_m_int_1=newp; Point(newp) = {x_[0]-gap_Core_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+p_Coil_2_m_ext_1=newp; Point(newp) = {x_[0]-(gap_Core_Coil_2+width_Coil_2), height_Coil_2/2, 0, c_Coil_2};
+
+l_Coil_2_m[]={};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_ext_0, p_Coil_2_m_ext_1};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_ext_1, p_Coil_2_m_int_1};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_int_1, p_Coil_2_m_int_0};
+
+l_Coil_2_m_Y0[]={};
+l_Coil_2_m_Y0[]+=newl; Line(newl) = {p_Coil_2_m_int_0, p_Coil_2_m_ext_0};
+
+ll_Coil_2_m=newll; Line Loop(newll) = {l_Coil_2_m[], l_Coil_2_m_Y0[]};
+s_Coil_2_m=news; Plane Surface(news) = {ll_Coil_2_m};
+
+Physical Surface("COIL_2_MINUS", COIL_2_MINUS) = {s_Coil_2_m};
+
+
+// AIR_WINDOW
+
+l_Air_Window_Y0[]={};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Leg_1_R_0, p_Coil_1_p_int_0};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Coil_1_p_ext_0, p_Coil_2_m_ext_0};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Coil_2_m_int_0, p_Leg_2_L_0};
+
+ll_Air_Window=newll; Line Loop(newll) = {-l_Core_In[], -l_Coil_1_p[], l_Coil_2_m[], l_Air_Window_Y0[]};
+s_Air_Window=news; Plane Surface(news) = {ll_Air_Window};
+
+Physical Surface("AIR_WINDOW", AIR_WINDOW) = {s_Air_Window};
+
+
+// AIR_EXT
+
+x_[]=Point{p_Leg_2_R_1};
+p_Air_Ext_1_R_0=newp; Point(newp) = {x_[0]+gap_Core_Box_X, 0, 0, c_Box};
+p_Air_Ext_1_R_1=newp; Point(newp) = {x_[0]+gap_Core_Box_X, x_[1]+gap_Core_Box_Y, 0, c_Box};
+
+x_[]=Point{p_Leg_1_L_1};
+p_Air_Ext_1_L_0=newp; Point(newp) = {x_[0]-gap_Core_Box_X, 0, 0, c_Box};
+p_Air_Ext_1_L_1=newp; Point(newp) = {x_[0]-gap_Core_Box_X, x_[1]+gap_Core_Box_Y, 0, c_Box};
+
+l_Air_Ext[]={};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_R_0, p_Air_Ext_1_R_1};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_R_1, p_Air_Ext_1_L_1};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_L_1, p_Air_Ext_1_L_0};
+
+
+l_Air_Ext_Y0[]={};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Leg_2_R_0, p_Coil_2_p_int_0};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Coil_2_p_ext_0, p_Air_Ext_1_R_0};
+
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Air_Ext_1_L_0, p_Coil_1_m_ext_0};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Coil_1_m_int_0, p_Leg_1_L_0};
+
+
+ll_Air_Ext=newll; Line Loop(newll) = {-l_Core_Out[], -l_Coil_2_p[], l_Coil_1_m[], l_Air_Ext[], l_Air_Ext_Y0[]};
+s_Air_Ext=news; Plane Surface(news) = {ll_Air_Ext};
+
+Physical Surface("AIR_EXT", AIR_EXT) = {s_Air_Ext};
+Physical Line("SUR_AIR_EXT", SUR_AIR_EXT) = {l_Air_Ext[]};
diff --git a/Magnetodynamics/transfo.pro b/Magnetodynamics/transfo.pro
new file mode 100644
index 0000000..5c1def3
--- /dev/null
+++ b/Magnetodynamics/transfo.pro
@@ -0,0 +1,260 @@
+/* -------------------------------------------------------------------
+ Tutorial 7b : magnetodyamic model of a single-phase transformer
+
+ Features:
+ - Time-domain and frequency-domain dynamical problems
+ - Use of a generic template formulation library
+ - Circuit coupling used as a black-box (see Tutorial 8 for more)
+
+ To compute the solution in a terminal:
+ getdp transfo -solve MagDyn_a_2D
+ getdp electromagnet -pos Map
+
+ To compute the solution interactively from the Gmsh GUI:
+ File > Open > transfo.pro
+ Run (button at the bottom of the left panel)
+ ------------------------------------------------------------------- */
+
+Include "transfo_common.pro";
+
+DefineConstant[
+ type_Conds = {2, Choices{1 = "Massive", 2 = "Coil"}, Highlight "Blue",
+ Name "Parameters/01Conductor type"}
+ type_Source = {2, Choices{1 = "Current", 2 = "Voltage"}, Highlight "Blue",
+ Name "Parameters/02Source type"}
+ type_Analysis = {1, Choices{1 = "Frequency-domain", 2 = "Time-domain"}, Highlight "Blue",
+ Name "Parameters/03Analysis type"}
+ Freq = {50, Min 0, Max 1e3, Step 1,
+ Name "Parameters/Frequency"}
+];
+
+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 */
+
+ VolCC_Mag = Region[{}]; // Non-conducting regions
+ VolC_Mag = Region[{}]; // Massive conductors
+ VolS_Mag = Region[{}]; // Stranded conductors, i.e., coils
+
+ // air physical groups
+ Air = Region[{AIR_WINDOW, AIR_EXT}];
+ VolCC_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];
+ VolCC_Mag += Region[Core];
+
+ Coil_1_P = Region[COIL_1_PLUS];
+ Coil_1_M = Region[COIL_1_MINUS];
+ 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 = Region[{Coil_2_P, Coil_2_M}];
+
+ Coils = Region[{Coil_1, Coil_2}];
+
+ If (type_Conds == 1)
+ VolC_Mag += Region[{Coils}];
+ ElseIf (type_Conds == 2)
+ VolS_Mag += Region[{Coils}];
+ VolCC_Mag += Region[{Coils}];
+ EndIf
+}
+
+
+Function {
+ mu0 = 4e-7*Pi;
+
+ mu[Air] = 1 * mu0;
+
+ mur_Core = 1;
+ mu[Core] = mur_Core * mu0;
+
+ mu[Coils] = 1 * mu0;
+ sigma[Coils] = 1e7;
+
+ // For a correct definition of the voltage
+ CoefGeo = thickness_Core;
+
+ // To be defined separately for each coil portion
+ Sc[Coil_1_P] = SurfaceArea[];
+ SignBranch[Coil_1_P] = 1; // To fix the convention of positive current (1:
+ // along Oz, -1: along -Oz)
+
+ Sc[Coil_1_M] = SurfaceArea[];
+ SignBranch[Coil_1_M] = -1;
+
+ Sc[Coil_2_P] = SurfaceArea[];
+ SignBranch[Coil_2_P] = 1;
+
+ Sc[Coil_2_M] = SurfaceArea[];
+ SignBranch[Coil_2_M] = -1;
+
+ // Number of turns (same for PLUS and MINUS portions) (half values because
+ // half coils are defined)
+ Ns[Coil_1] = 1;
+ Ns[Coil_2] = 1;
+
+ // Global definitions (nothing to change):
+
+ // Current density in each coil portion for a unit current (will be multiplied
+ // by the actual total current in the coil)
+ js0[Coils] = Ns[]/Sc[] * Vector[0,0,SignBranch[]];
+ CoefGeos[Coils] = SignBranch[] * CoefGeo;
+
+ // The reluctivity will be used
+ nu[] = 1/mu[];
+}
+
+If(type_Analysis == 1)
+ Flag_FrequencyDomain = 1;
+Else
+ Flag_FrequencyDomain = 0;
+EndIf
+
+If (type_Source == 1) // current
+
+ Flag_CircuitCoupling = 0;
+
+ElseIf (type_Source == 2) // voltage
+
+ // PLUS and MINUS coil portions to be connected in series, with applied
+ // voltage on the resulting branch
+ Flag_CircuitCoupling = 1;
+
+ // Here is the definition of the circuits on primary and secondary sides:
+ Group {
+ // Empty Groups to be filled
+ Resistance_Cir = Region[{}];
+ Inductance_Cir = Region[{}] ;
+ Capacitance_Cir = Region[{}] ;
+ SourceV_Cir = Region[{}]; // Voltage sources
+ SourceI_Cir = Region[{}]; // Current sources
+
+ // Primary side
+ E_in = Region[10001]; // arbitrary region number (not linked to the mesh)
+ SourceV_Cir += Region[{E_in}];
+
+ // Secondary side
+ R_out = Region[10101]; // arbitrary region number (not linked to the mesh)
+ Resistance_Cir += Region[{R_out}];
+ }
+
+ Function {
+ deg = Pi/180;
+ // Input RMS voltage (half of the voltage because of symmetry; half coils
+ // are defined)
+ val_E_in = 1.;
+ phase_E_in = 90 *deg; // Phase in radian (from phase in degree)
+ // High value for an open-circuit test; Low value for a short-circuit test;
+ // any value in-between for any charge
+ Resistance[R_out] = 1e6;
+ }
+
+ Constraint {
+
+ { Name Current_Cir ;
+ Case {
+ }
+ }
+
+ { Name Voltage_Cir ;
+ Case {
+ { Region E_in; Value val_E_in; TimeFunction F_Cos_wt_p[]{2*Pi*Freq, phase_E_in};}
+ }
+ }
+
+ { Name ElectricalCircuit ; Type Network ;
+ Case Circuit_1 {
+ // PLUS and MINUS coil portions to be connected in series, together with
+ // E_in (an additional resistance should be defined to represent the
+ // Coil_1 end-winding (not considered in the 2D model))
+ { Region E_in; Branch {1,2}; }
+
+ { Region Coil_1_P; Branch {2,3} ; }
+ { Region Coil_1_M; Branch {3,1} ; }
+ }
+ Case Circuit_2 {
+ // PLUS and MINUS coil portions to be connected in series, together with
+ // R_out (an additional resistance should be defined to represent the
+ // Coil_2 end-winding (not considered in the 2D model))
+ { Region R_out; Branch {1,2}; }
+
+ { Region Coil_2_P; Branch {2,3} ; }
+ { Region Coil_2_M; Branch {3,1} ; }
+ }
+ }
+
+ }
+
+EndIf
+
+Constraint {
+ { Name MagneticVectorPotential_2D;
+ Case {
+ { Region Sur_Air_Ext; Value 0; }
+ }
+ }
+ { Name Current_2D;
+ Case {
+ If (type_Source == 1)
+ // Current in each coil (same for PLUS and MINUS portions)
+ { Region Coil_1; Value 1; }
+ { Region Coil_2; Value 0; }
+ EndIf
+ }
+ }
+ { Name Voltage_2D;
+ Case {
+ }
+ }
+}
+
+Include "Lib_MagStaDyn_av_2D_Cir.pro";
+
+PostOperation {
+ { Name Map_a; NameOfPostProcessing MagDyn_a_2D;
+ Operation {
+ Print[ j, OnElementsOf Region[{VolC_Mag, VolS_Mag}], Format Gmsh, File "j.pos" ];
+ Print[ b, OnElementsOf Vol_Mag, Format Gmsh, File "b.pos" ];
+ Print[ az, OnElementsOf Vol_Mag, Format Gmsh, File "az.pos" ];
+
+ If (type_Source == 1) // current
+ // In text file UI.txt: voltage and current for each coil portion (note
+ // that the voltage is not equally distributed in PLUS and MINUS
+ // portions, which is the reason why we must apply the total voltage
+ // through a circuit -> type_Source == 2)
+ Echo[ "Coil_1_P", Format Table, File "UI.txt" ];
+ Print[ U, OnRegion Coil_1_P, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion Coil_1_P, Format FrequencyTable, File > "UI.txt"];
+ Echo[ "Coil_1_M", Format Table, File > "UI.txt" ];
+ Print[ U, OnRegion Coil_1_M, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion Coil_1_M, Format FrequencyTable, File > "UI.txt"];
+
+ Echo[ "Coil_2_P", Format Table, File > "UI.txt" ];
+ Print[ U, OnRegion Coil_2_P, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion Coil_2_P, Format FrequencyTable, File > "UI.txt"];
+ Echo[ "Coil_2_M", Format Table, File > "UI.txt" ];
+ Print[ U, OnRegion Coil_2_M, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion Coil_2_M, Format FrequencyTable, File > "UI.txt"];
+
+ ElseIf (type_Source == 2)
+ // In text file UI.txt: voltage and current of the primary coil (from E_in)
+ // (real and imaginary parts!)
+ Echo[ "E_in", Format Table, File "UI.txt" ];
+ Print[ U, OnRegion E_in, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion E_in, Format FrequencyTable, File > "UI.txt"];
+
+ // In text file UI.txt: voltage and current of the secondary coil (from R_out)
+ Echo[ "R_out", Format Table, File > "UI.txt" ];
+ Print[ U, OnRegion R_out, Format FrequencyTable, File > "UI.txt" ];
+ Print[ I, OnRegion R_out, Format FrequencyTable, File > "UI.txt"];
+ EndIf
+ }
+ }
+}
diff --git a/Magnetodynamics/transfo_common.pro b/Magnetodynamics/transfo_common.pro
new file mode 100644
index 0000000..a116935
--- /dev/null
+++ b/Magnetodynamics/transfo_common.pro
@@ -0,0 +1,49 @@
+
+// Dimensions
+
+width_Core = 1.;
+height_Core = 1.;
+
+// Thickness along Oz (to be considered for a correct definition of voltage)
+thickness_Core = 1.;
+
+width_Window = 0.5;
+height_Window = 0.5;
+
+width_Core_Leg = (width_Core-width_Window)/2.;
+
+width_Coil_1 = 0.10;
+height_Coil_1 = 0.25;
+gap_Core_Coil_1 = 0.05;
+
+width_Coil_2 = 0.10;
+height_Coil_2 = 0.25;
+gap_Core_Coil_2 = 0.05;
+
+gap_Core_Box_X = 1.;
+gap_Core_Box_Y = 1.;
+
+// Characteristic lenghts (for mesh sizes)
+
+s = 1;
+
+c_Core = width_Core_Leg/10. *s;
+
+c_Coil_1 = height_Coil_1/2/5 *s;
+c_Coil_2 = height_Coil_2/2/5 *s;
+
+c_Box = gap_Core_Box_X/6. *s;
+
+// Physical regions
+
+AIR_EXT = 1001;
+SUR_AIR_EXT = 1002;
+AIR_WINDOW = 1011;
+
+CORE = 1050;
+
+COIL_1_PLUS = 1101;
+COIL_1_MINUS = 1102;
+
+COIL_2_PLUS = 1201;
+COIL_2_MINUS = 1202;
diff --git a/Magnetostatics/electromagnet.pro b/Magnetostatics/electromagnet.pro
index d582938..629379f 100644
--- a/Magnetostatics/electromagnet.pro
+++ b/Magnetostatics/electromagnet.pro
@@ -8,7 +8,7 @@
To compute the solution in a terminal:
getdp electromagnet -solve MagSta_a
- getdp electromagnet -pos Map
+ getdp electromagnet -pos Map_a
To compute the solution interactively from the Gmsh GUI:
File > Open > electromagnet.pro
@@ -66,13 +66,13 @@ Group {
The abstract regions in this model have the following interpretation:
- Vol_Nu_Mag = region where the term [ nu[] * Dof{d a} , {d a} ] is assembled
- Vol_Js_Mag = region where the term [ - Dof{js} , {a} ] is assembled
- - Vol_Inf = region where the infinite ring geometric transformation is applied
+ - Vol_Inf_Mag = region where the infinite ring geometric transformation is applied
- Sur_Dir_Mag = Homogeneous Dirichlet part of the model's boundary;
- Sur_Neu_Mag = Homogeneous Neumann part of the model's boundary;
*/
Vol_Nu_Mag = Region[ {Air, AirInf, Core, Ind} ];
Vol_Js_Mag = Region[ Ind ];
- Vol_Inf = Region[ {AirInf} ];
+ Vol_Inf_Mag = Region[ {AirInf} ];
Sur_Dir_Mag = Region[ {Surface_bn0, Surface_Inf} ];
Sur_Neu_Mag = Region[ {Surface_ht0} ];
}
@@ -168,7 +168,7 @@ Include "electromagnet_common.pro";
Val_Rint = rInt; Val_Rext = rExt;
Jacobian {
{ Name Vol ;
- Case { { Region Vol_Inf ;
+ Case { { Region Vol_Inf_Mag ;
Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
{ Region All ; Jacobian Vol ; }
}
--
GitLab