On multiple crack identification by ultrasonic scanning

The present work develops an approach which reduces operator equations arising in the engineering problems to the problem of minimizing the discrepancy functional. For this minimization, an algorithm of random global search is proposed, which is allied to some genetic algorithms. The efficiency of the method is demonstrated by the solving problem of simultaneous identification of several linear cracks forming an array in an elastic medium by using the circular Ultrasonic scanning.


Introduction
Various optimization methods are widespread in the engineering science due to their efficiency when applied to both theoretical and practical problems, see, e.g., [1][2][3]. Classical regular methods are typically applied to optimize regular functionals, see [4]. If it is required to find the global minimum or maximum of a certain functional of a complex structure, then classical standard algorithms can hardly be used. However, some advanced algorithms of global random search demonstrate high efficiency just in applications to the problems with complicated target functionals, see, e.g., [5].
In the present paper, we demonstrate that many problems in engineering science can be reduced to certain classes of operator equations which can then be reduced to optimization problems. We construct an algorithm which permits efficient numerical treatment of the arising optimization problems and demonstrate the application of the proposed approach to the important problem of defect identification by using the Ultrasonic (US) scanning with an example, where the defects are represented by a finite array of linear cracks.
Many published works were currently dealing with inverse identification problems. This interest is caused by the importance of such investigations in many practical applications, in particular, as discussed above, in the US Nondestructive Testing. Chronologically, the first recognition methods were based on an approximation of weak wave interaction, namely on the theories of Born and Rytov [6]. The well-developed methods of acoustical tomography based on the Radon transformation are also related to the class of theories using the weak scattering hypothesis [7]. More strict diffraction theories are based on direct and inverse mathematical diffraction problems [8,9]. The present authors have published a series of works on the image identification of single defects [10,11]. Here we demonstrate that the optimization methods can also efficiently be applied to identify several defects simultaneously, on the basis of US methods.

Operator equations and optimization problems
Many problems in engineering science can mathematically be reduced to operator equations of the form where f is known but ϕ is unknown. In principle, Eq. (1) can be studied in a pair of Hilbert spaces, however for the applications discussed below it is sufficient to consider operator A acting from a finite-dimensional Euclidean space R N of dimension N to another Euclidean space R M of dimension M . To be more specific, we assume that D ϕ and D f are some simply connected convex bounded domains in R N and R M , respectively. Very often, the operator A is strongly nonlinear, and in such cases it is impossible to draw a conclusion about the existence and uniqueness of the solution of Eq. (1). However, in practice, a solution to Eq. (1) must be constructed even under the condition of such an uncertainty. Let us assume that there is a priori information that there exists a certain solution ϕ 0 to Eq. (1). Then the key point is to construct an algorithm for determining this solution or finding at least a good approximation to ϕ 0 . In the case where the operator A is continuous, even being nonlinear, an approximate solution ϕ * close to ϕ 0 corresponds to the respective right-hand side f * close to f . Let us introduce the following discrepancy functional: Obviously, it is nonnegative: Ω(ϕ) ≥ 0, and Ω(ϕ 0 ) = 0 since Aϕ 0 = f . Therefore, the minimization of the functional Ω(ϕ) may give a good approximation to the solution ϕ 0 . In the case where the operator A has a continuous inverse operator in the considered domain, the lower is the attained minimal value Ω min ≥ 0 of the functional, the more precise is the constructed approximation ϕ * ≈ ϕ 0 . It should be noted that, in nonlinear problems, the existence of the inverse operator cannot, as a rule, be proved as a theorem, but sometimes this fact is evident from the physical point of view. The important question regarding the uniqueness of the solution is beyond the present study.

Key properties of the applied algorithm
The functional (2) can be minimized by any classical optimization method [4]. However, the main restriction of regular iterative schemes is that they give only a local minimum of the respective functional. Under such conditions, it is not evident that a local minimum is simultaneously the global minimum of the functional. In fact, it is well known that, for nonlinear equations, such values of local minima may be too far from the desired value Ω = 0. The described property is demonstrated with an example of a one-dimensional functional which, in the one-dimensional case, is simply a normal function of the variable ϕ: This function is plotted in figure 1. This is a strongly oscillating continuous nonnegative function, which has several minimal zero values. These points arise every time when the function inside the modulus symbol in (3) changes its sign as the argument ϕ monotonically increases. If the analytical structure of function (3) is unknown, then one needs to apply an appropriate numerical algorithm for minimization. Obviously, if in the frames of any regular iteration process like the steepest descent method or Newton method, one starts from a certain initial point, then one arrives at a local minimum of the function at hand. Therefore, the global minimization cannot be resolved by such regular algorithms.
For this reason, in our numerical experiments, we use a version of the global random search method [5] contiguous to the one described in detail in [12]. This algorithm was developed to seek maxima, but it can also be used to seek minima. It is constructed so that it moves both up and downhill and, as the optimization process proceeds, focuses on the most promising area. As the first step, it randomly chooses a trial point within the step from the starting point selected by the user. The function is evaluated at this trial point, and its value is compared with its value at the initial point. In the minimization problem, all downhill moves are accepted and the algorithm continues from that trial point. The relationship between the initial value of Ω and the resulting step length is function dependent. This algorithm shows perfect convergence for many problems, also for our inverse identification problem, but unfortunately it sometimes arrives at a local extremum instead of the global one in the cases where there are a lot of global minima of the objective function. The algorithm applied is a slight modification of these ideas. It possesses the following two specific features: (1) random sampling of values in a neighborhood of the points at which the values of the functional are smaller happens more frequently than in a neighborhood of worse points, and (2) the domains in which random values of variables are chosen are gradually contracted to small neighborhoods of the points with smaller values of the functional. This technique demonstrates remarkable convergence for many known multidimensional functions of irregular behavior, like Rosenbrock's function, Rastrigin function, Himmelblau's function, and many other test functions. This turns out very efficient also in the engineering problem discussed below and dealing with crack identification by the Ultrasonic (US) scanning echo-method.

Formulation of the crack identification problem
Let a finite array of linear cracks be located in a linear isotropic homogeneous elastic medium inside a certain bounded domain; l n (n = 1, . . . , N c ) denotes the surface of the nth crack, and L = Nc n=1 l n denotes the full set of cracks. Let us choose the origin of the Cartesian coordinate system somewhere in this domain, close to the cluster of cracks. In order to identify the geometry of each crack in this array, we apply a circular scanning by an US sensor which operates in the echo-regime, see figure 2. The current position of the sensor is given by the Cartesian coordinates (R cos α, R sin α), and it is assumed that the amplitude of the back-scattered impulse is known for the full angular scanning: α ∈ (0, 2π).
Let us cite the governing equations in the problem under study. In the case of the so- called "anti-plane" (or "shear-stress," SH) problem, this stress is directed perpendicularly to the considered plane (x, y), so that the deformation mode is identical for all cross-sections z = const. Then in the fixed rectangular Cartesian coordinate system Oxyz, the displacement vector components areū = {0, 0, u z (x, y)}, and this determines the two nontrivial components of the stress tensor where µ is the shear elastic modulus. Under such conditions the equations of motion reduce to the single scalar Helmholtz equation for the function u z (x, y) where k is the wave number related to the transverse wave speed, the time-dependent factor e −iωt is omitted in all formulas, and c s is the transverse wave speed. The internal stress vectorT = {T x , T y , T z } over an arbitrary elemental area with the normaln = {n x , n y , 0} in the considered case of anti-plane deformation has the only nontrivial component Then the boundary conditions over the set of cracks L located in the medium should be satisfied over the faces of all defects which are free of load. Let us represent the full wave field as a sum of the incident and the scattered ones where α is the angle of incidence. It is obvious that the boundary condition for the scattered wave field implies that the crack faces are loaded by a certain tangential stress: In discrete formulation, in order to apply a numerical algorithm, all cracks are subdivided into I small elementary cracks. Then, by applying some known results of the classical potential theory, the diffraction problem can be reduced to a system of linear algebraic equations (SLAE) of the form I j=1 K lj g j = ik(n lx cos α + n ly sin α) exp[−i k(x l cos α + y l sin α)], l = 1, 2, . . . , I, where x ′′ l = (r j −r l ) ·t j , y ′′ l = (r j −r l ) ·n j .
Herer j is the radius vector of the jth elementary crack in the chosen Cartesian coordinate system,n j andτ j are the unit normal and unit tangential vectors to this elementary crack, ε j is its length, k is the wave number which is assumed to be fixed, and α is the angle of incidence with respect to the axis x. Besides, where H (1) 1 is the Hankel function. Once SLAE (9) is solved, i.e., all quantities g j are found, the amplitude of the back-scattered far-field wave can be calculated as A(α) = I j=1 g j tan[arccos(t jx cos α + t jy sin α)] sin k(t jx cos α + t jy sin α)ε j 2 . (14) Let the system under identification be an array of linear cracks l n , n = 1, . . . , N c . Each crack is defined by its central point with the Cartesian coordinates (a n , b n ), by its length ζ n , and by the angle of slope θ n (|θ n | ≤ π/2) with respect to the axis x. If the nth crack contains J n elementary cracks of length ε n = ζ n /J n , then the full set of elementary cracks is a union of ones when running over the given system of linear cracks: (q = 1, . . . , J n ; n = 1, . . . , N c ).
Coming to the formulated identification problem, let us estimate the total number of unknown parameters to be reconstructed. If the number of linear cracks is N c , then one has four unknown parameters for each of them: a n , b n , ζ n , θ n (n = 1, . . . , N c ). Therefore, the total number of the unknown parameters to be reconstructed is 4N c . In order to find all these unknowns, we construct an objective functional, and reduce the reconstruction to an optimization problem for this functional, as discussed in section 2. For this, we represent the SLAE (9) in operator form whose solution can be expressed in terms of the inverse matrix as Obviously, the right-hand side of Eq. (17) depends on all 4N c unknown parameters, as well as on the angle of incidence α: K −1 f = (K −1 f )(a n , b n , ζ n , θ n , α). Now, let us pass to the question: What is the measured information which can be used as the input data for the inverse identification problem? We assume that the scanning US sensor can measure the amplitude of the echo-impulse at a fixed distance R, in a far zone, for M positions of the sensor corresponding to the values of the irradiation angle α = α m (m = 1, . . . , M ). To be more specific, we assume that the values α m are uniformly distributed over the full circular interval (0, 2π). The registered values of the back-scattered amplitude form an array of the input data F m , (m = 1, . . . , M ). In our numerical experiments, the quantities F m may be taken as the amplitude A(α) calculated from the solution of the respective direct problem: F m = A(α m ), see Eq. (14). Then, by substituting (17) into (14), one comes to the system of nonlinear equations for the parameters a n , b n , ζ n , θ n written in discrete form

Numerical algorithm
The system of equations (18) can be resolved by minimizing the discrepancy functional, as discussed in section 2: min[Ω(a n , b n , ζ n , θ n )], Ω(a n , b n , ζ n , θ n ) = I j=1 (K −1 f )(a n , b n , ζ n , θ n , α m ) tan[arccos(t jx cos α m + t jy sin α m )] (K −1 f )(a n , b n , ζ n , θ n , α m ) tan[arccos(t jx cos α m + t jy sin α m )] The algorithm in such a form was realized in [13]. This shows a brilliant precision in reconstruction of the inclination angles θ and crack length ζ for each crack. The position of the cracks forming the array is identified incorrectly. This can easily be explained by the fact that the far-field scattered wave field, see expression (14), is free from the coordinates x j , y j and depends only on the angle of incidence. It is obvious that one cannot identify the real  position of a reflector by using only the information from the far-field zone. In order to improve the identification of the position of cracks, i.e., the parameters a and b, let us involve the information about the "time-of-flight" of the US impulse. The latter is always reflected in the echo-regime from the point which is the nearest to the sensor; hence this information gives the value of the most distant point on the cracks in each direction of incidence, i.e., the quantity d m , Therefore, it is reasonable to add an analogous discrepancy functional of the operator relation (20) to the functional (19): The final functional is thus constructed as a sum of expressions presented in equations (19) and (21) with appropriate weights. Numerous tests show that, to provide good precision for all geometric parameters under identification, the weight in front of functional (21) should be at least by an order of magnitude greater than the one for functional (19), For multiple crack array, the algorithm operates so that it selects the most likely geometry, sequentially for the number of cracks N c = 1, 2, 3, 4. Then, in the set of these four geometries, the algorithm chooses the best value of the functional Ω to come to the true solution including the recognition of the true number of cracks by itself.
Let us test the proposed method with some examples of the simultaneous identification of several cracks. If one applies US probes with the cyclic frequency f = 1 mHz and the wave speed in the medium is 6 km/s, then the wave length is λ = 6 × 10 3 /1 ×   Tables 1-4. All sizes are given in mm. The configuration of the system of cracks is shown in figure 3.

Conclusions
1. The number of cracks, parameter N c , for any configuration, as a rule, is identified correctly. In practice, this property takes place for both good and rough precision of the identification. 2. Following the conclusions of [13], one can see that, when ignoring the time-of-flight factor, i.e. with identification by minimizing functional (19) only and ignoring functional (21), the 9 1234567890 ''""  length of the crack (parameter ζ) and the inclination angle (parameter θ) are identified quite precisely. The Cartesian coordinates of the cracks cannot be correctly identified with this particular technique, which takes into account only the back-scattered amplitude of the US impulse. This is quite natural from the physical point of view when the impulse amplitude is measured in the far-field zone, where the value of the registered signal is insensitive to the real position of a defect. 3. It can also be seen from the presented Tables 1-4 that, by involving the time-of-flight information, one can attain more precise identification of cracks' location too. This is also physically natural, since the time-of-flight information is directly connected with the distance between a defect and the sensor. 4. It is very interesting to notice that even with an incorrect number of cracks, the identification, as a rule, contains the information about some true cracks. This can easily be seen on example of table 2, where the first crack among others, for "Three cracks" identification, has all four reconstructed parameters very close to the one for the first of the two real cracks. 5. It should be finally noted that in some extreme situations the identification cannot in principle be performed precisely. We mean the cases when a certain crack is hidden in an acoustic "shadow" of other cracks, like a crack located between a pair of two long parallel cracks.