[Getdp] Electrokinetic - GetDP vs. Analytical Solution
mail at jchn.de
mail at jchn.de
Tue Jul 31 16:26:35 CEST 2012
Hello Everyone,
since a couple of days, I am comparing the solution of GetDP with the analytical
solution of a simple earthing arrangement.
There are two hemispheres, both have their origin in point (1)={0,0,0,lc1}, the
earth electrode (inner hemisphere) has a radius of r0 = 1 m, the reference
ground (outer hemisphere) has a radius of r1 = 200 m.
Between these hemispheres a current of -10 kA flows, the potential of the inner
hemisphere is set to 0 V (yes, really). The soil resistivity is set to rho = 30
Ohm*m.
The voltage at the outer hemisphere can be calculated as follows: U = I * rho /
2 / Pi * ( 1/r0 – 1/r1 ) = 47508 V
With the attached files, I get a result from GetDP of U = 44666 V
Is it possible that there’s something wrong with the meshing? I think the
deviation is too big to explain it with the different "solving procedure"...
Greetings,
Jochen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/getdp/attachments/20120731/deb7c1f3/attachment.html>
-------------- next part --------------
Include "num.geo";
lc1=0.05;
lc2=5;
Point(1) = {0, 0, 0, lc1};
Point(2) = {0, -1, 0, lc1};
Point(3) = {1, 0, 0, lc1};
Point(4) = {-1, 0, 0, lc1};
Point(5) = {0, 0, 1, lc1};
Point(6) = {0, 0, -1, lc1};
Point(7) = {0, -200, 0, lc2};
Point(8) = {200, 0, 0, lc2};
Point(9) = {-200, 0, 0, lc2};
Point(10) = {0, 0, 200, lc2};
Point(11) = {0, 0, -200, lc2};
Circle(1) = {4, 1, 5};
Circle(2) = {5, 1, 3};
Circle(3) = {3, 1, 6};
Circle(4) = {6, 1, 4};
Circle(5) = {4, 1, 2};
Circle(6) = {6, 1, 2};
Circle(7) = {3, 1, 2};
Circle(8) = {5, 1, 2};
Circle(9) = {9, 1, 10};
Circle(10) = {10, 1, 8};
Circle(11) = {8, 1, 11};
Circle(12) = {11, 1, 9};
Circle(13) = {8, 1, 7};
Circle(14) = {9, 1, 7};
Circle(15) = {10, 1, 7};
Circle(16) = {11, 1, 7};
Line Loop(19) = {6, -7, 3};
Ruled Surface(20) = {19};
Line Loop(21) = {7, 2, -8};
Ruled Surface(22) = {21};
Line Loop(23) = {8, 1, -5};
Ruled Surface(24) = {23};
Line Loop(25) = {5, -6, 4};
Ruled Surface(26) = {25};
Line Loop(29) = {16, -13, 11};
Ruled Surface(30) = {29};
Line Loop(31) = {13, -15, 10};
Ruled Surface(32) = {31};
Line Loop(33) = {15, -14, 9};
Ruled Surface(34) = {33};
Line Loop(35) = {14, -16, 12};
Ruled Surface(36) = {35};
Line(39) = {11, 6};
Line(40) = {4, 9};
Line(41) = {10, 5};
Line(42) = {3, 8};
Line Loop(43) = {11, 39, -3, 42};
Ruled Surface(44) = {43};
Line Loop(45) = {42, -10, 41, 2};
Ruled Surface(46) = {45};
Line Loop(47) = {41, -1, 40, 9};
Ruled Surface(48) = {47};
Line Loop(49) = {40, -12, 39, 4};
Ruled Surface(50) = {49};
Surface Loop(51) = {46, 44, 30, 36, 34, 32, 48, 24, 22, 20, 26, 50};
Volume(52) = {51};
Physical Surface (Reference_Ground) = {30,32,34,36};
Physical Surface (Earth_Electrode) = {20,22,24,26};
Physical Volume (Soil) = {52};
-------------- next part --------------
Include "num.geo";
// Group Data
// ===========
//
Group{
Soil = #Soil;
Earth_Electrode = #Earth_Electrode;
Reference_Ground = #Reference_Ground;
Int_Vol = Region[{Soil}];
Elem_all = Region[{Soil, Earth_Electrode, Reference_Ground}];
}
// Function
// ============
//
Function{
rho[Soil] = 30; //Ohm-Meters
Reference_Potential = 0; //Volts
Earthing_Current = -10000; //Ampere
}
// Constraint
// ==========
//
Constraint{
{Name Ref_Grnd; Type Assign;
Case{
{Region Earth_Electrode; Value Reference_Potential; }
}
}
}
// Jacobian
// =========
//
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
{ Name JLin ;
Case {
{ Region All ; Jacobian Lin ; }
}
}
}
Include "Integration_Lib.pro"; // Same as in the Reference Manual
// FunctionSpace
// ==============
//
FunctionSpace {
{ Name Funspc_Name ; Type Form0 ;
BasisFunction {
{ Name sn ; NameOfCoef vn ; Function BF_Node ;
Support Elem_all ; Entity NodesOf[ All, Not Reference_Ground ] ;}
{ Name sf; NameOfCoef vf; Function BF_GroupOfNodes;
Support Elem_all ; Entity GroupsOfNodesOf[Reference_Ground]; }
}
GlobalQuantity {
{ Name GlobalCurrent; Type AliasOf; NameOfCoef vf; }
}
Constraint {
{ NameOfCoef vn ; EntityType NodesOf ;
NameOfConstraint Ref_Grnd ;}
}
}
}
// Formulation
// ============
//
Formulation{
{ Name Electrokinetics; Type FemEquation;
Quantity{
{ Name v; Type Local; NameOfSpace Funspc_Name;}
{ Name V; Type Global; NameOfSpace Funspc_Name [GlobalCurrent]; }
}
Equation {
Galerkin { [ 1/rho[]*Dof{Grad v} , {Grad v} ]; In Int_Vol ;
Jacobian JVol ; Integration GradGrad ; }
GlobalTerm { [ -Earthing_Current, {V} ]; In Reference_Ground ; }
}
}
}
// Resolution
// ===========
//
Resolution{
{Name EleKin_v;
System{
{Name Syst; NameOfFormulation Electrokinetics;}
}
Operation{
Generate[Syst]; Solve[Syst]; SaveSolution[Syst];
}
}
}
// Post Processing
// ================
//
PostProcessing{
{Name EleKin_v; NameOfFormulation Electrokinetics; NameOfSystem Syst;
Quantity{
{Name v; Value {Local{[{v}];In Elem_all;Jacobian JVol;}}}
}
}
}
// Post Operation
// ===============
//
PostOperation{
{Name Map_v; NameOfPostProcessing EleKin_v;
Operation{
Print[v, OnElementsOf Elem_all, File "Potential.pos"];
Print [ v, OnLine { {0,0,0} {200,0,0} } {2000}, Format Table, File "Potential_rad.txt"