CURVILINEAR VECTOR FINITE DIFFERENCE APPROACH TO THE COMPUTATION OF WAVEGUIDE MODES

We describe here a Vector Finite Difference approach to the evaluation of waveguide eigenvalues and modes for rectangular, circular and elliptical waveguides. The FD is applied using a 2D cartesian, polar and elliptical grid in the waveguide section. A suitable Taylor expansion of the vector mode function allows to take exactly into account the boundary condition. To prevent the raising of spurious modes, our FD approximation results in a constrained eigenvalue problem, that we solve using a decomposition method. This approach has been evaluated comparing our results to the analytical modes of rectangular and circular waveguide, and to known data for the elliptic case.


Introduction
In many applications, such as the analysis of waveguide junctions using mode matching [1], the knowledge of both eigenvalues and field distributions of waveguides modes [2] is required.The same type of information is also required in the analysis, using the method of moments (MoM), of thick-walled apertures [3].Indeed, these apertures can be considered as stub waveguides, and the modes of these guides are the natural basis functions for the MoM [4,5].
Apart from some simple geometries, mode computation cannot be done in closed form, so that suitable numerical techniques must be used.Until now, many different numerical techniques have been proposed.The most popular are the Finite Difference (FD) Method [6,7] and scalar and vector Finite Element Method (FEM).FD techniques, despite of their long history, are still very popular because of their simplicity and computational effectiveness.However, they compute the Hertz potential eigenfunctions, so that, since the basis functions of the MoM are the vector eigenfunctions, a numerical derivative is required, which can results in a reduced accuracy.
Aim of this work is to present the direct computation of mode vectors in a waveguide, using a finite difference (FD) approximation of the vector Helmholtz equation on a suitable discretization grid.Since we are mainly interested in using those modes in the MoM, the entire development will be expressed in term of equivalent magnetic surface currents.For each grid point, suitable second-order Taylor approximations allow to replace the continuous eigenfunction problem with a discrete one.This leads to a matrix eigenvector problem, when suitable conditions are added.These come out from the boundary conditions (which are included directly in the problem matrix), and the solenoidal or irrotational condition on mode vectors.This constrained eigenvalue problem can then be solved using linear algebra techniques [8].
We consider here rectangular, circular and elliptic [9] waveguides.In order to improve both the accuracy and the computational effectiveness, a discretization grid fitting exactly the waveguide boundary is chosen.Therefore, our FD approach has been implemented using cartesian, polar and elliptic frameworks.

Description of the Technique
Let us consider a waveguide.The TE mode vectors − → e are the eigenfunctions of the Helmholtz equation : where C is the contour of the waveguide section.If we introduce the (two-dimensional) magnetic current When substituted in (1) , after replacing and collecting terms we get: Figure 1: Vectors geometry with respect to the contour of the conductor.
The TE eigenvalue problem can therefore be rewritten, taking into account the contour geometry as in Fig. 1, as: In order to solve problem (5) with FD, both the unknown − → M and the equations of (5) are discretized, i.e., evaluated only on the points of a suitable grid.In this way, for each discretization point, we have a difference equation, corresponding to the first of (5), expressed in terms of the samples of − → M .This equation takes into account also the boundary condition on the waveguide contour.On the other hand, the requirement ∇ t × − → M = 0, suitably discretized, must be imposed separately.As a consequence, (5) becomes a constrained matrix eigenvalue problem.It is worth noting that the latter requirement could be incorporated into the eigenvalue equations.In this way a standard matrix eigenvalue problems follows, but it can contain spurious solutions, which are not present in our formulation.

Rectangular waveguides
Let us consider first a rectangular waveguide, whose discretization grid is shown in Fig. 2. Since this is the simpler structure, we discuss in detail, for this case, the passage from (5) to the discretized form.The discretization of the first of (5) depends on the grid point that we consider.For an internal point P , as in Fig. 3, we can use a second-order Taylor expansion as: We need also to discretize the second of (5).Using a first-order Taylor expansion (compare (6) ) we easily get: For a boundary point P , as in Fig. 4, we need a different approach, since D is not a sampling point for the current.We force in D the boundary condition M D,y = 0, and incorporate it into the FD matrix.Now M D,y can be obtained as a Taylor expansion: (11) which can be added to 1 2 M B,y given by ( 6), to get: and so for a boundary point P : For the x-component ∇ 2 t M P,x we need another grid point H (see Fig. 4), in which: From ( 14) and M B,x given by ( 6) we get: and the x-component of ( 9) is replaced by: In the same way, the condition (10) becomes: The FD matrix equivalent to ( 5) is then given by ( 9) (or (15) and ( 16) for a boundary point) with the constraint (10) or (17).
It is worth noting that the appproach described in this section can be used on more general waveguides, i.e., on all waveguides whose boundary is made of segments parallel to the cartesian axes.Let us consider now the discretization of (15) for a circular waveguide.In order to maintain the correct boundary condition, we use a polar discretization grid, as in Fig. 5, with spacing ∆r, ∆ϑ.
We start from the expansion of the vector Laplace operator in polar coordinates: ) For each internal point (whose geometry is quite similiar to the one of Fig. 3), using suitable Taylor approximations, as in the rectangular case, the derivatives in (18), in the sample point P , can be expresssed as: When substituted in (18) we get: For a boundary point, such as P in Fig. 6, the boundary condition is M D,r = 0 and this condition can be incorporated into the FD matrix, much in the same way as in the rectangular case.The result is that (19) is then replaced by: The constants of ( 19) and (20) are: The constraint: can be discretized in an internal point using a first-order Taylor expansion as: which becomes, in a boundary point: (23) It remains to consider the center of the circle.In this point it is not possible to use a Taylor expansion since it is a singular point for the polar framework.Therefore, we integrate the first of ( 5) on a circle S with radius ∆r/2 : Because of (18), the Laplace operator becomes : and substituting in (24): We apply the theorem of the gradient [10] to the l.h.s of (25) : From the divergence of − → M in polar coordinates we get, for the r.h.s of (26): The second integral of ( 28) and ( 29) can be evaluated by parts as: As a consequence, (25) becomes: − → M 0 is the value at the center of the circle, which must be expressed in Cartesian coordinates: To evaluate (32) we need M ϑ ( ∆r 2 , ϑ along C, and therefore outside the discretization grid.In order to include (32) into the FD matrix, we must express both in terms of − → M at the discretization points, one as an average, and the other as a finite difference.Since − → M 0 does not have polar components, M ϑ ( ∆r 2 , ϑ must be computed as (compare Fig. 7) : M r,q sin (ϑ q ) ∆ϑ (35) where M r,q = M r (∆r, q∆ϑ), and : [M ϑ,q + M 0,y cos ϑ q − M 0,x sin ϑ q ] cos ϑ q (37) where M ϑ,q = M ϑ (∆r, q∆ϑ).Summing (35) and (37) we get, for the x-component of the l.h.s. of (32): [2M r,q cos ϑ q + M ϑ,q sin ϑ q +M q 0,y sin ϑ q cos ϑ q − M q 0,x sin 2 ϑ q ] (38) In the same way, summing (34) and (36) we get, for the y-component of the l.h.s. of (32): [2M r,q sin ϑ q − M ϑ,q cos ϑ q −M q 0,y sin ϑ q cos ϑ q + M q 0,x cos 2 ϑ q ] (39)

Elliptic waveguides
In the same way as the circular case, in an elliptic waveguide, − → M is evaluated only on the points of an elliptic grid (see Fig. 8) with spacing ∆u, ∆v.

The expression of the Laplace vector operator in elliptic coordinates can be simplified if we let − →
is the common value of the scale factor, 2a f being the inter- and the v components is: Figure 9: Internal point of the elliptic cylindrical coordinates (grid T E).
For each internal grid point, as in Fig. 9, a second order Taylor approximation allows to discretize the Laplace operator (40,41) in terms of the current samples at the neighboring points.The resulting derivatives are as follows: The same simplification is obtained for the second equation ( 5) for the elliptic case, since letting − → A = h − → M makes this equation almost identical to the rectangular case: The singular points of the elliptical framework, either the foci, or the points on the inter-focal segment, require a different treatment.For a focus of the ellipse (Fig. 10), as in the circular case, we use the integral form of the eigenvalue equation, as in (25-27).The central term of (26) becomes The integrals are divided in 4 parts.We describe here in details only the evaluation of the part over ) and R = ( ∆u 2 , 0 ) , we have (45) and The same approach can be used for points on the interfocal segment, such as the one in Fig. 11 Figure 11: Point between foci.

Numerical Solution and Results
In all cases considered, we get a constrained eigenvalue problem as: where A is a (2n,2n) matrix, and C is (2n,m) with n > m and λ = −k 2 t .Of course, A is the discrete laplace operator, including the boundary condition, and C is the discrete form of the constrain ∇ t × − → M = 0. We can solve [11] (47) using the QR factorization: and R is a upper triangular matrix (2n,m).Substituting in second equation of (47) we get: which suggest to change the unknown as: so as (48) becomes: Replacing x = Q • y in the first of (47) and multiplying all members by Q T we obtain: where ) matrix.The unknown can be partitioned in two vector y = u v , where u and v are (n,1) vector.Then equations (50), (51) are rewritten in partitioned form as: where B ij are (n,n) matrices, and T 1 , is an (invertible) triangular matrix .The second equation of (52) becomes: Therefore the first of (53) becomes: and we need to extract the eigenvalues and eigenvectors of B 22 .This discretized eigenvalues problem must be solved by numerical routine and the full matrix routines of Matlab have been used.The waveguide modes can then be obtained, from the eigenvalues v of B 22 , as Using MATLAB routines, we have performed an extensive validation of the approach presented.Some results are shown in the following, while a wider collection of data can be found in [12].
Table 1 shows the validation test for a rectangular waveguide, whose size (in normalized units) is 3.35 × 1.65, discretized in 134 × 66 square discretization cells.Our data have been compared with the known analytical results [2], both for the eigenvalues and the mode vectors.The same comparison has been made in Table 2 for a circular waveguide, with a radius equal to 4 (in normalized units), and using discretization steps of ∆r = 0.0396, ∆ϑ = 1 o .In both tables k ta and k tpv are the eigenvalues computed analytically [2] and using our approach, e % is the percentage error between them, and the last column shows the RM S difference between analytical mode vectors, and those compured using (55).
From these data, it appears that our technique is able to give all mode data with a very small error, both on eigenvalues and on eigenvectors.The error on the former is no larger than 0.3 % for the rectangular guide and significantly smaller than this for the circular case.The RM S error on normalized modes is even smaller, being less than 5 • 10 −8 in both cases.
Table 3 shows the validation test for an elliptic waveguide, with a length of the minor axis (in normalized units) equal to 4, and different eccentricities e cx , discretized with ∆v = 1 o , and a different number N u of discretization ellipses.Though an analytical solution is available for elliptic waveguides [9], it effectiveness is very poor, and many numerical techniques has been proposed in the literature.As in many of those papers, we use, for our comparison, the eigenvalues k ta obtained from the cut-off wavelength data reported in [13].
From the data in Table 3, it is clear that, for all cases shown, the eigenvalues k tpv computed from our techniques are very accurate.

Conclusions
A Vector Finite Difference tecnique has been proposed for the computation of the modes and eigenvalues of rectangular, circular and elliptical waveguides.A suitable choice of the discretization grid, which exacly fits the waveguide boundary, allows to obtain a very effective and accurate procedure.The presented approach has been assessed against analytical modes of rectangular and circular waveguides, and data available in the open literature for elliptical waveguides.The agreement is excellent in all the cases.

Figure 2 :
Figure 2: TE grid for the rectangular case.

Figure 3 :
Figure 3: Internal Point of TE grid.

Figure 5 :
Figure 5: TE grid for a circular waveguide.

Figure 6 :
Figure 6: Boundary point of polar framework.

Figure 7 :
Figure 7: component of M between center and next point.

Figure 8 :
Figure 8: Geometry of the elliptic cylindrical coordinates.

Figure 10 :
Figure 10: Focus A of the ellipse.

Table 1 :
Validation test for a rectangular waveguide.

Table 2 :
Validation test for a circular waveguide.