Axisymmetric Electromagnetic Wave Propagation

Computation Using the Constrained Interpolation

Profile Scheme with Large Time Steps
Jean Porto, Paul-Quentin Elias

To cite this version:

Jean Porto, Paul-Quentin Elias. Axisymmetric Electromagnetic Wave Propagation Computation
Using the Constrained Interpolation Profile Scheme with Large Time Steps. IEEE Transac-
tions on Antennas and Propagation, Institute of Electrical and Electronics Engineers, In press,
�10.1109/TAP.2020.3044582�. �hal-03206748�

HAL Id: hal-03206748

https://hal.archives-ouvertes.fr/hal-03206748
Submitted on 23 Apr 2021

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Axisymmetric Electromagnetic Wave Propagation
Computation Using the Constrained Interpolation
Profile Scheme with Large Time Steps
Jean Porto, Paul-Quentin Elias

Abstract—Conventional explicit time-domain methods used [9] under the condition of identical cell size [10]. A proof
to solve Maxwell’s equations are reliable and robust but are of the weak stability of the scheme is presented in [11].
conditionally stable and often require the fulfillment of the The CIP scheme was combined with the Method of Char-
condition CFL ≤ 1. While this is acceptable for many ap-
plications, in some instances where the Maxwell’s equations acteristic (MoC) to get an accurate simulation of Maxwell’s
are solved alongside systems with slower propagation velocities, equations by Ogata et al. in 2006 [12]. Since then, several
explicit methods prove costly. This is the case for non-relativistic research developments have been published around the use of
electromagnetic Particle-In-Cell methods which are required to the CIP method for electromagnetics [13], [14]. Indeed, the
study plasma thrusters. Several algorithms have been proposed CIP method has several interesting features for the computa-
to retain a nearly explicit formulation using large time steps to
achieve higher CFL values. Among these is the semi-Lagrangian tion of electromagnetic fields in time domain, among them
Constrained Interpolation Profile method. While the ability of its higher order means lower dispersion and lesser points
this method to handle CFL > 1 has been demonstrated for required per wavelength. It can be used with variable cell
planar 2D–3D cases, this has not been done for 2D cases with sizes or subgridding techniques [15]. Finally, being a semi-
cylindrical symmetry. In this paper, a procedure is presented to Lagrangian technique, it can be used with larger CFL numbers
compute the electromagnetic wave propagation in 2D domains
with cylindrical symmetry using the Constrained Interpolation than classical explicit FD techniques.
Profile (CIP) method. The CIP scheme is extended for CFL ≥ 1 For the particular case of cylindrical systems, Tanaka et al.
cases, and a ghost node method is proposed to deal with the axis [13] proposed a method to deal with the radial terms with the
singularity and with the wall boundary condition. The results are CIP method. By considering that the computational domain
compared to the fields of a Hertzian dipole and with a coaxial can be seen as a medium with radially varying impedance,
cable, and they show a good agreement.
the characteristic variables can propagate radially with the
Index Terms—Constrained Interpolation Profile (CIP) method, semi-Lagrangian method taking into account the transmission
electromagnetic wave propagation, 2D cylindrical coordinate and reflection of waves due to the impedance variation. We
system, Courant-Friedrichs-Lewy condition (CFL) condition.
have tested this procedure and it has proved to be accurate
for the treatment of electromagnetic problems in a cylindrical
I. I NTRODUCTION configuration. However, the procedure described by Tanaka et
Finite-difference (FD) methods are one of the most extended al. [13] is only limited to CFL ≤ 1.0.
and successful techniques for computational electromagnetic The CIP scheme features, coupled to the capacity to get
problems [1], [2]. One of the advanced techniques among accurate results in cylindrical coordinate systems, could be
the finite-difference methods is a stable third-order accurate an interesting tool for the simulation of plasma propulsion
and non-dissipative scheme initially developed in the field of devices as they can often be represented by an axisymmetric
Computational Fluid Dynamics [3], known as the Constrained configuration [16]–[18]. An electromagnetic solver for cylin-
Interpolation Profile (CIP) method. It is a semi-Lagrangian drical simulations based on the CIP method would have the
scheme which circumvent the Courant-Friedrichs-Lewy (CFL) advantages of updating the fields at the same grid location,
stability condition [4], [5]. This feature allows computations which simplifies its integration to the Particle-In-Cell codes
with CFL values ≥ 1.0, as can be seen in [6] and [7] where used to model these low pressure plasma sources [19], [20].
the authors performed simulations using a CFL value of 2.6 In addition, the characteristic form used in the CIP method
in a Cartesian coordinate system. It considers not only the for electromagnetic computation offers an easy handling of
electromagnetic fields, but also their spatial derivatives, there- different types of boundary conditions such as a Perfectly
fore suppressing instabilities and providing lower numerical Matched Layer (PML) and conducting or dielectric surfaces.
dispersion even when using coarse grids and large time steps Nevertheless, with the methodology described in [13], it will
[8]. It has been shown that it provides higher accuracy than the still be limited by the CFL ≤ 1.0 condition.
conventional finite-difference time-domain (FDTD) method In this paper we extend the results of Tanaka. et al. [13] to
get a version of the 2D cylindrical CIP scheme that could be
Jean Porto is with Sorbonne Université, Observatoire de Paris, PSL Re- used with a CFL condition greater than one, thus providing a
search University, LERMA, CNRS UMR 8112, F-75005, Paris, France. mean to considerably speed-up the calculations. In addition,
Jean Porto and Paul-Quentin Elias are with DPHY, ONERA,
Université Paris Saclay, F-91123 Palaiseau – France (e-mail: we define a systematic strategy for the treatment of boundary
jean carlos.porto [email protected] ; [email protected]). conditions with this extended approach, using ghost nodes
outside the computational domain to take into account the
singularity at the symmetry axis in a cylindrical coordinate
system and conducting surfaces.
To investigate the accuracy of the proposed scheme, we
compared its results with those obtained using analytical ex-
pressions [21] to describe a Hertzian dipole. Another compari-
son seeking to test the stability of our approach was performed
with the steady-state propagation of electromagnetic waves
inside an open-ended coaxial wave guide. In that case the
time-domain simulation with the CIP method is compared with
the result using the COMSOL Multiphysics ™ [22] software
based on a frequency-domain finite-element analysis.
The paper is organized as follows. The sections II-A and
II-B recall the main features of the CIP method adapted for
2D cylindrical systems. Then, section II-C extends the method
proposed in [13] for cases where CFL ≥ 1. Section II-D
provides a simple method to account of the axis singularity
and perfectly conducting boundaries in the frame of the CIP
method with CFL ≥ 1. Finally, section III considers two tests
cases for the validation of the proposed method.

II. T HE CIP METHOD IN A TWO - DIMENSIONAL

CYLINDRICAL COORDINATE SYSTEM

A. Fundamentals of the CIP Scheme

Fig. 1. Fields update on the CIP method for the propagation of an electro-
Just as the finite-difference time-domain (FDTD) method, magnetic wave where ECIP and BCIP are the interpolated fields between
the CIP considers the electromagnetic fields on each grid the grid cells: a) CFL ≤ 1.0, b) 1.0 < CFL ≤ 2.0.
point. However, it also propagates the spatial derivatives values
for each field through an additional advection equation. This
feature produces less numerical dispersion or instabilities [8] are interpolated using a cubic polynomial function from the
because for a given propagating wave, its values between the values at the endpoints of the cell, generating what we will
grid points are interpolated using not only the wave informa- call FCIP and GCIP located at the square marker as seen on
tion at the grid but also its spatial derivatives. It allows to get Fig. 1. If we consider the right-propagating wave (+z) f and
a better approximation of the wave all over the computational g at grid points zi and zi−1 , the interpolated values between
domain and to maintain its original shape through the whole those grid points are given by:
simulation. Consider the following advection equation for a
one-dimensional problem : FCIP (z) =ai (z − zi )3 + bi (z − zi )2 + gi (z − zi ) + fi
∂F (z) (3)
∂f ∂f GCIP (z) = = 3ai (z − zi )2 + 2bi (z − zi ) + gi
+c =0 (1) ∂z
∂t ∂z
Where ai and bi are coefficients calculated based on the
It represents the propagation of the a wavefield f at a functions f and g values at the grid points as follows:
constant speed c. The CIP method solves not only eq. 1, but
also another differential equation for its spatial derivative along
its propagation direction, obtained by differentiating eq. 1: gi + gi−1 2(fi − fi−1 )
ai = +
(−∆z)2 (−∆z)3
∂g ∂g ∂f (4)
+c =0 ; g= (2) 3(fi−1 − fi ) 2gi + gi−1
∂t ∂z ∂z bi = −
(−∆z)2 (−∆z)
The advection phase may be better understood thanks to
Fig. 1. It represents the field update process on the CIP method B. Electromagnetic Fields Numerical Analysis in the Two-
where at each time step, the value at any given grid point is Dimensional Cylindrical Coordinate System
updated using the value the field had before traveling up to that The Maxwell-Faraday and Maxwell-Ampère equations in a
point the distance given by the time step and its propagating 2D axisymmetric cylindrical coordinate system can be written
speed (c∆t). It can be summarized as follows: consider a wave as follows if the medium is nonconducting:
propagating along the positive direction of the z-axis (+z): the
function f and its spatial derivative g at each time step (∆t) ∂Er ∂Bϕ 1
+ c2 = − Jr (5)
and at each grid point (zi ) are obtained by shifting its values ∂t ∂z 
from the position zi −c∆t, where c is the wave velocity. If the ∂Ez ∂Bϕ c2 1
required values to be shifted are not on the grid points, they − c2 = Bϕ − J z (6)
∂t ∂r r 
 ∂Er ∂Ez ∂Bϕ
− + =0 (7)
∂z ∂r ∂t
Where  is the permittivity of the vacuum, µ is its per-
meability, the electric field has a radial and a longitudinal
component, Er and Ez respectively, and the magnetic field
is oriented in the azimuthal direction Bϕ . This set of equation
can be expressed in matrix form as:

∂W ∂W ∂W
+Λ +Γ =S (8)
∂t ∂r ∂z

     
Bϕ 0 0 −1 0 1 0
W =  Er  , Λ =  0 0 0  , Γ = c2 0 0
Ez −c2 0 0 0 0 0
h iT
2 Fig. 2. Computational domain for the simulations in cylindrical coordinate
S = 0, − 1 Jr , cr Bϕ − 1 Jz systems. The colourless nodes represent ghost nodes outside the domain used
for the boundary condition at the symmetry axis Z.
We now applied the change of variable proposed by Tanaka
et al. [13] to get rid of the azimuthal magnetic field term
in the right hand side of eq. 8: Iϕ = µr Bϕ . The key idea • Along the radial direction (r-axis): Rp and Rm
here is to consider that the domain impedance is radially Z Z
varying. Around a radial grid point, Tanaka et al. assume Rp = Ez − Iϕ ; Rm = Ez + Iϕ (13)
rj rj
that the impedance is constant by part, thus removing the 1/r
dependence and transforming the medium as a graded-index ∂Rp ∂Ez Z ∂Iϕ ∂Rm ∂Ez Z ∂Iϕ
= − ; = + (14)
medium. The final form of our system can be written as: ∂r ∂r rj ∂r ∂r ∂r rj ∂r

Z      The subscript p (plus) stands for a wave propagating in

rj Iϕ 0 0 −c 0 c 0 the positive direction of its axis (right-going wave), and m
W =  Er  , Λ =  0 0 0  , Γ =  c 0 0 (minus) in the negative direction (left-going wave).
Ez −c 0 0 0 0 0
Once equation 10 is solved, we can treat the non-advection
T
S = 0, − 1 Jr , − 1 Jz term and obtain the field at the time step n+1 by adding the


pµ current source J as follows:

Where Z =  and rj is the radial distance from the
axis of symmetry to the grid point. We will treat the source ∆t
En+1 = E∗∗ − J (W∗∗ → Wn+1 ) (15)
term S as a non-advection term that will only be taken 
into account once the solution of the equation 8 is found
under the hypothesis of S = 0. To do so, thanks to the C. Wave Propagation with Multi-Layer Radial Scattering
change of variable previously introduced, we can now use the The computational domain used for our calculations in a
dimensional-splitting method [23]: 2D cylindrical coordinate configuration is presented in Fig. 2.
The longitudinal axis goes from zi=0 to zi=N z , and the radial
∂W ∂W axis from rj=0 to rj=N r . The empty circles are ghost nodes
+Λ =0 (Wn → W∗ ) (9)
∂t ∂r helpful for the boundary conditions definition.
∂W ∂W Thanks to the dimensional-splitting method presented in
+Γ =0 (W∗ → W∗∗ ) (10) equations 9 and 10, the problem can be treated along each
∂t ∂z
propagating axis. Regarding the longitudinal axis, the charac-
Where Wn is the field at time step n, W∗ represents the teristic and its derivatives described in equations 11 and 12
field after propagation along the radial direction, and W∗∗ is are advected along the positive and negative direction of the
the field after the propagation along the longitudinal direction. z-axis using the interpolation polynomial described in section
We will therefore express our system as the following set of II-A and as seen in the schematic example of Fig. 1.
propagating characteristics with its spatial derivatives along However, special treatment is required for the radial prop-
the propagation axis: agation because of the dependence of the Z/rj coefficient to
• Along the longitudinal direction (z-axis): Lp and Lm the radial distance from the longitudinal axis. Tanaka et al.
Z Z [13] has demonstrated that a radial scattering technique for
Lp = Er + Iϕ ; Lm = Er − Iϕ (11) the radial propagation produces accurate results. In order to
rj rj
extend this procedure for cases with CFL ≥ 1.0, let us begin
∂Lp ∂Er Z ∂Iϕ ∂Lm ∂Er Z ∂Iϕ by considering a right-going wave Rp propagating along the
= + ; = − (12) radial direction and a CFL ≤ 1.0. For a given node position
∂z ∂z rj ∂z ∂z ∂z rj ∂z
Fig. 3. Radial scattering procedure with CFL ≤ 1.0.

(zi , rj ) in the mesh, we can update its value at time step n+1
(Rpn+1j ) based on Rpnj−1 and Rpnj . The space between two
grid points (with radial distances to the symmetry axis rj and
rj−1 ) is considered as containing the interface between two
different medias (located at ∆r/2) each one described by a
fictitious impedance (Z/rj and Z/rj−1 ).
In Fig. 3, each medium is represented either as a grey zone
or as a hatched one. At the interface between different media,
as proposed by Okubo et al. in [24], the wave propagation Fig. 4. Multi-layer scattering of the advected characteristics in the radial
across this interface can be treated using the standard electro- direction. Field update for a right-going wave Rp using 1.0 < CFL ≤ 2.0.
magnetic interface conditions between two dielectrics. Thus,
each incoming wave will be divided into a ’R’-wave reflected can be generalized to greater CFL numbers. We will use an
back into the original media and a transmitted ’T’-wave that intermediate stage described as n + 12 , as seen in Fig. 4, where
propagates to the second medium. We will therefore call radial the required propagating characteristics will be calculated with
scattering the procedure used to obtain the wave RS which the following steps:
is the sum of a ’T’ wave coming from rj−1 and an ’R’-wave
coming from rj . The resulting wave RS can be thought of as 1) It starts with a radial scattering of the wave Rpnj−2
an equivalent wave coming from the left in the same medium. through the medias with impedance Z/rj−2 (hatched
This RS wave is used to perform the CIP propagation. This zone) and Z/rj−1 (grey zone) using equations 16. The
variable RS is shown in Fig. 3, and it is represented as the sum of the transmitted part of Rpnj−2 from the striped
sum of the two double-arrow lines in the grey zone. It can be zone to the grey one (solid line with double-arrow), and
obtained as follows for a right-going radial propagating wave: the reflected part of Rmnj−1 back into the grey zone
(dashed line with double-arrow) is called RSpj−2 .
2) This right-going characteristic RSpj−2 and the right-
RS =T+ Rpnj−1 + R− Rmnj going wave Rpnj−1 are both used to interpolate the
∂(RS) ∂Rpnj−1 ∂Rmnj (16) characteristic value between rj−2 and rj−1 using eq.
=T+ + R−
∂r ∂r ∂r 3. The result of this CIP interpolation is shifted to
With a transmission coefficient T+ = rj +rj−1
2r
and a reflec- the grid position rj−1 and it is an intermediate-stage
j−1 n+ 1
r −rj−1
tion coefficient R− = rjj +rj−1 . We check that T+ + R− = 1. characteristic called RSpj−12 .
Similarly, for a left-going radial propagating wave: 3) A radial scattering for the left-going wave Rmnj+1 using
equations 17 produces RSmj+1 . A CIP interpolation
between RSmj+1 and Rmnj , and a shift of the result of
RS =T− Rmnj+1 − R+ Rpnj
this interpolation to the position rj allows us to get the
∂(RS) ∂Rmnj+1 ∂Rpnj (17) n+ 1
=T− + R+ intermediate-stage left-going wave RSmj 2 .
∂r ∂r ∂r 4) Finally, once the two intermediate characteristics are
2r r −r n+ 1 n+ 1
Where T− = rj+1j+1 j+1 j
+rj and R+ = rj+1 +rj . As rj → ∞, obtained ( RSpj−12 and RSmj 2 ), the radial scattering
T → 1 and R → 0, meaning for example that for a right procedure is repeated with the right-going (RSpj−12 )
n+ 1
going wave RS → Rpj−1 , which is equivalent to the classical n+ 1
CIP method in a planar Cartesian coordinate system. and the left-going (RSmj 2 ) characteristics between
t = n + 12 and t = n + 1. The result of this scattering
Under the condition of CFL > 1.0, a multi-layer radial is shifted to rj and it gives us the update for Rpn+1
j .
scattering procedure is required in order to update the right- Note that two interpolations are performed in this procedure,
going characteristic Rpn+1
j . For this discussion we will restrict one of them is using RSpj−2 and Rpnj−1 , and the other
ourselves to the case CFL ≤ 2, but the method described below RSmj+1 and Rmnj . For a left-going wave, the entire pro-
cedure is symmetrical, starting from the rightmost point rj+2 fields: α+ and β+ for a conducting wall at r ≤ rj (eq.
and moving leftward to rj with the same multi-layer scattering. 19) receiving left-going waves Rm , and α− and β− for a
This procedure can be extended to higher CFL numbers. conducting wall at r ≥ rj (eq. 20) receiving right-going
First of all, the stencil changes as a function of the CFL waves Rp .
number. The higher the CFL, the greater the stencil extent • Perfect Conductor (PEC) wall parallel to the r-axis
spreading over as many grid cells as the rounded upper value To keep a zero tangential electric field on the conducting
of the CFL number. Having defined the stencil scope, the steps wall, therefore ensuring Er = 0, we can impose the fol-
1 and 2 are performed for the leftmost grid nodes, then the lowing conditions for any given real number k to produce
radial scattering procedure is repeated for RSp and RSm over a mirrored wavefield with respect to the conducting wall:
the remaining grid cells until the wave quantities reach rj . For
Er (z−k , rj ) = −Er (zk , rj )
example, for a CFL = 2.5, the radial scattering is performed
using a stencil scope going from rj−3 to rj+2 . ∂Er (z−k , rj ) ∂Er (zk , rj )
=
∂z ∂z (21)
Bϕ (z−k , rj ) = Bϕ (zk , rj )
D. Boundary Conditions
∂Bϕ (z−k , rj ) ∂Bϕ (zk , rj )
Different types of boundary conditions can be easily im- =−
posed with the use of ghost nodes as seen in Fig. 2. The ∂z ∂z
main idea is to update the values of the ghost nodes based Once the ghost nodes have been updated, the corresponding
on the waves reaching the boundary, and then proceed to a characteristics can be computed on those locations, and its
propagation back into the computational domain. Thus, no values can be propagated up to the inner nodes of the compu-
special stencil is required for the boundary nodes which uses tational domain applying the same multi-layer radial scattering
the same stencil as nodes in the domain. The rules for the procedure previously described.
ghost nodes update are the following:
• Symmetry axis
III. N UMERICAL RESULTS
Based on [25], the ghost nodes of the symmetry axis A. Hertzian Dipole
(smaller than rj=0 ) can be obtained by applying a simple We simulated the electromagnetic radiation of a Hertzian
transformation rule to the incoming wave: we multiply dipole with a 2D cylindrical coordinate system as seen in Fig.
by -1.0 any radial derivative, and also any radial or 2, by imposing a Gaussian function for the source term Jz
azimuthal component of a vector quantity (letting the along the symmetry axis with a width of σz = 2∆z and σr =
axial component unchanged). Under this rule, for any 4∆r where the mesh spacing ∆r = ∆z is fixed to 40 µm.
given real number k, we obtain: The dipole is located along the symmetry axis at r = 0 and
Ez (zi , r−k ) = Ez (zi , rk ) z = zd = 0.5∗Nz , as can be seen in Fig. 5, and it oscillates at
∂Ez (zi , r−k ) ∂Ez (zi , rk ) a frequency of 300 GHz. The expression for Jz is as follows
=− (with J0 = 2487.75 A.m−2 ):
∂r ∂r (18)
Bϕ (zi , r−k ) = −Bϕ (zi , rk )  2
(z − zd )2

∂Bϕ (zi , r−k ) ∂Bϕ (zi , rk ) r
= Jz (r, z, t) = J0 sin(ωt) exp − 2 − (22)
∂r ∂r σr σz2
• Perfect Conductor (PEC) wall parallel to the z-axis
Three different values of CFL were tested: CFL = 1.0,
Under the constraint of keeping a zero tangential electric
CFL = 1.5 and CFL = 2.0. The number of points per
field and normal magnetic field on a conducting wall lo-
wavelength was 25. The simulations results are compared with
cated at (zi , rj ), we can impose the following conditions:
the analytical values calculated using the equations presented
in Appendix A based on [21]. Figures 6 and 7 show a good
Rp (zi , rj−k ) = α+ Rm (zi , rj+k ) + β+ Rm (zi , rj ) (19) agreement among the three cases with a longitudinal view
Rm (zi , rj+k ) = α− Rp (zi , rj−k ) + β− Rp (zi , rj ) (20) along the line (z = 0 → Nz , r = 0.5 ∗ Nr ), and a radial
view along (z = 0.25 ∗ Nz , r = 0 → Nr ). This position for
the radial view was chosen over the alternative located at the
(rj + ∆r)(2rj − ∆r) middle point of the computational domain because the fields
α+ = −
(rj − ∆r)(2rj + ∆r) reach extremely high values in the zone close to the dipole.
2rj ∆r A radial view along the middle point of the computational
β+ = − domain would not be detailed enough to appreciate at the
(rj − ∆r)(2rj + ∆r)
(rj − ∆r)(2rj + ∆r) same time the two scales of the solution, close to the dipole
α− = − and faraway from it.
(rj + ∆r)(2rj − ∆r)
In Fig. 5, we can also notice the presence of an absorbing
2rj ∆r
β− = + boundary condition imposed in all the walls of the computa-
(rj + ∆r)(2rj − ∆r) tional domain (except for the symmetry axis). Those absorbing
Where the coefficients are obtained in order to ensure layers can be seen in Fig. 6 and 7 as a yellow zone. The
the Ez = 0 condition over the wall with the interpolated incoming waves into the Perfectly Matched Layer (PML)
 TABLE I
R ELATIVE ERRORS FOR THEH ERTZIAN DIPOLE SIMULATION .

Ez Er Bϕ
Radial direction 3.9% 12.7% 8.3%
Longitudinal direction 4.5% 12.8% 5.7%

show an attenuation until the fields completely drop to zero

as desired. The implementation of this boundary condition
was done based on [26] with a linearly increasing attenuation
parameter kf that depends on the distance from the PML
starting point, and it is zero outside that layer. This attenuation
parameter is treated as a non-advection term for equation 8,
and it is taken into account in the same way as the source
term S in eq. 15. For example, after the advection phase,
the following equation is applied to the right-propagating Lp
characteristic:

Ln+1 n+1
p(i,j) = Lp(i,j) (1.0 − ∆tkf ) (23)
−zo
Where kf = kfmax ( zzipml ). kfmax is the maximum value of
the attenuation parameter, zo is the starting point of the PML,
and zpml is the total length of the layer. The same procedure is
applied to the other outbound characteristics, adapting the kf
coefficient for each direction. One should also underline that Fig. 5. Computational domain for a) the Hertzian dipole, b) the coaxial cable.
the equations being in characteristic form, the implementation
of the PML is straightforward and does not require a modified
stencil.
A constant offset between the analytical results and the
Analytical CFL = 2.0 CFL = 1.5 CFL = 1.0
plots from the CIP simulations can be seen in Fig. 6 at 1
z = 0. An explanation to this offset comes from the fact that
the analytical expressions used for this comparison describe 0,5
Ez [V/m]

the radiation from the oscillation of two point charges in 0

free space separated by an infinitesimally small distance. The
numerical representation of this infinitesimal nature of the −0,5

Hertzian dipole is limited by the size of the mesh spacing. −1

Therefore, since a better match between the simulation and −1 −0,5 0 0,5 1
the theoretical results can be achieved with either a smaller 4
mesh spacing or a greater number of points per wavelength, 3
2
Bφ [nT]

the accuracy obtained with the previously described simulation

1
parameters is deemed sufficient for a proper comparison point. 0
The relative error along a given section compared to the −1
analytical result is given by: −2
−3
(F A − FiCIP )2 −1 −0,5 0 0,5 1
P
η(F ) = i Pi A 2 0,4
i (Fi )
0,2
Er [V/m]

Where F = Er , Ez , Bϕ , F A and F CIP are the analytical

and the simulation results, respectively. The results can be 0

seen in Table I. The overlapping of all the plots obtained with −0,2
the simulations using three different CFL values is a sign
of the robustness and precision of the multi-layer scattering −0,4
−1 −0,5 0 0,5 1
procedure even for larger CFL. Longitudinal axis [cm]

Fig. 6. Electromagnetic radiation of a Hertzian dipole in the longitudinal

B. Steady-State Open-Ended Coaxial Cable direction. Comparison among the exact results from an analytical model and
Another comparison was performed using the simulation different values of CFL (1.0, 1.5, and 2.0). Longitudinal electric field (top),
azimuthal magnetic field (middle), and radial electric field (bottom).
of an open-ended coaxial wave guide in vacuum, as shown
in Fig. 5. This test case has been chosen because it is
 Analytical CFL = 2.0 CFL = 1.5 CFL = 1.0 TABLE II
0,4 P OWER BUDGET FOR THE OPEN - ENDED COAXIAL CABLE .

0,2
Ez [V/m]

CIP Frequency Analytic

0
Domain [27]
Input power [W] 200.0 200.0 200.0
−0,2
Reflected power [W] 198.5 198.5 198.4
−0,4 Radiated power [W] 1.5 1.5 1.6
0 0,2 0,4 0,6 0,8 1
2
150
1 CIP COMSOL
Bφ [nT]

100

Er [kV/m]
0 50
0
−1
−50
−2
0 0,2 0,4 0,6 0,8 1 −100

0,4 −150
0 50 100 150 200
0,2 Longitudinal axis [mm]
Er [V/m]

100
0 CIP COMSOL
80

Er [kV/m]
−0,2 60

−0,4 40
0 0,2 0,8 1
Radial axis [cm] 20

Fig. 7. Electromagnetic radiation of a Hertzian dipole in the radial direction. 0

2 4 6 8 10 12
Comparison among the exact results from an analytical model and different Radial axis [mm]
values of CFL (1.0, 1.5, and 2.0). Longitudinal electric field (top), azimuthal
magnetic field (middle), and radial electric field (bottom). 0,6 CIP COMSOL
0,4
Bφ [mT]

0,2
0
representative of the coaxial configuration used in an ECR
−0,2
plasma thruster for which this model is developed [18], [20]. −0,4
The dimensions of the coaxial wave guide are: an inner radius −0,6
Rin = 1.15 mm, an outer radius of Rout = 13.75 mm, a 0 50 100 150 200
length L = 200 mm, and a rectangular computational domain Longitudinal axis [mm]
of z = 400 mm and r = 100 mm. A 2.45 GHz Transverse
Electro-Magnetic (TEM) input wave is imposed at the left- Fig. 8. Electromagnetic fields in a coaxial cable simulated using the CIP
hand wall. More precisely, the analytic solution for a right method and the COMSOL Multiphysics™ [22] software based on the finite-
propagating coaxial TEM mode was used to set the right- element method with a CFL = 1.0. The radial cross section is taken at
z = 80 mm, and the longitudinal cross section at the outer surface of the
propagating characteristics Lp and ∂z Lp in eq. 12 at the left inner conductor.
of the simulation domain boundary. The exiting left-going
characteristic Lm was left free. For the boundary conditions
we used a PML at the right-hand side and at the top right the reflected and radiated power are compared as well, as
corner, and a perfectly conducting surface boundary condition shown in Table II. In the limit of large wavelength compared
at the inner and outer radius. A constant mesh spacing of 0.23 to the coaxial diameter λ  Rout , an analytic formula for
mm was used. The simulation is run for 3800 iterations, which the radiated power exists as presented in Appendix B, and it
corresponds to more than two round trips of the incident wave, is included for comparison [27].
to establish a standing wave in the waveguide.
The standing wave pattern is retrieved for the CIP A good agreement is found between the two cases. The
simulation by taking the results once an amplitude peak relative error between the CIP radial electric field and the
in reached for the time-varying stationary wave. As a frequency domain result does not exceed 0.1%. The radi-
comparison, the electric field obtained using a Frequency ated power obtained with the CIP method also matches the
Domain electromagnetic solver from a commercial software is predicted value from the frequency-domain solver and the
compared to the wave pattern obtained with the CIP method. analytic formulation. This fact give us confidence about the
The same computational domain is used, using a triangulation stable nature of the solution obtained by our procedure.
of characteristic length of 0.2 mm. The electromagnetic Our time-domain solution successfully reproduced the steady
fields amplitudes are compared to the standing wave envelope state solution obtained with a frequency-domain finite-element
obtained from the CIP method, as shown in Fig. 8. In addition, method.
 IV. C ONCLUSION
k 4 (Rout
2 2 2
− Rin )
We have proposed a procedure for computing the propaga- Pout = Pin √  
Rout
tion of electromagnetic waves in axisymmetric 2D domains 12  ln Rin
using the Constrained Interpolation Profile scheme. The main
The formula is a function of the input power Pin , the inner
contribution of this development is the extension of the radial-
Rin and outer Rout radius, and the wave number k.
scattering algorithm for CFL conditions greater than 1.0 in
cylindrical configurations, as well as a simple treatment of
boundary conditions for the symmetry axis and conducting ACKNOWLEDGMENT
walls. For this work, we did not exceed a CFL condition of This work was made in the framework of the project
2.0, but in principle higher values may be achieved. In the MINOTOR that has received funding from the the European
future, our goal is to explore the CFL limit until which the Union’s Horizon 2020 research and innovation program under
procedure can be used based on a Von Neumann stability grant agreement No 730028.
analysis. In addition, the CIP scheme described above will
be applied to 2D cylindrical simulations of an ECR plasma R EFERENCES
source [20] where the source terms (charge density and plasma
[1] Shiraishi, K., and Matsuoka, T. “Wave Propagation Simulation Using
current) will be obtained using a Particle-In-cell code. In the CIP Method of Characteristic Equations.” Communications in Com-
this case where a strong coupling between the plasma and putational Physics, Vol. 3, No. 1, 2008, pp. 121–135.
the electromagnetic fields exists, the divergence conservation [2] Merino, M., Sanchez-Villar, A., Ahedo, E., Bonoli, P., Lee, J., Ram,
A., and Wright, J. ”Wave Propagation and Absorption in ECR Plasma
properties of the scheme are important and will be assessed. In Thrusters.” 35th International Electric Propulsion Conference, Atlanta,
particular, the need for a divergence cleaning method will be USA, 2017.
tested. If such were the case, one should note that the Method [3] T. Yabe and T. Aoki, A universal solver for hyperbolic equation by cubic-
polynomial interpolation, I. One-dimensional solver, Comput. Phys.
of Characteristic used in this approach enable a straightforward Commun., 66 (1991), 219-232.
extension to hyperbolic divergence cleaning methods [29]. [4] T. Yabe, F. Xiao and T. Utsumi, The constrained interpolation profile
method for multiphase analysis, J. Comput. Phys., 169 (2001), 556-593.
[5] P.K.Smolarkiewicz and J.A.Pudykiewicz,”A Class of Semi- Lagrangian
A PPENDIX A Approximations for Fluids,”J. Atmospheric Sciences, Vol.49, pp.2028-
E XPLICIT E XPRESSIONS FOR THE E LECTROMAGNETIC 2096, 1992.
[6] Nie, Y., Fu, K., and Lv, X. “CIP Method of Characteristics for the
F IELDS OF A H ERTZIAN D IPOLE Solution of Tide Wave Equations.” Advances in Mathematical Physics,
Vol. 2018, 2018. doi:10.1155/2018/3469534.
Explicit expressions for the electromagnetic fields of an [7] Tachioka, Y., Yasuda, Y., and Sakuma, T. “Application of the Con-
oscillating z-directed dipole p(t) = pẑ cos (ωt) were derived strained Interpolation Profile Method to Room Acoustic Problems:
and presented in [21] in spherical coordinates. Here r is the Examination of Boundary Modeling and Spatial/Time Discretization.”
Acoustical Science and Technology, Vol. 33, No. 1, 2012, pp. 21–32.
vector pointing to the observation point from the dipole and doi:10.1250/ast.33.21.
θ = ∠(r, z) is the angle with the dipole axis. [8] Kajita, K., Baba, Y., Nagaoka, N., and Ametani, A. “Computation of
Lightning Electromagnetic Pulses Using the Constrained Interpolation
   Profile Method.” Electric Power Systems Research, Vol. 115, 2014, pp.
sin (kr − ωt) sin θ 94–101. doi:10.1016/j.epsr.2014.02.017.
Hφ (r) = −k cos (kr − ωt) + pω [9] K. S. Yee, “Numerical solution of initial boundary value problems in-
r 4πr volving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas
  
cos (kr − ωt) 2 cos θ Propag., vol. AP-14, no. 4, pp. 302–307, May 1966.
Er (r) = k sin (kr − ωt) + p [10] K. Okubo and N. Takeuchi, “Analysis of an electromagnetic field created
r 4π0 r2 by line current using constrained interpolation profile method,” IEEE
  
cos (kr − ωt) sin θ Trans. Antennas Propag., vol.55, no.1, pp.111–119, Jan. 2007.
Eθ (r) = k sin (kr − ωt) + p [11] Tanaka Daiki. Stability Analysis of the CIP Scheme and its Applications
r 4π0 r2 in Fundamental Study of the Diffused Optical Tomography. PhD thesis,
pk 2 sin θ Kyoto University Graduate School of Informatics, 2014.
− cos (kr − ωt) [12] Ogata, Y., Yabe, T., and Odagaki, K. ”An Accurate Numerical Scheme
4π0 r for Maxwell Equations with CIP-Method of Characteristics”. Commu-
These expressions were transformed to a cylindrical coor- nications in Computational Physics. Vol 1. No. 2, 2006, pp. 311–335.
[13] Tanaka, Y., Baba, Y., Nagaoka, N., and Ametani, A. “Computation of
dinate system in order to perform the comparison with the Lightning Electromagnetic Pulses with the Constrained Interpolation
simulations results obtained using the CIP method. The next Profile Method in the 2-D Cylindrical Coordinate System.” IEEE Trans-
formulas are used for the radial and longitudinal electric field: actions on Electromagnetic Compatibility, Vol. 56, No. 6, 2014, pp.
1497–1505.
[14] Chakarothai, J., Chen, Q., and Sawaya, K. “Three-Dimensional Elec-
tromagnetic Scattering Analysis Using Constrained Interpolation Profile
Er (r, z) =Er (r) sin θ + Eθ (r) cos θ Method.” IEICE Transactions on Communications, Vol. E93–B, No. 10,
Ez (r, z) =Er (r) cos θ − Eθ (r) sin θ 2010, pp. 2619–2628.
[15] Kobayashi, S., Suzuki, Y., and Baba, Y. “Lightning Electromag-
netic Field Calculation Using the Constrained Interpolation Profile
A PPENDIX B Method with a Subgridding Technique.” IEEE Transactions on Elec-
tromagnetic Compatibility, Vol. 58, No. 5, 2016, pp. 1682–1685.
A NALYTICAL F ORMULA FOR THE R ADIATED P OWER FROM doi:10.1109/TEMC.2016.2575079.
AN O PEN -E NDED C OAXIAL C ABLE [16] Domı́nguez, Adrián & Cichocki, Filippo & Merino, Mario & Fajardo,
Pablo & Ahedo, Eduardo. (2018). Axisymmetric plasma plume charac-
In [27], we can find the following analytical expression for terization with 2D and 3D particle codes. Plasma Sources Science and
the radiated power Pout from an open-ended coaxial cable: Technology. 27. 10.1088/1361-6595/aae702.
[17] Szabo, James & Warner, Noah & Martinez-Sanchez, Manuel & Jean Porto received the B.S. degree in mechanical
Batishchev, Oleg. (2014). Full Particle-In-Cell Simulation Methodology engineering from the Universidad del Norte in Bar-
for Axisymmetric Hall Effect Thrusters. Journal of Propulsion and ranquilla, Colombia, and the M.S. degree specialized
Power. 30. 197-208. 10.2514/1.B34774. in aerospace engineering from the Arts et Métiers
[18] F. Cannat, T. Lafleur, J. Jarrige, P. Chabert, P.-Q. Elias, and D. Packan, Engineering School in Paris, France. He used to
“Optimization of a coaxial electron cyclotron resonance plasma thruster work as a simulation engineer using finite element
with an analytical model,” Phys. Plasmas, vol. 22, no. 5, p. 053503, analysis for the automotive industry. He is currently
2015. pursuing a Ph.D. degree in physics at the French
[19] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer National Aerospace Research Centre (ONERA) and
Simulation. Philadelphia, IOP Publishing, 1991. the Paris Observatory. His research interests include
[20] Jean C. Porto, Paul-Quentin Elias. ”Full-PIC Simulation of an ECR computational plasma physics for electric thrusters
Plasma Thruster with Magnetic Nozzle”. International Electric Propul- and numerical methods for partial differential equations.
sion Conference 2019, Sep. 2019, Vienne, Austria.
[21] Orfanidis, S. J. Chapter 15: “Radiation Fields”, in “Electromagnetic
Waves and Antennas.”, Vol. 2, No. Rutgers University, 2004, pp.
313–321. doi:10.1016/B978-075064947-6/50011-3.
[22] COMSOL Multiphysics® v. 5.4. www.comsol.com. COMSOL AB,
Stockholm, Sweden.
[23] Macnamara, S., and Strang, G. “Splitting Methods in Communication, Paul-Quentin Elias holds a M.S degree in applied
Imaging, Science, and Engineering.” 2016, pp. 1–21. doi:10.1007/978- physics from Ecole Centrale Paris, France (2003).
3-319-41589-5. He earned a PhD in Energetics in 2007, from the
[24] Okubo, K., Yoshida, Y., and Takeuchi, N. “Consideration of Treatment same institution. He is a Senior Research Scientist
of the Boundary between Different Media in Electromagnetic Field at the French National Aerospace Research Center
Analysis Using the Constrained Interpolation Profile Method.” IEEE (ONERA).
Transactions on Antennas and Propagation, Vol. 55, No. 2, 2007, pp.
485–489. doi:10.1109/TAP.2006.889984.
[25] Mohseni, K., and Colonius, T. “Numerical Treatment of Polar Coordi-
nate Singularities.” Journal of Computational Physics, Vol. 157, No. 2,
2000, pp. 787–795. doi:10.1006/jcph.1999.6382.
[26] Ishizuka, T., and Okubo, K. “Formulation and Examination of the Per-
fectly Matched Layer for Sound Field Analyses Using the Constrained
Interpolation Profile Method.” Acoustical Science and Technology, Vol.
34, No. 5, 2013, pp. 378–381. doi:10.1250/ast.34.378.
[27] Schelkunoff, S. A ”On Diffraction and Radiation of Electromagnetic
Waves” Physical Review, Vol. 56, pp. 308-316,1939.
[28] Ogata, Y., Yabe, T., Takizawa, K., and Ohkubo, T. “The Analysis of
Electromagnetic Waves Using CIP Scheme with Soroban Grid.” Compu-
tational Fluid Dynamics 2004, No. 3, 2006, pp. 141–146. doi:10.1007/3-
540-31801-1-16.
[29] C.-D. Munz, P. Ommes, R. Schneider, ”A three-dimensional finite-
volume solver for the Maxwell equations with divergence cleaning on
unstructured meshes”, Computer Physics Communications, Volume 130,
Issues 1–2,2000, pp. 83-117, doi:10.1016/S0010-4655(00)00045-X.
[30] Huang, W., Sloan, D. M. (1993). Pole condition for singular problems:
the pseudospectral approximation. Journal of Computational Physics,
107(2), 254-261.