EM Pulse Propagation Modeling for Tunnels by Three-Dimensional ADI-TDPE Method

The three-dimensional time domain parabolic equation (3D-TDPE) is derived to simulate the EM pulse propagation in straight, curved tunnels. Then the finite differential (FD) method is applied to accurately and flexibly describe the fine structures in the tunnel, such as cars. By using the alternating direction implicit (ADI) scheme to the spatial unknowns, the TDPE can be solved line by line in each transverse plane at any time step. More specifically, the computational efficiency can be improved significantly since a three-dimensional problem is changed into several one-dimensional problems. Moreover, both the Dirichlet and Neumann boundary conditions are used to estimate the transmission loss in tunnels. The rigorous numerical method, finite difference time domain (FDTD), is applied to verify the accuracy and efficiency of the proposed method. Numerical results are given to demonstrate that the proposed 3D-ADI-TDPE method can be used as an efficient tool to predict EM pulse propagation in tunnels including fine barriers.


I. INTRODUCTION
The accurate propagation modeling in tunnels for wide band plays an essential role in wireless communication systems. It is difficult to predict the EM pulse wave propagation in realistic environments since the rigorous solution of the Maxwell's equation is hard to be obtained for its insufferable computational resources. As a result, several asymptotic approaches were used to solve the wave propagation problems [1]- [3], namely geometric optics (GO), physical optics (PO), geometrical theory of diffraction (GTD) and empirical models. However, these techniques cannot provide accurate results since both the reflection and diffraction effects cannot be fully considered in the in-homogeneous atmosphere. The parabolic equation (PE), an approximation of wave equation, is extremely attractive for its stability and easily implementation. It should be noted that the effects of atmospheric refraction, diffraction, and reflection can be modeled by the PE method. However, it is time consuming to analyze broad-The associate editor coordinating the review of this manuscript and approving it for publication was Guido Lombardi . band EM pulse propagation problems with the frequencydomain methods. Therefore, the development of time-domain parabolic equation (TDPE) has a great significance.
Then it was introduced to analyze EM pulse propagation for irregular terrains [20], [21]. In the last decades, the 2D-TDPE was widely used to calculate the pulse propagation in the tunnels [22] and troposphere [23]. Based on the traditional 2D-TDPE method, some modified techniques are proposed to improve its accuracy, such as the two-way TDPE method [24], [25] and higher order TDPE method [26], [27]. More specifically, the two-way TDPE method can take both the forward and backward propagating effects into consideration. And the higher order TDPE method can provide accurate results at wider angles by employing high VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ order Padé expansion. All these studies were concerned on two-dimensional cases. More recently, Apaydin proposed a 3D finite-element and split-step Fourier-based two-way PE (3D-2W-SSPE) method to obtain the propagation in waveguide in frequency domain [28]. Then a hybrid TDPE/FDTD method was derived to model O2I radio wave propagation [29]. Moreover, we have investigated the 3D-TDPE to analyze the wide-band scattering from electrically large targets [30], [31]. All in all, the Fourier split-step (FSS) method is widely used to solve the propagation problems while the finite difference (FD) method is applied to the scattering problems. However, the FSS method cannot handle tunnels including barriers for its lack flexibility of complicated boundary modelling. In this paper, the 3D-TDPE is derived to model the EM pulse propagation in tunnels. Both the straight and curved tunnels can be taken into consideration by modifying the TDPE. The algorithm starts with applying FD method along the paraxial direction and time step, thus the calculation can be taken in a marching manner at each time step. Then the alternating direction implicit (ADI) scheme is used in each transverse plane. In this way, the calculation can be further accelerated since the unknowns can be calculated line by line. It should be noted that the FD method is used for better description of complicated boundary. At last, a series of numerical results are presented and compared with the finite difference time domain (FDTD) method.

A. STRAIGHT TUNNEL PROPAGATION MODEL OF 3D ADI-TDPE
The scalar three-dimensional wave equation can be derived as in which k is wave number, ψ is field component and n is the complex refractive index. By introducing the wave function u = e −ikx ψ(suppose the wave propagates along the paraxial direction x) to equation (1), then the following equation can be obtained The factorization method is applied in equation (2), then the wave equation can be rewritten as Then the forward parabolic wave can be obtained as Use the first order Taylor series to approximate the differential operator Q, which yields as Then the standard forward parabolic equation can be expressed as Define the Fourier transform (7) in whichF (k) = ∞ 0 E i e iωt dt is a spectrum function. s = ct − x describes the distance from the source to the observed point.
Using formula (7) to take the Fourier transform in (6) and consider free space satisfies n = 1, then the 3D TDPE can be rewritten as By using the central finite difference scheme, the 3D CN-TDPE can be derived as [30] 8 in which, m,l is the unknown transient fields of the mth transverse plane at lth time step, x, y, z are the spatial steps and ∇ y , ∇ z are the differential operators in space. At last, the 3D ADI-TDPE is used to enhance the computational efficiency. As shown in Fig. 1, the iteration diagram of the ADI scheme is given. The fields at the (m+1/2)th marching plane can be calculated column by column by the fields at the mth marching plane. Then the fields at the (m+1)th marching plane can be calculated row by row by the fields at the (m+1/2)th marching plane. By using the ADI scheme, the 3D ADI-TDPE can be derived as (11) in which, formula (10) represents the calculation from the mth transverse plane to (m +1/2)th plane, while formula (11) represents the calculation from the (m+1/2)th transverse plane to (m +1)th plane. p,q m,l is unknown transient field of mth transverse plane at lth time step at y = p y, z = q z.

B. CURVED TUNNEL PROPAGATION MODEL OF 3D ADI-TDPE
Actually, the tunnel is curved in real life. According to [13], the 3D ADI-FDPE in curved tunnels can be written as in which, a = i x 4k y 2 , b = i x 4k z 2 , δ y , δ z is the differential operator in y, z direction. A n p , B n p describe the curved tunnel and their expressions are given as It should be noted that when dealing with formula (12), an approximation is introduced. The formula (12) could be expanded as in which Assuming that the realistic curvature radius of tunnel is usually large enough. Thus the (q yk x) 2 can be ignored. Then the formula (17) can be approximated to Supposing the correction factors A, B, C are set as Then using the Fourier transform to equation (12) and (13), the equation can be rewritten in the following matrix form.
The boundary conditions of the tunnels are set as two types, namely Dirichlet boundary condition and Neumann boundary condition. The expressions for the Dirichlet and Neumann boundary conditions in time domain are written as follows = 0 · · · · · · Dirichlet ∂ /∂n = 0 · · · · · · Neumann (22) It should be noted that the boundary conditions of vehicles in tunnels are set as Dirichlet boundary condition. Moreover, the proposed method can be extended to realistic lossy tunnel or obstacle boundaries by using the impedance boundary conditions.

D. EXCITATION SOURCE
In this paper, all examples use modulated Gaussian pulse as excitation source. Assume that the center frequency of Gaussian pulse is f (GHz) and pulse width is τ (ns), the expression of the excitation source is given in which η = τ 3 and the unit of η, t are nanosecond.
However, the disperse in time-domain must satisfy Courant-Friedrichs-Lewy (CFL) condition, that is The source is located at the first transverse plane and boundary conditions are added at the surface of automobile and tunnels. It should be noted that the mesh size is usually set as one-tenth of a wavelength or less. As a result, the outline of the automobiles can be described with higher accuracy. By using the Fast Fourier Transform (FFT) algorithm, the field components distributions can be obtained within a certain bandwidth. However, the proposed 3D ADI-TDPE method can only model the forward propagation problems for electrically large convex objects.
in which k x = mπ x, k y = pπ y, k z = qπ z. Then substitute equation (25) into the former marching equation (20) with x = y = z. Moreover, suppose the electromagnetic wave mainly propagates along the x direction, such that k x ≈ k, k y ≈ 0, k z ≈ 0. Then equation (20) can be rewritten as The growth factor can be obtained Therefore, the proposed 3D ADI-TDPE for the former marching equation is unconditionally stable. Similarly, the latter marching equation can be proved to be unconditionally stable. To get the numerical dispersion error of the proposed method, the time-harmonic electric field can be expressed as E (x, y, z, t) = Re E 0 e −i ωt−k x x−k y y−k z z (28) By using the relationship between E and , substitute equation (28) into the former marching equation (20), then the following equation can be obtained cos χ +k s+ k −k x (1/2) x 1+ x s 2 y 2 sin 2k y y 2 x s 2 y 2 sin 2k y y 2 (29) in which χ = kl s + k −k x m x −k y p y −k z q z. Let k x ≈k cos θ ,k y ≈k sin θ cos ϕ,k z ≈k sin θ sin ϕ, then the equation (29) can be simplified to Assume that σ = k −k , in which k is the real wave number whilek is the calculated wave number. By ignoring the high-order error term of σ 2 , equation (30) can be rewritten as Similarly, the dispersion error of the latter marching equation can be obtained as It can be concluded that the dispersion error of the proposed method is approximate to zero when the wave propagates along the paraxial direction (θ ≈ 0).

III. NUMERICAL EXAMPLES
The computer of Intel Xeon E7-4850 CPU equipped with 8GB RAM is used to simulate all the numerical results. For the tunnel models, the modulated Gaussian pulse is used as the incident plane wave and it can be defined as where y = y − y 0 and z = z − z 0 . η y and η z are the parameter to limit the energy of Gaussian distribution, and can be written as η y = σ 2 y ln 10, η z = σ 2 z ln 10. Firstly, the time-domain response of a rectangular waveguide is simulated to validate the accuracy and efficiency of the proposed method. The geometry of the rectangular waveguide is shown in Fig. 2. Wave port is used as the excitation, and only the TE10 is considered. Meanwhile the observation point is set at (20m, 0m, 1.5m). As shown in Fig. 3, both the FDTD and ADI-TDPE methods are used to simulate the timedomain response of the rectangular waveguide. Here, the FDTD results are obtained by our own coding. FDTD is a fullwave numerical method of wave equation, which can be used as the benchmark. In this paper, the comparisons are made between the FDTD and the proposed 3D ADI-TDPE methods to demonstrate the accuracy of the proposed method. On the  other hand, the parabolic equation is an approximation of the wave equation, which can be solve in a marching manner along the paraxial direction. As a result, a 3D problem can be converted to a series of 2D problems along the paraxial direction and the computational resources can be reduced significantly. By using the ADI scheme, the unknowns in each transverse plane can be calculated line by line. In other words, a 3D problem can be converted to a series of 1D problems.    Therefore, both the CPU time and memory requirement can be saved significantly. There are three different mesh sizes are used for the FDTD method, namely 0.1m, 0.05m and  0.02m. More specifically, the results can be convergent for the mesh size of 0.02m. It can be seen that the proposed TDPE method with the mesh size of 0.1m can agree with the results of FDTD method for 0.02m. Since the dispersion error of the proposed 3D ADI-TDPE method is small enough when the wave propagated along the paraxial direction, the convergence of the time-domain response for the proposed method performs better than the FDTD method. Moreover, the computing resources are compared between these two methods. As shown in Table 1, both the CPU time and the memory requirement can be reduced significantly for the proposed ADI-TDPE method when compared with the FDTD method. Since the unknowns can be calculated line by line in each transverse plane at any time step by the proposed method, the proposed ADI-TDPE will achieve higher efficiency than the FDTD method. The time-domain reponse of magnetic field can be obtained by the solving the differential form of Maxwells' equation. As shown in Fig.4, the time domain response of magnetic field calculated by both the FDTD and VOLUME 8, 2020   in which ylength is the length of waveguide in y direction and its value is 4m. The range of y is [−2m, 2m]. Secondly, we consider a rectangular tunnel with the crosssection size of 6.0m×4.0m. A Gaussian source of excitation is set at the center of the first plane. The parameters in formula (33) are set as η y = 1.2m, η z = 0.8m, y 0 = 0m and z 0 = 2m. The center frequency of Gaussian pulse is 0.8GHz and pulse width is 5ns. Both the Dirichlet and Neumann boundary conditions are used. Fig.5 shows the field distributions at the frequencies of 462MHz, 698MHz and 1.2GHz for two kinds of boundary conditions.
In addition, we compare the field distribution between straight and curved tunnel based on the previous tunnel model. The model parameter is set similarly. Both the straight and curved cases are simulated by the proposed method and the curvature radius is set as 1000m. The center frequency of Gaussian pulse is 0.5GHz and pulse width is 10ns. As shown in Fig. 6, the time-domain responses of the straight and curved rectangular tunnels are compared at different distances. It can be seen that the time-domain response will change more obviously with the observation point locates farther. Moreover, the field components distributions of yoz plane at different VOLUME 8, 2020 frequencies are given in Fig. 7 and Fig. 8 for straight and curved cases in two boundary conditions. At last, the propagation factors (PF) of the proposed method at the rage of 50m are compared with the frequency-domain PE (FDPE) method in Fig. 9 (in Dirichlet boundary) and Fig.10 (in Neumann boundary). It can be seen that there is a good agreement between them. More specifically, the calculation should be taken for FDPE method at different sample frequencies.
As a result, the computational time will increase largely for wide band problems. It can be concluded that the proposed ADI-TDPE method can be used as an efficient tool to simulate EM pulse propagation in tunnels for wide band.
The propagation factor is calculated as follow PF = 20 log 10 (|ψ|) + 10 log 10 (d) in which ψ is the field component. d means the distance from the source to the observed cross section. At last, a more practical and complex tunnel is modeled. As shown in Fig. 11, a vaulted tunnel of Dirichlet boundary condition with three different kinds of conducting vehicles is analyzed. The curvature radius is set to be 600m. Meanwhile, Fig.12 shows the real sizes of three automobile models. The time-domain responses for tunnels with and without automobiles are given at different distances in Fig.13. Meanwhile, as shown in Fig.14, the field components distributions of yoz plane at different frequencies are given for tunnels with and without automobiles. The parameter of Gaussian source are set as η y = 2.5m, η z = 1.5m, y 0 = 0m and z 0 = 3m. It can be seen that the multipath effect of the wave propagation will be more obvious when vehicles exist in tunnel. In addition, there will be some differences in the time-domain far-field response. Influenced by the bending of the tunnel, the phenomenon of field accumulation occurs at the distances of 25m (see Fig. 13 (a)-(c)) and 50m (see Fig. 13 (d)-(f)). The wave reflection of the vehicles causes some ripples in the field distribution diagram. Since the observation position is exactly where the car is, the field distribution is composed of two parts (see Fig. 13 (g)-(i)). The first part is the field distribution in the car, where fields are equal to zero. For the rest of the region, the wave reflection is obvious near the cars. The fields will be newly accumulated at the distances of 55m (see Fig. 13 (j)-(l)). It can be seen that the proposed 3D ADI-TDPE method can be used to predict the time-domain response of the wave propagation in the complex curved tunnels with vehicles. However, the multipath effects can not be modeled by the proposed one-way PE method. It should be noted that the two-way parabolic equation (2W-PE) can be used to estimate the multi-path propagation with high accuracy, which can be considered in the future work.

IV. CONCLUSION
In this paper, a 3D ADI-TDPE method is derived to model the EM pulse propagation in tunnels with vehicles. By taking advantage of the ADI scheme, the complicated boundary of vehicles can be modeled more accurately and flexibly than the FSS method. Moreover, the ADI-TDPE method is extremely attractive for its high computational efficiency. Some numerical results for different boundary conditions are simulated to predict the EM pulse propagation in tunnels. It can be concluded that the proposed method has a good performance to model complex and practical tunnel environments.