[Getdp] given current => calculate the current density
Peter Kaufmann
mail at peter.in-berlin.de
Fri Apr 18 17:42:35 CEST 2008
Hello!
My problem sounds simple (and I think it have to be simple), but I did
not manage to solve it.
I want to calculate the current density in a 3D Volume. The current on
two surfaces are given (incoming and outgoing), the other surfaces are
isolated:
______
/_____/|
incoming | | | outgoing
=> | | | =>
current |____|/ current
My first attempt (testA.pro) was to set two Dirichlet Constraints
(incoming and outgoing current), use the equation of continuity and
calculate the current density direct via
Galerkin { [ Dof{j} , {d j} ] ;
In v_conductor ; Jacobian JVol ; Integration I1 ; }
... but I only get a current density on the two surfaces.
My second attempt was to calculate first the scalar potential and then
via PostProcessing the current density (testB.pro). The problem here is -
as far as I understand the problem - the Neumann Constraint on the
surfaces. My thoughts about this:
Generally
j[] = - sigma[] * {Grad{ phi }
holds. So on the surfaces I could write the Neumann Constraint using the
following Galerkin equation:
Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ;
In s_current_in ; Jacobian JSur ; Integration I1 ; }
... but running getdp I only get the error message "Wrong Initial
Constraints: remaining Dof(s) with non-fixed initial conditions".
Maybe some one can give me a hint or an example?
Is there a way to calculate the size of a surface (for example integrate
1 over the surface) in the "Function" section and use it in the
Constraint conditions?
Kind regards,
Peter Kaufmann
-------------- next part --------------
/*
test case for the calculation of a current density
works on a gold cube (1mm x 1 mm x 1mm)
*/
Group {
v_conductor = Region[38];
s_current_in = Region[40];
s_current_out = Region[39];
c_all = Region[{v_conductor, s_current_in, s_current_out}];
}
Function {
I = 1; // current in Ampere
f = 1; // surface in mm??
j0[ c_all ] = I/f; // incoming current density
resistivity_au = 2.44e-6; // electrical resistivity in Ohm*mm
conductance_au = 1/resistivity_au; // electrical conductance in 1/(Ohm*mm)
sigma [ c_all ] = conductance_au; // electrical conductance
}
Constraint {
// set the incoming current density
{ Name jIn ; Type Assign;
Case{
//{ Region s_current_in ; Value Vector[j0[], 0, 0] ; }
{ Region s_current_in ; Value Vector[j0[]/2, 0, 0] ; }
{ Region s_current_out ; Value Vector[j0[], 0, 0] ; }
}
}
}
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 6 ; }
}
}
}
}
}
FunctionSpace {
{ Name Jspace ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef tn ; Function BF_Node ;
Support Region[ v_conductor ] ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef tn ; EntityType NodesOf ; NameOfConstraint jIn; } // Dirichlet constraint
}
}
}
Formulation {
{ Name StaticJ ; Type FemEquation ;
Quantity {
{ Name j ; Type Local ; NameOfSpace Jspace ; }
}
Equation {
Galerkin { [ Dof{j} , {d j} ] ; In v_conductor ; Jacobian JVol ; Integration I1 ; }
}
}
}
Resolution {
{ Name StaticJSolution ;
System {
{ Name StaticJEquations ; NameOfFormulation StaticJ ; }
}
Operation {
Generate[StaticJEquations] ; Solve[StaticJEquations] ; SaveSolution[StaticJEquations] ;
}
}
}
PostProcessing {
{ Name StaticJPP; NameOfFormulation StaticJ ;
Quantity {
{ Name j ; Value { Local { [ {j} ]; In v_conductor ; Jacobian JVol ; } } }
}
}
}
PostOperation {
{ Name coils ; NameOfPostProcessing StaticJPP;
Operation {
Print[ j, OnElementsOf v_conductor, File "j.pos"] ;
}
}
}
-------------- next part --------------
/*
test case for the calculation of a current density
*/
Group {
v_conductor = Region[38];
s_current_in = Region[40];
s_current_out = Region[39];
c_all = Region[{v_conductor, s_current_in, s_current_out}];
}
Function {
I = 1; // current in Ampere
f_in = 2; // surface in mm??
f_out = 1; // surface in mm??
j[ s_current_in ] = Vector[I/f_in, 0, 0]; // incoming current density
j[ s_current_out ] = Vector[I/f_out, 0, 0]; // outgoing current density
resistivity_au = 2.44e-6; // electrical resistivity in Ohm*mm
conductance_au = 1/resistivity_au; // electrical conductance in 1/(Ohm*mm)
sigma [ c_all ] = conductance_au; // electrical conductance
}
Constraint {
// set the electrical potential on the drain side to zero (only the
// potential difference matters!)
{ Name phiScalarPotential ; Type Assign;
Case{
{ Region s_current_out ; Value 0.0 ; }
}
}
}
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 3 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 6 ; }
}
}
}
}
}
FunctionSpace {
{ Name phi ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef tn ; Function BF_Node ;
Support Region[ v_conductor ] ; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef tn ; EntityType NodesOf ; NameOfConstraint phiScalarPotential; } // Dirichlet constraint
}
}
}
Formulation {
{ Name StaticPotential ; Type FemEquation ;
Quantity {
{ Name phi ; Type Local ; NameOfSpace phi ; }
}
Equation {
Galerkin { [ sigma[] * Dof {Grad phi} , {Grad phi} ] ;
In v_conductor ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ; // Neumann Constraint
In s_current_in ; Jacobian JSur ; Integration I1 ; }
Galerkin { [ j[] + sigma[] * {Grad phi} , Dof {phi} ] ; // Neumann Constraint
In s_current_out ; Jacobian JSur ; Integration I1 ; }
}
}
}
Resolution {
{ Name StaticPotentialSolution ;
System {
{ Name StaticPotentialEquations ; NameOfFormulation StaticPotential ; }
}
Operation {
Generate[StaticPotentialEquations] ; Solve[StaticPotentialEquations] ; SaveSolution[StaticPotentialEquations] ;
}
}
}
PostProcessing {
{ Name StaticPotentialPP; NameOfFormulation StaticPotential ;
Quantity {
{ Name phi ; Value { Local { [ {phi} ]; In v_conductor ; Jacobian JVol ; } } }
{ Name j; Value { Local { [ -sigma[] * {Grad phi} ]; In v_conductor ; Jacobian JVol ; } } }
}
}
}
PostOperation {
{ Name coils ; NameOfPostProcessing StaticPotentialPP;
Operation {
Print[ phi, OnElementsOf v_conductor, File "phi.pos"] ;
Print[ j, OnElementsOf v_conductor, File "j.pos"] ;
}
}
}
-------------- next part --------------
a=1;
Point(1) = {0, -0.5, 0, 0.5*a};
Point(2) = {0, 1.5, 0, 0.5*a};
Point(3) = {1, 1, 0, 0.5*a};
Point(4) = {1, 0, 0, 0.5*a};
Line (1) = {1, 2};
Line (2) = {2, 3};
Line (3) = {3, 4};
Line (4) = {4, 1};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
Extrude {0, 0, 1} {
Surface{6};
}
Line Loop(29) = {3,22,-10,-18};
Plane Surface(30) = {29};
Line Loop(31) = {8,-14,-1,13};
Plane Surface(32) = {31};
Line Loop(33) = {9,10,11,8};
Plane Surface(34) = {33};
Plane Surface(35) = {5};
Surface Loop(36) = {28,23,6,19,27,15};
Volume(37) = {36};
Physical Volume(38) = {1}; // volume
Physical Surface(39) = {23}; // right surface
Physical Surface(40) = {15}; // left surface