[Solved] Electrostatic: Charge spere in vacuum

Numerical methods and mathematical models of Elmer
manu1905
Posts: 5
Joined: 10 Oct 2016, 21:13
Antispam: Yes

[Solved] Electrostatic: Charge spere in vacuum

Post by manu1905 »

Dear all,

I am quite new with FEM and Elmer. I wish to perform some electrostatic calculations in 3D with some boundaries at a fixed potential. Within the « structure », I would like to include punctual charges.

Since I don’t know how to simulate punctual charges in FEM, I was thinking to simply have a sphere, the surface of which would have a particular surface charge density. I read quite a lot to know which boundary conditions should be applied at the surface of the sphere. I read that a surface with a given charge density can be seen as a surface with a given flux. I have also read the direct correspondance g = sigma with g being the flux and sigma the surface charge density (but I am not sure anymore if it was for Elmer or not). However, I could not precisely find what would be the units of sigma.

I then tried to a simple model: A sphere in a box. The sphere represents the charge I am willing to model, and I calculate the potential in the box (vacuum). You can find the mesh attached together with the SIF file.
The problem is that when I check the output of the simulation (with paraview), the potential doesn't show the behavior I would expect for a charged sphere (phi = q/(4*pi*epsilon0)…
Since I am new, I can possibly do many things wrong ... (mesh, boundary condition). Could you tell me what I should check?

Thanks a lot.
Attachments
squareWithHole4.msh
Gmsh mesh
(233.53 KiB) Downloaded 382 times
case.sif
(1.66 KiB) Downloaded 383 times
Slide cutting the sphere (shows potential)
Slide cutting the sphere (shows potential)
output.png (51.59 KiB) Viewed 6784 times
Last edited by manu1905 on 14 Oct 2016, 19:30, edited 1 time in total.
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Electrostatic: Charge spere in vacuum

Post by mzenker »

Hi,

you can post your sif file so someone can have a look and see if something is obviously wrong.
Also you can check your coordinates. Elmer assumes m by default, so if your geometry is e.g. in mm, you need to say

Code: Select all

Coordinate Scaling = 0.001
in the simulation section. If you use ElmerGUi, there is a free text field where you can put that command.

HTH,

Matthias
manu1905
Posts: 5
Joined: 10 Oct 2016, 21:13
Antispam: Yes

Re: Electrostatic: Charge spere in vacuum

Post by manu1905 »

Hi,

Thanks for your fast answer.
The coordinate scaling will help me a lot for the future. I am using ElmerGUI. That should be indicated in Model/Setup ?

Anyway, I think the scaling is not the problem here. Independently of the scale, I would expect the potential to decay ~1/r as I go away from the sphere?
I included the sif file as an attachment, but here you can read it directly:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals = 1
  Timestepping Method = BDF
  BDF Order = 1
  Solver Input File = case.sif
  Post File = case.ep
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

Body 1
  Target Bodies(1) = 2
  Name = "Body 1"
  Equation = 1
  Material = 1
End

Solver 2
  Equation = Result Output
  Output Format = Vtk
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = result
  Binary Output = False
  Scalar Field 1 = Potential
  Exec Solver = After Simulation
End

Solver 1
  Equation = Electrostatics
  Variable = Potential
  Procedure = "StatElecSolve" "StatElecSolver"
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Direct
  Linear System Direct Method = Banded
End

Equation 1
  Name = "Electrostatic"
  Active Solvers(2) = 2 1
End

Material 1
  Name = "Vacuum"
  Relative Permittivity = 1
End

Boundary Condition 1
  Target Boundaries(1) = 1 
  Name = "Charge"
  Electric Infinity BC = True
  Electric Flux = 1E-20
End
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Electrostatic: Charge spere in vacuum

Post by raback »

Hi

You need at least two boundary conditions. One for the sphere surface and the other for the farfield. Now there is just one BC.

Also, make sure that the sphere is really properly joined to the mesh of the cube. Try with some Dirichlet BCs at start to really see this.

-Peter
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Electrostatic: Charge spere in vacuum

Post by mzenker »

manu1905 wrote: I included the sif file as an attachment
sorry, I overlooked that...

Matthias
manu1905
Posts: 5
Joined: 10 Oct 2016, 21:13
Antispam: Yes

Re: Electrostatic: Charge spere in vacuum

Post by manu1905 »

Hi,

Thanks for the advise. Actually, imposing a fixed voltage (Dirichlet) condition on the sphere's surface also leads to not coherent results.
I think, you're right, the problem may come from the surface of the sphere incorrectly joined with the square.

I tried to modify my mesh without success. I copy the gmsh code I used to generate the mesh in the case somebody has an idea how I can fix it...

Code: Select all

Mesh.CharacteristicLengthFromCurvature = 0.5;

Point(1) = {0, 0, 0, 2.0};
Point(2) = {0, 20, 0, 2.0};
Line(1) = {2, 1};
Extrude {20, 0, 0} {
  Line{1};
}
out[] = Extrude {0, 0, 20} {
  Surface{5};
};



centerChargeX = 10.0;
centerChargeY = 8.0;
centerChargeZ = 10.0;
radiusCharge = 1.0;
lcCharge = 0.5;

Point(31) = {centerChargeX, centerChargeY, centerChargeZ, lcCharge};
Point(32) = {centerChargeX+radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(33) = {centerChargeX, centerChargeY+radiusCharge, centerChargeZ, lcCharge};
Circle(31) = {32,31,33};
Point(34) = {centerChargeX-radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(35) = {centerChargeX, centerChargeY-radiusCharge, centerChargeZ, lcCharge};
Circle(32) = {33,31,34};
Circle(33) = {34,31,35};
Circle(34) = {35,31,32};
Point(36) = {centerChargeX, centerChargeY, centerChargeZ-radiusCharge, lcCharge};
Point(37) = {centerChargeX, centerChargeY, centerChargeZ+radiusCharge, lcCharge};
Circle(35) = {33,31,36};
Circle(36) = {36,31,35};
Circle(37) = {35,31,37};
Circle(38) = {37,31,33};
Circle(39) = {32,31,37};
Circle(40) = {37,31,34};
Circle(41) = {34,31,36};
Circle(42) = {36,31,32};
Line Loop(33) = {32,38,-40};
Ruled Surface(34) = {33};
Line Loop(35) = {40,33,37};
Ruled Surface(36) = {35};
Line Loop(37) = {-38,-39,31};
Ruled Surface(38) = {37};
Line Loop(39) = {-41,-32,35};
Ruled Surface(40) = {39};
Line Loop(41) = {-35,-42,-31};
Ruled Surface(42) = {41};
Line Loop(43) = {-33,41,36};
Ruled Surface(44) = {43};
Line Loop(45) = {-37,34,39};
Ruled Surface(46) = {45};
Line Loop(47) = {-34,42,-36};
Ruled Surface(48) = {47};
Surface Loop(49) = {48,46,36,34,40,44,42,38};
Volume(50) = {49};

//Compound Volume(51) = {1,50};


Physical Surface("Charge") = {48,46,36,34,40,44,42,38};
//Physical Volume("Vacuum") = {1,50};
Physical Surface("Extern") = {22, 27, 14, 18, 5, 26};
Physical Volume("Vacuum") = {1};

If the mesh would be ok, should I apply a Electric Infinity BC to the outer part of the square? I didn't really understand to which boundary this condition should be applied (I tried the sphere and the box with my wrong mesh).
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Electrostatic: Charge spere in vacuum

Post by mzenker »

Hi,

I had a quick look at your mesh. Using the clipping functionality iun gmsh, I have the impression that your two bodies overlap each other. Your sphere needs to be cut out of the cube, and its outer surface included in the cube mesh. To do that, when you have built all surfaces, you have to include the sphere outer surface as a hole in the cube volume - the gmsh GUI asks you for holes automatically, then just select the sphere outer face. Then you declare the sphere as second volume and mesh everything. When saving, either declare everything you need in the final mesh as physical, or choose "save all (ignore physical entities)".

HTH,

Matthias
manu1905
Posts: 5
Joined: 10 Oct 2016, 21:13
Antispam: Yes

Re: Electrostatic: Charge spere in vacuum

Post by manu1905 »

Hi Matthias,

Thanks a lot for your indications.
However, I don't fully understand the procedure?

So, if I start from the surfaces, I have the following geo file:

Code: Select all

Mesh.CharacteristicLengthFromCurvature = 0.5;

Point(1) = {0, 0, 0, 2.0};
Point(2) = {0, 20, 0, 2.0};
Line(1) = {2, 1};
Extrude {20, 0, 0} {
  Line{1};
}

centerChargeX = 10.0;
centerChargeY = 8.0;
centerChargeZ = 10.0;
radiusCharge = 1.0;
lcCharge = 0.5;

Point(31) = {centerChargeX, centerChargeY, centerChargeZ, lcCharge};
Point(32) = {centerChargeX+radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(33) = {centerChargeX, centerChargeY+radiusCharge, centerChargeZ, lcCharge};
Circle(31) = {32,31,33};
Point(34) = {centerChargeX-radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(35) = {centerChargeX, centerChargeY-radiusCharge, centerChargeZ, lcCharge};
Circle(32) = {33,31,34};
Circle(33) = {34,31,35};
Circle(34) = {35,31,32};
Point(36) = {centerChargeX, centerChargeY, centerChargeZ-radiusCharge, lcCharge};
Point(37) = {centerChargeX, centerChargeY, centerChargeZ+radiusCharge, lcCharge};
Circle(35) = {33,31,36};
Circle(36) = {36,31,35};
Circle(37) = {35,31,37};
Circle(38) = {37,31,33};
Circle(39) = {32,31,37};
Circle(40) = {37,31,34};
Circle(41) = {34,31,36};
Circle(42) = {36,31,32};
Line Loop(33) = {32,38,-40};
Ruled Surface(34) = {33};
Line Loop(35) = {40,33,37};
Ruled Surface(36) = {35};
Line Loop(37) = {-38,-39,31};
Ruled Surface(38) = {37};
Line Loop(39) = {-41,-32,35};
Ruled Surface(40) = {39};
Line Loop(41) = {-35,-42,-31};
Ruled Surface(42) = {41};
Line Loop(43) = {-33,41,36};
Ruled Surface(44) = {43};
Line Loop(45) = {-37,34,39};
Ruled Surface(46) = {45};
Line Loop(47) = {-34,42,-36};
Ruled Surface(48) = {47};
Then I don't fully understand your comment:
To do that, when you have built all surfaces, you have to include the sphere outer surface as a hole in the cube volume - the gmsh GUI asks you for holes automatically, then just select the sphere outer face.
At which moment should gmsh ask me automatically about holes?
I tried to perform an extrusion of the squared surface, but gmsh GUI didn't ask me anything...
Could you describe the procedure to include the sphere's outer surface as a hole in the cube volume?

Sorry for what seems to be very basic questions...

Regards,
Manuel
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Electrostatic: Charge spere in vacuum

Post by mzenker »

Hi,

just use the GUI to define your volume, it is rather intuitive. It's under geometry->elementary entities->add->volume. Select the outer faces, press "e", select the hole surface(s), press "e", and it should be done.
You can do that in the geo file also, but I don't remember the syntax by heart. There is documentation on the gmsh website for the geofile language.

HTH,

Matthias
manu1905
Posts: 5
Joined: 10 Oct 2016, 21:13
Antispam: Yes

Re: Electrostatic: Charge spere in vacuum

Post by manu1905 »

Hi,

Thanks a lot. That solved the problem.
Actually, the volume of the cube is created when I performed the extrusion. I had to delete the volume (using the GUI) and to add a new one where I could specify the holes.

In attachment I show the potential over a slice of the cube passing through the center of the sphere. This is done twice, onces with a fixed potential at the sphere and then for a defined surface charge.

I post the code here in the case that somebody would have the same problem.

geo file

Code: Select all

// Gmsh project created on Sat Oct  8 14:04:08 2016

Mesh.CharacteristicLengthFromCurvature = 0.5;

Point(1) = {0, 0, 0, 2.0};
Point(2) = {0, 20, 0, 2.0};
Line(1) = {2, 1};
Extrude {20, 0, 0} {
  Line{1};
}
//out[] = Extrude {0, 0, 20} {
  //Surface{5};
//};



centerChargeX = 10.0;
centerChargeY = 8.0;
centerChargeZ = 10.0;
radiusCharge = 1.0;
lcCharge = 0.5;

Point(31) = {centerChargeX, centerChargeY, centerChargeZ, lcCharge};
Point(32) = {centerChargeX+radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(33) = {centerChargeX, centerChargeY+radiusCharge, centerChargeZ, lcCharge};
Circle(31) = {32,31,33};
Point(34) = {centerChargeX-radiusCharge, centerChargeY, centerChargeZ, lcCharge};
Point(35) = {centerChargeX, centerChargeY-radiusCharge, centerChargeZ, lcCharge};
Circle(32) = {33,31,34};
Circle(33) = {34,31,35};
Circle(34) = {35,31,32};
Point(36) = {centerChargeX, centerChargeY, centerChargeZ-radiusCharge, lcCharge};
Point(37) = {centerChargeX, centerChargeY, centerChargeZ+radiusCharge, lcCharge};
Circle(35) = {33,31,36};
Circle(36) = {36,31,35};
Circle(37) = {35,31,37};
Circle(38) = {37,31,33};
Circle(39) = {32,31,37};
Circle(40) = {37,31,34};
Circle(41) = {34,31,36};
Circle(42) = {36,31,32};
Line Loop(33) = {32,38,-40};
Ruled Surface(34) = {33};
Line Loop(35) = {40,33,37};
Ruled Surface(36) = {35};
Line Loop(37) = {-38,-39,31};
Ruled Surface(38) = {37};
Line Loop(39) = {-41,-32,35};
Ruled Surface(40) = {39};
Line Loop(41) = {-35,-42,-31};
Ruled Surface(42) = {41};
Line Loop(43) = {-33,41,36};
Ruled Surface(44) = {43};
Line Loop(45) = {-37,34,39};
Ruled Surface(46) = {45};
Line Loop(47) = {-34,42,-36};
Ruled Surface(48) = {47};
Coherence;
Extrude {0, 0, 20} {
  Surface{5};
}
Surface Loop(71) = {69, 5, 57, 61, 65, 70};
Volume(72) = {71};
Surface Loop(73) = {34, 40, 44, 36, 46, 48, 42, 38};
Volume(74) = {71, 73};
Volume(75) = {73};
Delete {
  Volume{1, 75};
}
Delete {
  Volume{72};
}
Delete {
  Volume{74};
}
Volume(74) = {71, 73};
Physical Volume("Vacuum") = {74};
Delete {
  Surface{48, 36, 40, 44, 34, 42, 38, 46};
}
Delete {
  Surface{69, 65, 70, 57, 61, 5};
}
Physical Surface("Extern") = {70, 65, 69, 57, 61, 5};
Physical Surface("Charge") = {34, 38, 42, 40, 48, 44, 46, 36};
Sif:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals = 1
  Timestepping Method = BDF
  BDF Order = 1
  Solver Input File = case.sif
  Post File = case.ep
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

Body 1
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1
End

Solver 2
  Equation = Result Output
  Output Format = Vtu
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = result
  Scalar Field 1 = Potential
  Exec Solver = After Simulation
End

Solver 1
  Equation = Electrostatics
  Variable = Potential
  Procedure = "StatElecSolve" "StatElecSolver"
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Direct
  Linear System Direct Method = Banded
End

Equation 1
  Name = "Electrostatics"
  Active Solvers(2) = 2 1
End

Material 1
  Name = "Vacuum"
  Relative Permittivity = 1
End

Material 2
  Name = "Material 2"
End

Boundary Condition 1
  Target Boundaries(1) = 3 
  Name = "FixedPotential"
  Electric Flux = 1.27e-20
End

Boundary Condition 2
  Target Boundaries(1) = 2 
  Name = "extern"
  Electric Infinity BC = True
End
Again, thanks a lot for the help.
Attachments
Sphere with a surface charge density of 1.27e-20 C/m2
Sphere with a surface charge density of 1.27e-20 C/m2
sliceConstantCharge.png (53.28 KiB) Viewed 6723 times
Sphere at fixed potential (1V)
Sphere at fixed potential (1V)
sliceConstantPotential.png (49.27 KiB) Viewed 6723 times
Post Reply