Comparison of FDTD and TDWP Methods for Simulating Electromagnetic Wave Propagation above Terrains

Brazilian Microwave and Optoelectronics Society-SBMO received 26 Jan 2018; for review 09 Feb 2018; accepted 16Aug 2018 Brazilian Society of Electromagnetism-SBMag © 2018 SBMO/SBMag ISSN 2179-1074 Abstract— Researchers have always been in face of path loss calculation in different media for applications such as telecommunication link design. Wave propagation calculations in large spaces using the FDTD method is time-consuming and imposes a great computational burden. For this reason, to replace the classical FDTD method for wave propagation simulation and path loss calculation in large spaces, optimized methods, namely TDWP, have been provided. In this paper, the use of the TDWP method for wave propagation simulation and path loss calculation above a terrain is investigated. Longitudinal components of ground-waves are taken into account (including direct waves, ground reflections, and surface waves). Propagation space is longitudinally divided into smaller FDTD windows with finite length. The electromagnetic pulse travels through these windows from left to right to the desired point. But despite its capability in reducing computational burden and increasing processing speed, TDWP has lower precision in instantaneous field simulation and calculation of propagation coefficients, so that results obtained from FDTD and TDWP are clearly different. In this paper, some efficient methods are proposed, which yield an increase in method accuracy.


I. INTRODUCTION
Modern communication and radar systems need effective tools to foresee propagation.In many cellular communication and HF and VHF radar sites, high precision propagation simulation is vital.
The results are necessary for choosing proper places for ground stations and optimum height for transmitters and receivers, and also reducing the costs of implementing communication systems.One of the primary aims of propagation simulation is to forecast path loss ( ) with high precision.Path loss is defined as: Comparison of FDTD and TDWP Methods for Simulating Electromagnetic Wave Propagation above Terrains Mahmoud Fallah, Amir.H Nazeri, School of Electrical Engineering, Iran University of Science and Technology, Tehran, Iran performed using both TWDP and FDTD, and the results are compared.Finally, the simulation process is modified to reduce these differences.

II. TDWP METHOD
The wave propagation problem may be set in city streets with numerous buildings, in a vast desert, or in a mountainous area.The classical FDTD method faces shortcomings when simulating wave propagation along long paths.Suppose we would like to study wave front propagation with spatial Gaussian distribution (along the height axis) in a 10 km path and a height of 50 m.If we divide this space into a mesh of 10 cm cubic cells, there would be 5000000 cells in this space, not considering PML blocks located at the left, right, and top of the computational space.The propagated wave consists of three components, regardless of its polarization.Every time these three components are updated, the new values should be stored in matrices.So we need three matrices with 5 milli on elements each, all with 15 million elements updated in every time step.Considering the fact that the utilized program uses some of the computer's RAM and CPU cycles, it can be concluded that processing such huge amounts of data with personal computers is impossible.Even if sufficient RAM capacity is obtained, processing the data takes considerable time.
To overcome this problem, instead of updating all components for every computational cell at every moment, a rectangular window around the pulse can be presumed, containing most of the propagated pulse energy.In this manner, the computational space can be divided into smaller windows, as shown in fig. 1 and 2. These smaller rectangular windows must be surrounded by PML from the left, right, and top.PML blocks surrounding the computational space prevent reflection and refraction of the waves along the boundaries.The aim is to compute propagation loss along the wave's path and also along the height axis.To this end, we trigger the computational space with a short duration pulse.The trigger source is located at the beginning of the first window.When it begins to stimulate field components in the In fact, in this manner, a rectangular computational space with very great length is divided into rectangular computational spaces with smaller length, but the same height as the huge rectangle, so that the simulation can be readily done with a personal computer.

A. Wave Propagation Formulation
Although the physics of propagation is taking place in a three-dimensional space, for a symmetrical problem with respect to one of the Cartesian coordinates, we can use a two-dimensional FDTD approximation.Maxwell's equations and the constitutive parameters, in the form of vector differentials, in an isotropic space containing magnetic and dielectric materials, are as follows:

(n+1)th window nth window
Cells of the first column of the (n+1)th window, where field components of the last column of the preceding window are restored.
Cells of the last column of the nth window, used to store field components .
. ( 4 )   If a node of the magnetic field is placed within the setting of the conductor, it is assumed that the entire cell containing that node is placed inside the conductor.Thus, all electric field components of that cell (which include that node of the magnetic field) are fixed to zero.

IV. PATH LOSS CALCULATION
In this problem, the propagation space is a semi-confined space, limited by and vertically, and by and in the plane parallel to the ground.We would like to calculate the loss caused by the propagation medium, without considering the effect of the receiver and transmitter antenna.In other words, we would like to calculate the power lost in propagation, after the wave is transmitted by the transmitter antenna, before it is absorbed by the receiving antenna.
Therefore, in (10), , , and should be substituted with 1.In every given point, can be calculated in terms of field magnitude [2]: (10) Therefore, to compute the propagation loss, we need to calculate the electric field magnitude at every given point in the propagation space.
To be able to compare the results for different media, we need to eliminate the effect of the reduction of wave power due to the distance from the source via the inverse square law, which is the only loss associated with a lossless medium without obstacles.The propagation factor at every given point is defined as (for any specific frequency).Where and are the measured and reference electric field at the desired point, respectively [2].The reference field is the field propagated from the same source, not taking the terrain into account (fig.5).When the wave reaches the boundaries of the propagation space into free space, it must be ensured that it leaves the studied medium and will not be reflected back.To simulate the free space surrounding the propagation space, electromagnetic absorber layers must be used.These layers (fig.5) attenuate the wave, so that it doesn't return to the propagation medium after leaving it.One of the most appropriate types of such absorber layers, for simulations, is the "Perfectly Matched Layer" or PML for short [11].So in this method, we first perform the simulation for the terrain and store the derived values for both a specific height and different distances, or also a specific distance and different heights.In this step, the propagation medium is surrounded by the PML from the top, left, and right.In the next step, the simulation is repeated with the terrain removed and the propagation space covered by PML from the bottom as well.In this step, the values are stored in the same cells as the last simulation [12].For every scenario, two simulation steps should be performed.In the first, the fields are calculated in the presence of the terrain.In the second, the fields are calculated without regarding the terrain.
Finally, after both steps are done, by performing DFT analysis on the time-domain data available (in the desired bandwidth), for each desired frequency is derived.Usually, it would be favorable to calculate and plot this ratio in decibel units (2* ) for a specific distance or height.

A. Source Stimulation
To stimulate the fields of the computational space, we need to model a physical source.Because finite time distributions, which also have limited frequency bandwidth, provide better responses in FDTD simulations, the pulse shape (in the time domain) is usually selected as Gaussian.However, here we prefer to use the derivative of the Gaussian function, in order to prevent DC effects that can cause shifting of the wave propagation frequency in some applications.However, we should note that the desired pulse for our source should cover the range of frequencies determined by the bandwidth.The Gaussian function derivative used in our algorithm is as follows: Where  is the time shift of the function toward the positive direction of the time axis; which is equal Now which i and j the simulation applies to depend on the type of stimulation.The distribution, e.g.pulse time distribution, can also be in different forms.One of the most popular spatial distributions along the vertical axis is the Gaussian distribution; where its centricity is shown with source _ h .The statistical process of spatially distributing the pulse along the z axis is achieved by its variance.As it can be seen in frames 3 and 4 of fig. 9 (in the TDWP method), when the propagated wave reaches the end of the first window, it is restored in the beginning of the second window.But, as seen in frames 6 of fig.9, reflections from the right side triangle (situated in the second window) are  the results derived from FDTD [5].Fig 12, shows the comparison between the propagation loss results of two method, Optimized TDWP and FDTD, in Lines 3 and 4.

C. Analysis of Improving Factors of Optimized TDWP Method
In following some relations are defined which can help us to make a suitable analysis regarding the efficiency of proposed Technics, in order to improve traditional TDWP method.If we denote the propagation loss factor of Traditional and our Optimized TDWP method as ) (op TDW P PF and ) (re TDW P PF , respectively, the amount of generated Error in terms of y (axis) will be calculated for Traditional and Optimized methods as below, respectively.It should be emphasized that the one-window FDTD method is not applicable for long pathes.It is applied to this study just for a path of 100 m.

Fig 1 .
Fig 1. Dividing natural terrain into smaller computational regions

Fig 3 .
Fig 3. Storing and restoring field components in two consequent windows

. 4 ,
the configuration of the electric and magnetic field nodes are shown for TM polarization.In this figure, fields components within closed dash lines have the same index in the implementation of the algorithm.The nodes of the electric field, with lines to illustrate the direction of field components, are shown with solid circles; and the nodes of the magnetic field are shown with empty squares.Triangle shaped dashed lines represent a group of nodes that have identical array indices in the computer algorithm.For example, in the bottom left of the grid, all three nodes inside the triangle have diagram on the bottom left of the grid shows nodes related to each field and their indices for a typical group.In the diagram on the right side of the grid, the same nodes and their array indices have been shown when used in a computer program.In this polarization, none of the nodes of the fields were located in places with the correct values of discretized space.

Fig 4 .
Fig 4. Configuration of the nodes of electric and magnetic fields for TM polarization used in simulation

Fig 5 .
Fig 5.For every scenario, two simulation steps should be performed.In the first, the fields are calculated in the presence of the terrain.In the second, the fields are calculated without regarding the terrain.

to 4 Fig 6 .
Fig 6. a) The diagram of the Gaussian pulse derivative in terms of time for the bandwidth 200MHz , b) The diagram of the Fourier transform of the derivative Gaussian pulse.
of the spatial distribution center for the source from the ground and v  is the spatial variance of the source.In fig.7, the spatial distribution of the source with centricity of m source h 25 _  from the ground at a height of 50 m is shown for different variances.

Fig 7 .
Fig 7. Spatial distribution of the source with centricity of 0 25 Zm  from the ground at a height

Fig 9 .
Fig 9. Instantaneous images captured for comparing the FDTD and TDWP methods in the setup shown in fig.8

Fig 10 .
Fig 10.Comparison between propagation losses for lines , , , and , obtained from FDTD and TWDP

Fig11.
Fig11.Comparison of the instantaneous images taken, with field regeneration in the first column of the second window cells a) not limited b) limited.

Fig 12 .
Fig 12.Comparison of the propagation loss calculated using FDTD and modified TWDP for lines and , with field regeneration duration limited in the second window.
a Improving Factor for each height of y. improving factor of proposed method in one given line (like , we should first calculate the below Integrals for each line.
re to op  is defined as a Improving Factor for each vertical line x.similar way, for calculating the propagation loss factor along a horizontal line (like as improving factor for each length of x. improving factor of proposed method in one given line (like , we should first calculate the below Integrals for each line. op to re  is defined as a Improving Factor for each vertical line x.for attaining to a total improving factor in one window, relation 16 (or 20) should be employed over a given surface like