# Nondestructive Approach for Complex-Shaped Cracks in Concrete Structures by Electromagnetic Waves with FDTD Technique

Authors: Ummu S. Sener; Sebahattin Eker
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Concrete cracks have no specific shape and do not show linearity. Since the natural occurrences of concrete cracks make simulation identification difficult, rectangular step function and a dynamic geometry are used to define a concrete surface crack in the natural process. A novel interior crack expression is obtained by accepting the area between two curves as a crack filled by air in concrete and modeling this area like a Riemann integral domain. Taking the partition of this integral domain, the most realistic definition of the crack is made. Electromagnetic (EM) waves are utilized for numerical simulation after identifying the defects, cracks, rebars, and geometry of concrete. Three different simulation setups with complex geometries with two different surface cracks and one internal crack are simulated using a finite-difference time-domain (FDTD) method with Gaussian pulse wave excitation. Simulations are obtained using both transverse electric (TEz) waves and transverse magnetic (TMz) waves and the results are compared with each other. Air-dried concrete specimens are molded following simulation setups with surface cracks and measurements are made nondestructively with a Vivaldi antenna array in the frequency range of 0.4–4.0 GHz. The reflection and transmission coefficients are validated by comparing the data obtained using the measurement with the results obtained from numerical simulation.


---

## Body

## 1. Introduction

Early detection of damage and cracks in concrete structures and infrastructures has great importance for public safety economic saves. Reinforced concrete structures require routine, precise, and reliable monitoring; nondestructive testing (NDT) is carried out to periodically check the structural integrity of many concrete structures, such as bridge trays, plates, and tunnels [1–3]. In the inspection of concrete structures with nondestructive testing, the integrity of the structure is not compromised, and the serviceability of the structure is not prevented since this inspection can be done in a short time [4–8]. NDT techniques such as acoustic, thermal, optical microwave, and radar methods may be the available option in many infrastructures and buildings [9, 10]. Microwave and radar NDT method is a widely used method as no chemicals are used; it is treated with nonionizing waves and can be used in concrete, plastic, ceramic, and many composite materials [11–13]. Detecting gaps in materials that have been eroded after long-term use, detecting porosity in ceramics, measuring epoxy resin curing rates, determining the amount of moisture in dielectric materials, and examining dielectric materials such as composites are important application areas for microwave NDT [14–16]. In addition, evaluations of changes in composites and concrete structures resulting from environmental factors such as temperature and moisture, including the effect of curing, are studies that can be done with microwave NDT. In addition to these, many features such as element thickness and location of metallic reinforcement bars and channels, determination of basic structure properties, voids, honeycomb formation, stratification, cracking and moisture ingress control, corrosion of rebars, and rebar diameter can be determined by microwave radar technique [17–19].This technique is based on transferring the EM wave to the examination environment. Although the frequency of the EM wave used in civil engineering applications varies according to the application area, inspections are performed at frequencies in the 30 MHz to 300 GHz ranges. Some parts of the incident wave are transmitted, and the remaining part of the wave is reflected and scattered by the material under test [20]. The returning wave includes information on the materials with different EM properties or defects on the surface and subsurface. Studying the transmission and reflection of radar signals in concrete structures is synonymous with studying the propagation and scattering of EM waves in a dielectric medium, that is, solving Maxwell’s equations in a dielectric medium. Although rebar is one of the indispensables of reinforced concrete structures, it is quite natural that there are cracks and delamination in the aging concrete structures or defects and void in new structures. In the circumstances, concrete to be examined can contain different types of cracks, rebar, delamination, void, and layered structures, so the computation domain becomes a complex geometry and analytical field solutions to Maxwell’s equations are very difficult to find. Fortunately, numerical techniques are a popular way to understand better and interpret and simulate wave propagation in the concrete medium [21–24].Various computational methods like the method of moments (MoM) [25], the finite element method (FEM) [26], and the finite-difference time-domain (FDTD) method [27, 28] have been developed for the solution of Maxwell’s equations. There are also methods that reformulate electromagnetic governing equations by converting local differential equations into nonlocal integral form using variational principles [29, 30]. However, the MoM method does not have numerical dispersion errors that occur in FDTD and finite element methods; since the MoM method is used for linear problems and simple geometries, FDTD and finite element methods are more suitable for modeling complex shapes and inhomogeneous materials [26, 27].Among the various time-domain numerical techniques, the FDTD algorithm is well suited to the analysis of the interaction between concrete walls and EM waves, due to the simplicity of the algorithm, and there is nearly no limitation to the definition of the geometrical properties of the simulated material. The basis of the Yee algorithm is very robust and has great usability. Instead of solving the wave equation for only electric field and magnetic field, a solution is obtained for the electric and magnetic fields by time and space discretization using Maxwell’s curl equations combined with the Yee algorithm. As bothE and H information are used to solve the problem, a better result is obtained in the arbitrarily shaped geometries [27, 28]. Besides deciding the geometrical structure for simulation with FDTD, the computation domain needs to be restricted to prevent the reflection of the wave. Therefore, to truncate the FDTD computation domain and to absorb the waves that come to the boundary of the domain without any reflection, a finite-thickness perfectly matched layer (PML) is used [31–35].The dielectric properties of the construction structures are one of the most important parameters in microwave radar NDT [36–40]. The dielectric properties of the concrete structure change according to the moisture level [41]. However, it is also known that the electrical permittivity changes according to the frequency in microwave measurements [42]. Depending on the moisture level, concrete structures can be having dispersive medium properties. Although the frequency dependence of concrete materials is represented by models such as the extended one-pole Debye model, Cole-Cole model, and Jonscher model, the dispersive FDTD modeling of these models has many computational burdens because they contain fractional derivatives [43]. Using a quadratic complex rational function (QCRF), an accurate FDTD dispersive modeling of concrete materials with different moisture levels can be obtained at different frequency ranges [44]. Since the simulations and measurements in our study are made in 0.4–4.0 GHz frequency ranges, the QCRF-based dispersive FDTD model proposed in [44] is extended to the 0.4–4.0 GHz range by using the complex permittivity of concrete structures measured in these frequency ranges [42].Measurement methods that provide information about the propagation of EM waves according to the dielectric properties of any material are parallel plate capacitor technique, free space technique, transmission line technique, and resonator/oscillator technique. Using the vector network analyzer (VNA), which is a device operating at microwave frequencies developed to process the magnitude and phase of the reflected and transmitted waves from the network and a specific measurement setup in free space, the transmission and reflection coefficients of the tested specimen can be measured, connected with the measurement scheme [45]. The reflection or transmission coefficients of any material or device can be made using measurement systems as well as using numerical methods. In cases where the material properties varied continuously/discontinuously (transition from one half-space to the next), transmission and reflection coefficients can be calculated with the FDTD method using a sufficiently small spatial step size [46–49].In this paper, internal crack that cannot be detected from the outside and two different surface crack geometries are defined in concrete structures. In order to define the surface crack, rectangles of decreasing width and circles of the same radius moving one after another are used, respectively, in two different concrete samples. In the third concrete sample, crack type is defined differently from the void that cannot be detected from the surface and a very thin crack definition is made by using the Riemann integral region to express the concrete defects. FDTD modeling is performed in the center frequency ofL (1.0–2.0 GHz) and S (2.0–4.0 GHz) bands for prepared simulation setups which are suitable for the defined cracks. In the FDTD calculations, the dispersivity of the concrete samples is also taken into account. The structures which have surface crack are measured with Vivaldi antenna array. The reflection and transmission coefficients obtained from the measurement results and FDTD calculations are compared.


## 2. Research Methodology

### 2.1. Maxwell’s Curl Equations and FDTD

Maxwell’s equations express the propagation of the electromagnetic wave. Among Maxwell’s electromagnetic equations, the FDTD method uses two curl equations presented below in Cartesian coordinates for the linear, isotropic, and lossy medium:(1)∂H∂t=−1μ∇×E−1μMsource+σ∗H,(2)∂E∂t=1ε∇×H−1εJsource+σE,where E is the electric field (in volts per meter) and H is the magnetic field (in amperes per meter) with permeability μ, electric conductivity σ, equivalent magnetic loss σ∗, and permittivity ε. The curl equations given in equations (1)-(2) can be decomposed into components and six scalar equations are obtained. When components of these equations are reformulated, two decoupled groups of equations are obtained, and the decoupled equations are evaluated separately. Since the changes in shape or position in the third direction of the geometric models examined here are constant and the incident wave is uniform in the z-direction, the structure to be simulated using FDTD can be considered to be infinite in this direction; all derivatives with respect to z are equal to zero. Six scalar equations reduced to 2D can be divided into two groups. In the first group, all the electric field components are transverse to the reference dimension z; this group of equations is called TEz mode; in the second group, all the magnetic field components are transverse to the reference dimension z and this group of equations is called TMz mode. TEz mode contains Ex,Ey,Hz components and TMz mode contains Hx,Hy,Ezcomponents.These modes may be existing simultaneously in any medium without any interactions. These two modes are very different from each other physically. This is due to the positioning of theE and H field lines relative to the surface of the modeled structure. The TEz mode creates E-field lines in a plane perpendicular to the z-axis with infinite length, while the TMz mode sets up only E-field lines parallel to the z-axis. Since the inspection of concrete specimens is examined in both TEz mode and TMz mode, the equations expressing these modes are given as follows.In the Cartesian coordinate system, 2D TEz can be written as(3a)∂Ex∂t=1ε∂Hz∂y−Jsourcex+σEx,(3b)∂Ey∂t=1ε−∂Hz∂x−Jsourcey+σEy,(3c)∂Hz∂t=1μ∂Ex∂y−∂Ey∂x−Msourcez+σ∗Hz.TMz mode with respect to z in 2D or TMz mode is as follows:(4a)∂Hx∂t=1μ−∂Ez∂y−Msourcex+σ∗Hx,(4b)∂Hy∂t=1μ∂Ez∂x−Msourcey+σ∗Hy,(4c)∂Ez∂t=1ε∂Hy∂x−∂Hx∂y−Jsourcez+σEz.FDTD is a method that provides convergence to Maxwell’s differential equations in time by dividing the calculation region into stepwise rectangular elements based on Yee’s algorithm [27]. Applying the central finite-difference approximation to the spatial and time derivatives of equations (3a)–(3c) and (4a)–(4c) yields the 2D TEz and TMz mode discrete equations. An explicit solution for Maxwell’s curl equations is obtained using the direct numerical implementation [27, 28].The 2D TEz and TMz case represented in equations (3a)–(3c) and (4a)–(4c) is discretized by I×J sample points, also known as nodes, I is the number of points in the x-axis, and J is the number of points in the y-axis. The E and H fields are evaluated at each node in the computational domain. The location of the node in the computational domain is given by x,y=iΔx,jΔy, Δx and Δy are the space increment in the x-direction and in the y-direction, respectively. The i and j indices correspond to the spatial location of the x-direction and y-direction, respectively. In addition, let Ux,y,t be any function of space and time evaluated at a separate point on the grid and at a discrete point in time as Ux,y,t=UiΔx,jΔy,nΔt=Uni,j, where the n-index corresponds to the temporal step and Δt corresponds to the temporal step size. EM field quantities in time steps n+1 are determined using the previous values [27].


### 2.2. The Stability Condition and Incident Waves

The selection of the increment in the discretization of temporal and spatial is important for stability criterion. Unstable solutions, which are likely to occur in solutions of differential equations, are undesirable. If the numerical stability is not provided, the results will be unlimited and irregular increase during the time increment. Since the FDTD algorithm samples both electric field and magnetic field at separate points in time and space, the choice of sampling in time and sampling in space must comply with certain constraints to ensure solution stability. According to accuracy and stability criteria, the computational domain grid size and time increment are based on the electromagnetic properties and source excitation. The simulation code provides the requirement to have at least ten grid points for the shortest wavelength [27]. The stability of the FDTD formulation is achieved by using the Courant condition when the grid size is known in the time [27]. For a uniform discretization scheme Δx=Δy=Δz=Δ, the Courant stability criteria can be expressed more simply:(5)Δt≤Δcn.Sources are indispensable parts of FDTD simulations as they generate electric and magnetic fields and they are divided into two as far field and near field according to the type of problem to be solved. In both cases, a source excites electric and magnetic fields with a waveform as a function of time [28]. The waveform type is specially selected according to the nature of the problem, considering some of the limitations of the FDTD method. A source used in the FDTD grid can be set up simply by assigning a desired time function to certain E and H components in the space grid. This time function can be any sinusoidal or impulse function and is independent of anything in the model. The frequency range in which valid and accurate results can be obtained in FDTD simulation is determined by the frequency spectrum of the source waveform. The frequency spectrum of the source waveform used in FDTD simulation should include all frequencies used in the simulation [28].A Gaussian waveform is the most suitable waveform as it can contain all frequencies down to the highest frequency that depends on a cell size by a factor. In all simulations, the point source electromagnetic wave is generated. The point source in the 2D numerical simulation is assumed since the point source is a spherical wave. A Gaussian waveform can be written as a function of time as(6)gt=exp−t−t0τ01/2,where t0 states how long after the start of the simulation the pulse will be included in the system and τ0 is the parameter determining the pulse width. Since the product of time and time bandwidth is constant, as τ0 gets smaller, the frequency band widens. In the simulations, the source is included in the update equations over the Hz component in the TEz mode and Ezcomponent in the TMz mode at 1.5 GHz and 3.0 GHz frequencies.


### 2.3. Boundary Condition

The problem investigated in this study is governed by the Maxwell equations and it is solved by using the FDTD in a computational domain to understand the interaction of the EM wave and concrete structure. The problem space must be restricted at a certain place where it will be sufficient for interactions and terminated appropriately since it is not possible to simulate the infinite expansion of the problem space and the FDTD calculation domain must be finite. When the scattered and radiating fields reach the boundary, the fields can reflect into the problem space if the necessary boundary conditions are not applied. In numerical calculations of wave propagation, an artificial boundary must be defined in order to truncate the calculation space and prevent the reflection of the outgoing waves. Physically, this process is like the walls of an anechoic chamber. Usually, the boundary of the FDTD problem space is chosen in such a way that the scattered or radiating fields are absorbed when they reach the boundary. In order to model the real infinity, the boundary conditions applied after the problem space is terminated in a specific region are called Absorbing Boundary Conditions (ABCs).The most used ABCs are Mur type ABC, dispersive boundary condition (DBC), and PML. The PML developed by Berenger [31] has proven to be one of the strongest ABCs compared to other methods [32–35]. PML is a finite-thickness lossy material boundary layer that is perfectly matched to the solution space and surrounds the computation domain. The main advantage of PML over other ABCs is that it can absorb all waves of any frequency, polarization, and angle on it ideally with zero reflections. Another advantage of the PML is that the numerical computation domain can be composed of inhomogeneous, anisotropic, dispersive, or nonlinear mediums. Such mediums cannot be modeled analytic ABCs. Although there are various types of PML in the literature, the first PML application is applied to the problem space of FDTD by Berenger. PML is used in FDTD simulations to truncate the FDTD computation space. Berenger PML is based on the principle of obtaining the reflection coefficient as zero by using nonphysical splitting field components in the PML region with electrical and magnetic losses. The PML absorbing medium’s dielectric properties provide the condition given in(7)σε0=σ∗μ0,where σ is the electric conductivity and σ∗ is the magnetic conductivity. This equation ensures that the impedance of the absorbing medium is zero so that the waves from the computation space cannot be reflected in the computation domain after they hit the absorbing medium. In numerical simulation, the absorber medium has a few mesh cell thicknesses [27].


### 2.4. Dielectric Properties of Materials

The (relative to free space) dielectric property of a material is a complex parameter (i.e.,εr=εr′−jεr″) consisting of its permittivity is referred to as dielectric constant (real part) and its loss factor (imaginary part) [20]. The dielectric constant refers to how polarizable a medium is or how much energy it can store when exposed to an electric field. The loss factor represents the medium’s or material’s ability to absorb microwave energy. The dielectric constant of concrete varies depending on the moisture content, measurement temperature, measurement frequency, and water-to-cement ratio [41, 42]. The dielectric constants vary between 5 and 25 according to the wet, saturated, air-dried, and oven-dried conditions of the concrete in the frequency ranges of 0–20 GHz, including the frequency ranges in our study [42]. If the moisture level in the concrete is high, the change of the dielectric constant with frequency becomes essential, which means that the concrete is dispersive. Using the QCRF, FDTD dispersive modeling of concrete materials can be obtained [44]. Since simulations and measurements are made in the frequency range of 0.4–4.0 GHz, the algorithm given in [44] has been extended to this frequency range by using the electromagnetic properties of concrete provided in [42].


## 2.1. Maxwell’s Curl Equations and FDTD

Maxwell’s equations express the propagation of the electromagnetic wave. Among Maxwell’s electromagnetic equations, the FDTD method uses two curl equations presented below in Cartesian coordinates for the linear, isotropic, and lossy medium:(1)∂H∂t=−1μ∇×E−1μMsource+σ∗H,(2)∂E∂t=1ε∇×H−1εJsource+σE,where E is the electric field (in volts per meter) and H is the magnetic field (in amperes per meter) with permeability μ, electric conductivity σ, equivalent magnetic loss σ∗, and permittivity ε. The curl equations given in equations (1)-(2) can be decomposed into components and six scalar equations are obtained. When components of these equations are reformulated, two decoupled groups of equations are obtained, and the decoupled equations are evaluated separately. Since the changes in shape or position in the third direction of the geometric models examined here are constant and the incident wave is uniform in the z-direction, the structure to be simulated using FDTD can be considered to be infinite in this direction; all derivatives with respect to z are equal to zero. Six scalar equations reduced to 2D can be divided into two groups. In the first group, all the electric field components are transverse to the reference dimension z; this group of equations is called TEz mode; in the second group, all the magnetic field components are transverse to the reference dimension z and this group of equations is called TMz mode. TEz mode contains Ex,Ey,Hz components and TMz mode contains Hx,Hy,Ezcomponents.These modes may be existing simultaneously in any medium without any interactions. These two modes are very different from each other physically. This is due to the positioning of theE and H field lines relative to the surface of the modeled structure. The TEz mode creates E-field lines in a plane perpendicular to the z-axis with infinite length, while the TMz mode sets up only E-field lines parallel to the z-axis. Since the inspection of concrete specimens is examined in both TEz mode and TMz mode, the equations expressing these modes are given as follows.In the Cartesian coordinate system, 2D TEz can be written as(3a)∂Ex∂t=1ε∂Hz∂y−Jsourcex+σEx,(3b)∂Ey∂t=1ε−∂Hz∂x−Jsourcey+σEy,(3c)∂Hz∂t=1μ∂Ex∂y−∂Ey∂x−Msourcez+σ∗Hz.TMz mode with respect to z in 2D or TMz mode is as follows:(4a)∂Hx∂t=1μ−∂Ez∂y−Msourcex+σ∗Hx,(4b)∂Hy∂t=1μ∂Ez∂x−Msourcey+σ∗Hy,(4c)∂Ez∂t=1ε∂Hy∂x−∂Hx∂y−Jsourcez+σEz.FDTD is a method that provides convergence to Maxwell’s differential equations in time by dividing the calculation region into stepwise rectangular elements based on Yee’s algorithm [27]. Applying the central finite-difference approximation to the spatial and time derivatives of equations (3a)–(3c) and (4a)–(4c) yields the 2D TEz and TMz mode discrete equations. An explicit solution for Maxwell’s curl equations is obtained using the direct numerical implementation [27, 28].The 2D TEz and TMz case represented in equations (3a)–(3c) and (4a)–(4c) is discretized by I×J sample points, also known as nodes, I is the number of points in the x-axis, and J is the number of points in the y-axis. The E and H fields are evaluated at each node in the computational domain. The location of the node in the computational domain is given by x,y=iΔx,jΔy, Δx and Δy are the space increment in the x-direction and in the y-direction, respectively. The i and j indices correspond to the spatial location of the x-direction and y-direction, respectively. In addition, let Ux,y,t be any function of space and time evaluated at a separate point on the grid and at a discrete point in time as Ux,y,t=UiΔx,jΔy,nΔt=Uni,j, where the n-index corresponds to the temporal step and Δt corresponds to the temporal step size. EM field quantities in time steps n+1 are determined using the previous values [27].


## 2.2. The Stability Condition and Incident Waves

The selection of the increment in the discretization of temporal and spatial is important for stability criterion. Unstable solutions, which are likely to occur in solutions of differential equations, are undesirable. If the numerical stability is not provided, the results will be unlimited and irregular increase during the time increment. Since the FDTD algorithm samples both electric field and magnetic field at separate points in time and space, the choice of sampling in time and sampling in space must comply with certain constraints to ensure solution stability. According to accuracy and stability criteria, the computational domain grid size and time increment are based on the electromagnetic properties and source excitation. The simulation code provides the requirement to have at least ten grid points for the shortest wavelength [27]. The stability of the FDTD formulation is achieved by using the Courant condition when the grid size is known in the time [27]. For a uniform discretization scheme Δx=Δy=Δz=Δ, the Courant stability criteria can be expressed more simply:(5)Δt≤Δcn.Sources are indispensable parts of FDTD simulations as they generate electric and magnetic fields and they are divided into two as far field and near field according to the type of problem to be solved. In both cases, a source excites electric and magnetic fields with a waveform as a function of time [28]. The waveform type is specially selected according to the nature of the problem, considering some of the limitations of the FDTD method. A source used in the FDTD grid can be set up simply by assigning a desired time function to certain E and H components in the space grid. This time function can be any sinusoidal or impulse function and is independent of anything in the model. The frequency range in which valid and accurate results can be obtained in FDTD simulation is determined by the frequency spectrum of the source waveform. The frequency spectrum of the source waveform used in FDTD simulation should include all frequencies used in the simulation [28].A Gaussian waveform is the most suitable waveform as it can contain all frequencies down to the highest frequency that depends on a cell size by a factor. In all simulations, the point source electromagnetic wave is generated. The point source in the 2D numerical simulation is assumed since the point source is a spherical wave. A Gaussian waveform can be written as a function of time as(6)gt=exp−t−t0τ01/2,where t0 states how long after the start of the simulation the pulse will be included in the system and τ0 is the parameter determining the pulse width. Since the product of time and time bandwidth is constant, as τ0 gets smaller, the frequency band widens. In the simulations, the source is included in the update equations over the Hz component in the TEz mode and Ezcomponent in the TMz mode at 1.5 GHz and 3.0 GHz frequencies.


## 2.3. Boundary Condition

The problem investigated in this study is governed by the Maxwell equations and it is solved by using the FDTD in a computational domain to understand the interaction of the EM wave and concrete structure. The problem space must be restricted at a certain place where it will be sufficient for interactions and terminated appropriately since it is not possible to simulate the infinite expansion of the problem space and the FDTD calculation domain must be finite. When the scattered and radiating fields reach the boundary, the fields can reflect into the problem space if the necessary boundary conditions are not applied. In numerical calculations of wave propagation, an artificial boundary must be defined in order to truncate the calculation space and prevent the reflection of the outgoing waves. Physically, this process is like the walls of an anechoic chamber. Usually, the boundary of the FDTD problem space is chosen in such a way that the scattered or radiating fields are absorbed when they reach the boundary. In order to model the real infinity, the boundary conditions applied after the problem space is terminated in a specific region are called Absorbing Boundary Conditions (ABCs).The most used ABCs are Mur type ABC, dispersive boundary condition (DBC), and PML. The PML developed by Berenger [31] has proven to be one of the strongest ABCs compared to other methods [32–35]. PML is a finite-thickness lossy material boundary layer that is perfectly matched to the solution space and surrounds the computation domain. The main advantage of PML over other ABCs is that it can absorb all waves of any frequency, polarization, and angle on it ideally with zero reflections. Another advantage of the PML is that the numerical computation domain can be composed of inhomogeneous, anisotropic, dispersive, or nonlinear mediums. Such mediums cannot be modeled analytic ABCs. Although there are various types of PML in the literature, the first PML application is applied to the problem space of FDTD by Berenger. PML is used in FDTD simulations to truncate the FDTD computation space. Berenger PML is based on the principle of obtaining the reflection coefficient as zero by using nonphysical splitting field components in the PML region with electrical and magnetic losses. The PML absorbing medium’s dielectric properties provide the condition given in(7)σε0=σ∗μ0,where σ is the electric conductivity and σ∗ is the magnetic conductivity. This equation ensures that the impedance of the absorbing medium is zero so that the waves from the computation space cannot be reflected in the computation domain after they hit the absorbing medium. In numerical simulation, the absorber medium has a few mesh cell thicknesses [27].


## 2.4. Dielectric Properties of Materials

The (relative to free space) dielectric property of a material is a complex parameter (i.e.,εr=εr′−jεr″) consisting of its permittivity is referred to as dielectric constant (real part) and its loss factor (imaginary part) [20]. The dielectric constant refers to how polarizable a medium is or how much energy it can store when exposed to an electric field. The loss factor represents the medium’s or material’s ability to absorb microwave energy. The dielectric constant of concrete varies depending on the moisture content, measurement temperature, measurement frequency, and water-to-cement ratio [41, 42]. The dielectric constants vary between 5 and 25 according to the wet, saturated, air-dried, and oven-dried conditions of the concrete in the frequency ranges of 0–20 GHz, including the frequency ranges in our study [42]. If the moisture level in the concrete is high, the change of the dielectric constant with frequency becomes essential, which means that the concrete is dispersive. Using the QCRF, FDTD dispersive modeling of concrete materials can be obtained [44]. Since simulations and measurements are made in the frequency range of 0.4–4.0 GHz, the algorithm given in [44] has been extended to this frequency range by using the electromagnetic properties of concrete provided in [42].


## 3. Description of the Structures

Some of the physical models representing the types of cracks that could be found in any concrete structure are created. The simulation geometries of these crack types are explained in detail under subheadings. All the simulation setups are surrounded by the PML layer in order to truncate the computation domain and prevent the waves coming to the boundary from returning to the domain. Eight point sources are used at the top of the computation domain in all simulation setups.


### 3.1. Case A: Surface Crack Definition Using Sequential Rectangles

Rectangles, whose centers are on the same axis, are used for defining the surface crack. Rectangles of equal height but narrowing downwards are added end to end to form a stair shape. The simulation setup of this surface crack can be expressed as in Figure1. The dimensions of the numerical domain are 300 grids, in other words, 300-unit cells in the x-axis and 150 grids in the y-axis. For the spatial increment in the simulation domain, uniform discretization in space is appealed with Δx=Δy=0.0015m, so the computation domain physically corresponds to a rectangular region of 0.45m×0.225m. The computation domain contains four embedded rebars (radius of rebar is 8 mm).Figure 1
Surface crack model: Case A.


### 3.2. Case B: Surface Crack Definition Using Dynamic Geometry

Another surface crack definition is made by moving the same centered circles on triangle edges. The tangent point of one of the circles is the center of another the successive circle and circles have the same radius. The dimensions of the numerical domain are 300 grids on thex-axis and 150 grids on the y-axis. In reality, the computation domain corresponds to a rectangular region of 0.45m×0.225m and Δx=Δy=0.0015m. The experimental setup of Case B contains four embedded rebars (radius of rebar is 8 mm), expressed in Figure 2.Figure 2
Surface crack model: Case B.


### 3.3. Case C: Interior Crack Description Using Riemannian Integral Definition

Interior concrete crack definition is made taking advantage of Riemann integral domain, in other words, the domain that is under a curve. Firstly, the area under the curve is defined using the Riemann integral expression, and then this definition is extended to the region expression between two curves. The area between the two curves is considered as an air-filled crack in simulations. To simulate this crack, which is the air-filled area, the integral of the function assumed as a crack curve is used. The Riemann integral of a function requires the calculations of the area by dividing it into rectangles for this purpose. If one rectangle is used, the rudest approximation is achieved. But approximating the actual value of the area between two curves and making better refinements to the area, more rectangles should be added. These statements can be expressed mathematically as such if the number of subintervals is increased, a more realistic area definition is made.The subintervals are completed to rectangles from the bottom as in Figures3(a) and 3(c) or top as in Figures 3(b) and 3(d) after dividing the area under the function curve into subintervals. Figures 3(a) and 3(b) are the rudest approximation, Figures 3(c) and 3(d) are a little better, but Figures 4(a) and 4(b) are the best since calculating the area under the curve n piece of the rectangle is considered [50].Figure 3
Riemann integral definition: rude partition.
(a)(b)(c)(d)Figure 4
Riemann integral definition: thin partition.
(a)(b)As a result of these definitions, in Case C, the most realistic crack explanation that can be encountered in any concrete structure is made and this type of crack refers to a branching internal crack caused by various external factors. The crack expressed in this simulation model is an internal crack and cannot be physically seen from the outside. For Case C, the dimensions of the computational domain are 300 grids in thex-axis and 150 grids in the y-axis, and this physically corresponds to a rectangular region of 0.45m×0.225m and Δx=Δy=0.0015m. 0.045 m of the computation domain is air; the rest of the domain, in other words, thickness of the concrete, is 0.18 m. The simulation experiment for Case C is given in Figure 5.Figure 5
Interior crack model: Case C.The curve expressed by the function that represents the concrete crack or delamination physically is divided into subintervals. These subintervals have completed the rectangle from the top and the bottom. By doing so, rectangular subdomains are obtained, and then electric field and magnetic fields are computed at all these subdomains with FDTD.


## 3.1. Case A: Surface Crack Definition Using Sequential Rectangles

Rectangles, whose centers are on the same axis, are used for defining the surface crack. Rectangles of equal height but narrowing downwards are added end to end to form a stair shape. The simulation setup of this surface crack can be expressed as in Figure1. The dimensions of the numerical domain are 300 grids, in other words, 300-unit cells in the x-axis and 150 grids in the y-axis. For the spatial increment in the simulation domain, uniform discretization in space is appealed with Δx=Δy=0.0015m, so the computation domain physically corresponds to a rectangular region of 0.45m×0.225m. The computation domain contains four embedded rebars (radius of rebar is 8 mm).Figure 1
Surface crack model: Case A.


## 3.2. Case B: Surface Crack Definition Using Dynamic Geometry

Another surface crack definition is made by moving the same centered circles on triangle edges. The tangent point of one of the circles is the center of another the successive circle and circles have the same radius. The dimensions of the numerical domain are 300 grids on thex-axis and 150 grids on the y-axis. In reality, the computation domain corresponds to a rectangular region of 0.45m×0.225m and Δx=Δy=0.0015m. The experimental setup of Case B contains four embedded rebars (radius of rebar is 8 mm), expressed in Figure 2.Figure 2
Surface crack model: Case B.


## 3.3. Case C: Interior Crack Description Using Riemannian Integral Definition

Interior concrete crack definition is made taking advantage of Riemann integral domain, in other words, the domain that is under a curve. Firstly, the area under the curve is defined using the Riemann integral expression, and then this definition is extended to the region expression between two curves. The area between the two curves is considered as an air-filled crack in simulations. To simulate this crack, which is the air-filled area, the integral of the function assumed as a crack curve is used. The Riemann integral of a function requires the calculations of the area by dividing it into rectangles for this purpose. If one rectangle is used, the rudest approximation is achieved. But approximating the actual value of the area between two curves and making better refinements to the area, more rectangles should be added. These statements can be expressed mathematically as such if the number of subintervals is increased, a more realistic area definition is made.The subintervals are completed to rectangles from the bottom as in Figures3(a) and 3(c) or top as in Figures 3(b) and 3(d) after dividing the area under the function curve into subintervals. Figures 3(a) and 3(b) are the rudest approximation, Figures 3(c) and 3(d) are a little better, but Figures 4(a) and 4(b) are the best since calculating the area under the curve n piece of the rectangle is considered [50].Figure 3
Riemann integral definition: rude partition.
(a)(b)(c)(d)Figure 4
Riemann integral definition: thin partition.
(a)(b)As a result of these definitions, in Case C, the most realistic crack explanation that can be encountered in any concrete structure is made and this type of crack refers to a branching internal crack caused by various external factors. The crack expressed in this simulation model is an internal crack and cannot be physically seen from the outside. For Case C, the dimensions of the computational domain are 300 grids in thex-axis and 150 grids in the y-axis, and this physically corresponds to a rectangular region of 0.45m×0.225m and Δx=Δy=0.0015m. 0.045 m of the computation domain is air; the rest of the domain, in other words, thickness of the concrete, is 0.18 m. The simulation experiment for Case C is given in Figure 5.Figure 5
Interior crack model: Case C.The curve expressed by the function that represents the concrete crack or delamination physically is divided into subintervals. These subintervals have completed the rectangle from the top and the bottom. By doing so, rectangular subdomains are obtained, and then electric field and magnetic fields are computed at all these subdomains with FDTD.


## 4. Numerical Simulation Results

The numerical simulation experiments are excited by TE and TM waves, and all electric and magnetic components are calculated for Case A to Case C.Ex field distribution from simulations using TE wave and Ez field distributions from simulations using TM wave are shown and compared. In the presented cases, the maximum wave velocity is equal to the speed of light propagated along with the free space. The simulation results are obtained at 1.5 GHz and 3.0 GHz frequencies. Figures on how the field distributions changed in different time steps for each case are obtained, and the results in which cracks and rebars are best distinguished are presented. The operation frequencies for simulations are 1.5 GHz and 3.0 GHz and some of the results at specific times are given.Simulation results are given for Case A as follows.Four different simulation results are given in Figure6 for Case A. The thickness of the concrete and free space is 0.18 m and 0.045 m, respectively. The simulation results for 1.5 GHz colors are more distinct, so crack and rebar are seen clearly in Figure 6(a). When we compare the frequencies in the same time step for the Ex field distribution, the 1.5 GHz frequency is more suitable than the 3.0 GHz frequency. Figure 6(a) is a more reasonable result than Figure 6(c). Ez field distribution obtained from TM mode at 1.5 GHz and 3.0 GHz is given, respectively, in Figures 6(b) and 6(d). The variation of the waves propagating from the air to the concrete appears to be more apparent in TM mode. On the other hand, rebars can be easily viewed in any situation.Figure 6
Simulation result for Case A at 600-time step: (a)Ex field distribution for 1.5 GHz, (b) Ez field distribution for 1.5 GHz, (c) Ex field distribution for 3.0 GHz, and (d) Ez field distribution for 3.0 GHz.Simulation results are given for Case B as follows.An equal number of simulation results with Case A is also given for Case B. The thickness of the concrete and free space is 0.18 m and 0.045 m, respectively. Although the rebars can be displayed in all results, the best result is obtained in Figure7(a). The transmission of the EM wave from one medium to another can be observed better in the field distribution obtained from TM mode. Since the crack in Case B is closer to the surface than in Case A and the TM mode is more suitable for observing the transition of the wave from the air to the concrete, the surface crack appears to be more appropriate in Figure 7(b). Although the irons can be distinguished well in Figure 7(c), the surface crack is not obvious. In Figure 7(d), although crack can be seen, the appearance of the rebars is not as clear as other shapes.Figure 7
Simulation result for Case B at 600-time step: (a)Ex field distribution for 1.5 GHz, (b) Ez field distribution for 1.5 GHz, (c) Ex field distribution for 3.0 GHz, and (d) Ez field distribution for 3.0 GHz.Simulation results are given for Case C as follows.Six different simulation results are given for Case C. The thickness of the concrete and free space is 0.18 m and 0.045 m, respectively. In this experiment, rebar is not used to prevent the scattering of the wave in the structure so internal cracks are clearly visualized. The best result is in Figure8(a). In Figure 8(b), the crack continuing to the right can be observed better. From the results of Figures 8(c) and 8(d), it can be concluded that increasing frequency makes internal crack detection difficult.Figure 8
Simulation result for Case C at 700-time step: (a)Ex field distribution for 1.5 GHz, (b) Ez field distribution for 1.5 GHz, (c) Ex field distribution for 3.0 GHz, and (d) Ez field distribution for 3.0 GHz.In Figure9, 3D simulation results at 1.5 GHz for Case A, Case B, and Case C are given, respectively. Since the most suitable frequency for both crack and rebar detection is 1.5 GHz, it is preferred to provide 3D figures of simulations at this frequency. The locations of the sources can be seen clearly in all results. In Case A, because the crack is deeper, the wave can propagate without attenuation too much, hitting the rectangular-shaped edges and creating a reflection. This situation can be seen in Figure 9(a). While a similar situation can be observed for Case C as in Figure 9(e), it could not be observed for Case B in Figure 9(c). Since the EM wave is scattered from rebars, they appear prominently in Figures 9(a)–9(d). As we know from the 2D simulation results, the air-concrete separation can be determined more clearly in the Ez field distribution in Figures 9(b), 9(d), and 9(f).Figure 9
3D simulation result for all cases for 1.5 GHz: (a)Ex for Case A at 600-time step, (b) Ez for Case A at 600-time step, (c) Ex for Case B at 600-time step, (d) Ez for Case B at 600-time step, (e) Ex for Case C at 700-time step, and (f) Ez for Case C at 700-time step.


## 5. Measurement Results

One of the most valuable devices in the microwave NDT technique is the microwave network analyzer with one or more ports and each port can pass, absorb, and reflect electromagnetic energy. The advantages of this device are that it has the flexibility required to test the use of microwave NDT in specific applications; it contains microprocessors that provide calibration and customization for each application, which make it easy to use.The parameters that define the electrical properties of the linear electrical network stimulated by electrical signals are called scattering parameters. While scattering parameters are measured using open- or short-circuit conditions to the linear electrical network, S-parameters are measured in terms of power using matched loads and they are easier to use at higher signal frequencies [49]. S-parameters are the most common way of representing VNA measurements, and they can be measured directly with the VNA. The S-parameters obtained from the VNA measurement are usually in the form of a symmetric matrix with rows and columns equal to the number of ports. For the S-parameter Sij, the subscript i indicates the exciting port (input port) and the subscript j denotes the output port. Diagonal elements of the S-matrix represent reflection and off-diagonal elements of this matrix are used to explain transmission from port j to port i. For example, S11 is the reflection seen looking into port 1. S21 can be found by applying an incident wave at port 1 and measuring the outcoming wave at port 2. This is equivalent to the transmission from port 1 to port 2. Because of the symmetry, S21 is equivalent to S12. In order to calculate the reflection and scattering coefficients of our samples at microwave frequencies, VNA designed for simultaneous measurement of S-parameters of networks with more than two ports is used. The physical structures of the surface cracks given as Case A and Case B and measurement setup belong to these specimens as shown in Figures 10 and 11, respectively.Figure 10
Physical structure: (a) Case A; (b) Case B.
(a)(b)Concrete specimens are prepared in accordance with the simulation setup given in Case A and Case B. In Figure10(a) a defect is created on the concrete surface by using rectangles that narrow down width but remain constant in height and in Figure 10(b), using successive circles, a wider crack is created. Four rebars are used in both specimens.In Figure11, measurement setup is given. Measurement is made in an anechoic chamber using Vivaldi antenna array, VNA, computer, and specimen under test (SUT). Eight Vivaldi antennas are placed in contact with the defect layer of concrete. Since the concrete is placed in contact with eight antennas in the measurement setup, eight sources are used in the simulations. Measurements are made in the 0.4–4.0 GHz frequency range. The measurement and FDTD calculation results are presented in Figure 12.Figure 11
Measurement setup: (a) Case A; (b) Case B.
(a)(b)Figure 12
Measurement results: (a) S-parameters for Case A; (b) S-parameters for Case B; (c) reflection coefficient for Case A; (d) reflection coefficient for Case B, (e) transmission coefficient for Case A, and (f) transmission coefficient for Case B.The S11 and S12 parameters for Case A and Case B are presented in Figures 12(a) and 12(b), respectively. The reflection and transmission coefficients are calculated using the S-parameters obtained from the measurement by coefficient calculation formula [48]. In addition, the reflection and transmission coefficients are computed with the FDTD method [45], and the results are compared in Figures 12(c)–12(f). It is seen that the measurement results and the FDTD results are overlapped.


## 6. Conclusion

Three different types of cracks that may be encountered in a concrete structure are defined. Rectangles are used one after another, deepening and shrinking downward for defining a surface crack. Circles are utilized to identify defects due to a wide breakage of the concrete surface. An interior crack whose depth and structure cannot be determined from the surface is considered as the Riemann integral domain. This domain is defined as the region between two curves as in the integral expression and is divided into smaller subdomains. The crack is expressed by assuming that these subregions are filled with air. Since the region formed by these subregions does not have a smooth shape, the most realistic crack definition is made.Maxwell’s equations are discretized in space and time using the FDTD algorithm to simulate how the electromagnetic wave propagates in concrete which includes rebar and cracks defined above. In simulations, it is investigated which mode is more advantageous in which situations by comparing the field distributions from both TE mode and TM mode. While TM mode is more suitable for observing the propagation of the EM wave as it propagates from one medium to a different medium and to detect less deep cracks, TE mode gives better results in determining the shapes of the deep cracks and the positions of the rebars. Surface cracks, rebars, and internal cracks are simulated at 1.5 GHz and 3.0 GHz frequencies; 2D and 3D states of the field distributions are given. In the simulation, as the frequency increases, the scattering of the wave also increases, making it difficult to detect defects and rebars. For this reason, it has been observed that the crack and rebar detection simulation give better results around the 1.5 GHz frequency.In accordance with the simulation setup, two air-dried concrete specimens containing surface cracks and rebar are measured nondestructively using Vivaldi antenna array and VNA. The reflection and transmission coefficients expressed by the S-parameters obtained from the measurement and the reflection and transmission coefficients calculated from the FDTD are compared and validated. In concrete structures with different crack types, the coefficients calculated by FDTD and the coefficients obtained from the measurement are overlapped. The reflection and transmission coefficients for narrow-deep (Case A) and wide-shallow (Case B) surface cracks in air-dried concrete specimens molded using concrete with the same content are approximately similar. Since the experimental processes require cost and time, it has been observed in our study that it is possible to model dispersive media using FDTD and give reliable results.



---

*Source: 6624982-2021-06-01.xml*