Hello, Netgen/NGSolve community,
I am trying to model the distribution of electric potential around a point source of DC electric current in radially symmetric cylindrical coordinates.
The problem I am trying to solve is given by the following equation:
[tex]\nabla \cdot (\sigma \nabla u) = - \nabla \cdot j[/tex]
with the source term
[tex]\nabla \cdot j = I \delta (x-x_s) \delta (y-y_s) \delta (z-z_s)[/tex]
where sigma is the distribution of conductivity in the domain, u is the electrostatic potential, j is the electric current density, I is a current injected at the point source located at (x[sub]s[/sub], y[sub]s[/sub], z[sub]s[/sub]) and delta is the Dirac’s delta.
In radially symmetric cylindrical coordinates (r, z), with the source located at (0, 0) this equation can be rewritten as:
[tex]\frac{1}{ r} \frac{ \partial}{ \partial r} [ \sigma \frac{\partial u}{\partial r}] + \frac{ \partial}{ \partial z} [ \sigma \frac{\partial u}{\partial z}] = \frac {I \delta (r) \delta (z)}{2 \Pi r}[/tex]
Originally I was using GMSH and FEniCS, but over time I realised that switching to Netgen/NGSolve might be beneficial since model geometry I am dealing with is easier to describe in Netgen than in GMSH. In addition, this will allow me to skip mesh conversion from GMSH format to FEniCS format which is surprisingly time-consuming.
Bellow, I supplied the code I managed to create. In the final implementation the domain around the current source will be much more complex but for the sake of simplicity of an example and to allow comparison with the analytical solution let’s assume that the current source is located in a homogenous medium. At the end of the code I’ve added a comparison with an analytical solution:
[tex]U (r_I) = \frac {I} {4 \Pi r_I \sigma } [/tex]
where r[sub]I[/sub] is a distance between the source point and the point at which the potential value is measured.
The numerical and analytical results should match (taking aside drop of potential caused by the finite size of the domain), but unfortunately, they are pretty far apart in case of my Netgen/NGSolve implementation while in FEniCS similar code gives me correct results. Since I have no experience with Netgen/NGSolve I probably messed something up. I guess that the problem is associated with the source term. In FEniCS I used build-in function PointSource() to add point source in one point of the domain. To my knowledge similar function is not present in Netgen/NGSolve, therefore I added point source using a coefficient function with Dirac’s delta approximation.
I would be very grateful if someone will help track the source of error. In addition, I wonder if it is possible to assign maxh value to the point or project it as a coefficient function on the entire mesh. In all examples I was looking at the maxh value was assigned to lines. In my implementation, I quite crudely added two short lines on both sides of the source location and assing maxh value to them.
from ngsolve import *
from netgen.geom2d import SplineGeometry
import numpy as np
import math
def MakeGeometry(d_r, d_z, maxh):
geometry = SplineGeometry()
# point coordinates
points = [ (0, 0), (0, 0.1), (0, d_z), (d_r, d_z), (d_r, -d_z), (0, -d_z), (0, -0.1)]
pnums = [geometry.AppendPoint(*p) for p in points]
# start-point, end-point, boundary-condition, domain on left side, domain on right side, maxh:
lines = [ (0,1,1,0,1,0.01), (1,2,1,0,1,maxh), (2,3,2,0,1,maxh), (3,4,2,0,1,maxh), (4,5,2,0,1,maxh), (5,6,1,0,1,maxh), (6,0,1,0,1,0.01)]
for p1,p2,bc,left,right, mesh_size in lines:
geometry.Append( ["line", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right, maxh=mesh_size)
return geometry
# Parameters
I = 1 #electric current
rho = 100 #resistivity
d_r = 100 #radius of domain
d_z = 100 # 1/2 of hight of domain
maxh = 10 # mesh size
alpha = 1/10 # parameter in Dirac's delta aproximation
sigma = CoefficientFunction([1/rho]) #conductivity
dd = CoefficientFunction(1/(abs(alpha)*math.sqrt(np.pi))*exp(-((x**2 + y**2)**(1/2)/alpha)**2)) # Dirac's delta
model_geometry = MakeGeometry(d_r, d_z, maxh)
mesh = Mesh(model_geometry.GenerateMesh(maxh=maxh))
fes = H1(mesh, order=2, dirichlet=[2])
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=False)
a += 2*np.pi*grad(u)*grad(v)*x*sigma*dx
f = LinearForm(fes)
f += dd*I*v*dx
a.Assemble()
f.Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f.vec
Draw (gfu)
#Comparison of numerical and analytical solution
import matplotlib.pyplot as plt
x_long = [0.01*i for i in range(d_r*100)]
x_short = [0.01*i for i in range(d_r*10)]
u_numerical_long = [gfu(mesh(p, 0.0)) for p in x_long]
u_numerical_short = [gfu(mesh(p, 0.0)) for p in x_short]
fig, axs = plt.subplots(2, figsize=(20, 10))
fig.suptitle('Netgen/NGSolve')
axs[0].plot(x_long, u_numerical_long, linewidth=3, label='Numerical')
axs[0].plot(x_long, I/(1/rho*4*np.pi*np.array(x_long)), linewidth=3, linestyle='--', label='Analytical')
axs[0].set_yscale('log')
axs[0].set_ylim(0.0001, 1000)
axs[0].grid()
axs[0].legend(loc=0)
axs[0].set_xlabel("r")
axs[0].set_ylabel("U(r)")
axs[1].plot(x_short, u_numerical_short, linewidth=3, label='Numerical')
axs[1].plot(x_short, I/(1/rho*4*np.pi*np.array(x_short)), linewidth=3, linestyle='--', label='Analytical')
axs[1].set_yscale('log')
axs[0].set_ylim(0.0001, 1000)
axs[1].grid()
axs[1].legend(loc=0)
axs[1].set_xlabel("r")
axs[1].set_ylabel("U(r)")
plt.show()EDIT: Attachments are in zip format because I was unable to add attachments as individual images.
Attachment: attachments.zip