Formulation of the problem

To calculate an inhomogeneously filled circular shielded waveguide, it is proposed to use a modified Galerkin method [Reference Alsuyuti, Doha, Ezz-Eldien, Bayoumi and Baleanu22, Reference Duvigneau23], which is a variant of the general spectral method.

Let us consider the problem of the propagation of symmetric E and H-waves in a circular shielded waveguide with partial dielectric filling, the value of the dielectric constant of which arbitrarily depends on the radial coordinates ɛ(r, z, φ) = ɛ(r) (Fig. 1). The value of the magnetic permeability is assumed to be constant.

Fig. 1. The distribution function of the dielectric constant in the cross-section of the waveguide.

From Maxwell's equations we get:

(1)$$rotrotE = k_0^2 \varepsilon ( r, \;\varphi ) E.$$

Using the following expressions

$$grad( \psi ) = r\displaystyle{{d\psi } \over {dr}} + \varphi \displaystyle{1 \over r}\displaystyle{{d\psi } \over {d\varphi }} + z\displaystyle{{d\psi } \over {dz}}, \;$$

$$div( E ) = \displaystyle{{\partial E_r} \over {\partial r}} + \displaystyle{1 \over r}\displaystyle{{\partial E_\varphi } \over {\partial \varphi }} + \displaystyle{{\partial E_z} \over {\partial z}}, \;$$

$$rot( E) = r\left({\displaystyle{1 \over r}\displaystyle{{\partial E_z} \over {\partial \varphi }}-\displaystyle{{\partial E_\varphi } \over {\partial z}}} \right) + \varphi \left({\displaystyle{{\partial E_r} \over {\partial z}}-\displaystyle{{\partial E_z} \over {\partial r}}} \right) + z\left({\displaystyle{{\partial E_\varphi } \over {\partial r}} + \displaystyle{1 \over r}E_\varphi -\displaystyle{1 \over r}\displaystyle{{\partial E_r} \over {\partial \varphi }}} \right), \;$$

we write equation (1) for the field components in a cylindrical coordinate system:

(1a)$$\eqalign{rot {rot( E) } \vert _r & = \displaystyle{1 \over r}\displaystyle{{\partial ^2E_\varphi } \over {\partial r\partial \varphi }} + \displaystyle{{\partial ^2E_z} \over {\partial r\partial z}}-\displaystyle{1 \over {r^2}}\displaystyle{{\partial ^2E_r} \over {\partial \varphi ^2}} + \displaystyle{1 \over {r^2}}\displaystyle{{\partial E_\varphi } \over {\partial \varphi }}-\displaystyle{{\partial ^2E_r} \over {\partial z^2}} \cr & = k_0^2 \varepsilon ( r, \;\varphi ) E_r, \;}$$

(1b)$$\eqalign{rot {rot( E) } \vert _\varphi & = \displaystyle{1 \over {r^2}} E_\varphi + \displaystyle{1 \over r}\displaystyle{{\partial ^2E_r} \over {\partial r\partial \varphi }}-\displaystyle{1 \over {r^2}}\displaystyle{{\partial E_r} \over {\partial \varphi }} + \displaystyle{1 \over r}\displaystyle{{\partial ^2E_z} \over {\partial \varphi \partial z}}-\displaystyle{{\partial ^2E_\varphi } \over {\partial r^2}} \cr & -\displaystyle{1 \over {\partial \varphi}} \displaystyle{{\partial E_\varphi } \over {\partial r}}-\displaystyle{{\partial ^2E_\varphi } \over {\partial z^2}} = k_0^2 \varepsilon ( r, \;\varphi ) E_\varphi,}$$

(1c)$$\eqalign{rot {rot( E) } \vert _z & = \displaystyle{{\partial ^2E_r} \over {\partial r\partial z}} + \displaystyle{1 \over r}\displaystyle{{\partial E_r} \over {\partial z}} + \displaystyle{1 \over r}\displaystyle{{\partial ^2E_\varphi } \over {\partial \varphi \partial z}}-\displaystyle{{\partial ^2E_z} \over {\partial r^2}}-\displaystyle{1 \over r}\displaystyle{{\partial E_z} \over {\partial r}}-\displaystyle{1 \over {r^2}}\displaystyle{{\partial ^2E_z} \over {\partial \varphi ^2}} \cr & = k_0^2 \varepsilon ( r, \;\varphi ) E_z.}$$

We represent the wave fields of the guiding structure in the form of expansions in terms of eigenfunctions of the Dirichlet and Neumann boundary value problems for a uniformly filled circular waveguide. The connection between the components of the electric field, in accordance with the spectral method, is established through the coefficients of the series of expansions substituted in (1).

Symmetrical H-waves

In the absence of the angular dependence of the field, we assume ∂/∂φ = 0, E _r = 0, E _z = 0. In this case, equation (1) will be reduced to a single equation for the φ component of the electric field

$$\displaystyle{1 \over {r^2}}E_\varphi -\displaystyle{{\partial ^2E_\varphi } \over {\partial r^2}}-\displaystyle{1 \over r}\displaystyle{{\partial E_\varphi } \over {\partial r}}-\displaystyle{{\partial ^2E_\varphi } \over {\partial z^2}} = k_0^2 \varepsilon ( r ) E_\varphi.$$

Writing E _φ(r, φ, z) = E _φ(r, φ)e ^−iβz we obtain an equation for the transverse coordinate function

(2)$$\displaystyle{{\partial ^2E_\varphi } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial E_\varphi } \over {\partial r}}-\displaystyle{1 \over {r^2}}E_\varphi + ( {k_0^2 \varepsilon ( r ) -\beta^2} ) E_\varphi = 0, \;$$

where $\varepsilon ( r) = \left\{{\matrix{ {\varepsilon_1-{{\varepsilon_1-\varepsilon_2} \over {a^2}}r^2, \;r \le a} \hfill \cr {1, \;a \le r \le R.} \hfill \cr } } \right., \;$

ɛ − relative dielectric constant.

Assuming the dependence of the field on the longitudinal coordinate and time, we obtain equations for the components of the electric field:

(2a)$$\left\{{\matrix{ {\displaystyle{{\partial^2E_\varphi } \over {\partial a^2}} + ( {k_0^2 \varepsilon ( r, \;a) -\beta^2} ) E_\varphi -\displaystyle{{\partial^2E_a} \over {\partial r\partial a}} + i\beta \displaystyle{{\partial E_\theta } \over {\partial r}} = 0, \;} \hfill \cr {\displaystyle{{\partial^2E_a} \over {\partial r^2}} + ( {k_0^2 \varepsilon ( r, \;a) -\beta^2} ) E_a-\displaystyle{{\partial^2E_r} \over {\partial r\partial a}} + i\beta \displaystyle{{\partial E_\theta } \over {\partial a}} = 0, \;} \hfill \cr {\displaystyle{{\partial^2E_\theta } \over {\partial r^2}} + \displaystyle{{\partial^2E_\theta } \over {\partial a^2}} + k_0^2 \varepsilon ( r, \;a) E_\theta + i\beta \displaystyle{{\partial E_r} \over {\partial r}} + i\beta \displaystyle{{\partial E_a} \over {\partial a}} = 0.} \hfill \cr } } \right.$$

The solution to equation (2a) will be sought [Reference Duvigneau23] in the form:

(3)$$E_\varphi ( r) = \sum\limits_{n = 0}^N {b_nJ_1} ( \alpha _nr) , \;$$

where J ₁(α _nr) is the Bessel function of the 1st order, the coefficients α _n are determined taking into account the boundary condition E _φ(r = R) from equation J ₁(α _nR) = 0.

Substituting (3) into (2), we obtain

$$\eqalign{-\sum\limits_{n = 0}^N {b_n\left[{\displaystyle{{\partial^2J_1( \alpha_nr) } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial J_1( \alpha_nr) } \over {\partial r}}-\displaystyle{1 \over {r^2}}J_1( \alpha_nr) } \right]} = \cr = \sum\limits_{n = 0}^N {b_nk_0^2 \varepsilon ( r) J_1( \alpha _nr) -} \sum\limits_{n = 0}^N {b_n\beta ^2J_1( \alpha _nr).} } $$

Considering that

$$\displaystyle{{\partial ^2J_1( \alpha _nr) } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial J_1( \alpha _nr) } \over {\partial r}}-\displaystyle{1 \over {r^2}}J_1( \alpha _nr) = {-}\alpha _n^2 J_1( \alpha _nr) , \;$$

we get

(4)$$\sum\limits_{n = 0}^N {b_n( {\alpha_n^2 + \beta^2} ) J_1( \alpha _nr) } = \sum\limits_{n = 0}^N {b_nk_0^2 \varepsilon ( r) J_1} ( \alpha _nr).$$

Multiplying both sides of equation (4) by rJ ₁(α _qr) and integrating within r ∈ [0;R], we obtain the equation

(5)$$( {\alpha_q^2 + \beta^2} ) Q_qb_q = \sum\limits_{n = 0}^N {b_n} k_0^2 \int_0^R {\varepsilon ( r) r} J_1( \alpha _nr) J_1( \alpha _qr) dr.$$

Here we used the orthogonality condition for the Bessel functions:

$$\int_0^R r J_1( \alpha _nr) J_1( \alpha _qr) dr = \left\{{\matrix{ {Q_q, \;q = n} \hfill \cr {0, \;q\ne n} \hfill \cr } , \;} \right.$$

where $Q_q = 0.5R^2J_0^2 ( \alpha _nR) , \;$ which takes place, since in this case the Bessel functions are a solution to the homogeneous boundary value problem on the Bessel equation.

Equation (5) can be represented in matrix form:

(6)$$M\cdot b = T\cdot b, \;$$

where

$$M_{q, n} = \left\{\matrix{( \alpha_q^2 + \beta^2) Q_q, \;q = n, \;\hfill \cr 0, \;q\ne n, \;\hfill} \right.$$

$$T_{q, n} = k_0^2 \int_0^R \varepsilon ( r) rJ_1( \alpha _nr) J_1( \alpha _qr) dr.$$

Writing equation (6) in the form (M − T) ⋅ b = 0 and equating the determinant of matrix (M − T) to zero, we obtain the dispersion equation for symmetric H-waves propagating in a circular waveguide with an arbitrary dependence of ɛ on r:

(7)$$Det( \beta ) = M-T = 0.$$

Note that the matrix T does not depend on β, therefore, when solving the dispersion equation (7), it is calculated only once, which significantly reduces the search time for the roots of the dispersion equation. Note that, when deriving equations (6) and (7), no restrictions were imposed on the form of dependence ɛ(r), i.e. this method allows one to calculate symmetric H-waves with a completely arbitrary nature of the change in the dielectric constant along the transverse coordinate, while ɛ can also be a complex quantity, which allows, for example, calculating waveguides with a complex absorption distribution in the cross-section, that is, to solve non-self-adjoint boundary value problems, in which the identity of the differential operators of the direct and adjoint boundary value problems is not satisfied.

Symmetrical E-waves

For symmetric E-waves, we put

$$\displaystyle{\partial \over {\partial \varphi }} = 0, \;E_\varphi = 0, \;H_r = H_z = 0.$$

In this case, equation (1) transforms into a system of two equations:

$$\displaylines{i\beta \displaystyle{{\partial E_z} \over {\partial r}} + ( {k_0^2 \varepsilon ( r) -\beta^2} ) E_r = 0, \;\cr i\beta \displaystyle{{\partial E_r} \over {\partial r}} + i\beta \displaystyle{1 \over r}E_r + \displaystyle{{\partial ^2E_\varphi } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial E_z} \over {\partial r}} + k_0^2 \varepsilon ( r) E_z = 0.} $$

Introducing variable $\widetilde{{E_z}} = i\beta \cdot E_z, \;$ we arrive at the equations:

(8)$$\eqalign{\displaystyle{{\partial \widetilde{{E_z}}} \over {\partial r}} & + ( {k_0^2 \varepsilon ( r) -\beta^2} ) E_r = 0, \;\cr \displaystyle{{\partial ^2\widetilde{{E_z}}} \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial \widetilde{{E_z}}} \over {\partial _r}} & + k_0^2 \varepsilon ( r) \widetilde{{E_z}}-\beta ^2\left({\displaystyle{{\partial E_r} \over {\partial_r}} + \displaystyle{1 \over r}E_r} \right) = 0.} $$

The boundary conditions on an ideally conducting surface for the tangential and normal components of the electric field $E_\tau \left\vert {_s } = 0 {, \;{{\partial E_n} \over {\partial_n}}} \right\vert _n = 0$ [Reference Islamov and Ismibayli18], in this case lead to the equation

(9)$$E_z\vert {_{r = R} } = 0.$$

The components of the electric field in accordance with the spectral method will be sought in the form of autonomous expansions:

(10)$$\widetilde{{E_z}} = \sum\limits_{n = 0}^N {A_n} J_0( \alpha _nr) , \;E_r = \sum\limits_{m = 0}^N {B_mJ_1( \alpha _mr).} $$

Taking into account the first boundary condition (9), the wave numbers α _n are determined from equation J ₀(α _nR) = 0.

Substituting (10) into (8), we obtain a system of two functional equations:

(11)$$\eqalign{ -\sum\limits_{n = 0}^N {A_n} \alpha _nJ_1( \alpha _nr) + \sum\limits_{m = 0}^N {B_m( {k_0^2 \varepsilon ( r) -\beta^2} ) J_1( {\alpha_mr} ) = 0, \;} \cr \sum\limits_{n = 0}^N {A_n} \left({\displaystyle{{\partial^2J_0( {\alpha_nr} ) } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial J_0( {\alpha_nr} ) } \over {\partial r}} + k_0^2 \varepsilon ( r) J_0( {\alpha_nr} ) } \right)- \cr \quad -\beta ^2\sum\limits_{m = 0}^N {B_m\left({\displaystyle{{\partial J_1( {\alpha_mr} ) } \over {\partial r}} + \displaystyle{1 \over r}J_1( {\alpha_mr} ) } \right) = 0.} } $$

Taking into account the equalities

$$\displaystyle{{\partial ^2J_0( {\alpha_nr} ) } \over {\partial r^2}} + \displaystyle{1 \over r}\displaystyle{{\partial J_0( {\alpha_nr} ) } \over {\partial r}} = {-}\alpha _0^2 J( r) , \;$$

$$\displaystyle{{\partial J_1( {\alpha_mr} ) } \over {\partial r}} = \alpha _mJ_0( \alpha _mr) -\displaystyle{1 \over r}J_1( \alpha _mr) , \;$$

system (11) can be rewritten as

(12a)$$-\sum\limits_{n = 0}^N {A_n\alpha _nJ_1( \alpha _nr) + } \sum\limits_{m = 0}^N {B_m( {k_0^2 \varepsilon ( r) -\beta^2} ) J_1( \alpha _mr) = 0, \;} $$

(12b)$$-\sum\limits_{n = 0}^N {A_n( {k_0^2 \varepsilon ( r) -\alpha^2} ) J_0( \alpha _nr) -\beta ^2\sum\limits_{m = 0}^N {B_m\alpha _mJ_0( \alpha _mr) } = 0.} $$

Multiplying equation (12a) by rJ ₁(α _qr) = 0, equation (12b) by rJ ₀(α _qr) = 0 and integrating within r ∈ [0;R], we obtain the system of equations:

(13)$$\eqalign{& -A_q\alpha _qQ_q + k_0^2 \sum\limits_{m = 0}^N {B_m} \int_0^R r \varepsilon ( r) J_1( \alpha _mr) J_1( \alpha _qr) dr- \cr & -B_q\beta ^2Q_q = 0, \;\cr& \quad \quad k_0^2 \sum\limits_{n = 0}^N {A_n} \int_0^R r \varepsilon ( r) J_0( \alpha _nr) J_0( \alpha _qr) dr- \cr & \!\!\! \quad \quad -A_q\alpha _q^2 Q_q-B_q\beta ^2\alpha _qQ_q = 0.} $$

Here we used the orthogonality conditions for the Bessel functions

$$\eqalign{\int_0^R r J_0( \alpha _nr) J_0( \alpha _qr) dr & = \int_0^R r J_1( \alpha _nr) J_1( \alpha _qr) dr = \cr & = \left\{{\matrix{ {Q_q, \;Q = n} \hfill \cr {0, \;q\ne n} \hfill \cr } , \;} \right.} $$

where $Q_{_q } = ( {R^2/2} ) J_1^2 ( {\alpha_qR} )$, since in this case the Bessel functions are a solution to a homogeneous boundary value problem.

The system of equation (13) can be written in matrix form:

(14)$$\left[{\matrix{ {T^{( 0, 0) }T^{( 0, 1) }} \cr {T^{( 1, 0) }T^{( 1, 1) }} \cr } } \right]\cdot \left[{\matrix{ A \cr B \cr } } \right] = 0, \;$$

where

(15)$$\eqalign{& T_{q, m}^{( 0, 0) } = {-}\alpha _qQ_q\delta _{q, m}, \;\cr T_{q, m}^{( 0, 1) } = k_0^2 \int_0^R r & \varepsilon ( r) J_1( {\alpha_mr} ) J_1( {\alpha_qr} ) dr-\beta ^2Q_q\delta _{q, m}, \;\cr T_{q, m}^{( 1, 0) } = k_0^2 \int_0^R r & \varepsilon ( r) J_0( {\alpha_nr} ) J_0( {\alpha_qr} ) dr-\alpha _q^2 Q_q\delta _{q, m}, \;\cr & T_{q, m}^{( 1, 1) } = {-}\beta ^2\alpha _qQ_q\delta _{q, m}, \;} $$

δ _q,n −Kronecker symbol.

Equating the determinant of matrix equation (14) to zero, we obtain a dispersion equation describing the symmetric E-waves of a circular waveguide with an arbitrary radial dielectric filling.

Numerical implementation of algorithms

Two-layer shielded waveguide. As an example, we use equations (7) and (14) to calculate the simplest test structure – a circular waveguide with a homogeneous dielectric rod (i.e. ɛ(r) = ɛ = const, Fig. 2) and compare the results with the exact ones obtained by the classical method of partial regions.

Fig. 2. Circular waveguide with a dielectric rod.

The calculations were carried out for a waveguide with parameters: R = 20mm,  a = 10mm,  ɛ = 3, at a frequency of f = 10GHz.

The classical calculation method gives the following results: for symmetric H-waves β _H = 2376891/m, for symmetric E-waves β _E = 22755000/m.

The calculation of test structures using the proposed technique was carried out by substituting the function $\varepsilon ( r) = \left\{{\matrix{ {3, \;r \le a} \hfill \cr {1, \;a \le r \le R} \hfill \cr } } \right.$ into equations (6) and (14).

The convergence of solutions obtained by the modified Galerkin method for symmetric E and H-waves is shown in Table 1 and in Fig. 3.

Fig. 3. Convergence in integral characteristics.

Table 1. Calculation by the modified Galerkin method

From Table 1 and Fig. 3 it follows that the convergence of the modified Galerkin method is monotonic and occurs rather quickly (already at N = 5 the difference between longitudinal wave numbers does not exceed 1.5%).

Figure 3 also shows that in the case of symmetric H-waves, convergence occurs faster, which, apparently, is associated with the difference in the number of equations to be solved (one equation (2) for symmetric H-waves and two equation (8) for symmetric E-waves).

In Fig. 4, the dotted line shows the dependences of the field components H _z and E _φ on the coordinate r, calculated for the symmetric H-wave at N = 5.

Fig. 4. Field distribution of the first symmetric H-wave: dotted line – partial domain method, solid line – modified Galerkin method.

From the graphs shown in Fig. 4 that the field distributions calculated by two different methods practically coincide.

Thus, using the example of a test problem with an exact solution, a high accuracy, efficiency of the method, and fast convergence of the solution obtained using the modified Galerkin method are shown.

Calculation of a waveguide with a rod, the dielectric constant of which changes according to the parabolic law

Based on equation (15), the dispersion characteristics of symmetric E-waves propagating in a circular waveguide with partial dielectric filling, the permeability of which changes according to the parabolic law, are described by the equation:

$$\varepsilon ( r) = \left\{{\matrix{ {\varepsilon_1-\displaystyle{{\varepsilon_1-\varepsilon_2} \over {a^2}}r^2, \;r \le a} \hfill \cr {1, \;a \le r \le R.} \hfill \cr } } \right.$$

Substituting this expression in (15) and calculating the integrals (numerically or analytically), we obtain a solution to the dispersion problem. Note that for any calculation of the integrals from (15) is carried out only once, since they do not depend on either the frequency or the longitudinal wavenumber, and are determined only by the filling parameters. This is an unconditional advantage of this method, which makes it possible to significantly reduce the time for calculating the characteristics of the structure.

The results of calculating the dispersion characteristics of symmetric E-waves of a circular waveguide with a parabolic profile of the dielectric filling are shown in Fig. 5. Figure 6 shows the distribution of the Umov-Poynting vector over the cross-section of the waveguide, calculated for three modes at frequency f = 14 GHz (points 1, 2, 3 in Fig. 5).

Fig. 5. Dispersion characteristics of symmetric E-waves of a circular waveguide with a parabolic profile of dielectric filling.

Fig. 6. Distribution of power flux density symmetric E-waves at frequency f = 14 GHz.

Calculation of a waveguide with a rod, the dielectric constant of which varies linearly

Based on equation (7), the structure is calculated in the form of a circular waveguide with partial dielectric filling, the permeability of which varies linearly (Fig. 7) within r ∈ [0/a]. The calculations were carried out for a waveguide with parameters R = 20 mm,  a = 10 mm, ɛ(r) = ɛ ₁ − ((ɛ ₁ − ɛ ₂)/a)r, ɛ ₁ = 6, ɛ ₂ = 2 frequency f = 10 GHz.

Fig. 7. Dielectric constant function.

For comparison, the calculation of the same structure was performed with the representation of the linear profile of the dielectric constant in the form of a step approximation (Fig. 7) with the number of steps equal to 20. The results of the calculation of the field distribution obtained by solving the dispersion equation (7) are shown in Fig. 8. The results of calculating the field distribution, performed according to the proposed technique and using the partial domain method, coincide with the graphic accuracy.

Fig. 8. Results of calculating the distribution of the wave field H ₀₁ of a waveguide with a linear distribution of permeability.