Identification of brittle damage in a UHPC beam using ultrasonic transmission tomography

Ultrasonic Transmission Tomography (UTT) utilises the information contained in the ultrasonic signal after it has propagated through the medium under investigation. Measurements are performed along paths (so-called rays) between transmitting and receiving points, recording either the signal travel time or its attenuation coefficient at a given frequency, or, alternatively, the frequency derivative of the attenuation coefficient (Polakowski & Sikora, 2016). On this basis, after appropriate mathematical reconstruction, a tomographic image of the internal structure of the medium is obtained and indirectly represented by maps visualising the spatial distributions of ultrasonic wave velocity or attenuation coefficient. When measuring the travel time of an ultrasonic pulse, it should be noted that excitation of mechanical waves in the tested element may generate several wave types. In such cases, the analysis focuses on those wave modes that propagate fastest (e.g., bulk longitudinal waves in massive elements) or that dominate the recorded response (e.g., Lamb waves in plates (Perkowski et al., 2018)). Their average velocity can be estimated from the measured propagation time between transmitter and receiver, under the simplifying assumption that the fastest paths of sound propagation are straight.

To determine the spatial distribution of this quantity, an appropriate system of equations must be formulated (see, for example, Kak & Slaney, 1999). The reconstructed image is assumed to consists of a finite number of pixels. A computational grid composed of N rectangular cells of dimensions δ₁ × δ₂ (Fig. 3), is superimposed on the examined element, with each cell corresponding to a single pixel.

Fig. 3.

Scheme of a cell system, transmitting/receiving points, and rays in a plane area examined tomographically (own research)

The unknown function f(x, y) in the investigated domain is then approximated in a piecewise constant form by assuming that, within each cell, it takes a constant value f_(k)(k = 1,2, … , N).

Subsequently, if reconstruction of the material structure is based on measurements of the ultrasonic time-of-flight between transmitting points i = 1,2, … , M and receiving points j = 1,2, … , M, the measured travel time is related to the pulse propagation velocity and to the geometry of the path (ray) along which it propagates by the following expression: (1) Δtij=∫lilj1cldl, \Delta {t_{ij}} = \int_{{l_i}}^{{l_j}} {{1 \over {c\left( l \right)}}dl} , where: c

– longitudinal wave velocity as a function of position along the ray [m/s];

Consequently, in the present formulation the unknown function f is defined as the reciprocal of the wave velocity c. Moreover, recognising that only averaged values of f within the individual cells are to be determined, we can write: (2) Δtij=∑k=1Nlijkfk, \Delta {t_{ij}} = \sum\limits_{k = 1}^N {l_{ij}^{\left( k \right)}{f_{\left( k \right)}}}, where lijk l_{ij}^{\left( k \right)} – length of the segment of the ray l_ij = |l_il_j| contained within the k-th cell (if the ray does not intersect the cell k, then lijk=0)m l_{ij}^{\left( k \right)} = 0)\left[ m \right] .

The above expression defines a system of equations with respect to the unknowns f_(k), which is used to determine the distribution of ultrasonic wave propagation velocity in tomography. Upon solving this system, one obtains: (3) ck=1fk, {c_{\left( k \right)}} = {1 \over {{f_{\left( k \right)}}}}, where c_(k) – the average velocity (in the sense of the present method) of the identified wave mode in the k-th cell [m/s].

One limitation of the adopted formulation is that, for ultrasonic waves that may diffract when passing through regions with different acoustic impedances, the actual ray geometry is not known. Therefore, it is assumed, as a simplification, that the rays l_il_j are straight-line segments connecting points i and j. This assumption is one of the primary sources of error in the method, particularly when the investigated element contains subdomains whose acoustic impedance differs significantly from that of the surrounding material.

The number of rays l_il_j employed must be at least equal to the assumed number M. For notational convenience, system (2) is written in the following matrix form: (4) Ax=b, {\bf{Ax = b}}, where: A=[Aij]P×N {\bf{A}} = {[{A_{ij}}]_{P \times N}}

– coefficient matrix [m];

x=f1,f2,…,fN1×NT {\bf{x}} = \left[ {{f_{\left( 1 \right)}},{f_{\left( 2 \right)}}, \ldots ,{f_{\left( N \right)}}} \right]_{1 \times N}^T

– vector of unknowns, whose components are equal to the reciprocals of the ultrasonic wave propagation velocity in the individual cells [s/m];

b=Δt11,Δt12,…1×PT {\bf{b}} = \left[ {\Delta {t_{11}},\;\Delta {t_{12}}, \ldots } \right]_{1 \times P}^T

– vector of measured wave travel times along the individual rays [s];

In tomographic inverse problems, the linear system (4) is typically severely ill-conditioned, because the system matrix A is non-banded and exhibits a large condition number, i.e. small perturbations in the data vector b may produce large relative errors in the solution vector x. Consequently, the tomographic reconstruction problem should be solved using iterative algorithms, often in combination with regularisation techniques.

A fundamental tomographic image reconstruction method is the iterative projection algorithm introduced by the Polish mathematician Stefan Kaczmarz (1937). In 1970 it was rediscovered and popularised in medical imaging by Gordon and co-workers under the name Algebraic Reconstruction Technique (ART) (Gordon et al., 1970). This method was employed in the first X-ray computed tomography (CT) scanner developed by Hounsfield (Hounsfield, 1973).

On the basis of Kaczmarz’s algorithm, numerous variants have been developed. In the classical literature, three main families of algebraic methods are distinguished: ART, the Simultaneous Iterative Reconstruction Technique (SIRT) and the Simultaneous Algebraic Reconstruction Technique (SART), the latter combining advantageous features of both ART and SIRT (Hansen & Saxild-Hansen, 2012; Kak & Slaney, 1999).

Moreover, one of the methods used to solve the linear systems arising in tomography is the least-squares method with Tikhonov regularisation (Phillips, 1962). In Perkowski et al. (2018), reconstructions of the tested elements obtained using ART and Tikhonov-regularised least squares were presented. It was shown that the use of Tikhonov-regularised least squares yields physically realistic velocity maps, which directly translates into a more reliable assessment of the spatial distribution of material stiffness. For this reason, the present study adopts the least-squares method with Tikhonov regularisation.

The classical Least-Squares (LS) method provides an approximate solution to overdetermined systems, i.e. when the number of equations exceeds the number of unknowns. Its most common application is linear regression, which determines the parameters of the straight line that best fits (in the least-squares sense) a given set of points in the plane. In general terms, solving the system Ax = b in the LS sense consists in finding the minimum of an objective function defined as the sum of squared differences between the components of the vectors Ax and b: (5) Fx=Ax−b2=min, F\left( x \right) = {\left\| {{\bf{Ax}} - {\bf{b}}} \right\|^2} = min , which leads to the solution of the form: (6) x=ATA−1ATb. {\bf{x}} = {\left( {{{\bf{A}}^{\rm{T}}}{\bf{A}}} \right)^{ - 1}}{{\bf{A}}^T}{\bf{b}}.

If the system Ax = b is ill-conditioned, the stability of its solution can be improved by introducing regularisation. One of the most widely used approaches is Tikhonov regularisation, in which a so-called regularisation term is added to the objective function (5), and the minimum of the resulting functional is sought for an appropriately chosen regularisation parameter α, namely: (7) Fxα=12Axα−b2+12αRxα−x02=min F\left( {{{\boldsymbol{x}}_\alpha }} \right) = {1 \over 2}{\left\| {{\bf{A}}{{\bf{x}}_\alpha } - {\bf{b}}} \right\|^2} + {1 \over 2}\alpha {\left\| {{\bf{R}}\left( {{{\bf{x}}_\alpha } - {{\bf{x}}_0}} \right)} \right\|^2} = min which leads to the regularised normal equations of the form: (8) xα=ATA+α RTR−1ATb+α RTRx0, {{\bf{x}}_\alpha } = {\left( {{{\bf{A}}^{\rm{T}}}{\bf{A}} + \alpha \;{{\bf{R}}^{\rm{T}}}{\bf{R}}} \right)^{ - 1}}\left( {{{\bf{A}}^{\rm{T}}}{\bf{b}} + \alpha \;{{\bf{R}}^{\rm{T}}}{\bf{R}}{{\bf{x}}_0}} \right), where: A, b

– coefficient matrix and measurement data vector;

xα=f1α,f2α,…,fNα1×NT {{\bf{x}}_\alpha } = \left[ {{f_{\left( 1 \right)\alpha }},{f_{\left( 2 \right)\alpha }}, \ldots ,{f_{\left( N \right)\alpha }}} \right]_{1 \times N}^T

– vector of unknowns, whose components are equal to the reciprocals of the ultrasonic wave propagation velocity in the individual cells for a given value of α [s/m];

x0=cref−1,cref−1,…,cref−11×NT {{\bf{x}}_0} = \left[ {c_{ref}^{ - 1},c_{ref}^{ - 1}, \ldots ,c_{ref}^{ - 1}} \right]_{1 \times N}^T

– vector whose components are equal to the reciprocals of the reference ultrasonic wave propagation velocities c₀ [s/m];

R=[Rij]N×N {\bf{R}} = {[{R_{ij}}]_{N \times N}}

– regularisation matrix [−].

The structure of the regularisation matrix can be specified in various ways (Lampe & Voss, 2010; Polak & Mroczka, 2007; Voss, 2010). In the simplest formulation, the matrix R is taken to be the identity matrix. More generally, R can be used to impose desired properties on the regularised solution, provided that some a priori information about the solution is available. In the present study, in order to smooth the regularised solutions by reducing the magnitude of second derivatives in the velocity distribution maps constructed from x_α, the matrix R was chosen as a discrete second-order differential operator (Tatara, 2019).

The selection of an appropriate value of the regularisation parameter α, in the present study was based on the L-curve criterion (Lawson & Hanson, 1995; Miller, 1970). This is a heuristic parameter-choice rule, in which α is determined as a compromise between the residual term and the regularisation term of functional (7). In a log–log plot of the norm of the residual versus the norm of the regularised solution, the L-curve typically consists of two characteristic branches (Fig. 4). The first is a “flat” branch, where the error is dominated by the regularisation term (increasing α), and the second is a “steep” branch, where the error is governed by perturbations associated with the ill-conditioning of the system Ax = b (decreasing α). The optimal value of α is usually located in the vicinity of the corner of the L-curve (the point of maximum curvature), where the regularisation is already sufficient to stabilise the solution while still retaining most of the information contained in the residual term of functional (7). In some cases, more than one pronounced corner may appear on the L-curve (Skubalska-Rafajłowicz, 2011). In such cases, the corresponding solutions must be assessed from an additional perspective, for example by examining their physical plausibility or by analysing the data misfit, i.e. the magnitude of the discrepancy between the components of the vectors Ab and x.

Fig. 4.

Example of an L-curve (own research)

Ultrasonic testing was performed on a steel fiber-reinforced UHPC (Ultra-High Performance Concrete) beam with dimensions of 0.2 m × 0.24 m × 6.0 m (height × width × length). The exact dimensions of the beam under investigation are shown in Figure 5, and the experimental test setup is presented in Figure 6.

The beam was reinforced with conventional reinforcing bars and steel fibres of type DG 12.5/0.3-E304, with a length of 12.5 mm and a diameter of 0.30 mm. The steel fibre content in the beam was 2.1 %. The concrete had a compressive strength of approximately 238 MPa (Bońkowski, 2023). The longitudinal reinforcement consisted of four Ø12 bars at the bottom and four Ø12 bars at the top. The transverse reinforcement was provided by two-legged Ø8 stirrups with a centre-to-centre spacing of 150 mm.

Fig. 5.

Geometry of the analysed beam (dimensions in mm) (Bońkowski, 2020)

Fig. 6.

View of the experimental setup (own photo)

The beam was subjected to static loading in 13 cycles according to the scheme shown in Figure 7, with the maximum load being successively increased in each cycle. Measurements of the longitudinal wave travel time were carried out up to the 11th cycle. In the last two load cycles, severe signal scattering prevented a reliable determination of the travel times.

Fig. 7.

Force – displacement relationship for the UHPC beam. Table indicates maximum forces applied in the first eleven stages (own study based on Bońkowski, 2023)

Before the first loading (stage 0) and before each of the first 11 cycles (stages 1–11), the longitudinal wave travel time was measured for 14 transmitter-receiver point pairs (Fig. 8). The measurements were performed along a horizontal line of points spaced at 11 cm intervals, located 7 cm above the bottom edge of the beam, and were always taken after the beam had been completely unloaded.

Fig. 8.

a) Longitudinal view of the investigated element, b) Localisation of the measuring points in the UHPC beam (dimensions in mm) (own study)

Time-of-flight measurements were taken for up to two receiver points located diagonally with respect to each transmitting point. The longitudinal wave travel times between the predefined transmitter–receiver pairs were recorded using ultrasonic transducers mounted on opposite faces of the beam (Fig. 9). Transducers with a centre frequency of 150 kHz were used, and each measurement configuration was repeated three times.

Fig. 9.

Schematic of time-of-flight measurement for a single ultrasonic ray (Tatara, 2019)

Moreover, based on the measured travel times, additional “synthetic” time-of-flight values were generated by polynomial approximation at locations between the original measurement points. In particular, four equally spaced intermediate points were introduced between each pair of actual measurement positions.

The reconstructed velocity distribution maps, obtained by solving system (4) with a Tikhonov-regularised Least-Squares algorithm in which the regularisation matrix R implements a discrete second-order differential operator, are shown in Figures 10 and 11. The inversion was performed on a grid of 120 rectangular cells (5 in the vertical direction and 24 in the horizontal direction).

Fig. 10.

Tomographic reconstructions of longitudinal ultrasonic wave velocity distributions at successive stages: a) stage 0, b) stage 1, c) stage 2, d) stage 3, e) stage 4, f) stage 5 (own research)

To ensure that the number of equations in system (2) was not smaller than 120 (the number of unknown cell values), the tomographic images were computed using both the measured ultrasonic travel times and additional values obtained from polynomial interpolation between the measurement points arranged according to the layout shown in Figure 8.

Fig. 11.

Tomographic reconstructions of longitudinal ultrasonic wave velocity distributions at successive stages: a) stage 6, b) stage 7, c) stage 8, d) stage 9, e) stage 10, f) stage 11 (own research)

On the velocity maps shown in Figures 10 and 11, it can be seen that during the first two stages the beam exhibits slight stiffness degradation due to the onset of microcracking. From the 3rd stage onward, the reconstructed maps reveal the formation of localised damaged zones in the vicinity of the applied load, which become clearly pronounced after the 9th loading stage. Figure 12 illustrates the evolution of cracking in the investigated region over successive loading stages (up to and including stage 11), whereas Figure 13 shows photographs of the examined zone of the UHPC beam after completion of the ultrasonic measurements. Note that apart from ultrasonic tests, the beam was also tested with a dynamic actuator. This can explain the appearance of cracks in the upper part of the beam after initial load stages. As can be observed, visible cracks in the lower part of the UHPC beam appeared only after the 4th loading stage. The locations of the degraded regions identified on the reconstructed velocity maps are largely consistent with those revealed by visual inspection (Fig. 12).

Fig. 12.

Evolution of cracks in the investigated area of beam for each damage state: a) stage 1, b) stage 2, c) stage 3, d) stage 4, e) stage 5, f) stage 6, g) stage 7, h) stage 8, i) stage 9, j) stage 10, k) stage 11, source (own study based on Bońkowski, 2023)

Fig. 13.

View of the experimental setup (own photo)