Nonlinear Ultrasonic Imaging for Porosity Evaluation

The influence of porosity on the mechanical behaviour of composite laminates represents a complex problem that involves many variables. Therefore, the evaluation of the type and volume content of porosity in a composite specimen is important for quality control and for predicting material behaviour during service. A suitable way to evaluate the porosity content in composites is by using nonlinear ultrasonics because of their sensitivity to small cracks. The main objective of this research work is to present an imaging method for the porosity field in composites. Two nonlinear ultrasound techniques are proposed using backscattered signals acquired by a phased array system. The first method was based on the amplitude of the half-harmonic frequency components generated by microbubble reflections, while the second one involved the frequency derivative of the attenuation coefficient, which is proportional to the porosity content in the specimen. Two composite samples with induced porosity were considered in the experimental tests, and the results showed the high accuracy of both methods with respect to a classic C-scan baseline. The attenuation coefficient results showed high accuracy in defining bubble shapes in comparison with the half-harmonic technique when surface effects were neglected.


Introduction
In the last few decades, composite materials have been widely used in aeronautic and civil fields due to their high strength and lightweight characteristics. Despite these advantages, composites reinforced with continuous fibres can be affected by the presence of porosity, defined as the concentration of microscopic interfacial voids located in the matrix between plies of fibres and scattered through the thickness of the composite laminate [1]. The presence of porosity in composites may be due to several causes, such as the material constituents and the manufacturing process. In many cases, insufficient pressure during the curing is the most prevalent factor that affects the formation of porosity in a composite laminate [1]. It is well-known that porosity has detrimental effects on the strength of composite specimens, especially with regard to the material's inter-laminar shear strength [2][3][4][5], compressive strength [6], and bending strength [7], which are associated with matrix-dominated mechanical properties. It has also been shown that these voids result in a measurable decrease in fatigue life [7]. In general, the influence of voids on the mechanical behaviour of composite laminates represents a complex problem that involves a large number of variables, such as shape, size and location of voids; mechanical properties of fibre, matrix and interface; and mechanical loads present and their nature (static or cyclic) [8]. Moreover, it should take into account several environmental factors, such as temperature and moisture [4]. Therefore, the evaluation of the type and volume content of porosity in a composite specimen is essential in the prediction of material behaviour during the service of the specimen and in evaluating its useful life. Current research has focused on the reduction of porosity during composite manufacturing [9], with the assessment of porosity during service largely neglected. Thus, the focus of this work is the non-destructive evaluation of porosity during the service life of composite materials using an ultrasonic technique. Ultrasonic methods are currently the most widely used nondestructive techniques for the quantitative evaluation of porosity in composites because of their simplicity and high sensitivity to the presence of voids. Indeed, there are many previous theoretical and experimental works that explore the dependence of the ultrasonic wave speed and attenuation in composite laminates on the amount of porosity within that same specimen [10][11][12][13][14][15][16][17].
It should be noted that other non-destructive methods involving X-ray computed tomography [18] and thermography [19], or destructive methods based, for example, on acid digestion [20], are not considered in this paper.
Among all the ultrasound techniques, the so-called "linear" methods have shown a high level of accuracy for the detection and localisation of damage in composites; nevertheless, they lack sensitivity in the detection of micro-damage, such as porosity. In the latter case, the creation of nonlinear elastic effects, such as higher harmonics and subharmonics of the excitation frequency, can be used as a signature for micro-damage detection. Therefore, to overcome the limitation of the linear methods in porosity estimation, a nonlinear approach is followed in this paper.
Currently, there are limited methods available for the investigation of porosity in composite materials. These include traditional methods, such as optical microscopy, where cross-sectional samples are inspected under a microscope, and image-based algorithms are used to identify voids; this has limited application due to its destructive nature. X-ray and Neutron Imaging provide a non-destructive alternative where the inspected volumes can be reconstructed and analysed to identify porosity [21,22]. Finally, laser ultrasonic techniques can be used to induce wave propagation within the material, providing information regarding subsurface porosity [23][24][25]. These methods have some limitations regarding equipment cost, inspection speed, and size.
Early work done by Ashwell et al. [26] and Lauterborn [27] showed that nonlinear oscillations of gas bubbles forced into volume pulsations by an acoustic field in liquids might occur at frequencies equal to half of the driving frequency, the subharmonic. Furthermore, Leighton et al. [23] showed that it is possible to accurately size and resolve the bubbles in liquids by interrogating the subharmonic emitted by the bubble, which can be directly related to the resonance of these bubbles. It has also been shown that nonlinear responses occur at the subharmonic as well as at sum and difference frequencies when two frequencies are used to excite bubbles [28,29]. Research into this area of nonlinear ultrasonics is driven by the increased sensitivity to damage and damage evolution, predicted at a few orders of magnitude higher compared with linear techniques, which are only sensitive to gross defects rather than micro-damage [30][31][32]. Nonlinear ultrasound methods can be broken into a few general methods, which include single-frequency excitation, modulation (multiple-frequency mixing), and harmonic generation (including subharmonics) [33][34][35][36][37][38][39]. Among them, nonlinear elastic wave spectroscopy (NEWS) [40][41][42], nonlinear elastic wave modulation spectroscopy (NEWMS) [43][44][45], nonlinear imaging methods [46][47][48], and phase symmetry analysis (PSA) techniques [49,50] were used to explicitly interrogate the nonlinear behaviour of the medium and its effect on wave propagation. Many studies have shown that nonlinearities are generated in structures in the presence of ruptures, voids, cohesive bonds, and the opening/closing of surfaces or micro-cracks [51][52][53]. Material nonlinearity of damaged regions, in terms of nonlinear harmonic generation (due to ultrasound excitation), has been widely used over the past few years as a reliable signature for material flaw detection in both metallic and composite structures [46,54,55].
The aim of this research work was to evaluate induced porosity in two undamaged composite specimens using two different nonlinear ultrasound approaches. The first evaluates and correlates the production of subharmonic responses to porosity, while the second takes into account the attenuation coefficient of the medium. It is expected that gas trapped within resin resonates at specific frequencies to provide subharmonic content in the captured signal, which can then be used to analyse porosity. The first step in this approach is to provide evidence that subharmonic generation can be achieved with porosity, with future work focused on sizing. The two techniques showed different advantages and disadvantages; however, both were able to visualise levels of porosity within the tested composite samples. It should be noted that this work looks to incorporate multiple techniques to increase the probability of the detection of porosity rather than relying on a single method.
The outline of this research work is as follows: the first part of Section 2 describes the nonlinear ultrasound approach and presents both the theory and applications of a standard phased array system. The subharmonic and attenuation coefficient techniques are illustrated, respectively, in Sections 2.1 and 2.2. Section 3 shows the setup used to perform the experimental tests, while in Section 4, all the results are illustrated. The conclusions of this paper are reported in Section 5.

Nonlinear Ultrasound Approach
As aforementioned in Section 1, the best approach to extract porosity and void information in a composite specimen is by using nonlinear ultrasound techniques. Therefore, in this paper, two nonlinear ultrasound approaches were developed and investigated to produce a porosity map of the tested composite samples. These maps were compiled by the application of the developed theory on the backscattered signals sent and acquired by a phased array system.
Nonlinear ultrasound phased array techniques have been assessed by many authors. Ohara et al. [56] and Park et al. [57] have extensively developed, evaluated, and improved the detection of open and closed cracks in metallic structures using a subharmonic phased array, while Potter et al. [58] have developed a nonlinear array based on the traditional TFM. These methods provide nonlinear ultrasound information through the thickness of the sample, rely on short excitation signals, have been assessed on metallic structures, and do not generally address issues of equipment-based nonlinearities and porosity.
Ultrasonic phased array techniques arguably lead the field in terms of damage detection capabilities and suitability for composite material structures. Phased array systems generally use an array of transmitting and receiving piezoelectric elements, where the sequence of firing and capturing determines the method and accuracy. These systems generally use three methods for damage assessment and imaging, such as the plane/focus swept method, full matrix capture (FMC), and total focusing method (TFM), a post-process technique [59,60]. The geometry of a general phased array is highlighted in Figure 1, while Figure 2 presents the scanning technique used. These B-scan techniques can then be performed along the surface of the structure (y-axis) in order to generate a C-scan image.   The TFM method uses a post-processing algorithm which discretises the thickness of a sample into a grid. The signals from all the elements in the array (referring to that grid point) are then summed to focus on every point in the grid. The intensity of the image, , at any point in the scan, is given by: where: is the width of the aperture; and are distances in the z-direction (normal to array, depth) and x-direction (position along array) for transmitting and receiving elements (see Figure 1); ℎ , is the Hilbert transform of the time domain signal; The TFM method uses a post-processing algorithm which discretises the thickness of a sample into a grid. The signals from all the elements in the array (referring to that grid point) are then summed to focus on every point in the grid. The intensity of the image, I(x, z) at any point in the scan, is given by: where: D is the width of the aperture; z and x are distances in the z-direction (normal to array, depth) and x-direction (position along array) for transmitting (tx) and receiving elements (rx) (see Figure 1); h tx,rx is the Hilbert transform of the time domain signal; c 1 is the speed of sound in the medium; s is the step along the array, which sets the resolution of the B-scan.
If an ultrasound array is considered and set up as in Figure 2a,b below, multiple k transmit and receive positions (pulse-echo testing) will result in multiple time-domain response signals for each point tested (A-scans). Multiple A-scans along a line (x-direction) can be used to generate a B-scan with all points (x and y-direction), generating a C-scan. The acquired data were post-processed by the authors using MATLAB software code.
Therefore, the described phased array system can acquire the back-scattered signals caused by several material variations, such as porosity, interfacial resin-rich areas, variations in fibre density, waviness in the plies, and others [61].

Subharmonic Method
The resonance of a bubble in a liquid is a well-studied phenomenon, with early works focusing on ultrasonic cavitation of bubbles within liquid mediums, where the amplitude and pressure of an ultrasonic pressure wave has been shown to either grow or collapse these bubbles under pressure. The oscillation of such bubbles in a liquid medium can be described using highly non-linear equations, such as the Rayleigh-Plesset equation [62,63]: where R is the bubble radius, R 0 its equilibrium value, σ surface tension, ρ density, µ viscosity of the liquid, P 0 hydrostatic pressure, P(t) time-varying pressure component, P v vapour pressure within the bubble, and κ the polytropic index of the gas.
To approximate the resonant oscillations of an unforced bubble, viscosity and surface tension are assumed negligible, P(t) is set to zero [29], and the frequency is deemed [64]: ( Although these formulations refer to bubbles in liquid mediums, they give insight into the resonant behaviour of entrapped gases in mediums, but also the porosity within the resin. Thus, it is possible to exploit the behaviour and ability of ultrasonic-excited porosity to generate subharmonic responses, as the principle behind the mechanism remains the same. Furthermore, it is well-known that both micro-cracks and delamination, when excited by ultrasonic waves, can generate nonlinear material responses that can be analytically modelled using the classical nonlinear elasticity (CNE) theory [65]. Assuming a onedimensional longitudinal wave propagation along the x-direction, the elastodynamic wave equation, as shown by Equation (1), can be expressed as the power series of the strain as follows: where σ is the stress, λ and µ are, respectively, Lamé's first and second parameter, β and γ are the second and third order elastic coefficients, defined as the ratios between the second and third harmonic amplitudes and the square fundamental amplitude [39]. A normalised version of these parameters, which are only functions of the fundamental and high-order harmonic amplitudes, can be used as damage indicators because of their highest values at the damage location. The bispectral analysis is a valid alternative to the nonlinear coefficient's methods, as it can be used to measure both the magnitude and phase of the higher-order harmonic frequency (second and third-order components) [42]. Nevertheless, out of the presented nonlinear methods, many research works demonstrated the benefit of subharmonic analysis when microbubbles, therefore porosity, are considered [53,[65][66][67][68][69][70]. A generic amplitude spectrum of an acquired filtered backscattered signal can be presented as shown in Figure 3. Figure 3 shows the presence of both higher-order harmonics (at twice A2-2f1 and trice A3-3f1) and a half harmonic (f1/2), with the latter indicating the presence of microbubbles, such as porosity. The name is justified by the fact that this component is at half the fundamental frequency, which represents the frequency of the transmitted tone-burst signal (see Figure 4). Moreover, it is well-known that if a damaged sample is excited with a single frequency, the interaction between the generated waves and, for example, delamination produces contact acoustic nonlinearity in the form of higher harmonics, which provides a good measure for identifying the existence of delamination and determining their sizes in composites. Therefore, the spectrum presented in Figure 3 is due to the presence of both porosity and delamination. burst signal (see Figure 4). Moreover, it is well-known that if a damaged sample is excited with a single frequency, the interaction between the generated waves and, for example, delamination produces contact acoustic nonlinearity in the form of higher harmonics, which provides a good measure for identifying the existence of delamination and determining their sizes in composites. Therefore, the spectrum presented in Figure 3 is due to the presence of both porosity and delamination.   In the experimental tests of this research work (see Section 4), only the contribution of the fundamental and the half-harmonic frequency components were considered since only microbubbles were present. The reliability of half-harmonic amplitudes of the received backscattered signals was noticeable (the A-scans, see an example in Figure 5) for  In the experimental tests of this research work (see Section 4), only the contribution of the fundamental and the half-harmonic frequency components were considered since only microbubbles were present. The reliability of half-harmonic amplitudes of the received backscattered signals was noticeable (the A-scans, see an example in Figure 5) for the visualisation of the porosity field in the composite samples, identified by a 2D map, with each pixel representative of the extracted amplitude value of the half-harmonic component of a single backscattered signal.

Attenuation Coefficient Method
In Section 1, the relationship between the wave attenuation and the microscopic spheroidal dispersed voids within a composite sample was emphasised [15]. Therefore, a technique based on the attenuation coefficient was investigated and presented in this research work. This technique is based on a pulse-echo method achieved by a phased array system, as performed in the previous algorithm (see Section 2.1), with a 1-cycle tone burst as the transmitted signal (see Figure 4). If a single A-scan is considered and acquired, the amplitude spectrum of the ultrasonic pulse is calculated. It should be noted that the spectrum can be written as: where is the amplitude spectrum of the reference ultrasonic signal, is the thickness of the specimen, and represents the attenuation coefficient function depending on the frequency that can be determined from Equation (5) as follows: Once the attenuation coefficient function is obtained, it is possible to evaluate the porosity content in a composite sample since the proportionality between porosity and frequency are derivative of the attenuation coefficient function. Figure 6 shows an example of the attenuation coefficient calculated using Equation (6). The red line indicates the linear trend of the function, where the slope gradient is representative of the porosity details along a line from a single point of the specimen surface

Attenuation Coefficient Method
In Section 1, the relationship between the wave attenuation and the microscopic spheroidal dispersed voids within a composite sample was emphasised [15]. Therefore, a technique based on the attenuation coefficient was investigated and presented in this research work. This technique is based on a pulse-echo method achieved by a phased array system, as performed in the previous algorithm (see Section 2.1), with a 1-cycle tone burst as the transmitted signal (see Figure 4). If a single A-scan is considered and acquired, the amplitude spectrum S( f ) of the ultrasonic pulse is calculated. It should be noted that the spectrum can be written as: where S 0 ( f ) is the amplitude spectrum of the reference ultrasonic signal, H is the thickness of the specimen, and α( f ) represents the attenuation coefficient function depending on the frequency that can be determined from Equation (5) as follows: Once the attenuation coefficient function is obtained, it is possible to evaluate the porosity content in a composite sample since the proportionality between porosity and frequency are derivative of the attenuation coefficient function. Figure 6 shows an example of the attenuation coefficient calculated using Equation (6). The red line indicates the linear trend of the function, where the slope gradient is representative of the porosity details along a line from a single point of the specimen surface through the thickness. The application of the described algorithm to all the acquired A-scans leads to a 2D map showing the inner porosity field. Currently, the attenuation coefficient is generally used among other factors, such as frequency parameters, as a rapid tool to evaluate the structural quality of the composites [66].

Experimental Setup
In this study, two squares of composite samples, with dimensions of 200 mm × 200 mm × 9 mm and induced porosity due to manufacturing, were considered for the experimental tests (see Figure 7).

Experimental Setup
In this study, two squares of composite samples, with dimensions of 200 mm × 200 mm × 9 mm and induced porosity due to manufacturing, were considered for the experimental tests (see Figure 7). The phased array system (Diagnostic Sonar-FI Toolbox) described in Section 2 (see Figure 8) utilised a 5 MHz 128-element probe (Diagnostic Sonar) with an element pitch of 0.5257 mm and width of 0.5 mm. A stepped linear (or sequential) C-scan routine was used to evaluate the porosity in the two samples. The step linear routine used beams of 32 elements, equating to a total scan width of~58.8 mm. A half-cycle sine-wave tone burst with an amplitude of 80 V was used as the transmitting input signal. Due to the span of the array (~58.8 mm), four horizontal scanning stripes were necessary to cover the whole area of the two samples (see Figure 8). Raw data over the C-scan length were then post-processed according to the algorithms outlined in Sections 2.1 and 2.2.

Experimental Setup
In this study, two squares of composite samples, with dimensions of 200 mm × 200 mm × 9 mm and induced porosity due to manufacturing, were considered for the experimental tests (see Figure 7). The phased array system (Diagnostic Sonar-FI Toolbox) described in Section 2 (see Figure 8) utilised a 5 MHz 128-element probe (Diagnostic Sonar) with an element pitch of 0.5257 mm and width of 0.5 mm. A stepped linear (or sequential) C-scan routine was used to evaluate the porosity in the two samples. The step linear routine used beams of 32 elements, equating to a total scan width of ~58.8 mm. A half-cycle sine-wave tone burst with an amplitude of 80 V was used as the transmitting input signal. Due to the span of the

Results
This section illustrates the obtained results by using the algorithms described in Sections 2.1 and 2.2. The results from the two developed methods are compared with classic C-scan results, which are regarded as the baseline and called "C-map". The mentioned outputs consist of 2D maps of the three specimens, visualising the detected inner porosity field. The phased array probe was used to capture A-scans over the inspected area, and the subharmonic amplitude of each A-scan was used to build a 2D map, here named "Smap". Likewise, another 2D map was created from the derivative, in the frequency domain, of the attenuation coefficient at each scan point. This map output is named "Dmap".
In Figure 9, for each specimen, S-maps and D-maps are compared to the respective C-maps. The percentage P of the porous area with respect to the total area was calculated for each specimen as follows (see Table 1): where Abubble is the detected 2D area of a single bubble, whilst Atot is the total area of the specimen.

Results
This section illustrates the obtained results by using the algorithms described in Sections 2.1 and 2.2. The results from the two developed methods are compared with classic C-scan results, which are regarded as the baseline and called "C-map". The mentioned outputs consist of 2D maps of the three specimens, visualising the detected inner porosity field. The phased array probe was used to capture A-scans over the inspected area, and the subharmonic amplitude of each A-scan was used to build a 2D map, here named "S-map". Likewise, another 2D map was created from the derivative, in the frequency domain, of the attenuation coefficient at each scan point. This map output is named "D-map".
In Figure 9, for each specimen, S-maps and D-maps are compared to the respective C-maps. The percentage P of the porous area with respect to the total area was calculated for each specimen as follows (see Table 1): where A bubble is the detected 2D area of a single bubble, whilst A tot is the total area of the specimen.   According to Table 1, the S-maps showed more defects due to larger surface areas for identified porosity regions as well as the identification of smaller porosity, with an average of 405% higher detection than the C-Map. However, the D-maps provided a more uniform identification of each void, with an average of 194% higher detection than the C-Map. It should be noted that in each of the samples highlighted in Figure 9, there are areas of porosity that are clearly visible using one method but not using the other method. This reinforces the assumption that a combined approach, using multiple techniques with different advantages and disadvantages, increases the probability of detection.  According to Table 1, the S-maps showed more defects due to larger surface areas for identified porosity regions as well as the identification of smaller porosity, with an average of 405% higher detection than the C-Map. However, the D-maps provided a more uniform identification of each void, with an average of 194% higher detection than the C-Map. It should be noted that in each of the samples highlighted in Figure 9, there are areas of porosity that are clearly visible using one method but not using the other method. This reinforces the assumption that a combined approach, using multiple techniques with different advantages and disadvantages, increases the probability of detection.

Conclusions
This paper presented two nonlinear ultrasound techniques able to furnish a detailed 2D map of the porosity field within two undamaged samples. A pulse-echo method was performed by using a phased array system, able to transmit a tone-burst signal and receive the backscattered signal, representative of the sample content through the thickness. Data coming from all the sample points were post-processed, visualising the inner porosity by a reconstructed 2D map. A half-harmonic technique was implemented due to the microbubble peculiarity to emit a frequency equal to half the transmitted one [67][68][69][70]. A second method was based on the attenuation coefficient due to the relation between the wave attenuation and the presence of microscopic voids in a composite sample. Indeed, the frequency derivative of the attenuation coefficient is proportional to the porosity content in the specimen. The obtained results were compared, with a classic C-scan regarded as a baseline. All the obtained maps showed different levels and locations of inner porosity, with the S-Map and D-Map showing, on average, 405% and 194% more porosity in terms of area compared with the standard C-Map. While these results would suggest that the C-Map should be disregarded, this is not the case, as each of the methods provided useful and differing evaluations of the levels of porosity in the tested samples in terms of the location and size of the affected areas. Thus, it is suggested that the methods are used in tandem in order to maximise the probability of detection.
In order to maximise the probability of detection, further research is needed to combine the outputs from the various methods by developing algorithms to provide a single image of the levels of porosity.