Developing a leap-frog meshless methods with radial basis functions for modeling of electromagnetic concentrator

1.   Introduction

The research on electromagnetic metamaterials is a long-term dream of mankind and has potential in many interesting applications, which mainly include invisibility cloak, telecommunications to sensing, sub-wavelength imaging and radar technology. We can refer to ^[18,27] and their references. Before transformation optics ^[27,34] technology and metamaterials ^[5,6,7] were proposed, these were unimaginable. Therefore, transformation optics technology provides a shortcut for designing interesting electromagnetic devices with metamaterials. Because of the difficulty and large expense of repeating physical experiments to test the performance of specific electromagnetic devices, the numerical simulation of electromagnetic metamaterials occupies a very important position in the field of metamaterial research. In recent years, some researchers have proposed some effective and accurate numerical calculation methods in the field of metamaterials to solve the time-domain Maxwell equations for the metamaterials. The main numerical calculation methods are the meshless method ^[2,16,26,29], the finite-difference time-domain method ^[1,40] and the time-domain finite element method ^{[8,24,25,36,37]}. These are also several popular numerical methods for solving partial differential equation.

The numerical method of the RBF was first proposed by Kansa in 1990 ^[23]. It is a new numerical method for solving partial differential equation. This new method does not need to generate a specific grid, as one can just generate scatter points in the calculation area, so the program implementation is simpler than that for the traditional finite element method and finite difference method. The finite difference method has relatively high requirements for regular regions, and it is difficult to solve problems on complex regions. In the numerical calculation of the finite element method, the area to be solved is first divided into triangles or quadrilaterals, so as to form elements and meshes, which makes mesh division and processing to take up a lot of time. For example, the meshing of complex three-dimensional structures is more difficult, and the meshing of special structures is more difficult. Therefore, we need to study the gridless method, the key is the generation of grid points, which can overcome the dependence of traditional numerical methods on grids. But the RBF method is not theoretically justified in terms of stability and convergence, and it is not robust in terms of those free parameters of the RBF.

Since the introduction of the meshless RBF method, people have paid more and more attention to using the RBF to solve various partial differential equation over various domains ^{[4,19,20,21,30,31,32,39]}, including nonlinear Schr$\ddot{o}$dinger equation ^[11], nonlinear Green-Naghdi equation ^[10], shallow water equation ^[12], magnetohydrodynamics equation ^[9]. For more detailed information on the RBF meshless method, you can refer to articles ^[15], monographs ^[2,16,26]. At the same time, there are some good papers about the application of the radial basis functions method to study simple media using Maxwell equations that change with time ^{[13,14,29,35]}. Recently, Li and Nan used the leap-frog RBF meshless method to simulate the backward wave propagation phenomenon by solving the complex Maxwell equations in metamaterials ^[29]. Inspired by the success of the previous paper, my main goal is to further extend this RBF method to solve the more complex Maxwell equations and simulate the wave propagation phenomenon in electromagnetic wave concentrator devices designed by transformation optical technology. As far as I know, this is the first work on the radial basis functions method of electromagnetic concentrator devices.

The remaining chapters of this paper are arranged as follows. In Section 2, I introduce the governing equations for the electromagnetic concentrator model and its leap-frog scheme. In order to simulate the wave propagation behavior of the electromagnetic concentrator, I introduce the governing equations of a perfect matching layer model in Section 3. In addition, to practically implement the RBF algorithm, I then rewrite both the concentrator model and the PML model into just one system of PDEs in Section 3. In Section 4, I validate the effectiveness of the meshless RBF method for this complex model equation through numerical tests. Finally, I have summarized this paper in Section 5.

2.   Governing equations for the electromagnetic concentrator model and its leap-frog scheme

In order to better understand the propagation behavior of electromagnetic waves in various complex media, we need to numerically solve the following well-known Maxwell equations:

where $\boldsymbol{D}, \boldsymbol{E}, \boldsymbol{B}, \boldsymbol{H}$ represent the electric flux density, electric field, magnetic flux density, magnetic field, respectively. $\boldsymbol{J}_{s}, \rho$ represent the electric current and electric charge density, respectively. In order to ensure Maxwell equations are well posed, we need to couple the Maxwell equations with the following constitutive equations for complex media:

where $\epsilon_{0}$ and $\mu_{0}$ are the permittivity and permeability in vacuum, and $\epsilon$ and $\mu$ are the relative permittivity and permeability of the specific medium, respectively.

For cylindrical electromagnetic concentrators, the ideal material parameters for the polar coordinate system can be derived from the following corresponding coordinate transformations ^[33]:

This transformation compresses the circle $r\leq R_{2}$ into a smaller circle $r'\leq R_{1}$ and stretches the ring $R_{2}\leq r\leq R_{3}$ to $R_{1}\leq r'\leq R_{3}$. Because stretching and compression are mutually compensated, it leads to consistency between the original space and the transformed space. The purpose of this transformation is to concentrate energy in the internal region $r'\leq R_{1}$. In order to better understand the function of the electromagnetic device, a structural diagram corresponding to the above coordinate transformation (2.3) is shown in Figure 1.

According to the coordinate transformation (2.3), using the form invariance of Maxwell's equation, that is, the principle of transformation optics, we can derive the following complex time-dependent governing equations to simulate the wave propagation in the electromagnetic concentrator device is Eq (2.17) on page 141 of ^[36], rewrite as follows:

In the above set of equations, $\epsilon_{r} = \frac{r'+K_{1}}{r'}$, where $K_{1} = \frac{(R_{2}-R_{1})R_{3}}{R_{3}-R_{2}}$ and $r'\in [R_{1}, R_{3}]$. The parameters $\epsilon_{\max} = \frac{R_{3}-R_{2}}{R_{3}-R_{1}}, \mu_{\max} = \frac{R_{2}(R_{3}-R_{2})}{R_{1}(R_{3}-R_{1})}, \omega_{e}, \; \omega_{m}$ are used in the following Drude model:

The coefficient matrices $M_{ca}$ and $M_{cb}$ in (2.4) are as follows

The stability of the model equation (2.4) in the continuous case can be found in Reference ^[22].

Next, I construct the RBF meshless method for Eq (2.4). In order to better develop numerical format, We can re-express $M_{ca}$ and $M_{cb}$ as follows:

and re-express the model equation (2.4) in the following scalar form:

To solve this problem numerically, I assume that the boundary condition of (2.5) is the perfect conductor boundary:

where $\boldsymbol{n}$ is the outward normal vector of units on $\partial\Omega$. For the two-dimensional rectangular domain $\Omega = [a, b]\times[c, d]$, (2.6) is simplified to

For the discrete model equation (2.5) in domain $\Omega$, assuming that the region is $\Omega = [a, b]\times[c, d]$, I cover the domain $\bar{\Omega}$ by $N$ scattered points $Q_{i}: = (x_{i}, y_{i}), i = 1, \cdots, N$, where $\{Q_{i}\}_{i = 1}^{N_{I}}$ represents the internal points and $\{Q_{i}\}_{N_{I}+1}^{N}$ represents boundary points. Furthermore, we divide the time interval $[0, T]$ into $N_{t}$ uniform intervals, i.e., we have discrete times $t_{i} = i\tau, i = 0, 1, iii, N_{t}$, where the time step size $\tau = T /N_{t}$.

I can construct the following RBF meshless method: Given the proper initial conditions $E_{x}^{-\frac{3}{2}},$ $E_{x}^{-\frac{1}{2}}$, $E_{y}^{-\frac{3}{2}}, E_{y}^{-\frac{1}{2}}, D_{x}^{-\frac{3}{2}}, D_{x}^{-\frac{1}{2}}, D_{y}^{-\frac{3}{2}}, D_{y}^{-\frac{1}{2}},$ $H^{0}$, $H^{-1}$, at all interior points, for any $n\geq 0$, find $D_{x}^{n+\frac{1}{2}}$, $D_{y}^{n+\frac{1}{2}},$ $E_{x}^{n+\frac{1}{2}},$ $E_{y}^{n+\frac{1}{2}}, B^{n+1}, H^{n+1}$ such that

Regarding how to approximate the spatial derivatives $\frac{\partial H^{n}}{\partial y}, \frac{\partial H^{n}}{\partial x}, \frac{\partial E_{x}^{n+\frac{1}{2}}}{\partial y}, \frac{\partial E_{y}^{n+\frac{1}{2}}}{\partial x}$, I will discuss this in detail in Section 3.

3.   Governing equations of a perfect matching layer model and the RBF meshless method for coupled model

In order to simulate the wave propagation phenomenon in the electromagnetic concentrator, I need to introduce PML to simplify the unbounded region problem to a bounded region problem in order to absorb the external electromagnetic wave. Moreover, the PML is used to truncate the unbounded wave propagation problem to a bounded domain problem. In this paper, I use the 2-D transverse electric (TE) Ziolkowski PML model, which is complementary to the 2-D transverse magnetic Ziolkowski PML model given in ^[28]. Fro more detailed, I presente the expressions for the TE PML modeling equations we need, as follows:

Here $\boldsymbol{E}, H, \epsilon_{0}, \mu_{0}$ represent the same meaning as that in model equation (2.5). $\boldsymbol{J} = (J_{x}, J_{y})$ and $K$ represent the induced current and the induced magnetic current, respectively. Moreover $M_{1}$, $M_{2}$ and $M_{3}$ are $2\times2$ diagonal matrices $M_{1} = diag(\sigma_{y}-\sigma_{x}, \sigma_{x}-\sigma_{y})$, $M_{2} = diag(\sigma_{x}, \sigma_{y})$ and $M_{3} = diag((\sigma_{x}-\sigma_{y})\sigma_{x}, (\sigma_{y}-\sigma_{x})\sigma_{y})$, respectively, where $\sigma_{x}$ and $\sigma_{y}$ are the electrical conductivities in the $x$ and $y$ directions, respectively. The specific expressions for $\sigma_{x}$ and $\sigma_{y}$ are given later and they are non-negative functions. In order to numerically solve the PML governing equation (3.1) using the method in this paper, I rewrite it into the scalar form as follows:

To maintain the same discrete format as (2.8), I propose the following RBF method for the scalar form of the PML model equation (3.2). Given the proper initial conditions $E_{x}^{-\frac{1}{2}}, E_{y}^{-\frac{1}{2}}, H^{0}, J_{x}^{0}, J_{y}^{0}, K^{-\frac{1}{2}}$, for any $n\geq 0$, find $E_{x}^{n+\frac{1}{2}}, E_{y}^{n+\frac{1}{2}}, H^{n+1}$, $J_{x}^{n+1}, J_{y}^{n+1}, K^{n+\frac{1}{2}}$ such that

In the actual numerical implementation process, I need to couple the concentrator model and the PML model into a new system, and then apply the RBF method for this new system. The following RBF meshless method is obtained: at all configuration points $Q_{i}$ and all $n\in[0, N_{t}-1]$:

The corresponding coefficient expressions are as follows:

In order to numerically solve the governing equations (3.4)–(3.10) by the RBF meshless method, I represent the unknown solution functions $D_{x}^{n+\frac{1}{2}}, D_{y}^{n+\frac{1}{2}}, E_{x}^{n+\frac{1}{2}}, E_{y}^{n+\frac{1}{2}}, K^{n+\frac{1}{2}}, J_{x}^{n+1}, J_{y}^{n+1}, H^{n+1}$ as expansions of RBFs:

where $D_{xj}^{n+\frac{1}{2}}, \; D_{yj}^{n+\frac{1}{2}}, \; E_{xj}^{n+\frac{1}{2}}, \; E_{yj}^{n+\frac{1}{2}}, \; K_{j}^{n+\frac{1}{2}}, \; J_{xj}^{n+1}, \; J_{yj}^{n+1}$ and $H_{j}^{n+1}$ are the unknown coefficients to be determined. And the RBF is $\varphi_{j}(x, y): = \varphi(\sqrt{(x-x_{j})^{2}+(y-y_{j})^{2}})$ at the node $Q_{j}, j = 1, \cdots, N$. Note that I have a wide variety of choices for RBF. In the numerical tests presented in this paper, I use two basis functions, first, I focus on the globally smooth multiquadric (MQ) RBF:

where $\gamma$ is an odd integer and $\delta > 0$ is a free parameter. Then we consider the Gaussian RBF:

where $s > 0$ is a free parameter. The selection of the above parameters $\delta$, $\gamma$ and $s$ has a very large effect on the numerical simulation. Here, I want to point out that how to select an appropriate shape parameter is very challenging because it is very sensitive to the influence of the numerical results, although there are some good works in the literature ^[3,17,38]. For convenience, I denote $\partial_{y}H$ for the partial derivative of $H$ with respective to $y$. Similar symbols are used for the other partial derivatives. Now we can use the basis function to represent the partial derivative terms $\partial_{y}H^{n}, \; \partial_{x}H^{n}, \; \partial_{x}E_{y}^{n+\frac{1}{2}}, \; \partial_{y}E_{x}^{n+\frac{1}{2}}$ in (3.4)–(3.10).

The detailed implementation process for using the RBF meshless method to solve the governing equations (3.4)–(3.10) is as follows:

$\bf{Step\; 1.}$ First use the initial value $H^{0}$ to evaluate the coefficients $H_{j}^{0}$ by solving a linear system, then by using $(3.14)$ to calculate $\partial_{y} H^{0}$ and $\partial_{x} H^{0}$, final substitute the initial values $D_{x}^{-\frac{1}{2}}$, $D_{y}^{-\frac{1}{2}}$ and $\partial_{y} H^{0}$, $\partial_{x} H^{0}$ into $(3.4)$ to work out $D_{x}^{\frac{1}{2}}$ and $D_{y}^{\frac{1}{2}}$.

$\bf{Step\; 2.}$ Substitute the initial values $D_{x}^{-\frac{3}{2}}$, $D_{x}^{-\frac{1}{2}}$, $D_{y}^{-\frac{3}{2}}$, $D_{y}^{-\frac{1}{2}}$, $E_{x}^{-\frac{3}{2}}$, $E_{x}^{-\frac{1}{2}}$, $E_{y}^{-\frac{3}{2}}$, $E_{y}^{-\frac{1}{2}}$, $J_{x}^{0}$, $J_{y}^{0}$, as well as $\partial_{y} H^{0}$, $\partial_{x} H^{0}$, $D_{x}^{\frac{1}{2}}$ and $D_{y}^{\frac{1}{2}}$ in $\bf{Step\; 1}$ into $(3.5)$ and $(3.6)$ yields $E_{x}^{\frac{1}{2}}$ and $E_{y}^{\frac{1}{2}}$.

$\bf{Step\; 3.}$ Apply initial values $H^{0}$ and $K^{-\frac{1}{2}}$ to work out $K^{\frac{1}{2}}$ in $(3.7)$.

$\bf{Step\; 4.}$ Apply the obtained $E_{x}^{\frac{1}{2}}$ and $E_{y}^{\frac{1}{2}}$ in $\bf{Step\; 2}$ and the given initial values $J_{x}^{0}$, $J_{y}^{0}$ in $(3.8)$ and $(3.9)$, to respectively obtain $J_{x}^{1}$, $J_{y}^{1}$.

$\bf{Step\; 5.}$ Apply the obtained $E_{x}^{\frac{1}{2}}$ and $E_{y}^{\frac{1}{2}}$ in $\bf{Step\; 2}$ to obtain the coefficients $E_{xj}^{\frac{1}{2}}, E_{yj}^{\frac{1}{2}}$ by solving the same linear system as described in $\bf{Step\; 1}$ to respectively calculate $\partial_{x}E_{y}^{\frac{1}{2}}$ and $\partial_{y}E_{x}^{\frac{1}{2}}$, then substitute the obtained $K^{\frac{1}{2}}$ in $\bf{Step\; 3}$ and the given initialvalues $H^{0}, H^{-1}, B^{0}, B^{-1}$ and the just calculated $\partial_{x}E_{y}^{\frac{1}{2}}$ and $\partial_{y}E_{x}^{\frac{1}{2}}$ into $(3.10)$ to get $H^{1}$.

$\bf{Step\; 6.}$ Repeat the above steps $\bf{Step\; 1}-\bf{Step\; 5}$ until the last step.

4.   Numerical results

In this section, I will verify the effectiveness of the RBF method by adjusting the parameters of the RBF. For the numerical simulations, I set the physical area to be $\Omega = [-0.025, 0.025]\times[-0.025, 0.025]$. I chose the grid size to be $h = 8*10^{-4}$. The plane wave source we applied is $H = 0.1*\sin(\omega*t)$, where $\omega = 2\pi f$ and the operating frequency of the wave source $f = 40.0\; GHz$. set the PML to have a thickness that 12 times the grid size. In this test, the total time step number is $TN = 200$, and the splitting step of time $\tau = 10*10^{-13}$, that is, final simulation time $T = 0.2ns$. Some parameters are given below.

The functions $\sigma_{x}$ and $\sigma_{y}$ corresponding to the PML model are as follows.

I chose $c = -log(P_{0})(l+1)\frac{c_{v}}{2d}$, where $d = 12h$ is the thickness of the PML. $P_{0} = 1.0\times10^{-7}$ is the desired reflection error, and the polynomial degree was set as $l = 4$. Here, $c_{v} = 3\cdot10^{8}$ m/s is the speed at which the waves propagate in the vacuum.

4.1.   Example 1

In this example, I simulate an electromagnetic concentrator. The sketch of this device and sketch of the distribution of the collocation nodes is shown in Figure 2. The position of the wave source $x = -0.025$ m, $y$ range of values $[-0.025, 0.025]$ m, as shown by the red line in Figure 2. In order to clearly see the collocation points, I applied the the coarse grid $h = 16*10^{-4}$.

In order to see more clearly how the waves propagate in the concentrator, some snapshots of the electric field $E_{y}$ simulated by the MQ RBF with $\gamma = 1$, $\delta = 6*10^{-4}$ in Figure 3 and $\delta = 5*10^{-3}$ in Figure 4. The CPU times for the data in Figures 3 and 4 are 152.75 s and 149.37 s respectively for this 200 time step simulation.

The numerical simulation results show that the plane source wave propagates normally in the vacuum before entering the metamaterial region. After the wave passes through the metamaterial region, it slowly concentrates in the region $r'\leq R_{1}$. Then out of the metamaterial area, the wave propagates to the right as usual in the vacuum. Figures 3 and 4 clearly shows how the wave deforms in the metamaterial region. The results we obtained with the RBF meshless method produced a concentration phenomenon similar to that observed via the finite element method ^[36]. From Figures 3 and 4, we can see that the concentrated structure has almost no scattering at the interface of the medium, and has an electromagnetic concentration effect.

Here, I want to point out that the parameters $\gamma$ and $\delta$ in the MQ RBF are very sensitive to the simulation of the numerical results. Through numerical tests, I found that we must choose $\gamma = 1$, otherwise no matter what $\delta$ value is used, we will not get the desired reasonable results. Therefore, I assumed $\gamma = 1$ and have found that good numerical simulation results can be obtained for any $6*10^{-4}\leq\delta\leq5*10^{-3}$. The numerical results show that the image associated with $\delta = 6*10^{-4}$ is not as clear as the image produced when $\delta = 5*10^{-3}$. Within the range of parameter $\delta$ selection, the image becomes more and more clear as the value of the parameter increases.

4.2.   Example 2

In this example, I use another basis function, the Gaussian RBF, to simulate the wave propagation phenomenon of the electromagnetic concentrator. All other parameters are the same as in the previous example. Some snapshots of the electric field $E_{y}$ simulated by the Gaussian RBF with $s = 1.8*10^{6}$ in Figure 5 and $s = 4*10^{5}$ in Figure 6. The CPU times for the data in Figures 5 and 6 are 161.16 s and 159.50 s respectively for this 200 time step simulation.

The analysis of the numerical simulation results are the same as in Example 1. Here I point out the influence of the parameters of the Gaussian RBF on the numerical results. Through numerical tests, I found that the selection range for the only parameter $s$ in the Gaussian RBF is $[4*10^{5}, 1.8*10^{6}]$, other parameters can not simulate good results. The numerical results show that the image associated with $s = 1.8*10^{6}$ is not as clear as the image produced when $s = 4*10^{5}$. Within the range of parameter $s$ selection, the smaller the parameter, the clearer the image. Through the testing of two RBF, we found that the method is able to reproduce the same physical phenomena with the sophisticated finite element method of our previous work^[36].

5.   Conclusions

In this paper, I developed a leap-frog RBF meshless method for solving wave propagation problems in electromagnetic wave concentrator devices. The numerical results of MQ RBF and Gaussian RBF show that this RBF method is very effective when the correct free parameters are selected. How to systematically select the correct free parameters in RBF is still an open issue. In the future, I will continue to study the theoretical analysis of the RBF meshless method and how to select free parameters. In addition, I will also consider the application of the meshless method to simulate electromagnetic invisibility cloaks.

Acknowledgments

The author is very grateful to the anonymous referees for their helpful comments that improved the paper. Bin He's research was supported by the Gansu Provincial Educational Science and Technology Innovation Foundation, China (No: 2022A-040).

Conflict of interest

The author declares that he has no conflicts of interest.