Gérard Orjubin
email: orjubin@univ-mlv.fr

Elodie Richalot

Member, IEEE Stéphanie Mengué

Member, IEEE Man-Fai Wong

Odile Picon

On the FEM Modal Approach for Reverberation Chamber Analysis

Keywords: Cavity, electromagnetic analysis, finite element methods, large-scale eigenproblem, reverberation chamber

We review the difficulties of a modal approach when modeling a Reverberation Chamber (RC) by the Finite Element Method (FEM). The numerical challenge is due to the large scale problem involved by the over-dimensioned cavity. Moreover, the field singularity on the stirrer has to be captured by the FEM. First the following issues are discussed: existence of nullfrequency solutions, convergence rate for h and p adaption, and formulation type in E or H field. Then the modal analysis is compared to the classical harmonic one. A focus is put on the field singularity at the source point.

I. INTRODUCTION

he Reverberation Chamber (RC) is a useful tool employed for Electromagnetic Compatibility measurement. It consists of a closed metallic cavity containing a metallic stirrer. High field levels can be generated in this low loss enclosure using modest source power, and allow susceptibility tests on electronic devices. The field within a cavity being essentially composed of standing waves, the moving metallic stirrer contributes to randomize the field. At high frequencies, the field can be statistically uniform and isotropic over a working volume, for a sufficient perturbation induced by a stirrer displacement, generally a rotation. This statistical approach has been validated experimentally for high frequencies [START_REF] Ladbury | Evaluation of the NASA Langley Research Center Mode-Stirred Facility[END_REF].

As illustrated by Bruns and Vahldieck in [START_REF] Bruns | A Closer Look at Reverberation Chambers-3-D Simulation and Experimental Verification[END_REF], modeling is helpful to better understand the RC stirring process at low frequencies. This author gives a complete overview of the different techniques already used, most of them being based on a point by point harmonic approach. As an exception, the modal analysis of the RC was first investigated by Bunting in [START_REF] Bunting | A Two-Dimensional Finite-Element Analysis of Reverberation Chambers[END_REF] for the 2D case, the author pointing out the numerical difficulties. Since this paper publication (1999), many improvements have been made in term of software and hardware: the main purpose of the present paper is to show how state-of-the-art techniques allow a resolution of the 3D problem on a simple PC.

The attention must be drawn on the fact that we deliberately restrict the study to a very basic RC, despite the importance of small geometric details, such as RC doors, as illustrated in [START_REF] Bruns | A Closer Look at Reverberation Chambers-3-D Simulation and Experimental Verification[END_REF]. As another restriction of our study, the frequency range concerns the first modes of the RC. Although these facts deserve a special attention, our previous work [START_REF] Orjubin | Modal Representation of a Reverberation Chamber for a Parametric Analysis of the Field Uniformity[END_REF] showed interesting properties of the modal technique, such as a possibility to study the influence of the quality factor on the RC stirring process, and a capability of analyzing the modes perturbation induced by the stirrer rotation. In fact, the present paper is dedicated to the fundamental and the numerical aspects of the modal method.

Mode determination is used for other applications than RC: for instance the modeling of dielectric resonator by the Finite Integral Technique (FIT) [START_REF] Feigh | Eigenmode Simulation of Electromagnetic Resonator Cavities with Gyrotropic Materials[END_REF] or the analysis of particle accelerator cavity by the Finite Element Method (FEM) [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF]. We used this latter modeling method implemented in the Open Source Pyfemax program [24] and in commercial Femlab® V3.0 software.

The lossless RC eigenmodes are directly determined by the Maxwell equations discretization that leads to a large scale eigenproblem.

We first compare the performance of the solvers Implicitly Restarted Arnoldi Method (IRAM) and Jacobi-Davidson (JD). Then a specific problem arising with FEM is exposed, i.e. the existence of null-frequency modes. A Tree-Cotree technique [START_REF] Trapp | Eigenvalue Computation by Means of a Tree-Cotree Filtering Technique[END_REF] has been proposed to discard these modes, using a topological analysis of the discretization mesh. We check in Section II that the technique used in [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF] also finds its origin in a topological matrix. Then a numerical example shows the reduction of the CPU time.

In Section III, the h and p adaption issues are investigated. As the field may be singular on stirrer edges, the consequences of these singularities on the convergence rate are presented in Section IV.

The excitation of this enclosure creates another field singularity at the source point. Its contribution to the modal expansion is examined and a regularization technique is implemented in Section V.

II. MODE DETERMINATION

The basic geometric model of the RC used in this study is depicted in Fig. 1: A planar metallic stirrer (1.5m×0.75m) is placed in a (3.10m×2.47m×3.07m) metallic rectangular cavity. The stirrer is centered, placed at 0.87m below the roof and parallel to it. Only one stirrer position is analysed, for which the stirrer edges are parallel to the cavity walls. In the main part of this text, results are presented for the frequency band 50~100 MHz, which contains 5 eigenmodes when considering the empty stirrerless cavity. For the eigenmode analysis (Sections II~IV), the antenna modeling is useless.

A. Large scale eigenvalue problem

This lossless problem verifies the double rotational equation for the complex E field:

     Ω ∂ = ∧ Ω =       = on 0 in 2 r r r r r r r r E n E E c E t o R t o R λ ω . ( 1
)
where Ω designates the cavity inside.

Projecting the field on N global edge elements where the stiffness R and mass M matrices are defined by:

     = = ∫∫∫ ∫∫∫ Ω Ω dV w w m dV w t o R w t o R r k j k j k j k j r r r r r r . . .
This generalized eigenproblem leads to a direct determination of modes: each eigenvalue is linked to a resonant frequency, whereas the eigenvector corresponds to the mode cartography. R and M are real symmetric matrices whose dimension (number of DoF) is high, say N = 10 5 . In such a case, specific techniques exist to solve this large eigenproblem.

Based on a projected solution space, Implicitly Restarted Arnoldi Method (IRAM), Implicitly Restarted Lanczos (IRL) and Jacobi-Davidson (JD) algorithms are suited to the solving of (2), particularly because R and M are sparse matrices.

Looking for a small set of eigenmodes, one has to specify the eigenvalue target σ and the number of modes sought. In our case, 6 eigenmodes are sought close to the target σ = 2.

Femlab® features an IRAM solver, which is associated to a Shift-and-Invert spectral transformation: a factorization of the real symmetric matrix R-σ M is therefore necessary [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF], for which several factorization algorithms are available on Femlab software. After numerical tests, the Spooles factorization algorithm is retained, the performance of this direct method being such that the iterative techniques are useless.

On another hand, JD is a factorization-free algorithm, implemented in Pyfemax [START_REF] Geus | The Jacobi-Davidson algorithm for solving large sparse symmetric eigenvalue problems with application to the design of accelerator cavities[END_REF] [START_REF] Hill | Aperture excitation of electrically large, lossy cavities[END_REF]. That program also includes hierarchical second order H(Rot)-conform elements and a powerful Algebraic MultiGrid (AMG) technique.

The comparison of these two solvers is given in Fig. 2, for finding the 6 RC modes present in the 50~100 MHz band. The first conclusion is about the CPU time: 10 3 s for 10 5 DoF is fairly good when the platform is a mere 1.2 GHz PC. The second observation is that the algorithmic complexity makes the MultiGrid technique a must for larger problems, with a complexity of 1.3 to be compared to the complexity of 1.8 obtained with Femlab.

B. Structure of the discrete solution space

Taking the divergence of (1) leads to

0 ) ( 2 =         E Div c r ω . (3)
Non-null-frequency modes are solenoidal (null-divergence): more generally, the solutions can be either solenoidal (E sol ) either null-frequency irrotationnal (E irr ). These two kinds of solutions are orthogonal in the following sense:

( )

0 . 0 , = ⇔ = ⇔ ⊥ ∫∫∫ Ω dV E E E E E E sol irr ol s irr ol s irr r v r v r v (4)
The solution space H 0 (Rot) can be split in two infinite dimension spaces, a solenoidal and an irrotational space, as indicates the Helmholtz decomposition:

H 0 (Rot) = W sol ⊕ GradH 0 (5) where W sol = {u ∈ H 0 (Rot) / Div(u)=0} GradH 0 = {u ∈ H 0 (Rot) / Rot(u)=0}.
When using the so called edge elements introduced by Nédélec [START_REF] Nédélec | Mixed finite element for IR3[END_REF], the splitting is respected in a discrete sense: to get an idea about this fact, we recall some properties of these elements. If L i and L j are the Lagrange scalar functions associated to i and j points of a given tetrahedron, then the edge element is

i j j i ij L d a Gr L L d a Gr L w r r r - = . ( 6
)
These elements (named local as they are associated to a given tetrahedron) verify Div(w ij )=0, whereas the global elements w k permit to describe correctly the non-null-divergence solutions E sol . On the other hand, we have

1 14 13 12 ) ( L d a Gr w w w r r r r = + + - . (7)
This means that Grad(L 1 ) is an algebraic sum of the edge elements linked to node 1. Equation ( 7), concerning local elements L j and w ij , can be generalized to elements ϕ j and

w k : ] 1 [ 1 P j w g d a Gr k N k kj j ∈ ∀ = ∑ = r r ϕ (8)
where N is the edge number and P the node number.

The (g kj ) permit to define a topological sparse matrix G that can be considered as a discretized gradient. If

( ) N v v v v ..., , , 2 1 = r stands for ∑ = = P i i i v V 1 ϕ , then v r G corresponds to ∑ ∑∑ ∑∑ = = = = = = = P i i i P i k N k i i k N k k P i i i k d a Gr v w v g w v g 1 1 1 1 1 ϕ r r r . (9)
More properties can be presented: if Φ h is the scalar space of the nodal functions, [START_REF] Geus | The Jacobi-Davidson algorithm for solving large sparse symmetric eigenvalue problems with application to the design of accelerator cavities[END_REF] means that Grad(Φ h ) is included in the vectorial space W h of Nédélec functions. In fact, this inclusion property can be extended to a discrete decomposition [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF]:

W h = W h sol ⊕ Grad(Φ h ). ( 10
) Because of R definition (2), it is easy to check that RG = 0. (11)
This means that Grad(Φ h ) space corresponds to the R kernel, i.e. the null-frequency solutions e 0 , as R e 0 -λ M e 0 = 0 implies λMe 0 = 0 or λ =0 (M being positive defined). Thus W hsol is the space of researched non-null-frequency modes.

The relation [START_REF] Manges | A Generalized Tree-Cotree Gauge for Magnetic Field Computation[END_REF] can be interpreted in the more conceptual way using the differential form spaces W n , as suggested by Bossavit in 1988 [START_REF] Bossavit | A rationale for the "Edge-Elements" in 3-D fields computations[END_REF]. These spaces respect inclusion relations, that are outlined in the De Rham complex:

3 D 2 Ro 1 G 0 W W W W    →     →    →  . ( 12
)
The differentiation of W 0 (resp. W 1 ) via the discrete Grad (Resp. Rot) operator is included in the W 1 (Resp W 2 ) space, and so on. The Rot(Grad)=0 relation is similar to [START_REF] Manges | A Generalized Tree-Cotree Gauge for Magnetic Field Computation[END_REF]. In the case of trivial topology, W 1 (i.e. W h ) is a sum of two subspaces: Grad(W 0 ) (the kernel of discretised Rot operator) and the complementary space of null-divergence elements, as indicated by [START_REF] Bossavit | A rationale for the "Edge-Elements" in 3-D fields computations[END_REF]. Summing up, the use of Nédélec elements leads to a singular R matrix, and the R kernel corresponds to the nullfrequency solutions of (2). Taking into account the frontier conditions, the number of these non physical solutions is equal to the number of mesh interior points and can be estimated [START_REF] Manges | A Generalized Tree-Cotree Gauge for Magnetic Field Computation[END_REF] as 14% of the R dimension, N. Because N is as large as 10 5 , one understands the need to deal with these solutions.

C. Elimination of the Null-frequency solutions

) 13
As the null-frequency solutions e irr belong to the R kernel, i.e. to the G image [START_REF] Manges | A Generalized Tree-Cotree Gauge for Magnetic Field Computation[END_REF], we can write e irr = G v. The Morthogonality of the solutions of (2) yields:

( ) ( )

0 0 = ∀ ⇒ = ∀ sol T sol T e v v e v v r r r r r r T T C M G (14) G M C C T = = ⇒ with 0 sol e r . ( 15
)
Respecting [START_REF] Dibben | Nodal and Edge Element Analysis of Inhomogeneously Loaded 3D Cavity[END_REF], the resonant frequencies will be non-null: condition (15) can be implemented in the solver through various techniques: see [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF] for JD or [START_REF] Perepelitsa | Finite-Element Analysis of Arbitrary Shaped Cavity Resonators Using H 1 (curl) Elements[END_REF] for IRL.

It is interesting to note that (15) can be derived in another way, enforcing DivE = 0 with a Lagrange multiplier p. This is Kikuchi method, recalled in [START_REF] Arbenz | Solving Maxvell eigenvalue problems for accelerating cavities[END_REF], that leads to the weak variational statement: find E Є H 0 (Rot) and p Є H 0 1 such that

( ) ( ) ( ) ( )      = = + 0 , , , , 1 1 1 1 p d a Gr E E E E p d a Gr E t o R E t o R r r r r r v r r r r λ . ( 16
)
for any test functions E 1 and p 1 .

Discretizing E on edge elements (vector e of unknowns) and p on Lagrange elements (vector x), it comes

[ ] [ ] x e 0 0 0 M x e 0 C C R T       =       λ (17)
where ( )

j i j i d a Gr w c ϕ r r , = .
C is actually the same as given by [START_REF] Dibben | Nodal and Edge Element Analysis of Inhomogeneously Loaded 3D Cavity[END_REF], as shows ( 18):

( ) ∑ ∑ ∑ = =       = k k k i j k k i j k k k j k i j i m g w w g w g w c r r r r , , . (18)
As x is a Lagrange multiplier, expressing DivE=0 (i.e. nonnull-frequency) leads to the same condition [START_REF] Dibben | Nodal and Edge Element Analysis of Inhomogeneously Loaded 3D Cavity[END_REF]. Whatever the method is, [START_REF] Dibben | Nodal and Edge Element Analysis of Inhomogeneously Loaded 3D Cavity[END_REF] permits to eliminate the non-physical modes.

To illustrate the advantage of this restricted research, we used the Pyfemax program that includes an implementation of (15) in the JD direct solver, called SAUG, and compare the computation times with SAUG algorithm and without this option (direct solver). In both cases, the first steps of assembling and conditioning the matrix are unchanged. The computation time (s) for determining the modes in the 50~100 MHz band are given in Table I (mesh with 17 10 3 tetrahedra). For both dicretization orders, elimination of irrotational solutions by SAUG algorithm permits to save 10% off the solver computation time.

III. NODAL AND EDGE ELEMENTS COMPARISON

The nodal elements, also called Lagrange elements, are well-known to generate the so-called spurious modes. In fact, these modes appearance is due to the equation discretization. In other words, nodal elements can successfully be used to solve another equation than (1). This is done in this Section.

The solution space structure yields the existence of irrotational solutions (e.g. Gradφ). As the Gradφ discretization must be continuous (class C 0 ), this implies that the φ discretization must be of class C 1 . This condition can not be respected by using first order Lagrange elements for the φ discretization: this is the origin of the spurious solutions, as mentioned in [13, Sec. 7.2.1].

In the next paragraph is illustrated the penalty technique [START_REF] Geus | The Jacobi-Davidson algorithm for solving large sparse symmetric eigenvalue problems with application to the design of accelerator cavities[END_REF] that overcomes this problem: it is then possible to compare the precision of nodal and edge elements, and check that these latter should be preferred, as reported by Bardi et al. in [14]. Furthermore, the numerical advantages of second order elements are examined as an extention of a previous work [START_REF] Dibben | Nodal and Edge Element Analysis of Inhomogeneously Loaded 3D Cavity[END_REF].

A. Penalty method

This method consists in adding a penalty term to (1):

     Ω ∂ = = ∧ Ω = - - on 0 and 0 in 0 E Div E n E E Div d a Gr s E t o R t o R r r r v r r r r r r r λ (19)
Taking the divergence of [START_REF] Conciauro | Advanced Modal Analysis[END_REF] and defining e = DivE leads to:

     Ω ∂ = Ω = = ∆ - on 0 in e e e s e µ λ (20)
The smallest eigenvalue µ 0 of this scalar equation can be determined using the FEM with first order Lagrange elements.

If µ < µ 0 , the unique solution of ( 20) is trivial, i.e. e = DivE = 0. Thus, the irrotational solutions of ( 19) are eliminated if the penalty parameter s is taken to be higher than λ /µ 0 .

B. Application to the rectangular stirrerless cavity

For instance, we found µ 0 = 3.8. λ = 3.8 being the higher eigenvalue in the frequency band 50 to 100 MHz, it means that s >1 guarantees that the modes check DivE = 0 and are not spurious. To illustrate this fact, we look for the eigenmodes present in this band, knowing that there are only 5 physical modes in the stirrerless cavity. For various penalty parameter values, the number of found modes is reported in Fig. 3. The spectral pollution by spurious modes is avoided for this frequency band if s > 1. As a particularity, we note that s = 1 is a right value for determining the fundamental of a cavity and leads to a closed-form equivalent of ( 19) as: 0

r r r v = + ∆ E E λ . ( 21
)
An advantage of the study of the stirrerless rectangular cavity is that resonant frequencies are known analytically. This enables a convergence study, for a h-adaption scheme, using Femlab (Lagrange elements) and Pyfemax (edge elements) software: the results are presented in Fig. 4. where p is the element order, a is the regularity coefficient and h the mean mesh width. As the field inside the rectangular cavity is regular, the convergence rate is 2p. An important point is the good performance of the H(Rot)-conform (both first and second order edge elements) in term of convergence rate, and in term of numerical properties, as the R matrix is sparser with edge elements.

However, these results would be far different when modeling the RC. Firstly, it has been shown that the penalty method must be modified in case of re-entrant corner [START_REF] Costabel | Weighed Regularization of Maxwell Equations in Polyhedral Domains[END_REF]. Actually, the use of the simplistic [START_REF] Conciauro | Advanced Modal Analysis[END_REF] for a cavity including a stirrer gives no trustable result. A second remark concerns the regularity of the electric field: the stirrer creates a field singularity, that limits the convergence rate seen in ( 22) (as the regularity a is low). When extending this property to H(Rot)conform elements, it comes that the convergence is slower if the cavity contains a stirrer. This is investigated in Section IV.

IV. FIELD SINGULARITY ISSUES

A. h and p adaption

As already seen, the penalty method [START_REF] Conciauro | Advanced Modal Analysis[END_REF] can not be employed to model the cavity with the stirrer. This one strongly modifies the field and a high gradient may exist. For instance, the first eigenvector cartography, portrayed in Fig. 5, shows high edge effects at the stirrer vicinity. The monotonous convergence illustrates the advantage of order 2 over order 1 (p-adaption). It also proves the benefit of a mesh finely refined close to the stirrer, called "adapted mesh" in Fig. 6. Note that this a priori non-uniform mesh is not stricto sensu an adapted mesh to a given mode. In fact hp adaption technique is suited to driven frequency analysis or to the determination of a single eigenmode, but not to our study of a set of eigenmodes.

Another remark must be done concerning the modeling of singularities: instead of higher order regular (polynomial) elements, it seems advantageous to use singular elements at a lower order [START_REF] Graglia | Singular Higher Order Complete Vector Bases for Finite Methods[END_REF].

B. E and H formulation

Some interesting conclusions may be drawn from the theory of singularity. In the case of a plane wave diffraction by a infinite plane, Van Bladel [START_REF] Van Bladel | Singular Electromagnetic Fields and Sources[END_REF] has shown that the E field components are equivalent to 1 0

1 < <      ∝ ∝ - ν ν ν r E r E T Z , ( 23
)
where E Z and E T are respectively the components along the edge direction and transverse, r is the distance to the edge and ν is the singularity exponent. In case of diffraction on a flat sector, the author gives a comparison (Fig. 7) between the singularity exponent ν (for E) and τ (for H):

0 360° ν 1 0.5 δ τ 180°
δ Fig. 7. E and H singularity exponents for a flat sector, from [18, Fig. 5.17].

For a convex corner angle, the H singularity is better than the E singularity. This property of the diffracted waves seems to be valid for the study 1 of modes, as shows Table II in the case of the RC: The first conclusion is that the mesh has to be farely refined to get close results for E and H formulations, as Mesh3 corresponds to more than 2 10 5 DoF. Furthermore, this convergence study indicates a clear precision enhancement for the H formulation. Although these results only correspond to the 6 first modes of our RC, it is worth noting that it is consistent with theoretical aspects outlined in Fig. 7.

V. COMPARISON TO HARMONIC ANALYSIS

Once the modes are determined, it becomes possible, in a single post-treatment, to get the spectrum of the field when the RC is excited by an antenna. A comparison with a direct harmonic modeling is presented in this section.

A. Modal expansion

For an enclosure excited by a harmonic driven current J of pulsation ω, (1) becomes

     Ω ∂ = ∧ Ω         + - = on 0 in 2 0 r r r r r r r r E n E c J j E t o R t o R ω µ ω . (24)
The general solution of ( 24) is expanded as follows [19, §2.3]:

stat n M n n irr n n n sol n n n E E E E r r r r ∑ ∑ ∑ - = ∞ = ∞ = + + = 1 1 1 1 δ β α , (25)
where E n sol are the solenoidal solutions of (1), E n stat are the static solutions due to the RC topology 2 and E n irr are the gradient of solutions of Dirichlet problem:

     Ω ∂ = Ω         - = ∆ on 0 in 2 φ φ ω φ c . (26)
In the same way that (1) yields ω 0n sol , (26) gives ω 0n irr , parameters having no physical reality, since these modes are associated to 0Hz.

In the case of perfectly conducting walls, the coefficients of (25) are easily determined:

( ) ( ) ( )            = = - = stat n n n irr n n n sol n sol n n n E J j W E J j W E J j W r r r r r r , 2 1 , 2 1 , ) ( 2 1 2 2 ω δ ω β ω ω ω α (27)
where W n is the energy of the mode n.

The total field is imaginary and presents singularities at the solenoidal resonant frequencies. It is worth mentioning that the irrotational eigenvalues of (26) do not appear in these formulas: the singularities of irrotational solutions appear at 0Hz and these components are negligible near solenoidal resonant frequencies. It gives a reason to employ the common solenoidal expansion [START_REF] Orjubin | Modal Representation of a Reverberation Chamber for a Parametric Analysis of the Field Uniformity[END_REF][20] [START_REF] Cerri | Source Stirring Mode for Reverberation Chambers[END_REF] as:

sol n n n E E r r ∑ ∞ = ≅ 1 α . ( 28
)
Another interesting property is that the spectrum of an excited and closed cavity is no more discrete. To illustrate this point, we calculate the field (28), adopting from now the following: the measured point is P(1,1,1), and the source current J(1,-1,1) is located at source point S(0.4,2,2), for coordinates given in meters and current amplitude in Amperes par square meter. As seen in Fig. 8, the spectrum of an excited lossless cavity is no more a Dirac sum and this fact has been confirmed by a direct harmonic analysis. The modes are not only excited at resonant frequencies. When the coupling between the cavity and the internal source is taken into account, the cavity response becomes non-null between resonant frequencies. It is worth noting that in real life, this phenomenon is hidden by another one: As the RC is excited by an antenna coupled to the exterior, the RC cannot be considered as lossless.

B. Losses implementation

To suppress the frequency singularities, losses have to be modeled. A perturbation technique is employed, stating that resonant frequencies obtained in the lossless analysis (1) are substituted by complex resonant frequencies [START_REF] Conciauro | Advanced Modal Analysis[END_REF] as

      + = Q j f f sol n sol n 2 1 0 . ( 29
)
Thus, the coefficients α n are finite at the resonant frequencies, and for weak losses (Q >> 1), they are given close to the resonant frequency sol n f 0 by: ( )

sol n sol n sol n sol n sol n sol n n E J f f f Q j f f W f Q r r , 2 1 4 0 0 0 0 - + - ≅ π α .(30)
The frequency dependent term of (30) is a narrow band approximation of a second order oscillator response.

Instead of using a mode-dependant quality factor, which would represent losses originating only from Joule effect, we consider an averaged Q factor that takes into account losses from various origins. The main advantage is a capability to model antenna losses that are preponderant at low frequencies. To illustrate this fact, we first evaluate the Joule losses by the following approximation [START_REF] Hill | Aperture excitation of electrically large, lossy cavities[END_REF] δ

S V Q J 2 3 = , ( 31
)
where S and V are the cavity surface and volume, and δ the skin depth. Leaking power by antenna can be modeled [START_REF] Hill | Aperture excitation of electrically large, lossy cavities[END_REF] by

3 2 16       = c f V Q Ant π . ( 32
)
Finally the composite quality factor is defined [START_REF] Hill | Aperture excitation of electrically large, lossy cavities[END_REF] by

Ant J C Q Q Q 1 1 1 + = . ( 33
)
The evolution of the different quality factors versus the frequency is graphed in Fig. 9. Figure 9 shows clearly that the Joule losses can be neglected in the very low part of the spectrum: losses are mainly due to antenna leaking. Considering two antennas (transmission and reception), Q is set to a low value (e.g. Q = 60) for the 50~100 MHz band. At higher frequencies, around 260 MHz, Q is set to 1000.

In Fig. 10 a comparison is given for the 50~100 MHz band between the null-divergence expansion (28) and a harmonic analysis. In this later one, losses are simulated by adding surface impedances on the cavity walls, as for Joule losses. To permit a comparison between the results obtained using both methods, the surface impedance of the walls is chosen to yield a mean Q factor of 60 over the frequency band. A good agreement is shown, as well as the expected f 0n /2Q resonant frequency offset [19 p. 34]. The numerous harmonic simulations necessary to approximate the resonance peaks illustrate the advantage of the modal method in matters of computational effort 3 . The method is then applied to a higher frequency: 16 modes are sought around 260 MHz, and the electric field is expanded on these 16 modes using (28) and (30), with Q = 1000. The result is plotted in Fig. 11, as well as many harmonic determinations, wherein the Surface Impedance of the walls is chosen to yield a mean Q factor of 1000 over the frequency band. Besides this satisfactory comparison, two phenomena can be observed in Fig. 11. First the use of a same Q factor for all modes leads to acceptable differences for the peaks height. Secondly, the behavior at the frequency band edges gives an indication of the convergence speed of expansion (28). Figure 11 shows that the modal expansion on 16 modes is reliable in the 254~267 MHz band, thus at a 1 MHz from the edges of the frequency band on which the eigenvalues are calculated. This 1 MHz value corresponds to roughly 4f 0 /Q MHz, where f 0 = 260 MHz.

Figures 1112have illustrated and validated the use of modal expansion (28): a single eigenmode determination followed by a modal expansion provides a field description in a narrow frequency band, which would need many harmonic analysis. However, one can note that (28) is only an approximation of (25), since it does not include the irrotational terms: its validity is thus reserved for measure points located far from the source point.

C. Singularity extraction

The source term creates a field singularity, which cannot be expanded on the only solenoidal terms (28). To illustrate the difference between (25) and (28), we conduct a harmonic analysis at f = 76.62 MHz, which corresponds to the RC third resonance in Fig. 10. The result is given in Fig. 12. As the modes are well separated, the main term of the solenoidal expansion (28) corresponds to the third mode. α n=3 being real at the resonance, the field real component (Fig. 12) has the same pattern as TM 110 . In fact, a small imaginary field component does exist, due to the other modes contributions. For example, the real part of the electric field at point P is 36.7 and the imaginary part is 4.6. These values are deduced from both harmonic modeling and modal expansion (28). Thus (28) is a good approximation of (25) for a common interior point of the cavity, far from the antenna.

Actually, we notice on Fig. 13 a high imaginary 4 field region close to the antenna, and a fast decrease of this imaginary part when the distance to the antenna increases. This phenomenon can only be explained by irrotational terms in (25). The first idea would be to expand this irrotational part on the basis of solutions of (26). For this, we discretize the scalar equation (26) by Lagrange 2 elements and portray in Fig. 14 the first irrotational mode: Of course, these irrotational modes do not present any singularity near the source point: this explains why the convergence is known to be very slow [19, Ac. 2.8], i.e. the number of irrotational modes to take into account in the expansion (25) is very high. In place of expanding this irrotational part of E field, we illustrate the singularity extraction technique [START_REF] Borji | Rapid Calculation of the Green's Function in a Rectangular Enclosure With Application to Conductor Loaded Cavity Resonators[END_REF] in a simple way. The Green function singularity in the RC close to the antenna is the same as the singularity of free-space Green function. On a similar way, we just suppose that the field close to the antenna corresponds to the free-space field E 0 . This field is numerically approximated by a harmonic analysis, wherein Matched Impedance is implemented on cavity walls. Figure 15 depicts the result of the substraction of this free-space field E 0 to the total field observed in Fig. 13. As confirms Fig. 15, the field is regularized in the antenna vicinity, and can then be expanded on solenoidal modes.

VI. CONCLUSION

In this paper a Reverberation Chamber is modeled using the Finite Element Method (FEM) modal approach. For this 3D FEM model, the large number of Degrees of Freedom implies a restriction to the lower part of the spectrum. In order to analyze the RC close to the Lowest Useable Frequency (LUF), high performance solvers must be used. A comparison is given between Open Source Pyfemax and Femlab commercial code, and the issue of the null-frequency solutions is exploited to decrease the CPU time.

Second, the RC stirrer presence induces a field singularity that worsens the convergence rate: an investigation on the element kind and order proved the advantage of the second order edge elements. This modal approach is validated by a comparison to harmonic analysis, and a regularization technique is illustrated for the field in the source region.

.

Fig. 1 .

 1 Fig. 1. RC with its rectangular metallic stirrer.

  vector e of DoF. These elements are named global, as they are not related to a specific tetrahedron. Then, the Galerkin method yields the algebraic equation

Fig. 2 .

 2 Fig. 2. Computation time for determining the 6 modes present in the band 50 to 100 MHz. With Femlab, first order edge elements are used, and the solving is performed using IRAM solver and the Spoole factorization algorithm. With Pyfemax, hierarchical second order elements are used with an AMG technique, and the solving involves a JD algorithm.

From

  the definitions of M (2) and of the field orthogonality (4), one can write 0

Fig. 3 .

 3 Fig. 3. Influence of the penalty parameter on the number of modes in the 50-100 MHz band for the stirrerless cavity, using first order Lagrange elements.

Fig. 4 .

 4 Fig. 4. Mean precision for determining the resonant frequencies of the stirrerless cavity in 50 -100 MHz band.The nodal element convergence rate fits to the theory: the error over the eigenvalue determination is [13, Sec. 2.5.3]:

Fig. 5 .

 5 Fig. 5. Cartography of the mode 1 of the RC.

Figure 6

 6 Figure 6 represents the first resonant frequency obtained with Pyfemax for various meshes, using H(Rot)-conform elements of order 1 (Nédélec) and order 2:

Fig. 6 .

 6 Fig. 6. Convergence study of the RC first mode, using H(Rot)-conform elements.

Fig. 8 .

 8 Fig. 8. Norm(E) at point P obtained by modal method without loss implementation.

Fig. 9 .

 9 Fig. 9. Evolution of quality factor versus frequency, for the studied RC.

Fig. 10 .

 10 Fig. 10. Norm(E) at point P obtained by harmonic and modal methods, 50 ~100 MHz frequency band (Q # 60).

Fig. 11 .

 11 Fig. 11. Norm(E) at point P obtained by harmonic and modal method , close to 260 MHz (Q = 1000).

Fig. 12 .

 12 Fig. 12. Slice cartography of real(Ez) at resonant frequency f = 76.62 MHz obtained by harmonic analysis.

Fig. 13 .

 13 Fig. 13. Slice cartography of 20log[abs(Imag(Ez))] at resonant frequency f=76.62 MHz, obtained by harmonic analysis (notice the E scale in dB).

Fig. 14 .

 14 Fig. 14. Slice cartography of the first irrotational mode.

Fig. 15 .

 15 Fig. 15. Slice cartography of |Imag(Ez-E0z))| at resonant frequency f =76.62 MHz.

TABLE I COMPARISON

 I

		Order 1 elements	Order 2 elements
		2 10 4 DoF	1.1 10 5 DoF
	Assembling	11	17
	Conditionning SuperLu: 37	SuperLu: 37
	Direct solver	28	405
	SAUG	25	352

OF PYFEMAX SOLVER OPTION FOR ELIMINATION OF NULL-FREQUENCY MODES (CPU TIMES IN S.)

In fact, only one physical static mode exists among the irrotational modes found by the solver, as the stirrer presence yields a double connected frontier (M=2).

The matrices issued from harmonic analysis are complex symmetric, factorized by UMFPACK algorithm using Femlab.

The cavity excitation is modeled by a current distribution.