Variational form of curl equation issue

variational formulation
1

Hello,

I am trying to evaluate the electric field in a 2D conductor generated by a known external time varying magnetic field. The electric field verifies the two equations :

curl(\boldsymbol{E}) = - {\partial \boldsymbol{B} \over \partial t}

div(\boldsymbol E) = 0

The 2D domain \Omega is surrounded by a closed contour \partial \Omega as its boundary where the electric field is tangential to the contour :

\boldsymbol{E} \cdot \boldsymbol{n} = 0 on \partial \Omega

The variational formulation with a vector test function \boldsymbol{F} for the first equation would be :

\int_{\Omega} curl(\boldsymbol{E}) \cdot \boldsymbol{F} \,dx = - \int_{\Omega} {\partial \boldsymbol{B} \over \partial t} \cdot \boldsymbol{F} \,dx.

  1. It is a first order equation, so no integration by parts is necessary to decrease the constraint on the finite elements. And also I donā€™t know how a curl fits in an integration by partsā€¦

  2. \Omega is a 2D domain so E_z = 0 and E_x and E_y are independant of z.
    So curl(\boldsymbol{E}) lies only in the z direction and the variational form will be not null only for F_z. I suspect this is the reason why naively assembling the A matrix doesnā€™t work in the example belowā€¦

Expressing the variational form in the z axis results in a scalar equation. That would require a scalar test function while \boldsymbol{E} is still a 2D vector. The test and trial function would then not be in the same space leading to a non square linear system to solve.

The boundary condition also is a scalar equationā€¦

What is the proper way to formalize this problem to solve it efficiently with Fenics ?

from fenics import *

mesh = UnitSquareMesh(8,8)

V = VectorFunctionSpace(mesh, 'P', 1, dim=3)

bcs=[]

E = TrialFunction(V)
F = TestFunction(V)
B_dot = Constant((0.0, 0.0, -1.0))

a = dot(curl(E),F)*dx
L = dot(B_dot,F)*dx

A = assemble(a)

# Returns : Index out of bounds.
2

I would suggest that you formulate the problem in terms of the magnetic vector potential \mathbf{A}.
\nabla\times\nabla\times\mathbf{A} = \mu_0\mathbf{J},
where, using the Coulomb Gauge \nabla\cdot\mathbf{A} = 0 and low
frequency operation is assumed.
Electric field is then given by
\mathbf{E}= -\frac{\partial\mathbf{A}}{\partial t} in the absence of charge sources.
and the magnetic field is \mathbf{B}=\nabla\times\mathbf{A}. This tends to be the standard way of handling these sorts of problems.

The current \mathbf{J} will contain terms proportional to the electric field (conduction currents), i.e. \mathbf{J} = \sigma\mathbf{E} = -\sigma\frac{\partial\mathbf{A}}{\partial t}.

The functional expression in this case would follow the usual curl-curl formalism found in any reference treating the finite element solution of electromagnetic problems. If you are solving the time-harmonic problem, the time derivative can be replaced with j\omega.

3

Thank you for your kind help : I wish I can find a ā€œstandard wayā€ of handling this issue !

In this case though, \boldsymbol{B} is an ā€œexternalā€ known field. The current densities or magnetisations generating this field are outside of \Omega.

Equation

\nabla \times \nabla \times \boldsymbol{A} = \mu_0 \boldsymbol{J}

or

\nabla \times \boldsymbol{B} = \mu_0 \boldsymbol{J}

will just result in \boldsymbol{J}=0 as the current sources generating the known magnetic field are precisely outside \Omega.

(Note that \boldsymbol{J} \neq 0 in this case would mean that a static external magnetic field does generate currents in a static conductorā€¦ :wink:)

So, ok, \boldsymbol{E} = {\partial \boldsymbol{A} \over \partial t} (that is Faradayā€™s law, the one inducing electric field (and currents in conductors) from time varying magnetic field), but to find \boldsymbol{A}, I need to solve the equation :

\nabla \times \boldsymbol{A} = \boldsymbol{B}

which is equivalent to the initial problem.

So I still canā€™t see how to link this problem to the curl-curl formalismā€¦

Is it really a big deal to solve \nabla \times \boldsymbol{u} = \boldsymbol{f} in 2D using finite elements ?

4

ā€œ> Is it really a big deal to solve āˆ‡Ć—u=f in 2D using finite elements ?ā€

You really want to base a formulation on some type of variational form, which exhibits useful stationarity properties. I suggest you have a look at a finite element reference to see why this is. Silvester and Ferrariā€™s ā€œFinite Elements for Electrical Engineersā€ is a good place to start.

In brief, you generally want to start with the functional
F(\mathbf{A}) = \frac{1}{2}\iint\limits_{\Omega}\left\lbrace\nabla\times\mathbf{A}\cdot\nabla\times\mathbf{A}+\sigma\mathbf{A}\cdot\frac{\partial\mathbf{A}}{\partial t}\right\rbrace\, dS + B(\mathbf{A},\mathbf{B}^\mathrm{inc})_{\partial\Omega},
where B() is a boundary operator over the perimeter \partial\Omega where you treat the applied field similar to how you would treat a scattering problem.

If you really want to treat \nabla\times\mathbf{u}=\mathbf{f} directly, you could always form the least-squares residual,
R=\vert\!\vert\nabla\times\mathbf{u-f}\vert\!\vert^2
and find the solutions that minimize R. Note that this changes how you treat the boundary conditions, though.

5

Thank you for your kind response. Unfortunately I will not have access to the Silvester and Ferrariā€™s book for some timeā€¦ I had a hard time trying to figure out how to solve the least-squares residual variational formulation and didnā€™t find any example in the fenics tutorials. But there is this publication from Rickard Bergstrƶm hidden somewhere in the fenicsproject website. And I tried this :

from fenics import *

mesh = UnitSquareMesh(16, 16)

n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, 'P', 2)

u = TrialFunction(V)
v = TestFunction(V)

def curl2D(vector) :
    return vector[1].dx(0) - vector[0].dx(1)

B_dot = Constant(1.0)
Zero = Constant(0.0)

a1 = inner(curl2D(u), curl2D(v))*dx
a2 = inner(div(u), div(v))*dx
a3 = inner(dot(u, n), dot(v, n))*ds

l1 = inner(-B_dot, curl2D(v))*dx
l2 = inner(Zero, div(v))*dx
l3 = inner(Zero, dot(v, n))*ds

a = a1 + a2 + a3
l = l1 + l2 + l3

u = Function(V)
solve(a == l, u)

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
plt.subplot(2, 1, 1)

plot(u)

print(f"u(0, 0.5) = {u(0, 0.5)}")
print(f"u(0.5, 0.5) = {u(0.5, 0.5)}")
print(f"u(1, 0.5) = {u(1, 0.5)}")
print(f"u(0.5, 0) = {u(0.5, 0)}")
print(f"u(0.5, 1) = {u(0.5, 1)}")

du_dx = project(u.dx(0), V)
du_dy = project(u.dx(1), V)

print(f"du_dx(0.5, 0.5) = {du_dx(0.5, 0.5)}")
print(f"du_dy(0.5, 0.5) = {du_dy(0.5, 0.5)}")

print(f"du_dx(0, 0.5) = {du_dx(0, 0.5)}")
print(f"du_dy(0, 0.5) = {du_dy(0, 0.5)}")


plot_range = np.arange(.0, 1, 1/32)
np.append(plot_range, [1.0])

plt.subplot(2, 1, 2)
plt.plot(plot_range, [u(x, 0.5)[0] for x in plot_range], '.-', label = "Ex")
plt.plot(plot_range, [u(x, 0.5)[1] for x in plot_range], '.-', label = "Ey")
plt.plot(plot_range, [du_dy(x, 0.5)[0] for x in plot_range], '.--', label = "dEx/dy")
plt.plot(plot_range, [du_dx(x, 0.5)[1] for x in plot_range], '.--', label = "dEy/dx")
plt.plot(plot_range, [du_dy(x, 0.5)[1] for x in plot_range], '.--', label = "dEy/dy")
plt.title("At y = 0.5")
plt.xlabel('x')
plt.legend()

plt.show()

The results look acceptable to meā€¦
Why is this simplistic way of solving this equation not popular ?

1 Like
6

The least squares approach, under certain circumstances, suffers from nonphysical loss. Solutions suffer non-physical dissipation, particularly strong where re-entrant edges are present. I have been investigating this off and on, and it seems that the LS operators are not optimal from the a norm-equivalence point of view, making convergence dependent on a finer discretization (or higher order) than would be needed for the traditional Galerkin approach.

I do not claim to fully understand the reasons for this yet, so if anyone out there has something to say on this issue, I would be very interested as well.

1 Like