# Numerical solution to the 3D Static Maxwell equations in axisymmetric singular domains with arbitrary data

We propose a numerical method to solve the three-dimensional static Maxwell equations in a singular axisymmetric domain, generated by the rotation of a singular polygon around one of its sides. The mathematical tools and an in-depth study of the problem set in the meridian half-plane are exposed in <cit.>, <cit.>. Here, we derive a variational formulation and the corresponding approximation method. Numerical experiments are proposed, and show that the approach is able to capture the singular part of the solution. This article can also be viewed as a generalization of the Singular Complement Method to three-dimensional axisymmetric problems.

There are no comments yet.

## Authors

• 5 publications
• 1 publication
12/22/2020

### Numerical approximation of singular-degenerate parabolic stochastic PDEs

We study a general class of singular degenerate parabolic stochastic par...
02/10/2020

### Numerical solution to stress distribution of a hole with corners on infinite plane

In this paper, we consider the stress of a hole with the given fourfold ...
09/24/2019

### A sharp error estimate of piecewise polynomial collocation for nonlocal problems with weakly singular kernels

As is well known, using piecewise linear polynomial collocation (PLC) an...
05/18/2021

### Numerical Solution for an Inverse Variational Problem

In the present work, firstly, we use a minimax equality to prove the exi...
05/11/2021

### Preconvergence of the randomized extended Kaczmarz method

In this paper, we analyze the convergence behavior of the randomized ext...
08/07/2020

### A Probabilistic Numerical Extension of the Conjugate Gradient Method

We present a Conjugate Gradient (CG) implementation of the probabilistic...
11/20/2019

### Balanced truncation model reduction for 3D linear magneto-quasistatic field problems

We consider linear magneto-quasistatic field equations which arise in si...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

There is a need to simulate electromagnetic wave phenomena of increasing complexity, leading to the development of more general and efficient numerical methods. Indeed, a plethora of engineering problems requires to simulate numerically devices working with or within electromagnetic fields.

This article is part of the efforts made in the more general framework of boundary value problems with singularity in their solutions, caused by the presence of geometrical singularities, that is reentrant corners or edges on the boundary of a domain, or similarly by a change in the type of boundary conditions [28, 29]. From a more physical point of view, they are called singularities, since they can generate in their vicinity very strong fields that have to be taken into account, and are very difficult to compute. Moreover, as illustrated in [9], the inability to properly handle these singularities may have dramatic consequences on the physical phenomenon one wants to study.

In this context, several authors have proposed to use methods that “extract” the singular part of the solution near these singularities, or to apply mesh refinement toward these singularities, in the case of weak singularities (roughly speaking, that belong to a regular enough space like ). This allows to construct numerical methods that are able to catch the singular behavior of the solution. The non-matching grid approach is also an interesting alternative [12].

In this article, we are more specifically concerned with solving three-dimensional Maxwell’s equations, that are often used to describe the physics of engineering problems, in their static or time-dependent form, sometimes coupled with other equations (see an overview in [7]). Moreover, many structures that are to be modeled have a complex three-dimensional geometry that presents a surface with reentrant edges and/or corners, namely singularities.

There exist many methods to solve the Maxwell equations numerically [36]. One can mention the edge finite element method [37, 38], more recently, the class of discontinuous Galerkin method introduced by [33], or adaptive finite element method in two dimensions, as proposed in [17]. However, it is interesting for some applications to have a continuous approximation of the solutions, for instance when coupling the Maxwell equations with other equations, like Vlasov’s one [1], [10].

As it is well-known, when solving Maxwell’s equations in a non-convex and non-smooth domain with a continuous approximation, the discretized spaces are always included in a closed, strict subspace of the space of real solutions, see the seminal work of [15, 16] for theoretical justifications, and more recent developments by Costabel and Dauge (see among others [26]). Consequently, it is not possible to approximate the singular field and needs special treatment, even for static problems [24]. In this case, mesh refinement techniques fail. The Singular Complement method (SCM) [6, 9] addresses this problem by explicitly adding some singular complements to the space of solutions, see also [2, 3, 8].

Numerical solution of three-dimensional boundary value problems in non-convex domains is basically different from the two-dimensional case and is often more difficult. Among many existing methods, Fourier Finite Element Method is an efficient method for solving problems in three-dimensional prismatic or axisymmetric domains, even for other equations, see for instance [13] for Stokes equations. The method uses the Fourier expansion in one space direction associated to a finite element approach in the other two space dimensions, see, among others [18], [31], [29], [35], or [32] for interface problems.

In the present work, that extends the SCM to three-dimensional axisymmetric singular domains with arbitrary data, we consider a situation in which the three-dimensional (3D) Maxwell equations can be reformulated as two-dimensional (2D) models. This principle was also derived in [20, 21] for the Poisson equation in a prismatic or axisymmetric geometry. More precisely, the computational domain boils down to a subset of

, with respect to the cylindrical system of coordinates. Nevertheless, the electric and magnetic fields, and other vector quantities, still belong to

. Hence, the electromagnetic field is the solution to an infinite set of 2D equations, and as a result a set of 2D variational formulations, obtained by Fourier analysis.

This paper is organized as follows: in a first section, we recall the Maxwell equations and their formulation in an axisymmetric domain. Then we present the principle of the 2D space reduction, based on the use of a Fourier transform in

. This reduces 3D Maxwell’s equations to a series of 2D Maxwell’s equations, depending on the Fourier variable . This allows us to compute the 3D solution by solving several 2D problems, each one depending on . Even if the solution remains singular for each in the 2D domain, we will be able to decompose it into a regular and a singular part (see Section 4). The regular part belongs to a regular space and will be computed by a standard finite element method. The singular part, that belongs to a finite-dimensional subspace, will be handled following the same principle as in the SCM. This is the subject of Section 5. In the last Section, numerical examples are proposed to illustrate the feasibility of the method.

In the remainder of this paper, we write vector fields or spaces with boldface. Similarly, names of function spaces of scalar fields usually begin by an italic letter, whereas they begin by a bold letter for spaces of vector fields (for instance, or ).

## 2 Maxwell’s equations in an axisymmetric domain

### 2.1 The static Maxwell equations

Let be a bounded and simply connected Lipschitz domain in , its boundary, assumed for simplicity to be a connected boundary, and n the unit outward normal to . If we let and be, respectively, the speed of light and the dielectric permittivity, the time-dependent Maxwell equations in vacuum read in ,

 ∂E∂t−c2curlB=−1ε0J, (1) ∂B∂t+curlE=0, (2) divE=ρε0, (3) divB=0, (4)

where E is the electric field, B the magnetic flux density, and J the charge and current densities. These quantities depend on the space variable and on the time variable .

These equations are supplemented with appropriate boundary conditions. In this article, we assume that the boundary is a perfect conductor, so that the electromagnetic field satisfies

 E×n=0 and B⋅n=0 on the boundary Γ. (5)

Since we are interested in the static Maxwell equations, we consider problems and solutions that are time-independent, namely static equations. In other words, we assume that the explicit time-dependence of the electromagnetic field in Maxwell’s equations vanishes. With non-vanishing charge and current densities, this assumption yields there are two div-curl problems, depending on the boundary condition.

The first one is, for a mean zero value right-hand side in , such that and , and for a right-hand side in :
Find such that

 curlE=fE in Ω, (6) divE=gE in Ω, (7) E×n|Γ=0. (8)

The boundary condition on is imposed by the condition (8) (cf. [27]). In order to prove the existence and uniqueness of the solution E to (6)-(8), a possible way is to reformulate these equations as a saddle-point formulation, and to check that the Lagrange multiplier is equal to 0 (see [23], [7] Chap.6 for details). Assuming the connectivity of the boundary is required here, since the use of the saddle point approach needs to use a Friedrichs-type inequality of the form . Equivalently, Eqs. (6)-(8) can represent the stationary problem associated with Maxwell’s equations, namely the quasi-electrostatic problem. This amounts to assuming that the time-dependent parts are known, and that and .

With analogous notations, the second div-curl problem is, for a given in such that , and for a mean zero value in :
Find such that

 curlB=fB in Ω, (9) divB=gB in Ω, (10) B⋅n|Γ=0. (11)

Similarly, this can model the quasi-magnetostatic Maxwell’s equations by assuming and . The fact that has a mean zero value stems from (11). The existence and uniqueness of B can also be inferred by using a saddle-point approach. In both cases, the existence and uniqueness result can be achieved thanks to the Weber inequality [40], see details in [4] or in [7] Chap.6.

### 2.2 Formulation in an axisymmetric domain

Now we make the supplementary assumption that is an axisymmetric domain, limited by the surface of revolution . We denote by and their intersections with a meridian half-plane (see Figure 1). One has , where either when is a closed contour (i.e. does not contain the axis), or is the segment of the axis lying between the extremities of . The natural coordinates for this domain are the cylindrical coordinates , with the basis vectors . A meridian half-plane is defined by the equation constant, and are Cartesian coordinates in this half-plane.

However, although the domain is assumed to be axisymmetric, the symmetry of revolution is not assumed for the data. In these conditions, the problem can not be reduced to a two-dimensional one by assuming that , as made for instance in [6]. We continue here to deal with a three-dimensional problem.

In these conditions, one can obtain the expressions of the static Maxwell equations simply by replacing into (6-8) and (9-11) the operators and by their cylindrical counterparts in the cylindrical coordinates , with the basis vectors , defined by

 divu=1r∂∂r(rur)+1r∂uθ∂θ+∂uz∂z (12) curlu=(1r∂uz∂θ−∂uθ∂z)er+(∂ur∂z−∂uz∂r)eθ+1r(∂∂r(ruθ)−∂ur∂θ)ez (13)

Similarly, the gradient operator in cylindrical coordinates is defined by

### 2.3 Variational formulations in 3D

We now introduce the variational formulations of the problem, which can be applied independently of the (non) convexity of the domain . Let us define the function spaces, with classical notations: the usual norm and scalar product of are denoted by and respectively. We shall also need to use the following Sobolev spaces and norms

We introduce likewise

 H0(curl,Ω)={v∈% H(curl,Ω):v×n|Γ=0}

and

 H0(div,Ω)={v∈H(div,Ω):v⋅n|Γ=0}.

The electric and magnetic field naturally belongs respectively to the spaces

The spaces and are compactly embedded in [40], [25]. Consequently, when the boundary is connected, one can define an equivalent scalar product and norm on and as

 a(u,v):=(curlu,% curlv)+(divu,divv),∥u∥X=∥u∥Y:=a(u,u)1/2

In other words, the -norm is uniformly bounded by the X and the Y norm for elements of and respectively. This is the Weber inequality, that basically claims that in or in , the semi-norm is a norm equivalent to the canonical one.

We have now to derive the (augmented) variational formulations associated to these problems. Following a classical approach, we first take the dot product of equations (6) (resp: (9)) by , (resp: , ) and integrate over , then add the variational form of the divergence equation for E (resp: B). This gives the variational formulations:
Find such that:

 a(E,F)=(fE,curl% F)+(gE,divF),∀F∈X(Ω), (15)

and similarly, for the magnetic field,
Find such that:

 a(B,C)=(fB,curl% C)+(gB,divC),∀C∈Y(Ω). (16)

Existence, uniqueness and continuous dependence with respect to the data of these variational formulations follow from the application of usual techniques, see for instance [7].

## 3 Principle of two-dimensional space reduction

Since we consider non axisymmetric data, we can not perform to reduce the 3D space problem to a 2D space one. However, we will use that the domain is axisymmetric. The scalar and vector fields defined on will be characterized through their Fourier series in , the coefficients of which are functions defined on . Note that such a technique together with the Fourier-Laplace transform is also used for stability analysis of numerical schemes [19] solving Maxwell’s equations. Let us also emphasize that the time dependent part of the problem is not involved here. What is explained below is the principle of the two-dimensional space reduction. For this reason, we do not mention the time variable in the Fourier series, which can be easily added the case occurring. For instance, we will consider for a given function (resp: for a vector field ), the Fourier expansion

 w(r,θ,z)=1√2π∑k∈Zwk(r,z)eikθ,

resp.

 w(r,θ,z)=1√2π∑k∈Zwk(r,z)eikθ

and the truncated Fourier expansion of w at order

 w[N](r,θ,z)=1√2πN∑k=−Nwk(r,z)eikθ. (17)

We also consider the weighted Lebesgue space

 L21(ω):={w measurable on ω:∬ω|w(r,z)|2rdrdz<∞}

which is the space of Fourier coefficients (at all modes) of functions in .

Let us now examine the space of relevant Fourier coefficients for the electromagnetic fields. One easily checks that for , resp. such that , there holds:

while for , resp. :

 divw=1√2π∑k∈Zdivkwkeikθ resp. curlw=1√2π∑k∈Zcurlkwkeikθ.

Above, the operators for the mode are defined as:

As explained in [14], the regularity of the function and of the components of w in the ad hoc Sobolev spaces are characterized by the regularity of the Fourier components and , for .

As a consequence, a function v belongs to if and only if, for all , its Fourier coefficients belong to the space defined by

 X(k)(ω)={vk∈L21(ω),curlkvk∈L21(ω),divkvk∈L21(ω),vk×n∣γb=0}

with .

Similarly, we introduce the analogous of for the Fourier coefficients, namely

 Y(k)(ω)={vk∈L21(ω),curlkvk∈L21(ω),divkvk∈L21(ω),vk⋅n∣γb=0}

and .

An important property concerning these spaces is proved in [22], Prop.2.9:

###### Proposition 1.

The spaces and are independent of , for

As a consequence, it will be sufficient to compute the singular subspaces only for the modes , while the modes will be used to compute all the higher modes .

### 3.1 Variational formulation in 2D for each k

Our aim is now to apply this space dimension reduction to the 3D equations, to derive the corresponding 2D formulations satisfied by Fourier coefficients (, ) for each mode .

More precisely, we use the linearity of Maxwell’s equations (6-8) and (9-11) (or equivalently of their variational formulations) together with the orthogonality of the different Fourier modes in . This implies that the Fourier coefficients and of E and B are solutions to variational formulations similar to (15) and (16), with the operators and . Consequently, introducing the operator as follows

 ak(u,v)=(% curlku,curlkv)+(divku,divkv), (19)

we get that each mode is the solution to the problem:
find such that, for all :

 ak(Ek,F)=(fkE,curlkF)+(gkE,divkF), (20)

where and denote the Fourier coefficients of the right-hand sides and respectively, that depend only on .

Similarly, the magnetic field B being solution to (16), its Fourier coefficients satisfy the formulation, for each mode :
find such that, for all :

 ak(Bk,C)=(fkB,curlkC)+(gkB,divkC). (21)

Here again, and denote the Fourier coefficients of the right-hand sides involved in the equations of the magnetic field.

For an analysis of the truncation error of the Fourier expansion, we refer the interested reader to [22], [39]. Basically, the convergence of the truncated solution (see (17)) toward E is in for a given norm, the value of depending on the regularity of the right-hand sides . Similar results are available for .

## 4 Decomposition in regular/singular parts

Up to now, in the same spirit as in [35], we have reduced the 3D Maxwell equations to a series of 2D Maxwell equations, depending on the Fourier variable . Nevertheless, the two dimensional domain being singular (see Figure 1), we have now to deal with this singularity. The construction of the numerical method will be based on theoretical results proved in [4], for the case , and in [22] for the general case.

For our purpose, we first consider, for each Fourier mode , the weighted Sobolev space that contains functions such that . Then, we introduce the regularized spaces and subspaces of , defined by:

 XR(k)(ω):=X(k)(ω)∩H1(k)(ω) and YR(k)(ω):=Y(k)(ω)∩H1(k)(ω).

We then have the following property ([22] Lemma 6.2 and §6)

###### Proposition 2.

The regularized spaces and are closed, respectively, within and .

In these conditions, for a singular domain, the spaces of solution and can be decomposed in

 X(k)(ω)=XR(k)(ω)⊕XS(k)(ω) and Y(k)(ω)=YR(k)(ω)⊕YS(k)(ω).

The subspaces and are the spaces of solutions in case of a regular domain, whereas and are singular subspaces, equal to for a regular domain. As a consequence, the electromagnetic field solution and can also be decomposed into a regular and a singular part, says

 Ek=EkR+EkS,Bk=BkR+BkS. (22)

Moreover, these singular subspaces are of finite dimension, the dimension of which depending on the number of singularities in the domain .

Hence, one can compute a numerical approximation of and by a standard numerical method, for instance a -conforming finite element method. The difficulty coming from the singular parts and , we have now to derive a way to characterize these singular fields. For this purpose, we refer to the following property, see for instance [5], Prop. 3.2 of [6] or §7.1 of [22]. Let be the solution to the following equation, which involves a Legendre function: . Its value , and we have

###### Proposition 3.

The singular spaces and are of finite dimension, namely

• For

 dimYS(k)(ω):=NB= number of % reentrant edges, dimXS(k)(ω):=NE=NB+ number of % conical points with vertex angle >πβ.
• For

 dimYS(k)(ω):=NB=dimXS(k)(ω):=NE= number of reentrant edges.

By introducing now and the basis of and for a given Fourier mode , we get that the singular parts of the Maxwell’s equations solution can be decomposed into

 EkS=NE∑j=1kjExkS,j%andBkS=NB∑j=1kjBxkS,j,

where and are constants we will have to determine. This will be detailed in Section 6.

We present now the characterization of the singular basis and . For simplicity, in what follows, we will assume that there is only one singularity, that is , and we will drop the index . In these conditions, we are looking for the equations satisfied by and . Following [4] §5.1 and §5.2, [6] §3.1 or [22] §7.1, we obtain that they can be characterized via their variational formulation. Indeed, they are solution to the following homogeneous formulations

• For the space , the basis solves:

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ak(xkS,F)=0,∀F∈XR(k)(ω)xkS×n|γb=0,xkS⋅n|γa=0, (23)
• For the space , the basis solves:

 {ak(ykS,C)=0,∀C∈YR(k)(ω)ykS⋅n|γ=0,γ:=γb∪γa (24)

As a consequence of Prop.1, we readily get that the spaces and are satisfying

 XS(k)(ω)=XS(2)(ω),YS(k)(ω)=YS(2)(ω), for |k|≥2.

As recalled above, this stabilization property has a fundamental consequence on the numerical method, based on the decomposition of and in a regular and singular subspace: it will be sufficient to compute the singular basis and only for and not for all , the modes also serving as a non-orthogonal complement for the modes . More details and numerical illustrations will be given in Section 6. Another choice would be to derive, for each , a mode-specific orthogonal basis. The interested reader will find a comparison for the Poisson equation in [34].

## 5 Computation of singular basis

### 5.1 The case of xkS

We now present the numerical method to compute the singular part . The basis solves the singular and homogeneous problem (23) in . Consequently, if we try to solve it with a standard finite element method, we will get a zero solution.

To overcome this difficulty, we rather use that the so-called principal part of the singularity S - the part that makes singular - does not depend on the Fourier mode . In these conditions, can be decomposed into

 xkS=xkS,reg+S

where denotes the regular part of , that can be computed by a classical finite element method. Note also that in the electric case, as recalled in Prop.3, there exist two kinds of geometrical singularities S:

1. for all , the one that exists in the neighborhood of a reentrant edge of , that is a reentrant corner of , ( in Figure 2-left), that we will denote ,

2. only for , the conical singularity ( in Figure 2-right), that we will denote , that exists in the neighborhood of a conical vertex with an aperture greater than the limit vertex angle for (so that ).

Figure 2 shows the notations associated to these singularities. In particular, denotes the local polar coordinates centered at the reentrant edge , the corresponding angle being called , . For the conical point , are the local polar coordinates centered at this point, with the origin of on the -axis.

In these conditions, the principal part can be written as , whereas the principal part at the conical point (if any) can be expressed as .

As above, denotes also here the Legendre function of index , where is the index such that . The expressions of the and in the basis are given by

 Se=−raαρα−1⎛⎜⎝sin((α−1)ϕ−ϕ0)0cos((α−1)ϕ−ϕ0)⎞⎟⎠, (25)
 Sc=νρν−1⎛⎜⎝Pν(cosϕ)cosϕ−P1ν(cosϕ)sinϕ0Pν(cosϕ)sinϕ+P1ν(cosϕ)cosϕ⎞⎟⎠, (26)

Remark that the term in (25) has been introduced to impose the boundary condition on the axis , and can be viewed as a particular cut-off function. Moreover, the singularity appearing only for , it is related to a “fully axisymmetric” case, which was already treated in [6], and will not be considered in the following (included numerical examples of Section 6).

Under these circumstances, we can compute by solving the following variational formulation

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ak(xkS,reg,v)=−ak(Se,v)∀v∈XR(k)(ω),xkS,reg×n|γb=−Se×n|γb,xkS,reg⋅n|γa=0,

The right-hand side of this equation can be computed analytically, by using the following expressions of and , involved in , see (19):

 curlkSe=αaρα−1⎛⎜⎝−ikcos((α−1)ϕ−ϕ0)cos((α−1)ϕ−ϕ0)iksin((α−1)ϕ−ϕ0)⎞⎟⎠,divkSe=−2αaρα−1sin((α−1)ϕ−ϕ0).
###### Remark 1.

By construction, see for instance Eqs.(23)-(24), the singular complement is orthogonal with respect to the form . Now, for the numerical implementation, it is also possible to orthonormalise the basis , for each singularity , and to compute the basis vectors which are orthogonal to one another and to the regular space (for ). The same is true for the magnetic case.
Now, concerning the principal part of the singularity S, it is the same, whatever the singular basis that we consider, the orthogonalisation process modifying only the regular part of , and not S. Computing such an orthogonal complement requires an additional computational effort on the one hand. On the other hand, the variational formulations (20)-(21) are easier to solve because they contains fewer terms, some canceling due to orthogonality.

### 5.2 The case of ykS

Let us turn our attention to the computational method for the singular part . This time, the basis solves the system of equations (24). This problem is singular and homogeneous in : if we try to solve it with a standard finite element approach, we will get a zero solution.

To overcome this difficulty, we rather use, as for the case of , that the principal part of the singularity S, i.e. the part that makes singular, does not depend on the Fourier mode . In these conditions, can be decomposed into

 ykS=ykS,reg+S

where denotes the regular part of , that can be computed by a classical finite element method. The expression of S in the basis is given by

 S=−raαρα−1⎛⎜⎝cos((α−1)ϕ−ϕ0)0−sin((α−1)ϕ−ϕ0)⎞⎟⎠, (27)

where the term in (27) is useful to impose the boundary condition on the axis , and is, here again, a particular cut-off function. Note also that in the magnetic case, there is no singularity due to the presence of conical vertex.

In these conditions, the function will be computed by solving the following variational formulation

 ⎧⎨⎩ak(ykS,reg,v)=−ak(S,v)∀v∈YR(k)(ω),ykS,reg⋅n|γ=−S⋅n|γ. (28)

The right-hand side of this equation is computed analytically by using the following expressions of and , involved in , see (19):

 curlkS=αaρα−1⎛⎜⎝iksin ((α−1)ϕ−ϕ0)−sin((α−1)ϕ−ϕ0)ikcos((α−1)ϕ−ϕ0)⎞⎟⎠,divkS=−2αaρα−1cos((α−1)ϕ−ϕ0).

For the practical purpose of the computation, it is useful to express the bilinear form , depending on the values of . Indeed, recall that, in our approach, we will have to compute the singular basis and only for . Performing a simple integration by parts shows that

 ak(u,v) = a0(um,vm)+k2(umr,vmr)+(curluθ,curlvθ)+k2(uθr,vθr) (29) + ıkB(u,v)+ıkC(u,v),

where denotes the operator for (namely in the “fully” axisymmetric case), , the vector of a scalar field being defined by

 curlw:=−∂zwer+r−1∂r(rw)ez.

In addition, the two bilinear forms and are defined by

 B(u,v):=∫γb(um⋅n)¯vθ−uθ(¯vm⋅n)dγ,

and

 C(u,v):=∫∫ω2(uθ¯vr−