Analysis of Acoustic Emission Activity during Progressive Failure in Heterogeneous Materials: Experimental and Numerical Investigation

: This work focuses on an experimental and numerical investigation into monitoring damage in a cube-shaped concrete specimen under compression. Experimental monitoring uses acoustic emission (AE) signals acquired by two independent measurement apparatuses, and the same damage process is numerically simulated with the lattice discrete element method (LDEM). The results from the experiment and simulation are then compared in terms of their failure load, ﬁnal conﬁgurations, and the evolution of global parameters based on AE signals, such as the b -value coefﬁcient and the natural time approach. It is concluded that the results from the AE analysis present a signiﬁcant sensitivity to the characteristics of the acquisition systems. However, natural time methods are more robust for determining such differences, indicating the same general tendency for all three data sets.


Introduction
Acoustic emission (AE) methods are a promising approach to monitoring the damage process in structural systems, and thus have found numerous applications in various engineering branches.Briefly, AE methods refer to the occurrence of micro-fractures within a structure inducing elastic waves, which propagate throughout the system and can be detected with suitable sensors, usually located on external surfaces [1].In the present work, the waves emitted by the sources are called events, and the responses captured by the sensors are called signals.
In this paper, the damage process in a cube of concrete is studied through AE methods.The cube is subjected to compression until it reaches failure, with global results presented in terms of the failure load and final configuration of the specimen.Two independent systems perform the corresponding AE data acquisition: eight sensors connected to Lunitek [2] and two AE sensors connected to Mistras [3] conditioning devices.The damage process is also simulated using a discrete element method (DEM) version.Information from this simulation and both acquisition systems support the investigation of damage evolution through different AE global parameters, such as the b-value, initially used by Ritcher [4] and applied in several articles, most notably by Carpinteri et al. [5], and a natural time (NT) analysis was proposed by Varotsos et al. [6].
In the following, we present the essential aspects of the primary analysis tools used here: the b-value and NT techniques for analyzing AE data and the DEM version used in the simulations, originally proposed by [7].
b-value analysis: this damage assessment criterion is provided by the statistical analysis of the amplitude distribution of AE signals generated by growing microcracks, as shown schematically in Figure 1.This distribution of amplitudes (A) conforms to the Gutenberg-Richter (GR) law [4], N(≥ A) ∝ A −b , where N is the number of AE signals with amplitude ≥ A. The GR law's exponent b, the so-called b-value, changes with the different stages of damage growth: while the initially dominant microcracking generates large numbers of low-amplitude AE signals, the subsequent macrocracking generates fewer signals of higher amplitude.As a result, the b-value progressively reduces as the specimen approaches impending failure.
Appl.Sci.2022, 12, x FOR PEER REVIEW 2 of 20 In the following, we present the essential aspects of the primary analysis tools used here: the b-value and NT techniques for analyzing AE data and the DEM version used in the simulations, originally proposed by [7].
b-value analysis: this damage assessment criterion is provided by the statistical analysis of the amplitude distribution of AE signals generated by growing microcracks, as shown schematically in Figure 1.This distribution of amplitudes ( A ) conforms to the Gutenberg-Richter (GR) law [4], ( ) , where N is the number of AE signals with amplitude A ≥ .The GR law's exponent b, the so-called b-value, changes with the different stages of damage growth: while the initially dominant microcracking generates large numbers of low-amplitude AE signals, the subsequent macrocracking generates fewer signals of higher amplitude.As a result, the b-value progressively reduces as the specimen approaches impending failure.Pradhan et al. [8] offers an alternative interpretation of the b-value in connection with the so-called bundle model, extensively explored by [9], which consists of a set of parallel-connected bars with a random strength distribution that are progressively loaded until model failure is reached.The histogram of burst distributions (i.e., the number of simultaneously broken bars) and their magnitudes have a characteristic slope in the bi-logarithmic dominium, which matches the b-value in AE signals from physical systems.Thus, b-values calculated from a simulated model can be used as a precursor to critical stages in similar physical systems.
Natural time (NT) analysis: This approach, proposed by Varotsos et al. [25], is proven to be ideal for enhancing the characteristics of the studied signals, reducing uncertainties, and optimizing the extraction of information, especially in systems approaching criticality.Relevant examples of this method's usefulness appear in a large variety of fields.For specific applications in the analysis of AE signals, see [26][27][28][29].
For a time series comprising N events, the k -th event's natural time k χ is de- fined as that event's normalization with respect to the total number of considered events, i.e., k k N χ = , serving as an index for its occurrence as shown in Figure 2. From its definition, it is clear that . In addition to k χ , we also define k Q as a desired measure to characterize the signals.Its specific physical meaning is arbitrary, but energy and amplitude are common choices.Its corresponding normalized value is given by serves as a basis for NT analysis.Pradhan et al. [8] offers an alternative interpretation of the b-value in connection with the so-called bundle model, extensively explored by [9], which consists of a set of parallel-connected bars with a random strength distribution that are progressively loaded until model failure is reached.The histogram of burst distributions (i.e., the number of simultaneously broken bars) and their magnitudes have a characteristic slope in the bilogarithmic dominium, which matches the b-value in AE signals from physical systems.Thus, b-values calculated from a simulated model can be used as a precursor to critical stages in similar physical systems.
Natural time (NT) analysis: This approach, proposed by Varotsos et al. [25], is proven to be ideal for enhancing the characteristics of the studied signals, reducing uncertainties, and optimizing the extraction of information, especially in systems approaching criticality.Relevant examples of this method's usefulness appear in a large variety of fields.For specific applications in the analysis of AE signals, see [26][27][28][29].
For a time series comprising N events, the k-th event's natural time χ k is defined as that event's normalization with respect to the total number of considered events, i.e., χ k = k/N, serving as an index for its occurrence as shown in Figure 2. From its definition, it is clear that χ k ∈ [0, 1].In addition to χ k , we also define Q k as a desired measure to characterize the signals.Its specific physical meaning is arbitrary, but energy and amplitude are common choices.Its corresponding normalized value is given by The pair (χ k , p k ) serves as a basis for NT analysis.
A given natural time χ k 's variance κ 1 is calculated as: κ 1 varies when a new k-th event occurs, once natural time χ k changes from k/N to k/(N + 1) and p k changes to Q k /∑ N+1 i=1 Q i .As shown in [25], such changes may identify a dynamical system approaching a critical point.κ 1 is often called the system's order parameter.A given natural time k χ 's variance 1 κ is calculated as: 1 κ varies when a new k -th event occurs, once natural time k χ changes from k N to ( ) As shown in [25], such changes may iden- tify a dynamical system approaching a critical point. 1 κ is often called the system's order parameter.
Additionally, the entropy in natural time S can be defined as: Notice that, as observed by [30], this parameter is only based on numerical signal data and does not indicate a thermodynamic relation.Just like the variance 1 κ , the en- tropy S is also a dynamic parameter depending on the sequential order of events [31].
The third parameter in NT analysis is the entropy under time reversal rev S .It is computed by reversing the order of the energy data ( k Q ), recalculating S in Equation (2) with this reversed series, and subtracting S 's original value.Finally, Varotsos et al. [6], propose yet another order parameter, describing the "average" distance D between the normalized power spectrum ( ) ( ) Additionally, the entropy in natural time S can be defined as: Notice that, as observed by [30], this parameter is only based on numerical signal data and does not indicate a thermodynamic relation.Just like the variance κ 1 , the entropy S is also a dynamic parameter depending on the sequential order of events [31].
The third parameter in NT analysis is the entropy under time reversal S rev .It is computed by reversing the order of the energy data (Q k ), recalculating S in Equation (2) with this reversed series, and subtracting S's original value.
Both S and S rev become smaller than the entropy of uniform noise S u = 0.0966.
Lattice discrete element method (LDEM): As LDEM can simulate crack propagation and evolution, it captures internal changes in the model that resemble the emission of AE signals during the damage process.This method was used by Riera [7] to solve structural dynamics problems in an orthotropic continuum, with stiffness coefficients derived for a linearly elastic solid by Nayfeh and Hefzy [37], and by using the explicit direct numerical integration of the equations of motion.The LDEM was also employed to represent the behavior of quasi-brittle specimens in [10,38], showing a good correlation with numerical methods and experimental results, not only in terms of final configuration and global parameters, but also in the simulation of acoustic emission.
The method consists of representing a solid by an arrangement of masses with three degrees of freedom, interconnected by massless elements carrying only axial loads.Each nodal mass's three degrees of freedom correspond to the nodal displacements in the coordinate directions.For the adopted cubic array presented in Figure 3a The uniaxial constitutive law adopted for the LDEM bars is based on Hilleborg's [41] proposal for quasi-brittle materials, which adopts a triangular-element constitutive relationship, Figure 3b.Note that failure in compression occurs by indirect tension.The element axial force F depends on the axial strain ε .An equivalent fracture area * i A of the i -th element is defined to ensure that the energies dissipated by fracture both in the continuum and in its discrete representation are equivalent, considering a cubic sample of dimensions L L L × × [10].A basic parameter of the bilinear constitutive relationship shown in Figure 3b is given by: where p ε denotes the strain at the peak of the bilinear law, f G is the fracture energy and eq d is a material characteristic length, also considered a material property, discussed in detail by Taylor [42].The i -th element's ultimate strain r ε , at which an element subjected to tension has no residual carrying capacity, is defined by: It is observed that r ε depends both on the material properties and the discretization level.More details about the LDEM formulation employed here and how to calibrate it for quasi-brittle structures appear in previous contributions by the authors, such as [43].
The non-homogeneous nature of engineering materials such as concrete is accounted for in the LDEM using two approaches.The first approach defines the material's The initial axial stiffness (Figure 3b) of normal and diagonal elements are EA n /L n and EA d /L d respectively, where E is the Young's modulus of the material.The cubic array adopted herein results in the exact representation of an isotropic continuum for ν = 0.25.On the other hand, shear differences arise for ν = 0.25, as detailed in [39].Such differences are very small and can be neglected if 0.20 ≤ ν ≤ 0.30.
A system of equations of motion is obtained by applying Newton's second law to each node.The material's internal damping is assumed to depend on nodal masses and corresponding velocities.The equations of motion are thus uncoupled and may be integrated in the time domain using an explicit finite difference scheme.The convergence of the solution using LDEM for both linear elastic problems and elastic stability problems was verified by different authors, as reported in [40].
The uniaxial constitutive law adopted for the LDEM bars is based on Hilleborg's [41] proposal for quasi-brittle materials, which adopts a triangular-element constitutive relationship, Figure 3b.Note that failure in compression occurs by indirect tension.The element axial force F depends on the axial strain ε.An equivalent fracture area A * i of the i-th element is defined to ensure that the energies dissipated by fracture both in the continuum and in its discrete representation are equivalent, considering a cubic sample of dimensions L × L × L [10].A basic parameter of the bilinear constitutive relationship shown in Figure 3b is given by: where ε p denotes the strain at the peak of the bilinear law, G f is the fracture energy and d eq is a material characteristic length, also considered a material property, discussed in detail by Taylor [42].The i-th element's ultimate strain ε r , at which an element subjected to tension has no residual carrying capacity, is defined by: It is observed that ε r depends both on the material properties and the discretization level.More details about the LDEM formulation employed here and how to calibrate it for quasi-brittle structures appear in previous contributions by the authors, such as [43].
The non-homogeneous nature of engineering materials such as concrete is accounted for in the LDEM using two approaches.The first approach defines the material's fracture energy G f as a 3D random field with a Weibull probability distribution, which also defines ε p for each i-th element used in the simulations (see Equation ( 5)) according to its location (considering the barycenter of the bar) [44].The correlation lengths L cx , L cy and L cz along the corresponding axes x, y and z determine the fracture energy's linear spatial correlation.The second is to introduce small mesh perturbations in the regular cubic LDEM arrangement through small initial dislocations of the nodal points.This approach leads to significant improvements in predicting the responses from structures under compressive loading when instability effects are not negligible [45,46].By numerical experimentation, considering a perfectly cubic nodal array as the initial state, the coefficient of variation, CV p , for the normally distributed perturbations, was chosen as 2.5%.
One should note that the focus of this study is not on the experimental procedure per se, but on the sensitivity of each analysis method on the data sets from the different sources, especially regarding their ability to identify the system's entry into a critical damage state.

Experimental Analysis 2.1. Test Description
Acoustic emission data were obtained from a standard axial compression test of a concrete sample, using a static hydraulic test machine [47] with a capacity of 500 kN.The concrete specimen, presented in Figure 4, was positioned directly between the press's steel plates, without Teflon layers.The compressive force steadily increased at a constant load rate of 1.5 kN/s until reaching the peak load.After that, the load was first reduced by 80% and then totally lifted.
state, the coefficient of variation, p CV , for the normally distributed perturbations, was chosen as 2.5%.One should note that the focus of this study is not on the experimental procedure per se, but on the sensitivity of each analysis method on the data sets from the different sources, especially regarding their ability to identify the system's entry into a critical damage state.

Test Description
Acoustic emission data were obtained from a standard axial compression test of a concrete sample, using a static hydraulic test machine [47] with a capacity of 500 kN.The concrete specimen, presented in Figure 4, was positioned directly between the press's steel plates, without Teflon layers.The compressive force steadily increased at a constant load rate of 1.5 kN/s until reaching the peak load.After that, the load was first reduced by 80% and then totally lifted.The acquisition of AE signals was made through two independent multichannel systems in parallel: "AEmission" by Lunitek [2] with 8 piezoelectric sensors (model LT18-003-PRD-00-R0 with a frequency range of 15-625 kHz) and "MicroSHM" by Mistras [3] with 2 piezoelectric sensors (model R15A-150 KHZ general-purpose AE sensor with a frequency range of 20-600 kHz).The piezoelectric transducers were fixed to the free faces of the cube with epoxy resin (cold glue X60) and aluminum adhesive tape.Both acquisition systems use a 5 MHz sampling frequency and store data in parametric form.The recorded parameters for each detected signal are:

•
Signal start time = instant of the first reading that exceeds the detection threshold: 40 dB (corresponding to 100 µV) for the Mistras [48], it can be approximated by the envelope area of the triangle formed by the signal peak amplitude and its duration.

Results and Discussions
Figures 5 and 6 show the AE results from the compression test using the Lunitek and Mistras systems, respectively.All plots used the normalized time t/t * , where t * = 3000 s is the time when the load reaches 80% of its peak value.The AE data presented after t/t * = 1 were caused by the acquisition systems continuing to run for some time after the compression test ended, but such data were not considered in the following analyses.
Figures 5a and 6a show the load against the normalized time.Notably, once again, the amplitude threshold changes with the acquisition system in each case.The maximum load P max = 4117 kN occurs at t/t * = 0.927, corresponding to a compressive stress of 45.75 MPa.
Figure 5b,c and Figure 6b,c illustrate the signal distribution in terms of average frequency (AF) in (kHz) and signal duration in (ms).The frequency distribution is similar for both acquisition systems, but the Lunitek system captured longer-duration events after t/t * =0.8, which did not appear in the Mistras system.
The instantaneous and accumulated number of signals are presented in Figures 5d and 6d.The accumulated number of counts appears with the same general tendency in both plots.For the interval roughly between t/t * = 0.2 and t/t * = 0.6, when the load reaches 70% of the peak load, the number of counts grows with an approximately uniform rate.After that, the count growth becomes pronouncedly nonlinear, which indicates an unstable evolution in the damage process.The expressive differences in total counts for each acquisition system come from the different threshold values.That is also the likely reason for the differences in accumulated AE energy, signal energy distribution, and signal duration time.
The evolution of the b-value during the test appears in Figures 5e and 6e.The test period was divided into 14 equal intervals, and the b-value was calculated for each of them.The same tendency occurs for both data sets: after presenting relatively low variations for about 60% of the total time, the coefficient suddenly decreases to within the range [1, 1.5], indicating imminent instability in the structure.
The behavior indicated through the b-value is coherent with the natural time (NT) results, Figures 5f and 6f, which also present the following auxiliary horizontal lines: κ 1 = 0.070 (dashed red), < D > limit = 10 −2 (dashed green), and entropy limit S u ≈ 0.0966 (dashed blue).The corresponding analysis was performed according to the same methods given in [49], taking the energy component (Q k ) as equal to A max 1.5 , where A max is the maximum amplitude of the signal (in µV), an approach also commonly found in seismology studies [50].The calculation of the term A max 1.5 was performed hit-by-hit, i.e., every time a new signal was identified, it was included in the NT analysis, leading to the rescaling of the (χ k ,Q k ) time series and the recalculation of κ 1 , S, S rev and < D >.For this analysis, imminent instability was indicated when all four parameters reached their indicated critical values.In Figures 5f and 6f, this occurs for the region highlighted with a gray background, which is within the interval where the b-value lies in the [1, 1.5] range.However, the unstable region identified with data from the Mistras system is broader than that identified from the Lunitek system.Note that for t/t * near 0.1, a narrow critical region was identified with the Mistras system but not with the Lunitek system.
Figures 5g and 6g present the information condensed by the parameter < D >, i.e., the comparison between expressions (3) and ( 4), representing the initial parts of the theoretical spectrum and the real one.The comparison is carried out for three different regions indicated in Figures 5f and 6f with black circles.It is clear that the two spectra are close when the region analyzed is inside the gray region, that is, when the critical state is reached.
Although the b-value and NT analyses yield results mostly in agreement for both data sets, some significant discrepancies are noticeable in terms of (i) differences in total counts for each acquisition system, the (ii) accumulated AE energy, and (iii) signal duration time.As we discuss now, such discrepancies are most likely attributable to differences in the following factors for each acquisition system: threshold levels, sensors characteristics, and the sensibility of signal-conditioning circuits.The influence of the different threshold levels can be somewhat directly addressed here.As for the other factors, a deeper understanding can be reached only through further analyses in future work.
Figure 7 represents the duration of the signals (in µs) recorded by the two acquisition systems versus their corresponding number of oscillations (counts or threshold crossings), which calculates the signals' frequencies.Although the nominal operation ranges for the two systems are very similar (20 kHz-600 kHz for Mistras, 15 kHz-625 kHz for Lunitek), the frequency ranges of captured signals are very different: between 5 and 94 kHz for Mistras, and between 5 and 45 kHz for Lunitek.Furthermore, the accumulated counts with the first system are much larger than with the second one.
Finally, both systems registered similar counts for large-and medium-amplitude signals (80-110 dB and 60-80 dB, respectively), whereas significant differences appeared for signals between 40 and 60 dB, which is precisely the range where the threshold differences occur: 40 dB for Mistras; 49 dB for Lunitek.Thus, although the b-value and NT parameters tend to be the same for both data sets, their calculation is affected by the number of available samples.
Figure 8 presents a comparison between the Lunitek (grey background) and Mistras (purple background) data sets in terms of their statistical distribution of peak amplitudes (db), average frequency (AF) (kHz), and signal duration (µs).The histograms' information is complemented in Table 1, containing their corresponding average values ( .), standard deviations (σ) and coefficients of variation (CV).Once again, the threshold difference clearly affects the shapes of all three distributions, especially the peak amplitude.It should also be noted that the higher signal duration in the case of the Lunitek results appears in the tail of the histogram presented in Figure 8c.
The AF distribution for both acquisition systems is similar with a sharper peak for the Lunitek system and a flatter profile for the Mistras system.The different threshold levels are certainly an important cause for this effect, as illustrated in Figure 7.With a lower threshold, the Mistras system registers a broader spectrum of frequencies because smaller cracks tend generate AE signals with a lower amplitude and higher frequencies (or equivalently, lower duration times).
As preliminary conjectures, one can also suggest other causes for the observed differences in the frequency distribution.For instance, the sensors from each system differ in mass and geometric dimensions, which affects their resonance characteristics.More specifically, since the Mistras sensor has a larger mass (105 g compared to 52 g for Lunitek), it is also expected to present a higher damping ratio, which makes it difficult to record different events that are close to each other because they would appear to be a single event in the sensor's response.Such challenges will be investigated in future work.s is the time when the load reaches 80% of its peak value.The AE data presented after * / t t = 1 were caused by the acquisition systems continuing to run for some time after the compression test ended, but such data were not considered in the following analyses.Finally, both systems registered similar counts for large-and medium-amplitude signals (80-110 dB and 60-80 dB, respectively), whereas significant differences appeared for signals between 40 and 60 dB, which is precisely the range where the threshold differences occur: 40 dB for Mistras; 49 dB for Lunitek.Thus, although the b-value and NT parameters tend to be the same for both data sets, their calculation is affected by the number of available samples.
Figure 8 presents a comparison between the Lunitek (grey background) and Mistras (purple background) data sets in terms of their statistical distribution of peak amplitudes (db), average frequency ( AF ) (kHz), and signal duration (μs).The histograms' information is complemented in Table 1, containing their corresponding average values ( .), standard deviations (σ ) and coefficients of variation ( CV ).Once again, the threshold difference clearly affects the shapes of all three distributions, especially the peak ampli- tude.It should also be noted that the higher signal duration in the case of the Lunitek results appears in the tail of the histogram presented in Figure 8c.The AF distribution for both acquisition systems is similar with a sharper peak for the Lunitek system and a flatter profile for the Mistras system.The different threshold levels are certainly an important cause for this effect, as illustrated in Figure 7.With a lower threshold, the Mistras system registers a broader spectrum of frequencies because smaller cracks tend generate AE signals with a lower amplitude and higher frequencies (or equivalently, lower duration times).

Amplitude (db)
As preliminary conjectures, one can also suggest other causes for the observed differences in the frequency distribution.For instance, the sensors from each system differ  We now investigate the effect of the sensor location, as illustrated in Figure 9, which depicts the distribution of AE sources along the z coordinate of the sample (see Figure 4) according to the localization methodology proposed in [1].The AE activity is preponderant along the positive semi-axis in both plots.When an elastic wave hits a piezoelectric AE sensor coupled to the surface of a material under damage, the surface displacement at the nanometer scale is converted into an electric signal.The wave mode and frequency content of the AE signal arriving at the sensor are affected by the geometry and the wave propagation properties of the material between source and sensor.An AE signal is detected if the arrival of the AE wave causes an electric signal in the AE sensor whose amplitude exceeds the selected detection threshold.
Transient AE signals are analyzed and reduced to a set of parametric features which are recorded by the data acquisition system.Essential AE signal features that are extracted and stored are arrival time at the sensor, maximum peak amplitude, signal duration, and energy.Most AE signal features depend on the selected measurement parameters, e.g., the detection threshold.Therefore, it is not possible to speak in absolute terms of the accuracy of experimental results in tests with AE.Too many parameters influence the obtained outcomes: the signal attenuation properties of the analyzed material, environmental noise, sensor sensitivity, and the predetermined detection threshold, in particular.In this way, when comparisons are made between the energy of AE signals detected by the piezoelectric sensors and those actually emitted by the materials under stress, they can only be made with arbitrary units or by normalizing the current values to the maximum values obtained.
Only by taking into account the problem of the AE source localization in various materials can an efficient study on the accuracy of results from AE tests be conducted.
In the literature, there are various studies on this problem; see for example [51][52][53].From these studies, it is generally clear that the accurate positioning of an AE event is realized by adjusting the number of sensors and improving the picking accuracy of the first arrival of the signal.As a matter of fact, the highest accuracy can be achieved for damage sources localized among the sensor distributions.Distances of sources far from the sensor network are less accurate.When an elastic wave hits a piezoelectric AE sensor coupled to the surface of a material under damage, the surface displacement at the nanometer scale is converted into an electric signal.The wave mode and frequency content of the AE signal arriving at the sensor are affected by the geometry and the wave propagation properties of the material between source and sensor.An AE signal is detected if the arrival of the AE wave causes an electric signal in the AE sensor whose amplitude exceeds the selected detection threshold.
Transient AE signals are analyzed and reduced to a set of parametric features which are recorded by the data acquisition system.Essential AE signal features that are extracted and stored are arrival time at the sensor, maximum peak amplitude, signal duration, and energy.Most AE signal features depend on the selected measurement parameters, e.g., the detection threshold.Therefore, it is not possible to speak in absolute terms of the accuracy of experimental results in tests with AE.Too many parameters influence the obtained outcomes: the signal attenuation properties of the analyzed material, environmental noise, sensor sensitivity, and the predetermined detection threshold, in particular.In this way, when comparisons are made between the energy of AE signals detected by the piezoelectric sensors and those actually emitted by the materials under stress, they can only be made with arbitrary units or by normalizing the current values to the maximum values obtained.
Only by taking into account the problem of the AE source localization in various materials can an efficient study on the accuracy of results from AE tests be conducted.
In the literature, there are various studies on this problem; see for example [51][52][53].From these studies, it is generally clear that the accurate positioning of an AE event is realized by adjusting the number of sensors and improving the picking accuracy of the first arrival of the signal.As a matter of fact, the highest accuracy can be achieved for damage sources localized among the sensor distributions.Distances of sources far from the sensor network are less accurate.
Onset time determination of an AE signal is another important factor in this analysis.A widely used approach, based on the Akaike information criterion (AIC), is confirmed to provide more reliable onset time determination of AE signals, see [53].Moreover, the inhomogeneity of the test object may lead to errors in the ultrasonic wave travelling model.Different-scale fluctuations and structure variations in the composite structure can result in random variations of the propagation velocity and systematic errors.An example of the difference in the accuracy of the experimental results on AEs is provided in Figure 9 which depicts the distribution of AE sources' localization along the z direction of the cube sample.Most importantly, the AE activity is preponderant along the positive semi-axis in both plots obtained by Lunitek and Mistras systems.
Therefore, the main goal of the present article, rather than focusing on the Lunitek and Mistras systems' accuracy, is to compare the two systems' outcomes.This comparison is important because the Lunitek and Mistras systems present different configurations and, consequently, different outputs, but in both cases, the global parameter evolution indicates the same trend, which is the greatest advantage of the AE technique.
Figure 10 presents a comparison between the signals captured over the two opposite faces of the cube where the sensors for both acquisition systems were placed together: at z = −15 cm (Lunitek L1 and L2, Mistras M1), and z = 15 cm (Lunitek L5 and L6, Mistras M2), as depicted in Figure 4b.A significant difference in signals counts for the Mistras sensors indicates that the damage developed more intensely close to the sensor M2.Such a difference is not as clear for the Lunitek system.This discrepancy is likely due to the location of the sensors on the cube's faces: the Mistras sensors are near the center of each face, whereas the Lunitek sensors are close to the specimen vertices.
Onset time determination of an AE signal is another important factor in this analysis.A widely used approach, based on the Akaike information criterion (AIC), is confirmed to provide more reliable onset time determination of AE signals, see [53].Moreover, the inhomogeneity of the test object may lead to errors in the ultrasonic wave travelling model.Different-scale fluctuations and structure variations in the composite structure can result in random variations of the propagation velocity and systematic errors.An example of the difference in the accuracy of the experimental results on AEs is provided in Figure 9 which depicts the distribution of AE sources' localization along the z direction of the cube sample.Most importantly, the AE activity is preponderant along the positive semi-axis in both plots obtained by Lunitek and Mistras systems.
Therefore, the main goal of the present article, rather than focusing on the Lunitek and Mistras systems' accuracy, is to compare the two systems' outcomes.This comparison is important because the Lunitek and Mistras systems present different configurations and, consequently, different outputs, but in both cases, the global parameter evolution indicates the same trend, which is the greatest advantage of the AE technique.
Figure 10 presents a comparison between the signals captured over the two opposite faces of the cube where the sensors for both acquisition systems were placed together: at z = −15 cm (Lunitek L1 and L2, Mistras M1), and z = 15 cm (Lunitek L5 and L6, Mis- tras M2), as depicted in Figure 4b.A significant difference in signals counts for the Mistras sensors indicates that the damage developed more intensely close to the sensor M2.Such a difference is not as clear for the Lunitek system.This discrepancy is likely due to the location of the sensors on the cube's faces: the Mistras sensors are near the center of each face, whereas the Lunitek sensors are close to the specimen vertices.Aside from using the experimental data sets from both acquisition systems, the damage analysis methods were also applied to the simulation results, as described in Section 3.

Model Description and Calibration
The simulated LDEM model consists of a cube with 300 mm edges, composed of cubic modules with edge length L = 4 mm, resulting in 75 modules along each axis.The model's vertical displacement is fixed at the solid's lower face (the horizontal plane at y = −150 mm), which corresponds to the area in contact with the hydraulic press's lower Aside from using the experimental data sets from both acquisition systems, the damage analysis methods were also applied to the simulation results, as described in Section 3.

Model Description and Calibration
The simulated LDEM model consists of a cube with 300 mm edges, composed of cubic modules with edge length L = 4 mm, resulting in 75 modules along each axis.The model's vertical displacement is fixed at the solid's lower face (the horizontal plane at y = −150 mm), which corresponds to the area in contact with the hydraulic press's lower steel plate in the experiments, while its upper face (the horizontal plane at y = +150 mm) suffers a prescribed vertical displacement to simulate the experimental load.This arrangement is illustrated in Figure 11a.
steel plate in the experiments, while its upper face (the horizontal plane at y = +150 mm) suffers a prescribed vertical displacement to simulate the experimental load.This arrangement is illustrated in Figure 11a.MPa, which could be reduced by about 50% due to the interaction between the structure's boundary conditions and the random distribution of the fracture energy.On the other hand, the compressive strength reached during the experiment was 45.75 MPa.In general, quasi-brittle materials tend to present failure global tensile stresses close to 10 times lower than the global compressive stress.Then, for this example, the tensile stress is about 5 MPa, which is compatible with the upper estimated here.
Regarding the 3D random field for determining f G , the adopted correlation length L for all three directions.The coefficient of variation, Gf CV , was taken as 50%.For further information about the methodology followed in the calibration, see [43].

Simulation Results and Discussion
Figure 11a presents both numerical and experimental results in terms of failure configurations.The numerical crack patterns are compatible with those from the experiment, Figure 11a, especially along the vertical direction.
Further comparisons between numerical and experimental data are shown in Figure 12a in terms of applied load versus time, normalized to the instant when the load returns to 80% of the peak load after that peak is reached.This return period is shorter for the LDEM because it is impossible to directly prescribe the displacement in the experiment, as the hydraulic-press-manipulated variable is the pressure it applies to the specimen.However, the load peak is very similar for both cases, differing by less than 10%.The material had the following mechanical properties: mass density ρ =2400 kg/m 3 , Young's modulus E = 40 GPa, and fracture energy G f = 200 N/m.Considering d eq = 0.045 m and using expression (5), the constitutive law's critical strain ε p (Figure 3b) is: Thus, the expected upper bound for the global tensile failure stress is Eε p = 13.3MPa, which could be reduced by about 50% due to the interaction between the structure's boundary conditions and the random distribution of the fracture energy.On the other hand, the compressive strength reached during the experiment was 45.75 MPa.In general, quasi-brittle materials tend to present failure global tensile stresses close to 10 times lower than the global compressive stress.Then, for this example, the tensile stress is about 5 MPa, which is compatible with the upper estimated here.
Regarding the 3D random field for determining G f , the adopted correlation length (L c ) was L/2 for all three directions.The coefficient of variation, CV G f , was taken as 50%.For further information about the methodology followed in the calibration, see [43].

Simulation Results and Discussion
Figure 11a presents both numerical and experimental results in terms of failure configurations.The numerical crack patterns are compatible with those from the experiment, Figure 11a, especially along the vertical direction.
Further comparisons between numerical and experimental data are shown in Figure 12a in terms of applied load versus time, normalized to the instant when the load returns to 80% of the peak load after that peak is reached.This return period is shorter for the LDEM because it is impossible to directly prescribe the displacement in the experiment, as the hydraulic-press-manipulated variable is the pressure it applies to the specimen.However, the load peak is very similar for both cases, differing by less than 10%.
The energy balance computed during the simulation appears in Figure 12b, where the kinetic energy is multiplied by 50 to facilitate visualization.Since this energy component is at significantly lower levels than those for the other two, the quasi-static condition is reached in the simulation.
The energy balance computed during the simulation appears in Figure 12b, where the kinetic energy is multiplied by 50 to facilitate visualization.Since this energy component is at significantly lower levels than those for the other two, the quasi-static condition is reached in the simulation.The acoustic emission activity during the simulation was calculated through two methods: (i) computing the increment kinetic energy [10], k E Δ , and (ii) using a virtual sensor [38], shown at the position S1 in Figure 11a.The increment kinetic energy is defined by Iturrioz et al. [10] between two successive integration times as: where k E is kinetic energy in the discretized time during the simulation.Since k E during a given event is proportional to that event's AE energy, ( ) to the time derivative of ( ) is also proportional to the signal's peak distribution.The acquisition of the AE waves in the simulation was performed by a virtual sensor placed on a node on the simulated cube's surface.Such a sensor registers the node's acceleration perpendicular to the surface where it stands, thus emulating the effect of the piezoelectric accelerometers used in the AE experiment.
Figure 13 shows a comparison of the acoustic emission activity using the two approaches.The vertical axis is normalized with respect to the highest signal captured before the specimen reached the peak load.Notice that the k E Δ distribution in Figure 13a is measured directly in the event source, whereas the acceleration wave must travel through the body to reach the virtual AE sensor.The identification of the same event through both methods is shown in detail in Figure 13c.The acoustic emission activity during the simulation was calculated through two methods: (i) computing the increment kinetic energy [10], ∆E k , and (ii) using a virtual sensor [38], shown at the position S1 in Figure 11a.The increment kinetic energy is defined by Iturrioz et al. [10] between two successive integration times as: where E k is kinetic energy in the discretized time during the simulation.Since E k during a given event is proportional to that event's AE energy, ∆E k (t i ) is proportional to the time derivative of E k (t i ), implying that ∆E k (t i ) is also proportional to the signal's peak distribution.The acquisition of the AE waves in the simulation was performed by a virtual sensor placed on a node on the simulated cube's surface.Such a sensor registers the node's acceleration perpendicular to the surface where it stands, thus emulating the effect of the piezoelectric accelerometers used in the AE experiment.
Figure 13 shows a comparison of the acoustic emission activity using the two approaches.The vertical axis is normalized with respect to the highest signal captured before the specimen reached the peak load.Notice that the ∆E k distribution in Figure 13a is measured directly in the event source, whereas the acceleration wave must travel through the body to reach the virtual AE sensor.The identification of the same event through both methods is shown in detail in Figure 13c.
As an AE parameter precursor of the failure of the numerically analyzed structure, only the natural time was used.Using the computed ∆E k (Figure 13a) as the measure of AE activity, the calculation of the natural time parameters followed the same procedure used with the experimental data, that is, using hit-by-hit and considering Q k to be equal to the event's AE energy.
Figure 14 shows the evolution of the NT parameters during the simulation.The region with critical behavior appears between normalized instants, 0.81 and 0.92, where the four parameters κ 1 , S, S rev and < D > reach the characteristics limits indicated in Section 1.As an AE parameter precursor of the failure of the numerically analyzed structure, only the natural time was used.Using the computed k E Δ (Figure 13a) as the measure of AE activity, the calculation of the natural time parameters followed the same procedure used with the experimental data, that is, using hit-by-hit and considering k Q to be equal to the event's AE energy.

Conclusions
This work analyzed acoustic emission (AE) data from a damage test performed on a cubic concrete specimen.The damage development within the structure was studied in terms of two AE-based analysis methods: b-value and natural time.These techniques were applied to three different data sets: two from the experimental test, registered independently with different acquisition systems, and one from simulating the same test using a version of the lattice discrete element method (LDEM).It was possible to conclude that:

•
The differences in data sets due to each acquisition apparatus have significant effects on the conclusions from AE analysis methods, such as the b-value.The main reason for this likely lies in the amplitude-frequency relationship in the emission of acous-

Conclusions
This work analyzed acoustic emission (AE) data from a damage test performed on a cubic concrete specimen.The damage development within the structure was studied in terms of two AE-based analysis methods: b-value and natural time.These techniques were applied to three different data sets: two from the experimental test, registered independently with different acquisition systems, and one from simulating the same test using a version of the lattice discrete element method (LDEM).It was possible to conclude that: and Technological Development (CNPq), the Coordination for the Improvement of Higher Level of Education Personnel (CAPES) and the sponsorship guaranteed with basic research funds provided by Politecnico di Torino, Italy for their financial aids in this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.b-value parameter extraction from an event set.

Figure 1 .
Figure 1.b-value parameter extraction from an event set.

Figure 2 .
Figure 2. The NT analysis method: (a) a typical AE time series; (b) corresponding NT representation.

Figure 2 .
Figure 2. The NT analysis method: (a) a typical AE time series; (b) corresponding NT representation.

20 Figure 3 .
Figure 3. LDEM: (a) Basic cubic module employed in the discretization and (b) bilinear constitutive law attributed to the elements.

Figure 3 .
Figure 3. LDEM: (a) Basic cubic module employed in the discretization and (b) bilinear constitutive law attributed to the elements.

Figure 4 .
Figure 4. (a) Compression test scheme with sample dimensions 0.3 × 0.3 × 0.3 m and, (b) AE sensors' locations.The coordinate system's origin is at the cube's geometric center.The acquisition of AE signals was made through two independent multichannel systems in parallel: "AEmission" by Lunitek[2] with 8 piezoelectric sensors (model LT18-003-PRD-00-R0 with a frequency range of 15-625 kHz) and "MicroSHM" by Mistras[3] with 2 piezoelectric sensors (model R15A-150 KHZ general-purpose AE sensor with a frequency range of 20-600 kHz).The piezoelectric transducers were fixed to the free faces of the cube with epoxy resin (cold glue X60) and aluminum adhesive tape.Both acquisition systems use a 5 MHz sampling frequency and store data in parametric form.

Figure 4 .
Figure 4. (a) Compression test scheme with sample dimensions 0.3 × 0.3 × 0.3 m and, (b) AE sensors' locations.The coordinate system's origin is at the cube's geometric center.

Figures 5
Figures 5 and 6 show the AE results from the compression test using the Lunitek and Mistras systems, respectively.All plots used the normalized time

Figure 5 .
Figure 5. b-value and NT analysis, data from the Lunitek system: (a) compressive-load and AE signal amplitudes, (b) average frequency (AF), (c) signal duration, (d) instantaneous and accumulated signal counts and normalized AE energy (E AE ), (e) b-value evolution and (f) natural time analysis.The horizontal axis is the same for all figures.(g) Evolution of Π(φ) with AE data from the Lunitek system for the true coincidence instants (black circles in (f)) in normalized time equal to (from left to right): 0.7492, 0.7877, and 0.8506.

Figure 6 .
Figure 6.b-value and NT analysis, data from the Mistras system: (a) compressive-load and AE signal amplitudes, (b) average frequency ( AF ), (c) signal duration, (d) instantaneous and accu- mulated signal counts and normalized AE energy ( AE E ), (e) b-value evolution, and (f) natural time analysis.The horizontal axis is the same for all figures.(g) Evolution of ( ) φ Π

Figure 6 .
Figure 6.b-value and NT analysis, data from the Mistras system: (a) compressive-load and AE signal amplitudes, (b) average frequency (AF), (c) signal duration, (d) instantaneous and accumulated signal counts and normalized AE energy (E AE ), (e) b-value evolution, and (f) natural time analysis.The horizontal axis is the same for all figures.(g) Evolution of Π(φ) with AE data from the Mistras system for the true coincidence instants (black circles in (f)) in normalized time equal to (from left to right): 0.698, 0.7601, 0.8621.

Figure 7
Figure 7 represents the duration of the signals (in μs) recorded by the two acquisition systems versus their corresponding number of oscillations (counts or threshold crossings), which calculates the signals' frequencies.Although the nominal operation ranges for the two systems are very similar (20 kHz-600 kHz for Mistras, 15 kHz-625 kHz for Lunitek), the frequency ranges of captured signals are very different: between 5 and 94 kHz for Mistras, and between 5 and 45 kHz for Lunitek.Furthermore, the accumulated counts with the first system are much larger than with the second one.

Figure 7 .
Figure 7. Relation between counts versus duration for: the (a) Lunitek and, (b) Mistras system.Amplitude in db is present in the colorbar.

Figure 7 .
Figure 7. Relation between counts versus duration for: the (a) Lunitek and, (b) Mistras system.Amplitude in db is present in the colorbar.

Figure 11 .
Figure 11.(a) Model discretization, boundary conditions, and simulated final rupture configuration.The red color represents fractured bars, and blue indicates partially damaged ones.Undamaged bars are omitted for enhanced visualization.(b) Experimental final rupture configuration.

Figure 11 .
Figure 11.(a) Model discretization, boundary conditions, and simulated final rupture configuration.The red color represents fractured bars, and blue indicates partially damaged ones.Undamaged bars are omitted for enhanced visualization.(b) Experimental final rupture configuration.

Figure 12 .
Figure 12.(a) Comparison between experimental and numerical results in terms of load versus normalized time.(b) Energy balance presented during the simulation.

Figure 12 .
Figure 12.(a) Comparison between experimental and numerical results in terms of load versus normalized time.(b) Energy balance presented during the simulation.

Figure 13 .
Figure 13.AE activity registered during the simulations: (a) using the increment of kinetic energy, k E Δ , (b) with a virtual acoustic emission sensor.(c) Details of two typical signals for the ap- proaches (a,b).
Figure 14 shows the evolution of the NT parameters during the simulation.The region with critical behavior appears between normalized instants, 0.81 and 0.92, where the four parameters 1 κ , S , rev S and D < > reach the characteristics limits indicated in Section 1.

Figure 13 .
Figure 13.AE activity registered during the simulations: (a) using the increment of kinetic energy, ∆E k , (b) with a virtual acoustic emission sensor.(c) Details of two typical signals for the approaches (a,b).

Figure 14 .
Figure 14.Natural time analysis versus the normalized time for AE activity registered in the LDEM.

Figure 14 .
Figure 14.Natural time analysis versus the normalized time for AE activity registered in the LDEM.
system, and 49 dB (280 µV) for Lunitek; • Peak amplitude, expressed in dB (A dB = 20log10(V max /1 µV)); • Number of oscillations (counts) = the number of intermediate crossings of the threshold by the signal.This measure is also commonly used to estimate the AE signal's average frequency (AF) through the counts/duration relationship (signal duration ≡ end time-start time); • AE energy (E AE ) = integral of the waveform.As indicated by RILEM

Table 1 .
Coefficient of variation ( CV ), average values ( . ) and standard deviations (σ ) of the AE signal parameters.

Table 1 .
Coefficient of variation (CV), average values ( . ) and standard deviations (σ ) of the AE signal parameters.