diff --git a/CircuitCoupling/RLC_circuit.geo b/CircuitCoupling/RLC_circuit.geo new file mode 100644 index 0000000000000000000000000000000000000000..2bb92992bfbf2bbaba55bbda503144ed469bf6e9 --- /dev/null +++ b/CircuitCoupling/RLC_circuit.geo @@ -0,0 +1,33 @@ +// Just 3 cubes + +Include "RLC_circuit_common.pro"; + +lc = 0.25; // Characteristic length + +Point(1) = {0, 0, 0, lc}; +Extrude {1, 0, 0} { + Point{1}; +} +Extrude {0, 1, 0} { + Line{1}; +} +Extrude {0, 0, 1} { + Surface{5}; +} +Translate {2, 0, 0} { + Duplicata { Volume{1}; } +} +Translate {4, 0, 0} { + Duplicata { Volume{1}; } +} + +Physical Surface(Cube1Top) = {27}; +Physical Surface(Cube1Bottom) = {5}; +Physical Surface(Cube2Top) = {34}; +Physical Surface(Cube2Bottom) = {29}; +Physical Surface(Cube3Top) = {61}; +Physical Surface(Cube3Bottom) = {56}; + +Physical Volume(Cube1) = {1}; +Physical Volume(Cube2) = {28}; +Physical Volume(Cube3) = {55}; diff --git a/CircuitCoupling/RLC_circuit.jpg b/CircuitCoupling/RLC_circuit.jpg new file mode 100644 index 0000000000000000000000000000000000000000..d29f1b86a53fab8a4eee8ba2144030aad967e946 Binary files /dev/null and b/CircuitCoupling/RLC_circuit.jpg differ diff --git a/CircuitCoupling/RLC_circuit.pro b/CircuitCoupling/RLC_circuit.pro new file mode 100644 index 0000000000000000000000000000000000000000..ff8fae3a5c8e077c4eab22ec76155b56537ae8d9 --- /dev/null +++ b/CircuitCoupling/RLC_circuit.pro @@ -0,0 +1,365 @@ +/* ------------------------------------------------------------------- + Tutorial 8b : circuit coupling - RLC circuit + + Features: + - Electrokinetic formulation coupled with RLC circuit + - Transient resolution + + To compute the solution in a terminal: + getdp RLC_circuit -solve dynamic -pos Cube_Top_Values_ASCII + + To compute the solution interactively from the Gmsh GUI: + File > Open > RLC_circuit.pro + Run (button at the bottom of the left panel) + ------------------------------------------------------------------- */ + +/* + This example shows how to implement a circuit with discrete resistors, + capacitors, inductors, voltage- and current-sources in a transient finite + element simulation. There are three separate cubes, each with a certain + conductivity. The bottom of all cubes is kept at 0V via the constraint + Voltage_3D. The cubes are connected according the `RLC_circuit.jpg` schematic. +*/ + +Include "RLC_circuit_common.pro"; + +// FEM domain group data +Group { + // Surfaces + Top1 = Region[Cube1Top]; + Bottom1 = Region[Cube1Bottom]; + Top2 = Region[Cube2Top]; + Bottom2 = Region[Cube2Bottom]; + Top3 = Region[Cube3Top]; + Bottom3 = Region[Cube3Bottom]; + + // Volumes + Cube1 = Region[Cube1]; + Cube2 = Region[Cube2]; + Cube3 = Region[Cube3]; + + // FEM Electrical regions + // * all volumes + Vol_Ele = Region[{Cube1, Cube2, Cube3}]; + // * all surfaces connected to the lumped element circuit + Sur_Ele = Region[{Top1, Bottom1, Top2, Bottom2, Top3, Bottom3}]; + // * total electrical computation domain + VolWithSur_Ele = Region[{Vol_Ele, Sur_Ele}]; +} + +// Circuit group data - Fictive regions for the circuit elements +Group { + + // Sources + CurrentSource1 = Region[{3001}]; + VoltageSource1 = Region[{3002}]; + + // Resistors + R1 = Region[{4001}]; + R2 = Region[{4002}]; + + // Inductors + L1 = Region[{4003}]; + + // Capacitors + C1 = Region[{4004}]; + + Resistance_Cir = Region[{R1, R2}]; // All resistors + Inductance_Cir = Region[{L1}]; // All inductors + Capacitance_Cir = Region[{C1}]; // All capacitors + SourceI_Cir = Region[{CurrentSource1}]; // All current sources + SourceV_Cir = Region[{VoltageSource1}]; // All voltage sources + + // Complete circuit domain containing all circuit elements + Domain_Cir = Region[{Resistance_Cir, Inductance_Cir, Capacitance_Cir, + SourceI_Cir, SourceV_Cir}]; +} + +Function { + // Simulation parameters + tStop = 10.0e-6; // Stop time [s] + tStep = 100e-9; // Time step [s] + nStrobe= 1; // Strobe periode for saving the result + + SaveFct[] = !($TimeStep % nStrobe); // Only each nStrobe-th time-step is saved + + // FEM domain function data + // ------------------------ + + Vbottom = 0.0; // Absolute voltage at the bottom of all cubes + + // Geometry of a cube + w = 1.0; // Width [m] + l = 1.0; // Length [m] + h = 1.0; // Height [m] + + // Resistance + Rcube1 = 10.0; // Resistance of cube 1 [Ohm] + Rcube2 = 40.0; // Resistance of cube 2 [Ohm] + Rcube3 = 20.0; // Resistance of cube 3 [Ohm] + + // Specific electrical conductivity [S*m/m^2] + kappa[Cube1] = h / (w * l * Rcube1); + kappa[Cube2] = h / (w * l * Rcube2); + kappa[Cube3] = h / (w * l * Rcube3); + + // Circuit domain function data + // ---------------------------- + + // Resistors + R[R1] = 30.0; // Resistance of R1 [Ohm] + R[R2] = 40.0; // Resistance of R2 [Ohm] + + // Inductors + L[L1] = 100.0e-6; // Inductance of L1 [H] + + // Capacitors + C[C1] = 250.0e-9; // Capacitance of C1 [F] + + // Note: All voltages and currents are assigned / initialized in the + // Constraint-block +} + +// FEM domain constraint data +Constraint { + { Name Current_3D ; Type Assign ; + Case { + } + } + { Name Voltage_3D ; Type Assign ; + Case { + { Region Bottom1; Type Assign; Value Vbottom; } + { Region Bottom2; Type Assign; Value Vbottom; } + { Region Bottom3; Type Assign; Value Vbottom; } + } + } +} + +// Circuit domain constraint data +Constraint { + { Name Current_Cir ; + Case { + { Region CurrentSource1; Type Assign; Value 1.0; } // CurrentSource1 has 1.0 A + { Region L1; Type Init; Value 1.0; } // Initial current of L1 is 1.0 A + } + } + { Name Voltage_Cir ; + Case { + { Region VoltageSource1; Type Assign; Value 80.0; } // VoltageSource1 has 80.0 V + { Region C1; Type Init; Value 0.0; } // Initial voltage of C1 = 0.0 V + } + } + + { Name ElectricalCircuit ; Type Network ; + Case Circuit1 { + // Circuit for cube 1: + { Region VoltageSource1; Branch{ 1, 10};} + { Region R1; Branch{ 1, 2};} + // The voltage between nodes 2 and 3 corresponds to the absolute voltage + // of the surface Top1: + { Region Top1; Branch{ 2, 3};} + // Note the reverse order of the nodes in the next branch: it's for + // subtracting the bottom voltage. + { Region Bottom1; Branch{ 10, 3};} + // With the two lines above, the voltage between node 2 and 10 corresponds + // to the voltage drop over the cube regardless what the absolute voltages + // are. The current between node 2 and 3 corresponds to the absolute + // current flowing through Top1. The same is valid for Bottom1. Due to the + // reverse order here the current has the opposite sign. This is correct + // because the current flowing IN Top1 is flowing OUT of Bottom1. Note + // that in this particular circuit it is actually not necessary to include + // the bottom faces in the netlist, because all bottoms are kept at 0 V + // via constraint Voltage_3D... But then all Bottom faces have to be + // excluded also in the region Sur_Ele! + + // Circuit for cube 2: + { Region L1; Branch{ 4, 10};} + { Region R2; Branch{ 4, 10};} + { Region Top2; Branch{ 4, 5};} + { Region Bottom2; Branch{10, 5};} + + // Circuit for cube 3 + { Region CurrentSource1; Branch{ 6, 10};} + { Region C1; Branch{ 6, 10};} + { Region Top3; Branch{ 6, 7};} + { Region Bottom3; Branch{10, 7};} + } + } +} + +Jacobian { + { Name JVol ; + Case { + { Region Region[{Vol_Ele}]; Jacobian Vol; } + { Region Region[{Sur_Ele}]; Jacobian Sur; } + } + } +} + +Integration { + { Name I1 ; + Case { + { Type Gauss ; + Case { + { GeoElement Point ; NumberOfPoints 1 ; } + { GeoElement Line ; NumberOfPoints 5 ; } + { GeoElement Triangle ; NumberOfPoints 6 ; } + { GeoElement Quadrangle ; NumberOfPoints 7 ; } + { GeoElement Tetrahedron ; NumberOfPoints 15 ; } + { GeoElement Hexahedron ; NumberOfPoints 34 ; } + { GeoElement Prism ; NumberOfPoints 21 ; } + } + } + } + } +} + +FunctionSpace { + // Function space for the FEM domain + { Name Hgrad_V; Type Form0; + BasisFunction { + // All nodes but the grouped nodes of the surfaces for connecting the + // circuitry + { Name sn; NameOfCoef vn; Function BF_Node; + Support VolWithSur_Ele; Entity NodesOf[ All, Not Sur_Ele ]; } + // Grouped nodes: Each surface of Sur_Ele has just one single voltage + { Name sf; NameOfCoef vf; Function BF_GroupOfNodes; + Support VolWithSur_Ele; Entity GroupsOfNodesOf[ Sur_Ele ]; } + } + GlobalQuantity { + { Name U; Type AliasOf; NameOfCoef vf; } + { Name I; Type AssociatedWith; NameOfCoef vf; } + } + Constraint { + { NameOfCoef vn; EntityType NodesOf ; NameOfConstraint Voltage_3D; } + { NameOfCoef U; EntityType GroupsOfNodesOf ; NameOfConstraint Voltage_3D; } + { NameOfCoef I; EntityType GroupsOfNodesOf ; NameOfConstraint Current_3D; } + } + } + + // Function space for the circuit domain + { Name Hregion_Cir; Type Scalar; + BasisFunction { + { Name sn; 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 ; } + } + } + +} + + +Formulation { + { Name dynamic; Type FemEquation; + Quantity { + { Name V; Type Local; NameOfSpace Hgrad_V; } + { Name U; Type Global; NameOfSpace Hgrad_V [U]; } + { Name I; Type Global; NameOfSpace Hgrad_V [I]; } + { Name Uz; Type Global; NameOfSpace Hregion_Cir [Uz]; } + { Name Iz; Type Global; NameOfSpace Hregion_Cir [Iz]; } + } + Equation { + // FEM domain + Galerkin { [ kappa[] * Dof{d V}, {d V} ]; In Vol_Ele; + Integration I1; Jacobian JVol; } + GlobalTerm{ [ Dof{I}, {U} ]; In Sur_Ele; } + + // Circuit related terms + // Resistance equation + GlobalTerm{ [ Dof{Uz}, {Iz} ]; In Resistance_Cir; } + GlobalTerm{ [ R[] * Dof{Iz}, {Iz} ]; In Resistance_Cir; } + + // Inductance equation + GlobalTerm{ [ Dof{Uz}, {Iz} ]; In Inductance_Cir; } + GlobalTerm{ DtDof[ L[] * Dof{Iz}, {Iz} ]; In Inductance_Cir; } + + // Capacitance equation + GlobalTerm{ [ Dof{Iz}, {Iz} ]; In Capacitance_Cir; } + GlobalTerm{ DtDof[ C[] * Dof{Uz}, {Iz} ]; In Capacitance_Cir; } + + // Inserting the network + GlobalEquation { Type Network; NameOfConstraint ElectricalCircuit; + { Node {I}; Loop {U}; Equation {I}; In Sur_Ele; } + { Node {Iz}; Loop {Uz}; Equation {Uz}; In Domain_Cir; } + } + } + } +} + +Resolution { + { Name dynamic; + System { + { Name A; NameOfFormulation dynamic; } + } + Operation { + InitSolution[A]; + TimeLoopTheta[0.0, tStop, tStep, 1.0]{ + Generate[A]; + Solve[A]; + Test[SaveFct[]] { SaveSolution[A]; } + } + } + } +} + +PostProcessing { + { Name Ele; NameOfFormulation dynamic; + Quantity { + { Name V; Value{ Local{ [ {V} ]; In Vol_Ele; Jacobian JVol;} } } + { Name Jv; Value{ Local{ [ -kappa[]*{d V} ]; In Vol_Ele; Jacobian JVol;} } } + { Name J; Value{ Local{ [ Norm[kappa[]*{d V}] ]; In Vol_Ele; Jacobian JVol;} } } + { Name U; Value { Term { [ {U} ]; In Sur_Ele; } } } + // The minus sign here is for getting positive currents at the top of the + // cubes. This is because incomming currents are negative. + { Name I; Value { Term { [ -{I} ]; In Sur_Ele; } } } + } + } +} + +PostOperation { + // Absolute voltage everywhere [V] + { Name V; NameOfPostProcessing Ele; + Operation { + Print[ V, OnElementsOf Vol_Ele, TimeLegend, File "Result_V.pos"]; + } + } + // Current through the top surfaces of the cubes [A] + { Name I_Top; NameOfPostProcessing Ele; + Operation { + Print[ I, OnElementsOf Region[{Top1, Top2, Top3}], TimeLegend, + File "Result_I_Top.pos"]; + } + } + // Current density vectors everywhere [A/m^2] + { Name J_vectors; NameOfPostProcessing Ele; + Operation { + Print[ Jv, OnElementsOf Vol_Ele, TimeLegend, File "Result_J_vectors.pos"]; + } + } + // Magnitude of current density everywhere [A/m^2] + { Name J_magnitude; NameOfPostProcessing Ele; + Operation { + Print[ J, OnElementsOf Vol_Ele, TimeLegend, File "Result_J.pos"]; + } + } + // Store the results in an ASCII file + { Name Cube_Top_Values_ASCII; NameOfPostProcessing Ele; + Operation { + Print[U, OnRegion Top1, File "Result_Cube1TopValues.txt", Format TimeTable]; + Print[U, OnRegion Top2, File "Result_Cube2TopValues.txt", Format TimeTable]; + Print[U, OnRegion Top3, File "Result_Cube3TopValues.txt", Format TimeTable]; + + Print[I, OnRegion Top1, File > "Result_Cube1TopValues.txt", Format TimeTable]; + Print[I, OnRegion Top2, File > "Result_Cube2TopValues.txt", Format TimeTable]; + Print[I, OnRegion Top3, File > "Result_Cube3TopValues.txt", Format TimeTable]; + } + } + +} diff --git a/CircuitCoupling/RLC_circuit_common.pro b/CircuitCoupling/RLC_circuit_common.pro new file mode 100644 index 0000000000000000000000000000000000000000..65360ad102bfbb4f23c065b9da1180e09e9f5a44 --- /dev/null +++ b/CircuitCoupling/RLC_circuit_common.pro @@ -0,0 +1,12 @@ +// Physical numbers for volumes +Cube1 = 1001; +Cube2 = 1002; +Cube3 = 1003; + +// Physical numbers for surfaces +Cube1Top = 2001; +Cube1Bottom = 2002; +Cube2Top = 2003; +Cube2Bottom = 2004; +Cube3Top = 2005; +Cube3Bottom = 2006; diff --git a/CircuitCoupling/R_circuit.geo b/CircuitCoupling/R_circuit.geo index 8a5b6f29daa0dbfc9fcddfe18494f921cc38f210..77d983849190f948560f7007829dd1aacf9936a5 100644 --- a/CircuitCoupling/R_circuit.geo +++ b/CircuitCoupling/R_circuit.geo @@ -2,7 +2,7 @@ Include "R_circuit_common.pro"; -lc = 0.25; // Cjaracteristic length +lc = 0.25; // Cjaracteristic length Point(1) = {0, 0, 0, lc}; Extrude {1, 0, 0} { @@ -24,14 +24,14 @@ Translate {2, 4, 0} { Duplicata { Volume{1}; } } -Physical Surface(Cube1Top) = {27}; -Physical Surface(Cube1Bottom) = {5}; -Physical Surface(Cube2Top) = {34}; -Physical Surface(Cube2Bottom) = {29}; -Physical Surface(Cube3Top) = {61}; -Physical Surface(Cube3Bottom) = {56}; -Physical Surface(Cube4Top) = {88}; -Physical Surface(Cube4Bottom) = {83}; +Physical Surface(Cube1Top) = {27}; +Physical Surface(Cube1Bottom) = {5}; +Physical Surface(Cube2Top) = {34}; +Physical Surface(Cube2Bottom) = {29}; +Physical Surface(Cube3Top) = {61}; +Physical Surface(Cube3Bottom) = {56}; +Physical Surface(Cube4Top) = {88}; +Physical Surface(Cube4Bottom) = {83}; Physical Volume(Cube1) = {1}; Physical Volume(Cube2) = {28}; diff --git a/CircuitCoupling/R_circuit.pro b/CircuitCoupling/R_circuit.pro index 02afd17565dd48928a4f44917735dee2dadd8bfb..8ec804bd7b9774f4bc47b47cf88bcb1e1b41a0f0 100644 --- a/CircuitCoupling/R_circuit.pro +++ b/CircuitCoupling/R_circuit.pro @@ -7,19 +7,21 @@ - Implementation of a netlist To compute the solution in a terminal: - getdp R_circuit -solve Jstatic -pos Cube_Top_Values_as_ASCII_File + getdp R_circuit -solve static -pos Cube_Top_Values_ASCII To compute the solution interactively from the Gmsh GUI: File > Open > R_circuit.pro Run (button at the bottom of the left panel) ------------------------------------------------------------------- */ -// This example shows how to implement a circuit with discrete elements in a -// finite element simulation. Only lumped resistors, voltage- and -// current-sources are used. There are four separate cubes, each with a certain -// conductivity. The bottom of all cubes is kept at a constant voltage via the -// constraint Voltage_3D. The top of the cubes are connected with a circuit -// according the `R_circuit.jpg' schematic. +/* + This example shows how to implement a circuit with discrete elements in a + finite element simulation. Only lumped resistors, voltage- and current-sources + are used. There are four separate cubes, each with a certain conductivity. The + bottom of all cubes is kept at a constant voltage via the constraint + Voltage_3D. The top of the cubes are connected with a circuit according the + `R_circuit.jpg' schematic. +*/ Include "R_circuit_common.pro"; @@ -252,7 +254,7 @@ FunctionSpace { } Formulation { - { Name Jstatic; Type FemEquation; + { Name static; Type FemEquation; Quantity { { Name V; Type Local; NameOfSpace Hgrad_V; } { Name U; Type Global; NameOfSpace Hgrad_V [U]; } @@ -280,21 +282,20 @@ Formulation { } Resolution { - { Name Jstatic; + { Name static; System { - { Name Jstatic; NameOfFormulation Jstatic; } + { Name A; NameOfFormulation static; } } Operation { - InitSolution[Jstatic]; - Generate[Jstatic]; - Solve[Jstatic]; - SaveSolution[Jstatic]; + Generate[A]; + Solve[A]; + SaveSolution[A]; } } } PostProcessing { - { Name Ele; NameOfFormulation Jstatic; + { Name Ele; NameOfFormulation static; Quantity { { Name V; Value{ Local{ [ {V} ]; In Vol_Ele; Jacobian JVol;} } } { Name Jv; Value{ Local{ [ -kappa[]*{d V} ]; In Vol_Ele; Jacobian JVol;} } } @@ -325,7 +326,7 @@ PostOperation { } } // Store the results in an ASCII file - { Name Cube_Top_Values_as_ASCII_File; NameOfPostProcessing Ele; + { Name Cube_Top_Values_ASCII; NameOfPostProcessing Ele; Operation { Print[U, OnRegion Top1, File "Result_CubeTopValues.txt", Format SimpleTable]; Print[U, OnRegion Top2, File > "Result_CubeTopValues.txt", Format SimpleTable]; diff --git a/CircuitCoupling/R_circuit_common.pro b/CircuitCoupling/R_circuit_common.pro index bb6a5d3bac6def106a8ba8714293dd2ac8927171..3384474c3fb76d58ace0d025a9f233af22e40ba7 100644 --- a/CircuitCoupling/R_circuit_common.pro +++ b/CircuitCoupling/R_circuit_common.pro @@ -1,15 +1,15 @@ // Physical numbers for volumes -Cube1 = 1001; -Cube2 = 1002; -Cube3 = 1003; -Cube4 = 1004; +Cube1 = 1001; +Cube2 = 1002; +Cube3 = 1003; +Cube4 = 1004; // Physical numbers for surfaces -Cube1Top = 2001; -Cube1Bottom = 2002; -Cube2Top = 2003; -Cube2Bottom = 2004; -Cube3Top = 2005; -Cube3Bottom = 2006; -Cube4Top = 2007; -Cube4Bottom = 2008; \ No newline at end of file +Cube1Top = 2001; +Cube1Bottom = 2002; +Cube2Top = 2003; +Cube2Bottom = 2004; +Cube3Top = 2005; +Cube3Bottom = 2006; +Cube4Top = 2007; +Cube4Bottom = 2008; diff --git a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro index fa9c1dc018b310a7897bf1e0325373a8e73dbc1f..3a05f854ddfb08d274977ec821ddf39ab2d24ee6 100644 --- a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro +++ b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro @@ -39,10 +39,10 @@ Group { DefineGroup[ SourceV_Cir, // voltage sources SourceI_Cir, // current sources - Resistance_Cir, // resistances (linear) - Inductance_Cir, // inductances - Capacitance_Cir, // capacitances - Diode_Cir // diodes (nonlinear resistances) + Resistance_Cir, // resistors (linear) + Inductance_Cir, // inductors + Capacitance_Cir, // capacitors + Diode_Cir // diodes (treated as nonlinear resistors) ]; EndIf } diff --git a/Magnetodynamics/transfo.pro b/Magnetodynamics/transfo.pro index 27e6bb867398e02aa43a5b56338ceba6de1e0d02..b3b1ece3e5c659a0754241f4f0605244c3b24f19 100644 --- a/Magnetodynamics/transfo.pro +++ b/Magnetodynamics/transfo.pro @@ -171,7 +171,7 @@ ElseIf (type_Source == 2) // voltage { 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 + // 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,2}; } @@ -180,7 +180,7 @@ ElseIf (type_Source == 2) // voltage } 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 + // 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}; }