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 nonmatching grid approach is also an interesting alternative [12].
In this article, we are more specifically concerned with solving threedimensional Maxwell’s equations, that are often used
to describe the physics of engineering problems, in their static or timedependent form, sometimes coupled with other equations (see an overview in [7]). Moreover, many structures that are to be modeled have a complex threedimensional 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 wellknown, when solving Maxwell’s equations in a nonconvex and nonsmooth 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 threedimensional boundary value problems in nonconvex domains is basically different from the twodimensional case and is often more difficult. Among many existing methods, Fourier Finite Element Method is an efficient method for solving problems in threedimensional 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 threedimensional axisymmetric singular domains with arbitrary data, we consider a situation in which the threedimensional (3D) Maxwell equations can be reformulated as twodimensional (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 finitedimensional 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 timedependent Maxwell equations in vacuum read in ,
(1)  
(2)  
(3)  
(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
(5) 
Since we are interested in the static Maxwell equations, we consider problems and solutions that are timeindependent, namely static equations. In other words, we assume that the explicit timedependence of the electromagnetic field in Maxwell’s equations vanishes. With nonvanishing charge and current densities, this assumption yields there are two divcurl problems, depending on the boundary condition.
The first one is, for a mean zero value righthand side in , such that and , and for a righthand side in :
Find such that
(6)  
(7)  
(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 saddlepoint 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 Friedrichstype inequality of the form . Equivalently, Eqs. (6)(8) can represent the stationary problem associated with Maxwell’s equations, namely the quasielectrostatic problem. This amounts to assuming that the timedependent parts are known, and that and .
With analogous notations, the second divcurl problem is, for a given in such that , and for a mean zero value in :
Find such that
(9)  
(10)  
(11) 
Similarly, this can model the quasimagnetostatic 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 saddlepoint 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 halfplane (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 halfplane is defined by the equation constant, and are Cartesian coordinates in this halfplane.
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 twodimensional one by assuming that , as made for instance in [6]. We continue here to deal with a threedimensional problem.
In these conditions, one can obtain the expressions of the static Maxwell equations simply by replacing into (68) and (911) the operators and by their cylindrical counterparts in the cylindrical coordinates , with the basis vectors , defined by
(12)  
(13) 
Similarly, the gradient operator in cylindrical coordinates is defined by
(14) 
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
and
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
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 seminorm 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:
(15) 
and similarly, for the magnetic field,
Find such that:
(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 twodimensional 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 FourierLaplace 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 twodimensional 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
resp.
and the truncated Fourier expansion of w at order
(17) 
We also consider the weighted Lebesgue space
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. :
Above, the operators for the mode are defined as:
(18) 
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
with .
Similarly, we introduce the analogous of for the Fourier coefficients, namely
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
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 (68) and (911) (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
(19) 
we get that each mode is the solution to the problem:
find such that, for all :
(20) 
where and denote the Fourier coefficients of the righthand 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 :
(21) 
Here again, and denote the Fourier coefficients of the righthand 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 righthand 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:
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
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
(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

For
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
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:
(23) 
For the space , the basis solves:
(24)
As a consequence of Prop.1, we readily get that the spaces and are satisfying
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 nonorthogonal complement for the modes . More details and numerical illustrations will be given in Section 6. Another choice would be to derive, for each , a modespecific 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
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 socalled 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
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:

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

only for , the conical singularity ( in Figure 2right), 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
(25) 
(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 cutoff 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
The righthand side of this equation can be computed analytically, by using the following expressions of and , involved in , see (19):
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
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
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
(27) 
where the term in (27) is useful to impose the boundary condition on the axis , and is, here again, a particular cutoff 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
(28) 
The righthand side of this equation is computed analytically by using the following expressions of and , involved in , see (19):
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
(29)  
where denotes the operator for (namely in the “fully” axisymmetric case), , the vector of a scalar field being defined by
In addition, the two bilinear forms and are defined by
and
Comments
There are no comments yet.