[Getdp] Help integrating

Daniel R Morris dmorris at sparkchaser.com
Sun Jan 17 20:47:48 CET 2010


Hello!  What a wonderful tool GetDP is - if only I knew how to use it!

 

I am trying to solve an electrical conduction problem, a simplified version
of which I present here.  I have a conducting bar, with two contacts held at
fixed potential.  I have no trouble calculating scalar voltage, E field and
current density.  I am struggling, though, with integrating the current
density over a cross-section of the bar to measure total current.  I have
set up the bar 1cm dia and 10cm long, with half on either side of the X-Z
plane.  I can integrate 1 over the cross section at the X-Z plane and
correctly derive the cross-sectional area of the bar (after correcting with
2pi geometry factor for global quantity inAxiSymmetric problem).  When I try
to integrate current, I get 0 as a result - because E field isn't defined on
the line element I use as my cross-section cut-plane.  I am not sure how to
resolve this - either how to get E field (thus, current) defined on the
cut-plane or define the integration differently.  Any help is greatly
appreciated - thank you!

 

Daniel Morris

 

 

 

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

  Bar.geo

  Geometry of simple conducting bar

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

 

Include "BarPar.pro";

 

/* Conducting bar dimensions */

BarDia = 0.010;    // Bar diameter 1.0cm

BarLen = 0.100;    // Bar length   10.0cm

 

/* Mesh size parameters */

MeshScale = 1.0/10.0;    // General scaling parameter

lmEnd = (BarDia/2)*MeshScale;

lmMid = (BarDia/2)*MeshScale;

 

/* Names for physical objects */

CONDUCTOR = 1000;

ELECPOS   = 1001;

ELECNEG   = 1002;

CUTPLANE  = 1010;

 

 

Point(1) = {          0,  BarLen/2,  0,  lmEnd };

Point(2) = {   BarDia/2,  BarLen/2,  0,  lmEnd };

Point(3) = {   BarDia/2,         0,  0,  lmMid };

Point(4) = {   BarDia/2, -BarLen/2,  0,  lmEnd };

Point(5) = {          0, -BarLen/2,  0,  lmEnd };

Point(6) = {          0,         0,  0,  lmMid };

 

Line(101) = {1, 2};

Line(102) = {2, 3};

Line(103) = {3, 4};

Line(104) = {4, 5};

Line(105) = {5, 6};

Line(106) = {6, 1};

 

Line(110) = {6, 3};

 

Line Loop(200) = {101, 102, 103, 104, 105, 106};

Plane Surface(300) = {200};   // Conductor

 

 

/**** Physical Surfaces and Lines ****/

// Physical Surface used to represent the Conductor

Physical Surface(CONDUCTOR) = {300};

 

// Physical Line used to represent the Positive Electrode

Physical Line(ELECPOS) = {101};

 

// Physical Line used to represent the Negative Electrode

Physical Line(ELECNEG) = {104};

 

// Physical Line used to determine cross sectional area and total current

Physical Line(CUTPLANE) = {110};

 

 

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

  BarPar.pro 

  Parameters describing conducting bar

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

 

/* Conducting bar dimensions */

BarDia = 0.010;    // Bar diameter 1.0cm

BarLen = 0.100;    // Bar length   10.0cm

 

/* Conducting bar physical properties */

RhoBar = 0.70;    // Conductivity of material [Ohm-m]

 

/* Mesh size parameters */

MeshScale = 1.0/10.0;    // General scaling parameter

lmEnd = (BarDia/2)*MeshScale;

lmMid = (BarDia/2)*MeshScale;

 

/* Factor for integration in axisymmetric geometry */

CoeffGeo = 2.0*Pi;

 

/* Names for physical objects */

CONDUCTOR = 1000;

ELECPOS   = 1001;

ELECNEG   = 1002;

CUTPLANE  = 1010;

 

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

   Bar.pro

      This file solves DC Conduction in a simple round bar

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

 

Include "BarPar.pro";

 

Group {

   Conductor     = Region[CONDUCTOR];

   ElectrodePos  = Region[ELECPOS];

   ElectrodeNeg  = Region[ELECNEG];

   CutPlane      = Region[CUTPLANE];

   Conductors    = Region[{Conductor, ElectrodePos, ElectrodeNeg}];

   SigmaDefined  = Region[{Conductor,CutPlane}];

}

 

 

Function {

  Sigma[ SigmaDefined ]      = 1.0/RhoBar;   // Material conductivity

}

 

 

Constraint {

   { Name VoltScalarPotential; Type Assign;

      Case{

         { Region ElectrodePos;     Value 1.0; }

         { Region ElectrodeNeg;     Value 0.0; }

      }

   }

}

 

 

 

Jacobian {

  { Name JVol; Case { { Region All; Jacobian VolAxi; } } }

  { Name JSur; Case { { Region All; Jacobian SurAxi; } } }

}

 

 

Integration {

  { Name GaussInt;

    Case { { Type Gauss;

             Case { { GeoElement Triangle    ; NumberOfPoints 4; }

                    { GeoElement Line        ; NumberOfPoints 4; }

                    { GeoElement Quadrangle  ; NumberOfPoints 4; }

                    { GeoElement Tetrahedron ; NumberOfPoints 4; }

                    { GeoElement Hexahedron  ; NumberOfPoints 6; }

                    { GeoElement Prism       ; NumberOfPoints 9; } }

           }

         }

   }

}

 

 

FunctionSpace {

  { Name Vgrad; Type Form0;

    BasisFunction {

      { Name sn; NameOfCoef tn; Function BF_Node;

        Support Region[{Conductors}]; Entity NodesOf[ All ]; }

    }

    Constraint {

      { NameOfCoef tn; EntityType NodesOf; NameOfConstraint
VoltScalarPotential; }

    }

  }

}

 

 

Formulation {

  { Name Voltage_Steady ; Type FemEquation ;

    Quantity { 

      { Name Volt; Type Local; NameOfSpace Vgrad; }

    }

    Equation {

      Galerkin { [ Sigma[] *  Dof {d Volt},  {d Volt} ];  

                 In Conductors; Jacobian JVol; Integration GaussInt; }

    }

  }

}

 

 

Resolution {

  { Name Voltage_Steady ;

    System {

      { Name AA ; NameOfFormulation Voltage_Steady ; }

  }

    Operation { 

      Generate[AA]; Solve[AA]; SaveSolution[AA]; 

    }

  }

}

 

 

 

PostProcessing {

  { Name  VoltSteady; NameOfFormulation Voltage_Steady ;

    Quantity {

      { Name Volt;   Value { Local { [ {Volt} ]; In Conductor; Jacobian
JVol; } } }

      { Name eField; Value { Local { [ -{Grad Volt} ]; In Conductor;
Jacobian JVol; } } }

      { Name eNorm;  Value { Local { [Norm[ -{Grad Volt} ]]; In Conductor;
Jacobian JVol; } } }

      { Name jNorm;  Value { Local { [Norm[ -Sigma[] * {Grad Volt} ]]; In
Conductor; Jacobian JVol; } } }

 

/*--------------------------------  This integration correctly calculates
the cross-sectional area of the bar ---------------------------*/

      { Name surf;   Value { Integral { [1*CoeffGeo]; In CutPlane;
Integration GaussInt; Jacobian JSur; } } }

 

/*-------------------------------- This integration fails to correctly
calculate total current in the bar: eField (Grad Volt) not defined in
Cutplane so result is zero  -------------*/

      { Name cur;    Value { Integral { [-{Grad
Volt}*Vector[0,1,0]*CoeffGeo*Sigma[]]; In CutPlane; Integration GaussInt;
Jacobian JSur; } } }

    }

  }

}

 

 

 

PostOperation {

 

  { Name VoltSteady ; NameOfPostProcessing VoltSteady;

    Operation {

        Print[ Volt,   OnElementsOf Conductors, File "BarVolt.pos"];

        Print[ eField, OnElementsOf Conductors, File "BarEField.pos"];

        Print[ eNorm,  OnElementsOf Conductors, File "BarEFieldNorm.pos"];

        Print[ jNorm,  OnElementsOf Conductors, File "BarcurrDensity.pos"];

 

        Print[ surf[CutPlane], OnGlobal, File "surf.txt", Format Table];
// Cross-sectional area

        Print[ cur[CutPlane], OnGlobal, File "current.txt", Format Table];
// current

    }

  }

}

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/getdp/attachments/20100117/8d6e7d30/attachment.html>