[Getdp] Electrokinetic - GetDP vs. Analytical Solution

mail at jchn.de mail at jchn.de
Thu Aug 16 11:15:07 CEST 2012


michael.asam at infineon.com <mailto:michael.asam at infineon.com> hat am 6. August
2012 um 11:55 geschrieben:

> Hi Jochen,
>
> sorry for the late response.
> The deviation is due to an inadequate mesh-size because the characteristic
> length varies linearly between the two electrodes, whereas the E field is
> proportional to r^-2.
> So the best is to use a field in Gmsh to specify the mesh size by
> an equation, or in case of more complex geometries do adaptive meshing.
>
> Please find attached a modification of your example for adaptive meshing.
> As your geometry is axisymmetrical I've reduced it to 2D.
> Just generate a mesh with Gmsh as usual and run GetDP. Do the post-
> processing "GenerateAdaptionFile" and plot it in Gmsh.
> There do a right click in the main window on the according view number
> and choose "Apply As Background Mesh" and then mesh the geometry
> once again.
> With the adapted mesh size the result should be better.
> For 3D the procedure should be similar.
>
> I hope this is of some help.
>
> Kind regards,
> Michael
>

Dear Michael,

thanks for your help!
Using your adaptive meshing example (2D), deviation from the analytical result
can be decreased from 0,05% with linear variation of mesh size to 0,006% with
applied adaptive meshing.
But when I transferred this procedure to my 3D problem, the deviation stayed at
approx. 5%. What I had to do, was decreasing the outer radius from 200 m to 50 m
and modeling just one quarter of the hemisphere to get more elements per volume
(see attached files). with approx. 600000 elements, I get a deviation of 1,03%,
which can be decreased by adaptive meshing to 0,94%.
Now my RAM is on its limit (32 Bit Win7). Are there any other ways to increase
the accuracy without increasing the element size?

Greetings,

Jochen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/getdp/attachments/20120816/93839e58/attachment.html>
-------------- next part --------------
// Gmsh project created on Wed Aug 15 09:11:36 2012

Earth_Electrode		=	1001	;
Reference_Ground	=	1002	;
Soil				=	1003	;

lc1	=	0.01;
lc2	=	1.1;

r1	=	1;
r2	=	50;

Point(1)	=	{0, 0, 0, lc1};
Point(2)	=	{-r1, 0, 0, lc1};
Point(3)	=	{0, -r1, 0, lc1};
Point(4)	=	{0, 0, -r1, lc1};
Point(5)	=	{-r2, 0, 0, lc2};
Point(6)	=	{0, -r2, 0, lc2};
Point(7)	=	{0, 0, -r2, lc2};

Line(1)	=	{2,5};
Line(2)	=	{3,6};
Line(3)	=	{4,7};

Circle(4)	=	{2,1,4};
Circle(5)	=	{5,1,7};
Circle(6)	=	{5,1,6};
Circle(7)	=	{7,1,6};
Circle(8)	=	{2,1,3};
Circle(9)	=	{4,1,3};

Line Loop(10) = {2, -6, -1, 8};
Ruled Surface(11) = {10};
Line Loop(12) = {1, 5, -3, -4};
Ruled Surface(13) = {12};
Line Loop(14) = {7, -2, -9, 3};
Ruled Surface(15) = {14};
Line Loop(16) = {7, -6, 5};
Ruled Surface(17) = {16};
Line Loop(18) = {8, -9, -4};
Ruled Surface(19) = {18};

Surface Loop(20) = {15, 17, 13, 11, 19};
Volume(21) = {20};

Physical Surface (Reference_Ground) = {17};
Physical Surface (Earth_Electrode) = {19};
Physical Volume (Soil) = {21};

-------------- next part --------------

Earth_Electrode		=	1001	;
Reference_Ground	=	1002	;
Soil				=	1003	;

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

Radius  =  50;  // Meters

Reference_Potential =	   0; //Volts
Earthing_Current = 	  -10000; //Ampere
Earthing_Current_Density = 	  0.25 * Earthing_Current / ( 0.25 * 2 * Pi * Radius^2 ) ; //Ampere / m^2

U_analytic[] 		= -Earthing_Current * 30 /  2 / Pi * ( 1/( Sqrt [ X[]^2 + Y[]^2 + Z[]^2 ] ) - 1 / 1 ) ;
E_analytic[] 		= -Earthing_Current * 30 /  2 / Pi * ( 1/( X[]^2 + Y[]^2 + Z[]^2 ) ) ;
}


// 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 ; }
    }
  }
}


// Integration
// ============
//
Integration {
  { Name GradGrad ;
    Case { {Type Gauss ;
            Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
                   { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
                   { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
                   { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
                   { GeoElement Prism       ; NumberOfPoints  9 ; } }
           }
         }
  }
}


// FunctionSpace
// ==============
// 
FunctionSpace {
  { Name Funspc_Name ; Type Form0 ;
    BasisFunction {
      { Name sn ; NameOfCoef vn ; Function BF_Node ;
        Support Elem_all ; Entity NodesOf[ All ] ;} 
	}
	Constraint {
      { NameOfCoef vn ; EntityType NodesOf ;
        NameOfConstraint Ref_Grnd ;}
    }
  }
}


// Formulation
// ============
// 
Formulation{
  { Name Electrokinetics; Type FemEquation;
    Quantity{
      { Name U; Type Local; NameOfSpace Funspc_Name;}
	}
    Equation {
      Galerkin { [ 1/rho[]*Dof{Grad U} , {Grad U} ]; In Int_Vol ;
                  Jacobian JVol ; Integration GradGrad ; }
	  Galerkin { [ -Earthing_Current_Density, {U} ]; In Reference_Ground ;
				  Jacobian JSur ; Integration GradGrad ; }
     }

  }
}


// Resolution
// ===========
//
Resolution{
  {Name EleKin_U;
    System{
      {Name Syst; NameOfFormulation Electrokinetics;}
    }
    Operation{
      Generate[Syst]; Solve[Syst]; SaveSolution[Syst];
    }
  }
}


// Post Processing
// ================
//
PostProcessing{
  {Name EleKin_U; NameOfFormulation Electrokinetics; NameOfSystem Syst;
    Quantity{
      {Name U;				Value {Local{[{U}];					In Elem_all;	Jacobian JVol;}}}
	  {Name Norm_U;			Value {Local{[Norm[{U}]];			In Elem_all;	Jacobian JVol;}}}
	  {Name U_analytic;     Value {Local{[ U_analytic[] ];      In Elem_all; 	Jacobian JVol;}}}
	  {Name E;				Value {Local{[{d U}];				In Elem_all;	Jacobian JVol;}}}
	  {Name Norm_E;			Value {Local{[Norm[{d U}]];			In Elem_all;	Jacobian JVol;}}}
	  {Name E_analytic;     Value {Local{[ E_analytic[] ];      In Elem_all; 	Jacobian JVol;}}}
	  {Name Error_E;  		Value {Local{[ Norm[ Norm[{d U}] / E_analytic[] - 1.0] * 100 ];  In Elem_all; Jacobian JVol;}}}
	  
	  }
  }
}


// Post Operation
// ===============
// 
PostOperation{
  {Name Map_U; NameOfPostProcessing EleKin_U;
    Operation{
      Print[ U, OnElementsOf Int_Vol, File "Result_U.pos"];
      // Print[ U, OnLine { {0,0,0} {200,0,0} } {2000}, Format Table, File "Result_U_rad.txt" ];
	  Print[ U_analytic, OnElementsOf Int_Vol, File "Result_U_analytic.pos"];
	  }
  }
  
    {Name Map_Error; NameOfPostProcessing EleKin_U;
    Operation{
      Print[ Error_E, OnElementsOf Int_Vol, File "Result_Error_E.pos"];
	  }
  }
  
  { Name GenerateAdaptationFile; NameOfPostProcessing EleKin_U; 
    Operation {	  
	  Print[ Norm_E,OnElementsOf Int_Vol, File "NewMeshSizeField.pos", Adapt H2, Dimension 3, Target 500000 ] ;
	  // Print[ Norm_E, OnElementsOf Int_Vol, File "Result_Norm_U.pos"];
    }
  }
}