[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>