3D Mesh - multiple material [gmsh] [electrostatic] [FreeCAD]

mesh
1

Good day,

I am trying to do some electrostatic field calculations, considering a sphere within a sphere within a sphere case, (let’s just say that a copper ball is isolated with a thick PTFE layer, which is bounded by a far-field condition (air)).

The procedure which I am using to create a mesh:

  1. Firstly I am creating three individual models on FreeCAD, the copper surface (a small sphere), the PTFE area (a sphere without an area bounded by the copper surface, and the Far-field area (a sphere without an area bounded by the copper surface and PTFE area).
  2. Export the .step files.
  3. Import .step files to gmsh and merge them.
  4. Define surface numbering.
  5. Create mesh.

I would upload the mesh, however, the size of the .msh file is too big.

The code which I am using for my potential and electrical field calculations:

import fenics as fn
import numpy as np

epsilon_0 = 8.85e-12

Import the mesh, identify the subdomains and boundaries

and then make a function space and a dx element.

mesh = fn.Mesh(‘Test23.xml’)
subdomains = fn.MeshFunction(“size_t”, mesh, ‘Test23_physical_region.xml’)
boundaries = fn.MeshFunction(‘size_t’, mesh, ‘Test23_facet_region.xml’)
dx = fn.Measure(‘dx’, domain=mesh, subdomain_data=subdomains)
V0 = fn.FunctionSpace(mesh, ‘DG’, 0)
V = fn.FunctionSpace(mesh, ‘P’, 1)

class permittivity(fn.UserExpression):
def init(self, markers, **kwargs):
self.markers = markers
super().init(**kwargs)

def eval_cell(self, values, x, cell):
    if self.markers[cell.index] == 4:
        values[0] = 5
    else:
        values[0] = 1

eps = permittivity(subdomains, degree=1)

Set the boundary conditions using the boundary numbers

specified in Gmsh.

top = fn.DirichletBC(V, fn.Constant(1000), boundaries, 3)
bottom = fn.DirichletBC(V, fn.Constant(3000), boundaries, 1)
bcs =[top, bottom]

Solve the Poisson equation with the source set to 0

(the Laplace equation) via the L = fn.Constant(‘0’)… line

u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant(0) * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs)

Find the electric field by taking the negative of the

gradient of the electric potential. Project this onto

the function space for evaluation.

electric_field = fn.project(-fn.grad(u))

output the results to two separate files

potentialFile = fn.File(‘output6/potential.pvd’)
potentialFile << u
e_fieldfile = fn.File(‘output6/e_field.pvd’)
e_fieldfile << electric_field

The results of the calculations:

image
image932×721 44.7 KB

I have expected that the potential would linearly decrease up to the point of outer boundary conditions. The PTFE area has a defined 0V potential even though it is not applied in the python script.

My question is why the electrical potential is not distributed equally through the whole model? Is it the model representation, the mesh or the script?

Best wishes,
Ted

1 Like