*4.1. Calculation Example 1*

The primary purpose of the example is to show to what extent a change in stiffness in a localized, elastically degraded concrete zone can cause the paths of the fastest sound wave propagation to deflect. Another objective is to propose, on the basis of the performed analysis, a useful way to reduce the inaccuracies in tomographic imaging that may arise from this.

The example presents Dijkstra's method of calculating the paths and times of the fastest longitudinal ultrasound wave propagation between the assumed transmitting/receiving points in the longitudinal section of the concrete beam model. One damaged zone was assumed to be formed in the beam in its entire cross-section of 0.4 m × 0.4 m (these dimensions were arbitrarily selected as typical, but this does not change the general conclusions which were finally drawn on this basis). The damaged area is located in the middle of a beam section of 1.2 m in length. It was assumed that the brittle damage changes the distribution of wave propagation velocity along its length in the *x* function according to the results shown in Figure 6b. It means that a damage level is uniform in a given cross-section of the beam. Consideration was limited to cases where a minimum ratio of *E*D/*E*<sup>0</sup> in the damaged area was 0.9, 0.8, 0.6, or 0.2 with the development of the defect, and, due to the geometry of the task, imaging was performed on the section of the beam along its central vertical plane of symmetry. In order to simplify the calculations, the impact of any reinforcement on the results is not taken into account in the example. It should be obviously noted that this influence can be important for practical research and must not be ignored in the case of higher reinforcement ratios, in particular for the sound wave paths along the reinforcing bars. It is recommended to avoid such measurement situations as much as possible (e.g., Reference [1]). For this reason, it should also be emphasized that the conclusions resulting from the presented calculation examples can be used for quantitative analyses directly only for beams with a low reinforcement ratio and arrangement of propagation paths between the primary reinforcing bars and not along them and arms of stirrups, as well. In other situations, the impact of reinforcement should be taken into account—e.g., using appropriate correction factors, reducing the measured wave velocities so that they correspond to propagation conditions in non-reinforced concrete [1].

The assumed scheme of the beam with transmitting/receiving points is shown in Figure 8. The transmitting/receiving points were selected in the system of opposite points distant from each other every Δ<sup>p</sup> = 6.25 mm. Using Dijkstra's algorithm, wave propagation times and paths were determined between the opposite points and with shifting of the receiving points relative to the transmitting points right or left so that the angle of inclination of the straight line between them to the axis *x* was 45◦ or 135◦, respectively. In Figure 8, the tomographic rays corresponding to the wave propagation paths analyzed are also marked. Furthermore, in order to provide a convenient quantitative interpretation, the numerical results are presented in a dimensionless form, normalizing them to values that would have been obtained in the absence of the damage (i.e., in material characterized by longitudinal wave velocity equal to *c*L0).

**Figure 8.** Scheme of the beam model with the assumed damage distribution and the system of transmitting/receiving points and rays.

First, it was checked how the distance between the nodes of the graph (Δ<sup>n</sup> according to Figure 7) in Dijkstra's method influences the convergence of results. The graph nodes were placed in the transmitting/receiving points and inside the imaged longitudinal section in the intersection points of the orthogonal grid with the step of Δ<sup>n</sup> vertically and horizontally. In order to reduce the calculation time, only the damaged area (in which *E*D/*E*<sup>0</sup> < 1) was covered by the internal nodes. The weights of the graph edges were determined from relation (21) by integrating the speed distributions shown in Figure 6b with a step of 1 mm using the rectangular method. Figure 9 shows paths of the fastest propagation illustratively only between transmitting/receiving points distant by 5 cm from each other, min(*E*D/*E*0) = 0.2 and Δ<sup>n</sup> ≈ 100 mm, 50 mm, 25 mm, 12.5 mm, 6.25 mm, or 3.125 mm. Adoption of Δ<sup>n</sup> in this way caused that the number of graph nodes in the width of the damaged area changed in extreme cases from 3 to 70. The width of this area can be read from the diagrams in Figures 5 and 6. Table 1 shows the maximum relative changes of the determined wave propagation times on the fastest paths *t*path,*<sup>i</sup>* (s) and their lengths *L*path,*<sup>i</sup>* (m) as the step of the graph nodes decreases where *i* is the path index. When interpreting the shape of paths in Figure 9, it should be noted that if there is more than one path with the same wave propagation time due to the selection of points of graph nodes and symmetry of the problem geometry then the algorithm used shows the first one found among them. On the other hand, when analyzing the convergence of the solution (Table 1), it can be stated that, with the decrease in the step of the graph node grid, its relative changes become less and less significant and converge to zero. Due to the time needed to solve one such task and the owned computer hardware, further calculations were carried out using the solutions obtained at Δ<sup>n</sup> = 3.125 mm (PC with processor 3.06 GHz and RAM 8 GB allowing to determine data for one path on average in approximately 4 min at Δ<sup>n</sup> = 3.125 mm).

**Figure 9.** The shape of the paths of the fastest propagation between the selected transmitting/receiving points in the longitudinal section of the beam with the area damaged at min(*E*D/*E*0) = 0.2 and Δ<sup>n</sup> ≈ (**a**) 100 mm, (**b**) 50 mm, (**c**) 25 mm, (**d**) 12.5 mm, (**e**) 6.25 mm, and (**f**) 3.125 mm.

Figure 10 collectively shows the shape of selected paths of the fastest wave propagation depending on the degree of damage. It is evident that as the minimum ratio *E*D/*E*<sup>0</sup> decreases the wave undergoes increasing deflection in the damaged area. At the same time, it is obvious that the time needed to travel such a path by the wave is shorter and shorter than the one that would be calculated with a simplistic assumption of its propagation in a straight ray. The scale of this difference is illustrated in Figure 11 showing how the propagation times change along the fastest paths and rays connecting consecutive opposite transmitting/receiving points and points lying diagonally at an angle of 45◦ to each other in relation to the axis of the beam. Figure 11 does not show the variability of times along the paths between points lying diagonally to each other at an angle of 135◦ relative to the axis *x* because it is the same as for those lying at an angle of 45◦ taking into account the symmetry. In order to make the diagrams more readable, they are also presented in the form of continuous curves. In the case of a damaged zone with a course perpendicular to the axis of the beam analysed in the example, it can be observed that greater differences in wave propagation times along the real paths and calculated on the assumption of their straightness occur for those that connect the opposite points and are increasing as the damage increases. These differences are much smaller for paths connecting the points diagonally. This observation may be used to reduce the inaccuracy of tomographic calculations due to the fact

that *t*path,*<sup>i</sup>* (in fact, the times determined from measurements) are inserted in real situations in place of *t*ray,*<sup>i</sup>* in the system of Equation (3) where it is assumed in a simplistic way that *t*ray measured,*<sup>i</sup>* = *t*path,*i*. Therefore, in the case of propagation times only for opposite points, an approximate relationship can be proposed between *t*path,*<sup>i</sup>* and *t*ray,*i*, i.e.,:

$$t\_{\text{ray},i} \approx t\_{\text{ray},\text{approx},i} = t\_{\text{path},i} + \beta \left( t\_{\text{path},i} - \min \left( t\_{\text{path},i} \right)\_{i=1,2,\dots,l} \right) \tag{22}$$

where: *t*ray approx,*i*—approximated propagation times assuming a wave passage over the ray (s), β—non-negative factor [-].

**Table 1.** Maximum relative changes in path length and wave propagation times when changing Δ<sup>n</sup> in different defect evolution phases.


**Figure 10.** The shape of paths of the fastest propagation between the selected transmitting/receiving points in the longitudinal section of the beam depending on the degree of elastic degradation in the damaged zone min(*E*D/*E*0) = (**a**) 0.9, (**b**) 0.8, (**c**) 0.6, and (**d**) 0.2 (for Δ<sup>n</sup> = 3.125 mm).

**Figure 11.** Times of wave propagation *t*path,*<sup>i</sup>* on the fastest paths, and *t*ray,*<sup>i</sup>* over the straight rays between the transmitting/receiving points in the different phases of the defect evolution: min(*E*D/*E*0) = (**a**) 0.9, (**b**) 0.8, (**c**) 0.6, and (**d**) 0.2. The results are presented as a function of the position of the transmitting points: a black line for the fastest paths determined by Dijkstra's algorithm (for Δ<sup>n</sup> = 3.125 mm); a grey line for the paths established as straight rays. The diagrams on the left refer to the paths connecting the opposite points and those on the right to the points lying diagonally at an angle of 45◦ to each other in relation to the axis *x*.

Thus, the factor β scales propagation times from real paths by elongating them in relation to the minimum value in the direction of positive values so as to bring their course as close as possible to the times that would correspond to the passage of a wave over the straight ray. As a result, this will clearly result in a reduction of calculation errors. Analyzing the graphs shown in Figure 11, this factor can be determined, e.g., from the least squares method. However, in a situation of actual measurement, due to the lack of such data, this is impossible. Therefore, the authors propose to apply the following heuristic approach to the determination of β. Approximate solution of the system of Equation (3) can be alternatively obtained using Tikhonov's method [51,52] minimizing the function of the mean square error of this solution (the so-called residual term *F*res) with additional consideration of the regularisation term *F*reg (e.g., Reference [24]). Using this approach in the measurement situation, when in Equation (3) we replace the appropriate components in the vector **b** *t*ray measured,*<sup>i</sup>* = *t*path,*<sup>i</sup>* across *t*ray measured,*<sup>i</sup>* = *t*ray appox,*i*(β), we will obtain:

$$\begin{aligned} F\_{\text{res}} &= \left. \frac{1}{2} \| \mathbf{A} \mathbf{x}\_{a,\emptyset} - \mathbf{b}(\beta) \| ^2 \right., \ F\_{\text{reg}} = \left. \frac{1}{2} a \| \mathbf{R} \left( \mathbf{x}\_{a,\emptyset} - \mathbf{x}\_0 \right) \| ^2 \right. \\\ F\_{\text{res}} \{ \mathbf{x}\_{a,\emptyset} \} + F\_{\text{reg}} \{ \mathbf{x}\_{a,\emptyset} \} &= \min \rightarrow \mathbf{x}\_{a,\emptyset} = \left( \mathbf{A}^T \mathbf{A} + a \mathbf{R}^T \mathbf{R} \right)^{-1} \left( \mathbf{A}^T \mathbf{b}(\beta) + a \mathbf{R}^T \mathbf{R} \mathbf{x}\_0 \right) . \end{aligned} \tag{23}$$

where: α—regularization parameter [-], **R**—regularization matrix (m), **x**α,β—solution of the problem for given values α and β (s/m), **x**0—the point in the vicinity of which a regularized solution is sought (s/m). Matrix **R** was assumed in further consideration as a unitary one multiplied by 1 m for the sake of the unit account, and **x**<sup>0</sup> as the vector according to Equation (6). In the proposed method of determination of β, it is assumed that the less data contained in the vector **b** will be affected by measurement errors due to the selection of β, the smaller the values of the term *F*reg in the case of an optimal solution to the problem. It is, therefore, proposed to determine the optimal value of the coefficient β = βopt from the relationships:

$$\beta\_{\text{opt}} = \max\left(g(z)\right) \text{ for } z \ge 0,\tag{24}$$

where

$$\log(z) = \arg\min \left( F\_{\text{reg}} \left( \mathbf{x}\_{\alpha=z, \beta>0} \right) \right) \text{ for } z \ge 0. \tag{25}$$

Selection of βopt as the maximum value of function (25), in turn, has the following meaning: determine what is the maximum possible increase in the propagation time of the wave in the vector **b** using formula (22) with a minimum effect of regularization so that the original information contained in system of Equation (3) is lost as little as possible at the same time. Using formulas (22)–(25), the following values of βopt were obtained as in Table 2 for the data from the graphs in Figure 11. Illustratively, the variability of the wave propagation time *t*ray approx,*<sup>i</sup>* modified according to relation (22) between opposite transmitting/receiving points at β = βopt and min(*E*D/*E*0) = 0.2 is shown in Figure 12. At the same time, it was compared with *t*ray, *<sup>i</sup>* and *t*path,*i*, where *t*ray approx,*<sup>i</sup>* has a much more similar course to that of *t*ray,*<sup>i</sup>* than *t*path,*i*.

**Table 2.** βopt calculated on the basis of data from Figure 11 in case of the opposite transmitting/receiving points in different phases of defect growth.

**Figure 12.** Propagation times *t*path,*i*, *t*ray,*<sup>i</sup>* and *t*ray approx,*<sup>i</sup>* between the opposite transmitting/receiving points for a defect with min(*E*D/*E*0) = 0.2. *t*ray approx,*<sup>i</sup>* were calculated according to relation (22) with β = βopt.

In order to assess the correctness of the proposed approximation method, appropriate tomographic reconstructions of the longitudinal wave velocity maps in the beam model section from Figure 8 are presented in Figure 13. The reconstructions were determined by the randomized Kaczmarz method in accordance with the information presented in point 2, where *c*L0 was adopted as *c*L ref. The results are shown only for the central area of the beam separated by a red dashed line through which all types of rays passed due to their inclination. Wave velocity distributions in the damaged area according to Figure 6b were assumed with minimal ratio *E*D/*E*<sup>0</sup> equal to 0.9, 0.8, 0.6, or 0.2. The resolution of δ<sup>1</sup> × δ<sup>2</sup> = 3.33 cm × 3.33 cm was applied, as well as the number and arrangement of rays, as shown in Figure 8. For comparison purposes, in addition to the reconstruction with the use of times *t*ray approx,*<sup>i</sup>* according to relation (22) and *t*path,*i*, the original map is also shown on the basis of data from Figure 6b, but in the resolution used, i.e., in each cell, its average speed value is taken from the formula:

$$c\_{\text{L. mean},k} = \frac{\int\_{A\_k} c\_{\text{L}} \, dx \, dy}{A\_k} \, \tag{26}$$

where: *Ak*—rectangular area occupied by the *k*-th cell (m2). Integral (26) was calculated in an approximate way using the rectangular method, taking the step of integration after *x* and *y*, respectively as δ1/10 and δ2/10. Therefore, relative errors were calculated for each reconstruction: global—mean square *e*<sup>g</sup> = √( *<sup>K</sup> <sup>k</sup>*=1(*c*L mean,*<sup>k</sup>* <sup>−</sup> *<sup>c</sup>*L tom,*k*)2/ *<sup>K</sup> <sup>k</sup>*=<sup>1</sup> *<sup>c</sup>*<sup>2</sup> L mean,*k* ) and local—maximal *e*l max = max((*c*L mean,*<sup>k</sup>* − *c*L tom,*k*)/*c*L mean,*k*)*k*=1,2,...,*K*, where *c*L tom,*<sup>k</sup>* was determined from system of Equation (3). Because *c*L mean,*<sup>k</sup>* is not known in a real measuring situation, it is also necessary to introduce a measure which will allow to estimate the reconstruction errors due to the assumption of straight-lined rays without this information.

**Figure 13.** Originally assumed distributions of wave velocity in the longitudinal section of the beam in different stages of defect evolution (**a**) and their tomographic reconstructions according to Equation (3): (**b**) calculated using the propagation times *t*ray,*<sup>i</sup>* = *t*ray approx,*<sup>i</sup>* for the paths connecting the opposite points and *t*ray,*<sup>i</sup>* = *t*path,*<sup>i</sup>* for the paths connecting the points diagonally, (**c**) calculated using only the propagation times *t*ray,*<sup>i</sup>* = *t*path,*<sup>i</sup>* (red dotted lines—explanation in the text).

In this paper, it is proposed to calculate index *d*cL ray max: the maximal relative difference between *c*L ray measured,*<sup>i</sup>* (wave velocities averaged over the *i*-th ray based in real directly on the measured propagation times) and *c*L ray approx,*<sup>i</sup>* (wave velocities averaged over the *i*-th ray calculated, respectively, on the basis of approximated propagation times according to relation (22)). This leads to the formula:

$$d\_{\text{cl. ray max}} = \max\left(\frac{\text{\textdegree\textquotesingle}\_{\text{Lray measured},i} - \text{\textquotesingle}\_{\text{Lray approach},i}}{\text{\textquotesingle}\_{\text{Lray measured},i}}\right)\_{i=1,2,\dots,l} = \max\left(1 - \frac{t\_{\text{path},i}}{t\_{\text{ray approach},i}}\right)\_{i=1,2,\dots,l} \tag{27}$$

A summary of *eg*, *e*l max, and *d*cL ray max is presented in Table 3, depending on the adopted calculation strategy, while for *d*cL ray max , the exact value of it is given in a comparative manner (if calculated with the use of times of propagation *t*ray,*<sup>i</sup>* instead of *t*ray approx,*i*).

**Table 3.** Global—mean square (*e*g) and local—maximal (*e*l max) relative tomographic reconstruction error using the propagation times for the opposite points in Equation (3) *t*ray,*<sup>i</sup>* = *t*ray approx,*<sup>i</sup>* or *t*ray,*<sup>i</sup>* = *t*path,*i*. Maximal relative difference in average wave velocities over the rays *d*cL ray max due to the use of *t*ray,*<sup>i</sup>* = *t*ray approx,*<sup>i</sup>* or *t*ray,*<sup>i</sup>* = *t*path,*i*.


The analysis of the presented results shows that the reconstruction error increases with the degree of development of localized brittle damage. If propagation times *t*path,*<sup>i</sup>* are used in vector **b** of Equation (3) (as would be obtained directly from the measurement in the real situation), then, in the analyzed case, the relative decrease of Young's modulus in the damaged zone from 10% to 80% is connected with the increase of *e*<sup>g</sup> from approximately 0.002 to 0.09. In parallel, *e*l max changes from 0.01 to 0.70. Introduction to the vector **b** propagation times for the opposite transmitting/receiving points according to formula (22) with the optimum value of β allows reducing the relative reconstruction error of both *e*g, as well as *e*l max, for ranges from approximately 0.001 to 0.03 and from 0.007 to 0.19, respectively. The use of formula (22) also allows for a more accurate evaluation of the shape and spatial range of the defect (Figure 13), i.e., regardless of the degree of damage, its original shape is visualized correctly. Using only propagation times *t*path,*<sup>i</sup>* in the calculations allows to strictly assess the shape of the analyzed defect only when Young's modulus drops to about 20%. On the other hand, the value of *d*cL ray max also allows a rough estimation of the differences in the reconstructed wave velocity maps that may occur due to the adoption of the fastest straight-line propagation paths in calculations (although, taking into account the limitations that result from its definition according to Equation (27)). Its advantage in turn is that it can be determined in research only on the basis of measurement data. It starts to increase noticeably when Young's modulus drops in the defected zone above 20%.
