Electrostatic problem

mathematics
1

Hi guys, I’m from Brazil and here no one know how to use FEniCS. I’m passing through a problem that probably you can help. I must simulate a little nanometric peace of a big metalic plane (on a potencial V) with some imperfections. Well, here is my idea: make a box that bottom surface is my irregular plane (V Dirichilet boundary condition), top surface is a simple plane (theoretically calculated V’ Dirichilet boundary condition) and side surfaces are simple planes (dV/dn = 0 vonn Neumann boundary condition). When I simulate on FEniCS I find values that no make sense (for example: my potencial energy muste be x 10^-10 but the simulation gives me y10^-32). I think could be a scale problem. In my geometry “one unity” = 100nm. I have already changed my constants, but it wasn’t work.

Could someone help me?

Here is my idea:

my_problem
my_problem2208×3243 192 KB

And here is my code:

from fenics import *
from dolfin import *
import meshio
import numpy


# MESH CONVERSION (.msh ---> .xdmf)

# read the 3D Gmsh file:
msh = meshio.read("fenics_test.msh")

# write the Gmsh file in xdmf format:
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))

# write surface data:
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

# write volume data:
meshio.write("cf.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}, cell_data={"tetra": {"name_to_read": msh.cell_data["tetra"]["gmsh:physical"]}}))

# create mesh:
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

# create boundary subdomain facet function:
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# create subdomain cell function:
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)



# SIMULATION

# Function space:
S_1 = FunctionSpace(mesh, 'P', 1)
VS_1 = VectorFunctionSpace(mesh, 'P', 1)

# Boundary conditions:
bc_bottom = DirichletBC(S_1, Constant(1000.0), mf, 1)   # 1: bottom boundary tag
bc_top = DirichletBC(S_1, Constant(904.5340), mf, 2)    # 2: top boundary tag

bcs = [bc_bottom, bc_top]

# Variational form:
dx = Measure('dx', domain=mesh, subdomain_data=cf)
phi = TrialFunction(S_1)
v = TestFunction(S_1)
f = Constant(0.0)
a = inner(grad(phi), grad(v))*dx
L = inner(f, v)*dx

# Solution:
phi = Function(S_1)                      # electric potential [V]
solve(a == L, phi, bcs)
e_0 = 8.854*10**(-12)                    # permittivity of free space
pi = 3.1416
a_O2 = 1.562*10**(-30)*4*pi*e_0          # O2 atomic polarizability
a_N2 = 1.710*10**(-30)*4*pi*e_0          # N2 atomic polarizability
E = Function(VS_1)
E.assign(project((-(grad(phi))), VS_1))        # electric field [V/m]
U = Function(S_1)
U.assign(project((-(a_O2*inner(E, E))), S_1))  # potential energy [J]
F = Function(VS_1)
F.assign(project((-(grad(U))), VS_1))          # force field [N]


# Check points
print(phi(0, 0, 5.001))
print(phi(0, 0, 50))
print(E(0, 0, 5.001))
print(U(0, 0, 5.001))
print(F(0, 0, 5.001))


# Plot:
XDMFFile('potencial.xdmf').write(phi)
XDMFFile('electric_field.xdmf').write(E)
XDMFFile('energy.xdmf').write(U)
XDMFFile('force_field.xdmf').write(F)

I really don’t know if I’m doing something wrong on my code or another thing. I think my geometry is ok, the potencial looks good on my domain.

3

You should not project the gradient of a CG-1 function into a CG-1 space, but a DG-0 space.

To avoid Machine precision errors, i would recommend not multiplying by a_O2, and rather scale your expected result by this constant.

1 Like