Open access peer-reviewed chapter - ONLINE FIRST

Signal Propagation in Soil Medium: A Two Dimensional Finite Element Procedure

By Frank Kataka Banaseka, Kofi Sarpong Adu-Manu, Godfred Yaw Koi-Akrofi and Selasie Aformaley Brown

Submitted: March 26th 2021Reviewed: July 8th 2021Published: August 4th 2021

DOI: 10.5772/intechopen.99333

Downloaded: 29

Abstract

A two-Dimensional Finite Element Method of electromagnetic (EM) wave propagation through the soil is presented in this chapter. The chapter employs a boundary value problem (BVP) to solve the Helmholtz time-harmonic electromagnetic model. An infinitely large dielectric object of an arbitrary cross-section is considered for scattering from a dielectric medium and illuminated by an incident wave. Since the domain extends to infinity, an artificial boundary, a perfectly matched layer (PML) is used to truncate the computational domain. The incident field, the scattered field, and the total field in terms of the z-component are expressed for the transverse magnetic (TM) and transverse electric (TE) modes. The radar cross-section (RCS), as a function of several other parameters, such as operating frequency, polarization, illumination angle, observation angle, geometry, and material properties of the medium, is computed to describe how a scatterer reflects an electromagnetic wave in a given direction. Simulation results obtained from MATLAB for the scattered field, the total field, and the radar cross-section are presented for three soil types – sand, loam, and clay.

Keywords

  • electromagnetic signal in soil
  • finite element method
  • radar cross-section
  • underground wireless communication

1. Introduction

Signal propagation underground and in soil medium constitute the backbone of the Internet of underground things (IoUT), which power many applications such as precision agriculture, border monitoring for intrusion detection, pipeline monitoring, etc. [1, 2, 3]. In the recent past, the study of signal propagation in soil medium for underground wireless communication has focused mainly on empirical techniques [4, 5, 6, 7, 8]. The most commonly used modeling techniques for implementing IoUT include electromagnetic waves, magnetic induction, and acoustic waves [9, 10]. Electromagnetic field analysis is performed in this chapter using numerical modeling with the finite element method to examine the signal strength in the soil medium.

The process of numerical modeling of how electromagnetic fields propagate and interact with physical objects and the environment is usually referred to as computational electromagnetics (CEM), numerical electromagnetics. The primary motivation of this process is to develop efficient approximations to Maxwell’s equations through numerical schemes for cases where closed-form analytical solutions of Maxwell’s equation cannot be obtained due to the complexity of geometries, material parameters, and boundary conditions. Therefore, several real-life problems that are not analytically computable, such as electromagnetic scattering, antenna radiation, electromagnetic wave propagation, electromagnetic compatibility, etc., can effectively be solved by numerical techniques. The mathematical model of the electromagnetic problem is usually obtained in terms of partial differential equations, integral equations, or integro-differential equations derived from Maxwell’s equations and a set of a priori constraints of the problem such as boundary and initial conditions material parameters and geometry. The problem is ideally defined on an infinite-dimensional function space. Numerical methods apply a discretization to the continuum to reduce infinite degrees of freedom to a finite degree of freedom. In other words, the solution of an infinite dimension-dimensional problem is projected into a finite-dimensional space. Hence, the solution to the problem becomes amenable on a digital computer. The main philosophy in most of the numerical methods is to apply the divide-and-conquer strategy. The idea is to divide an intractable continuous problem into smaller pieces (divide), express the solution over each small piece (conquer), and then combine the piecewise solutions to obtain a global solution. In this chapter, FEM will be applied to the two-dimensional boundary value problem in EM wave propagation through the soil to evaluate the signal strength of the wave propagation in soil. This evaluation will be based on the incidence angle of the transmitted wave. The radar cross-section of the scatterer will be used to evaluate the direction of the wave.

Advertisement

2. Time-harmonic EM model

When deriving the wave equation for electromagnetic wave propagation through a vacuum or a dielectric medium such as soil or pure water, the free charge density ρf=0and the free current density Jf=0. However, in the case of conductors like seawater or metals, we cannot control the flow of charges and, in general, Jfis certainly not equal to zero. With this, Maxwell’s equations for linear media assume the form [11].

i·E=1ερfiiixE=Btii·B=0ivxB=μσE+μεEtE1

Apply the curl to (iii) and (iv), and we obtain modified wave equations for the electric field Eand the magnetic field B:

2u=με2ut2+μσutE2

where urepresents the scalar component of the electric or magnetic field, εis the permittivity of the medium μ=μrμ0, is the magnetic permeability (μris the relative permeability of the soil μr=1, for non-magnetic soil, μ0=4πx107N/A2is the magnetic constant in vacuum), and σis the electric conductivity of the soil medium. Eq. (2) admits the plane wave solution

uxt=ueikxωtE3

where ω=2πfis the angular frequency, tis the time, and kis the complex wavenumber or propagation constant, which can be derived from Eq. (3). The boundary value problem (BVP) used to solve the time-harmonic electromagnetic problem in 2-D, can be expressed in its generic form as

xpxutypyuy+qu=fforxyΩE4
u=u0onΓD:DirichletBoundaryConditionE5
pxutâx+pyuyây.n̂=βonΓN:NeumannBoundaryConditionE6
pxutâx+pyuyây.n̂+αu=βonΓM:MixedBoundaryConditionE7

where uxyis the unknown function to be determined, and pxxypyxyqxyfxyare given functions. u0α, and βare given functions in boundary conditions (BCs); ΓDΓNand ΓMrefers to boundaries where Dirichlet, Neumann, and mixed BCs are imposed, respectively; n̂=âxnx+âynyis the unit vector normal to the boundary in the outward direction.

The weak form is used to provide the finite element solution. The weak form is written as [12, 13, 14].

Ωpxuxψx+pyuyψydsΓψpxuynx+pyuynydl+ΩqψudsΩψfds=0E8

where ψxyis the weight function.

The weak form is applied in each element domain, and element matrices are formed by expressing the unknown function as a weighted sum of nodal shape functions. The sum of line integrals of two neighboring elements cancels out while combining the element matrices. Therefore, the line integral can be omitted for interior elements and should be considered only for elements adjacent to the boundary. Due to its special form, the line integral makes easier the imposition of mixed types of BCs. The mesh generation will be discussed first, and the shape functions will be given, and then the finite element solution of the time-harmonic problem of electromagnetic scattering in a dielectric medium (soil) will be presented.

3. The Helmholtz time harmonic model

In 2-D, the geometry and boundary conditions do not vary along an axis (say the z-axis). Hence fields can be represented as a superposition of fields of two orthogonal polarizations using the linearity property. Any field can be decomposed into transverse magnetic (TM) and transverse electric (TE) parts for the z-variable. In the TM case (horizontal polarization), only the z-component of the electric field E=âzEzxyand the x-y-components magnetic field H=âxHxxy+âyHy(xy)exist. However, in the TE case (vertical polarization), which is a dual of TM, only the z-component of the magnetic field H=âzHzxyand the x-y-components electric field E=âxExxy+âyEy(xy)exist.

The material tensors in this case are defined as

εrc=εrcxxεrcxy0εrcyxεrcyy000εrczzwithεrcsub=εrcxxεrcxyεrcyxεrcyy,E9
μr=μrxxμrxy0μryxμcyy000μrzzwithμrsub=μrxxμrxyμryxμryy,E10

The generalized homogeneous Helmholtz equation in TM and TE case can be written, respectively, in the following forms [15].

ΛμEz+k02εrczzEz=0,E11
ΛεHz+k02μrzzHz=0,E12

Where

Λμ=1μrsubTμrsubT,E13
Λε=1εrcsubTεrcsubT,E14

Where the superscript T indicates the transpose and .refers to the determinant of the corresponding submatrix; k0is the free-space wavenumber; εrc=εr/ωε0is the complex relative permittivity εr, and μrare the relative permittivity and permeability, respectively; and σis the conductivity. For isotropic mediums, Eq. (8) becomes

μr1Ez+k02εrcEz=0,E15

Or

x1μrEzx+x1μrEzy+k02εrcEz=0,E16

The same way Eq. (11) becomes

εrc1Hz+k02μrHz=0,E17

Or

x1εrcHzx+x1εrcHzy+k02μrHz=0,E18

4. Scattering from a dielectric medium (soil)

In [16], a new wave number model is proposed with the combination of the Peplinski principle and multiple scattering from particles in the soil medium. The new wave number is used in the computation of the path loss. In another recent work, sensitivity analysis of the Ku-band scattering coefficient to soil moisture was performed under single-polarized, dual-polarized, and dual-angular combinations [17]. Similarly, a model of parabolic equations for reflection and refraction in an environment with an obstacle where the area is decomposed into two different domains. The discrete mixed Fourier transformation is used to compute the field strength in the upper subdomain, and the finite difference method is used to calculate the field strength in the lower subdomain [18].

In our case, an infinitely large dielectric object of an arbitrary cross-section is considered and illuminated by an incident wave that is not a function of z. An illustration of the scattering problem for TM mode is shown in Figure 1 where a general FEM is depicted. The problem is defined in TE mode by replacing the electric field with a magnetic field.

Figure 1.

Electromagnetic scattering in soil: FEM modeling with PML.

4.1 The perfectly matched layer (PML)

Since the domain extends to infinity, an artificial boundary or layer is used to truncate the computational domain. That is an absorbing boundary condition (ABC) or perfectly matched layer (PML) [15].

Locally-Conformal PML (LC-PML) is a powerful PML method whose implementation is straightforward. It is implemented by replacing the real coordinates ρ=xyinside the PML region with their corresponding complex coordinates ρ˜=x˜y˜[15]. The computational domain is the union of the PML region, the free-space region, and the scatterer region.

Mesh generation is the process of representing the domain of interest as a collection of elements. The two most commonly used elements in 2-D problems are triangular and quadrilateral. Figure 2 shows the mesh that is formed by triangular elements for a rectangular domain of the boundary value problem (BVP). The figure was generated in MATLAB.

Figure 2.

Mesh formed by triangular elements for the rectangular domain of the BVP.

For a computational domain with such curved boundaries, discretization errors will occur due to the inability to capture the exact geometry of the domain. Triangular elements are preferred due to their simplicity and the possibility of developing algorithms for automatic triangulation for a computational domain such as Delaunay triangulation [19, 20]. Although fewer elements are needed when quadrilateral elements are used, Triangular elements are well-suited for complex geometries and cause fewer numerical dispersion errors. Furthermore, the calculation of element matrices is easier in triangular elements. Discretization error might inevitably occur unless sufficient mesh density is used. This happens regardless of whether triangular or quadrilateral elements are used for meshing. One way to overcome this problem is to refine the mesh. Mesh refinement might be needed especially if the geometry has a curved boundary or some corners, sharp edges, small features, or discontinuities. This might be important especially if linear interpolation functions are used. Another approach that increases accuracy is to use high-order elements at the expense of increased computational load. To achieve this, we use extra nodes within an element and use high-order interpolation functions. Figure 3 shows a simple mesh with six linear triangular elements having eight nodes. Local node numbers must follow the anticlockwise orientation in all elements to guarantee that the area of each element is obtained as a positive quantity. During the mesh generation, certain data arrays must be created [19].

Figure 3.

Mesh of a 2-D domain using linear triangular elements.

An element connectivity matrix of size Mx3, where M is the number of elements vectors of node coordinates (xand y), each of size Nx1, where N is the number of nodes representation of the given functions (px, py, q, and f) in each element (each of which is of size Mx1 Arrays containing special nodes and/or element. The connectivity matrix belonging to Figure 3 is given in Table 1.

Element (e)Node 1Node 2Node 3
1128
2238
3378
4347
5467
6456

Table 1.

Structure of element connectivity for triangular mesh.

The Delaunay function is used to create the connectivity matrix and automatically enumerates the nodes. The triangular mesh is generated for a rectangular domain. The element size (els) = Δh=λ0/5, Δh=λ0/10, Δh=λ0/20and the like. The higher the denominator, the finer the element. Figure 4 shows the mesh generated in MATLAB for els = Δh=λ0/15.

Figure 4.

Triangular mesh generated withΔh=λ0/15.

4.2 The scattered field formulation

Fields in the presence of the scatterer can be decomposed into two parts:

  • The first one is the incident field which is produced without the scatterer and

  • The second one is the scattered field produced by an equal amount of current induced on the scatterer or the surface enclosing the scatterer.

We first assume the TM polarization case, where the incident field Einc=âzEzinc, and the scattered field Escat=âzEzscat, are expressed in terms of the z-component, the z-component of the total field can be expressed as Ez=Ezscat+Ezinc. For an isotropic case, the total field satisfies the differential Eq. [9].

x1μrEzscatxx1μrEzscatyk02εrcEzscat=fxy,E19

Where the source term is given by

fxy=x1μrEzincx+x1μrEzincy+k02εrcEzinc,E20

This differential equation is a special form of Eq. (4a), where u=Ezscatxy, px=py=p=1μand q=k02εrc.

The incident field can be arbitrary and usually chosen as a uniform plane wave since the incident field sources are sufficiently far away from the object. The incident field with the unit magnitude is given by

Einc=âzEzinc=âzexpjkxcosφinc+ysinφinc,E21

where φincis the angle of incidence for the x-axis in cartesian coordinates.

Similar calculations can be performed in TE polarization mode by replacing the electric field with the magnetic field and permittivity with permeability. Hence, the differential equation in terms of the scattered magnetic field is given by

x1εrcHzscatxx1εrcHzscatyk02μrHzscat=fxy,E22

Where the source term is given by

fxy=x1εrcHzincx+x1εrcHzincy+k02μrHzinc,E23

This differential equation is a special form of Eq. (4a), where u=Hzscatxy, px=py=p=1εrcand q=k02μ.

For dielectric objects, the right-hand side of the matrix equation can be obtained by using the source term f(x, y) in each element. An alternative simpler approach is that the source term in Eq. (20) is just the differential operator applied to the known incident field within the object. The left-hand side of Eq. (19) is the same differential operator being applied to the unknown scattered field, which yields the left-hand side of the matrix equation Au, where urefers to the vector of nodal values of the scattered field. Therefore, Au=Auincmust be satisfied within the scatterer region. After forming the global matrix, as usual, the right-hand side vector can be modified by just multiplying the incident field vector by the global matrix, i.e. b=Auinc, only for entries bcorresponding to the nodes lying inside the object.

4.3 Radar cross section

The radar cross-section (RCS) of the scatterer is perhaps the most critical parameter that must be evaluated in the post-processing phase of FEM for the electromagnetic scattering problem [9]. RCS is the reflection of the scattering electromagnetic ware at an incident angle in a particular direction. In other words, it is the area capturing that amount of scattered power produced at the receiver in an isotropic medium. This is a density that is equal to that scattered by the actual target. It is a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and material properties of the object. In 2-D, it is mathematically defined as

σ2D=limr2πρufarscat2uinc2E24

Where ufarscat=1ρfφis the scattered electric or magnetic field at far-zone observed in a given direction, when ρsatisfies the inequality ρ2D2/λ, where Dis the largest dimension of the scatterer and λis the wavelength.

If the incident and observation directions are the same, the RCS is called monostatic or backscatter RCS; otherwise, it is referred to as bistatic RCS. In 2-D, RCS is usually normalized for λ(wavelength) or m (meter). The unit for σ2D/λis dBw, and for σ2D/mis dBm.

5. Scattering inside the dielectric medium

For each such element, the far-zone field is computed and then superposed over all elements, as follows:

Efarscat=ηk8πρej3π4ejkρi=1KJziejkxmicosφ+ymisinφΔli,E25

where Kis the number of elements adjacent to the boundary, xmiymithe i-th segment midpoint, and Δlithe i-th segment length.

Since FEM is formulated in terms of the scattered electric field, computation of the derivatives of the scattered field requires additional effort. This can simply be achieved by using the weighted sum of derivatives of shape functions in terms of nodal scattered fields.

5.1 Scattering inside the dielectric medium, TM case

For the dielectric object, the TM case, the interior part of the object should be included. Since the magnetic current density becomes nonzero, the integral containing the magnetic current density should be evaluated. The magnetic current density has x- and y-components, which is defined as M=âxMx+âyMy. Hence, it shows that ρ̂×M=cosφMysinφMxâz. The components of the magnetic current density can be determined in terms of the scattered electric field as follows:

M=E×n̂=nyEzâx+nxEzây,E26

Finally, the scattered electric field can be obtained as follows:

Efarscat=k8πρej3π4ejkρi=1Knxcosφ+nysinφEziejkxmicosφ+ymisinφΔli+ηk8πρej3π4ejkρi=1KJziejkxmicosφ+ymisinφΔliE27

Here, Ezican be determined as the average of nodal field values connected to the boundary.

5.2 Scattering inside the dielectric medium, TE case

For dielectric object, TE case: The magnetic current density is nonzero and has only z-component, M=âzMzbecause ρ̂×ρ̂×M=Mzâz. The scattered magnetic field is determined as:

Hfarscat=k8πρej3π4ejkρi=1KJxisinφJyicosφejkxmicosφ+ymisinφΔli+ηk8πρej3π4ejkρi=1KMziejkxmicosφ+ymisinφΔliE28
Advertisement

6. Simulation results and discussion

The operating frequency used in the simulation is 300 MHz, and the free-space wavelength is λ0=1m. The mesh size is approximatively Δh=λ0/40. The dielectric medium is a rectangular computational domain, the relative permittivity for sandy soil εrc1=4+i0.1, the relative permittivity for loamy soil εrc2=12+i1.8, and the relative permittivity for clay soil εrc3=25.3+i5.7. The angle of incidence is 0o, that is, the incident plane wave travels in the ydirection. The wavelength inside the dielectric medium (soil) decreases in terms of the wavelength in free-space and relative permittivity, i.e. λ=λ0/εrc, and hence, mesh size Δh=λ0/40resolution corresponds to λεrc/40resolution inside the dielectric. Therefore, smaller elements are used in dielectric media to preserve the level of accuracy. All simulation results are obtained using MATLAB.

Figure 5 shows scattering from a dielectric medium for TE mode, with sandy soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe a minimal electric field scattering and an intense total electric field inside the dielectric medium of sandy soil. This is because this medium is porous and allows for better signal propagation.

Figure 5.

Scattering from a dielectric medium (Sandy soil): (a) scattered field (TE) (b) Total field (TE).

Figure 6 shows scattering from a dielectric medium for TM mode, with sandy soil characteristics. In this case, we observe a minimal magnetic field scattering in (a) and an intense total magnetic field in (b) inside the dielectric medium of sandy soil.

Figure 6.

Scattering from a dielectric medium (Sandy soil): (a) scattered field (TM) (b) Total field (TM).

Figure 7 shows scattering from a dielectric medium for TE mode, with loamy soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe an intense electric field scattering and a minimal total electric field inside the dielectric medium of loamy soil. This is because this medium is less porous and presents some challenges in signal propagation.

Figure 7.

Scattering from a dielectric medium (loamy soil): (a) scattered field (TE) (b) Total field (TE).

Similarly, Figure 8 shows scattering from a dielectric medium for TM mode, with loamy soil characteristics. In this case, we observe a high magnetic field scattering in (a) and a low total magnetic field in (b) inside the dielectric medium of loamy soil.

Figure 8.

Scattering from a dielectric medium (loamy soil): (a) scattered field (TM) (b) Total field (TM).

Figure 9 shows scattering from a dielectric medium for TE mode, with clay soil characteristics. In (a) and (b), the scattered and total fields are shown respectively. We observe very high electric field scattering and a low total electric field inside the dielectric medium of clay soil. This is because this medium is non-porous and presents a very poor signal propagation.

Figure 9.

Scattering from a dielectric medium (clay soil): (a) scattered field (TE) (b) Total field (TE).

Similarly, Figure 10 shows scattering from a dielectric medium for TM mode, with clay soil characteristics. In this case, we observe a very high magnetic field scattering in (a) and a low total magnetic field in (b) inside the dielectric medium of clay soil.

Figure 10.

Scattering from a dielectric medium (clay soil): (a) scattered field (TM) (b) Total field (TM).

Figures 1113 show the bistatic RCS profiles to describe how scatterers reflect the incident electromagnetic wave in a given direction. This is the area intercepting that amount of power which, when scattered in the soil medium, produces at the receiver a density which is equal to that scattered by the actual target. The RCS is a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and properties of the soil medium. It is shown in the TE modes (a) and TM modes (b) for sandy soil, loamy soil, and clay soil, respectively.

Figure 11.

Radar cross section in Sandy soil for (a) TE mode and (b) TM mode.

Figure 12.

Radar cross section in loam soil for (a) TE mode and (b) TM mode.

Figure 13.

Radar cross section in clay soil for (a) TE mode and (b) TM mode.

7. Conclusion

In a two-Dimensional Finite Element Analysis of EM wave Propagation through the soil, a boundary value problem (BVP) used to solve the time-harmonic electromagnetic problem in 2-D, has been expressed in its generic form. In TM and TE cases, the Helmholtz model has considered an infinitely large dielectric object of an arbitrary cross-section for scattering from a dielectric medium and illuminated by an incident wave. Since the domain extends to infinity, an artificial boundary, an absorbing boundary condition (ABC), or a perfectly matched layer (PML), has been used to truncate the computational domain. The incident field, the scattered field, and the total field in terms of the z-component are expressed for the TM and TE modes. The radar cross-section (RCS), a function of several parameters, such as operation frequency, polarization, illumination angle, observation angle, geometry, and material properties of the medium, has been computed to describe how a scatterer reflects an incident electromagnetic wave in a given direction. Simulation results for the scattered field, the total field, have been presented for soil types, and the radar cross-section for different element refinements have also been presented.

DOWNLOAD FOR FREE

chapter PDF

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Frank Kataka Banaseka, Kofi Sarpong Adu-Manu, Godfred Yaw Koi-Akrofi and Selasie Aformaley Brown (August 4th 2021). Signal Propagation in Soil Medium: A Two Dimensional Finite Element Procedure [Online First], IntechOpen, DOI: 10.5772/intechopen.99333. Available from:

chapter statistics

29total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us