Three dimensional indoor electromagnetic wave propagation modeling using a hybrid method

(1)

Master Thesis

Three dimensional indoor electromagnetic wave propagation modeling using a hybrid method

Alberto González Escudero

Stockholm, Sweden 2011

XR-EE-ETK 2011:011

(2)

Estudio de propagaci´ on electromagn´etica en interiores en 3D mediante una t´ecnica de

an´ alisis h´ıbrida.

Autor: Alberto Gonz´ alez Escudero

albgonzalez@gmail.com

Tutor: Martin Norgren

mnorgren@kth.se

School of Electrical Engineering Electromagnetic Engineering KTH Royal institute of Technology

Coordinadora Erasmus:

Birgitta Beskow birbes@kth.se

5 de Octubre de 2011

Resumen

El objetivo de este proyecto ha sido el dise˜ nar un método para pre- decir en tres dimensiones la propagación de ondas electromagnéticas en interiores.

Para ello se ha usado un método h´ıbrido, combinando técnicas quasi- anal´ıticas con otras p´ uramente numéricas.

En el m´etodo presentado, las paredes fueron modeladas como Con-

ductor El´ectrico Perfecto (CEP), de modo que los pasillos se compor-

tarsen como gu´ıas de onda.

(3)

Las variaciones en el contorno han sido modeladas mediante Ajuste Modal, mientras que la presencia de giros, puertas o ventanas ha sido modelado usando Diferencias Finitas en el Dominio del Tiempo.

Ajuste Modal es un m´etodo que se basa en encontrar el acoplo entre los distintos modos a ambos lados de una uni´ on entre gu´ıas.

Estos modos son cada una de las soluciones para la ecuacion de Helmholtz con las condiciones de frontera impuestas por el contorno de la gu´ıa.

Este método es capaz de resolver casos en los que se mantenga la simetr´ıa de traslación, y será utilizado para resolver situaciones como puertas entre habitaciones, o estrechamientos en pasillos.

El m´etodo de las Diferencias Finitas en el Dominio del Tiempo (DFDT) se basa en discretizar un dominio, e iterar las ecuaciones rota- cionales de Maxwell (leyes de Amp`ere y Faraday) para ver la evoluci´ on de los campos en el tiempo dentro de este. En el caso de que los cam- pos sean excitados por un pulso, la respuesta en frecuencia del sistema puede ser hallada para un gran numero de frecuencias con una sola simulaci´ on.

La celda usada para crear el enrejado ha sido la celda de Yee, la cual implicitamente satisface las condiciones de divergencia. El los l´ımites del dominio se han situado una condici´on de contorno CEP.

En este proyecto s´olo se han tratado dominios con forma cuboide, en donde cada una de las caras equivaldr´ıa a bien un muro, ventana, o un puerto. Las ventanas han sido modeladas aadiendo una capa perfec- tamente adaptada (CPA), que evita reflexiones al introducir p´erdidas sin variar la impedancia propia del medio.

Para el caso de puertos, también ha sido a˜ nadida una CPA. En cada instante de tiempo, en uno de ellos la excitación es introducida, modulada por una fuente pulsada, para luego ser extraidos los campos eléctrico y magnético en cada uno de ellos, y de este modo estudiar la respuesta del sistema.

Este proceso es llevado a cabo para cada modo en cada gu´ıa.

La capa adaptada permite obtener f´acilmente la matrix de parmet-

ros S, debido a que permite muestear correctamente las ondas de por-

tencia incidentes y reflejadas.

(4)

Esta matriz de parametros S permite describir completamente un sistema el´ectrico lineal. Cada una de las secciones tendr´a asociadas submatrices S donde cada uno de los modos estar´ an representados.

El concatenado de matrices, una vez que los par´ ametros de cada secci´ on han sido extraidos, es inmediato, y permite obtener una ma- trix global que define completamente el sistema.

Los resultados del Ajuste Modal muestran que en los tres casos considerados (cambio de plano E, cambio de plano H y doble salto) coinciden conlos resultados obtenidos mediante HFSS.

Para el caso de los resultados de DFDT, primeramente se compro- baron los resultados obtenidos para una secci´ on de gu´ıa vac´ıa, donde es posible obtener una soluci´on anal´ıtica, para diferente n´ umero de divisiones por longitud de onda. Los resutados fueron satisfactorios, y se procedi´ o a comparar los resultados obtenidos del siguiente caso de inter´es, cuando existe un giro en la gu´ıa.

Cuando hay un giro presente, los resultados no son satisfactorios.

Para el modo T E

¹⁰

se puede ver que ambos resultados son similares, pero presentando m´ ultiples errores. Para modos superiores, los resulta- dos se ven r´ apidamente degradados, llegando al caso de ser totalmente diferentes de los resultados obtenidos mediante HFSS, adem´as de pre- sentar multitud de artefactos num´ericos.

El ´ ultimo caso estudiado ha sido una ventana situada frente a un edificio, con la intenci´ on de modelar c´omo son las reflexiones produci- das por este. para ello se ha situado un puerto en una de las caras, una placa de CEP en la cara opuesta del cuboide, mientras las otras cuatro caras han sido cubiertas con una capa perfectamente adaptada, obteniendose buenos resultados.

Finalmente, los resquisitos de tiempo han sido analizados para el

caso de las Diferencias Finitas, mostrando que el crecimiento en tiem-

po de proceso no s´olo depende del n´ umero de divisiones por longitud

de onda, si no tambi´en del n´ umero de modos considerados, debido a

que las multiplicaciones matriciales necesarias para el muestreo de los

campos el´ectrico y magn´etico igualan o superan las operaciones nece-

sarias para actualizar los campos en el dominio.

(5)

La conclusi´ on obtenida en este proyecto ha sido que el m´etodo no es viable, y que es necesaria una nueva implementaci´on de las DFTF para logar que lo sea, ya que ni el tiempo de proceso se ha visto reduci- do de una manera significativa, ni los resultados son lo suficientemente exactos.

Palabras Clave : Ajuste Modal, Diferencias Finitas en el Dominio del

Tiempo, FDTD, Propagaci´on en Interiores, Scattering, Par´ametros S,

Gu´ıa de Onda, UPML

(6)

Abstract

The objective of this project has been to develop a method to predict three dimensional indoor electromagnetic propagation.For this purpose, a hybrid method has been introduced, consisting in modeling the hallways as waveguides with Perfect Electric Conductor walls, and analyzing the propagation from a modal point of view.

Two different methods were used in order to model changes in the guide, Modal Match- ing for changes in width and height, and Finite Differences in Time Domain (FDTD) for the cases in which the translation symmetry is not maintained.

Cuboids were used for the shape of the FDTD sections. The Yee cell was used for discretization, and Uniaxial Perfectly Matched Layers to model windows and the different ports were the excitation is inserted, and the fields are extracted. This approach allows reflection-free ports, so they can be considered adapted, and therefore a modal scattering matrix representation becomes possible.

The results for the Modal Matching subproblem were satisfactory, as opposed to those obtained by FDTD. For the case of a turn, the results for the mode T E10 were accept- able, despite presenting systematic errors and numerical artifacts, which got worse for the superior modes. The time requirements did not show a dramatic improvement, due to the number of matrix multiplications needed to sample the the fields in the ports located on the faces of the cuboid after every FDTD iteration.

The conclusion obtained is that a new implementation of the FDTD subproblem is needed in order to achieve the objectives of accuracy and computing speed.

Keywords: Modal Matching, Finite Differences in Time Domain, FDTD, Indoor Wave Propagation, Scattering, S-matrix, Waveguide, Uniaxial Perfectly Matched Layer

(7)

To my friends, for helping me become the the person I am today.

To my girlfriend, for standing by my side when I felt lost.

But specially to my family, for supporting me, and giving me the best they could for all this years.

(8)

List of Figures

2.1 ZT M variation with the frequency. . . 15

2.2 ZT E variation with the frequency. . . 16

2.3 Rectangular guide. . . 17

2.4 Dispersion diagram for the 20 first modes for a square guide (a = b) . . . 21

2.5 Dispersion diagram for the 20 first modes for a rectangular guide (a = 2b). . . 21

3.1 Frontal and lateral view of a junction. . . 23

3.2 Intersection (double step) between two guides. . . 28

3.3 Example of a cascade connection. . . 31

4.1 The Yee cell configuration. . . 36

4.2 Example of Yee’s latice. The arrows represent the fields as in Figure 4.1 . . . 41

4.3 Two field configurations used in Modal Matching and in the Yee Cell. . . 47

4.4 Diagram of the process of extracting the scattering parameters from a FDTD section. 50 5.1 FDTD configuration for a uniform waveguide section . . . 54

5.2 Evolution of a simulation for a uniform section. . . 56

5.3 Voltage and current waves for the mode T E10 in a two port uniform section . . . 57

5.4 S21 parameter comparison for a uniform section . . . 58

5.5 Numerical S11 for a uniform section . . . 59

5.6 Absolute error for different divisions per wavelength. . . 59

5.7 FDTD configuration of a turn. . . 61

5.8 Turn excited with the mode T E10 . . . 64

5.9 Turn excited with the mode T E20 . . . 67

5.10 Voltage and current waves for a turn . . . 68

5.11 Absolute value of S11for the mode T E10 prior to concatenate. . . 69

5.12 Absolute value of S11for the mode T E10. . . 70

5.13 Absolute value of S21for the mode T E10. . . 70

5.14 Absolute value of S11for the mode T E20 in a turn. . . 71

5.15 Absolute value of S21for the mode T E20 in a turn. . . 72

5.16 A general configuration of a turn. . . 73

5.17 Example simulation for the mode T E10 facing a reflective plate in free space. . . . 75

5.18 Voltage and current waves of a port facing a reflective plate. . . 76

5.19 S11 for the mode T E10 for a reflective plate in free space calculated using FDTD. 77 5.20 S11 for the mode T E10 for a reflective plate in free space provided by COMSOL. . 77

5.21 CPU time distribution with different amount of modes . . . 79

5.22 CPU time needs variation with the modes and divisions per wavelenght . . . 80

5.23 CPU time needs represented as a 2D-plot . . . 80

5.24 E-plane step. . . 81 4

(11)

List of Figures 5

5.25 Scattering parameters in the big guide for the mode T E10 for E-plane step . . . . 82

5.26 Transmission coefficients for the mode T E10for the E-plane step. . . 82

5.27 Scattering parameters in the big guide for the mode T E20 for the E-plane step. . . 83

5.28 Transmission coefficients for the mode T E20for the E-plane step. . . 83

5.29 S-parameters for the mode T E10 in the small guide for the E-plane step. . . 84

5.30 Transmission coefficients for the mode T E10in the small guide for E-plane step. . 84

5.31 A H-plane step configuration. . . 85

5.32 H-plane’s results for the big guide’s mode T E10 . . . 86

5.33 S-parameters S1,1 to S1,10 for a double step with square concentric guides. . . 87

5.34 S-parameters S1,21 to S1,25 for a double step with square concentric guides. . . 87

5.35 S-parameters S21,21 to S21,25 for a double step with square concentric guides. . . . 88

5.36 S-parameters S21,1 to S21,10 for a double step with square concentric guides. . . . 88 5.37 S-parameters S1,1 to S1,10 for a double step configuration with off-centered guides. 89 5.38 S-parameters S1,21 to S1,25 for a double step configuration with off-centered guides. 89 5.39 S-parameters S21,21 to S21,25for a double step configuration with off-centered guides. 90 5.40 S-parameters S21,1 to S21,10 for a double step configuration with off-centered guides. 90 5.41 S-parameters S4,1 to S4,10 for a double step configuration with off-centered guides. 91 5.42 S-parameters S5,1 to S5,10 for a double step configuration with off-centered guides. 91

(12)

(13)

Chapter 1 Introduction

1.1 Motivation

The need for faster indoor communications rises every day. As the complexity of the systems increases (802.11n, HDSPA+, other MIMO systems), the positioning becomes critical.

Shadowing and reflections from walls or interference with other equipment can occur if the placement is not correct, and an increase in the range and throughput can be achieved with a good positioning.

In order to find the optimal locations, the environment must be well known, which requires empirical modeling, a task that can be costly, and only applicable to existing structures.

To model the medium approximations are needed, and even though in most of the cases an analytical solution can not be found.

It is here where the analysis techniques come into play. This computational techniques try to simulate the behavior of the environment in different ways, and can be classified into two different groups:

• Purely numerical techniques: This group is characterized for its versatility, meaning that it has almost no geometric nor electric limitations, but needs a big amount of computational power. In this group it is possible to find the method of moments (MoM), the finite element method (FEM), and the finite-difference time-domain method (FDTD) among many others [1].

• Quasi-analytical techniques: This methods have a strong analytical content, and thus, can save a lot of processing time, but they are only valid for special cases or situations, where the analytical approach is valid. The Modal Matching Method belongs to this group, and as it will be shown on, Its simplicity leads to orders of magnitude in time saving compared to other numerical methods.

To use a purely numerical technique would be desirable, because whatever the configuration, it will be able to find a solution. But in big scenarios, the computational requirements can be huge.

This is the case with buildings, where the floor planes span from a few tens of square meters to hundreds or even thousands of them.

7

(14)

8 CHAPTER 1. INTRODUCTION

If a building with 1000m² of floor plate, 3 meters of height, at 1GHz and 10 divisions per wavelength is going to be simulated with an FDTD method, the memory requirements for one of the field elements will require an array with ten million elements, meaning that if double precision in wanted, the amount of memory needed to store this array will be close to 80 GB.

Another 5 components will have to be stored, and then operations among this matrices done until the propagation is complete, which can mean to operate them a few ens of thousands of times.

This is impossible current computer power available, at least at a reasonable price.

This is why new methods are needed, and that is why an hybrid method is going to be used: a numerical/quasi-analytical technique to increase the speed of the simulations for indoor simulations in large areas.

The hallways will be modeled as ideal rectangular waveguides, meaning that all the wave propagation through the corridors is done analytically, and when any change in height or with is present will be modeled by using the Modal Matching Method.

When a change in the direction appears, or a window is present, the FTDT section comes into play. As the volume to be modeled is orders of magnitude smaller, the FDTD method is a viable alternative, even though several sweeps shall be done in order to combine the results with those obtained from the quasi-analytical technique

1.2 Organization and purpose

The main goal in this thesis is to achieve a fully working three dimensional hybrid method, and analyze the results comparing them to commercial software, as well as try to improve the simulation speeds without degrading the results.

Finite Differences in Time Domain and Modal Matching will be combined, and hopefully, it will be possible to use the strengths of each method to increase the overall efficiency of the code.

In the second chapter the theoretical bases for this project will be introduced. Modes, guided propagation and the rectangular waveguide are concepts fundamental for understanding the two methods used, and more importantly, the approach used to solve the problem.

In the third chapter, the Modal Matching Method is explained. Modal Matching will be used to solve any discontinuity that does not imply a change in the direction of propagation, that means, change in height or width in the hallway.

The General Scattering Matrix is also introduced here, as well as how to cascade elements defined by the S-matrix.

In the fourth chapter, the FTDT method is explained, including the boundaries problem and the code stability requirements, as well as the scheme used to obtain the frequency-dependent scattering parameters from a time simulation.

This will be used to compute any situation that MMM can not, which includes changes in the direction of propagation, windows, etc..

In the fifth chapter, the results are commented, including the errors produced by the different steps used in the method, and compared with commercial software.

(15)

1.2. ORGANIZATION AND PURPOSE 9

The computer needs will also be analyzed and compared with other methods.

The last chapter reports the conclusions and gives and idea about future work in relation with this project.

(16)

(17)

Chapter 2 Theoretical background

In this chapter the theory beyond guided propagation will be studied, what modes are, and finally, focusing on the particular case of the rectangular waveguide.

The general problem under study is formed by homogeneous media surrounded by perfect electric conductor. Then, after this general case, move to the problem with a single homogeneous medium, analyze the general solutions, and then particularize again for the case of the rectangular waveguide.

2.1 Guided propagation

In the frequency domain, the electromagnetic problem is composed by i homogeneous media, characterized by the magnitudes ǫi and µi, sourceless and surrounded by a perfect electric conductor. This problem is governed by the equations:

∇ · ~Di(~r, ω) = 0 (2.1)

∇ · ~Bi(~r, ω) = 0 (2.2)

∇ × ~Ei(~r, ω) = −jωµⁱH~i(~r, ω) (2.3)

∇ × ~Hi(~r, ω)j = jωǫiE~i(~r, ω) (2.4) By calculating the curl on 2.3

∇{∇ × ~Ei(~r, ω) = −jωµiH~i(~r, ω)} ←→ ∇ × ∇ ~Ei(~r, ω) = ∇ × (−jωµiH~i(~r, ω)) And having

∇ × ∇ × ~Ei= ∇∇ · ~Ei− ∆ ~Ei

Lead us to

∇∇ ~Ei

| {z }

0

−∆ ~Ei = ∇ × (−jωµⁱH~i) (2.5)

Replacing ∇ × ~Hi= jωǫiE~i in 2.5, and defining

γ0i² = −ωµⁱǫi (2.6)

It is possible to write

11

(18)

12 CHAPTER 2. THEORETICAL BACKGROUND

∆ ~Ei− γ0i²E~i= 0 (2.7)

And carrying out a similar process with 2.4:

∆ ~Hi− γ0i²H~i= 0 (2.8)

As the problem has translation symmetry, the nabla operator becomes

∇ ≡ ∇t+ ∇z

The electric and magnetic fields can also decompose into longitudinal and transversal components:

E~i(u1, u2, z) = ~Eti(u1, u2, z) + ~Ezi(u1, u2, z) (2.9) H~i(u1, u2, z) = ~Hti(u1, u2, z) + ~Hzi(u1, u2, z) (2.10) This decomposition allow us to split 2.7:

∆ ~Ei− γ0i²E~i= 0

(∆ ~Eti− γ0i²E~ti= 0

∆ ~Ezi− γ²0iE~zi= 0 (2.11) As ∆ ~Ezi = (∆tEzi+ ∂²Ezi

∂x²)~z is it possible to write:

∆tEzi+∂²Ezi

∂x² − γ0i²Ezi= 0 (2.12)

Applying variable splitting Ezi(u1, u2, z) = FEi(u1, u2)Zi(z) and substituting:

∆tFEi

FEi

+Z_i^′′

Zi − γ²0i= 0





 Z_i^′′

Zi = γ² with γ_ci² = γ_0i² − γ²

∆tFEi− γci²FEi= 0

(2.13)

Following the notation used, γ should have been named γi, but as the fields have to satisfy the boundary conditions for every z, the constant has to be the same in all of them, so the use of subindex is unnecessary.

The solution to the longitudinal electric field in every dielectric region has the form:

Ezi(u1, u2, z) = FEi(u1, u2)e^±γz (2.14) being FEithe solution to the equation

∆FEi− γci²FEi= 0 (2.15)

and following a parallel process for the magnetic field

Hzi(u1, u2, z) = FHi(u1, u2)e^±^γz (2.16) being FHi the solution to the equation

∆FHi− γci²FHi= 0 (2.17)

which are the Helmholtz equation for the E and H field, and γcia complex constant so far unknown.

(19)

2.2. SOLUTIONS TO THE EQUATION 13

2.2 Solutions to the equation

The next step is to find the families of solutions to the Helmholtz equation when there is only one homogeneous medium present. The solutions to this equation can be classified according to the presence of not of ~Ez or ~Hz.

TEM modes

This case correspond to the trivial solution of 2.17 and 2.15. In this case ~Ez= 0 and ~Hz= 0, which means that the fields ~E and ~H has to be completely contained in the perpendicular plains to the z axis. This solutions must fulfill:

∇ × ~Et= −jωµ ~Ht

(∇t× ~Et= 0

∇^z× ~Et= −jωµ ~Ht

∇ × ~Ht= jωǫ ~Et

(∇t× ~Ht= 0

∇z× ~Ht= jωǫ ~Et

(2.18)

By multiplying the last equation by −jωµ and using the second one:

∂²E~t

∂z² − γ0²E~t= 0 (2.19)

Which forces a solution with the form:

E~T EM = ~FEt(u1, u2)e⁻^γ⁰^z

Once the solution for the transversal electric field has been found, the next step is to find the transversal magnetic field. The use of the following relation:

∇z× ~Et= −jωµ ~Ht= −γ0zˆ× ~Et (2.20) Leads to the solution for ~Ht:

H~t= γ0

jωµzˆ× ~Et=zˆ× ~Et

pµ/ǫ = zˆ× ~Et

η (2.21)

This last result also allow us to define the impedance of the TEM mode as:

ZT EM =zˆ× ~Et

H~t

= η (2.22)

Where η in the intrinsic impedance of the medium.

This family of modes will not be of our interest, as it requires at least two conductors to occur.

(20)

TM modes

This is the first family of modes that will be of use in our work. In this case, ~Hz = 0 (the ~H field will be contained in perpendicular plains to ˆz) and ~Ezwill be:

E~z= FE(u1, u2)e⁻^γzzˆ (2.23)

Being FE solution to the equation 2.15 with γ_c²= γ₀²− γ². Having

(∇ × ~E =−jωµ ~H

∇ × ~H = jωǫ ~E →

(∇t× ~Ez+ ∇z× ~Et= −jωµ ~Ht

∇^z× ~Ht= jωǫ ~Et

(2.24)

And multiplying the first equation by jωǫ and using the second:

jωǫ∇t× ~Ez+ ∇z× (∇z× ~Ht)

| {z }

−∂²H~t

∂z²

= −γ²0H~t (2.25)

Rewriting and rearranging terms:

∂²H~t

∂z² − γ0²H~t= jωǫ∇^t× ~Ez (2.26) The solution to the homogeneous case (TEM) has already been studied, so focusing on the TM solution:

H =~ −jωǫ γc²

∇tEz× ˆz (2.27)

In an analogous process it is possible to obtain can obtain:

E~t= γ

γc²∇^tEz and rewrite H~t= jωǫ

γ ˆz× ~Et (2.28) Once the transversal fields have been found, an impedance for the mode TM can be defined as the relation between them:

ZT M =zˆ× ~Et

H~t

= γ

jωǫ (2.29)

Note that this impedance will lack of the physical meaning that had on the TEM case, where the impedance was a relation between the voltage and current waves.

(21)

2.2. SOLUTIONS TO THE EQUATION 15

fc

ZT M

f imag(ZT M)

real(ZT M)

Z _{T M}

0

pµ/ǫ

−∞

Figure 2.1: ZT M variation with the frequency.

TE Modes

This is the second family of modes that will be used.

Following an analogous approach to the previous case, the ~E field will lack of ~Ez component, and ~Hz will be on the form:

H~z= FH(u1, u2) e⁻^γzzˆ (2.30) where FH(u1, u2) is solution to the equation 2.17 with γc² = γ0²− γ. The equations that determine the relations between the fields are:

∇z× ~Et= −jωµ ~Ht (2.31)

∇t× ~Hz+ ∇z× ~Ht= jωǫ ~Et (2.32) By multiplying by −jωµ the equation 2.32, and using 2.31:

∂²E~t

∂z² − γ⁰= −jωµ∇^t× ~Hz (2.33)

Discarding again the TEM solution:

Et= jωµ

γ²c ∇tHz× ˆz (2.34)

H~t= γ

γc²∇^tHz= γ

jωµzˆ× ~Et (2.35)

(22)

The impedance of the TE mode can be defined as the relation between the transversal fields:

ZT E =zˆ× ~Et

H~t

=jωµ γ

fc

ZT E

f imag(ZT E) real(ZT E)

Z _{T E}

0

pµ/ǫ

∞

Figure 2.2: ZT E variation with the frequency.

(23)

2.3. PARTICULAR CASE: THE RECTANGULAR GUIDE 17

2.3 Particular case: The rectangular guide

y

z x b

a

Figure 2.3: Rectangular guide.

The rectangular guide is one particular case of a system with translation symmetry.

In this case, the transversal section will correspond with a rectangle, filled with an homogeneous dielectric, and the walls perfect electric conductor.

The equation for the modes in the rectangular guide is:

∆tF− γ²F = 0 (2.36)

Where F can be FH or FE, and the boundary conditions are:

• ∂FH

∂n

c

= 0, a Neumann boundary condition in the surface of the perfect electric conductor for the modes TE.

• FE

c= 0, a Dirichlet boundary condition in the surface of the PEC for the modes TM Being c the PEC boundary and n the normal to its surface.

In the rectangular case, the boundary will be given by:

(24)

x = 0, 0≤ y ≤ b x = a, 0≤ y ≤ b

y = 0, 0≤ x ≤ a y = b, 0≤ x ≤ a

The function F (x, y) can be expressed as F (x, y) = X(x)Y (y) as the rectangular domain is delimited by two boundary conditions per axis.

The Helmholtz equation:

∆tF− γc²F = ∂²F

∂x² + ∂F

∂y² − γ²c = 0 (2.37)

Can be rewritten using the separation of variables as:

∂ XY

∂x² +∂ XY

∂y² − γc²XY = 0 (2.38)

Or rearranging the terms:

X^′′

X +Y^′′

Y − γc²= 0 (2.39)

Resulting in one equation for each variable where −γc²= kx²+ k²y: X^′′(x)

X(x) = −kx²

Y^′′(y) Y (y) = −k²y

The general solutions to this two equations are:

X(x) = A sin (kxx) + B cos (kxx) (2.40a) Y (y) = C sin (kyy) + D cos (kyy) (2.40b)

Once the solutions have been found, the boundary conditions shall be applied in order to particularize for the TM and TE cases, determining the value of the constants A, B, C, D, kx

and ky.

For the TE modes the following conditions have to be met:











∂FH

∂x

x=0= kx(A cos (0) − B sin(0)) Y = 0 → A = 0

∂FH

∂x

_x=a= kx(A cos (kxa)− B sin(k^xa)) Y = 0 → ky= mπ

a , ∀m ∈ Z











∂FH

∂y

y=0= X ky(C cos (0) − D sin(0)) = 0 → C = 0

∂FH

∂y

y=b= X ky(C cos (kyb)− D sin(kyb)) = 0 → ky= nπ

b ,∀n ∈ Z

(25)

Which leads to the TE solution:

FH= P cos mπx a

cos nπy b

(2.41)

The expression for the longitudinal magnetic field:

Hz,mn= Pmncos mπx a

cos nπy b

e⁻^γ^mn^z (2.42)

In an equivalent way, the boundary conditions for the TM modes are:





 FE

x=0= (A sin (0) − B cos(0)) Y = 0 → B = 0 FE

x=a = (A sin (kxa)− B cos(kxa)) Y = 0 → ky= mπ

a , ∀m ∈ Z





 FE

y=0= Xky(C sin (0) − D cos(0)) = 0 → D = 0 FE

_y=b= Xky(C sin (kyb)− D cos(kyb)) = 0 → ky =nπ

b , ∀n ∈ Z Leading to the TM solution, which is:

FH = Q sin mπx a

sin nπy b

(2.43)

And the associated longitudinal electric field:

Ez,mn= Qmnsin mπx a

sin nπy b

e^−γ^mn^z (2.44)

In both TE and TM cases, the constant γmn will be given by:

γmn= r

−ω²µǫ + mπ a

2

+ mπ a

2

(2.45)

Once the longitudinal fields have been found, obtaining the complete solution is pretty straight-forward. In the case of TE, by using the solution from 2.42 in 2.34 and 2.35 gives the total fields inside the guide

TE^z Ex= −ky

kc

√ǫ1mǫn

√2 a b

√Z cos (kxx) sin (kyy)e^−γz

Ey = kx

kc

√ 1 ǫmǫn

√2 a b

√Z sin (kxx) cos (kyy)e⁻^γz

Hx= −kx

kc

√ 1 ǫmǫn

√2 a b

√Y sin (kxx) cos (kyy) e⁻^γz

Hy = −ky

kc

√ǫ1mǫn

√2 a b

√Y cos (kxx) sin (kyy) e^−γz

Hz= −kc

γ

√ 1 ǫmǫn

√2 a b

√Y cos (kxx) cos (kyy) e⁻^γz

(26)

Where ǫmand ǫn equal 2 when m = 0 or n = 0 respectively, and 1 otherwise.

For the TM case, the solution to be used in order to find the transversal fields is 2.44. By inserting it in the equation 2.28, it is possible to obtain:

TM^z

Ex= kx

kc

√2 a b

√Z cos (kxx) sin (kyy)e⁻^γz

Ey = ky

kc

√2 a b

√Z sin (kxx) cos (kyy)e^−γz

Ez= −kc

γ

√2 a b

√Z sin (kxx) sin (kyy)e⁻^γz

Hx= −ky

kc

√2 a b

√Y sin (kxx) cos (kyy)e⁻^γz

Hy = kx

kc

√2 a b

√Y cos (kxx) sin (kyy)e^−γz

The values of the constants used to describe the fields are:

kx= mπ a

2

ky= nπ b

2

k²c = kx²+ k²y= −γc,mn²

γmn=p

k²c− k²0= r

−ω²ǫµ + mπ a

2

+ nπ b

2

ZT E =jωµ

γ ZT M = γ

jωǫ

The indexes m and n represent the variation in the x and y axis. None of them can be zero for the TM modes, while for the TE modes only one of them can be zero at the same time.

γ can be expressed as an attenuation constant, α, plus a phase constant, β.

Analyzing the expression for γ it is possible to see that below a certain frequency, β will be zero, with gives γ = α, so the mode will not propagate, it will just be attenuated. When above this frequency, β is non zero, meaning that the mode will propagate without attenuation.

This behavior can be seen in the figures 2.2 and 2.1 for the TE and TM modes.

This frequency it is called ’cutoff frequency’ and it is calculated as:

fc=

r mπ a

2

+ nπ b

2

2π√

µǫ (2.46)

Any mode above its frequency will be called propagating mode. If below, it will be named evanescent mode.

This evanescent modes will decay in amplitude exponentially as the waves move in the z axis, so their contribution will be significant only in the region of space where there are excited.

(27)

Figure 2.4: Dispersion diagram for the 20 first modes for a square guide (a = b)

Figure 2.5: Dispersion diagram for the 20 first modes for a rectangular guide (a = 2b).

(28)

In the previous two images it is shown the dispersion diagram (also called Brillouin diagram) of two common guides: the square guide and the rectangular guide with a = 2b. In the images can be seen the appearance of propagating modes as the frequency is increased. The divisions in the horizontal axis indicate multiples of the cutoff frequency of the dominant mode.

2.4 Conclusions

In this chapter the concept of guided propagation has been introduced, the basics beyond the modal propagation, and one particular case, the rectangular guide.

Many other guide configurations are possible, like circular, micro-strip and coaxial, among many others. But due to the similarity between the shape of rectangular guides and real hallways, the study will be limited to this particular case.

All this concepts are going to be useful in the two following chapters, in order to understand the basics of Modal Matching, and the modal excitation and extraction in FDTD.

(29)

Chapter 3 Modal Matching

3.1 Introduction

The modal matching technique is a computer method that analyzes the coupling among modes in a junction between two uniform guide sections.

The main idea is to match the total fields at both sides of the junction so that the power is maintained. Then, deduce the amplitudes of the modes at the output in terms of the mode spectrum at the input.

The power of the method relies on that the relation among the modes can be expressed as the components of a scattering matrix, and cascaded to form an overall scattering matrix when multiple junctions are present.

This method is widely used to design horn antennas and microwave filters and cavities, among others problems where translation symmetry is present and the modes can be analytically known.

3.2 General formulation

S ₀ S _a

Region I Region II

zˆ ˆz

Figure 3.1: Frontal and lateral view of a junction.

23

(30)

24 CHAPTER 3. MODAL MATCHING

The total field in a waveguide can be represented as an infinite weighted sum of the modes in it. The modes, as eigenfunctions, form a complete orthogonal system (see [2], or [3] for a more detailed explanation), meaning that no matter what field is inside the guide, it can be totally decomposed into the modal base, and thus represented as a linear combination of them.

Every mode has an associated field, and this field is composed by waves traveling forward or backward through the longitudinal axis of the guide. Their amplitudes will be defined by the boundary conditions in the transversal surfaces and the beginning and end of it, depending on the excitation and load conditions.

If the guide under study is filled with an homogeneous dielectric, the field will only be composed of TE and TM modes, and the transversal fields can be expressed as :

E~t^I = X∞ i=1

[aÎie^−γÎⁱ^z+ bÎie^+γⁱÎ^z]~eÎi (3.1)

H~t^I = X∞ i=1

[aÎie⁻^γⁱÎ^z− bÎie^+γÎⁱ^z]~hÎi (3.2)

Where aiis the complex amplitude of the incident mode, biis the complex amplitude of the reflected mode, γi is the propagation constant of the mode i and ~ei and ~hi are the transversal electric and magnetic fields of the mode i, all of them in the region I.

At the other side of the junction, in an equivalent way, the fields can be written as:

E~_t^II = X∞ j=1

[aÎI_j e⁻^γ^jÎI^z+ bÎI_j e^+γ^jÎI^z]~e_jÎI (3.3)

H~t^II = X∞ j=1

[aÎIj e^−γÎI^j ^z− bÎIj e^+γ^jÎI^z]~hjÎI (3.4)

Where i and j refers to both TE and TM, and the sums include all the modes, propagating and evanescent.

If the z-origin is set at the junction, as z = 0 all the exponentials become equal to unity, so the fields can be expressed in a neater way as:

E~t^I = X∞ i=1

[aÎi + bÎi] ~eiÎ E~tÎI = X∞ j=1

[aÎIj + bÎIj ] ~ejÎI

H~t^I = X∞ i=1

[aÎi − bÎi] ~hiÎ H~tÎI = X∞ j=1

[aÎIj − bÎIj ] ~hjÎI

(3.5)

Two boundary conditions must be fulfilled at the junction: the continuity of the tangential components of the fields, and the nullity of the transversal electrical field in the surface of the perfect electric conductor:

(31)

3.2. GENERAL FORMULATION 25

E~t^I =

(E~t^II S0

0 Sa− S0

(3.6)

H~t^I = ~Ht^II , S0 (3.7)

Or, by using 3.5 to reexpress the the conditions:

X∞ i=1

[aÎi + bÎi] ~eiÎ = (P∞

j=1[aÎI_j + bÎI_j ] ~e_jÎI S0

0 Sa− S0

(3.8a) X∞

i=1

[aÎ_i − bÎi] ~hÎ_i = X∞ j=1

[aÎI_j − bÎIj ] ~hÎI_j (3.8b)

In this equation, the fields ~e and ~h are supposed to be known, and the relation among amplitudes, both in amplitude and phase, to be unknown.

This equations are not algebraic, as the fields take a different value depending of the spot in the surface of the junction. That is why its resolution it is not direct, and a numeric method is needed.

To convert this equations into an algebraic form, it is needed do the cross product in 3.8a by the magnetic field in the mode m from the guide I and integrate through the surface SI = Sa. Then, take 3.8b and multiply it by the electric field of the mode n in the guide II and integrate trough the surface SII = S0.

It must also be considered that the electric field is annulled in the surface Sa− S⁰.

The orthogonality property of the modes will be used. As the modes are eigenfunctions, the orthogonality is implicit, but a complete proof can be found in [3]:

Z

Sa

(~ei^I× ~h^Im) · ~ds = 0 ∀i 6= m Z

S0

(~en^II× ~h^IIj ) · ~ds = 0 ∀j 6= n

(3.9)

Rewriting the equations in 3.8 using the process described above:

X∞ i=1

(a^Ii + b^Ii) Z

Sa

(~ei^I× ~h^Im) · ~ds = X∞ j=1

(a^IIj + b^IIj ) Z

S0

(~ej^II× ~h^Im) · ~ds X∞

i=1

(a^Ii − b^Ii) Z

Sa

(~en^II× ~h^Ii) · ~ds = X∞ j=1

(−a^IIj + b^IIj ) Z

S0

(~en^II× ~h^IIj ) · ~ds

That applying orthogonality (3.9) will result in:

(32)

(a^Im+ b^Im) Z

Sa

(~em^I × ~h^Im) · ~ds = X∞ j=1

(a^IIj + b^IIj ) Z

S0

(~ej^II× ~h^Im) · ~ds X∞

i=1

(a^Ii − b^Ii) Z

Sa

(~enÎI× ~hÎi) · ~ds = (−aÎIn + bÎIn ) Z

S0

(~en^II× ~h^IIn ) · ~ds

(3.10)

In practice, this infinite sums shall be truncated to be calculated as a finite sum. The number of elements in each side of the equation can be different. From now on, m = 1, . . . , N1 and n = 1, . . . , N2 will be used, were N1 and N2 are the number of terms in each sum.

The equations in 3.10 can be written if a clearer way by making some definitions:

Z

Sa

(~e^Im× ~h^Im) · ~ds = δ1m

Z

S0

(~e^IIj × ~h^Im) · ~ds = Xjm

Z

S0

(~e^IIn × ~h^IIn ) · ~ds = δ2n

Z

Sa

(~e^IIn × ~h^Ii) · ~ds = Xni

(3.11)

And thus, 3.10 can be rewritten as:

(a^Im+ b^Im)δ1m= X∞ j=1

(a^IIj + b^IIj )Xjm

X∞ i=1

(aÎi − bÎi)Xni= (−aÎIn + bÎIn )δ2n

(3.12)

If the modes are not normalized, δ1m and δ2n will not be equal to unity, and thus the S matrix does not retain all of its properties. To normalize the matrix they will be divided by

√δ1mand√ δ2n.

Now, by writing the equations in a matrix form and truncating to N1and N2 terms:

∆1(aÎ+ bÎ) = X^T(aÎI + bÎI) (3.13a) X(aÎ− bÎ) = ∆2(−aÎI + bÎI) (3.13b) Where the column vectors are defined as:

aÎ = [aÎ1aÎ2 . . . aÎN1]^t bÎ = [bÎ1bÎ2 . . . bÎN1]^t aÎI = [aÎI₁ aÎI₂ . . . aÎI_N1]^t

bÎI = [bÎI1 bÎI2 . . . bÎIN1]^t

(3.14)

(33)

3.3. THE SCATTERING MATRIX 27

And the matrices:

∆1= diagonal[δ1m]

∆2= diagonal[δ2n] X = [Xnm]

(3.15)

The two first matrices have size N1 and N2 and the third one has a size N2× N1. The last step is to use the cross product relation:

~hi = 1 Zi

(ˆz× ~eⁱ)

~ei = Zi(~hi× ˆz) And thus

(~ei× ~hj) · ˆz = 1

Zi(~ei· ~ej) = Zi(~hi· ~hj) (3.16) where Zi is the impedance of the ith-mode.

Now it is possible to rewrite the terms as:

Xnm= Ym^I

Z

S0

(~en^II· ~em^I) ds

δ1m= Ym^I

Z

S0

(~emÎ · ~emÎ) ds δ2n= YnÎI

Z

S0

(~en^II· ~en^II) ds

(3.17)

3.3 The Scattering Matrix

The scattering matrix describes the junction as a relation between the incident and reflected amplitudes for each mode at each side.

It is called generalized because it also takes into account the vanishing modes, while the stan- dard scattering matrix used to describe circuits only takes into account the propagating modes.

The relation is:

b^I b^II

=

S11 S12

S21 S22

·

a^I a^II

(3.18)

The S parameters can be found by operating 3.13. S11can be found by clearing b^II in 3.13b and substituting in 3.13a. Applying the conditions for calculating S11(a^II = 0, the second port is perfectly matched), the parameters can be calculated. A similar approach can be done to calculate the rest of the parameters.

(34)

S11= [X X^T + I]⁻¹[X X^T − I] (3.19)

S12= 2[X X^T + I]⁻¹X (3.20)

S21= X^T{I − [XX^T + 1]⁻¹[XX^T− 1]} = X^T{I − S¹¹} (3.21) S22= I − 2X^T[XX^T + I]⁻¹X = I− X^TS12 (3.22)

3.4 The Mutual Coupling Matrix

a1

b1

a2

b2

S0

Sa

Figure 3.2: Intersection (double step) between two guides.

The next step is to see how to calculate the coupling matrix particularizing for the case of two square waveguides.

For this purpose, the problem treated will be the most generic one, a double step junction, even though it is possible to improve the performance significantly by considering two particular cases: the E-plane discontinuity, a discontinuity with only variation in height, and the H-plane discontinuity, in which the width changes [4].

(35)

3.4. THE MUTUAL COUPLING MATRIX 29

The coupling between two modes can be expressed as:

Xpq= Yp^I

Z

S0

(~eq^II· ~ep^I) ds (3.23)

In the double step case, the most generic situation, all the modes shall be considered.

The transversal electric field in the joint for the TE case is:

T Emn^z → ~e =−ky

kc

2√

√ Z ǫmǫn

√abcos (kxx) sin (kyy)ˆx + kx

kc

2√

√ Z ǫmǫn

√absin (kxx) cos (kyy)ˆy

Or rearranging elements:

~e = 2√ Z kc√

ǫmǫn

√ab

(−kycos (kxx) sin (kyy)ˆx) + (kxsin (kxx) cos (kyy)ˆy)

(3.24)

With Z = jωµ

γ for this family of modes.

For the TM modes, the expression is similar:

T Mmn^z → ~e = kx

kc

2√

√Z

ab cos (kxx) sin (kyy)ˆx +ky

kc

2√

√ Z

ab sin (kxx) cos (kyy)ˆy And again, rearranging elements:

~e = 2√ Z kc

√ab

(kxcos (kxx) sin (kyy)ˆx) + (kysin (kxx) cos (kyy)ˆy)

(3.25)

With Z = γ

jωǫ in this case.

Being all the parameters above defined as:

kx1=m1π

a1 ky1= n1π

b1 kx2=m2π

a2 ky2= n2π b2

kc1² =

m1π a1

2

+

n1π b1

2

kc2² =

m2π a2

2

+

n2π b2

2

γ1=p

(−ωµǫ)²+ kc1 γ2=p

(−ωµǫ)²+ kc2

(3.26)

And ǫm and ǫn as defined in 2.3 for every TE.

Once the expressions for the different modes has been stated, the products can be classified as follows:

(36)

TE with TE :

The mode p will be T Em^z1n1 in the guide I and the mode q will be the mode T Em^z2n2 in the guide II.

The coupling expression will be:

Xpq= 4√

Z2

√Z1kc1kc2√ǫp1ǫp2√ǫq1ǫq2√ a1b1√

a2b2

(ky1ky2)i1i2+ (kx1kx2)i3i4

(3.27)

TM with TE :

In this case, the mode p will be T Mm^z1n1 in the guide I and the mode q will be the mode T Em^z2n2 in the guide II.

And coupling expression will be:

Xpq= 4√

Z2

√Z1kc1kc2√ǫq1ǫq2√ a1b1√

a2b2 (−kx1ky2)i1i2+ (ky1kx2)i3i4

(3.28)

TE with TM :

As in analogous way as in the two previous cases, the mode p will be T Mm^z1n1 in the guide I and the mode q will be the mode T E^zm2n2 in the guide II:

Xpq= 4√

Z2

√Z1kc1kc2√ǫp1ǫp2√ a1b1√

a2b2 (−ky1kx2)i1i2+ (kx1ky2)i3i4

(3.29)

TM with TM :

The last scenario is the coupling between two TM modes.The mode p will be T Mm^z1n1 in the guide I and the mode q will be the mode T Em^z2n2 in the guide II:

Xpq= 4√

Z2

√Z1kc1kc2√ a1b1√

a2b2 (kx1kx2)i1i2+ (ky1ky2)i3i4

(3.30)

In all this expressions, i1, i2, i3and i4represent the analytical integrals through the junction.

i1= Z x1

x0

cos (kx1x) cos (kx2(x − x0))dx

i2= Z y1

y0

sin (ky1y) cos (ky2(y − y0))dy

i3= Z x1

x0

sin (kx1x) cos (kx2(x − x0))dx

i4= Z y1

y0

cos (ky1y) cos (ky2(y − y0))dy

Where x0 and y0 denote the beginning of the area where the coupling shall be computed (S0) and x1 = x0 + a1and y1 = y0 + b1 then end of it.

(37)

3.5. CASCADING MATRICES 31

Finite length section

The S-matrix of a waveguide section of length L is:

S =

0 D

D 0

(3.31)

Were D = Diag(e⁻^jγL) and γ is the propagation constant of each of the modes.

This diagonal matrix will include propagating and evanescent modes because if the distance to the next discontinuity is small, the contribution of this non-propagating modes might be considerable.

In the case of propagating modes, the elements in D will have the shape of a complex exponential, while the non propagating will be corresponded by negative exponentials.

3.5 Cascading Matrices

S _L ₁₁ S _L ₁₂ S _L ₂₁ S _L ₂₂

S _R ₁₁ S _R ₁₂ S _R ₂₁ S _R ₂₂

S _T ₁₁ S _T ₁₂

S _T ₂₁ S _T ₂₂

a

^I

a

^{I I}

a

^{I I I}

a

^{I V}

a

^I

a

^{I V}

b

^I

b

^{I V}

b

^{I V}

b

^{I I I}

b

^{I I}

b

^I

Figure 3.3: Example of a cascade connection.

Once the coupling matrix of a junction has been found, it is needed to know how to cascade it with the rest of the elements, that can be uniform waveguide sections or other elements described by its S-parameters (junctions or sections whose scattering parameters have been found

Three dimensional indoor electromagnetic wave propagation modeling using a hybrid method

Master Thesis