This paper considers the reconstruction of a defect in a two-dimensional waveguide using a derivative-based optimization approach. The propagation of the waves is simulated by the Scaled Boundary Finite Element Method (SBFEM) that builds on a semi-analytical approach. The simulated data is then fitted to a given set of data describing the reflection of a defect to be reconstructed. For this purpose, we apply an iteratively regularized Gauss-Newton method in combination with algorithmic differentiation to provide the required derivative information accurately and efficiently. We present numerical results for three different kinds of defects, namely a crack, a delamination, and a corrosion. These examples show that the parameterization of the defect can be reconstructed efficiently and robustly in the presence of noise.
Ultrasonic guided waves, especially Lamb-waves, are used in Non-Destructive Testing (NDT) and Structural Health Monitoring (SHM) to identify defects. Here the long testing range of the guided wave that can be covered is one of the main advantages over conventional scanning techniques. However, this advantage also yields new challenges. The multi-modal nature of guided waves with different wave velocities and the large distance between the defect and the sensor makes the localisation and characterisation of defects more difficult. On the other side, Lamb-waves provide the flexibility to choose the most appropriate mode, which is the most sensitive to a specific type of defect.
Most models for crack identification assume that the data containing the waves scattered by the defect are available at certain sensor positions. This data can be an electric signal, the displacement, or the velocity field depending on the type of sensor. We will refer to this data as measurement data, even though for many publications including this one, it is generated artificially by appropriate simulations. There are two main types of methods to detect a defect in a waveguide. One approach is based on imaging to localise and characterise the defects. Another class of algorithms fits a damage model directly with the measured data. The resulting damage model can then be further analyzed, e.g., using fracture mechanics techniques. This also supports other goals like NDT and SHM investigations and may answer the question of whether the damage is severe leading to a failure of the structure.
Most imaging algorithms assume a sensor network around the defect and consider a surface view of the waveguide as sketched in Figure 1. Often these approaches are called guided wave tomography or migration in the seismic community. The measurement data is then propagated backwards from the sensors on the surface. For one backward propagation method, each pixel of the surface is analyzed by the direct sound path between this pixel and transmitter/receiver-pairs [ihn2008pitch]. This approach relies on the assumption that there is only one dominant wave mode to use its group velocity for the backward propagation. Other methods do not need a dominating mode, for example, if a time-reversal approach is utilized as backward propagation method [mori2019damage]. There are many competing approaches for the backward propagation method. Since a full review is out of the scope of this introduction, we refer the reader to [zhao2011ultrasonic] and [neubeck2020matrix] for a comparison of the various methods. If the backward propagation method is defined by the transmitter/receiver-path, it depends strongly on the number of sensors and their position [zhao2011ultrasonic]. In other cases, the resolution of these approaches is directly linked to the physical properties of the propagating waves, for example, the wavelength [mori2019damage]. Most tomography algorithms show only a projection on the surface of the waveguide that does not show any depth of the defect.
As an alternative to the above imaging algorithms, the localization and characterization of defects can be viewed as an inverse problem. Here, one tries to find the set of parameters of a damage model that yield the smallest distance of the simulated data to the measurement data. These approaches consist of three main elements. First, one needs a forward model that generates simulated data representing the damage for the current set of parameters. Second, one uses an objective function that compares the measurement and simulation. Third, an optimization algorithm updates the parameter values to minimize the objective function, i.e., obtain a better fit of simulation and measurement. There are very different approaches with regard to the simulation method and the damage model, the objective function and optimization as discussed briefly in the next paragraphs.
For a discrete damage model, the extended finite element method (XFEM) is often used, because of its possibility to define mesh independent cracks [rabinovich2007xfem, rabinovich2009crack, jung2014modeling, agathos2018multiple, livani2018identification]. Rabinovich et al. [rabinovich2009crack]
used the XFEM in the frequency domain to reconstruct a single straight crack. The objective function is based on the least square error between measurement and simulation data. This objective function oscillates and has several local minima. A genetic algorithm optimizes this objective function to overcome the difficulty of getting trapped in a local minimizer. The same authors improved the results by switching to the time-domain[rabinovich2009crack] and changing the objective function. The objective function in the time-domain is based on the ”arrival time” of the reflection, which is often used in classic defect detection in an ultrasonic test. Using the arrival time, the objective function oscillates less, but there are still local minima. Livani et al. [livani2018identification]
investigated multiple defects with particle swarm optimization. Finally, three dimensional problems with multiple cracks are solved by Agathos et al.[agathos2018multiple]. They used a two-step GA approach for elliptical cracks. Jung et al. [jung2014modeling] consider curved cracks and inclusions together with a gradient descent as the optimization method. Multiple optimization runs are performed with different starting points representing different cracks to circumvent local minima.
A time-reversal approach can also be used to optimize defect parameters. If a sensor signal is reversed in time and transmitted back in the domain, it will be refocused at the source. Amitt et al. [amitt2014time] defined a objective function based on these refocusing of the wave in a specific part of a membrane111The word membrane should indicate that the Helmholtz equation is solved instead of the elastic wave equation.. Unfortunately, this objective function is oscillating [amitt2014time]. A full scan of the parameter space that defined the crack and its position was performed for optimization. [lopatin2017computational] considers an experimental validation of this approach.
Seidle et al. [seidl2016iterative] used another time-reversal approach for a membrane. Here, an integral damage model is used which lowers the stiffness for certain elements. The connection between XFEM and a stiffness-reduction can be found in [broumand2021inverse]. This approach, also known as Full Wavefield Inversion, allows many possible defect geometries and is not limiting their number, but the resolution of the defect depends on the computational grid. A highly refined computational grid can lead to very time-consuming calculations, especially for three-dimensional problems. A theoretical and experimental investigation of a similar approach can be found, e.g., in Rao et al. [rao2016guided, rao2017investigation, rao2020multi], where an acoustic inversion is used for elastic wave propagation to lower the computational demand. A change in both stiffness and density is considered in [rao2020multi].
The approaches mentioned above exploit the surface of the waveguide, whereas also a cross-section view may serve as alternative description. A plane strain simplification can be used to derive a two-dimensional problem - see Figure 1. Here, it is important to stress that these alternative dimensions of the cross-sectional view enable the analysis of the depth of the defect and the mode conversion due to defects with a certain depth. Wu et al. [wu2002inverse] analyzed composite plates with the strip element method. The model allows investigating horizontal and vertical cracks in the cross-section. Semi-analytical methods such as the strip element method have the advantage that the forward calculation is very efficient. This simulation approach is combined with a genetic algorithm to perform the optimization.
Gravenkamp proposed another semi-analytical method, the Scaled Boundary Finite Element Method (SBFEM), for simulating plate-like structures in the context of NDT and SHM [gravenkamp2012simulation]. The SBFEM was first developed by Song and Wolf to compute bounded and unbounded domains with the same approach [song2000scaled, wolf2000scaled]
. The domain is divided into super elements. Quasi-polar coordinates, called scaled boundary coordinates, construct this first type of super element. For problems in the frequency domain, the super element has only degrees of freedom on the boundary, which leads to a reduction of one dimension in the computational cost. Later, Gravenkamp et al. used another coordinate transformation to construct super elements with a constant cross-section for bounded and unbounded domains in two-dimensions[gravenkamp2012numerical, gravenkamp2015simulation]. These super elements are very efficient for simulating waveguides. The approach was extended to curved waveguides, three-dimensional-waveguides and waveguides with fluid/structure interaction [krome2017semi, krome2017prismatic, wasmer2018fluid]. In addition to the simulation of unbounded domains, a discrete crack model is an essential development in SBFEM [song2018review]. Quite recently, the SBFEM was used in a shape optimization scheme for a horn speaker and meta-materials [khajah2021shape] illustrating the flexibility of this formulation in the meshing process and and with respect to infinite boundary conditions.
In this contribution, we propose a general method for reconstructing a single discrete defect model inside a cross-section model of a waveguide. The SBFEM is used as a very efficient forward model for the simulation. The defect reconstruction is formulated as an optimization problem and a derivative-based optimization algorithm is used to solve it in a systematic and efficient way. The technique of Algorithmic Differentiation (AD, [GrWa08, Nau12]) is used to provide the required gradient information exactly and with low computational costs. We analyse the resulting approach with respect to properties of the objective function, the reconstruction quality and the robustness against noise.
summarizes the optimization scheme. This includes a discussion about the objective function and an estimate for the initial value. Three numerical examples are presented in section4 showing the versatility for different damage geometries. Conclusions are given in section 5. The appendix contains additional information on the implementation.
2 SBFEM forward model
This section introduces the forward model to simulate the guided waves. Throughout, we consider only the two-dimensional case. The extension to three dimensions is the subject of our current research. As mentioned in the introduction, the damage is incorporated directly into the waveguide model such that specific parameters model the damage. The forward model incorporates all important aspects of an ultrasound test on a waveguide. However, some aspects can only be designed for a specific experiment. Hence, certain simplifying assumptions have to be made to achieve greater general validity.
The first of these assumptions is that the defects are in a specific area within the waveguide modeled as an infinite domain. It is assumed that the vertical edges are far enough away from this area of interest so that the reflections from the vertical edges of the waveguide do not interfere with the reflections from the damage.
The second assumption is that the excitation by a sensor can be modeled by a traction force
on one part of the boundary. This traction is mapped into the frequency domain by the Laplace-Transformation. In contrast to the Fourier-Transformation, the imaginary part of the Laplace-Transformation should weaken the resonance and wrap-around effects of the numerical model. In the discrete-time version, this procedure, known as Exponential Window Method[kausel2017advanced], is given by
where denotes the complex conjugated of the previous term. The complex angular frequency is defined by a uniform frequency step and a small parameter . Note that only the real part of the complex angular frequencies changes. In practice, this can be simply derived by multiplying with the window function and then use the Fast-Fourier-Transform (FFT) on the product, where
is the resulting frequency step from the time vector and. For evaluation in the time-domain, the displacement is first transformed by the inverse FFT, and then the inverse window function is multiplied.
For all examples, the time dependent part of the excitation is a sin-modulated Gaussian pulse, i.e.,
with the center frequency , a time shift of . Figure 1(a) shows the pulse in the time-domain, while 1(b) illustrates the spectrum. Additionally, the dispersion curves of the waveguide are presented. This figure shows that the relevant spectrum of the excitation lies in the area, where only the two fundamental modes are propagating.
For each discrete angular frequency step, the SBFEM is used to approximate the displacement in the equations of linear elastodynamics
with the displacement , the linear stress and the density .
As for many other numerical methods, when using the SBFEM, the PDE in Equation (3) is converted into systems of linear equations
with the global, dynamic stiffness matrix that depends on the angular frequency, the nodal displacement vector and the nodal traction vector . To compute the global dynamic stiffness matrix , the domain has to be sub-divided into super elements, i.e., . Figure 3 shows such a subdivision into super elements. For each super element , parts of the boundary are meshed with finite elements that are used to compute a local stiffness matrix . These local stiffness matrices are then assembled into the global matrix, while the global nodal traction vector is assembled directly from the finite elements. The finite element approximation uses spectral shape functions based on the Gauss-Lobatto points. Figure 3 depicts elements with shape functions of degree with 5 nodes (marked with ). The corner points of the finite elements are dependent on a set of parameters that is later modified by the optimization process.
The prismatic super elements are used to approximate the extended parts of the waveguide. Figure 3 shows examples with the super elements , , and . These super elements require a constant cross-section. The cross-section is then meshed with finite elements and scaled in -direction. If the super element is finite, as in Figure 3, the finite element meshes must coincide on both sides. Note that the length in -direction does not influence the computational cost because there are only nodes at the cross-sections. The low number of nodes leads to a very efficient calculation for long parts of the domain. The outer boundaries parallel to the -direction, illustrated by thin lines in Figure 3, have traction-free conditions. Such super elements can also approximate semi-infinite parts of the model that are used to bypass reflections, see, e.g., and in Figure 3. For the derivation of the local stiffness matrix of the super element, see [gravenkamp2015simulation, gravenkamp2018efficient]
. The main part of the computation to set up the local stiffness matrix is an eigenvalue problem. A direct formula[van2007computation] computes the derivative of the eigenvalue problem inside the differentiated algorithm when applying AD.
The super element that forms a star-convex polygon is always required when the prismatic super element can no longer define the geometry or boundary conditions. Here, a finite element line mesh is scaled toward a single point. The single point is called the scaling center. In general, the position of the scaling center is also dependent on the parameter set to be determined by the optimization. The super element only needs nodes on the boundary yielding a high level of effectiveness. An additional special feature is that the super element also provides a simple crack model because a double node can introduce it. The crack is illustrated by in Figure 3, where a red node marks the double node. The red line is a traction-free boundary due to the double node. Despite the crack tip singularity, one still observes exponential convergence under -refinement [bulling2019comparison, tschoke2018numerical, bulling2019high] using such super elements. For the derivation of the local stiffness matrix of the super element, we followed [chen2014high]. The derivation of the local stiffness matrix is quite involved. A continued-fraction-based approach builds on one algebraic Riccati equation and several algebraic Lyapunov equations [song1997scaled, bazyar2008continued]. For this work, both algebraic equations are solved by an eigenvalue decomposition inside the SBFEM algorithm. The solution for the Riccati equation can be found in [song1997scaled], while the Lyapunov equation is solved in the appendix A.
It is worth mentioning that previous investigations have shown that polygonal super elements are quite robust for small and large finite element sizes inside the same mesh [xing2018scaled]. This robustness is important if geometric optimization changes the element sizes on the fly.
3 The Adapted Optimization Approach
The Inverse Problem
In general, the procedure described here to identify the material defect properties is an indirect method. The aim is to solve the inverse problem
where is given by physical measurements performed to detect the defect and denotes the Forward Operator. It is given by a simulation of the same quantity assuming that the material defect is characterized by some vector . The space and its elements represent one possible parameterization of the considered defect. For the examples considered in this article the space is scaled either to or . This scaling to the same range of values for all parameters supports the solution of the inverse problem. Solving the optimization problem (5), one obtains a parameterization for the material defect which best approximates the measurement data.
In practice, inverse problems are usually hard to solve for several reasons. Firstly, the resulting optimization problem need not (and generally does not) have a unique solution which may not depend continuously on the measurement data . Such inverse problems are called ill-posed problems and are pervailant throughout almost all applications of parameter identification problems. In order to solve ill-posed problems one is usually well-advised to deploy custom-tailored regularization techniques.
On the other hand, a different issue that may arise is that the error function in (5) may be hard to optimize as it may have adverse optimization properties such as not being sufficiently smooth or incorporating many only locally optimal points. In this case, it may be advantageous to modify the forward operator and the measurements so that the objective function is better suited for optimization purposes.
A simple approach uses the response data measured directly at a point. The -component is more interesting as the -component because of its practical relevance. For example, a single laser doppler vibrometer measures the -component of displacement or velocity, depending on the type of measurement method. We will assume that the displacement is measured. However, since both quantities have very similar characteristics, all results are transferable to the velocity. For all numerical examples, the first geometric parameter is the global position of the defect. Figure 3(a) shows in gray the y-displacement response of the model in Figure 3 at the point for selected parameters of , where the other parameters of the crack are kept constant at . Additionally the envelope of the signal is shown in black or dark green. Figure 3(b) depicts the objective function based on the response, where the green signal in Figure 3(a) is used as measurement data . However, as the waveguide is excited via a sinusoidal pulse, the measured response also incorporates the oscillating behavior. In fact, it can be seen that the objective function (5) also inherits these features, see Figure 4(a) for an illustration, which is turned zoom of the blue box in Figure 3(b). We note that the resulting objective function has many only locally optimal points near the global optima. Hence, many optimization methods would have difficulties solving this optimization problem to global minimality.
The authors also experimented with different variations of the data such as using the spectrum of the response. However, the problem caused by many local minima remained. Many test revealed that the objective function with the best properties from an optimization point of view was the calculation of the envelope of the response via a Hilbert-transform, see Figure 3(c) and 4(b). Analog to the response, Figure 3(c) shows the objective function for all crack positions between the green envelope in Figure 3(b) with a forward operator which included the calculation of the envelope.
In addition to the global minimum, Figure 3(c) shows two local minima. These local minima are associated with the mode conversion. The excitation, to the left of the blue line, is a pulse that corresponds to the S0-mode. However, the S0-mode and the mode converted A0-mode are reflected. If the S0-wavepackage of the forward simulation is at the position of A0-wavepackage of the measurement data, this leads to a local minimum. Compare the green line with the line above in Figure 3(a). Our approach to circumvent the local minima is to get a sufficient initial guess, which lies near the global minimum. The details of the initial guess can be found below.
For the rest of the paper, the measurement data is the envelope of the y-displacement at a single point, while the forward operator includes the calculation of the envelope.
In order to use gradient-based optimization methods for the reconstruction of material defects, derivatives of solutions of the partial differential equation are required. These could be approximated via finite differences (e.g.[jung2014modeling]) or computed by setting up and solving the adjoint equation of the partial differential equation considered (e.g. [seidl2016iterative, rao2016guided]). We tested the use of finite differences and found the approximation to be insufficient for the optimization procedure. Furthermore, in the development phase of this project we experimented with different boundary conditions. Hence, an adaptation of the adjoint equation as a time consuming and error-prone process would have been necessary in order to compute derivatives via the adjoint approach.
Therefore, we apply algorithmic differentiation, also called automatic differentiation (AD), which provides derivative information of arbitrary order for a given code segment. The derivatives are provably accurate up to working accuracy due to the fact that the computer program evaluating the function is broken down into a sequence of elementary evaluations upon which the chain rule of calculus is systematically applied. There are numerous AD tools available for a wide range of programming languages, see the web site www.autodiff.org for an overview of tools and references. Since our forward operator described in detail in the previous section is implemented in MATLAB, we apply the AD-tool ADiMat[adimat] to compute derivatives of solutions of the partial differential equations. Some mild modifications to the code were necessary such as computing derivatives of the Lyapunov equation which originally was not available in ADiMat.
The Resulting Optimization Procedure
Even though we made some effort to reformulate the objective function (5) so that it is easier to optimize, we cannot eliminate all issues entirely. Hence, in this subsection, we present the optimization algorithms that we deployed.
A first simple approach for the solution of the optimization problem (5) is to use standard optimization software naively. This can involve derivative-based approaches such as IPOPT [ipopt], WORHP [worhp], or the optimization routine fmincon
provided by Matlab. Alternatively, one may use derivative-free approaches such as coordinate descent, genetic optimization or other heuristic-based optimization methods like particle-swarm methods.
Derivative-based optimization algorithms have the great benefit of fast and provable convergence but generally only converge locally. With derivative-free methods often there is the hope associated that these methods converge globally. However, global convergence cannot be proven. Moreover, the whole optimization procedure may require tens of thousands of function evaluations. The required number of evaluations is of particular significance as the simulation in our case involves solving a partial differential equation; thus one objective function evaluation is already quite expensive – solving the partial differential equations tens of thousands of times is simply impractical.
In light of these facts, the approach presented here is a two-step approach: We first generate a refined initial guess for all variables. After a good initial point is established, we deploy a gradient-based optimization with fast convergence properties.
Initial guess for the global position
As mentioned above, the first geometric parameter is the global position of the defect for all numerical examples. To get an initial guess for the defect position , we consider the time of flight to get a physical position inside the waveguide. This physical position is then mapped linearly into the parameter space. The physical position of the defect can be estimated by
where is the time of flight, is the group velocity at the center frequency of the excited mode, and is the group velocity of the reflected mode with the highest amplitude. The computation of the time of flight is based on cross-correlating the excitation pulse, which is separated by a defined threshold from the reflections, with the rest of the signal. If the envelope of the cross-correlation is maximal, the argument is an approximation of the time of flight, i.e.,
where is the signal, is the envelope function, is the cross-correlation operator and is the Iverson bracket222The Iverson bracket is one if the statement inside is true and zero otherwise.. Figure 3(a) shows the threshold as a blue line, and the excitation pulse is on the left side while the rest of the signal is on the right side.
Refined initial guess for all parameters
In order to obtain a refined initial guess for the gradient-based optimization, we can use the approximation Equation (6). The approximation for the first coordinate given there is fairly accurate and is located within the valley of the objective function, see Figure 4(b). However, with a simple line-search in the vicinity of the valley, the accuracy of the initial guess can be greatly improved. This simple line-search requires very few function evaluations, and the computational cost is insignificant. The accuracy of the first coordinate is important as inaccurate values lead to locally optimal values (especially in the presence of noise) for the other coordinates with which optimization methods have difficulty escaping. The other coordinates are determined by selecting the vector with the smallest objective function from random vectors, where our default value of is 10. As the line-search is computationally very cheap, we also refine the other coordinates in the same fashion, beginning with the least critical or problematic coordinates. Due to many only locally optimal points, see Figure 10(a), the amount of randomly selected points is increased to 100 for corrosion-type defects. The whole procedure is sketched in Algorithm 1.
We experimented with different optimization methods in order to solve the optimization problem Equation (5). One difficulty we encountered was that the objective function has many locally optimal points that are difficult to escape from. However, optimization approaches aiming at a more global scale, such as genetic algorithms, or particle swarm methods are infeasible due to the computation effort caused by one evaluation of the forward operator.
We were able overcome these issues by applying an Iteratively Regularized Gauss-Newton Method (IRGNM), see e.g. [kaltenbacherbuch]. This method modifies the well known Gauss-Newton approach for nonlinear least-squares fitting problems by adding terms that make sure that matrices occurring in the Newton-step are invertible and adds a Thikonov-type term so that the overall solution stays near the initial guess. The overall procedure is given in Algorithm 2.
4 Numerical Experiments
In this section, we present several numerical examples to show the versatility of the approach. Each model considers another type of defect. The first example is about classic crack identification. The second example shows a similar set-up but tries to find a horizontal crack or rather a delamination defect. The third type of defect is a model of a corrosion defect.
All waveguides in the following examples are made out of steel. The isotropic material parameters are given in Table 1, and for all waveguides a plane strain is assumed. All meshes use shape functions of degree , i.e., elements with five nodes. The necessary degree is determined in advance by investigating the critical geometric parameter. For example, the smallest and the largest crack length are the critical geometric parameter for the first model. For all our models, a simulation with a degree of leads to a difference in the signals below compared to a simulation with . It is worth noting that there are numerous convergence studies for the SBFEM [tschoke2018numerical, gravenkamp2020mass, bulling2019comparison, chen2014high] and the code is thoroughly validated using commercial FEM tools.
4.1 Waveguide with a Crack
|global x-position of|
|x-position of the crack tip with 0 as the midpoint of|
|y-position of the crack tip with 0 as the midpoint of|
This example considers a steel plate with a single straight crack. Figure 3 shows an overview of the model. Note that the -axis is interrupted to fit the model into the figure. The first section of this figure is for the excitation. Arrows in Figure 3 show the area and direction where spatially constant traction is applied. At these parts, the excitation takes place by a normal traction force
where is the outer normal vector. The traction is a model for the excitation of a double transducer set-up [ihn2008pitch]. A double transducer set-up can produce both fundamental modes by applying a force in the same or opposite direction on the plate. For this traction, the S0-mode is excited, so the group velocity in equation (6) is S0-group velocity . Here, the S0-excitation is appropriate because the S0-mode leads to larger reflections for these kinds of cracks.
The second section of Figure 3 shows the evaluation point . The envelope of the -displacement is evaluated as a time series. As described above, the -component is of practical relevance as it can be measured using a single laser doppler vibrometer. The highest reflections for the -component are always associated with the A0-mode (see Fig. 3(a)), so the reflection velocity is the A0-group velocity at the center frequency - see Figure 1(b). Figure 4 shows the -displacement at point in gray and the envelope as a black curve over the signal of different crack positions.
The third section contains the defect. A double node is inserted into a continued-fraction-based super element. The emerging inner traction-free boundary is marked in red in Figure 3. Note that the continued-fraction approach handles the singular stress at the end of this crack model. In total, the model has degrees of freedom. This low number is possible because the semi-analytical approach approximates the large undamaged part of the waveguide. At this center frequency, the wavelength for the S0 mode is around . With other simulation methods, this leads to a much higher computational effort.
There are three parameters for the optimization. These parameters change the defect and its location. The first parameter defines the global position of the crack. This change in position is achieved by shifting the super element along the -axis and changing the length of the super element accordingly. The second and third parameters control the crack length and angle by moving the crack tip inside the super element relative to the boundary. The maximal and minimal values for the parameters are listed in Table 2. The crack length is between and ca. . Here the parameters are restricted to allow a numerically stable solution. Similar parameter restrictions can be found in other publications like [agathos2018multiple].
The objective function dependent on and with regard to the optimization problem (5) is displayed in Fig. 5(a). Here, is kept constant at , i.e., only orthogonal cracks are considered. The artificial target measurement corresponds to the midpoint of the parameter space. The surface is relatively smooth and has only a few local minima. The global minimum at has clearly the lowest value. Using the result of Algorithm 1 as refined initial guess, the defect properties can be obtained by applying the IRGN method Algorithm 2 without any need for using a more in-depth initial guess search such as applying genetic algorithms.
Another important aspect of an inverse method, especially in view of the measurement data, is its performance under noise. To analyze the performance, we create artificial measurement data by a single forward simulation with a parameter vector and add random noise to this target signal. The noise is defined in terms of the excitation pulse. The same threshold as for the time of flight is used to identify the excitation pulse. In Figure 4, the blue line shows this threshold. The excitation pulse is to the left of the blue line. A noise level ofof the maximum amplitude of the excitation pulse has been added to the target signal. The envelope of the target signal is calculated afterward. Figure 5(b) shows the reconstruction error per parameter for the midpoint of the parameter space. Here, the first parameter, the global position of the crack, shows the smallest error. A lower error in the first parameter is also expected because the range in the physical space is much larger compared to the other parameters – see Table 2. Translating the error in the physical space, the error stays under for a noise level under , which is a satisfactory level of precision for most applications. The error in the second parameter, on the other hand, shows lower robustness to noise. This parameter, which is linked to the crack angle, is likely to produce incorrect results. However, the more critical third parameter, which is more related to the crack length, has a significantly lower reconstruction error. In general, this investigation is only valid for the current target parameter ; as the crack length gets smaller with the third parameter, it is also expected that error increases due to the negative influence of the noise. In conclusion, there are also physical boundaries to an appropriate parameter space depending on the noise inside the experimental data.
Figure 6(a) shows the convergence of the reconstruction error with a small noise level of . This noise level is introduced to prevent an ”Inverse Crime” as both the target data and the reconstruction algorithm use the same simulation method. We observe an asymptotic rate that coincides with the rate of the regulation parameter in the Algorithm 2.
To further validate the inverse method, ten different target parameters are randomly drawn. Noise with a standard deviation of is added to the corresponding signals, and the parameters are reconstructed. Table 3 lists the reconstruction error for these points. For all reconstructions, the error is below .
4.2 Waveguide with a Delamination
|global x-position of|
|delamination length inside|
This example considers a steel plate with a single straight horizontal crack. This horizontal crack can also be seen as a delamination model. Figure 8 shows an overview of the model. Note that once again, the -axis is interrupted at some places to fit the model into the figure. The two colors inside the figure indicate that these defects are often present if there are two different layers. However, the material of these layers in our example is the same.
where is the second unit vector, excites an A0-mode. Hence, in equation (6) is the A0-group velocity at the center frequency - see Figure 1(b). An A0-excitation is chosen because it will lead to larger reflections for a delamination.
The second section in Figure 8 highlights again the evaluation point , where the envelope of the -component of the displacement is taken as the data. As for the previous example, the highest reflections for the -component are always associated with the A0-mode, so the reflection velocity is the A0-group velocity .
The third section of Fig. 8 contains the defect. The delamination is modeled by inserting double nodes in two super elements. The semi-analytical solution of the SBFEM once again captures the points with singular stress.
There are two parameters for the optimization. The first parameter denotes the global position of the delamination. The second parameter changes the length of the delamination. The maximal and minimal values for the parameters are listed in Table 4.
Figure 8(a) depicts the objective function dependent on and . As in the previous example, the target signal is generated with the midpoint of the parameter space for . Figure 8(a) and Figure 5(a) show a similar shape, but the two side channels with local minima are not present in this example. These side channels are missing because delaminations in the center of a waveguide do not lead to mode conversion for the A0-mode.
For this example, we could set the regularization parameter and achieve convergence within five iteration steps. In Figure 8(b), the error for reconstruction with noise is given. Both parameters show high robustness against the noise. Additionally, Table 5 lists the reconstruction error for ten different randomly drawn target parameters with a noise level of . The error indicates a high level of precision.
In total, the model has degrees of freedom. A single forward calculation is performed in under seconds on a modern desktop computer, while the reconstruction is computed in under minutes, including all preprocessing steps. The time is averaged over 100 reconstructions.
4.3 Waveguide with a Corrosion Defect
|global x-position of|
|corrosion length on the surface|
This example considers a steel plate with a simple model for a corrosion damage. The assumption is that parts of the steel are missing, and a tapering has formed in the waveguide. Figure 10 shows an overview of the model. Note that again the -axis is interrupted. The traction is the same as in first example (Equation (8)), which leads to an S0-mode excitation. The S0-mode shows greater sensitivity than the A0-mode to changes in thickness. Therefore, the S0-excitation is used. The two re-entrant corners [bulling2019comparison] in sub-domains and are located at the scaling center to handle the possible singular stresses for a deep corrosion. The model has degrees of freedom in total.
There are three parameters to be optimized. Once more, the first parameter defines the global position of the defect. The second parameter defines the length of the corrosion by changing the size of the super element . The third parameter denotes the depth of the corrosion - see Figure 10. The maximal and minimal values for the parameters are listed in Table 6. Figure 10(a) shows the objective function for a constant corrosion depth. The surface is more complicated due to the two main reflection points at the scaling centers, leading to mode conversion. As mentioned earlier, we had to increase the number of random vectors drawn for the initial guess because of the more complicated shape of the objective function. In Figure 10(b) the error for reconstruction with noise is given. Additionally, Figure 6(b) shows the convergence behavior for the reconstruction error at the midpoint of the parameter space. Despite the more complicated behavior of the objective function, the parameters show robustness against noise. For all ten tests, the reconstruction leads to small errors as listed in Table 7.
In the paper, we showed that a semi-analytical waveguide model given by SBFEM can be combined successfully with derivative-based optimization to reconstruct defects of different nature in a two-dimensional waveguide, even in the presence of noise. For this purpose, SBFEM is implemented in MATLAB, coupled with the AD tool ADiMat and embedded in an iterative solution procedure. The reconstruction can be performed on a standard desktop computer within a very moderate time such that our approach offers considerable runtime advantages compared to alternative methods. Three numerical tests were conducted to illustrate our approach: The first example covers a crack identification, the second one a delamination, and the third one a corrosion. In all three cases, the waveguides were made out of steel.
Future work will be dedicated to the extension to the three-dimensional case and anisotropic material behavior. The future work will allow the use of the proposed approach for non-destructive testing and structural health monitoring in real-life applications like carbon-enforced structures.
The authors gratefully acknowledge the German Research Foundation for funding (DFG project number 428590437).
Appendix A Implementation Details
This appendix shows how the Lyapunov equation can be solved with the help of an eigenvalue decomposition. The Lyapunov equation is
where is the conjugate transpose. Define the right eigenvalue decomposition of
with the right eigenvector matrixand the diagonal eigenvalue matrix . Note that this leads to the left eigenvalue decomposition of
After pre-multiplying with and post-multiplying with , Equation (10) leads to a new Lyapunov equation with a diagonal coefficient matrix
The matrix can be computed element-wise by
The solution matrix of the original Lyapunov equation can be derived by the inverse of Equation (14), i.e., .