diff --git a/Electrostatics/microstrip.geo b/Electrostatics/microstrip.geo
index 04ecbf515d566ca77ed14566b1be84e24bb52706..51dbe5d123487014117677e94377c8638b9d0abf 100644
--- a/Electrostatics/microstrip.geo
+++ b/Electrostatics/microstrip.geo
@@ -15,7 +15,8 @@ xBox = w/2. * 6. ; yBox = h * 12. ;
/* Definition of parameters for local mesh dimensions */
-s = 1. ;
+DefineConstant[ s = {1., Name "Parameters/Global mesh size factor"} ] ;
+
p0 = h / 10. * s ;
pLine0 = w/2. / 10. * s ; pLine1 = w/2. / 50. * s ;
pxBox = xBox / 10. * s ; pyBox = yBox / 8. * s ;
diff --git a/Magnetodynamics/electromagnet.pro b/Magnetodynamics/electromagnet.pro
index a693c6469079269d9d3156d4688c4e958b661945..22fdf163c5e7c63fa1a275e947bb8606eddd19a2 100644
--- a/Magnetodynamics/electromagnet.pro
+++ b/Magnetodynamics/electromagnet.pro
@@ -4,14 +4,18 @@
Features:
- Use of a template formulation library
- Identical to Tutorial 2 for a static current source
- - Frequency-domain solution (phasor) for a dynamic current source
+ - Frequency-domain or time-domain solution for a time-dependent
+ current source
To compute the static solution in a terminal:
getdp electromagnet -solve Magnetostatics2D_a -pos Map_a
- To compute the time-harmonic dynamic solution in a terminal:
+ To compute the frequency-domain solution in a terminal:
getdp electromagnet -solve Magnetodynamics2D_av -pos Map_a
+ To compute the time-dependent dynamic solution in a terminal:
+ getdp electromagnet -setnumber TimeDomain 1 -solve Magnetodynamics2D_av -pos Map_a
+
To compute the solution interactively from the Gmsh GUI:
File > Open > electromagnet.pro
You may choose the Resolution in the left panel:
@@ -42,10 +46,27 @@ Function {
DefineConstant[
murCore = {100, Name "Model parameters/Mur core"},
Current = {0.01, Name "Model parameters/Current"},
- frequency = {1, Name "Model parameters/Frequency"}
+ TimeDomain = {0, Choices{0 = "Frequency-domain", 1 = "Time-domain"},
+ Name "Model parameters/03Analysis type"}
+ frequency = {50, Visible !TimeDomain,
+ Name "Model parameters/Frequency" }
];
- Freq = frequency;
+ If(TimeDomain)
+ // Fix parameters from the "Lib_Magnetodynamics2D_av_Cir.pro" template:
+ Flag_FrequencyDomain = 0;
+ TimeInit = 0; // start simulation at time = 0s
+ TimeFinal = 20e-3; // stop simulation at time = 20 ms
+ DeltaTime = 1e-3; // use time steps equal to 1 ms
+ // Define the time modulation of the current source, i.e. a linear ramp from
+ // 0 to 1 until 10 ms, then a constant value of 1:
+ myModulation[] = ($Time < TimeFinal / 2) ? (2 / TimeFinal * $Time) : 1;
+ Else
+ // Fix parameters from the "Lib_Magnetodynamics2D_av_Cir.pro" template:
+ Flag_FrequencyDomain = 1;
+ Freq = frequency;
+ EndIf
+
mu0 = 4.e-7 * Pi;
nu[ Region[{Air, Ind, AirInf}] ] = 1. / mu0;
nu[ Core ] = 1. / (murCore * mu0);
@@ -55,10 +76,13 @@ Function {
Ns[ Ind ] = 1000 ; // number of turns in coil
Sc[ Ind ] = SurfaceArea[] ; // area of coil cross section
+
// 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];
- CoefGeos[] = 1;
+ js0[ Ind ] = Ns[] / Sc[] * Vector[0, 0, -1];
+
+ // For a correct definition of the voltage:
+ CoefGeos[] = 1; // planar model, 1 meter thick
}
Constraint {
@@ -70,8 +94,13 @@ Constraint {
}
{ Name Current_2D;
Case {
- // represents the phasor amplitude (peak to peak value) for a dynamic analysis
- { Region Ind; Value Current; }
+ If(Flag_FrequencyDomain)
+ // Amplitude of the phasor is set to "Current"
+ { Region Ind; Value Current; }
+ Else
+ // Time-dependent value is set to "Current * myModulation[]"
+ { Region Ind; Value Current; TimeFunction myModulation[]; }
+ EndIf
}
}
{ Name Voltage_2D;
@@ -87,8 +116,8 @@ PostOperation {
{ Name Map_a; NameOfPostProcessing Magnetodynamics2D_av;
Operation {
Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
- Print[ b, OnElementsOf Vol_Mag, File "b.pos" , HarmonicToTime 20];
- Print[ j, OnElementsOf Vol_Mag, File "j.pos", HarmonicToTime 20];
+ Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
+ Print[ j, OnElementsOf Vol_Mag, File "j.pos" ];
}
}
}
diff --git a/Magnetodynamics/transfo.pro b/Magnetodynamics/transfo.pro
index 026b91dc2a0f9ebf56becea298364419fa5f53ae..abeeed7b0be6f393b0cb342b7b428d85fae3d588 100644
--- a/Magnetodynamics/transfo.pro
+++ b/Magnetodynamics/transfo.pro
@@ -3,7 +3,6 @@
Features:
- Use of a generic template formulation library
- - Frequency- and time-domain dynamic solutions
- Circuit coupling used as a black-box (see Tutorial 8 for details)
To compute the solution in a terminal:
@@ -17,13 +16,12 @@
Include "transfo_common.pro";
DefineConstant[
- type_Conds = {2, Choices{1 = "Massive", 2 = "Coil"}, Highlight "Blue",
+ // The "Massive" option is non-physical, but it's interesting to highlight the
+ // fact that both massive and stranded conductors can be handled in the same
+ // way as far as circuit-coupling is concerned
+ ConductorType = {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,
+ Freq = {1, Min 0, Max 1e3, Step 1,
Name "Parameters/Frequency"}
];
@@ -43,147 +41,98 @@ Group {
// Abstract regions that will be used in the "Lib_Magnetodynamics2D_av_Cir.pro"
// template file included below;
Vol_Mag = Region[{Air, Core, Coils}]; // full magnetic domain
- If (type_Conds == 1)
+ If (ConductorType == 1)
Vol_C_Mag = Region[{Coils}]; // massive conductors
- ElseIf (type_Conds == 2)
+ ElseIf (ConductorType == 2)
Vol_S_Mag = Region[{Coils}]; // stranded conductors (coils)
EndIf
}
Function {
mu0 = 4e-7*Pi;
-
+ mur_Core = 1000;
mu[Air] = 1 * mu0;
-
- mur_Core = 100;
mu[Core] = mur_Core * mu0;
-
mu[Coils] = 1 * mu0;
- sigma[Coils] = 1e7;
-
- // For a correct definition of the voltage
- CoefGeo = thickness_Core;
+ nu[] = 1 / mu[];
- // 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)
+ sigma[Coils] = 1e7;
- Sc[Coil_1_M] = SurfaceArea[];
+ // To be defined separately for each coil portion, to fix the convention of
+ // positive current (1: along Oz, -1: along -Oz)
+ SignBranch[Coil_1_P] = 1;
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;
+ If(ConductorType == 2)
+ // Number of turns in the coils (same for PLUS and MINUS portions) - half
+ // values because half coils are defined geometrically due to symmetry!
+ Ns[Coil_1] = 100;
+ Ns[Coil_2] = 200;
+
+ // To be defined separately for each coil portion:
+ Sc[Coil_1_P] = SurfaceArea[];
+ Sc[Coil_1_M] = SurfaceArea[];
+ Sc[Coil_2_P] = SurfaceArea[];
+ Sc[Coil_2_M] = SurfaceArea[];
+
+ // Current density in each coil portion for a unit current (will be
+ // multiplied by the actual total current in the coil), in the case of
+ // stranded conductors
+ js0[Coils] = Ns[] / Sc[] * Vector[0, 0, SignBranch[]];
+ EndIf
- // The reluctivity will be used
- nu[] = 1/mu[];
+ // For a correct definition of the voltage
+ CoefGeos[Coils] = SignBranch[] * thickness_Core;
}
-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}];
- R_in = Region[10002]; // arbitrary region number (not linked to the mesh)
- Resistance_Cir += Region[{R_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;
+// We will use a circuit coupling to connect the PLUS and MINUS portions of the
+// coil in series, for both the primary and secondary. We will then apply a
+// voltage source (with a small resistance in series) on the resulting branch on
+// the primary, and connect a resistive load on the resulting branch on the
+// secondary.
+Flag_CircuitCoupling = 1;
- // End-winding primary winding resistance for more realistic primary coil
- // model
- Resistance[R_in] = 1;
- }
-
- 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}; }
- }
- }
+// Note that the voltage will not be equally distributed in the PLUS and MINUS
+// parts, which is the reason why we must apply the total voltage through a
+// circuit -- and we cannot simply use a current source like in Tutorial 7a.
- { Name ElectricalCircuit ; Type Network ;
- Case Circuit_1 {
- // PLUS and MINUS coil portions to be connected in series, together with
- // E_in (an additional resistor should be defined to represent the
- // Coil_1 end-winding (not considered in the 2D model))
- { Region E_in; Branch {1,4}; }
- { Region R_in; Branch {4,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 resistor 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} ; }
- }
- }
-
- }
+// Here is the definition of the circuits on primary and secondary sides:
+Group {
+ // Empty Groups to be filled
+ Resistance_Cir = Region[{}]; // all resistances
+ Inductance_Cir = Region[{}] ; // all inductances
+ Capacitance_Cir = Region[{}] ; // all capacitances
+ SourceV_Cir = Region[{}]; // all voltage sources
+ SourceI_Cir = Region[{}]; // all current sources
+
+ // Primary side
+ E_in = Region[10001]; // arbitrary region number (not linked to the mesh)
+ SourceV_Cir += Region[{E_in}];
+ R_in = Region[10002]; // arbitrary region number (not linked to the mesh)
+ Resistance_Cir += Region[{R_in}];
+
+ // Secondary side
+ R_out = Region[10101]; // arbitrary region number (not linked to the mesh)
+ Resistance_Cir += Region[{R_out}];
+}
-EndIf
+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;
+
+ // End-winding primary winding resistance for more realistic primary coil
+ // model
+ Resistance[R_in] = 1e-3;
+}
Constraint {
{ Name MagneticVectorPotential_2D;
@@ -193,17 +142,46 @@ Constraint {
}
{ Name Current_2D;
Case {
- If (type_Source == 1)
- // Current in each coil (same for PLUS and MINUS portions)
- { Region Coil_1; Value 1; TimeFunction F_Sin_wt_p[]{2*Pi*Freq, 0}; }
- { Region Coil_2; Value 0; }
- EndIf
}
}
{ Name Voltage_2D;
Case {
}
}
+ { Name Current_Cir ;
+ Case {
+ }
+ }
+ { Name Voltage_Cir ;
+ Case {
+ { Region E_in; Value val_E_in;
+ // F_Cos_wt_p[] is a built-in function with two parameters (w and p),
+ // that can be used to evaluate cos(w * t + p) in both frequency- and
+ // time-domain
+ 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 resistor is defined in series to represent the
+ // Coil_1 end-winding, which is not considered in the 2D model.
+ { Region E_in; Branch {1,4}; }
+ { Region R_in; Branch {4,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 resistor could be defined to represent the Coil_2
+ // end-winding - but we can directly add it to R_out as well)
+ { Region R_out; Branch {1,2}; }
+
+ { Region Coil_2_P; Branch {2,3} ; }
+ { Region Coil_2_M; Branch {3,1} ; }
+ }
+ }
}
Include "Lib_Magnetodynamics2D_av_Cir.pro";
@@ -214,40 +192,19 @@ PostOperation {
Print[ j, OnElementsOf Region[{Vol_C_Mag, Vol_S_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_Analysis == 1) // frequency domain
- 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!)
+ If (Flag_FrequencyDomain)
+ // 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)
+ // 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
- EndIf
}
}
}
diff --git a/Magnetostatics/electromagnet.pro b/Magnetostatics/electromagnet.pro
index 4261455468d6a29f9867713205cf8285d2d48b20..082881750ad4801906ddd20c4ace6aec5f6079cd 100644
--- a/Magnetostatics/electromagnet.pro
+++ b/Magnetostatics/electromagnet.pro
@@ -218,10 +218,10 @@ Integration {
/* The function space "Hregion_j_Mag_2D" provides one basis function, and hence
one degree of freedom, per physical region in the abstract region
- "Vol_S_Mag". The constraint "SourceCurrentDensityZ" fixes all these dofs, so
- the FunctionSpace "Hregion_j_Mag_2D" is fully fixed and has no FE unknowns.
- One could thus have replaced it by a simple function and the Integral term
- below in the weak formulation could have been written as
+ "Vol_S_Mag". The constraint "SourceCurrentDensityZ" fixes all these dofs (a
+ single one!), so the FunctionSpace "Hregion_j_Mag_2D" is fully fixed and has
+ no FE unknowns. One could thus have replaced it by a simple function and the
+ Integral term below in the weak formulation could have been written as
Integral { [ Vector[ 0,0,-js_fct[] ] , {a} ];
In Vol_S_Mag; Jacobian Vol; Integration Int; }