[Getdp] Magnetodynamic 3D - Electromagnetic circuit element

Florin CIUPRINA florin at elmat.pub.ro
Wed Mar 26 13:43:02 CET 2003


Hello,

I am still modeling the L-shape conductor (3D) but now in the magnetodynamic
case (time harmonic and transient). I tried to use the a-v formulation with 
edge elements for "a"
and nodal elements with floating scalar potential on the terminals for "v".
The boundary conditions I tried to impose are a = 0 on the whole boundary, 
v = 1 on one terminal and
v = 0 on the other one. Since edge elements are used for a, I did not 
include the penalty term in the equation.
I obtained the solution, BUT  with a = 0 in all the domain which is not 
right (I think). The solution is quite the one obtained
in the electrical conduction formulation. I think that something is wrong 
in how  I tell GetDP the problem.
I appreciate very much if you would help me to solve this problem. (see 
problem files attached)

Thank you in advance!
Florin

Florin Ciuprina, Ph.D.
Associate Professor
----------------------------------------------------------------------------
POLITEHNICA University of Bucharest
Electrical Engineering Faculty
Electrical Materials Laboratory
Spl.Independentei 313, 77206, Bucharest, Romania
Tel:  +40.21.402 92 91
Fax: +40.21.410 43 55
e-mail: florin at elmat.pub.ro
----------------------------------------------------------------------------  
-------------- next part --------------
/* -------------------------------------------------------------------
   File "L3Dind.pro"

   This file defines the problem dependent data structures for
   an Magnetodynamic 3D problem.


   ------------------------------------------------------------------- */

Group {
   
  /* Let's start by defining the interface (i.e. Elementary groups)
     between the mesh file and GetDP (no mesh object is defined, so 
     the default mesh will be assumed to be in GMSH format and located 
     in "L3D.msh") */

  Conductor = Region[100];
  Ground = Region[111]; Hot = Region[110]; 
  Surf_Lat = Region[112];

  /* We can then define a global group (used in "MagDyn_av.pro",
     the file containing the function spaces and formulations) */

   DomainC_Mag = Region[{Conductor}]; 
   Surf = Region[{Surf_Lat, Ground, Hot}];
   Electrode_Mag = Region[{Hot}] ; // only used with global quantities


}

Function {

  /* The relative permeability (needed in the formulation) is piecewise
     defined in Elementary groups */

   sigma[Conductor] = 5.e7;
   mur[Conductor] = 1.;
   mu0 = 4.e-7 * Pi;
   nu[] = 1. / (mu0 * mur[]);
   Freq1 = 1.e1;
   Freq2 = 1.e2;
   Freq3 = 1.e3;
   Freq4 = 1.e4;
   Freq5 = 1.e5;

   Tc = 1.; 
   Mag_Time0 = 0.; Mag_TimeMax = 1.*Tc; Mag_DTime[] = Tc/20.; 
   Mag_Theta[] = 1.; 


}

Constraint {

/*
A 1-form function space containing vector potentials can be associated with a gauge condition, 
which can be defined as a constraint, e.g. a zero value is fixed for all circulations along edges 
of a tree (EdgesOfTreeIn) built in the mesh (Domain), having to be complete on certain 
boundaries (StartingOn Surf)
*/

{ Name GaugeCondition_a_Mag_3D; Type Assign;
    Case {
      { Region DomainC_Mag ; SubRegion Surf; Value 0.; }
    }
}

/* UNUSED
{ Name MagneticVectorPotential_3D; 
    Case {
      { Region Surf; Value 0.; }
    }
}
*/


  /* Now, some Dirichlet conditions are defined. The name 
     'ElectricScalarPotential' refers to the constraint name given in
     the function space */

  { Name ElectricScalarPotential; Type Assign;
    Case {
      { Region Ground; Value 0.; }
    }
  }
  { Name GlobalElectricPotential ;
    Case {
      { Region Hot ; Value 1.; }
    }
  }

}

/* The formulation used and its tools, considered as being 
   in a black box, can now be included */

Include "Jacobian_Lib.pro"
Include "Integration_Lib.pro"
Include "MagDyn_av.pro"

/* Finally, we can define some operations to output results */

PostOperation Map_vf_v UsingPost MagDyn_avf {
  Print[ v, OnElementsOf DomainC_Mag, File "Map_v.pos" ] ;
}

PostOperation Map_vf_j UsingPost MagDyn_avf {
  Print[ j, OnElementsOf DomainC_Mag, File "Map_j.pos" ] ;
}

PostOperation Map_vf_h UsingPost MagDyn_avf {
  Print[ h, OnElementsOf DomainC_Mag, File "Map_h.pos" ] ;
}

PostOperation Map_vf_a UsingPost MagDyn_avf {
  Print[ a, OnElementsOf DomainC_Mag, File "Map_a.pos" ] ;
}
PostOperation Map_vf_jg UsingPost MagDyn_avf {
  Print [ j, OnPlane {{0,1.5,0}{1,1.5,0}{0,1.5,1}} {10,10}, File "Cut_j" ];
}

PostOperation Val_V UsingPost MagDyn_avf {
  Print[ V, OnRegion Electrode_Mag, Format Table, File "VIZ.txt" ] ;
  Print[ I, OnRegion Electrode_Mag, Format Table, File >> "VIZ.txt"  ] ;
  Print[ Z, OnRegion Electrode_Mag, Format Table, File >> "VIZ.txt"   ] ;
}

// For test: post-calculation of a global quantity
// - can be used if the GlobalTerm with I is not included in the Formulation).
// - a function of type vf (in a SubSpace) must be defined.
// - OnRegion is for BasisFunction entity number.
// - the whole domain (DomainC_Mag) is the domain of integration.
PostOperation Ipos UsingPost MagDyn_avf {
  Print[ Ipos[DomainC_Mag], OnRegion Hot, Format Table, File "Ipos.txt" ] ; 
}
PostOperation Power UsingPost MagDyn_avf {
  Print [  roj2[Conductor], OnGlobal, Format Table, File "Power.txt" ];
}
 
-------------- next part --------------
/* -------------------------------------------------------------------
   File "L3D.geo"

   This file is the geometrical description used by GMSH to produce
   the file "L3D.msh".
   ------------------------------------------------------------------- */

/* Definition of some parameters for geometrical dimensions, i.e.
    w (width of 'Line')*/

w = 1.e-3 ;
h = 2.e-3 ; 

/* Definition of parameters for local mesh dimensions */

s = 1. ;
p0 = w / 3. * s ;
 

/* Definition of gemetrical points */

Point(1) = { 0, 0, 0, p0} ;
Point(2) = { w, 0, 0, p0} ;
Point(3) = { 2*w, 0, 0, p0} ;
Point(4) = { 2*w, w, 0, p0} ;
Point(5) = { w, w, 0, p0} ;
Point(6) = { 0, w, 0, p0} ;
Point(7) = { 0, 2*w, 0, p0} ;
Point(8) = { w, 2*w, 0, p0} ;

Point(21) = { 0, 0, h, p0} ;
Point(22) = { w, 0, h, p0} ;
Point(23) = { 2*w, 0, h, p0} ;
Point(24) = { 2*w, w, h, p0} ;
Point(25) = { w, w, h, p0} ;
Point(26) = { 0, w, h, p0} ;
Point(27) = { 0, 2*w, h, p0} ;
Point(28) = { w, 2*w, h, p0} ;



/* Definition of gemetrical lines */

/* L in the plane z = 0*/
Line(1) = {1,2};   Line(2) = {2,5};  Line(3) = {5,6};
Line(4) = {6,1};   Line(5) = {2,3};  Line(6) = {3,4}; 
Line(7) = {4,5};   Line(8) = {5,8};  Line(9) = {6,7};
Line(10) = {7,8}; 

/* L in the plane z = h*/
Line(21) = {21,22};   Line(22) = {22,25};  Line(23) = {25,26};
Line(24) = {26,21};   Line(25) = {22,23};  Line(26) = {23,24}; 
Line(27) = {24,25};   Line(28) = {25,28};  Line(29) = {26,27};
Line(30) = {27,28}; 

/* Lines between L-s*/ 
Line(41) = {1, 21};   Line(42) = {2, 22};  
Line(43) = {3, 23};   Line(44) = {4, 24}; 
Line(45) = {5, 25};   Line(46) = {6, 26};  
Line(47) = {7, 27};   Line(48) = {8, 28};  
 
/* Definition of geometrical surfaces */

/* Squares on L at z = 0 */
Line Loop(11) = {1,2,3,4};   Plane Surface(12) = {11};
Line Loop(13) = {5,6,7,-2};   Plane Surface(14) = {13};
Line Loop(15) = {-3,8,-10,-9};   Plane Surface(16) = {15};

/* Squares on L at z = h */
Line Loop(31) = {21,22,23,24};   Plane Surface(32) = {31};
Line Loop(33) = {25,26,27,-22};   Plane Surface(34) = {33};
Line Loop(35) = {-23,28,-30,-29};   Plane Surface(36) = {35};

/* Squares between L-s  */
Line Loop(51) = {1, 42, -21, -41}; Plane Surface(52) = {51};
Line Loop(53) = {5, 43, -25, -42}; Plane Surface(54) = {53};
Line Loop(55) = {-3, 45, 23, -46}; Plane Surface(56) = {55};
Line Loop(57) = {-7, 44, 27, -45}; Plane Surface(58) = {57};
Line Loop(59) = {10, 48, -30, -47}; Plane Surface(60) = {59};

Line Loop(61) = {-9, 46, 29, -47}; Plane Surface(62) = {61};
Line Loop(63) = {4, 41, -24, -46}; Plane Surface(64) = {63};
Line Loop(65) = {-8, 45, 28, -48}; Plane Surface(66) = {65};
Line Loop(67) = {-2, 42, 22, -45}; Plane Surface(68) = {67};
Line Loop(69) = {-6, 43, 26, -44}; Plane Surface(70) = {69};

/* Definition of geometrical volumes */

Surface Loop(101) = {12, 32, 52, 68, 56, 64}; Volume(102) = {101};
Surface Loop(103) = {14, 34, 68, 54, 70, 58}; Volume(104) = {103};
Surface Loop(105) = {16, 36, 56, 66, 60, 62}; Volume(106) = {105};

/* Definition of Physical entities (volumes, surfaces). The Physical
   entities tell GMSH the elements and their associated region numbers
   to save in the file 'L3D.msh'. */


Physical Surface (110) = {60} ;   /* Hot */
Physical Surface (111) = {70} ;   /* Ground */
Physical Surface (112) = {12,14,16,32,34,36,62,64,52,54,58,66} ;  /* Surf_Lat */
Physical Volume (100) = {102, 104, 106}; /* Mag + Conductor */


-------------- next part --------------
/* -------------------------------------------------------------------
   File "MagDyn_v.pro"

   Magnetodynamic - Magnetic vector potential and electric  scalar
   potential a-v formulation
   with floating electric potential
   ------------------------------------------------------------------- 

   I N P U T
   ---------

   Global Groups :  (Extension '_Mag' is for Magnetic problem)
   -------------
    Domain_Mag               Whole magnetic domain
    DomainC_Mag              Conducting regions
    Electrode_Mag            Electrodes (surfaces)
    Surf                     Lateral surfaces (Surf_Lat) + Electrodes

   Function :
   --------
   mur[]                   Relative permeability
   sigma[]                 Conductivity


   Constraint :
   ----------
    MagneticVectorPotential  Fixed magnetic vector potential    
                             (boundary condition)
    ElectricScalarPotential  Fixed electric scalar potential
                             (classical boundary condition)
    GlobalElectricPotential  Fixed global electric scalar potential
                             (for perfect electrodes)
    GlobalElectricCurrent    Fixed global electric current
                             (for perfect electrodes)


/* O U T P U T
   -----------

  PostQuantities :
  --------------
   v : Electric scalar potential
   a : Magnetic vector potential
   h : Magnetic field
   b : Magnetic flux density
   j : Electric current density
   V : Voltage
   I : Current
   Z : Impedance

*/
/* --------------------------------------------------------------------------*/

Group {
  DefineGroup[ Domain_Mag, DomainC_Mag, Electrode_Mag, Surf ] ;
}

Function {
  DefineFunction[ mur ];
  DefineFunction[ sigma ];
  DefineFunction[ nu ];
  DefineVariable[ Freq1, Freq2, Freq3, Freq4, Freq5 ];
  DefineVariable[Mag_Time0, Mag_TimeMax, Mag_DTime, Mag_Theta];
}

/* --------------------------------------------------------------------------*/

/* Edge finite element space with gauge condition */

FunctionSpace {
  { Name Hcurl_a_Gauge; Type Form1;
    BasisFunction {
     // a = a  s  
      //     e  e  
      { Name se; NameOfCoef ae; Function BF_Edge;
        Support DomainC_Mag ; Entity EdgesOfTreeIn[DomainC_Mag, StartingOn Surf]; }
    }
    Constraint {
     //gauge condition ae = 0 on tree edges and on some boundaries
      { NameOfCoef ae;
        EntityType EdgesOfTreeIn; EntitySubType StartingOn;
        NameOfConstraint GaugeCondition_a_Mag_3D; }
    }
  }


/* Edge finite element space */
/* UNUSED
FunctionSpace {
  { Name Hcurl_a_Gauge; Type Form1;
    BasisFunction {
     // a = a  s  
      //     e  e  
      { Name se; NameOfCoef ae; Function BF_Edge;
        Support DomainC_Mag ; Entity EdgesOf[All]; }
    }
    Constraint {
     //condition ae = 0 on some boundaries
      { NameOfCoef ae; EntityType EdgesOf; 
        NameOfConstraint MagneticVectorPotential_3D; }
    }
  }
*/

  { Name Hgrad_vf_Mag ; Type Form0 ;
    BasisFunction {
      // v = v  s  + v    s
      //      n  n    c,k  c,k
      { Name sn ; NameOfCoef vn ; Function BF_Node ;
        Support DomainC_Mag ; Entity NodesOf[ All, Not Electrode_Mag ] ; }
      { Name sck ;NameOfCoef vck ; Function BF_GroupOfNodes ;
        Support DomainC_Mag ; Entity GroupsOfNodesOf[ Electrode_Mag ] ; }
    }
    SubSpace { // only for a PostOperation
      { Name vf ; NameOfBasisFunction sck ; }
    }
    GlobalQuantity {
      { Name V ; Type AliasOf        ; NameOfCoef vck ; }
      { Name I ; Type AssociatedWith ; NameOfCoef vck ; }
    }
    Constraint {
      { NameOfCoef vn ;
        EntityType NodesOf ; NameOfConstraint ElectricScalarPotential ; }

      { NameOfCoef V ;
        EntityType GroupsOfNodesOf ; NameOfConstraint GlobalElectricPotential ; }
      { NameOfCoef I ;
        EntityType GroupsOfNodesOf ; NameOfConstraint GlobalElectricCurrent ; }
    }
  }
}


Formulation {
  { Name MagDyn_avf; Type FemEquation ;
    Quantity {
      { Name a ; Type Local  ; NameOfSpace Hcurl_a_Gauge ; }
      { Name v ; Type Local  ; NameOfSpace Hgrad_vf_Mag ; }
      { Name I ; Type Global ; NameOfSpace Hgrad_vf_Mag [I] ; }
      { Name V ; Type Global ; NameOfSpace Hgrad_vf_Mag [V] ; }
      { Name vf; Type Local  ; NameOfSpace Hgrad_vf_Mag [vf] ; }
    }
    Equation {
      Galerkin { [ nu[] * Dof{d a} , {d a} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }

      Galerkin { DtDof [ sigma[] * Dof{a} , {a} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      Galerkin { [ sigma[] * Dof{d v} , {a} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      
      Galerkin { DtDof [ sigma[] * Dof{a} , {d v} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      Galerkin { [ sigma[] * Dof{d v} , {d v} ]; In DomainC_Mag;
                 Jacobian Vol; Integration CurlCurl; }
      GlobalTerm { [ Dof{I} , {V} ]; In Electrode_Mag; }
    }
  }

}


Resolution {

/* Frequency resolution */
  { Name MagDyn_vf;
    System {
      { Name Sys_MagDyn; NameOfFormulation MagDyn_avf;
        Type ComplexValue; Frequency {Freq1, Freq2, Freq3, Freq4, Freq5}; }
    }
    Operation { 
      Generate[Sys_MagDyn]; Solve[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
    }
  }

/* Time domain resolution */
  { Name MagDyn_avf_Time;
    System {
      { Name Sys_MagDyn; NameOfFormulation MagDyn_avf; }
    }
    Operation {
      InitSolution[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
      TimeLoopTheta[Mag_Time0, Mag_TimeMax, Mag_DTime[], Mag_Theta[]] {
        Generate[Sys_MagDyn]; Solve[Sys_MagDyn]; SaveSolution[Sys_MagDyn]; 
      }
    }
  }
}



PostProcessing {
  { Name MagDyn_avf ; NameOfFormulation MagDyn_avf ;
    PostQuantity {
      { Name v ; Value { Term { [ {v} ]              ; In DomainC_Mag ; } } }
      { Name h ; Value { Term { [ nu[] * {d a} ]     ; In DomainC_Mag ; } } }
      { Name a ; Value { Term { [ {a} ]              ; In DomainC_Mag ; } } }
      { Name j ; Value { Term { [ - sigma[]*(Dt[{a}]+ {d v}) ]; In DomainC_Mag ; } } }
      { Name V ; Value { Term { [ {V} ]              ; In Electrode_Mag ; } } }
      { Name I ; Value { Term { [ {I} ]              ; In Electrode_Mag ; } } }
      { Name Z ; Value { Term { [ -{V}/{I} ]         ; In Electrode_Mag ; } } }

      { Name Ipos ;
        Value {
	  Integral {
	    [ - sigma[] * {d v} * BF{d vf} ] ;
	    In DomainC_Mag ; Jacobian Vol ; Integration CurlCurl ;
	  }
	}
      }

      { Name roj2;
        Value { 
          Integral { [ sigma[]*SquNorm[Dt[{a}]+{d v}] ]; In DomainC_Mag; 
                  Jacobian Vol; Integration CurlCurl;} 
        } 
      }



    }
  }
}


/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/