**3. Localized Elastic Degradation in Concrete—Crack Model in Concrete According to the Damage Mechanics**

In concrete, due to its brittleness, the evolution of damage occurs in particular at the action of tensile stresses [27–31]. The process begins with the formation of microcracks which, with further increase in stress, grow into localized damage zones and cracks visible to the naked eye. For this reason, the presence of such zones at the final stage can be easily detected visually, as they contain cracks of the order of tenths of a millimeter in width. On the other hand, they are not visible at an early stage of development, and what is important, the initial microcracks usually do not occur in random places of concrete structure. When the cementitious material structure develops, the high and low internal cohesion zones can be distinguished. Especially, the places with low cohesion are the ones where cracks begin to grow because, in such places, even a small amount of released energy causes their propagation. They form mainly around the micro-pores or near the phase separation surfaces [11]. From the analyses presented in this respect within the damage mechanics, it is known that this process obviously leads to elastic degradation of concrete, which means a local decrease in material stiffness [27–31,49]. This phenomenon gives grounds for detection and control of this type of damage by ultrasound tomography if the spatial distribution of propagation velocity of a selected type of mechanical waves (e.g., Reference [24,33]) is assessed within such a framework. That is why, in the article, the preliminary calculations were oriented on determining the spatial distribution of concrete stiffness change in the case of localized damage under tension. The necessary information was obtained in this way in order to analyze the propagation of ultrasound waves in a localized elastically degraded concrete zone which evolution leads finally to forming macro-cracks. The calculations were performed using the assumptions of one of the most popular models in this respect, introduced by Chaboche [27] and Mazars [28], which takes into account the weakening of the material due to the isotropic damage accumulation and which was developed by Pijaudier-Cabot [29–31] in non-local terms. In a uniaxial tensile state, the model assumes the following relationship between stress and deformation (using the principle of strain equivalence):

$$
\sigma = (1 - D)E\_0 \varepsilon\_\prime \tag{8}
$$

where: σ—tensile normal stress (Pa); ε—normal strain [-] *E*0—Young's modulus in undamaged material (Pa); *D*—damage parameter [-]. Due to thermodynamic limitations of the process, the following must be satisfied: if *fD* <sup>=</sup> <sup>ε</sup> <sup>−</sup> <sup>κ</sup> <sup>=</sup> 0 and .

$$\begin{aligned} \text{if } f\_D = \overline{\varepsilon} - \kappa = 0 \text{ and } \overline{\varepsilon} > 0, \text{ then} \\ D = 1 - (1 - A\_{\mathfrak{t}}) \frac{\mathbb{x}\_0}{\mathbb{x}} - A\_{\mathfrak{t}} \exp(-B\_{\mathfrak{t}}(\kappa - \mathbb{x}\_0)) \text{ and } \dot{\kappa} = \dot{\overline{\varepsilon}}, \end{aligned} \tag{9}$$

$$\begin{aligned} \text{if } f\_{\rm D} = \stackrel{\cdot}{\overline{\varepsilon}} - \kappa < 0 \text{ or } \stackrel{\cdot}{\overline{\varepsilon}} \le 0, \text{ then} \\ \dot{D} = 0 \text{ and } \dot{\kappa} = 0, \end{aligned} \tag{10}$$

where: *A*t, *B*t—material parameters [-]; *f*D—load function [-]; κ—variable describing the process of material weakening [-]; κ0—initial value of the variable κ [-]; ε—non-local equivalent tensile strain [-]. To simplify the problem, if we assume that pure tension occurs in a concrete element with its axis described by variable *x* (m), the non-local equivalent strain at point *x* will be defined as:

$$\overline{\boldsymbol{\varepsilon}}(\mathbf{x}) = \frac{1}{V\_{\mathbf{r}}} \int\_{s\_{(1)}}^{s\_{(2)}} \overline{\boldsymbol{\varepsilon}}(\mathbf{s}) \, \boldsymbol{\psi}(\mathbf{x} - \mathbf{s}) A(\mathbf{s}) \, d\mathbf{s},\tag{11}$$

$$V\_{\mathbf{r}}(\mathbf{x}) = \int\_{\mathbf{s}\_{(1)}}^{\mathbf{s}\_{(2)}} \psi(\mathbf{x} - \mathbf{s}) A(\mathbf{s}) \, d\mathbf{s} \, \tag{12}$$

$$\psi(\mathbf{x} - \mathbf{s}) = \exp\left(\frac{-4\left(\mathbf{x} - \mathbf{s}\right)^2}{l\_\mathbf{c}^2}\right) \tag{13}$$

$$\overline{\varepsilon} = \sqrt{\sum\_{i=1}^{3} \langle \varepsilon\_i \rangle^2},\tag{14}$$

where: .ε—local equivalent tensile strain [-]; <sup>ε</sup>*i*—local principal strain [-]; *<sup>V</sup>*r—representative volume of the material (m3), *s*—variable describing the position along axis *x* (m); *s*(1), *s*(2)—starting and ending point of the considered element (m); *A*—cross-sectional area of the element (m2); *l*c—the characteristic dimension of the non-locality or the so-called internal length (m); ψ—weight function [-]; ...—Macauley's operator. Variability of function ψ is adopted in the model identical to the normal distribution with standard deviation *<sup>l</sup>*<sup>c</sup> 2 √ 2 . It also means that the range of the non-locality practically ceases to be significant above the distance *<sup>l</sup>*<sup>c</sup> because <sup>ψ</sup>(*<sup>x</sup>* <sup>−</sup> *<sup>s</sup>* = *<sup>l</sup>*c) = 1.83·10−<sup>3</sup> (Figure 2). In addition, in the case of uniaxial tension, .<sup>ε</sup> <sup>=</sup> <sup>ε</sup> where <sup>ε</sup> is the normal strain along the axis of the beam, <sup>κ</sup><sup>0</sup> will be equal to this deformation at the moment of initiation of the damage evolution, and κ<sup>0</sup> can reach a maximum value of *f*t/*E*<sup>0</sup> where *f*<sup>t</sup> is the tensile strength. In the incremental version needed for the numerical analysis of the problem, the stress and damage evolution Equations (8)–(10) take the following forms:

$$\begin{cases} f\_{\rm D} = \overline{\varepsilon} - \kappa = 0 \text{ and } \Delta \overline{\varepsilon} > 0, \text{ then} \\ \Delta \sigma = (1 - D)E\_0 \Delta \varepsilon - \Delta D E\_0 \varepsilon \\ \Delta D = \frac{\partial D}{\partial \overline{\varepsilon}} \Delta \overline{\varepsilon} = \frac{\partial D}{\partial \overline{\varepsilon}} \frac{1}{V\_{\rm r}} \int\_{\kappa\_{(1)}}^{\kappa\_{(2)}} \Delta \overline{\varepsilon}(s) \, \psi(x - s) A(s) \, ds \\ \qquad \qquad \qquad \Delta \kappa = \Delta \overline{\varepsilon} \end{cases} \tag{15}$$

$$\begin{cases} f\_{\rm D} = \overline{\epsilon} - \kappa < 0 \text{ or } \Delta \overline{\epsilon} \le 0, \text{ then} \\ \begin{cases} \Delta \sigma = (1 - D) E\_0 \Delta \varepsilon \\ \Delta D = 0 \\ \Delta \kappa = 0 \end{cases} \end{cases} \tag{16}$$

where: Δ ...—the finite increment of a given quantity.

**Figure 2.** Distribution of the weight function ψ for *x* = 0.

A separate issue related to the use of the discussed model is the adoption of appropriate material parameters that will enable the most accurate capture of the real behavior of concrete. The authors in paper [31] give the approximate typical ranges of these parameters in case of concrete of moderate strength, i.e., *<sup>E</sup>*<sup>0</sup> = 30–40 GPa, *<sup>A</sup>*<sup>t</sup> = 0.7–1.2, *<sup>B</sup>*<sup>t</sup> = 104–5·104, <sup>κ</sup><sup>0</sup> = <sup>10</sup><sup>−</sup>4, and *<sup>l</sup>*<sup>c</sup> = 3–5 *<sup>d</sup>*a max where *<sup>d</sup>*a max is the maximum size of the aggregate. However, so far, there are no exhaustive items in the literature devoted to research and formulation of automatic optimization calculation techniques allowing for precise estimation of all parameters mentioned above for a given type of concrete due to the length *l*c. It is assessed using mainly the scale effect and the model is calibrated by a manual trail-and-error method (e.g., Reference [31,49,50]). For this purpose, a suitable coefficient inverse problem is formulated in this article which uses illustratively the data from [34]. These tests concern the uniaxial tensile testing of a series of "dog bone" shaped concrete specimens of mean compressive strength *f*cm,cube = 50 MPa and *d*a max = 8 mm. The scheme of specimens is shown in Figure 3, and the experimental load-elongation dependencies (*P*exp-*u*exp) of a selected series of specimens of an overall length of 30 cm within a range of up to 50 μm are shown in Figure 4a. This relation was obtained on the basis of digitalization of graphs presented in Reference [34], and due to their quality by entering the curve in the middle area formed by the envelope of all several curves measured on the specimens of the length of 30 cm. For calculations, data from this series of samples were selected from all 7.5 cm, 15 cm, 30 cm, 60 cm, 120 cm, and 240 cm length series, which were tested in Reference [34], because in their case the relatively highest number of measurements was obtained with possibly small scatter of measured tensile strength. On the other hand, samples with a length of 7.5 cm were characterized by the largest scatter of measurement results on the basis of which the authors of [34] concluded that the width of their smallest section is smaller than the length *l*<sup>c</sup> and it can be larger. For this reason, its value was suggested as 6–7 *d*a max. The elongation measurement base was 0.6 times the notch length in the specimen (Figure 3). It should be noted at this point that a very interesting analysis of the development of elastically degraded zones based on the data from work [34] was also presented in Reference [50] but without formulating a coefficient inverse problem.

The tests included in Reference [34] were also selected for analysis due to the fact that the shape of the specimens enables, in tension, the location of the brittle damage zone in their middle area, leading to the formation of a single macro-crack when the load capacity is exhausted. In the considerations, for the sake of simplification, a bar in tension of variable cross-section was used as a model of specimens, according to their geometrical features. This choice was dictated by the need to carry out the calculation procedures described below in an acceptable time frame due to the computer hardware available (PC with a processor 3.06 GHz and RAM 8 GB—time of solving one task approximately 2 min). The boundary problem analyzed was solved using the finite element method (FEM) in an incremental version using physical Equations (15),(16). Two-node bar finite elements with 2 degrees of freedom and constant and averaged over the length physical and geometric features (with one integration point in the middle of the element) were used. One hundred and twenty-one finite elements were adopted on the axis of the specimen of the length of 30 cm. In addition, the load increments of the specimen were assumed to be continuously elongated during the calculation, which meant switching the force sign at the transition to the weakening phase. The force increment, in turn, was selected in each step so that the increments in elongation of the measuring base and normal strain in each of the elements would not exceed, respectively, the values of 0.15 <sup>μ</sup>m and 1.25 <sup>×</sup> 10<sup>−</sup>6. Taking into account the above mentioned conditions of calculations allowed to obtain a satisfactory convergence of the solution. Further increase of the number of elements and decrease of permissible increments of strain and displacements did not significantly affect the result (Figure 4b).

**Figure 4.** (**a**) Comparison of relations *P*exp-*u*exp and *P*-*u* using 121 finite elements. (**b**) Comparison of relations *P*-*u* using 121 and 241 finite elements.

First, *E*<sup>0</sup> was estimated by adjusting the slope of the initial straight-line section to the measuring point *P*exp, *u*exp = [20 kN, 5 μm] in the load-elongation relation (*P*-*u*) obtained from the model in the elastic range. The calculus procedure in this case was realized using ordered domain search. Hence, it was determined that *E*<sup>0</sup> = 38.29 GPa. Other parameters, i.e., *A*t, *B*t, κ0, and *l*<sup>c</sup> were estimated by minimizing the objective function:

$$F(\mathbf{p}) = \sum\_{i=1}^{n} \left( P\_i(\mathbf{p}) - P\_{\text{exp},i} \right)^2 \tag{17}$$

where: *P*exp,*i*—specimen loads from the experiment (N) determined in the weakening phase (Figure 4a) for elongations at equidistant intervals starting from *u*exp = 8.25 μm (for *P*exp = *P*max = 30.72 kN—i.e., maximum load) to *u*exp = 50 μm (for *P*exp = 8.18 kN); *Pi*—equivalents *P*exp,*<sup>i</sup>* determined on the basis of the assumed model (N); **p** = [*p*1, *p*2, *p*3, *p*4]—vector of variables corresponding to the parameters searched for, i.e., *A*t, *B*t, κ0, and *l*c, respectively. On that basis, it was assumed that:

$$\arg\min F(\mathbf{p}) = [A\_{\mathbf{t}\prime}, B\_{\mathbf{t}\prime} \kappa\_0, l\_{\mathbf{t}}].\tag{18}$$

The calculations in formula (17) assumed *n* = 31. The minimization of the implicit function (17) was carried out in three stages. In the first stage, a genetic algorithm was used on a population of 20 parameter vectors **p** in 20 selections. The vector components of the first population were drawn in preselected intervals: *<sup>p</sup>*<sup>1</sup> of [0.7, 1.2], *<sup>p</sup>*<sup>2</sup> of [104, 5 <sup>×</sup> 104], *<sup>p</sup>*<sup>3</sup> of [5 <sup>×</sup> 10−4, *<sup>P</sup>*max/(*E*0*A*min)] and *<sup>p</sup>*<sup>4</sup> of [3 *d*a max, 11 *d*a max], where *A*min [m2] is the minimum cross-sectional area of the specimen at the center of its length; hence, *P*max/*A*min = *f*t. It should be noted that the upper limit of interval selected for *p*<sup>4</sup> is greater than that given in Reference [31,34]. The authors assumed finally the value 11 *d*a max as the first smallest one for which the genetic algorithm did not estimate placing the minimum point of *F* right next to the upper limit of *p*4. Subsequent generations of the population were created as follows: **p** with the lowest value of *F* passed unconditionally to the next generation, the next 10 new vectors **p** were obtained by arithmetic crossover of randomly selected **p** from the group of the first 14 vectors of the old generation after their ordering from the lowest to the highest value of *F*. The set of the last nine new vectors **p** were drawn in the same way as in the first generation. This procedure was repeated independently five times. In the second stage, the same procedure as in the first stage was followed, however, changing the draw intervals of vector **p** components. The boundaries of them were defined from 0.8 to 1.2 of the the values of the **p** vector components for which the lowest value of *F* was obtained in stage one. In the third and final stage, an orderly search of the domain of acceptable solutions was performed in the vicinity of the point defined by the vector **p** components with the lowest *F*-value found in the second stage. The search interval was selected from 0.99 to 1.01 of the values of the parameters of this vector. If the minimum *F* in this area was at the boundary of any interval, the procedure was repeated where the point with the current minimum value *F* determined the midpoint of the intervals

in the next step. In this way, the following values of the parameters of the concrete damage evolution equations were obtained: *<sup>A</sup>*<sup>t</sup> = 8.01 <sup>×</sup> 10−1, *<sup>B</sup>*<sup>t</sup> = 1.95 <sup>×</sup> 104, <sup>κ</sup><sup>0</sup> = 5.64 <sup>×</sup> 10−<sup>5</sup> and *<sup>l</sup>*<sup>c</sup> = 112 mm while the objective function *<sup>F</sup>* reached the value 6.07 <sup>×</sup> <sup>10</sup><sup>6</sup> N2. This outcome corresponds to the global mean square relative error √( *<sup>n</sup> <sup>i</sup>*=<sup>1</sup> (*Pi*(**p**) − *P*exp,*i*) 2 / *<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *P*exp,*<sup>i</sup>* <sup>2</sup>) = 0.027 which expresses matching the experimental curve to the theoretical one in the weakening phase. In Figure 4a, the experimental curve and model one for the determined parameters were compared. Figure 4b also shows the course of the model curve after increasing the number of elements to 241 and reducing the permissible increments of measurement base elongation and normal strain to the value of 0.075 <sup>μ</sup>m and 6.25 <sup>×</sup> <sup>10</sup><sup>−</sup>7, respectively, to confirm the correctness of the calculations from the point of view of ensuring the convergence of the solution. In this case, the mean square relative difference between the solutions in the weakening phase with dense and rare discretization amounted to √( *<sup>n</sup> <sup>i</sup>*=<sup>1</sup> (*P*241,*<sup>i</sup>* − *P*121,*i*) 2 / *<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *P*121,*<sup>i</sup>* <sup>2</sup> ) = 7.6 <sup>×</sup> <sup>10</sup>−<sup>3</sup> where in the subscripts the number of elements used is given. On the other hand, Figure 5a shows the distribution of the damage parameter *D* at its different maximum values along the axis of the specimen at its middle section during the growth of the localized damage. It is also an equivalent way of modeling the development of macro-crack formation from the point of view of damage mechanics [30,49]. The *D*-distributions with a proportional 2-fold increase of all dimensions of the specimen model (Figure 5b) with the same values of parameters determining the accuracy of the solution were also calculated in a comparative way. A very similar variability of individual distributions was obtained, while the maximum width of the equivalent macro-crack zone reached a value approximately equal to 2*l*c. This confirms the correctness of the presented method of modeling the localized damage evolution in the case of tension in concrete. The presented considerations also justify the adoption of a specific resolution in tomographic imaging by means of wave velocity maps, i.e., dimensions δ<sup>1</sup> and δ2, according to Figure 1. In an extreme case, they should not be greater than 2*l*c, however, in order to ensure adequate image sharpness and precise identification of the degree of damage, they should be adequately smaller. Taking into account the variability of the distribution *D*, it is reasonable that δ<sup>1</sup> and δ<sup>2</sup> were assumed in the interval of approximately *l*c/3 − *l*c/5. The presented result also shows that the knowledge of the value of *l*<sup>c</sup> is crucial for the correct assessment of the distribution of damage around the crack in ultrasound tomography, and, on the other hand, it is the tomography that can be used for the direct assessment of this value.

**Figure 5.** Calculated distributions of *D* during the development of the damaged zone in the middle of the specimen: (**a**) model of the specimen *d* × *h* = 20 cm × 30 cm (**b**) specimen model with twice the total width and length increased *d* × *h* = 40 cm × 60 cm (*d*, *h* according to Figure 3).

Given that the growth *D* according to relation (8) causes the decrease of Young's modulus of the material, i.e.,: *<sup>E</sup>*<sup>D</sup>

$$\frac{E\_D}{E\_0} = 1 - D\_\prime \tag{19}$$

where: *E*D—tangential Young's modulus (Pa) of the damaged material during inactive growth of the damage, it is in its micro-cracked zones that the velocity of ultrasound waves decreases. Hence, based on the simplified isotropic damage evolution model for concrete developed in Reference [27–31] for three-dimensional stress state and assuming negligible change in material density during the evolution of brittle damage, the following estimated relationships can be given:

$$\frac{c\_{\rm LD}}{c\_{\rm L0}} = \sqrt{\frac{E\_{\rm D}}{E\_0}} \text{ or } \frac{E\_{\rm D}}{E\_0} = \left(\frac{c\_{\rm LD}}{c\_{\rm L0}}\right)^2,\tag{20}$$

where: *c*L0, *c*LD—the velocity of the longitudinal wave in the virgin and damaged material, respectively, (m/s). Based on relations (19), (20), variabilities of Young's modulus and longitudinal wave velocity are shown in Figure 6 which correspond to the damage parameter distributions from Figure 5b. They will be used in the computational examples presented in point 4 that illustrate the problem of tomographic detection of cracks in concrete members at various stages of their evolution.

**Figure 6.** Relative changes of Young's modulus (**a**) and longitudinal wave velocities (**b**) for distributions of *D* from Figure 5b.

It needs to be highlighted that the numerical results presented in this point can be directly used only in the case of pure tension in concrete. However, they can be also generalized usefully for the case of bent reinforced concrete (RC) beams when estimating the damage level in concrete of such beams with the following simplifying assumptions: the planar cross-sections of RC beam remain planar after loading (as in Bernoulli–Euler theory), failure of the concrete is determined by the normal cross-sectional stresses, and the average normal strains along beam axis are equal in the bonded reinforcing longitudinal bars and surrounding concrete. These assumptions are also commonly used in the design of RC beams, for example, as described in the EN-1992-1-1:2004 standard. Under the mentioned conditions, the results shown in Figure 6 can be also applied to a damage estimation of individual fibers of RC beam. The authors used the numerical results in this way for damage level estimation of experimentally tested RC beams in point 5.
