Three-dimensional analysis of radiowave propagation over regular and irregular perfectly conducting terrain

ISTANBUL TECHNICAL UNIVERSITY  GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY

Ph.D. THESIS

JANUARY 2018

THREE-DIMENSIONAL ANALYSIS OF RADIOWAVE PROPAGATION OVER REGULAR AND IRREGULAR PERFECTLY CONDUCTING TERRAIN

Zeina EL AHDAB

Department of Electronics and Communications Engineering Telecommunications Engineering Program

Electronics and Communications Engineering Department Telecommunications Engineering Program

FEBRUARY 2018

ISTANBUL TECHNICAL UNIVERSITY  GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY

THREE-DIMENSIONAL ANALYSIS OF RADIOWAVE PROPAGATION OVER REGULAR AND IRREGULAR PERFECTLY CONDUCTING TERRAIN

Ph.D. THESIS Zeina EL AHDAB

(504112308)

Elektronik ve Haberleşme Mühendisliği Anabilim Dalı Telekomünikasyon Mühendisliği Programı

ŞUBAT 2018

İSTANBUL TEKNİK ÜNİVERSİTESİ  FEN BİLİMLERİ ENSTİTÜSÜ

MÜKEMMEL İLETKEN ENGEBELİ VE ENGEBESİZ ARAZİ ÜZERİNDEKİ ÜÇ BOYUTLU RADYO DALGA YAYILIMI ANALİZİ

DOKTORA TEZİ Zeina EL AHDAB

(504112308)

v

Thesis Advisor : Prof. Dr. Funda AKLEMAN ... ISTANBUL Technical University

Jury Members : Ass. Prof. Dr. Güneş Zeynep KARABULUT KURT

Istanbul Technical University ... Doğuş University

Ass. Prof. Dr. Özgür ÖZDEMİR

Istanbul Technical University ...

Prof. Dr. Levent SEVGİ ... Okan Üniversitesi

Prof. Dr. Prof. Dr. Ercan TOPUZ

Doğuş Üniversitesi ... Zeina EL AHDAB, a Ph.D. student of İTU Graduate School of Science Engineering and Technology student ID 504112308, successfully defended the thesis/dissertation

entitled “THREE-DIMENSIONAL ANALYSIS OF RADIOWAVE

PROPAGATION OVER REGULAR AND IRREGULAR PERFECTLY

CONDUCTING TERRAIN”, which she prepared after fulfilling the requirements specified in the associated legislations, before the jury whose signatures are below.

Date of Submission : 11 January 2018 Date of Defense : 19 February 2018

vii

To my family my backbone, to all the teachers who enlightened my path, to Humanity...

ix FOREWORD

This thesis is written as completion to the Telecommunications Engineering PhD program at the Faculty of Electronics and Communications. I was enrolled in the program in September 2011 and started the extensive research on the topic in September 2013. I have experienced this period as very instructive and enriching. On a quite fundamental level this thesis would not have been possible without the support of the following: Ali Kuşçu Science and Technology Scholarship fund, TÜBITAK 2215 program, and the support of the Scientific Research Projects Coordination Unit at Istanbul Technical University (BAP-ITU).

I am deeply grateful to my PhD advisor Prof. Dr. Funda Akleman for her supervision, advice and key scientific insights guiding me to the completion of this thesis. Whenever I lost track on how to proceed she offered helpful guidance and general support. I am lucky to have met and worked with such an amazingly giving person. I also wish to thank Prof. Dr. Ercan Topuz and Ass. Prof. Dr. Güneş Zeynep Karabulut Kurt for their constructive criticism as members of the steering committee of the thesis. During this path, I fully realized the amount of support I have received from my parents over the past decades which allowed me to pursue this research career path. Their strong appreciation for science, enlightenment and knowledge, their fostering of my interests in engineering and physics and their continued support not only during these past years made this thesis possible. Thanking them is not enough.

xi TABLE OF CONTENTS Page FOREWORD ... ix TABLE OF CONTENTS ... xi ABBREVIATIONS ... xiii SYMBOLS ... xv

LIST OF TABLES ... xvii

LIST OF FIGURES ... xix

SUMMARY ... xxiii ÖZET ... xxv 1. INTRODUCTION ... 1 1.1 Significance of Thesis ... 1 1.2 Purpose of Thesis ... 2 1.3 Literature Review ... 3 1.4 Hypothesis ... 5

2. PARABOLIC EQUATION METHOD FOR 3-D VECTOR WAVE PROPAGATION PREDICTION ... 7

2.1 Formulation of the 3-D Vector Parabolic Equation ... 7

2.2 Numerical Evaluation of the 3-D VPE ... 10

2.2.1 Initialization of the algorithm ... 10

2.2.2 3-D VPE algorithm steps ... 13

2.3 Numerical Simulations and Validation of the Algorithm for Simple Basic Propagation Spaces ... 14

2.3.1 Free space simulations ... 14

2.3.2 Propagation over a flat terrain ... 15

2.3.3 Scattering from a PEC knife-edge ... 16

2.3.4 Scattering from multiple PEC knife-edge ... 17

2.3.5 Propagation over a hilly nonuniform terrain ... 18

2.4 Concusion ... 21

3. 3-D ANALYSIS OF NONHOMOGENEOUS ATMOSPHERE PROPAGATION ... 23

3.1 The Formulation ... 23

3.2 Numerical Results ... 25

3.2.1 2-D Simulations for z-varying refractivity profiles ... 26

3.2.2 3-D Simulations for z- and zy-varying refractivity profiles ... 31

3.3 Conclusion ... 40

4. BIDIRECTIONAL VECTOR PARABOLIC EQUATION ... 43

4.1 General Bidirectional Approach ... 43

4.1.1 Formulation ... 44

4.1.2 Backward propagating wave and the bidirectional 3-D algorithm ... 45

4.2 Hz Induction and Computation from zx-Discontinuity ... 46

4.3 Numerical Results ... 49

4.4 Conclusion ... 59

xii

5.1 FDTD-PE Approach and its Implementation ... 61

5.1.1 The FDTD region ... 62

5.1.2 The PE region ... 62

5.1.3 The transition between the FDTD and PE regions ... 63

5.2 Numerical Results ... 64

5.3 Computational and Numerical Efficiency ... 77

5.4 Conclusion ... 78

6. CONCLUSIONS AND RECOMMENDATIONS ... 81

REFERENCES ... 83

APPENDICES ... 89

xiii ABBREVIATIONS

2-D : Two-Dimensional 3-D : Three-Dimensional FD : Finite-Difference FFT : Fast Fourier Transform FEM : Finite Element Method

FDTD : Finite-Difference Time-Domain GUI : Graphical User Interface

MoM : Method of Moments PE : Parabolic Equation PF : Propagation Factor

PEC : Perfect Electric Conductor RCWA : Rigorous coupled wave analysıs RF : Radio Frequency

SSPE : Split-Step Parabolic Equation

SS-VPE : Split-Step Vector Parabolic Equation TE : Transverse Electric

TM : Transverse Magnetic VPE : Vectro Parabolic Equation

xv SYMBOLS

𝐀

⃗⃗ : Magnetic vector potential 𝑨̃̃ : 2-D Fourier transform of the 𝐀⃗⃗ 𝐄⃗ : Electric field vector

𝐅 : Electric vector potential 𝐇

⃗⃗ : Magnetic field vector

𝐉 , 𝐌⃗⃗⃗ : Electric and magnetic current dencities 𝐤_𝟎 : Wave number of free space

𝒌_𝒙, 𝒌_𝑦, 𝒌_𝒛 : Wave number components 𝐧̂ : unit normal vector

𝒏 : reflection coefficient

𝐓𝐄𝐳, 𝐓𝐌𝐳 : Transvers Electric and Magnetic modes

∆𝐬 : Surface impedance

𝛈_𝟎 : Intrinsic impedance 𝛆_𝟎 : Permittivity of free space

Г : Surface impedance reflection coefficient 𝛍_𝟎 : Permeability of free space

𝛔_𝐳 : Standard deviation 𝛚 : angular frequency

xvii LIST OF TABLES

Page Table 5.1 : FDTD Computational Parameters... 67 Table 5.2 : PE Computational Parameters. ... 67 Table 5.3 : FDTD-PE vs FDTD Computational speed. ... 78

xix LIST OF FIGURES

Page Figure 2.1 : Simple 3-D Geometry of the VPE. ... 10 Figure 2.2 : 3-D Gaussian distribution of a source placed at a height z = 2 meters. 11 Figure 2.3 : Gaussian distribution with respect to height (z-direction) with a

transmitter positioned at z = 2 meters. ... 12 Figure 2.4 : 3-D Gaussian distribution of a 2-D source placed at a height z = 2

meters. ... 13 Figure 2.5 : Validation of the SSPE algorithm: Propagation in free space. ... 15 Figure 2.6 : Propagation over perfectly conducting flat terrain. ... 16 Figure 2.7 : Geometry representing propagation over flat earth and an absorbing

knife-edge. ... 16 Figure 2.8 : Propagation over flat earth and an absorbing knife-edge, a) Variations

with respect to the vertical displacement of the receiver.b) Variations with respect to later displacement of the receiver. ... 17 Figure 2.9 : Multiple PEC knife-edge model. ... 18 Figure 2.10 : 120 knife-edges.Vertical variations of the Propagation Factor at a

range x = 25 m, where reciever height of 0 coresponds to the kife-edge height. ... 18 Figure 2.11 : Propagation over a sinusoidal hilly terrain modeled by absorbing knife-edges: a) geometry b) numerical results with different values of the hill width. ... 19 Figure 2.12 : Propagation over a trapezoidal hilly terrain modeled by absorbing

knife-edges: a) geometry b) numerical results with different values of the hill width. ... 20 Figure 3.1 : Refractivity Profiles a) Linear Profile b) Surface-based duct c)

Evaporation duct d) Elevated duct. ... 26 Figure 3.2 : Propagation in a Homogeneous atmosphere. ... 27 Figure 3.3 : Homogeneous atmosphere: 2-D simulations. ... 27 Figure 3.4 : Propagation in an Inhomogeneous atmosphere (Simple Refractivity

Profile). ... 28 Figure 3.5 : Inhomogeneous atmosphere: Simple Refractivity Profile... 28 Figure 3.6 : Inhomogeneous atmosphere: Surface Duct, transmitter at 20m height. 29 Figure 3.7 : Inhomogeneous atmosphere: Surface Duct, transmitter at 100m height.

29

Figure 3.8 : Propagation in an Inhomogeneous atmosphere (Elevated Duct). ... 30 Figure 3.9 : Inhomogeneous atmosphere: Elevated Duct. ... 30 Figure 3.10 : Effect of the inhomogeneous atmosphere on the propagation of the

wave. ... 31 Figure 3.11 : Height variations of the electric field at x = 1000 m and y = 0 for

homogeneous atmosphere and three inhomogeneous atmosphere

xx

Figure 3.12 : Range variations of the electric field at z = 9 m and y = 0 for homogeneous atmosphere and three inhomogeneous atmosphere

profiles. ... 32 Figure 3.13 : Range variations of the electric field at z = 180 m and y = 0 for

homogeneous atmosphere and three inhomogeneous atmosphere

profiles. ... 33 Figure 3.14 : 2-D Refractivity variations : a) 3-D plot Linearly increasing variations in z- and linearly decreasing in y -, b) its corresponding 2-D plot, c) 3-D plot Linearly increasing variations in z- and linearly increasing in y -, d) its corresponding 2-D plot. ... 34 Figure 3.15 : Height variations at x = 1000 m and y = 0 for linear increasing z-

profile and three different y- refractivity profiles. ... 35 Figure 3.16 : Range variations at z = 9 m and y = 0 for linear increasing z- profile

and three different y- refractivity profiles. ... 35 Figure 3.17 : Range variations at z = 180 m and y = 0 for linear increasing z- profile and three different y- refractivity profiles. ... 35 Figure 3.18 : 2-D Refractivity variations : a) 3-D plot surface duct profile in z- and

linearly decreasing in y -, b) its corresponding 2-D plot, c) 3-D plot surface duct in z- and linearly increasing in y -, d) its corresponding 2-D plot. ... 36 Figure 3.19 : Height variations at x = 1000 m and y = 0 for surface duct z- profile

and three different y- refractivity profiles. ... 37 Figure 3.20 : Range variations at z = 9 m and y = 0 for surface duct z- profile and

three different y- refractivity profiles. ... 37 Figure 3.21 : Range variations at z = 180 m and y = 0 for surface duct z- profile and

three different y- refractivity profiles. ... 38 Figure 3.22 : 2-D Refractivity variations : a) 3-D plot elevated duct profile in z- and

linearly decreasing in y -, b) its corresponding 2-D plot, c) 3-D plot elevated duct in z- and linearly increasing in y -, d) its corresponding 2-D plot. ... 38 Figure 3.23 : Height variations at x = 1000 m and y = 0 for elevated duct z- profile

and three different y- refractivity profiles. ... 39 Figure 3.24 : Range variations at z = 9 m and y = 0 for elevated duct z- profile and

three different y- refractivity profiles. ... 39 Figure 3.25 : Range variations at z = 180 m and y = 0 for elevated duct z- profile

and three different y- refractivity profiles. ... 40 Figure 4.1 : Geometry of the two-way problem with two conducting edge. ... 46 Figure 4.2 : Geometry of inclide surface with discontinuity in the zx- plane. ... 47 Figure 4.3 : Single knife-edge, 3-D-VPE simulated to model a 2-D space. (a)

Forward wave. (b) Bidirectional wave. ... 50 Figure 4.4 : Single knife-edge. 2-D versus 3-D total field values between source and knife-edge. Observation at 25 m from source (xobs = 25 m). ... 51 Figure 4.5 : Single knife-edge. 2-D vs 3-D total field values at a height of z = 3m, 2

m below transmitting antenna. ... 51 Figure 4.6 : Single knife-edges. 2-D vs 3-D total field values at a height of z = 9 m.

52

Figure 4.7 : Two knife-edges at 50 and 60 m from source. 3-D-VPE simulated to model a 2-D space. (a) Forward wave. (b) Bidirectional wave. ... 52

xxi

Figure 4.8 : Two knife-edges at 50 and 60 m from source. Height variations between source and first knife-edge. Receiver placed at 30 m from source (xrx =

30 m). ... 53 Figure 4.9 : Two knife-edges at 50 and 60 m from source. Receiver placed between

the two knife-edges at 55 m from source (xrx = 55 m). ... 53

Figure 4.10: Two knife-edges at 50 and 60 m from source. Range variations. Receiver placed at a height of z = 10 m, corresponding to the height of the second edge. ... 54 Figure 4.11 : Single knife edge of width 2λ. PF variations with range for a lateral

displacement y = 1 m and a height of z = 2 m (2λ/3). ... 55 Figure 4.12 : Single knife edge of width 𝜆 and 2𝜆; Lateral variation of the PF for a

receiver at 𝑥 = 3 𝑚 (𝜆) of range and at a height of 𝑧 = 1 𝑚 (𝜆/3) for two knife-edge widths. ... 56 Figure 4.13 : Two knife-edges at 5and 6 meters both of width 2λ; Propagation

Factor variations with range for a lateral displacement y=1m and a height of z = 2m (2λ/3). ... 56 Figure 4.14 : Total Field comparaisons for 2-D and two different 3-D models: the

first with a 2-D source and the second with a D source modeling a 3-D beam. ... 57 Figure 4.15 : Orthogonal triangular prism. a) 3-D Geometry, b) xz-plane. ... 58 Figure 4.16 : Orthogonal triangular prism. The ratio variation between induced Hz

and total Ez field with respect to range. ... 59 Figure 5.1 : 3-D Geometry of the FDTD-PE model: FDTD region, PE region and

transition plane. ... 61 Figure 5.2 : PE initialization matrix: the 𝑁𝑦 𝑥 𝑁𝑧 data matrix (in blue) is extracted

from the 3-D FDTD region along with the PML data (white frame) at 𝑥𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛. Zero padding is then applied to expand its dimensions to conform to the PE space dimensions 𝑁𝑦𝑃𝐸𝑥 𝑁𝑧𝑃𝐸. ... 64 Figure 5.3 : Geometries of the simulation models with the corresponding numerical

paramenters and space partitioning. a) finite edge b) two knife-edges . c) Infinite wall with open window. ... 65 Figure 5.4 : PE Initialization matrix of the electric field (without the zero tappering

region) for Geometry 1: a) Ex component b) Ey component c) Ez component. ... 68 Figure 5.5 : Geometry 1: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥)

variations of field components at 𝑦 = 0. a) FDTD Ex b) PE Ex c) FDTD Ey d) PE Ey. e) FDTD Ez f) PE Ez ... 69 Figure 5.6 : One knife-edge: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 1.5 𝑚 b) z = 3 𝑚. . 70 Figure 5.7 : Two knife-edges: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 1.5 𝑚 b) z = 3 𝑚. . 71 Figure 5.8 : Infinite wall with window: FDTD-PE vs FDTD range plots of electric

field componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 10 c𝑚 b) z = 75c𝑚. ... 72 Figure 5.9 : Geometry 4: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥)

variations of field components at 𝑦 = 0. a) FDTD Ex b) FDTD Ey c) FDTD Ez d) PE Ex e) PE Ey f) PE Ez. ... 74 Figure 5.10 : Cube with window: FDTD-PE vs FDTD range plots of electric field

xxii

Figure 5.11 : Geometry 5: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥) variations of the magnetic field component 𝐻𝑧 at 𝑦 = 0. a) FDTD b) PE. ... 76 Figure 5.12 : FDTD-PE vs FDTD range plots of the magnetic field component 𝐻𝑧, at 𝑦 = 0m for a) 𝑧 = 0.5m b) z = 4𝑚. ... 76 Figure A.1 : FDTD transition matrices at x=7m for a) Ex, b) Ey and c) Ez ... 91 Figure A.2 : Geometry 1: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 0.5 𝑚 b) z = 2 𝑚 ... 91 Figure A.3: Geometry 1: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥)

variations of field components at 𝑦 = 0. a) FDTD Ex b) PE Ex c) FDTD Ey d) PE Ey e) FDTD Ez f) PE Ez. ... 92 Figure B.1: FDTD transition matrices at x=7m for a) Ex, b) Ey and c) Ez ... 93 Figure B.2 : Geometry 2: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 1 𝑚 b) z = 2 𝑚. ... 93 Figure B.3 : Geometry 2: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥)

variations of field components at 𝑦 = 0. a) FDTD Ex b) PE Ex c) FDTD Ey d) PE Ey e) FDTD Ez f) PE Ez. ... 94 Figure C.1 : FDTD transition matrices at x=7m for a) Ex, b) Ey and c) Ez ... 95 Figure C.2 : Geometry 3: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 1 𝑚 b) z = 2 𝑚 ... 95 Figure C.3 : Geometry 3: FDTD-PE vs FDTD plots of range and height (𝑧. 𝑥)

variations of field components at 𝑦 = 0. a) FDTD Ex b) PE Ex c) FDTD Ey d) PE Ey e) FDTD Ez f) PE Ez. ... 96 Figure D.1 : FDTD transition matrices at x=7m for a) Ex, b) Ey and c) Ez ... 97 Figure D.2 : Geometry 4: FDTD-PE vs FDTD range plots of electric field

componenst Ex, Ey and Ez at 𝑦 = 0m for a) 𝑧 = 0.5 𝑚 b) z = 3 𝑚 ... 97 Figure E.1 : FDTD transition matrices at x=6.5 m for a) Ez and b) Hz ... 98 Figure E.2 : Geometry 5 FDTD-PE vs FDTD plots of range and height (z.x)

variations of field components at y = 0. a) FDTD Ez b) FDTD-PE Ez . 98

Figure E.3: Geometry 5: FDTD-PE vs FDTD range plots of electric field

xxiii

THREE-DIMENSIONAL ANALYSIS OF RADIOWAVE PROPAGATION OVER REGULAR AND IRREGULAR PERFECTLY CONDUCTING

TERRAIN SUMMARY

The wireless transmission of information from one location to another or to multiple other locations in space has become an omnipresent need with the evolution of technology. One major method of wireless transmission over the earth’s surface is radio wave transmission. It is at the same time a very challenging analysis problem and has triggered many studies during the last century. Radiowaves are separated according to the atmospheric layer in which they propagate. Radiowave propagation analysis of ground waves depends on both the characteristics of the lower layers of the atmosphere and the characteristics of the terrain over which the waves propagate. One method used for the analysis of long distance propagation is the Parabolic Equation method. It generally consists of splitting the wave equation into two paraxial terms and provides highly accurate results for radio wave propagation over irregular terrain and plane surfaces. It is also used for propagation analysis in urban areas with buildings and multiple scatterers. This method has been applied to various propagation spaces mostly two dimensional with few three dimensional. Knowing that a full wave analysis implies a three dimensional modeling of the propagation space, it is interesting to study the performance of the PE method in this type of spaces. Knowing that each method has its own limitations, it would be of interest to attempt to merge this method along with another numerical method in order to form a robust and efficient hybrid algorithm.

First, the three-dimensional vector parabolic equation method (3D-VPE) is analyzed and applied to different geometries. The vector fields are expressed in terms of electric and magnetic vector potentials, which are initiated at the source and then marched in range by applying the split-step Fourier algorithm. Adequate boundary conditions are applied, which include perfectly conducting (PEC) or lossy obstacles and terrain. The method is tested for simple PEC obstacles over a given range. It is then applied to irregular terrain geometries. The 3-D results are compared to those obtained from 2-D simulations in order to see the effect neglecting the third dimension in the latter numerical simulations. This study of the 3D-VPE permits to see its implementation and performance in 3-D propagation problems.

Second, a 3D-VPE based approach is used to analyze groundwave propagation in an atmosphere with non-homogeneous electromagnetic properties. The inhomogeneity in the atmosphere imply the formation of ducts that result in modifying the propagation direction and range of the waves. One mechanism of wave propagation in the atmosphere is ducting which consists of the bending of waves due to tropospheric refraction and can be strong enough to cause signals to propagate while following the surface of the Earth for both flat and curved earth models. This effect is mostly observed at high frequencies at VHF and UHF and thus the behavior of a signal in an atmosphere with ducts depends directly on the frequency of the propagating signal itself. Ducting is more commonly observed in certain areas than others depending on

xxiv

parameters such as the temperature distribution of the atmosphere, its humidity, the concentration of certain atmospherical particles, etc. In this work, the wave is considered to propagate in an inhomogeneous atmosphere and over a regular terrain. At each marching step, the variations of the index of refraction of the atmosphere are added to the propagator that marches the wave in range. The formulation is explained in details, and the numerical simulations and results are presented where 3-D PE results are compared to 2-D results for different type of atmospherical ducts.

After analysis of the 3D-VPE, and the addition of the variability of the refraction index of the atmosphere, a novel 3-D bidirectional solution to the parabolic approximation of the wave equation is investigated. The backward propagating wave is integrated to the classical parabolic equation approach, which only includes the forward propagating wave. The addition of the backward propagating wave extends the region of validity of the VPE to include regions between the source and the obstacles and/or between multiple obstacles. The propagation over flat terrain in the presence of knife-edges is considered as well as over irregular terrain consisting of hills modeled by the succession of knife-edges. At each knife-edge, appropriate boundary conditions are enforced, and the wave is partly reflected in the backward direction. The wave is marched in both directions by using the split-step algorithm. Different tests are conducted in order to analyze and validate the results obtained by the proposed algorithm. The algorithm is also tested for a geometry inducing coupling between the electric and magnetic fields. All results are compared to those obtained from The Finite-Difference Time-Domain (FDTD) method. The bidirectional algorithm shows a high degree of agreement in the regions targeted by the backward scattering of the wave. The main advantange of the proposed algorithm over other full wave methods, such as the FDTD, is that the PE based algorithm is much faster and requires less computational resources.

Finally, a novel 3-D hybrid model is presented for the prediction of propagation of radiowaves in the presence of complex obstacles. The FDTD method is implemented to model the source and obstacle region and the wave is then propagated by using the Parabolic Equation (PE) method based model. The source region contains all the complex obstacles and takes into consideration all the possible wave interaction in that region. The transition region between the FDTD domain and the PE domain is carefully treated in order to avoid numerical incompatibilities and to realize the time/frequency domain transition. The method is applied for different scenarios; starting with relatively simple obstacles like knife-edges for the verification and validation of the approach and then including more complex obstacles to the propagation medium. The results are compared to those obtained from full FDTD implementations of the same models. The computational speed and memory consumption of the hybrid algorithm are compared to those of the FDTD and to the bidirectional 3D-VPE algorithm. The hybrid algorithm, in addition to providing results conforming to those obtained by FDTD, is much more efficient computationally and much faster. For geometries wherer the bidirectional algorithm can be implemented, it remains the more efficient between the considered three methods.

In conclusion, this thesis includes radiowave propagation analysis approaches which are built on the PE approximation to the wave equaiton. Two novel approach are proposed, the bidirectional 3-D PE approach as well as a novel hybrid FDTD-PE algorithm. They are implemented to simulate different propagation spaces. The main goal to obtain precise electromagnetic field prediction results by the means of fast and efficient numerical approaches is achieved.

xxv

MÜKEMMEL İLETKEN ENGEBELİ VE ENGEBESİZ ARAZİ ÜZERİNDEKİ ÜÇ BOYUTLU RADYO DALGA YAYILIMI ANALİZİ

ÖZET

Gelişen teknolojinin bir sonucu olarak, bir noktadan başka bir noktaya veya birden fazla noktaya kablosuz bilgi aktarımı yapılması her alanda var olan bir ihtiyaç haline gelmiştir. Radyo dalgası iletimi, yerküre üzerinden sağlanan kablosuz iletişime imkân tanıyan yöntemlerin en başında gelmektedir ve aynı zamanda geçtiğimiz yüzyıl boyunca birçok çalışmaya konu edilmiş oldukça zorlu bir analiz problemidir. Radyo dalgaları yayılabildikleri atmosfer katmanlarına göre sınıflandırılırlar. Zemin dalgalarının, radyo dalgası yayılım analizi hem atmosferin alt katmanlarının özelliklerine hem de dalgaların yayıldığı arazinin özelliklerine bağlıdır. Uzun mesafe yayılımının analizi için kullanılan bir yöntem ‘Parabolik Denklem (PE)’ yöntemidir. Bu yöntem dalga denkleminin iki paraksiyel terime ayrılmasından oluşmaktadır ve düzgün yüzeylerdeki ve düzensiz arazilerdeki radyo dalgası yayılımı için çok yüksek doğruluğa sahip sonuçlar sağlamaktadır. Bu teknik aynı zamanda, binaların ve başka birçok farklı saçıcıların bulunduğu kentsel yerleşim bölgelerindeki radyo dalgası yayılımı analizinde de kullanılmaktadır. Bahsedilen metot, çoğunlukla iki boyutlu farklı yayılım uzaylarına uygulanmış ve sayıca çok az olacak şekilde üç boyutlu yayılım uzayı için gerçeklenmiştir. Tam dalga analizi için yayılım uzayının üç boyutlu modellemesinin uygun olacağı bilindiğinden, bu tür uzaylarda PE yönteminin de performansını incelemek ilgi çekicidir. Her yöntemin kendi sınırlamaları olduğu düşünüldüğünde, sağlam ve verimli bir hibrit algoritma oluşturma amacıyla bu yöntemi bir başka sayısal yöntemle birleştirmeye çalışmak da ilginç olacaktır.

İlk olarak, üç boyutlu vektör parabolik denklem yöntemi (3D-VPE) analiz edilmiş ve farklı geometrilere uygulanmıştır. İlk simülasyonda, radyo dalgasının boş uzayda yayılması modellenmiştir. İkinci adımda ise, yayılımın gerçekleştiği mükemmel iletken düz yeryüzü eklenmiş ve etkisini analiz edilmiştir. Üçüncü aşamada ise, yeryüzünün üzerine tek bir mükemmel iletken bıçak sırtı eklenmiş ve yayılmakta olan dalganın bıçak sırtının arkasına saçılmasının etkisi gözlemlenmiştir. Bir sonraki adımda, çoklu bıçak sırtların radyo dalgalarının yayılımı üzerindeki etkisi göz önüne alınmıştır. Son olarak, dağlık arazi gibi düzgün olmayan mukemmel iletken bir yeryüzü üzerindeki yayılım, düz bir arazinin üzerindeki farklı yüksekliklere ve genişliklere sahip tepeleri ardışık mükemmel iletken bıçak sırtları olarak modellenerek simüle edilimiştir. Vektör alanları, kaynak noktasında başlayan ve daha sonra adım adım Fourier algoritması (SSF) uygulanmasıyla tanımlanan aralıkta ilerletilmiş olan elektrik ve manyetik vektör potansiyelleri cinsinden ifade edilmiştir. Mükemmel iletken engellerin olduğu arazi durumları için gereken sınır koşulları uygulanmıştır. Belirtilen yöntem, belirli bir aralıkta basit mükemmel iletken engelleri için test edilmiş ve daha sonra düzensiz arazi geometrileri için uygulanmıştır. Buradan elde edilen iki boyutlu sonuçlar ile üç boyutlu sonuçların karşılaştırılması neticesinde, ilkinde üçüncü boyutun ihmal edilmesinin oluşturduğu etkinin anlaşılması sağlanmıştır. Bu 3D-VPE

xxvi

çalışması, yöntemin uygulamalarının görülmesi fırsatını yaratmış ve üç boyutlu problemlerdeki performansının incelenmesine imkân tanımıştır.

İkinci olarak, 3D-VPE temelli yaklaşım homojen olmayan elektromanyetik ortam özelliğine sahip atmosferdeki zemin dalgaları yayılımını analiz etmek için kullanılmıştır. Atmosferdeki homojen olmayan ortam oluklar oluşmasına neden olacağından, dalgaların yayılım yönünü ve aralığını değiştirecektir. Atmosferdeki dalga yayılımının bir mekanizması, troposferik kırılma nedeniyle dalgaların bükülmesinden oluşan kanaldır ve sinyallerin yassı ve kıvrımlı toprak modelleri için Dünya yüzeyini izlerken yayılmasına neden olacak kadar güçlü olabilir. Bu etki çoğunlukla VHF ve UHF gibi yüksek frekanslarda gözlemlenir ve bu nedenle atmosferdeki kanallarda bir sinyalin davranışı doğrudan yayılım sinyalinin frekansına bağlıdır. Kanal etkileri atmosferin sıcaklık dağılımı, nemi, belli atmosferik parçacıkların konsantrasyonu vb. gibi parametrelere bağlı olarak bazı bölgelerde diğerlerinden daha sık görülür. Yapılan çalışmada dalgaların homojen olmayan atmosfer ortamında ve düzenli bir arazide yayıldığı düşünülmektedir. Her bir yürütme adımında, atmosferin kırıcılık indeksindeki değişim, dalganın yürütüldüğü aralıktaki yayılıma eklenmiştir. Yöntemin formülasyonu detaylarıyla birlikte anlatılmış ve sayısal benzetimler ile elde edilen üç boyutlu PE yönteminin sonuçları ile atmosferik kanalların farklı çeşitlerindeki iki boyutlu sonuçları karşılaştırılarak sonuçlar sunulmuştur. Farklı kırılma indisi türlerini temsil eden dört kırılma profili: standart atmosfer, yükseltilmiş oluk, yüzey oluüu ve buharlaşma oluğu modellenmiştir. Bu oluklar, dalganın farklı yollarla yayılmasını etkiler. Bazıları yüzeye yakın yayılım üzerinde çok az etkili olurken bazıları yüzeye yakın yayılım için yayılım menzilini büyük ölçüde genişletebilir. Yükseklikle Değişken kırılma indısının kırılmanın üç boyutlu algoritmasına eklenmesi, uzun mesafelerde yer dalgalarının üç boyutlu yayılımı ile ilgili daha fazla bilgi sağlar. Bu analiz, atmosferin elektromanyetik özelliklerinin göz önüne alınması gereken bölgelerde alan tahmini yapılması gereken modellerde PE yönteminin doğruluğunu arttırır.

Üç boyutlu vektör parabolik denklem yöntemi analiz edildikten ve atmosferin kırılma indisinin değişimi eklendikten sonra dalga denklemine uygulanan parabolik yaklaşımın yeni bir çift yönlü çözümü araştırılmıştır. Sadece ileri dalga yayılımını içeren klasik parabolik denklemine geri yansıyan dalgaların etkileri de dahil edilmiştir. Potansiyel fonksiyonları kaynaktan başlatılır ve daha sonra Split Step Fourier algoritması uygulanarak menzil yönünde ilerletilir. İlerletilen potansiyel değerleri, yayılımın ana yönüne dik olan düzlemdeki değerlerdir. İlerleyen düzlem bıçak sırtıyla karşılaştığında, bıçak kenarına uygun sınır koşulları potansiyel değerlerine uygulanmiştir. Gelen dalga ileri yönde ilerlemeye devam ederken, yansıyan dalga geriye doğru ilerletilmiştir. Sınır yüzeyinde, yansıyan potansiyel değerler geriye doğru ilerleyen dalga için başlangıç değerleri haline gelir.Geriye doğru yansıyan dalgaların eklenmesi kaynak ile engeller arasındaki ve/veya birden fazla engel arasındaki bölgeleri içerecek şekilde VPE yönteminin geçerlilik bölgesinin genişlemesini sağlamıştır. Yayılımın bıçak sırtı olduğu düz arazi üzerinden olabileceği düşünüldüğü gibi, bıçak sırtları ardışık olarak sıralanmasıyla modellenen tepelerden oluşan düzensiz arazinin üzerinden gerçekleşebileceği göz önünde bulundurulmuştur. Her keskin kenarda uygun sınır koşulları uygulanarak dalganın bir kısmının geriye doğru yansıması sağlanmış ve böylece dalga her iki yöne doğru adım yürütülmüştür. Önerilen algoritmanın denenmesi amacıyla farklı testler uygulanmış ve geçerli sonuçlar elde edilmiştir. İlk olarak, tek bir bıçak sırtı içeren düz bir arazinin üzerindeki yayılım düşünülmüştur ve daha sonra iki bıçak sırtı içerecek şekilde genişletilmiştir.

xxvii

Benzer bir yaklaşım kullanılarak, düz olmayan arazideki yayılımın analizi gerçekleştirilir. Arazi değişken genişlik ve uzunluğa sahip olan bıçak sırtı dizisi olarak modelenmiştir. Söz edilen algoritma ayrıca elektrik ve manyetik alanlar arasındaki kuplajı içeren geometrilerde de test edilmiştir. Elde edilen bütün sonuçlar aynı zamanda Zamanda Sonlu Farklar yöntemiyle de (The Finite-Difference Time-Domain method - FDTD) hesaplanmış ve sonuçların karşılaştırılmıştır. 3D-VPE algoritmasının doğruluğunu kontrol ettikten sonra, farklı üç boyutlu modeller için iki ve üç boyutlu simülasyonların arasındaki farklar gösterilmiştir. İki ve üç boyutlu modelleri karşılaştırabilmek için yalnızca üçüncü boyuttaki uzunlukların (yani, yanal boyut) sınırlı olarak değil, aynı zamanda elektrik kaynağının özellikleri ve boyutları buna göre değiştirilmelidir. Üç farklı üç boyutlu senaryo göz önüne alınmıştır. İlk ikisinde, yanal olarak değişmeyen bir kaynak ve yan genişliği sonlu olan bir bıçak sırtı, üçüncü senaryoda ise üç boyutlu bir işima örüntüsüne sahip bir kaynak ve sonsuz genişliğe sahip bir bıçak sırtı modellenmiştir. İki boyutlu simülasyon yoluyla elde edilen sonuçlar ile yapılan bu karşılaştırmalar, tartışılan parametrelerin alan öngörüleri üzerindeki etkisi ve hangi ölçüde iki boyutlu modelin kesin sonuçlar sağlamaığı belirlenmiştir. Çift yönlü algoritma, dalganın geriye doğru saçılmasıyla hedeflenen bölgelerde yüksek bir eşleşme göstermiştir. Önerilen PE temelli algoritmanın en büyük avantajı FDTD tekniği gibi diğer tam dalga yöntemlerine göre çok daha hızlı olması ve hesaplama maliyetinin çok daha az olmasıdır.

Son olarak, karmaşık engellerin olduğu durumlarda radyo dalgası yayılımını öngörmek için yeni ve özgün bir üç boyutlu hibrit model geliştirilmiştir. Kaynak ve engellerin olduğu alandan oluşan modele FDTD yöntemi uygulanmış ve daha sonra PE temelli model kullanılarak yayılımı sağlanmıştır. Kaynak ile birlikte tüm engeller ve karmaşıklıklar FDTD yoluyla tam dalga analizinin gerçekleştirildiği bölgede bulunmaktadır. Verilen modele göre, sayısal uyuşmazlıklardan kaçınmak için, FDTD ve PE bölgeleri arasındaki birleştirme konumu olarak uygum bir mesafe dikkatlice seçilir. Parabolik denklemin başlangıç vektörü olarak ayarlanan iki boyutlu matris daha sonra SSPE algoritmasına aktarılır ve alıcıya kadar ilerletilir. Alanlar, FDTD alanında tanımlanan bir Gauss zaman kaynağı tarafından başlatılır. Zaman domeninin alanları FDTD-PE bağlantı arayüzünde biriktirilir ve PE uzayının başlangıç vektörlerini elde etmek için Fourier Dönüşümü kullanılarak spektral alana dönüştürülür. FDTD ve PE bölgelerinin ızgara aralığı ve hesaplama alanının boyutu gibi sayısal parametrelerin seçimi, sonuçların doğruluğunda hassas bir rol oynar ve daha ayrıntılı olarak ele alınır. Sonuç olarak, kaynaktan alıcıya bütün menzildeki alan değerleri, elde edilir. Anlatılan yöntem bir çok farklı senaryoya uygulanmıştır. Öncelikle, göreceli olarak daha basit olan bıçak sırtı engellerin olduğu durumlara uygulanarak yöntemin geçerliliği doğrulanmıştır ve algoritma daha sonra yayılım ortamında daha karmaşık engellerin olduğu senaryolara uygulanmıştır. Sonuçlar aynı senaryolara uygulanan FDTD yönteminin sonuçları ile karşılaştırılmıştır. Hibrit yönteminin hesaplama hızı ve hafıza kullanım maliyeti FDTD yöntemininki ve çift yönlü 3D-VPE yöntemininki ile kıyaslanmıştır. Hibrit algoritmasının uygulanması sonucunda, çok daha hızlı bir şekilde ve daha verimli hafıza kullanımıyla FDTD yönteminin sonuçlarını elde edebildiği görülmüştür. Çift yönlü algoritmanın uygulanabildiği geometrilerde, üç yöntemin arasındaki en verimli yöntem olduğu görülmüştür.

Sonuç olarak bu tez, dalga denklemine uygulanan PE yaklaşımı üzerine kurulmuş olan radyo dalgası yayılım analizi yaklaşımlarını içermektedir. İki özgün yaklaşık yöntem olan, çift yönlü 3-D PE yaklaşımı ve FDTD-PE hibrit yaklaşımı önerilmektedir. Bu

xxviii

yöntemlerin farklı yayılım uzaylarında benzetimleri gerçekleştirilmiştir. Temel hedef olarak görülen elektromanyetik alanların kesin değerlerinin, çok hızlı ve etkin sayısal yöntemlerle yaklaşık olarak elde edilmesi işlemi yapılan çalışmalar sonucunda başarılmıştır.

1 1. INTRODUCTION

1.1 Significance of Thesis

The work in this thesis concentrates on ground waves whose characteristics depend on the lower layers of the atmosphere and the characteristic of the ground over which the propagation takes place. Whether over long rural distances characterized by inhomogeneous terrain, or over short urban distances, few are the cases in which a simple direct line-of-sight is present between the transmitter and the receiver. Such environments are dynamic, and the path between the transmitter and the receiver can considerably vary from a direct line-of-sight to one that undergoes multiple scatterings due to obstacles such as buildings and mountains. This implies that the performance of the radio transmission system directly depends on the propagation environment, which causes severe limitations to the transmission itself. This dependence on the propagation environment makes the determination of ground wave propagation characteristics complex between two points, and leads to site-specific analysis based on statistical measurements performed in the transmission area where the communication system is intended to operate. Although these methods can provide accurate solutions, their accuracy remain dependent on the physical and electromagnetic properties of the environment, which limits their application to specific geometries or spaces. Hence, they remain incomplete in the sense of providing a general solution to the wave problem. Here rises the need to develop more efficient methods and algorithms that could provide valid and reliable solutions while trying to analyze radio wave propagation models. Approaches that provide solutions to general propagation problems can be categorized based on the techniques they use. Some use analytical approximations to the wave equation; others use numerical approaches while other methods include both numerical and analytical approximations, whether 2-D or 3-D.

2 1.2 Purpose of Thesis

A three-dimensional deterministic analysis of radio wave propagation is important for both national commercial and military applications. The aim of this work is to bring together and implement a new 3-D hybrid numerical method, which can be generally applied to the analysis of wave propagation over non-homogeneous terrain. Due to the complexity of the problem, applications of 3-D methods remain few in the literature. In order to develop a solid 3-D hybrid method, existent numerical methods are analyzed in 3-D space in order to determine their regions of validity and their numerical efficiency. The work is initialized by analyzing 3-D Parabolic Equation (PE) based approaches, and developing them in order to extend their region of application. The main goal of this thesis is to design an algorithm to predict signal levels for electromagnetic wave propagation over non-homogeneous terrain including atmospheric effects, which requires the solution of Maxwell’s equations in a 3-D medium. Since a full-wave 3-D analytic solution has not yet appeared in literature, appropriate numerical techniques are applied together with hybridization procedures. The main constraint in such numerical problems depends on the ratio between the computation space and the wavelength of the propagating waves, which is relatively high for long-range propagation. Therefore, statistical models are usually used in conventional radio frequency (RF) planning. However, models that are more realistic are required especially for emerging technologies. These models can be used as pre-design analysis tools instead of site measurements that are cumbersome and expensive. The proposed algorithm will enable the analysis of radiowave propagation within a certain frequency range by using geographical and atmospheric data.

The transition from a 2-D model to a more realistic 3-D model helps to improve the accuracy of a given method but at the same time requires higher computational sources and increases time cost. For instance, previously developed 2-D PE algorithms do not take into account the laterally propagating waves which results in an incomplete model [1], whereas 3-D PE methods based on the split-step algorithm (Split-Step Parabolic Equation: SSPE) have been developed and provide clearly more accurate results [2]. However, the SSPE algorithm, which is used to solve the parabolic wave equation, does not take into account the backscattered waves whose effect becomes un-neglectable when the receiver is positioned between the transmitter and a given obstacle or when the propagation takes place over a back-scattering irregular terrain.

3

Hence, in this work, the addition of the backscattered waves to the 3-D SSPE algorithm will be considered and is expected to improve the results obtained for regions that are not completely modeled by the one-way SSPE approach.

Other aspects to be considered in the propagation of ground waves is the atmosphere in which the propagation takes place and specifically its variations with respect to the refraction index. These variations affect directly the propagating wave, thus should be taken into consideration in the model to be analyzed.

As a result, the aim of this study is to merge efficiently two methods in one efficient, accurate and widely applicable hybrid method that will have high impact on the current radio wave propagation analysis methods.

1.3 Literature Review

Many two-dimensional (2-D) implementations of various numerical methods used in the analysis of radio-wave propagation for both indoor and outdoor spaces can be found in the literature. In these studies, numerical methods such as Finite-Difference (FD) [3], Method of Moments (MoM) [4,5], Finite Element Method (FEM) [6], Finite-Difference Time-Domain (FDTD) [7], Integral Equation [8] and the Tabulated Interaction Method [9] have been successfully implemented in both closed and open 2-D spaces for which deterministic results are provided. However, the implementation of these methods leads to the formation of large matrices that need to be solved, which requires high processing sources even for 2-D numerical spaces. This need for high processing sources increases exponentially in 3-D numerical spaces. In order to overcome this limitation, methods based on analytical approximations such as Ray Tracing methods [10] and the Parabolic Equation method (PE) [1] are used and result in numerically less demanding algorithms. Under appropriate conditions, these methods provide satisfactory results when implemented to 2-D numerical spaces. For instance, the PE method provides results that are valid when applied under certain approximations. It is used for long distance propagation and propagation in urban areas with buildings and multiple scatterers [2]. It consists of one of the methods that are often considered when combining efficiency and far field computations. A lot of 2-D PE based analysis can be found in the literature where the propagation in the third dimension is neglected [11-16]. Few are the formulations found in the literature, which model full 3-D propagation medium [17-19]. Still, these 3-D models are based either

4

on scalar PE formulation or on the implementation of noncoupling boundary conditions. More general implementations of the PE method where vector fields are treated on the boundaries are presented in [20,2]. In the first work, the PE is numerically discretized by using the FD technique, whereas in the second work the Fourier transform based split-step (SS) algorithm is used to march the waves in range. 3-D finite difference-based PE is also successfully applied to model transmission in closed environments consisting of tunnels of different geometries in [21,22].

However, each of the mentioned method has its own limitations and regions of validity that directly depend on the geometry of the problem. For instance, a method can provide excellent results in a specific region while providing unreliable results in another region. In some cases, these methods can be numerically highly demanding and result in non-realizable models. An example to that is the transition from a closed space to an open space while analyzing the propagation of electromagnetic waves [23]. In practice, this problem is equivalent to the determination of the propagation characteristics of an electromagnetic wave in a topologically different terrain outside a given structure (e.g. building) while the transmitter is inside it. In such situations, the development of hybrid algorithms that are appropriate to the region of interest becomes meaningful and gains importance. In the literature, various applications of hybrid methods are found like the FDTD/Rigorous coupled wave analysıs (RCWA) [23], Ray-Optics/ParFlow[24], Parabolic Equation/ Equivalent Current [25], PE/MoM [26] PE/Integral Equation method [27] and PE/Gaussian Beam [28]. The majority of these implementations are two-dimensional (2-D) and the fact that they cannot provide an exact model of the real problem limits their applications. A 3-D implementation of a hybrid Finite Difference (FD)/PE method is represented in [29] for propagation in a closed space, such as straight and curved tunnels. An interesting hybrid approach coupling the FDTD method with the PE method for near field-far field analysis is presented in [30,31], where it is applied to predict 2-D sound propagation in the presence of complex obstacles and wind profiles. The Latter corresponds to solving the acoustical wave equation governing the propagation of sound and the main difference is that the entities computed are 2-D scalar whereas the electric and magnetic fields are 3-D vector entities.

5 1.4 Hypothesis

The hypothesis on which the thesis is based on can be stated as follow, “The numerical analysis of radio wave propagation problems an inhomogeneous atmosphere and over inhomogeneous terrain can be highly improved in terms of accuracy and speed by means of a 3-D hybrid numerical algorithm uniting full wave analysis and approximation based analysis”.

The modeling of a radiowave propagation problem is closest to reality when realized in 3-D. This modeling takes into consideration the propagation of waves that are not accounted for in 2-D modeling (e.g. lateral waves), and hence provides more realistic and accurate results. As stated in the previous section, numerical methods have weaknesses in terms of performance (i.e. speed of computation and accuracy) depending on the space in which they are implemented, which makes the use of one single method in 3-D analysis restraining. In order to overcome this, a novel hybrid numerical algorithm that combines multiple methods will be able to provide a better solid performance that benefits from the strength of a method in one region and overcome its weaknesses by using another method in other regions. The proposed method will add to the literature a new approach with improved performance over already existing methods, which can be implemented while solving radiowave propagation problems over inhomogeneous terrain, and thus will help in the development of new simulation tools.

7

2. PARABOLIC EQUATION METHOD FOR 3-D VECTOR WAVE PROPAGATION PREDICTION

In this chapter, the detailed formulation of the solution of the 3-D wave equation using the vector parabolic approach is presented. Several simulations of propagation models including different terrain types and simple PEC obstacles are conducted. Comparisons between 3-D and 2-D results are shown in order to test the validity of the 3-D algorithm and the effects of the lateral wave diffractions on the 2-D simulation results. The results obtain from this work were presented in [32].

2.1 Formulation of the 3-D Vector Parabolic Equation

Considering the propagation in a homogeneous atmosphere where the index of refraction is that of free space, and supposing a time dependency of 𝑒−𝑗𝜔𝑡, the electric and magnetic fields can be expressed in terms of two magnetic and electric vector potentials, respectively 𝐴 , and 𝐹 [33], as

𝐸⃗ = 𝑗𝜔 (𝐴 + 1 𝑘02∇ (∇. 𝐴 )) − 1 𝜀0(∇ x 𝐹 ), (2.1) 𝐻⃗⃗ = 1 𝜇₀(∇ x 𝐴 ) + 𝑗𝜔 (𝐹 + 1 𝑘₀2∇ (∇. 𝐹 ) ), (2.2)

where 𝜀₀ and 𝜇₀ are the permittivity and permeability of free space respectively, 𝑘₀ is the free space wavenumber.

Supposing that the fields are excited by vertically polarized (z-oriented) sources at a given range 𝑥₀, they can be expressed as,

𝐸𝑧𝐴 = − 𝑗𝜔 𝑘2 (𝑘2+ 𝜕2 𝜕𝑧2) 𝐴𝑧 = 𝐸𝑧 , (2.3) 𝐻_𝑧𝐹 _{= −}𝑗𝜔 𝑘2 (𝑘2+ 𝜕2 𝜕𝑧2) 𝐹𝑧 = 𝐻𝑧 . (2.4)

8

The total electric and magnetic fields are decomposed into two sets of three TMz and

three TEz components [2] and are given as,

𝐸_𝑥𝐴 =𝑗𝜂0 𝑘₀ 𝜕2_𝐴 𝑧 𝜕𝑥𝜕𝑧 , (2.5) 𝐻_𝑥𝐴 ₌𝜕𝐴𝑧 𝜕𝑦 , (2.6) 𝐸𝑦𝐴 = 𝑗𝜂₀ 𝑘0 𝜕2𝐴_𝑧 𝜕𝑦𝜕𝑧 , (2.7) 𝐻_𝑦𝐴 = −𝜕𝐴𝑧 𝜕𝑥 , (2.8) 𝐸_𝑥𝐹 _{= −}𝜕𝐹𝑧 𝜕𝑦 , (2.9) 𝐻_𝑥𝐹 = 𝑗 𝜂₀𝑘₀ 𝜕2_𝐹 𝑧 𝜕𝑥𝜕𝑧 , (2.10) 𝐸_𝑦𝐹 ₌𝜕𝐹𝑧 𝜕𝑥 , (2.11) 𝐻_𝑦𝐹 = 𝑗 𝜂₀𝑘₀ 𝜕2_𝐹 𝑧 𝜕𝑦𝜕𝑧 . (2.12)

Using this representation of the fields, the propagation problem can be solved in terms of the potentials, who have solutions that satisfy the following equations as in [2],

∇2𝐴 + 𝑘₀2𝐴 = −𝜇0𝐽 , (2.13)

∇2𝐹 + 𝑘₀2𝐹 = −𝜀0𝑀⃗⃗ , (2.14) where the magnetic current source 𝑀⃗⃗⃗⃗ = 0, and electric current source 𝐽 = 𝐽_𝑧𝑒 _𝑧 at the source and zero elsewhere.

Applying the wide-angle parabolic equation approximation to the (2.13) and (2.14) in a source-free space gives the following equations as in [2],

9 𝜕𝐴_𝑧 𝜕𝑥 − 𝑗𝑘0√1 + 1 𝑘₀2( 𝜕2 𝜕𝑦2+ 𝜕2 𝜕𝑧2) 𝐴𝑧 = 0 , (2.15) 𝜕𝐹_𝑧 𝜕𝑥 − 𝑗𝑘0√1 + 1 𝑘₀2( 𝜕2 𝜕𝑦2+ 𝜕2 𝜕𝑧2) 𝐹𝑧 = 0 , (2.16)

which represent the forward propagating electric and magnetic potentials.

In addition, the potentials should satisfy the boundary conditions on the boundary of the obstacles. Surface impedance boundary conditions are applied and are expressed as in [2] as,

𝑛̂ × 𝐸⃗ = ∆_𝑠𝜂₀𝑛̂ × (𝑛̂ × 𝐻⃗⃗ ), _(2.17) where 𝑛̂ is the unit normal vector, ∆_𝑠 is the surface impedance of the boundary and 𝜂₀ is the intrinsic impedance of the medium.

This results in the following set of boundary conditions for both the PEC ground and vertical surfaces as [2],

∂𝐴

𝜕𝑛 = 0, 𝐹 = 0

for the ground, _(2.18)

𝐴 = 0, ∂𝐹 𝜕𝑛 = 0

for vertical surfaces. _(2.19)

While applying these boundary conditions to the wave equation in a perfectly conducting medium, the solution of the problem is quite straight forward and the potentials can be easily computed at each marching step. Nevertheless, the case is not so when we consider lossy obstacles, i.e. dielectric, and the solution of the problem becomes complex as it consists of solving complex matrices which require a high level of computation. Details of these formulations are shown in [2].

10

The following Section 2.2 includes a detailed presentation of the Split-Step Vector Parabolic Equation (SS-VPE) algorithm as well as results obtained from different propagation models.

2.2 Numerical Evaluation of the 3-D VPE

The PE method consists of approximating the wave equation in a way that models the emitted wave from a given source as a conical beam propagating around a given axis or direction [34]. In order to solve the parabolic wave equation for a simple geometry model with two screens as shown in Figure 2.1, the initial vector potential generated by the source is propagated in range after being initialized at 𝒙 = 𝟎 . The potential at the next range step 𝒙 = 𝒙 + ∆𝒙 is computed by using the Split-Step Fourier approach.

Figure 2.1 : Simple 3-D Geometry of the VPE.

Details of the algorithmical steps are presented in section 2.2.1. Results of different tests and simulations are presented and explained in section 2.2.2.

The first step is to develop the 3-D VPE algorithm based on the SSPE approach that solves the 3-D VPE and to start testing it for simple propagation environments. The algorithm is based on the vectorial formulation presented in section 2.1 and on the approach suggested in by Janaswamy in [2].

2.2.1 Initialization of the algorithm

The source excitation is chosen as a vertically polarized current source with a 3-D source with a 1-D Gaussian aperture distribution at central lateral position 𝑦 = 0 as shown in Figure 2.2 and Figure 2.3 and given by [2],

Source

Forward wave

x z _y

11 𝑔(𝑦, 𝑧) = 1

σ_z√2𝜋𝑒

−(𝑧−𝐻𝑡)2

2𝜎𝑧2 . 𝛿(𝑦), (2.20)

where 𝐻_𝑡 represents the height of the transmitting antenna, and 𝜎_𝑧 is the source standard deviation used to set the 3 dB elevation beam width of the transmitting antenna. Furthermore, the initial magnetic vector potential due to an electric current of amplitude 𝐼₀ can be expressed as given in [2] and [35] as,

𝐴(0+, 𝑦, 𝑧) ≈ 𝑗𝜇0𝐼0𝑙

2 𝑔(𝑦, 𝑧) ,

(2.21)

where 𝑙 is the short dipole length and 𝜇0 is the permittivity of free space. As for the

initial electric vector potential value 𝐹(0+_{, 𝑦, 𝑧), it is set to zero because no magnetic}

current is set to excite the source [16].

Figure 2.2 : 1-D Gaussian aperture distribution of a 3-D source at central lateral position 𝑦 = 0 placed at a height 𝑧 = 2 meters.

12

Figure 2.3 : Gaussian distribution with respect to height (z-direction) with a transmitter positioned at 𝑧 = 2 meters.

It is to be noted that, as will be explained in the following section, in order to apply surface boundary conditions on the flat terrain considered as a perfectly electric conducting (PEC) surface, the image theory is applied, implying that the fields and their corresponding vector potentials are represented in terms of odd and even parts. This infers that the Gaussian source representation is equally split into an odd component and an even component.

In this thesis, in examples that include 2-D/3-D comparisons of propagation models, the source in the 3-D algorithms is initialized with a 2-D source in a 3-D space. This is realized by setting the initial height distribution as in Figure 2.3, and then repeating it in the lateral dimension, which is the dimension that is not modeled in the 2-D space. This repetition in the lateral dimension permits the elimination of source incompatibilities between the 2-D and the 3-D simulations, which enables the analysis of the effects of 3-D obstacles with identical sources. This process results a 3-D source (2-D distribution at 𝑥 = 0) shown in Figure 2.4 in contrast to that in Figure 2.2.

13

Figure 2.4 : 1-D Gaussian aperture distribution of a 3-D source placed at a height z = 2 meters with no lateral variation.

2.2.2 3-D VPE algorithm steps

After initializing the field 𝐴(0+, 𝑦, 𝑧) at 𝑥 = 0+, field values at step 𝑥( 0++ ∆𝑥) are computed by applying the propagator 𝑒𝑗𝑘𝑥∆𝑥 _{to the Fourier transformed} field 𝐴̃̃(0+, 𝑘𝑦, 𝑘𝑧). Hence, obtaining the field values at the successive location 𝑥 +

∆𝑥 by an inverse Fourier transform. As can be seen in equation 2.22, the double integration corresponds to an inverse Fourier transform where 𝑘_𝑦 and 𝑘_𝑧 are the Fourier domain variables and is given by,

𝐴_𝑧𝑓(𝑥 + ∆𝑥, 𝑦, 𝑧) = 1 4𝜋2∬ 𝑒𝑗𝑘𝑥∆𝑥𝐴̃̃𝑧 𝑓 (𝑥, 𝑘_𝑦, 𝑘_𝑧) ∞ −∞ 𝑒 𝑗(𝑘𝑦𝑦+𝑘𝑧𝑧)_𝑑𝑘 𝑦𝑑𝑘𝑧 , (2.22) where 𝑘_𝑥 = √𝑘₀2_{− 𝑘}

𝑦2− 𝑘𝑧2 is the wave number in the x direction and 𝐴̃̃ is the 2-D

Fourier transform of the potential 𝐴 and is given by,

𝐴̃̃(𝑥, 𝑘_𝑦, 𝑘_𝑧) = ∫ ∫∞ 𝐴(𝑥, 𝑦′_{, 𝑧}′_)𝑒−𝑖(𝑘𝑦𝑦′)_{. [𝑒}−𝑖(𝑘𝑧𝑧′)₊ 𝑦′_=−∞ ∞ 𝑧′₌₀ Г(𝑘_𝑧)𝑒𝑖(𝑘𝑧𝑧′)_]𝑑𝑘 𝑦𝑑𝑘𝑧 , (2.23)

where Г(𝑘𝑧) is the reflection coefficient of an impedance surface and is given for 𝑇𝑀𝑧

14 Г(𝑘𝑧) = { 𝑘𝑧− ∆𝑠𝑘0 𝑘𝑧+ ∆𝑠𝑘0 for TMz, ∆_𝑠𝑘_𝑧− 𝑘₀ ∆𝑠𝑘𝑧+ 𝑘0 for TEz. (2.24)

The value of the reflection coefficient is appropriately chosen for 𝑨 and 𝑭 from (2.24). This is the series of operations to be applied in order to obtain field values at any new range position, given that there are no obstacles and that the wave propagates in free space. Supposing that the wave propagates over a perfectly conducting terrain, and that electrical obstacles are present, proper boundary conditions should apply. The wave can be propagating in an urban environment where the obstacles can be buildings, houses, or other edifices or in a rural environment where there can be hills, mountains and trees. It is supposed that these obstacles can be modeled as a succession of perfectly absorbing knife-edges.

2.3 Numerical Simulations and Validation of the Algorithm for Simple Basic Propagation Spaces

2.3.1 Free space simulations

The first simulation is that of the propagation of a radio wave in free space and the results are shown in Figure 2.5. A Cartesian coordinate system is used, where x is the direction of propagation, y is the lateral direction and z is the vertical direction. The transmitter height is 50 meters, and the receiver height is 2.1 meters. The operating frequency is 1 GHz and the y and z Fourier parameters are 𝑁𝑦 = 4096 and 𝑁𝑧 = 2048. The source is set as in Figure 2.2. The performance of the VPE algorithm is tested against the analytical free space field representation in terms of the magnetic potential 𝐴₀𝑒 given in [2], [35] as

𝐸_𝑧0 = 𝑖𝜂₀𝑘₀(1 − 𝑠𝑖𝑛2_{𝛽𝑐𝑜𝑠}2_𝛼)𝐴 0

𝑒 _(2.25)

where the magnetic potential is represented in terms of the Fourier transform of the Gaussian aperture distribution function 𝑔̃̃(𝑘_𝑦0, 𝑘_𝑧0) as in [2],

15 𝐴₀𝑒(𝑥, 𝑦, 𝑧) = 𝐼0𝑙

𝑒𝑖𝑘0𝑅

4𝜋𝑅 𝑔̃̃(𝑘𝑦0, 𝑘𝑧0) ,

(2.25)

with 𝑅 the distance between the transmitter and the receiver.

From the results obtained, it can be clearly observed that the PE method provides consistent results up to 200 meters of lateral displacement from the main direction of propagation that corresponds to the maximum propagation angle 𝛽 of ~21 degrees that the PE can accurately model [2].

Figure 2.5 : Validation of the SSPE algorithm: Propagation in free space. 2.3.2 Propagation over a flat terrain

The second step was to add a perfectly conducting terrain over which the propagation takes place and to analyze its effect. The results are shown in Figure 2.6, where the propagation factor is shown for variable receiver heights at a distance of 1000 meters from the transmitter. The operating frequency is 1GHz and the transmitter height is of 30 meters. The source is Gaussian with characteristics shown in Figure 2.3 and Figure 2.4. The Fourier parameters are 𝑁𝑦 = 𝑁𝑧 = 1024. The grid displacement in the 𝑦 and 𝑧 directions is 𝛿𝑦 = 𝛿𝑧 = 𝜆. The results are compared to those obtained from a two-dimensional analysis of the same propagation space. The results shown in Figure 2.6 shown that identical results are obtained from both simulations.

16

Figure 2.6 : Propagation over perfectly conducting flat terrain. 2.3.3 Scattering from a PEC knife-edge

The third step was to include a single PEC knife-edge over the terrain and to observe the effect of the scattering of the propagating wave behind the edge. The knife-edge’s height and width are both 50 meters, its distance from the transmitter is 125 meters, and the receiver is placed at a distance of 500 meters from the transmitter and at a height 𝐻𝑡 = 60 𝑚𝑒𝑡𝑒𝑟𝑠. The geometry of the problem is shown in Figure 2.7.

Figure 2.7 : Geometry representing propagation over flat earth and an absorbing knife-edge.

The results shown in Figure 2.8 a) represent the vertical variations of the Propagation Factor (PF) for a receiver placed at 500 meters from the transmitter and a variable height from 0 to 80 meters. The results shown in Figure 2.8 b) represent the horizontal variations of the PF for the same transmitter and a receiver placed at a height of 2.1 meters while varying its horizontal position (y direction) between 0 and 150 meters.

17 a)

b)

Figure 2.8 : Propagation over flat earth and an absorbing knife-edge, a) Variations with respect to the vertical displacement of the receiver.b) Variations with respect to

later displacement of the receiver.

The propagation factor is computed by normalizing the field values obtained in the presence of the knife-edge, to the free space field values.

2.3.4 Scattering from multiple PEC knife-edge

In the following step, the effect of multiple absorbing knife-edges on the propagation of radio waves was considered. A 120 knife-edge with a separation distance of 50 meters on a total distance of 6 km were placed over the terrain as shown in Figure 2.9. The transmitter is positioned at a height of 125 meters. The source has the same

18

parameters as the case of a single knife-edge. The effect of the scattering of the wave due to the presence of these knife-edges on the propagation factor at a distance of 6000 meters and a variable height is shown in Figure 2.10.

Figure 2.9 : Multiple PEC knife-edge model.

Figure 2.10 : 120 knife-edges.Vertical variations of the Propagation Factor at a range x = 25 m, where reciever height of 0 coresponds to the kife-edge height. 2.3.5 Propagation over a hilly nonuniform terrain

As a next step, the propagation over a nonuniform terrain such as a hilly terrain is simulated by modeling the hills as successive absorbing knife-edges of different heights and width over a flat terrain. Two different geometries of hills have been considered: the first with a sinusoidal profile Figure 2.11 a) and the second with a trapezoidal profile Figure 2.12 b). The geometry of the hills is represented in two dimensions. Results for different hill widths are shown in Figure 2.11 b) and Figure 2.12 b).

19 a)

b)

Figure 2.11 : Propagation over a sinusoidal hilly terrain modeled by absorbing knife-edges: a) geometry b) numerical results with different values of the hill width. The effect of the width of the hill over the field values behind the knife-edge is clear and this is mainly due to the lateral waves propagating on the lateral sides of the hill. For hill width of 100 m, the hill can be considered as electrically infinite in the lateral direction, which permits it to be modeled in a 2-D space.

The geometries that have been considered show that 2-D simulations would not take into consideration the lateral waves diffracted by the vertical edges of the obstacles which leads to a misestimation of the field value beyond the obstacle. Nevertheless, the width of the hills does not seem to affect the value of the filed at a far distance, as can be seen from Figure 11 b) and Figure 12 b) where the field seems to converge to

20

the same value at 6 km beyond the hill. This implies that 2-D models would be quiet satisfying in such cases.

a)

b)

Figure 2.12 : Propagation over a trapezoidal hilly terrain modeled by absorbing knife-edges: a) geometry b) numerical results with different values of the hill width. On another hand, the computational cost of the method depends on the size of the obstacles and the size of the spectral sampling, which varies from geometry to geometry. The duration of the 3-D simulation is on average ~10 times that of 2-D.

21 2.4 Concusion

In this chapter, the split-step Fourier technique used to implement the Vector Parabolic Equation proposed by Janaswamy in [2] is introduced and its key formulation steps were presented. 2-D and 3-D simulations are compared for a set of models starting from free space to more complex terrain models including single and multiple knife-edges and hilly terrain. The main goal of these comparisons is to test the 3-D VPE algorithm that has been elaborated based on the given formulation and compared it to 2-D PE simulations. On an other hand, these comparisons give a preliminary idea about the geometries that require full 3-D simulation and those who can accurately be modeled by 2-D simulations. In general, the simulations and tests that are realized in this chapter allow to build and test on a simple level the 3-D VPE algorithm that is the backbone of the work in the following chapters.