Gauss’ Law of Electrostatics#

Here we use the discretize package to solve for the electric potential (ϕ) and electric fields (e) in 2D that result from a static charge distribution. Starting with Gauss’ law and Faraday’s law:

E=ρϵ0×E=0E=ϕs.t.ϕ|Ω=0

where σ is the charge density and ϵ0 is the permittivity of free space. We will consider the case where there is both a positive and a negative charge of equal magnitude within our domain. Thus:

ρ=ρ0[δ(r+)δ(r)]

To solve this problem numerically, we use the weak formulation; that is, we take the inner product of each equation with an appropriate test function. Where ψ is a scalar test function and f is a vector test function:

Ωψ(E)dV=1ϵ0ΩψρdVΩfEdV=Ωf(ϕ)dV

In the case of Gauss’ law, we have a volume integral containing the Dirac delta function, thus:

Ωψ(E)dV=1ϵ0ψq

where q represents an integrated charge density. By applying the finite volume approach to this expression we obtain:

ψTMcDe=1ϵ0ψTq

where q denotes the total enclosed charge for each cell. Thus qi=ρ0 for the cell containing the positive charge and qi=ρ0 for the cell containing the negative charge. It is zero for every other cell.

ψ and q live at cell centers and e lives on cell faces. D is the discrete divergence operator. Mc is an inner product matrix for cell centered quantities.

For the second weak form equation, we make use of the divergence theorem as follows:

ΩfEdV=Ωf(ϕ)dV=1ϵ0Ω(fϕ)dV+1ϵ0Ω(f)ϕdV=1ϵ0Ωn(fϕ)da+1ϵ0Ω(f)ϕdV=0+1ϵ0Ω(f)ϕdV

where the surface integral is zero due to the boundary conditions we imposed. Evaluating this expression according to the finite volume approach we obtain:

fTMfe=fTDTMcϕ

where f lives on cell faces and Mf is the inner product matrix for quantities that live on cell faces. By canceling terms and combining the set of discrete equations we obtain:

[McDMf1DTMc]ϕ=1ϵ0q

from which we can solve for ϕ. The electric field can be obtained by computing:

e=Mf1DTMcϕ

Import Packages#

Here we import the packages required for this tutorial.

from discretize import TensorMesh
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import numpy as np
from discretize.utils import sdiag

Solving the Problem#

# Create a tensor mesh
h = np.ones(75)
mesh = TensorMesh([h, h], "CC")

# Create system
DIV = mesh.face_divergence  # Faces to cell centers divergence
Mf_inv = mesh.get_face_inner_product(invert_matrix=True)
Mc = sdiag(mesh.cell_volumes)
A = Mc * DIV * Mf_inv * DIV.T * Mc

# Define RHS (charge distributions at cell centers)
xycc = mesh.gridCC
kneg = (xycc[:, 0] == -10) & (xycc[:, 1] == 0)  # -ve charge distr. at (-10, 0)
kpos = (xycc[:, 0] == 10) & (xycc[:, 1] == 0)  # +ve charge distr. at (10, 0)

rho = np.zeros(mesh.nC)
rho[kneg] = -1
rho[kpos] = 1

# LU factorization and solve
phi = spsolve(A, rho)

# Compute electric fields
E = Mf_inv * DIV.T * Mc * phi

# Plotting
fig = plt.figure(figsize=(14, 4))

ax1 = fig.add_subplot(131)
mesh.plot_image(rho, v_type="CC", ax=ax1)
ax1.set_title("Charge Density")

ax2 = fig.add_subplot(132)
mesh.plot_image(phi, v_type="CC", ax=ax2)
ax2.set_title("Electric Potential")

ax3 = fig.add_subplot(133)
mesh.plot_image(
    E, ax=ax3, v_type="F", view="vec", stream_opts={"color": "w", "density": 1.0}
)
ax3.set_title("Electric Fields")
Charge Density, Electric Potential, Electric Fields
Text(0.5, 1.0, 'Electric Fields')

Total running time of the script: (0 minutes 1.010 seconds)

Gallery generated by Sphinx-Gallery