[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-0001.htm>
```

More information about the getdp mailing list