A sphere in a point charge field.
In general, the LaplaceBIE allows finding the field on the body's surface if the potential in the absence of an object(s) is known. (Though one can not add objects one by one to get the solution directly. One, however, should be able to make an iterative scheme in the lines of multipole method.) To demonstrate that, let's consider a field on the dielectric sphere (radius r1) due to a point charge at distance ζ away. The solution inside the sphere is given as series [Straton page 221]:
where $L_n$ is Legendre polynomial and $\theta$ angle from line connecting point charge and spehre.
using Jacobi
function ψt(cosθ,ϵ,ζ,r1)
s = 0
for n in 0:25
s += (2*n + 1)/(n*ϵ + n + 1)*r1^n/ζ^(n+1)*legendre(cosθ,n)
end
return s
end
ψt (generic function with 1 method)
The normal derivative can also be given for the sphere which coincides with derivative in the radial direction:
function ∇ψn(cosθ,ϵ,ζ,r1)
s = 0
for n in 1:25
s += n*(2*n + 1)/(n*ϵ + n + 1)*r1^(n-1)/ζ^(n+1) * legendre(cosθ,n)
end
return s
end
∇ψn (generic function with 1 method)
First, we can define the exterior potential in the absence of objects:
using LinearAlgebra
ζ = 1.2
r1 = 1
y = [0,0,ζ]
ψ0(x) = 1/norm(x.-y)
∇ψ0(x) = -(x.-y)/norm(x.-y)^3
∇ψ0 (generic function with 1 method)
Let's now place a sphere in the field with relative permittivity $\epsilon=10$.
ϵ = 10
include("sphere.jl")
msh = unitsphere(2)
vertices, faces = msh.vertices, msh.faces
n = vertices
162-element Array{GeometryTypes.Point{3,Float64},1}:
[-0.525731, 0.850651, 0.0]
[0.525731, 0.850651, 0.0]
[-0.525731, -0.850651, 0.0]
[0.525731, -0.850651, 0.0]
[0.0, -0.525731, 0.850651]
[0.0, 0.525731, 0.850651]
[0.0, -0.525731, -0.850651]
[0.0, 0.525731, -0.850651]
[0.850651, 0.0, -0.525731]
[0.850651, 0.0, 0.525731]
⋮
[0.259892, 0.433889, -0.862668]
[0.262866, 0.16246, -0.951057]
[0.262866, -0.16246, -0.951057]
[0.961938, 0.0, 0.273267]
[0.951057, 0.262866, 0.16246]
[0.862668, 0.259892, -0.433889]
[0.951057, 0.262866, -0.16246]
[0.69378, 0.702046, 0.160622]
[0.850651, 0.525731, 0.0]
With LaplaceBIE to calculate the surface field, we execute these three lines:
using LaplaceBIE
ψ = surfacepotential(vertices,n,faces,ϵ,ψ0);
P∇ψ = tangentderivatives(vertices,n,faces,ψ);
n∇ψ = normalderivatives(vertices,n,faces,P∇ψ,ϵ,∇ψ0);
162-element Array{Float64,1}:
-0.06205534122187404
-0.06205534122187409
-0.062055341221873776
-0.0620553412218737
0.09154529388081165
0.09154529388081104
-0.06365584711404176
-0.06365584711404139
-0.06367958215679825
-0.04611694351733581
⋮
-0.06373122123157463
-0.06363392295695919
-0.06363392295695924
-0.058521276281101015
-0.06092044651491289
-0.06376068401697485
-0.06308197836685234
-0.06102685443985198
-0.061876924424271076
Now to compare with analytics, we use azimuthal symmetry, which in numerics was chosen as the x-axis:
using Winston
cosθ = range(-1,1,length=100)
cosθs = [x[3]/1 for x in vertices]
162-element Array{Float64,1}:
0.0
0.0
0.0
0.0
0.85065080835204
0.85065080835204
-0.85065080835204
-0.85065080835204
-0.5257311121191336
0.5257311121191336
⋮
-0.8626684804161862
-0.9510565162951536
-0.9510565162951536
0.2732665289126717
0.16245984811645317
-0.4338885645526948
-0.16245984811645317
0.16062203564002314
0.0
Potential on the surface vs analytics:
scatter(cosθs,ψ)
oplot(cosθ,ψt.(cosθ,ϵ,ζ,r1) )
savefig("potential.svg")
Normal derivatives on the surface vs analytics:
scatter(cosθs,n∇ψ)
oplot(cosθ,∇ψn.(cosθ,ϵ,ζ,r1) )
savefig("nderivatives.svg")