**5. Numerical analysis of multimodal waveguide: Model of confined environments**

So far, we have considered the modeling of guiding structures with metallic walls. However, tunnels have dielectric walls and other materials within them (trains, cars, objects, etc.). They also involve boundary conditions which can be time or frequency dependent. The limitation of the computational domain within dielectric materials is usually achieved by using perfectly matched layer (PML) for transmission-line method (TLM), finite-difference time domain (FDTD), or finite element method (FEM) implementations. However, this solution is too expensive in terms of complexity and consequently, computation cost when dealing with very large electrical objects. Since we are interested only in fields within the hollow region of the tunnel, appropriate boundary conditions, namely, the surface impedance boundary conditions (SIBCs), constitute the mathematical artifice to limit the computational domain.

### **5.1. Assumptions for tunnel wall modeling: SIBC concept**

To apply the SIBC concept, fields in the walls region are assumed to be time exponential decaying, source free, and with no reflections returning to the interface air and tunnels walls. For practical purposes, the wall thickness *t* should be larger than the skin depth *δ*, (*t* >> *δ*). A thickness of 2*δ* would ensure an amplitude of 1.8% of the launched wave [30]. Another important requirement is the smallest radius of curvature on a surface *Rmin*, which must be larger than the skin depth (*Rmin* >> *δ*). Consequently, the interface can be considered locally flat and plane-wave concept can be applied. Considering a numerical method like TLM, using the common criteria of *λ*/10 for the maximum mesh size, waves can be treated locally as plane waves. Moreover, at this resolution level, the staircase mesh is a good approximation of curved surfaces.

### **5.2. SIBC implementation in time domain for tunnel wall modeling**

The implementation of the SIBC concept was first introduced for good conductors. An overview of the implementations of these techniques in time domain can be found in [5, 23, 27]. Tunnel walls are typically modeled by a lossy dielectric with relative permittivity *<sup>r</sup>* <sup>=</sup> <sup>10</sup> <sup>−</sup> 12 and conductivity *<sup>σ</sup>* <sup>=</sup> 0.01 <sup>−</sup> 0.02*Sm*−1[9]. Maloney in [21] first introduced the SIBC for lossy dielectric in FDTD by using an expansion of the impedance model proposed by [5] for good conductors. Finally, [27] proposed a way to improve its implement in FDTD.

In contrast with good conductors, the SIBC for lossy dielectric varies with the angle of incidence and polarization of incident waves. If one considers Figure 10, surface impedances for vertical and horizontal polarizations are given by [24]

$$Z\_{Y}\left(\mathbf{s}\right) = \frac{E\_{\mathbf{x}}^{'+}\left(\mathbf{x}, \mathbf{0}, \mathbf{s}\right)}{H\_{\mathbf{y}}^{'+}\left(\mathbf{x}, \mathbf{0}, \mathbf{s}\right)} = \eta \frac{\cos\left(\theta\_{l}\right)}{\mathbf{s} + \sigma/\epsilon} \sqrt{s^{2} + 2as} \tag{62}$$

$$Z\_H(\mathbf{s}) = -\frac{E\_y^{'+}\left(\mathbf{x}, \mathbf{0}, \mathbf{s}\right)}{H\_{\mathbf{x}}^{'+}\left(\mathbf{x}, \mathbf{0}, \mathbf{s}\right)} = \eta \frac{\mathbf{s}}{\cos\left(\theta\_l\right)\sqrt{s^2 + 2as}}\tag{63}$$

where

illustrate this point, let us consider the calculation of the propagated error on *α*. The total differential of a function *z* = *f* (*x*1, *x*2, ..., *xN*), where *xi* is the *i*-th independent variable, is

the term *dz* denotes the error in *z*. The uncertainty in the calculation of the attenuation constant for waveguides can be found by considering the differential of the attenuation and phase constants, *α* and *β*, respectively. The differential of the attenuation and phase constants

> *∂α ∂ω*

> *∂β ∂ω*

Finally, by dividing (59) by (60), the differential of the attenuation constant in terms of the

As frequency increases, the first partial derivative on the right-hand side of (61) tends to zero and the second to the speed of light in free space, *c*0. Hence, for high frequencies, the error ∆*α* is bounded. For low frequencies, the first derivate tends to infinity and the second to a finite value, resulting in a larger error. Indeed, in Figure 9, the highest difference between both approaches occurs at cutoff frequencies of the modes and decreases far from these frequencies. Apart from this relative error on *α*, the comparison between numerical results and theory for both propagation constant *β* and the attenuation constant *α* confirms the suitability of the 2.5D TLM to characterize modal parameters in waveguiding structures.

 *∂ω ∂β* 

 *∂α ∂ω*

**5. Numerical analysis of multimodal waveguide: Model of confined**

So far, we have considered the modeling of guiding structures with metallic walls. However, tunnels have dielectric walls and other materials within them (trains, cars, objects, etc.). They also involve boundary conditions which can be time or frequency dependent. The limitation of the computational domain within dielectric materials is usually achieved by using perfectly matched layer (PML) for transmission-line method (TLM), finite-difference time domain (FDTD), or finite element method (FEM) implementations. However, this solution is too expensive in terms of complexity and consequently, computation cost when dealing with very large electrical objects. Since we are interested only in fields within the hollow region of the tunnel, appropriate boundary conditions, namely, the surface

∆*x*<sup>2</sup> + ... +

 *∂ f ∂xN* ∆*ω*, (59)

∆*ω*. (60)

∆*β* (61)

∆*xN* (58)

 *∂ f ∂x*<sup>2</sup> 

∆*α* =

∆*β* =

∆*α* =

given by

can be computed by

258 Advanced Electromagnetic Waves

is obtained.

**environments**

*dz* =

differential of the phase constant is given by

 *∂ f ∂x*<sup>1</sup> 

∆*x*<sup>1</sup> +

$$\mathfrak{a}\_{\parallel} = \frac{\sigma}{2\varepsilon \cos^2(\theta\_l)} = \frac{\sigma}{2\varepsilon\_0 \varepsilon\_r \cos^2(\theta\_l)}$$

$$\cos\left(\theta\_{l}\right) = \sqrt{1 - \frac{v^{2}}{c^{2}}\sin^{2}\left(\theta\_{l}\right)} = \sqrt{1 - \frac{\sin^{2}\left(\theta\_{l}\right)}{\mu\_{r}\epsilon\_{r}}}\tag{64}$$

$$\upsilon \quad = \qquad \qquad \frac{1}{\sqrt{\mu\epsilon}} = \frac{1}{\sqrt{\mu\_{0}\mu\_{r}\epsilon\_{0}\epsilon\_{r}}}$$

**Figure 10.** Geometry and coordinate system for a TE and TM polarized waves incident onto a dispersive half space

*Z* (*s*) is the surface impedance; *s* is the complex frequency; *η* is the intrinsic impedance of the lossy dielectric; , *µ*, and *σ* its constitutive parameters; and *θ<sup>i</sup>* and *θ<sup>t</sup>* are the angles of incidence and refraction, respectively. In tunnel environments, the incident angle *θ<sup>i</sup>* for the waves is unknown. Moreover, multiple reflections of waves with different angles of incidence are involved. To the best of the authors' knowledge, no work has been reported for SIBC under these conditions. A procedure to implement it in any time-domain method as well as its validation is presented in [2, 4]. Here, a brief explanation of its implementation is described.

### **5.3. SIBC formulation for tunnel wall modeling**

The SIBC is inherently a frequency-domain concept. When it is transformed into the time domain for use in methods such as TLM, it is replaced by a convolution integral. Maloney in [21] proposed a recursive formula to overcome the computational difficulties for its calculation. To implement the SIBC in TLM, the reflection coefficient of (4) has to be calculated. This coefficient depends on the polarization, frequency, and angle of incidence. Thus, their influences have to be studied.

For tunnel walls, the relative permittivity *<sup>r</sup>* >> 1. Then, for these values in (64), the term *cos*(*θt*) in (Eqs.62, 63) tends to one, obtaining the two angle-independent expressions of the surface impedance for both polarizations:

$$Z\_V\left(\mathbf{s}\right) = Z\_{||pol-TM} = \frac{\eta}{\mathbf{s} + \sigma/\epsilon} \sqrt{\mathbf{s}^2 + 2\mathbf{a}\mathbf{s}}\tag{65}$$

$$Z\_H\left(\mathbf{s}\right) = Z\_{\perp pol-TE} = \eta \frac{\mathbf{s}}{\sqrt{\mathbf{s}^2 + 2\mathbf{a}\mathbf{s}}}\tag{66}$$

By using a *rational polynomial approximation* [12], the resulting coefficient of (4) can be calculated for both polarizations into the form

$$P\left(\mathbf{s}\right) \approx \tilde{P}\left(\mathbf{s}\right) = \frac{N\left(\mathbf{s}\right)}{D\left(\mathbf{s}\right)} = \Pi\_{i=0}^{M} \left[\frac{\mathbf{s} - \tilde{\xi}\_{i}}{\mathbf{s} - \tilde{\zeta}\_{i}}\right] \approx \frac{V'\left(\mathbf{s}\right)}{V^{i}\left(\mathbf{s}\right)'}\tag{67}$$

where *ξ<sup>i</sup>* and *ζ<sup>i</sup>* correspond to the poles and zeros of *P*˜ (*s*) and *si* are the sample points. Consider discrete values of a function *P* (*S*) at points *Si*. One can find two polynomials *N* (*s*) and *D* (*s*) of order *M*, such that *P*˜ (*s*) = *N* (*s*) /*D* (*s*) best approximates *P* (*s*), in the least square sense.

The next step is to obtain (67) in *z*-domain. The *bilinear z-transform* can be employed, obtaining the rational polynomial expression

$$P\left(z\right) = \Pi\_{i=0}^{M} \left[\frac{z - \tilde{\xi}\_{i}^{\prime}}{z - \tilde{\xi}\_{i}^{\prime}}\right] \tag{68}$$

Finally, the discrete-time state variable representation of (68) can be obtained, which gives us the required relationship between the reflected and incident voltages necessary to implement the boundary in TLM. A canonical representation for the state-space and output equations is given, respectively, by

**Figure 10.** Geometry and coordinate system for a TE and TM polarized waves incident onto a dispersive half space

described.

260 Advanced Electromagnetic Waves

**5.3. SIBC formulation for tunnel wall modeling**

Thus, their influences have to be studied.

surface impedance for both polarizations:

calculated for both polarizations into the form

*Z* (*s*) is the surface impedance; *s* is the complex frequency; *η* is the intrinsic impedance of the lossy dielectric; , *µ*, and *σ* its constitutive parameters; and *θ<sup>i</sup>* and *θ<sup>t</sup>* are the angles of incidence and refraction, respectively. In tunnel environments, the incident angle *θ<sup>i</sup>* for the waves is unknown. Moreover, multiple reflections of waves with different angles of incidence are involved. To the best of the authors' knowledge, no work has been reported for SIBC under these conditions. A procedure to implement it in any time-domain method as well as its validation is presented in [2, 4]. Here, a brief explanation of its implementation is

The SIBC is inherently a frequency-domain concept. When it is transformed into the time domain for use in methods such as TLM, it is replaced by a convolution integral. Maloney in [21] proposed a recursive formula to overcome the computational difficulties for its calculation. To implement the SIBC in TLM, the reflection coefficient of (4) has to be calculated. This coefficient depends on the polarization, frequency, and angle of incidence.

For tunnel walls, the relative permittivity *<sup>r</sup>* >> 1. Then, for these values in (64), the term *cos*(*θt*) in (Eqs.62, 63) tends to one, obtaining the two angle-independent expressions of the

By using a *rational polynomial approximation* [12], the resulting coefficient of (4) can be

*s* + *σ*/

*s* <sup>√</sup>*s*<sup>2</sup> <sup>+</sup> <sup>2</sup>*α<sup>s</sup>*

*s*<sup>2</sup> + 2*αs* (65)

(66)

*ZV* (*s*) <sup>=</sup> *<sup>Z</sup>*||*pol*−*TM* <sup>=</sup> *<sup>η</sup>*

*ZH* (*s*) = *<sup>Z</sup>*⊥*pol*−*TE* = *<sup>η</sup>*

$$
\begin{bmatrix} \mathbf{x}\_1 \\ \vdots \\ \mathbf{x}\_M \end{bmatrix} = \mathbf{z}^{-1} \begin{bmatrix} f\_1 & 0 & \cdots & 0 \\ 0 & f\_2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & f\_M \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \vdots \\ \mathbf{x}\_M \end{bmatrix} + \mathbf{z}^{-1} \begin{bmatrix} \mathcal{G} 1 \\ \mathcal{G} 2 \\ \vdots \\ \mathcal{G} M \end{bmatrix} V^i \tag{69}
$$

$$V^r = \begin{bmatrix} h\_1 \ h\_2 \ \cdots \ h\_M \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \vdots \\ \mathbf{x}\_M \end{bmatrix} + z^{-1} j\_1 V^i \tag{70}$$

This procedure guarantees a good agreement with theory for a relative permittivity *<sup>r</sup>* >> 1, i.e., it applies to roadway and railway tunnel modeling.

### **5.4. Wave Propagation in a Rectangular Tunnel with Lossy Dielectric Walls Using the SIBC**

The wave propagation in a rectangular tunnel with transverse dimension *w* = 7 m × *h* = 5 m was studied. This tunnel was simulated using the TLM for guiding structures and time-domain SIBC implementations presented previously. The tunnel walls were modeled by using a lossy dielectric with *<sup>r</sup>* = 12 and *σ* = 0.02 Sm−1. The tunnel is shown in Figure 11.

In a rectangular hollow dielectric waveguide, the modes are classified into two families: the first ones with most of its electric field polarized in the *y*-direction, designed as *E<sup>y</sup> nm*, and the other ones with most of its electric field polarized in the *x*-direction, designed as *E<sup>x</sup> nm*.

**Figure 11.** Left: Rectangular tunnel with dimensions w = 7 m ≈ 70 *λ* and w = 5 m ≈ 50 *λ*. Right: Equivalent 2.5D problem terminated by SIBC (dashed line: TLM staircase section approximation)

Using this notation, the figure on the left of Figure 12 shows a comparison in the range of frequency from 0.4 GHz to 2 GHz of the dispersion curves for the modes *E<sup>x</sup>* <sup>11</sup> and *<sup>E</sup><sup>y</sup>* 11 between our procedure and the commonly used Marcatilli's theory for dielectric rectangular waveguides (see Section 2.2). Some very good agreement was obtained for both modes.

**Figure 12.** Dispersion curve for: a) *E<sup>y</sup>* <sup>11</sup> mode and b) *E<sup>x</sup>* <sup>11</sup> mode of the rectangular tunnel calculated in the frequency range from 0.4 GHz to 2 GHz with Marcatilli's theory and our procedure, b) dispersion diagram for a rectangular tunnel in the frequency range from 0.2 GHz to 2 GHz using our numerical procedure [4]

As it can be seen, the dispersion curves for the modes tend to the limit *β* = *ω*/*c*, and up to 3 GHz, they are propagating modes. They propagate with different phase constants and the propagation is strongly multimodal, e.g., around 120 modes at 3 GHz. The figure on the right of Figure 12 shows the comparison between both approaches in terms of attenuation constant. One can notice some discrepancies in the low region of the spectrum. One should stress that (54) is valid for *α* << *β*, and the approximation fails close to the cutoff frequency of the modes (see Section 4.1). Moreover, discrepancies can be explained by the fact that Marcatilli's theory neglects fields at the corners. As the frequency increases, fields are more concentrated in the core region, reducing this error. In Figure 12, the difference between the attenuation constants diminishes as the frequency increases, confirming this fact.

Figure 13 shows the field profiles for the *E<sup>y</sup>* <sup>11</sup> and *<sup>E</sup><sup>x</sup>* <sup>11</sup> modes in the waveguide calculated with Marcatilli's theory. Results obtained in Figure 14 were with the 2.5D TLM approach for guiding structures at 1 GHz. The structure was excited with a Gaussian pulse and injecting the corresponding phase constant of the mode, *β* ≈ 20 in Figure 12. Good agreement was found for the field profiles of both modes.

**Figure 13.** Field configuration in dBV/m of the *E<sup>y</sup>* <sup>11</sup> and *E<sup>x</sup>* <sup>11</sup> modes calculated with Marcatilli's theory [18, 19].

**Figure 14.** Field configuration in dBV/m of the *E<sup>y</sup>* <sup>11</sup> and *E<sup>x</sup>* <sup>11</sup> modes calculated with our procedure [4].

To conclude, it is important to mention that the application of this method can be extended to any uniform waveguiding structure, and also, it can be used to optimally determine antenna field specifications and positioning in confined or waveguide environments as presented in [3].

## **6. Conclusion**

**Figure 11.** Left: Rectangular tunnel with dimensions w = 7 m ≈ 70 *λ* and w = 5 m ≈ 50 *λ*. Right: Equivalent 2.5D

Using this notation, the figure on the left of Figure 12 shows a comparison in the range

between our procedure and the commonly used Marcatilli's theory for dielectric rectangular waveguides (see Section 2.2). Some very good agreement was obtained for both modes.

range from 0.4 GHz to 2 GHz with Marcatilli's theory and our procedure, b) dispersion diagram for a rectangular

As it can be seen, the dispersion curves for the modes tend to the limit *β* = *ω*/*c*, and up to 3 GHz, they are propagating modes. They propagate with different phase constants and the propagation is strongly multimodal, e.g., around 120 modes at 3 GHz. The figure on the right of Figure 12 shows the comparison between both approaches in terms of attenuation constant. One can notice some discrepancies in the low region of the spectrum. One should stress that (54) is valid for *α* << *β*, and the approximation fails close to the cutoff frequency of the modes (see Section 4.1). Moreover, discrepancies can be explained by the fact that Marcatilli's theory neglects fields at the corners. As the frequency increases, fields are more concentrated in the core region, reducing this error. In Figure 12, the difference between the

<sup>11</sup> and *<sup>E</sup><sup>x</sup>*

with Marcatilli's theory. Results obtained in Figure 14 were with the 2.5D TLM approach for guiding structures at 1 GHz. The structure was excited with a Gaussian pulse and injecting

<sup>11</sup> and *<sup>E</sup><sup>y</sup>*

<sup>11</sup> mode of the rectangular tunnel calculated in the frequency

<sup>11</sup> modes in the waveguide calculated

11

of frequency from 0.4 GHz to 2 GHz of the dispersion curves for the modes *E<sup>x</sup>*

<sup>11</sup> mode and b) *E<sup>x</sup>*

attenuation constants diminishes as the frequency increases, confirming this fact.

tunnel in the frequency range from 0.2 GHz to 2 GHz using our numerical procedure [4]

problem terminated by SIBC (dashed line: TLM staircase section approximation)

**Figure 12.** Dispersion curve for: a) *E<sup>y</sup>*

262 Advanced Electromagnetic Waves

Figure 13 shows the field profiles for the *E<sup>y</sup>*

The TLM (*transmission-line matrix*) to characterize modes in guided structures was proposed and described. The 2.5D TLM approach is useful for finding the dispersion diagram, amplitude, and attenuation constant of modes propagating in uniform lossy guiding structures with arbitrary cross section. Its validation in a metallic-rectangular waveguide was presented. This confirms that this approach is well suited for numerical implementation.

Then, the modeling of the wave propagation in lossy waveguides as mine, roadway, or railway tunnels was studied. A procedure to implement the concept of SIBC (*surface* *impedance boundary condition*) to simulate the lossy dielectric walls was described. This concept allows us to avoid the meshing of the tunnel walls by replacing the wall medium by a frequency-dependent reflection coefficient. Thus, field values beyond the air-wall interface are not required for mode parameter computation. However, the presence of the dielectric beyond the interface is accounted for by the SIBC. The SIBC concept yields satisfactory results as long as *ε<sup>r</sup>* >> 1, for arbitrary conductivity. The illustration of the capabilities of this method in a rectangular tunnel shows that our procedure is in agreement with theory for the mode amplitudes, attenuation, and phase constants for the *E<sup>y</sup>* <sup>11</sup> and *<sup>E</sup><sup>x</sup>* <sup>11</sup> modes when compared with Marcatilli's model.

The 2.5D TLM approach for guiding structures and the SIBC concept might be also applied to study the radio-wave propagation for any cross section. The results in terms of amplitudes and dispersion curves and attenuation constants will allow to study different types of structures and the sensibility of these parameters with the polarization, excitation, shape, size, or frequency.
