[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