# Propagation Modeling of Point Source Excited Magnetoinductive Waves Based on a New Plane Wave Expansion Approach

Authors: Feng Liu; Zhijun Zhang; Lianqing Zhu; Yajun Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The signal fading in wireless underground sensor networks (WUSNs), which is caused by lossy media such as soil and sand, can be reduced by applying technology of magnetoinductive (MI) propagation. This technology can effectively establish a communication at very low frequency (VLF). In contrast to the previous studies in the literature, which mostly focus on the propagation of plane waves, we propose a new approach based on the plane wave expansion (PWE) to model the near field MI waves. The proposed approach is based on excitation of a point source, which is a common case in a practical WUSN. The frequent usage of square lattice MI structure is investigated. To verify the mathematical derivation, the simulation of time domain based on the fourth-order Runge-Kutta (RK) method is carried out. Simulation results show that the new model can provide a precise prediction to the MI wave’s propagation, with the computation load being one-tenth of that of the time domain simulation. The characteristics of the propagation of the MI waves are presented and discussed. Finally, the reflection on the edge of the MI structure is reduced by analysing the terminal matching conditions and calculating a method for matching impedances.


---

## Body

## 1. Introduction

Wireless underground sensor networks (WUSNs) become an active research in monitoring underground environments such as mineral contents, oil leakage, and formation stress [1–4]. Since the propagation medium of underground environments is not air and contains soil, rocks, and water, wireless communication in WUSN faces severe challenges caused by the high path loss.The classical wireless signal propagation for electromagnetic (EM) waves is not effective in working in underground environments [5–7]. Since the waves with low frequency EM suffer from lower fading in loss media, this can be helpful to set the WUSN working at the VLF band similar to most of the telecommunication systems [8–10]. However, large size antennae and high drive voltage are required to transmit EM waves in VLF band. This approach is usually impracticable in underground environments. The emergence of MI propagation technology provides a practical solution to this problem.The MI technology has been extensively studied in the metamaterials field [11–13]. In [14], Shamonina et al. applied the concept of MI waves for the loop current waves in one, two, and three dimensions. Syms et al. studied the positive and negative refraction of MI waves in two dimensions [15]. In addition, they applied the concept of MI waveguide with the quantitative analysis to several kinds of waveguide devices [16]. The MI waveguide was introduced by Sun and Akyildiz in [17] to the WUSNs field as a new proposed technology in the physical layer of sensor networks. The studies in [18–21] show that the MI waveguides have the capabilities of transmitting powers. Therefore, the MI waveguides are useful to recharge the WUSN nodes. In conclusion, the current studies of the MI structures mostly focused on the plane wave propagation. However, in a real WUSN, signals are always transmitted by a single loop, and the wave excited by the point source cannot be considered as a plane wave.Since the frequency of MI networks only depends on the resonance frequency of the split ring resonators, this can be exploited to adjust the load capacitor and the number of windings to set the MI structure working at the VLF band, keeping the transmit antenna at a low size level. At the same time, the resonant capacitor will provide auxiliary power potential to drive the loop coils, which reduce the required driving voltage to a very low level.Since the propagation model and wave patterns provide fundamental basis for link budget, topology controlling, and the design of communication protocol, the aim of this paper is to show the characteristics of the near field propagation of the MI waves that are excited by a point source. The frequent usage of square lattice MI structure is investigated. The simulation of time domain based on the fourth-order Runge-Kutta (RK) method is carried out to verify the mathematical derivation.The rest of this paper is organized as follows. Section2 presents the system model, including the dispersion equations of MI structure. The wave propagation analysis based on the plane wave expansion method is illustrated in Section 3. The simulation results and the verification using a time domain simulation are given in Section 4. The discussions are presented in Section 5, and Section 6 concludes this paper.


## 2. System Model

In classical WUSNs, nodes are distributed in monitoring areas. A square lattice with the same lattice constant “a” can be considered a simple topological structure and suitable for WUSNs. As shown in Figure 1, every node in WUSNs is equipped with a capacitive loaded loop antenna.Figure 1
Schematic representation of a capacitively loaded wire loop.In order to keep the isotropy, only the coplanar configuration of the network is discussed in this paper which is shown in Figure2.Figure 2
Geometry of the loops (m is the column number and n is the row number).For the array of MI loops, only the mutual effects between nearest neighbours need to be considered [22]. Thus, in a square lattice WUSN, each node can be considered to have relationships only with four neighbours.In the coplanar loops, as shown in Figure3, the flux is excited by a loop itself which has opposite direction compared with that excited by other loops. Thus, the mutual inductance is negative and denoted by -M. Figure 4 shows a simple equivalent circuit for one loop in the two-dimensional (2D) wave as shown in Figure 1. In Figure 4, each resonator is modelled by the connection of an inductance (L), a capacitance (C), and four mutual inductances of -M, which are coupled from its nearest neighbors. In addition, in each loop, Kirchhoff’s voltage law is applied considering the relations with four neighbour loops’ currents (Im,n).Figure 3
Mutual inductance between two coplanar loops. (ϕs is the self-induced flux and ϕm is the mutual induced flux.)Figure 4
Equivalent circuit model of the loop in the lattice array. (M is the mutual inductance, L is the self-inductance, and C is the capacitor load.)In the lossless case, we obtain(1)um,n(t)=Ldim,n(t)dt+1C∫0tim,n(t)dt-M(dim-1,n(t)dt+dim+1,n(t)dtlllllllllllllll+dim,n-1(t)dt+dim,n+1(t)dt),
where um,n(t) is the voltage source of the node at nth row and mth column, im,n(t) is the antenna current, L and C are the inductance and capacitance of the antenna, respectively, and -M is the mutual inductance.In the frequency domain, (1) can be rewritten as
(2)U=ZIm,n-ZM(Im-1,n+Im+1,n+Im,n-1+Im,n+1),
where Z=jωL+1/jωC is the loop impedance and the coupling term is defined as -ZM=-jωM. When the voltage source is vanished, we assume a travelling-wave solution for the simplest solution as
(3)Im,n=Iejω0te-jk→r→,
where k→=ux→+vy→ is the wave vector and r→=xx→+yy→ is the position vector. Substituting (3) into (2), we obtain the constraint equation
(4)ZIe-j(mu+nv)a-ZmIe-j(mu+nv)a×(ejau+e-jau+ejav+e-jav)=0.Equation (4) can be rewritten as
(5)cos⁡au+cos⁡av=(ω2-ω02)/ω2κ,
where ω0=1/LC is the resonance frequency and κ=2M/L is the coupling coefficient. Since ((ω2-ω02)/ω2)/κ<2 and κ is a very small value in WUSNs, the bandwidth will be very narrow. Figure 5 shows some curves of wave vectors. Considering the monochrome signal at the resonance frequency ω0, (5) can be rewritten as
(6)cos⁡au+cos⁡av=0.Figure 5
The propagation vector curves.The resonant frequency wave vector curves are four straight lines that are joined end to end forming a square frame. In Figure5, the curve of the k→ at lower frequencies is located inside the square frame of the resonance frequency curve, with higher frequency curves outside of the enclosure.


## 3. Propagation of MI Waves

### 3.1. The Plane Wave Expansion of MI Waves Excited by Point Source

We can consider theIe-jk→r→ in (3) as a space signal whose Fourier transformation is known as wavenumber spectrum. As it is known, certain wave vectors are allowed in a 2D MI waveguide. In addition, we can describe the phenomena in other words that the MI waveguide is a space system and produces different response to the MI wave amplitude with different propagation constants; that is, MI waveguide is a system of space signal. It can be verified that the MI waveguide space system is a linear space invariant system (LSI) and is mathematically similar to a linear time invariant system (LTI).Since the space parameters of the loop are the coordinates of the circle centre, a single loop exciting source can be considered as a point source and expressed with a pulse function. The wavenumber spectrum of point source is 1. As shown in Section2, the MI waveguide is an ideal filter of wave vector with unit pass-band gain and the stop-band gain of 0, and there is a uniform system gain for every possible k. As proven in the current studies, the possible form of wave vectors is four separate continuous curves in the loop array of square lattice. Thus, the space distribution of the loop current excited by a single loop source (point source) represents a single frequency signal and can be written using the inverse Fourier transformation as
(7)I(x,y,ω)=A∫L(ω)ej(ux+vy)dl,
where L(ω) is the curve that is composed of the dots with coordinates (u,v) as shown in (5). A in (7) can be obtained considering the node at the origin of the space coordinates, which can be expressed by
(8)A=I(0,0,ω)∫L(ω)dl,
where I(0,0,ω) is the amplitude of the loop current at position of x=0 and y=0. In addition, I(0,0,ω) could be a complex number.Considering the current source is a periodic signal with a period ofT0, the loop current can be expressed by
(9)I0=∑q=-∞+∞(Fqejq(2π/T0)t).The current amplitude at any coordinates can be obtained from (6), (7), and (8) as
(10)I(x,y,t)=∑q=-∞+∞{[Fq∫L(q(2π/T0))dl∫L(q(2π/T0))ej(ux+xy)dl]lllllllllllllllllllllllllllll×ejq(2π/T0)t[Fq∫L(q(2π/T0))dl∫L(q(2π/T0))ej(ux+xy)dl]}.WhenT0→∞, we can get the current at position (x,y) excited by any nonperiodic point source:
(11)I(x,y,t)=∫-∞+∞{Ix=0,y=0(ω)∫L(ω)ej(ux+vy)dl∫L(ω)dlejωt}dω.Because the propagation of vector curve,L(ω), is symmetrical in both the X and Y coordinates, (9) and (10) are rewritten as
(12)I(x,y,t)=∫-∞+∞{4Ix=0,y=0(ω)∫u1u21+v′2dulllllllllllllllllll×∫u1u21+v′2cos⁡(ux)cos⁡(vy)duejωt4Ix=0,y=0(ω)∫u1u21+v′2du}dω,
where v=arccos((ω2-ω02)/ω2κ-cos⁡(au))/a as shown in (5).


### 3.2. The Reflection and Matching of Regular Boundaries

In most applications, two kinds of regular edges are commonly used as shown in Figure6.Figure 6
The line edge of a 2D MI waveguide.
(a)
Side edge
(b)
HypotenusesIn Figure6, the loops in the shadow are the actual waveguide. For the convenience of analysis, the virtual loops are drawn with dotted lines and attached to the edge with zeros currents as the boundary conditions.In an infinite 2D MI waveguide, an identical boundary condition is a mirror current source at the position symmetrical to the actual source for the straight line through the centers of attached loops. The MI waves with the reflection of the edge in the actual areas of the 2D MI waveguide could be considered as the superposition of the waves. In addition, these waves are excited by the two point sources.As shown in Figure7, in a finite 2D MI waveguide with more than one edges, the theory of geometrical optics can help us to find all of the mirror sources. In addition, this is useful to obtain the wave propagation that is resulted by the superposition of all PWE from the sources.Figure 7
The multimirror source.In some applications, the reflection is very harmful to the signal propagation. Thus, the edges should be designed to absorb any incident powers and the network should be terminally matched.According to the situation in Figure6(a), we first assume that there exists an infinite 2D MI waveguide. Furthermore, we remove half of the loops to form a line edge. As mentioned previously, the second loop (2 as shown in Figure 6(a)) has effect only on the first loop (1 as shown in Figure 6(a)). To keep everything normal in the first loop and remove the second loop, we should connect the impedance in series form as the second loop is still there. Thus, (2) can be rewritten in (13) with substituting ZM by the matching impedance Zmatching,
(13)U=ZIm,n-ZM(Im-1,n+Im,n-1+Im,n+1)-ZMIm+1,n=ZIm,n-ZM(Im-1,n+Im,n-1+Im,n+1)+ZmatchingIm,n,
where Im,n is not zero and Zmatching is obtained as
(14)Zmatching=-ZMIm+1,nIm,n.Substituting (12) into (14) and considering the Fourier transform formula, we obtain
(15)Zmatching=-jωM∫u1u21+v′2cos⁡[u(m+1)a]cos⁡(vna)du∫u1u21+v′2cos⁡(uma)cos⁡(vna)du=-jωX(ω,m,n),
where X(ω,m,n) is a real function of angular frequency and loop position. Thus, the matching device would be an inductor or a capacitor and its value should be calculated for every position and frequency.Note that it is impossible to realize perfect terminal impedance matching in a 2D MI waveguide. However, an approximate matching may be acceptable for a narrow band telecommunication with the consideration of the matching impedance of the resonance frequency.


## 3.1. The Plane Wave Expansion of MI Waves Excited by Point Source

We can consider theIe-jk→r→ in (3) as a space signal whose Fourier transformation is known as wavenumber spectrum. As it is known, certain wave vectors are allowed in a 2D MI waveguide. In addition, we can describe the phenomena in other words that the MI waveguide is a space system and produces different response to the MI wave amplitude with different propagation constants; that is, MI waveguide is a system of space signal. It can be verified that the MI waveguide space system is a linear space invariant system (LSI) and is mathematically similar to a linear time invariant system (LTI).Since the space parameters of the loop are the coordinates of the circle centre, a single loop exciting source can be considered as a point source and expressed with a pulse function. The wavenumber spectrum of point source is 1. As shown in Section2, the MI waveguide is an ideal filter of wave vector with unit pass-band gain and the stop-band gain of 0, and there is a uniform system gain for every possible k. As proven in the current studies, the possible form of wave vectors is four separate continuous curves in the loop array of square lattice. Thus, the space distribution of the loop current excited by a single loop source (point source) represents a single frequency signal and can be written using the inverse Fourier transformation as
(7)I(x,y,ω)=A∫L(ω)ej(ux+vy)dl,
where L(ω) is the curve that is composed of the dots with coordinates (u,v) as shown in (5). A in (7) can be obtained considering the node at the origin of the space coordinates, which can be expressed by
(8)A=I(0,0,ω)∫L(ω)dl,
where I(0,0,ω) is the amplitude of the loop current at position of x=0 and y=0. In addition, I(0,0,ω) could be a complex number.Considering the current source is a periodic signal with a period ofT0, the loop current can be expressed by
(9)I0=∑q=-∞+∞(Fqejq(2π/T0)t).The current amplitude at any coordinates can be obtained from (6), (7), and (8) as
(10)I(x,y,t)=∑q=-∞+∞{[Fq∫L(q(2π/T0))dl∫L(q(2π/T0))ej(ux+xy)dl]lllllllllllllllllllllllllllll×ejq(2π/T0)t[Fq∫L(q(2π/T0))dl∫L(q(2π/T0))ej(ux+xy)dl]}.WhenT0→∞, we can get the current at position (x,y) excited by any nonperiodic point source:
(11)I(x,y,t)=∫-∞+∞{Ix=0,y=0(ω)∫L(ω)ej(ux+vy)dl∫L(ω)dlejωt}dω.Because the propagation of vector curve,L(ω), is symmetrical in both the X and Y coordinates, (9) and (10) are rewritten as
(12)I(x,y,t)=∫-∞+∞{4Ix=0,y=0(ω)∫u1u21+v′2dulllllllllllllllllll×∫u1u21+v′2cos⁡(ux)cos⁡(vy)duejωt4Ix=0,y=0(ω)∫u1u21+v′2du}dω,
where v=arccos((ω2-ω02)/ω2κ-cos⁡(au))/a as shown in (5).


## 3.2. The Reflection and Matching of Regular Boundaries

In most applications, two kinds of regular edges are commonly used as shown in Figure6.Figure 6
The line edge of a 2D MI waveguide.
(a)
Side edge
(b)
HypotenusesIn Figure6, the loops in the shadow are the actual waveguide. For the convenience of analysis, the virtual loops are drawn with dotted lines and attached to the edge with zeros currents as the boundary conditions.In an infinite 2D MI waveguide, an identical boundary condition is a mirror current source at the position symmetrical to the actual source for the straight line through the centers of attached loops. The MI waves with the reflection of the edge in the actual areas of the 2D MI waveguide could be considered as the superposition of the waves. In addition, these waves are excited by the two point sources.As shown in Figure7, in a finite 2D MI waveguide with more than one edges, the theory of geometrical optics can help us to find all of the mirror sources. In addition, this is useful to obtain the wave propagation that is resulted by the superposition of all PWE from the sources.Figure 7
The multimirror source.In some applications, the reflection is very harmful to the signal propagation. Thus, the edges should be designed to absorb any incident powers and the network should be terminally matched.According to the situation in Figure6(a), we first assume that there exists an infinite 2D MI waveguide. Furthermore, we remove half of the loops to form a line edge. As mentioned previously, the second loop (2 as shown in Figure 6(a)) has effect only on the first loop (1 as shown in Figure 6(a)). To keep everything normal in the first loop and remove the second loop, we should connect the impedance in series form as the second loop is still there. Thus, (2) can be rewritten in (13) with substituting ZM by the matching impedance Zmatching,
(13)U=ZIm,n-ZM(Im-1,n+Im,n-1+Im,n+1)-ZMIm+1,n=ZIm,n-ZM(Im-1,n+Im,n-1+Im,n+1)+ZmatchingIm,n,
where Im,n is not zero and Zmatching is obtained as
(14)Zmatching=-ZMIm+1,nIm,n.Substituting (12) into (14) and considering the Fourier transform formula, we obtain
(15)Zmatching=-jωM∫u1u21+v′2cos⁡[u(m+1)a]cos⁡(vna)du∫u1u21+v′2cos⁡(uma)cos⁡(vna)du=-jωX(ω,m,n),
where X(ω,m,n) is a real function of angular frequency and loop position. Thus, the matching device would be an inductor or a capacitor and its value should be calculated for every position and frequency.Note that it is impossible to realize perfect terminal impedance matching in a 2D MI waveguide. However, an approximate matching may be acceptable for a narrow band telecommunication with the consideration of the matching impedance of the resonance frequency.


## 4. Numerical Result and Time Domain Simulation

Since it is difficult to obtain the analytical solution of the PWE in close forms, the numerical evaluation is carried out. Table1 shows the conditions of numerical analysis.Table 1
Numerical analysis conditions.
Item
Value
1
Shape
Circle
2
Loop radius
0.1 m
3
Turns
20
4
Loop space
1 m
5
Operating freq.
10 KHzWe consider node array with 121 by 121. The source is located in the array centre. The time domain simulation is carried out to verify the numerical results. As shown in (1) and using a loop at the mth column and nth row, the differential equations can be set up with the form as
(16)ddtum,n(t)=Ld2dt2im,n(t)+1Cim,n(t)-M[d2dt2im-1,n(t)+d2dt2im+1,n(t)llllllllllllll+d2dt2im,n-1(t)+d2dt2im,n+1(t)].For the loops at the edge and corner, some of the items are removed. Matching impedance is also inserted according to (15).Assuming thatdim,n(t)/dt=wm,n(t), (16) is rewritten as
(17)dim,n(t)dt=wm,n(t),ddtum,n(t)=Lddtwm,n(t)+1Cim,n(t)-M[dwm-1,n(t)dt+dwm+1,n(t)dtllllllllllllll+dwm,n-1(t)dt+dwm,n+1(t)dt].The second-order differential equation can be transformed into two equations with one order. Thus, the whole number of differential equations in this test is 29282. With the fourth-order RK method, we can obtain the time domain numerical solution. To limit the bandwidth of the signal, a rising cosine window is used to modulate the carrier as(18)u(t)=0.01*[1-cos⁡(ωst)]*sin(ω0t),0<t<2πωs,
where ωs=10π and ω0=20000π.Both PWE and simulation results reveal that the narrow MI wave packet propagates in four main directions. Figure8 shows a snapshot of the propagation at the time of 0.9 seconds.Figure 8
Snapshot of the propagation of MI narrow wave packet excited by a point source at the 61th row and the 61th column (the grey level represents the amplitude of the current).
(a)
The PWE result
(b)
The time domain simulationWe redraw the current amplitude of the loops along the diagonal line in Figure9. As can be seen in Figures 8 and 9, the PWE results match with the results of time domain simulations, except that there are more fluctuations in the PWE results. In addition, we see that due to the symmetry of the loop array, the propagation of the MI wave is also symmetrical at 61th row and 61th column, which divides the loop array into four identical parts. Figure 8 illustrates that the propagation speeds and the wave patterns of the wave packets (obtained with two different methods) are almost the same along the diagonal line. The wave packet patterns also show that there is serious dispersion with the signal transmission.Figure 9
The current amplitude along the diagonal line.Figure8 shows that the MI wave packet propagates along several different paths. This is clearly seen in Figure 10, where the amplitudes of the wave front at 0.9 s are given for different azimuth angles form 0 to 90 degrees. In Figure 10, the PWE results contain more chaotic ups and downs than the time domain simulation results. However, the outlines of the pattern are coincident.Figure 10
The wave front amplitude.
(a)
The PWE result
(b)
The time domain simulation


## 5. Discussion

The PWE method is based on the plane wave superposition. Thus, at any time, the space is already full of steady plane waves. As can be seen in Figures8(a) and 8(b), before arriving of the narrow wave packets, there are plenty of wave distributions in the spare area. Moreover, in the real and negative times, the future wave packets introduce fluctuations at presence of time and area (without physical meaning). This explains why the PWE results contain more ripples than the simulations of time domain.The simulation of time domain is based on solving the differential equations and the Runge-Kutta method. These are commonly used as tools for this purpose. The loop array of 121 by 121 includes solving 29282 differential equations and introducing a huge amount of calculations. This occupies a great deal of computing resources. Therefore, in this paper, the running time takes about 70 hours to accomplish the MI wave simulation using a computer with a 3.4 GHz 4-Core CPU and 32 GB memory. However, it takes only lower running time than half an hour to accomplish identical simulation with the PWE method. Therefore, even though the PWE method is only an approximation of the real propagation, it is a very efficient method to analyse the MI waves.To explain the formation of the pattern of the MI wave propagation, we review the periodic structure about the bandwidth and the energy flow. Using (5), we can obtain the bandwidth of different propagation directions.The normalized bandwidth of the simulated loop array is given in Figure11.Figure 11
The bandwidth of different propagation directions.As can be seen in Figure11, the diagonal direction has the maximum bandwidth.By considering the center node(m,n) and its four neighbors in Figure 2, the neighbors’ power which is transmitted and received by the center node is expressed as
(19)Plefttocenter=ZMIm-1,nIm,n*,Pcentertoright=ZMIm,nIm+1,n*,Pcentertoup=ZMIm,nIm,n+1*,Pdowntocenter=ZMIm,n-1Im,n*.Substituting (3) into (19) and assuming that I=Im,n, we can obtain the dispersion equation as
(20)Pcentertoright=-ωMI2sin(ua),Pcentertoup=-ωMI2sin(va),
where Plefttocenter is the power flow in x direction and is marked with Px and the Pdowntocenter is marked with Py; thus (20) can be rewritten as
(21)Px=-ωMI2sin(ua),Py=-ωMI2sin(va).The result in (21) is also deduced through the group speed calculation. Figure 12 shows the normalized energy propagating speed of the simulated loop array.Figure 12
The energy flow speed.As illustrated in Figures11 and 12, the main part of the signal power is bounded into the diagonal directions and this is called the main propagation direction (MPD). In other directions, the propagation is more complicated. To illustrate this issue, an example in Figure 13 is presented.Figure 13
The established process of APD.In Figure13, a narrow wave packet, called W1, departs from the source into the direction of D1 and the energy flowing direction will be D3 towards the MPD. The secondary wave packet, called W2, which is excited by W1, will eventually turn to the direction of D2. This is due to the energy flow and the new propagation path which is parallel to the MPD. This path is called auxiliary propagating direction (APD). As shown in Figure 8, the wave packet propagates along the MPD directions; thus new secondary waves will be excited and resulted in several trajectories of discrete parallel wave propagation.


## 6. Conclusion

The numerical results and time domain simulations are presented. The results and simulations show that the proposed PWE method is efficient and accurate to model the MI wave propagation in the near field of a point source. It can speed up the analysis dozens of times. From the theoretical derivations and simulations, we can conclude that the MI wave propagation in a 2D resonance array is not similar to that of an EM wave in some continuous mediums. The MI waves, which are excited by a point source, mainly propagate along the diagonal directions (MPDs). Secondary waves form multipropagation trajectories (APDs) are parallel to the MPDs. Since the main power is transmitted along the MPDs, the 2D array should employ a switching and relay mechanism when setup of the communication between two nodes is not in the same MPD. For future works, the effects of other topologies, the communication protocols, and beam forming technologies will be analysed in the MI 2D waveguides.



---

*Source: 951718-2015-03-04.xml*