[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"];
}
}
}