Time Reverse Modeling of Damage Detection in Underwater Concrete Beams Using Piezoelectric Intelligent Modules

Underwater cracks in concrete structures are often difficult to detect due to their complexity of the service environment. With numerical and experimental analysis of concrete beams immersed in water, an active monitoring system, based on a cement-based piezoelectric intelligent module array (CPIMA), was developed to locate and quantify the underwater cracks. Time reversal (TR) of the stress wave field is accomplished to focus on the crack area through the concrete beam specimen by the system. First, a piezoelectric actuator is applied to emit the initial propagating wave, which can be reflected, attenuated, and diffracted by the crack, transmitted through water filled in the crack, as well as diffracted by the coarse aggregates. To extract the damage waveforms associated with the crack and analyze the robust time-reversal invariance under the high-order multiple scattering effect, a pair of homogeneous and heterogeneous forward finite element (FE) models is established. Then, the damage waveforms are time-reversed and re-propagated in the inverse numerical model, where an optimal refocusing is achieved on the crack that behaves as an acoustic source. Finally, the damage area is obtained in the form of the stacked energy distribution of each time step. The focus results are represented by cloud images and compared with root-mean-square deviation (RMSD) values. Numerical simulation and experiments show that this method can identify and quantify underwater cracks effectively.


Introduction
As a porous solid colloid, concrete usually suffers from micro cracks resulting from pore water, pore air, and chemical admixtures in the early stages of construction. For underwater concrete structures, such as dams, pipelines, and bridge piers, these micro cracks may propagate to macro cracks under cyclic environmental loads, creep, and shrinkage [1,2]. To evaluate the damage degree of concrete structures and arrange remedial measures, it is often necessary to locate and quantify the depth of these cracks. In order to detect underwater cracks, the traditional optical image detection techniques have been used for decades. However, due to the loss of information (scattering, absorption, and convolution effects) in underwater conditions, the quality of underwater images is usually poor and the crack characteristics are difficult to be extracted. Other technologies, such as electrical impedance tomography [3], ground-penetrating radar testing [4], and vibration-based technique [5,6], are time, energy, and financially consuming, and they are easily affected by the noise caused by the complex underwater service environment and the heterogeneity of concrete. both the piezoelectric module and wave propagating numerical models to analyze the relationship between RMSD and different cracks with different sizes, positions, and directions. The results of both the simulation and experiment show the sensitivity of the wave energy attenuation method to the presence of concrete cracks. Unlike the scattered waves, however, the attenuated waves have seldom been used for crack imaging in previous studies since they it is hard to capture them directly.
To image underwater cracks, the time-reversal method (TRM) was herein introduced [28]. TRM is used to converge waves synchrony at the source (and possibly reflected or scattered by the target of interest) [29], which is widely applied in underwater communication, medical imaging, and nondestructive testing [30][31][32]. The beauty of this method is that it can retrace complex multipaths of damage signals through a high-order scattering medium, which effectively improves the signal-to-noise ratio of underwater concrete problems [29]. Prada C. et al. [33] proposed a decomposition method of the time reversal operator to achieve accurate energy focusing on the position of multiple point targets. Based on this, Ulrich T. J. et al. applied the time-inelastic nonlinear diagnosis (TREND) technique to locate damages at different depths in a concrete specimen [34,35]. Kocur G. K. [36] applied TRM and signal-based acoustic emission (AE) technology to locate concrete cracks. This method can use the waveform of received signals to locate the AE activity from cracks, and neglect the pre-selection process. Although TRM has been widely used to identify damage in concrete structures, it has rarely been reported to locate underwater concrete cracks. The accuracy of traditional diffracted wave-based TRM can suffer from being influenced by the transmitted wave through water in the crack. First, the transmitted wave through water creates direct wave paths from the actuator to the sensors, which can be retraced by backpropagating waves to focus on the actuator, rather than the crack. Second, as a wave-based method, the accuracy of TRM depends on the sharpness and strength of the wave source, which can be weakened by the interference of the transmitted wave through water. To overcome this problem, the damage-related waves (including reflected, attenuated, and diffracted waves) are proposed to make the crack a stronger signal source than the crack tip as a diffractor. The multiple paths of different wave components can be retraced by TRM; in this way, a constructive interference can be achieved on the wave source due to time-reversal invariance and the reciprocity theorem. In addition, the proposed damage-related wave considers the transmitted waves through water as a part of the attenuated wave component, in which way the negative effect of water in the crack is reduced.
In this paper, a TRM-based imaging technique is developed to detect underwater concrete cracks with the help of a piezoelectric-based monitoring system. The application of TRM is improved by using the damage-related waves and the crack information is better used. The feasibility and robustness of the proposed method is studied by the forward and inverse heterogeneous FE simulation, respectively. In the experiment, a group of concrete specimens with artificial cracks were immersed in water and tested by the CPIMA-based active monitoring system. Finally, the imaging results are discussed and compared with the energy-based RMSD results.

Piezoelectric Intelligent Module
Piezoelectric intelligent modules are the basic monitoring components in this method. Since the piezoelectric patches are composed of a thin brittle material, they are encapsulated in cement as intelligent modules to enhance the durability and the compatibility with concrete [37]. In this paper, the intelligent module is made of an epoxy cement protection of 30 mm (diameter) × 30 mm (height), with a circular lead zirconate titanate (PZT) patch of a diameter of 20 mm (diameter) × 1 mm (thickness) inside (see Figure 1c). The production process is as follows: first, the PZT patch was placed in a self-made mold (see Figure 1a). Then, the epoxy cement was stirred evenly and poured into the mold. After, the mold was continuously vibrated and vacuumized to reduce the porosity of the cement. Finally, the obtained sample was placed in the standard curing chamber for 7 days, with the control temperature at (20 ± 1) • C and the relative humidity ≥90%.

Wave-Based SHM Method
Because of the direct and inverse piezoelectric effect, the piezoelectric intelligent modules can be used as sensors, which can record real-time stress waves as voltage signals, and can also function as actuators to stimulate modulated signals. For the self-made piezoelectric intelligent module, the stimulated signal would propagate in the form of both P-waves and S-waves through solid materials. The P-wave can travel across a long distance along the polarization direction. Since the good orientation of the P-wave can help to form more effective propagating paths between CPIMA, the Pwave is herein simulated. If there is a crack in this propagation path, the traveling wave will be blocked, reflected, and diffracted. According to Huygens' principle, all these influences behave as secondary wave fronts, which originate from the crack, see Figure 2. Concrete is a typical heterogeneous material since it is composed of aggregates, mortar, and cement matrix. When waves are propagating through such heterogeneous mediums, noises are generated due to refraction, diffraction, and multiple scattering. The random spatial distributed physical composition of concrete also leads to the inhomogeneities of sound speed, which can distort and redirect the traveling wave. Previous studies have shown how an optimal focusing is achieved through inhomogeneous media, using the time-reversal invariance of wave equations. The TR focusing process acts as a spatial-temporal matching filter [38]. Its self-adaptive character can compensate for the distortion caused by the propagation through heterogeneous media.
The aggregate of concrete can be approximately equivalent to a scatterer, which results in information loss of the signals whose wavelengths are smaller than the classical diffraction limit (diameter-to-wavelength ratio of 0.5). Given the P-wave speed of 4000 m/s, when the governing frequency is set as 150 kHz, the calculated wavelength is 26.7 mm. Considering a maximum concrete aggregate size dmax = 5mm (diameter-to-wavelength ratio of 0.19 < 0.5), low scattering herein is expected. The numerical simulation also shows the feasibility of robust focusing through inhomogeneous media when the propagating wavefield generated by the crack is dominant over the noise.

Wave-Based SHM Method
Because of the direct and inverse piezoelectric effect, the piezoelectric intelligent modules can be used as sensors, which can record real-time stress waves as voltage signals, and can also function as actuators to stimulate modulated signals. For the self-made piezoelectric intelligent module, the stimulated signal would propagate in the form of both P-waves and S-waves through solid materials. The P-wave can travel across a long distance along the polarization direction. Since the good orientation of the P-wave can help to form more effective propagating paths between CPIMA, the P-wave is herein simulated. If there is a crack in this propagation path, the traveling wave will be blocked, reflected, and diffracted. According to Huygens' principle, all these influences behave as secondary wave fronts, which originate from the crack, see Figure 2.

Wave-Based SHM Method
Because of the direct and inverse piezoelectric effect, the piezoelectric intelligent modules can be used as sensors, which can record real-time stress waves as voltage signals, and can also function as actuators to stimulate modulated signals. For the self-made piezoelectric intelligent module, the stimulated signal would propagate in the form of both P-waves and S-waves through solid materials. The P-wave can travel across a long distance along the polarization direction. Since the good orientation of the P-wave can help to form more effective propagating paths between CPIMA, the Pwave is herein simulated. If there is a crack in this propagation path, the traveling wave will be blocked, reflected, and diffracted. According to Huygens' principle, all these influences behave as secondary wave fronts, which originate from the crack, see Figure 2. Concrete is a typical heterogeneous material since it is composed of aggregates, mortar, and cement matrix. When waves are propagating through such heterogeneous mediums, noises are generated due to refraction, diffraction, and multiple scattering. The random spatial distributed physical composition of concrete also leads to the inhomogeneities of sound speed, which can distort and redirect the traveling wave. Previous studies have shown how an optimal focusing is achieved through inhomogeneous media, using the time-reversal invariance of wave equations. The TR focusing process acts as a spatial-temporal matching filter [38]. Its self-adaptive character can compensate for the distortion caused by the propagation through heterogeneous media.
The aggregate of concrete can be approximately equivalent to a scatterer, which results in information loss of the signals whose wavelengths are smaller than the classical diffraction limit (diameter-to-wavelength ratio of 0.5). Given the P-wave speed of 4000 m/s, when the governing frequency is set as 150 kHz, the calculated wavelength is 26.7 mm. Considering a maximum concrete aggregate size dmax = 5mm (diameter-to-wavelength ratio of 0.19 < 0.5), low scattering herein is expected. The numerical simulation also shows the feasibility of robust focusing through inhomogeneous media when the propagating wavefield generated by the crack is dominant over the noise. Concrete is a typical heterogeneous material since it is composed of aggregates, mortar, and cement matrix. When waves are propagating through such heterogeneous mediums, noises are generated due to refraction, diffraction, and multiple scattering. The random spatial distributed physical composition of concrete also leads to the inhomogeneities of sound speed, which can distort and redirect the traveling wave. Previous studies have shown how an optimal focusing is achieved through inhomogeneous media, using the time-reversal invariance of wave equations. The TR focusing process acts as a spatial-temporal matching filter [38]. Its self-adaptive character can compensate for the distortion caused by the propagation through heterogeneous media.
The aggregate of concrete can be approximately equivalent to a scatterer, which results in information loss of the signals whose wavelengths are smaller than the classical diffraction limit (diameter-to-wavelength ratio of 0.5). Given the P-wave speed of 4000 m/s, when the governing frequency is set as 150 kHz, the calculated wavelength is 26.7 mm. Considering a maximum concrete aggregate size d max = 5 mm (diameter-to-wavelength ratio of 0.19 < 0.5), low scattering herein is expected. The numerical simulation also shows the feasibility of robust focusing through inhomogeneous media when the propagating wavefield generated by the crack is dominant over the noise.

Principle of the Time Reversal Method
The principle of TRM in a concrete beam basically contains two parts: forward simulation and inverse simulation. The forward simulation aims to extract the damage-related signal to be focused, which is used as the input signal in the following reversal simulation. As shown in Figure 3, a transducer T and n receivers R 1 , R 2 , ..., R n are embedded in the test specimen, with a crack divided into m point-like areas D 1 , D 2 , ..., D m . The propagation wave is emitted by a transducer T, then the damage-related wave is extracted by differentiating the wave field of the damaged beam from the healthy one, which originates from the crack. The inverse simulation is a re-emission process. The damage-related signal is time reversed and imported into the numerical model to re-emit from receivers R 1 , R 2 , ..., R n . As a result, the damage-related signals are focused on each point-like area D i at different time steps. Then, the whole crack area is obtained by accumulating the focus results by time.

Principle of the Time Reversal Method
The principle of TRM in a concrete beam basically contains two parts: forward simulation and inverse simulation. The forward simulation aims to extract the damage-related signal to be focused, which is used as the input signal in the following reversal simulation. As shown in Figure 3, a transducer T and n receivers R1, R2,...,Rn are embedded in the test specimen, with a crack divided into m point-like areas D1, D2,...,Dm. The propagation wave is emitted by a transducer T, then the damagerelated wave is extracted by differentiating the wave field of the damaged beam from the healthy one, which originates from the crack. The inverse simulation is a re-emission process. The damagerelated signal is time reversed and imported into the numerical model to re-emit from receivers R1, R2,...,Rn. As a result, the damage-related signals are focused on each point-like area Di at different time steps. Then, the whole crack area is obtained by accumulating the focus results by time.
where ⊗ is the convolution operator. H (t, Di→Rj) is the transfer function from the damage point Di to the receiver Rj.
If the signal is reversed and re-emitted in the time domain, the received signal at any point P in the control domain is: For the wave source where P = Di, according to the symmetry of the transfer function that H (t, Rj→Di) = H (t, Di→Rj), a constructive temporal coherence is obtained as follows: Equations (2) and (3) can be expressed in the frequency domain: where * denotes the complex conjugate of the function, and ω is the angular frequency. yP (ω), yDi (ω), and xi (ω) are the Fourier transforms of the signals yP (t), yDi (t), and xi (t), respectively. Let the whole crack Di (I = 1,2,...,m) emit the damage signals and the receiver array Rj (j = 1,2,...,n) re-emit the input signals. The output signal at a given damage area Dk (k = 1,2,...,m) can be expressed as: The input signal x i (t) is generated from damage point D i , then the received signal at R j in the time domain is: where ⊗ is the convolution operator. H (t, D i →R j ) is the transfer function from the damage point D i to the receiver R j . If the signal is reversed and re-emitted in the time domain, the received signal at any point P in the control domain is: For the wave source where P = D i , according to the symmetry of the transfer function that H (t, R j →D i ) = H (t, D i →R j ), a constructive temporal coherence is obtained as follows: Equations (2) and (3) can be expressed in the frequency domain: where * denotes the complex conjugate of the function, and ω is the angular frequency. y P (ω), y Di (ω), and x i (ω) are the Fourier transforms of the signals y P (t), y Di (t), and x i (t), respectively. Let the whole crack D i (I = 1, 2, ..., m) emit the damage signals and the receiver array R j (j = 1, 2, ..., n) re-emit the input signals. The output signal at a given damage area D k (k = 1, 2, ..., m) can be expressed as: From the above formula, we can find that for all the damage signals diverging from source D k , their reversed waves can only converge in the original source D k among all the complex paths, while a poor spatial coherence is obtained in the other area due to the interference. This phenomenon can also be explained by the time-reversal invariance of the wave equations.

Wavelet Packet Decomposition and Damage Index
Wavelet packet analysis is an effective tool in signal processing, which can decompose a time-domain signal into different frequency components and evaluate the sensitivity of different frequency bands to damage. In this paper, the wavelet packet is introduced to calculate the energy of the signal, and the damage index is obtained to analyze the energy attenuation of cracks of different depths, and then to be compared with the imaging results.
The sensor signal S is decomposed into 2 n signal sets as {X 1 , X 2 , . . . , X 2 n }, where each signal can be presented as where j is the frequency range (j = 1, 2, . . . , 2 n ), m is the number of samples contained in the signal, and n is the level of wavelet packet decomposition, where n = 3 is herein assumed. The energy of each signal is calculated as: The energy of the received signal of the healthy beam E h and of the damage beam E d can be expressed as: The energy-based damage index can be presented in the form of the root-mean-square deviation (RMSD) as: Sensors 2020, 20, 7318 7 of 19

Property and Modeling Principle of Piezoelectric Patch
For linear elastic materials, the mechanical and electrical behavior is governed by the following piezoelectric constitutive equations in the matrix form: where σ and ε are the mechanical stress and strain tensor, D is the vector of electric displacement, C E is the elastic stiffness matrix, E is the electric field intensity vector, e is the piezoelectric stress factor matrix, and ξ ε is the dielectric constant matrix. When the piezoelectric patch is exposed to the coupled stress-electric field, the equilibrium equation of the piezoelectric material is derived using the virtual displacement principle. In FEM, displacements and electric potentials are interpolated by nodal quantities, then the equilibrium equation is derived in the following matrix form: where K uu , K ϕϕ , and K uu are the mechanical stiffness, dielectric constant, and piezoelectric coupling matrix, respectively. F u is the mechanical force and ρ ϕ is the electric charge volume density. Due to the high sensitivity and low cost, the PZT-5A patch is selected as the core transducer of a CPIM, the dielectric property of which is defined by the permittivity values D 11 of 424, D 22 of 424, and D 33 (in polar direction) of 450. The mechanical property is presented by the stiffness matrix as follows: The tie constraint surface contact is used to define the contact between the concrete and the piezoelectric patch. More information of the piezoelectric patch model is presented in [26].

Random Aggregates Model
To describe the heterogeneity of concrete components (such as the geometry, elastic properties, and spatial distribution), a random aggregate model is established by the secondary development of ABAQUS software based on Python. The random aggregate structure consists of two parts: the randomly distributed concrete aggregates and mortar filling the remaining area. The model is generated by the following steps: Step 1-estimate the particle size distribution. The random sampling principle of the Monte Carlo method is applied to obtain the particle size distribution in the simulation, following the grading curve with a maximum grain size of 5 mm.
Step 2-determine the geometry of each aggregate. First, a circular particle size is generated with the size from step 1. Then, it is cut into a polygon, whose side length and inner angle are controlled by the irregularity parameter.
Step 3-release the aggregates into the beam. Before placing the aggregates, the beam should be meshed into structured grids. The grid size is considered to be 1mm to achieve a good balance of image resolution and computing cost. The aggregates are then randomly mapped to the grid in the order from large ones to small ones. It is also necessary to check the potential problems, such as aggregate overlap, boundary intersection, and gaps between cement faces. If the above problems occur, the aggregates should be re-dropped.
An obtained random aggregates model is presented in Figure 4. of image resolution and computing cost. The aggregates are then randomly mapped to the grid in the order from large ones to small ones. It is also necessary to check the potential problems, such as aggregate overlap, boundary intersection, and gaps between cement faces. If the above problems occur, the aggregates should be re-dropped.
An obtained random aggregates model is presented in Figure 4.

Modeling of Wave Propagation
The dynamic explicit algorithm adopts the time difference scheme of the dynamic equation. It is effective to solve non-linear dynamic problems, such as wave propagation modeling in an inhomogeneous medium, due to the fact that it does not need to iterate or to solve the tangent stiffness directly. The stability of the explicit model depends on the time increment Δt, which should be small enough to satisfy the Neumann criterion. The quadrangle unit is herein used and the element size has a great influence on the calculation accuracy and computing cost. Based on the simulation theory, the number of finite elements n per wavelength should not be less than 10 [39]. Supposing that the wavelength is 26.7 mm, the structured element size is chosen as 1 mm (26.7 elements per wavelength). The forward mesh method is adopted, and the time increment is set as 0.2 μs. To reflect the change of the wave field caused by crack rather than meshes, the elements remain the same for the damaged and healthy beam, while the crack and PZTs are modeled by assigning corresponding materials to prescribed finite elements. The material parameters are shown in Table 1.  For the healthy state shown in Figure 5a, the wave fronts travel horizontally and undisturbed, and the corresponding displacements at nodes of CPIMA are used as the baseline. In Figure 5b, we can see the obvious attenuation of the propagating wave caused by the crack, as well as the slight "leakage" of the propagating wave caused by water filling in the crack as a media. The comprehensive influence of the crack on the wave field can be seen in Figure 5c, which is a snapshot

Modeling of Wave Propagation
The dynamic explicit algorithm adopts the time difference scheme of the dynamic equation. It is effective to solve non-linear dynamic problems, such as wave propagation modeling in an inhomogeneous medium, due to the fact that it does not need to iterate or to solve the tangent stiffness directly. The stability of the explicit model depends on the time increment ∆t, which should be small enough to satisfy the Neumann criterion. The quadrangle unit is herein used and the element size has a great influence on the calculation accuracy and computing cost. Based on the simulation theory, the number of finite elements n per wavelength should not be less than 10 [39]. Supposing that the wavelength is 26.7 mm, the structured element size is chosen as 1 mm (26.7 elements per wavelength). The forward mesh method is adopted, and the time increment is set as 0.2 µs. To reflect the change of the wave field caused by crack rather than meshes, the elements remain the same for the damaged and healthy beam, while the crack and PZTs are modeled by assigning corresponding materials to prescribed finite elements. The material parameters are shown in Table 1.  diffuses around from the crack tip. To capture all the signal components, the modules are suggested to be as more as possible and embedded along the boundary of the monitoring region to capture waves in all directions from the crack.
It can be seen from Figure 6a-c that the inhomogeneities of the media can distort and scatter the wave field. However, this distortion cannot redirect the wave fronts and the scattering effect is not strong enough to overwhelm the wave path from the wave source (crack) to the piezoelectric sensors. The wave reversibility to the crack is valid and hence TRM is adaptive for concrete media.   to be as more as possible and embedded along the boundary of the monitoring region to capture waves in all directions from the crack.

Inverse Simulation and Imaging Conditions
It can be seen from Figure 6a-c that the inhomogeneities of the media can distort and scatter the wave field. However, this distortion cannot redirect the wave fronts and the scattering effect is not strong enough to overwhelm the wave path from the wave source (crack) to the piezoelectric sensors. The wave reversibility to the crack is valid and hence TRM is adaptive for concrete media.   For the healthy state shown in Figure 5a, the wave fronts travel horizontally and undisturbed, and the corresponding displacements at nodes of CPIMA are used as the baseline. In Figure 5b, we can see the obvious attenuation of the propagating wave caused by the crack, as well as the slight "leakage" of the propagating wave caused by water filling in the crack as a media. The comprehensive influence of the crack on the wave field can be seen in Figure 5c, which is a snapshot of the difference of the displacement field of Figure 5a,b. By differencing, the damage signal is extracted as if it is generated by the crack, including three components: (1) the attenuated wave, which travels right;

Inverse Simulation and Imaging Conditions
(2) the reflected wave, which travels left; and (3) the diffraction wave, which diffuses around from the crack tip. To capture all the signal components, the modules are suggested to be as more as possible and embedded along the boundary of the monitoring region to capture waves in all directions from the crack.
It can be seen from Figure 6a-c that the inhomogeneities of the media can distort and scatter the wave field. However, this distortion cannot redirect the wave fronts and the scattering effect is not strong enough to overwhelm the wave path from the wave source (crack) to the piezoelectric sensors. The wave reversibility to the crack is valid and hence TRM is adaptive for concrete media.

Inverse Simulation and Imaging Conditions
The inversion model simulates the focusing process of the backward wave field. Crack 1 was introduced as an example to the random aggregates FE model, whose location is shown in Table 2. Three piezoelectric actuators of the array are selected to re-emit the damage signal after the time reverse process. Based on the elastic wave theory, the energy density can be separated by two parts: P-wave (rotation free part) and S-wave (isochore part). Hence, three imaging conditions can be developed using P-wave energy, S-wave energy, and total energy density. For the good orientation, the P-wave energy-based imaging conditions are chosen in the inverse simulation. Figure 7a-f are snapshots of the typical focusing moments of the first P-wave fronts from 54.07 to 61.125 µs, from which we can see the P-wave fronts interfere constructively at different damage areas in different time steps. Separation of the S-wave and P-wave is also observed due to the difference in the velocity and direction. In some snapshots, the focus areas do not match the crack accurately; however, the accuracy can be improved by accumulating the energy values of displacement during the focusing period.

Experimental Setup
A total of three concrete specimens were studied, with the strength of the cement of 32.5 MPa and a mixture ratio of water, cement, fine sand, and stones of 0.44:1:1.42:3.17. The diameter of coarse aggregates was 3 mm~5 mm, and the size of the specimens was 500 mm × 150 mm × 150 mm. Four piezoelectric modules (PZT1~4) were fixed in each concrete beam before pouring. All piezoelectric modules were located on the neutral surface of the specimen (see Figure 8), with a maximum horizontal distance of 340 mm and a vertical distance of 50 mm.
After the curing period of 28 days, each concrete specimen was first tested in the healthy state to get reference signals by recording the received signals, which were propagated from the piezoelectric actuator and captured by the remaining three piezoelectric sensors. Note that the actuator was chosen in turn from PZT1 to PZT4, hence there were four groups of reference signals To display the energy field, the energy density of the P-wave at the position x is calculated as [40]: and normalized energy density is obtained: where u (x, T − t) is the displacement at location x at time step T − t. T is the total focusing time period. Ω is the detection domain. E(x) is the average value of the energy E(x) in the monitoring region. The different image conditions (energy calculation) have a great influence on the imaging results. In a previous study, there are two typical energy calculation methods based on the maximum energy value and the accumulated the energy by time [32,36]. The maximum energy is obtained by combining the maximum energy density of each node during the entire focusing period together, or recording the energy distribution when the overall maximum energy density appears. Considering the complex service environment and the background noise of underwater conditions, the accumulation energy method is applied to improve the robustness of the imaging method. For CPIMA with n piezoelectric modules, n tests are conducted for the detection of each crack. Each piezoelectric module is chosen as an actuator in turn from the first to the last, while the rest function as sensors. Hence, the average of n groups of energy density is used as the image conditionÊ avg (x).

Experimental Setup
A total of three concrete specimens were studied, with the strength of the cement of 32.5 MPa and a mixture ratio of water, cement, fine sand, and stones of 0.44:1:1.42:3.17. The diameter of coarse aggregates was 3 mm~5 mm, and the size of the specimens was 500 mm × 150 mm × 150 mm. Four piezoelectric modules (PZT1~4) were fixed in each concrete beam before pouring. All piezoelectric modules were located on the neutral surface of the specimen (see Figure 8), with a maximum horizontal distance of 340 mm and a vertical distance of 50 mm.

Experimental Setup
A total of three concrete specimens were studied, with the strength of the cement of 32.5 MPa and a mixture ratio of water, cement, fine sand, and stones of 0.44:1:1.42:3.17. The diameter of coarse aggregates was 3 mm~5 mm, and the size of the specimens was 500 mm × 150 mm × 150 mm. Four piezoelectric modules (PZT1~4) were fixed in each concrete beam before pouring. All piezoelectric modules were located on the neutral surface of the specimen (see Figure 8), with a maximum horizontal distance of 340 mm and a vertical distance of 50 mm.
After the curing period of 28 days, each concrete specimen was first tested in the healthy state to get reference signals by recording the received signals, which were propagated from the piezoelectric actuator and captured by the remaining three piezoelectric sensors. Note that the actuator was chosen in turn from PZT1 to PZT4, hence there were four groups of reference signals for each specimen. Then, each of the three specimens was cut into an artificial crack of crack 1, crack 2, and crack 3, respectively. The crack setup is shown in Table 2 and Figure 8. The damaged signals could be obtained by testing the cracked concrete specimen in the same manner. During the test, the specimens were immersed in water. The Agilent 33522A arbitrary waveform signal generator was employed to send signal to the piezoelectric actuators. A digital filter was applied to filter the low-frequency AC noise. Signals were recorded by an Agilent DSO7034B oscilloscope with a bandwidth of 350 MHz and a sampling frequency of 3 MHz. The experimental setup is shown in Figure 9. After the curing period of 28 days, each concrete specimen was first tested in the healthy state to get reference signals by recording the received signals, which were propagated from the piezoelectric actuator and captured by the remaining three piezoelectric sensors. Note that the actuator was chosen in turn from PZT1 to PZT4, hence there were four groups of reference signals for each specimen. Then, each of the three specimens was cut into an artificial crack of crack 1, crack 2, and crack 3, respectively. The crack setup is shown in Table 2 and Figure 8. The damaged signals could be obtained by testing the cracked concrete specimen in the same manner. During the test, the specimens were immersed in water.
The Agilent 33522A arbitrary waveform signal generator was employed to send signal to the piezoelectric actuators. A digital filter was applied to filter the low-frequency AC noise. Signals were recorded by an Agilent DSO7034B oscilloscope with a bandwidth of 350 MHz and a sampling frequency of 3 MHz. The experimental setup is shown in Figure 9. A 5-cycle 150 kHz Hanning-windowed tone burst, with a maximum amplitude of 10 V, was applied as the excitation signal. In order to reduce the influence of the reflected waves from the external boundary, the first arrival wave of the response signal was studied. The experimental data used were filtered by wavelet analyzer to extract signals of the dominant frequency of 150 kHz. The normalized excitation and response signal are shown in Figure 10a. The focus period was 0.24 ms, and the forward period and inverse period were both 0.12 ms, respectively. As shown in Figure 10b, the red line presents the received signal of the healthy beam, which was used as a reference. The green line is the received signal of the damaged beam. The blue line is the difference between the received signal with the crack and the reference signal, which presents the crack-related signal. The blue dash line is the time-reversed crack-related signal, which is used as the input signal in the inverse model.  A 5-cycle 150 kHz Hanning-windowed tone burst, with a maximum amplitude of 10 V, was applied as the excitation signal. In order to reduce the influence of the reflected waves from the external boundary, the first arrival wave of the response signal was studied. The experimental data used were filtered by wavelet analyzer to extract signals of the dominant frequency of 150 kHz. The normalized excitation and response signal are shown in Figure 10a. The focus period was 0.24 ms, and the forward period and inverse period were both 0.12 ms, respectively. As shown in Figure 10b, the red line presents the received signal of the healthy beam, which was used as a reference. The green line is the received signal of the damaged beam. The blue line is the difference between the received signal with the crack and the reference signal, which presents the crack-related signal. The blue dash line is the time-reversed crack-related signal, which is used as the input signal in the inverse model. A 5-cycle 150 kHz Hanning-windowed tone burst, with a maximum amplitude of 10 V, was applied as the excitation signal. In order to reduce the influence of the reflected waves from the external boundary, the first arrival wave of the response signal was studied. The experimental data used were filtered by wavelet analyzer to extract signals of the dominant frequency of 150 kHz. The normalized excitation and response signal are shown in Figure 10a. The focus period was 0.24 ms, and the forward period and inverse period were both 0.12 ms, respectively. As shown in Figure 10b, the red line presents the received signal of the healthy beam, which was used as a reference. The green line is the received signal of the damaged beam. The blue line is the difference between the received signal with the crack and the reference signal, which presents the crack-related signal. The blue dash line is the time-reversed crack-related signal, which is used as the input signal in the inverse model.
(a) (b) Figure 10. Extraction of the damage signal: (a) Excitation signal and response signal; (b) received signal of the beam with/without crack, their differential and reversed signal. Figure 10. Extraction of the damage signal: (a) Excitation signal and response signal; (b) received signal of the beam with/without crack, their differential and reversed signal.

Experimental Procedures
The experiment consists of three steps: (1) Forward process: Each of the three healthy beams were first placed underwater and tested to get reference signals. The depth of the water was about 0.15 m to precisely immerse the beam and avoid the effect of water pressure. After the test, the beam was carried out of the water to prevent water seepage in the concrete. Then, each beam was cut by a sewing machine into crack 1, crack 2, and crack 3, respectively. The cracked beams were immersed in water and tested again to get the damaged signals; thus, damage-related signals were extracted. The entire test process was finished within one week to avoid the effect of temperature change; (2) Inversion process: The first arrival wave of the time damage-related signals was reversed in the time domain and was imported as input signals in the inverse model. The obtained displacement and energy field were recorded; (3) imaging process: The energy value of each node was superimposed and converted to a true-color colorful cloud image after normalization to visualize the crack, see Figure 11.

Experimental Procedures
The experiment consists of three steps: (1) Forward process: Each of the three healthy beams were first placed underwater and tested to get reference signals. The depth of the water was about 0.15 m to precisely immerse the beam and avoid the effect of water pressure. After the test, the beam was carried out of the water to prevent water seepage in the concrete. Then, each beam was cut by a sewing machine into crack 1, crack 2, and crack 3, respectively. The cracked beams were immersed in water and tested again to get the damaged signals; thus, damage-related signals were extracted. The entire test process was finished within one week to avoid the effect of temperature change; (2) Inversion process: The first arrival wave of the time damage-related signals was reversed in the time domain and was imported as input signals in the inverse model. The obtained displacement and energy field were recorded; (3) imaging process: The energy value of each node was superimposed and converted to a true-color colorful cloud image after normalization to visualize the crack, see Figure 11. In order to verify the effectiveness of the proposed method and its robustness to noise, the above process was simulated by FE models as a comparison, and the imaging results are shown in Figure  12. In order to verify the effectiveness of the proposed method and its robustness to noise, the above process was simulated by FE models as a comparison, and the imaging results are shown in Figure 12

Results and Discussion
The imaging results include colors from dark blue to red consistent with the average energy densityÊ avg (x) from low to high, see Figure 11. Due to the reflection of the beam surface, the crack opening is usually intensively focused, represented by the area where the maximum energy density appears. A threshold level of 30% of the normalized maximum energy density was set to extract the damage areas. Considering that the crack shape is rectangle, the quantitative identification results of the crack can be measured by the crack opening position and crack depth, see Table 3.

Results and Discussion
The imaging results include colors from dark blue to red consistent with the average energy density  ( ) avg E x from low to high, see Figure 11. Due to the reflection of the beam surface, the crack opening is usually intensively focused, represented by the area where the maximum energy density appears. A threshold level of 30% of the normalized maximum energy density was set to extract the damage areas. Considering that the crack shape is rectangle, the quantitative identification results of the crack can be measured by the crack opening position and crack depth, see Table 3.    Figure 11 shows that underwater cracks, with different locations and depths, are imaged by the proposed TRM-based detection method, demonstrating that the extracted damage-related waves can carry the crack information and can be focused in the damage area through heterogeneous media (concrete). Additionally, the effect of the transmitted waves through water-filled-in cracks has been considered as a part of the attenuated wave component of the damage-related waves, hence the function of the crack as a wave source will not be influenced by the water. Except crack 3 in the experiment group, all the cracks in both the FE simulation group and experiment group are basically covered by the identified damage areas with shapes consistent with the direction of the actual cracks.
In the FE simulation group, note that the width of the identified damage area is wider near the crack opening with higher energy density, while the opposite phenomenon can be observed near the crack tip. This rule is more obvious for cracks with a shorter depth. For crack 2, crack 1, and crack 3, the width of the identified damage areas on the upper edge is 30, 25, and 18 mm, respectively. This is because for surface cracks with a shorter depth, a higher proportion of the crack-initiated waves is reflected by the surface of the concrete beam during propagation, making the surface a strong wave source comparable to the wave strength emitted from the crack. For crack 3, the identification of a deep surface crack is not seriously influenced by the reflected waves from the upper concrete surface; therefore, the energy density is equally distributed and the identified width is always about 18 mm along the vertical direction. However, an artifact appears about 20 mm below the crack tip. There are two possible reasons to explain this: The strong scattering of an exceptionally large concrete aggregate, or the information loss. From the random aggregates model (see Figure 4), no concrete aggregate with an exceptional large size is observed near the crack tip area, hence the possibility of strong scattering is excluded. Another reason is the information loss, which is caused by the layout and the number of piezoelectric modules employed. Due to the limited size and the shape of the beam specimen, all the four piezoelectric modules are embedded close to the middle line of the concrete beam, and their polar directions are parallel to the middle line of the beam. Therefore, only the information of waves propagating almost horizontally is available (see Figure 7). In addition, the useful information is limited by the number of sensor modules since the damage-related information is available at the sensors' location only. In the authors' opinion, this problem can be solved by placing the piezoelectric modules symmetrically on the boundary of the monitoring area and put them in multiple directions.
For the experimental imaging results, all three cracks are successfully detected, and the imaging results of crack 1 and crack 2 have good consistency with the actual ones. Crack 3 can also be identified from the image; however, the energy density of is unevenly distributed in the damage area due to the environmental effect, such as electromagnetic interference. Table 3 shows that the crack opening position in the horizonal direction is accurately identified in both the FE simulation and experiment groups. Except crack 3 in the experiment group (11 mm), all the errors of the crack opening position are less than 4 mm. Note that the wavelength of the P-wave used herein is about 26.7 mm (150 kHz frequency), and the errors of the crack opening position are much smaller than the wavelength. This indicates the possibility that the resolution of the wave-based method can be increased by TRM focus. In both the FE simulation and experiment groups, the relative errors of crack depth identification are under 12%, except the crack 3 in the FE simulation group (37.8%). However, the error of crack 3 is caused by the artifact. If the effect of artifact is eliminated, the relative error and error would be 3.8% and 19 mm, respectively. To compare the errors of vertical depth with that of the horizontal position, the errors of the vertical depth are nearly one order larger than that of the horizontal position, showing that the lateral resolution is higher than the longitudinal resolution. Considering the wavelength (26.7 mm), the error of the crack depth is acceptable.
The energy attenuation-based detection method and RMSD damage index are very established in concrete crack detection. For underwater conditions, all the RMSD values are above 0.2, indicating the sensitivity to the existence of concrete cracks. Further, for crack 2 and crack 3 in the same position, the RMSD value is changed consistently with the increase of the crack depth. This shows that the energy-based RMSD damage index is an effective and reliable approach to the underwater concrete crack detection and can be a reference to TRM. However, the TRM imaging result is able to provide much more information about the crack than RMSD.
In FE simulation, the fineness of the mesh influences the resolution of the model, not only for shaping the aggregates and cementing surface but also for simulating the propagating wave fronts. However, since the higher element density will consume additional computing resources and the self-focusing character of TRM does not depend on the wave paths, extra fine mesh is unnecessary when Neumann criterion is satisfied. The element size of 1 mm is conformance to the resolution requirements of the presented frequency (100, 150, and 200 kHz) and the scale of the applied specimen (meter scale).
Based on the wave-based detection theory, the dominant frequency is one of the determining factors of detection resolution. For the frequency of 150 kHz, the wavelength (26.7 mm, centimeter scale) is one order larger than the actual crack width (1 mm, millimeter scale). As a result, the width of the identified crack area is usually lager than the actual one. In theory, higher frequency means higher resolution but will lead to more serious attenuation (material damping) and scattering in concrete. Concrete can be regarded as a low-pass filter, and for frequencies over 400 kHz, which is believed in [41] to be the cut-off frequency in concrete, the elastic wave is extremely difficult to propagate due to the high damping of concrete. In many concrete damage detection experiments, the commonly used maximum size of concrete aggregate is less than 20 mm. When dominating frequencies are above 200 kHz, the wavelengths (below 20 mm) are close to the maximum size of concrete aggregates and this will lead to serious scattering. Therefore, the influence of frequency on the resolution is complex, and the increase of frequency is not equivalent to the improvement of accuracy in practice. It is necessary to evaluate the performance of the imaging method at different frequencies.
To analyze the influence of frequency on the accuracy of the proposed method, crack 1, crack 2, and crack 3 were identified by FE simulation with the dominating frequencies of 100 and 200 kHz. The imaging results are shown in Figure 12 and the quantitative crack imaging results of the crack opening position and crack depth are shown in Table 4.  Figure 12 shows that all the cracks in both the 100-kHz group and 200-kHz group can be detected and basically covered by the identified damage areas. In the 100-kHz group, the widths of the identified damage areas are significantly larger than that of the 200-kHz group, showing that the imaging resolution is directly related to the dominating frequency. In crack 3 of both groups, the artifacts are also observed near the crack tips. Table 4 shows that the crack opening position is accurately identified in both groups with errors (below 6 mm) substantially less than the wavelengths (20-40 mm). To compare the accuracy of the two groups, it can be observed that the accuracy of both groups is similar with respect to the horizontal opening position and vertical depth. Further, by comparing the quantitative imaging results in Table 3 and in the FE simulation group in Table 4, it can be concluded that similar accuracy is achieved at frequencies of 100, 150, and 200 kHz. This shows that the robustness to the frequency is in the presented frequency range, and an appropriate excitation frequency can be selected based on the specific working conditions, such as the scale of the specimen, maximum size of the concrete aggregate, and resonant frequency of PZT transducers.

Conclusions
In this study, a CPIMA monitoring system and TRM-based imaging technique were applied to detect and image underwater cracks of concrete structures. To overcome the challenge of water filled in the crack and to improve the application of TRM, the concept of damage-related waves was proposed. A heterogeneous FE numerical model, including the piezoelectric constitutive model and the wave propagation model characterized by displacement, was established. This model was used to simulate both the forward process and inverse process of TR focusing. The forward process depicts the wave propagation through healthy and damaged underwater concrete beams, showing the effectiveness of the proposed damage-related wave based TRM focus in underwater conditions and its robustness to heterogeneous material (concrete).
The inverse process simulates the backpropagating wave field, then the focusing results of different cracks were obtained. Three underwater concrete cracks, with different locations and depths, were imaged by the proposed TRM-based detection method by both FE simulation and experiment. The FE simulation results show the feasibility of the proposed method in heterogeneous media to detect underwater cracks. The reflection of the concrete surface near the crack opening has some influence on the shape of the identified areas but has little influence on the accuracy of the identification of crack opening. The influence of the dominant frequency on the resolution is complex, and for the presented different frequencies, the accuracy of the imaging results is similar. The experimental imaging results show that the underwater surface cracks in concrete beams are imaged by limited piezoelectric modules with an acceptable error, indicating the application potential in underwater crack detection of concrete structures. Various effects (water pressure, temperature, and seepage etc.) on the performance of piezoelectric modules and received signals need to be further studied before being adapted to in situ monitoring for underwater concrete structures.