Defect reconstruction is essential in non-destructive testing and structural health monitoring with guided ultrasonic waves. This paper presents an algorithm for reconstructing notches in steel plates, which can be seen as artificial defects representing cracks by comparing measured results with those from a simulation model. The model contains a parameterized notch, and its geometrical parameters are to be reconstructed. While the algorithm is formulated and presented in a general notation, a special case of guided wave propagation is used to investigate one of the simplest possible simulation models that discretizes only the cross section of the steel plate. An efficient simulation model of the plate cross section is obtained by the semi-analytical scaled boundary finite element method. The reconstruction algorithm applied is gradient-based, and algorithmic differentiation calculates the gradient. The dedicated experimental setup excites nearly plane wave fronts propagating orthogonal to the notch. A scanning laser Doppler vibrometer records the velocity field at certain points on the plate surface as input to the reconstruction algorithm. Using two plates with notches of different depths, it is demonstrated that accurate geometry reconstruction is possible.

Characterizing damage in shell-like structures, such as plates, laminates, or pipes, is an important issue in non-destructive testing and structural health monitoring. Methods based on ultrasonic guided waves are appealing for damage characterization because the guided wave propagation allows the inspection of large areas of the structure from only a few sensor positions. However, determining the damage size using these methods can be challenging, as the various dispersive modes of guided waves do not always permit a straightforward correlation between reflections from damage and damage size.1–4 In addition, the wave packet reflected from the defect must be distinguished from the natural wave propagation in the pristine structure for many methods. Therefore, many state-of-the-art methods are based on residual signals that compare a measurement of the current, possibly damaged, condition with a baseline measurement of the pristine condition.5–7 In order to use these methods, the baseline measurement of the pristine condition must be available, and it must be guaranteed that the pristine condition is free of damage. This requirement can be restrictive, and baseline-free alternatives are desirable.

A possible alternative and baseline-free approach to solve this problem is an algorithm that compares the measured signals with those generated by a simulation model with damage of a specific size. If the simulated signal is close to the measured signal, the damage in the simulation model is a reconstruction of the actual damage.

The main objective of this paper is to validate an algorithm that is an extension of the concept published in Ref. 8, where the algorithm is tested with artificial data from simulation. The algorithm is able to reconstruct a defect defined by certain geometric parameters. The proposed algorithm is derived by solving an inverse problem formulated as an optimization problem. A forward model is iteratively updated until the output matches the measured signals as closely as possible.

The concept of the algorithm has two key aspects for efficient inversion. The first key aspect of the concept is the application of an efficient forward model. The forward model is based on the semi-analytic scaled boundary finite element method (SBFEM), first introduced by Song and Wolf,9–11 which results in a comparatively small number of degrees of freedom in the final system of equations. The second key aspect of the concept is the gradient-based optimization procedure. The gradient is computed by applying algorithmic differentiation (AD)12,13 to the forward model. By computing the gradient using AD, the optimization is fast and accurate.

This paper discusses how the algorithm can be extended to incorporate data from experiments and situations where the excitation has some undetermined properties. Here, the reconstruction of rectangular notches in steel plates using planar ultrasonic guided waves is demonstrated. Investigations of notches as a surrogate for natural cracks in guided wave studies have a long tradition.1–4,14,15 While the description of the proposed algorithm can be easily adapted to general problems, this artificial defect is used because it can be manufactured with high precision. The wave propagation and the type of notches investigated in the current work can also be mathematically reduced to a two-dimensional (2D) problem by considering only the cross section of the plate. The cross-sectional model is well studied in the guided wave literature (e.g., Refs. 1 and 2) and is a simple case with important features of guided wave propagation, such as dispersion and mode conversion. Reducing the mathematical model to the cross section often allows much faster computation, and the problem is more limited in the number of parameters. Consequently, specimens where the waves propagate in only one direction are often used for the first steps of testing an inverse method (e.g., Ref. 16). On the other hand, the experiment must be designed to fit the mathematical model of propagating waves. This paper, therefore, discusses the challenges of validating 2D models through experimentation.

Another favorable condition for the modeling of the experimental setup is a scanning laser Doppler vibrometer (LDV). This complex measurement system is currently required to make a qualitative comparison between the wave shape of the experiment and a simulation model. The advantage is that the model requires no calibration parameters other than the material parameters of the plate. A measurement with a LDV facilitates modeling. In contrast to a measurement with piezoelectric sensors, a LDV measurement can be treated as a simple point evaluation of the velocity, which does not affect the ultrasonic wave itself. Note that the performance of piezoelectric sensors can be highly dependent on the adhesive layer between the sensor and the plate.17,18 In Ref. 17, for example, significant changes are observed when the adhesive layer thickness varies in the range from 10 μm to 40 μm. In addition, the performance of the sensor-glue combination is temperature dependent (e.g., Ref. 19), although this temperature dependence is only significant for metal plates if there are large temperature fluctuations. Similar concerns are discussed in the qualitative comparison between the experiment and a simulation in Ref. 20. In order to create a simple simulation model that does not require any calibration parameters other than the material parameters of the plate, the excitation sensor and the adhesive layer are not directly modeled, but are replaced by a transfer function-weighted traction. This transfer function-weighted traction is the major extension of the algorithm required for the current experimental validation, because it automatically adjusts the ratio of the guided wave modes excited by the sensor to the measurement. The adjustment ensures independence from the exact parameters of the sensor and the adhesive layer. It should be emphasized that the validation of 2D models is a crucial step towards a general defect reconstruction method, and this example shows that geometric optimization is possible under laboratory conditions.

As an alternative to geometric defect optimization in the context of model-based reconstruction, defects can be modeled as material changes. This approach is known as full waveform inversion (FWI) and was pioneered for seismic waves.21 References 22–24 are examples of FWI for guided waves that include experimental validation. A similar idea can be found in Ref. 25, where wave velocities define the model parameters. Although these results are very promising, the methods require a numerical model that allows for a material change at each point where a defect might exist, making these models computationally expensive. In contrast, the method presented here takes advantage of the fact that large parts of the structure can be calculated semi-analytically by assuming a homogeneous material model, resulting in more efficient models.

As an alternative to gradient-based optimization of defect parameters, Ref. 16 presents statistical optimization of surrogate cracks in beam-like structures, while Ref. 26 considers pipes. The presented methods are very versatile and allow the reconstruction of multiple cracks, but they require a large number of simulations to obtain the statistical values. Both references consider one-dimensional (1D) models, where the cracks are approximated by special elements.

This paper is organized as follows. Section II summarizes the measurement setup for the two investigated plates. Section III defines the defect reconstruction algorithm and shows which extensions to previous work had to be made to adapt to the experiment. Section IV analyzes the performance of the algorithm and shows the reconstruction results. The conclusion is given in Sec. V. The  Appendix summarizes the details of the SBFEM derivation.

The reconstruction is studied on two 350 mm × 350 mm × 2 mm steel plates having a notch near the center. The notches are 0.5 mm wide and 300 mm long, while the notch depths differ with 0.4 mm and 0.8 mm, respectively. Figure 1 shows photos of the plates and an overview of the experimental setup, while Fig. 2 presents the schematic layout of the plates. The positioning of the origin of the coordinate system is arbitrary, but for convenience, the origin was chosen to be at the left edge of the notch. This makes deviations of the reconstruction from the nominal notch position trivial to calculate.

FIG. 1.

(Color online) Photos in the laboratory.

FIG. 1.

(Color online) Photos in the laboratory.

Close modal
FIG. 2.

(Color online) Overview of the LDV measurement of the plate with the original measurement grid of the case with the notch depth of 0.8 mm.

FIG. 2.

(Color online) Overview of the LDV measurement of the plate with the original measurement grid of the case with the notch depth of 0.8 mm.

Close modal

As mentioned in the Introduction, the reconstruction algorithm is based on a 2D cross-sectional model. However, a cross-sectional model simulates only plane wave fronts, which are constant orthogonal to the direction of wave propagation. Therefore, the experimental setup should favor the propagation of plane waves to achieve the best possible agreement between the simulation model and the measurement.

In summary, the data collection process can be divided into the following steps, which will be discussed in more detail below:

  1. The excitation signal is defined on a personal computer (PC).

  2. Based on step 1, the electrical input signal for the sensor is generated and amplified.

  3. The sensor uses the electrical signal to excite the ultrasonic waves inside the plate.

  4. A LDV measures the ultrasonic waves across a scanning area (Fig. 1).

  5. The measurement is post-processed to create the input for the reconstruction algorithm.

Step 1: Data acquisition begins with the definition of the normalized excitation signal s. Here, the normalized excitation signal is a sinusoidal Gaussian pulse, i.e., s ( t ) = sin ( 2 π f c t ) exp ( 0.5 ( t t c ) 2 / f c 2 ) with a center frequency, fc, of 500 kHz and a time shift, tc, of 5 f c 1 = 10 μ s. The choice of a Gaussian pulse is intended to promote agreement between experiment and simulation, since numerical models for smooth excitations achieve higher accuracy with the same number of degrees of freedom. Figure 3(a) shows the normalized excitation signal in the time domain, while Fig. 3(b) shows the amplitude in the frequency domain with renormalization. In addition, Fig. 3(b) shows the wavelengths of the S0, A0, and S1 modes of the plate. For the corresponding spectrum, only the S0 and A0 modes are excitable. The dispersion relationship between frequency and wavelength is calculated by the SBFEM.27,28

FIG. 3.

Time and frequency dependency of the excitation signal.

FIG. 3.

Time and frequency dependency of the excitation signal.

Close modal

Steps 2 and 3: The normalized excitation signal s is converted into an electrical signal using the arbitrary waveform generator (33521A; Agilent/dataTec GmbH, Reutlingen, Germany) and then amplified by a high-voltage amplifier (HVA-400-A; Ciprian, Grenoble, France). The amplified signal is sent to a macro fiber composite (MFC) sensor (M8503P2MFC; Smart Material GmbH, Dresden, Germany). The MFC sensor has an active piezoelectric area of 3 mm × 85 mm, and a total area of 7 mm × 94.5 mm. The long extension of the sensor in the y-direction compared to the wavelength of the propagating modes (about λ ̂ A 0 = 4.6 mm and λ ̂ S 0 = 10.4 mm) produces a nearly planar wavefront and thus to the wave type that also fits the model assumption of a plane strain state. The velocities in the y-direction are decoupled from the velocities in the x- and z-directions for plane waves propagating in the x-direction. Due to the choice of the sensor, shear horizontal modes associated with oscillations in the y-direction are negligibly excited, so only velocities in the x- and z-directions are considered.

Step 4: A scanning LDV (PSV-500-3D-HV; Polytec GmbH, Waldbronn, Germany) records the resulting ultrasonic waves. Due to the high reflectivity of the steel plate, it was necessary to treat the scanning area with white, non-permanent paint (ARDROX 9D1B; BASF, Ludwigshafen, Germany) to ensure a defused reflection of the laser (see Fig. 1). The scanning LDV measures the velocity field in all three directions in a predefined grid. For our application, the y-velocity field is neglected. To increase the signal-to-noise ratio, the signal is averaged over 1000 times for each measurement point. For this feasibility study, it was decided to test the algorithms on a very high data quality, in most cases 200 replicates are enough for metals (e.g., Ref. 20).

As mentioned above, the model simulates planar waves. Therefore, an evaluation of the planar waves from the experimental data is necessary for a good fit. Figure 4 shows that the excitation already produces a nearly planar wave. Nevertheless, some small changes in the y-direction can be observed. As a filter for waves with a wave vector not parallel to the x axis, the average of the velocities along the y axis is taken. In addition, Fig. 4 shows that the sensor setup excites the S0 and A0 modes simultaneously, which can be distinguished by their wavelength.

FIG. 4.

(Color online) Normalized velocity fields at t = 20 μ s for the original measurement grid (Fig. 2) on the plate with a notch depth of 0.8 mm.

FIG. 4.

(Color online) Normalized velocity fields at t = 20 μ s for the original measurement grid (Fig. 2) on the plate with a notch depth of 0.8 mm.

Close modal

Step 5: Since the automatic focusing of the lasers and the alignment of the plate lead to slight deviations from an ideal uniform grid, the x- and z-velocity fields are linearly interpolated to a smaller uniform grid. As Fig. 2 shows, the uniform grid has almost the same grid width and is only used to correct the angles and to reduce the number of measurement points. If the exact excitation is known, a measurement at a single x-coordinate is theoretically sufficient to reconstruct the defect parameters8 (compare also Ref. 16) However, since in our case the exact excitation is not sufficiently determined, we need several measuring points with different x-coordinates to additionally calculate the transfer function [Eqs. (8) and (9)] of the sensor. The reduced number of measuring points still allows a sufficiently accurate derivation of the transfer function, but also shows that the measuring points can be placed at a significant distance from the defect. The large size of the original scan grid is only required for Fig. 4. In addition, the new uniform grid enables averaging of the data in the y-direction, which, as mentioned above, comes closer to a plane strain signal. This results in a preliminary array V ̃ meas of size 2 × 11 × 9 × 782: here, the first dimension distinguishes between the x- and z-velocities, the second dimension is for the change in the x-coordinate, the third dimension is for the change in the y-coordinate, and the fourth dimension is for the change in time.

Finally, the data are normalized to the largest amplitude. In summary, this pre-processing of the preliminary array V ̃ meas is
V meas = nrm ( mean y ( V ̃ meas ) ) 2 × 11 × 1 × 782 ,
where nrm ( A ) = A / max ( vec ( env t ( A ) ) ) is the normalization of an arbitrary array A to its largest amplitude. The normalization is defined by an envelope operator over the time axis, env t ( A ), through a discrete Hilbert transform and a reshaping of a multidimensional array A into a vector, vec ( A ). To give an impression of the experimental data, Fig. 5 shows the mean velocity fields stored in the array V meas for both plates. The data for both plates were recorded using the same steps. The only difference is that the original measuring grid is slightly different because the scanning area is defined by hand in the LVD software. In Fig. 5, mainly the time interval between 33 μ s and 50 μ s shows the reflected waves coming from the notches.
FIG. 5.

(Color online) Normalized mean velocity fields for the interpolated uniform grid (Fig. 2) stored in V meas.

FIG. 5.

(Color online) Normalized mean velocity fields for the interpolated uniform grid (Fig. 2) stored in V meas.

Close modal
The defect reconstruction is formulated as an inverse problem. The geometric defect parameters q must be estimated so that the simulation model best fits the experimental data. The final reconstruction parameters are defined as
q min = arg min q Q F ( q ) y meas 2
with the forward operator F : q y sim K × 1 (see Sec. III A), a vector, y meas K × 1, containing the experimental data, and the parameter space Q.
Our previous research has shown that similar inverse problems have significantly fewer local minima if the difference between the experimental data and the forward operator is based on the envelope of the signals.8 The problem that an objective function based on the L2-norm of direct ultrasound signals has many local minima is also known as cycle skipping and has been discussed in detail in the FWI literature.29–31 Here, the most recent approaches rely on an objective function based on optimal transport.30,31 While an objective function based on optimal transport is an interesting alternative to using the envelope, it requires more computation time and has an additional parameter to be selected. In our previous study,8 using the envelope was sufficient for the optimization. However, our previous studies with an objective function based on envelope comparison evaluated a single displacement signal, while the current work considers the average velocity at multiple points. Nevertheless, we came to the same conclusion that this is an objective function with fewer local minima than the L2-norm of the direct signals, and the objective function is sufficient for optimization. Consequently, the vector y is defined by calculating the envelope over the time axis, and the array is reshaped into a vector, i.e.,
y meas = vec ( env t ( V meas ) ) ,
(1)
where V meas is the multi-dimensional array described in Sec. II, and thus, K = 2 × 11 × 782 = 17 204.

The forward operator mainly calculates the wave propagation in a cross-sectional model of a steel plate. The main idea is to use a semi-analytical model based on the SBFEM, which leads to a comparatively small number of degrees of freedom in the final system of equations. The SBFEM approximates the solution of the wave equations by dividing the domain into smaller subdomains called super-elements. The solution approach for each super-element is derived by a coordinate transformation that uses a finite element approximation of parts of the boundary. The current work uses the (complex) frequency domain to solve the elastic wave equation. For the frequency domain, the finite element boundary approximation contains all degrees of freedom of the approximation, which leads to an order of magnitude reduction in the degrees of freedom compared to a pure finite element approximation of the entire domain. In this work, two types of super-elements are considered: on the one hand, the high-order approach for polygonal elements,32 and on the other hand, the waveguide super-elements.33,34

As mentioned above, the problem is solved in the (complex) frequency domain to obtain a model with very few degrees of freedom. The exponential window method (EWM)8,35 is used to compute a finite domain in the frequency domain. The EWM is based on the discrete Laplace transformation, DLT ζ ( s ) = FFT ( s exp ( ζ t ) ), and its inverse, DLT ζ 1 ( s ̂ ) = FFT 1 ( s ̂ ) exp ( ζ t ), where FFT is the fast Fourier transform over the time dimension, s and s ̂ are arbitrary signal arrays in the time and frequency domains, t is the uniform time vector of the measurement, ζ = 0.5 ω Δ is the name-giving exponential factor, ω Δ is the resulting angular frequency step due to the sampling frequency of the time vector t, exp is the entry-wise exponential function, and is the entry-wise product along the time dimension. In addition to using these transforms, the wave equation in the EWM is derived for discrete complex angular frequencies: ω = ( ω Δ i ζ ) . Note that only the real part of the complex angular frequency changes with the index, while the imaginary part remains constant.

In our and many other experimental setups, there is no validated piezoelectric model of the sensor used. The current approach fits not only the damage to the experimental data but also the excitation model to account for this difficulty. For this fitting, the experimental data are transformed into the frequency domain
V ̂ meas = DLT ζ ( V meas ) ,
(2)
where V meas is the multi-dimensional array described in Sec. II. The excitation model is based on two tractions acting on the plate boundary at the sensor region with unknown amplitudes and phases. Therefore, the sensor itself is not modeled, but only its force. For each complex angular frequency ω in a certain frequency range with relevant amplitudes of the experimental data, the cross-sectional plate model (Fig. 6) approximates the two mean velocities v ¯ j , j { 1 , 2 }, one for each basic type of mean sensor traction,
τ ¯ 1 = ( 1 , 0 ) T and τ ¯ 2 = ( 0 , 1 ) T ,
(3)
and then, the two results v ¯ j are weighted by a transfer function to fit the experimental data in the frequency domain [see Eq. (9)]. The range between 10 kHz and 1.5 MHz is selected as the relevant frequency range for the current application. This range contains 75 discrete frequencies.
FIG. 6.

(Color online) Cross-sectional SBFEM model of the centerline (Fig. 2) of the plate with q = ( 0 mm , 0.8 mm ) T.

FIG. 6.

(Color online) Cross-sectional SBFEM model of the centerline (Fig. 2) of the plate with q = ( 0 mm , 0.8 mm ) T.

Close modal
For simplicity, we omit the index dependence of for the complex angular frequency ω = ω in the following. For the cross section at the centerline (Fig. 2), the mean velocity fields are v ¯ j = ( v ¯ j x , v ¯ j z ) T = ( i ω ) u ¯ j and are based on the elastic wave equations
( i ω ) 2 ρ u ¯ j ( x ) = · σ ( u ¯ j ( x ) ) , x Ω ,
(4)
σ ( u ¯ j ( x ) ) n = τ ¯ j ( x ) , x Γ sensor ,
(5)
σ ( u ¯ j ( x ) ) n = 0 , x Γ free ,
(6)
where ρ is the density, u ¯ j = ( u ¯ j x , u ¯ j z ) T the mean displacement, σ the mean linear stress operator, n the outer normal vector, and τ ¯ j the mean traction and the cross section coordinates are x = ( x , z ) T. The linear stress tensor σ is defined by Hooke's law, i.e., σ ( u ¯ j ) = C : ε ( u ¯ j ), and the linear strain-displacement relation, ε ( u ¯ j ) = 0.5 ( ( u ¯ j ) T + ( u ¯ j ) ), where C is the fourth-order elasticity tensor with the plane strain assumption given by the material parameters. The material parameters for the plates studied are density, ρ = 7 . 9 g cm 3, Young's modulus, E = 200 GPa, and Poisson's ratio, ν = 0.3.
The cross section at the centerline is the domain Ω with a parameterized notch, which can be expressed mathematically as
Ω = ( 180 mm , 170 mm ) × ( 0 mm , 2 mm ) ( q 1 , q 1 + 0.5 mm ) × ( 0 mm , q 2 ) ¯ ,
where q1 is the x-position of the notch, and q2 is the depth of the notch. The tractions are applied to the sensor area Γ sensor = [ 73 mm , 70 mm ] × { 0 }, while the rest of the boundary, Γ free = Ω Γ sensor , is traction-free.

For the SBFEM forward model (Fig. 6), the domain is subdivided into super-elements. For the current forward model, two different types of super-elements are used to calculate a local stiffness matrix, S ̂ b. The approximation approaches for both types, the waveguide super-element and the polygonal super-element, are briefly summarized in the  Appendix. The waveguide super-element is used for the long, undisturbed parts of the cross-section model, as it only requires degrees of freedom in the z-direction and can be of any length (see the  Appendix). The polygonal super-element, on the other hand, models all detailed areas of the waveguide. See the detailed views in Fig. 6. It is worth noting that no grid refinement is required for the stress singularities in the notch corners, as these are handled by the semi-analytical nature of the approximation approach (see the  Appendix).

The local stiffness matrices, S ̂ b, are assembled to a global stiffness matrix, S ̂. To approximate Eqs. (4)–(6), the two traction terms [Eq. (3)] must also be taken into account. Since the boundary on which the forces act is discretized with finite boundary elements, the nodal versions of the tractions F are represented analogously to the finite element method (FEM) by an N × 2 matrix, where N is the number of degrees of freedom, and each column is associated with a traction term in Eq. (3). Note that the nodal tractions F are frequency independent. A linear system of equations must be solved for each frequency:
S ̂ U ̂ = F with S ̂ N × N and U ̂ , F N × 2 .
(7)
For the model in Fig. 6, the number of degrees of freedom N is 346.
As mentioned above, the concept of a transfer function is utilized to generate the simulated velocity. Since the simulation model is constructed such that there are nodes at the measurement points, the velocities V ̂ sim for the two traction terms [Eq. (3)] are calculated by V ̂ sim = ( i ω ) A U ̂ with the assignment matrix
( A ) m n = { 1 if the n th simulation degree corresponds to the m th to the m th measurement degree , 0 otherwise , M × N .
The two columns of velocities are afterward fitted to the measurement data
v ̂ meas = vec ( V ̂ meas | ω ) M × 1 ,
where V ̂ meas | ω is the part of the array in Eq. (2) for the current frequency ω, with M = 2 × 11 = 22, by the transfer functions
h ̂ = ( V ̂ sim H V ̂ sim + β I ) 1 V ̂ sim H v ̂ meas .
(8)
A small regularization parameter, β = 1 × 1−9, is introduced, because the invertibility of V ̂ sim H V ̂ sim is not guaranteed. The final simulated velocities for the current frequency at the measurement points are
v ̂ sim = V ̂ sim h ̂ .
(9)
For all relevant frequencies, the simulation results are collected, reshaped, and transferred back into the time domain by the inverse Laplace transformation DLT ζ 1 in an array, V sim. Analogously to Eq. (1), the final output of the forward operator is
y sim = vec ( env t ( V sim ) ) .

Note that in the case under investigation, the sensor and the defect are located on different sides of the measurement grid. The derived transfer function can only reproduce the waves propagating in the positive x-direction and not the signals reflected by the defect propagating in the negative x-direction. This observation allows the transfer function and the defect parameters to be correctly determined at the same time. To create a similar situation where the position of the defect can be arbitrarily located, the LDV measurements could be performed on both sides of the sensor. Note that a defect between the sensor and the measurement points could lead to a mathematical ambiguity.

As mentioned in the Introduction, the forward operator in Ref. 8 is mainly extended by the calculation of the transfer function. The number of vectors given by the linear system of Eqs. (7) increases by a factor of 2, but no new square matrices corresponding to all degrees of freedom must be calculated. The 2 × 2 linear system of Eq. (8) and the matrix vector product in Eq. (9) add negligible computational complexity.

A simplified version of the forward operator is shown in Algorithm 1 with the matrices for each the super-element from the  Appendix. This version is intended as a guide for a first implementation. However, many steps can be taken to speed up this version. For example, quantities that are independent of the geometric parameters q, such as F, A, ω , and V ̂ meas, can be calculated only once and treated as additional input values. In our implementation, the frequency-independent matrices, e.g., M 0 , E 0 , E 1 , E 2, are stored for all super-elements that are independent of the geometric parameters q. These super-elements are to the left of x = 50 mm.

The reconstruction algorithm utilizes gradient-based optimization strategies. In contrast to typical derivative-free optimization methods, gradient-based optimization methods have provable convergence properties such as guaranteed local convergence with high convergence speed. However, such algorithms determine only local minima. Since the solution of the parameter estimation problem considered here is given by the global minimizer, for gradient-based methods it is important to provide a good initial guess which is in proximity of the globally optimal point. We achieve this goal by computing the error between simulation and measurement at 100 random points uniformly distributed throughout the parameter space Q. The point with the smallest error is selected as an initial guess for the optimization procedure.

As an optimization approach, we selected an iteratively regularized Gauss–Newton method (IRGNM),36 as shown in Algorithm 2.

In each iteration n, the IRGNM minimizes a linearization of the original nonlinear least-squares problem (the error between simulation and measurement) with an additional Thikonov regularization term that vanishes as n . The method is well understood and analyzed and is predominantly used in ill-posed parameter identification problems. Classical theory also guarantees convergence results in the presence of measurement errors. Different from in prior cases,8 when using synthetic data instead of real measurement data, it was not necessary to utilize regularization. Hence, we set α n = 0 for the current work.

This method is highly dependent on accurate and cheaply available gradient information. Derivatives could be computed numerically by applying a finite difference scheme to the forward operator. However, finite difference approximations of the derivative are usually very error prone and may be very costly to compute if the function to be differentiated has many variables. In the past, we have experienced some issues with accuracy in this context and thus chose to use a different approach: As the whole simulation code is available to us, we apply AD.12,13 AD as used here is a technique to automatically convert the code of programmed functions with input parameters into code that calculates both the function and the derivative of the function with respect to a given subset of input parameters. The basic principle of the conversion is to decompose the code into its basic commands and insert new lines of code that calculate the derivation with respect to the given parameters of the basic commands using the chain rule. We selected the AD tool adimat (see Ref. 37) in forward mode as it can differentiate simulations written in matlab. Furthermore, adimat can handle some technical programming issues such as data structures, etc., that other AD tools for matlab cannot. Overall, an examination of the computed derivative values showed that the derivatives were within negligible tolerance of expected derivative values. A matlab implementation of an algorithm similar to Algorithm 1 is then translated by adimat into a version that computes F ( q ) and the gradient F ( q ) with respect to the defect parameters q.

In this section, we present the performance of our reconstruction algorithm for the two studied plates with different notch depths. We will refer to the plate with a nominal notch depth of 0.4 mm as plate I, while the plate with a nominal notch depth of 0.8 mm will be referred to as plate II. The reconstruction is performed for a single parameter space that includes both notch depths. The position of the notch q1 is in the range between −40 mm and + 40 mm, which means that the distance is more than 17 A0-wavelengths. The notch depth is between 0.1 mm and 1.1 mm. The notch width is set to 0.5 mm, the nominal width of both plates. The notch width is not of much interest because the notch is only used as a replacement geometry for a crack. Normally, cracks are assumed to be infinitesimally wide, so no width needs to be reconstructed. It is essential that the numerical model is stable for the entire parameter space, since the initial values of the experiments are randomly drawn from this parameter space, as described in Sec. III B.

As a preliminary investigation to the optimization, the objective function is examined. This additional investigation only provides insight into the optimization problem but is not necessary for the actual reconstruction. Figure 7 shows the objective function for both plates over the entire examined parameter space. The colored contour plots are based on a parameter grid of size 161 × 21. Both plates show a complex shape for the objective function with several local minima. However, the global minimum is located near the nominal parameters q * = ( 0 mm , 0.4 mm ) T and q * = ( 0 mm , 0.8 mm ) T, respectively. There are several reasons why the global minimum of the objective function is not exactly at the nominal parameters. First, the measurement naturally contains a certain amount of noise. Also, the assumptions for the numerical model may not be accurate enough. It should be emphasized that this is only a 2D cross-sectional model. In addition, the plates were machined with a certain tolerance so that the actual values differ from the nominal parameters. To illustrate the machining-related variance, the statistical values of several gauging points at which the notch is assessed with a measuring microscope (MS4-BM01; Uhl Technische Mikroskopie, Aßlar, Germany) are listed in Table I. The distances were measured at 14 equally distributed points between y = 130 mm and y = + 130 mm along the notches (see Fig. 2). From the minimum and maximum depth values, it is evident that the bottoms of the notches are not completely flat, but there is no clear trend from one side of the notch to the other. On average, the notches are deeper than the nominal values. For the sake of simplicity, and since q1 varies over the length of the notch, the nominal parameter is usually referred to in the following. Furthermore, the MFC sensors were manually aligned to the predefined position and glued on by hand. The measuring points were aligned visually using the camera system of the LDV, so a more precise determination than the nominal position of q 1 cannot be given. Nevertheless, Fig. 7 shows that the reconstruction parameters in the area with a low objective function value are close to the nominal parameters.

FIG. 7.

(Color online) Contour plot of the objective function with the nominal parameters q * and the reconstruction parameters q min.

FIG. 7.

(Color online) Contour plot of the objective function with the nominal parameters q * and the reconstruction parameters q min.

Close modal
TABLE I.

Statistical data on the width and depth of the two notches, gauged with a measuring microscope.a

Plate Width (mm) Depth (mm)
Minimum Mean Maximum Minimum Mean Maximum
Plate I  0.498  0.504  0.507  0.408  0.427  0.450 
Plate II  0.496  0.500  0.505  0.763  0.809  0.850 
Plate Width (mm) Depth (mm)
Minimum Mean Maximum Minimum Mean Maximum
Plate I  0.498  0.504  0.507  0.408  0.427  0.450 
Plate II  0.496  0.500  0.505  0.763  0.809  0.850 
a

Gauged at 14 positions between y = 130 mm and y = + 130 mm (see Fig. 2).

The final reconstruction with the corresponding parameters q min (Table II) will be examined in the following. The initial parameter vectors are q 0 = ( 0.24 mm , 0.34 mm ) and q 0 = ( 0.07 mm , 0.59 mm ) for plates I and II, respectively. The final reconstruction is achieved after 17 and 13 iterations, where ϵ in Algorithm 2 is set to 1 × 1−5. Figure 8 shows the spectrum of the calculated transfer functions h ̂ i [Eq. (8)]. It can be seen that both tractions make a significant contribution to the signal. Moreover, the maximum value of the transfer functions is close to the center frequency of 500 kHz of the input signal s ̂. However, there are also qualitative differences between the two cases. While Fig. 8(a) shows an excellent agreement in the overall shape of the spectra between the transfer functions h ̂ i and the input signal s ̂, a clear difference in shape can be observed in Fig. 8(b). The transfer functions represent not only the different MFC sensors, but also the adhesive bond to the plate. The fact that the adhesive bond is handmade could explain the differences, despite the same design of the sensors and exactly the same signal generator and amplifier. However, the difference shows that calculating the transfer functions makes the optimization more robust to variations in sensor performance.

TABLE II.

Reconstruction parameters and reference values.

Plate Parameters Measured mean depth (mm)
Nominal Reconstruction
Plate I  q * = ( 0 mm , 0.4 mm ) T  q min = ( 0.09 mm , + 0.44 mm ) T  q 2 meas = 0.427 mm 
Plate II  q * = ( 0 mm , 0.8 mm ) T  q min = ( + 0.03 mm , + 0.77 mm ) T  q 2 meas = 0.809 mm 
Plate Parameters Measured mean depth (mm)
Nominal Reconstruction
Plate I  q * = ( 0 mm , 0.4 mm ) T  q min = ( 0.09 mm , + 0.44 mm ) T  q 2 meas = 0.427 mm 
Plate II  q * = ( 0 mm , 0.8 mm ) T  q min = ( + 0.03 mm , + 0.77 mm ) T  q 2 meas = 0.809 mm 
Algorithm 1.

Forward operator F ( q ).

   Input: geometric parameter q, time vector t, the measurement array V meas, regularization parameter β
Block of all frequency-independent quantities: 
1: compute ζ from t 
2: V ̂ meas = DLT ζ ( V meas ) 
3: allocate V ̂ sim, A, F 
4: compute A 
5: for all finite elements e do 
6:    compute F e 
7:    add F e to F 
8: for all super-elements do 
9:    compute and store M 0 , E 0 , E 1 , E 2 depending on the type of the super-element 
10:   if the current super-element is a polygonal super-element then 
11:     for m = 1,…,M do 
12:        compute and store S 0 ( m ) , S 1 ( m ) , X ( m + 1 ) 
   Block of all frequency-dependent quantities: 
13: determine a relevant frequency range ω from V ̂ meas 
14: for all ω in ω do 
15:    allocate S ̂ 
16:    for all super-elements do 
17:       compute S ̂ b depending on the type of the super-element 
18:       add S ̂ b to S ̂ 
19:     U ̂ = S ̂ \ F 
20:     V ̂ sim = ( i ω ) A U ̂ 
21:     v ̂ meas = vec ( V ̂ meas | ω ) 
22:     h ̂ = ( V ̂ sim H V ̂ sim + β I ) \ ( V ̂ sim H v ̂ meas ) 
23:     v ̂ sim = V ̂ sim h ̂ 
24:    sort v ̂ sim into V ̂ sim 
   Transform back to time domain and final steps: 
25: V sim = DLT ζ 1 ( V ̂ sim ) 
26: y sim = vec ( env t ( V sim ) ) 
   Input: geometric parameter q, time vector t, the measurement array V meas, regularization parameter β
Block of all frequency-independent quantities: 
1: compute ζ from t 
2: V ̂ meas = DLT ζ ( V meas ) 
3: allocate V ̂ sim, A, F 
4: compute A 
5: for all finite elements e do 
6:    compute F e 
7:    add F e to F 
8: for all super-elements do 
9:    compute and store M 0 , E 0 , E 1 , E 2 depending on the type of the super-element 
10:   if the current super-element is a polygonal super-element then 
11:     for m = 1,…,M do 
12:        compute and store S 0 ( m ) , S 1 ( m ) , X ( m + 1 ) 
   Block of all frequency-dependent quantities: 
13: determine a relevant frequency range ω from V ̂ meas 
14: for all ω in ω do 
15:    allocate S ̂ 
16:    for all super-elements do 
17:       compute S ̂ b depending on the type of the super-element 
18:       add S ̂ b to S ̂ 
19:     U ̂ = S ̂ \ F 
20:     V ̂ sim = ( i ω ) A U ̂ 
21:     v ̂ meas = vec ( V ̂ meas | ω ) 
22:     h ̂ = ( V ̂ sim H V ̂ sim + β I ) \ ( V ̂ sim H v ̂ meas ) 
23:     v ̂ sim = V ̂ sim h ̂ 
24:    sort v ̂ sim into V ̂ sim 
   Transform back to time domain and final steps: 
25: V sim = DLT ζ 1 ( V ̂ sim ) 
26: y sim = vec ( env t ( V sim ) ) 
ALGORITHM 2.

IRGNM.

  Input: initial value q 0, regularization parameter sequence ( α n ) n 0, maximal number of iterations N
1: for n = 0 , , N do 
2:     q n + 1 = q n + ( F ( q n ) H F ( q n ) + α n I ) 1 ( F ( q n ) H ( y meas F ( q n ) ) + α n ( q 0 q n ) ) 
3:    if q n + 1 q n < ϵ, then 
4:       q min = q n + 1 
5:      STOP 
  Input: initial value q 0, regularization parameter sequence ( α n ) n 0, maximal number of iterations N
1: for n = 0 , , N do 
2:     q n + 1 = q n + ( F ( q n ) H F ( q n ) + α n I ) 1 ( F ( q n ) H ( y meas F ( q n ) ) + α n ( q 0 q n ) ) 
3:    if q n + 1 q n < ϵ, then 
4:       q min = q n + 1 
5:      STOP 
FIG. 8.

(Color online) Transfer functions of the model for the reconstruction parameter q min in the frequency domain. Colors correspond to those in Fig. 6.

FIG. 8.

(Color online) Transfer functions of the model for the reconstruction parameter q min in the frequency domain. Colors correspond to those in Fig. 6.

Close modal

The optimization procedure should achieve the best possible agreement between simulation and measurement at the measurement points. To give an impression of the accuracy of the model, Fig. 9 shows a direct comparison of the signals. Although the inverse method only uses the envelopes to compute the objective function, the phase of the signal from the simulation also matches the measured signal. The matching phase can be explained by the calculation of the transfer function, where the phase information is adjusted.

FIG. 9.

(Color online) Comparison of z-velocities at point x = ( 60 mm , 0 mm ) T between forward models for reconstruction parameter q min (sim.) and the measurement (meas.).

FIG. 9.

(Color online) Comparison of z-velocities at point x = ( 60 mm , 0 mm ) T between forward models for reconstruction parameter q min (sim.) and the measurement (meas.).

Close modal

For the evaluation of the reconstruction, we refer to Fig. 10. In this figure, the SBFEM mesh of the reconstruction is plotted against the nominal geometry outlined by a dashed line. The figures show only a small section around the notch. It should be noted that the actual parameter space is much larger. The figure shows a typically numbered x axis at the bottom, while an x axis measured with the A0-wavelength of the center frequency fc is shown at the top. This additional axis should make it possible to estimate the error for other plate dimensions, since the wavelength is the characteristic length of the problem. The errors for the x-position | q 1 * q 1 min | of the notch are + 0.09 mm and + 0.03 mm for plates I and II, respectively. In relation, the nominal errors for the notch depth | q 2 * q 2 min | are + 0.04 mm and + 0.03 mm, while the errors to average measured notch depth | q 2 meas q 2 min | are 0.01 mm and 0.04 mm.

FIG. 10.

(Color online) Reconstruction model based on q min compared to the outline of the nominal notch based on q *, shown with a dashed line.

FIG. 10.

(Color online) Reconstruction model based on q min compared to the outline of the nominal notch based on q *, shown with a dashed line.

Close modal

An algorithm has been presented that enables defect reconstruction using guided ultrasonic waves by solving an inverse optimization problem. The output of a simulation model, which contains a parametric representation of the defect, is intended to approximate the experimental data as closely as possible. For this purpose, a gradient-based method is used, where the computation of the derivatives according to the defect parameters is based on AD. A second aspect of the presented work is the highly efficient ABFEM, which requires very small, final linear systems of equations to approximate the linear waveguide equation.

This article presents the successful reconstruction of two notches in steel plates based on data from experiments. Although the reconstruction algorithm is formulated in a very general way, the initial focus was on the simplest case of wave propagation, where the ultrasonic wave can be assumed to be planar, and the problem can be reduced to a 2D cross-sectional model. An error between nominal geometry and reconstruction on the order of 0.06 mm was achieved. Thus, it could be shown that, first, a geometric reconstruction of defects is possible and, second, a cross-sectional model is capable of describing the wave field with sufficient accuracy for a reconstruction, given a suitable experimental setup.

The first important finding of this work is that the envelopes of the recorded data provide sufficient information while leading to an objective function that allows optimization with a gradient method. It should be noted that the current sensor setup leads to mixed mode excitation, and both modes are converted into each other at the defect. In addition, the modes are subject to dispersion. Nevertheless, the notch parameters are correctly determined by the approach because these effects are taken into account within the model.

The second finding is that another linear fit must be introduced for the sensor's force distribution. This allows reconstruction even for an ultrasound source that is not fully characterized, where the ultrasound source is determined by both the sensor and the adhesive layer.17,18 However, real-world applications would benefit from a sensor model that includes the adhesive layer and allows quantitative comparison between simulation and measurement with a small error—eliminating the need for LDV measurement.

Future work will extend the reconstruction algorithm to three-dimensional (3D) models of wave propagation. Reconstructing a crack in three dimensions introduces the y position of the crack, the angle of the crack in the xy plane, and the crack length as important parameters. This leads to new challenges. For example, the space of possible initial values is at least five-dimensional (5D), so a random draw is unlikely to come close to the true crack parameters. In addition, if there is only one source of ultrasonic waves, the crack may be orthogonal to the wavefront, resulting in very small reflections. This makes a single source technique unlikely to detect all cracks.

The authors gratefully acknowledge the German Research Foundation for funding (DFG Project No. 428590437). In addition, the authors would also like to thank Yevgeniya Lugovtsova and Tobias Homann for their support with the LDV measurement and Andreas Knöppchen for the data obtained with the measuring microscope.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This appendix briefly summarizes the derivation of the two types of super-elements used and introduces the notation of the matrices for Algorithm 1.

1. The waveguide super-element
The first type of super-element approximates the displacement in the undisturbed parts of the plate's cross section. Gravenkamp et al. developed this super-element for waveguides;33 therefore, it will be referred to as waveguide super-element in the following. Figure 11(a) shows an example of a super-element, where γ ( η ) is a 1D finite element mesh of the left and right edges in the z-direction with an element-wise parameterization coordinate, η. The coordinates x and xr are the positions of the left and right boundaries, respectively. Using the coordinate transformation
( x z ) = ( ξ + x γ ( η ) ) ,
with ξ between 0 and the super-element length L = x r x , a product ansatz,
u ¯ ( x , z ) = N ( η ) u ̂ b ( ξ ) ,
(A1)
is assumed, where N ( η ) is the matrix of the vector-valued shape functions defined by the finite element mesh and u ̂ b ( ξ ) is the vector of nodal displacements.33,34 While in general a variety of different shape functions can be adopted,38 for this work, spectral shape functions of degree p = 6 (with 7 nodes) based on the Gauss-Lobatto integration points are used. Substituting the product ansatz Eq. (A1) as the test and ansatz functions into the weak form of Eq. (4) and integrating by parts, we obtain the ordinary differential equation (ODE), called the SBFE equation in displacement,
E 0 ξ ξ u ̂ b ( ξ ) + ( E 1 T E 1 ) ξ u ̂ b ( ξ ) E 2 u ̂ b ( ξ ) + ω 2 M 0 u ̂ b ( ξ ) = 0 ,
(A2)
where the coefficient matrices E i and M 0 are calculated similarly to the FEM by integrating the shape functions and their derivatives. The ODE Eq. (A2) leads to a quadratic eigenvalue problem,
E 0 ( i k ̂ ) 2 ψ ̂ + ( E 1 T E 1 ) ( i k ̂ ) ψ ̂ E 2 ψ ̂ + ω 2 M 0 ψ ̂ = 0 ,
where the eigenvalue ( i k ̂ ) is the product between the wavenumber of a mode and the imaginary unit, and the eigenvector ψ ̂ is the nodal mode shape vector. The solution for the nodal displacement has the form
u ̂ b ( ξ ) = Ψ ̂ diag ( exp ( i k ξ ) ) c ̂ + Ψ ̂ r diag ( exp ( i k r ( ξ L ) ) ) c ̂ r ,
where Ψ ̂ , Ψ ̂ r are the eigenvector matrices with columns ψ ̂, k ̂ , k ̂ r are vectors with entries k ̂, and A = diag ( a ) is the diagonal matrix with entries according to the vector a. The coefficient vectors c ̂ , c ̂ r are the unknowns of the wave field for this super-element and must be determined by the boundary values. The index indicates the waves propagating from left to right, while the index r indicates the waves propagating from right to left. A rearrangement of the equation for the displacements u ̂ b ( ξ ) and an equation for the internal nodal forces ( E 0 ξ + E 1 T ) u ̂ b ( ξ ) with ξ = 0 and ξ = L then allows determination of the local stiffness matrix S ̂ b for this super-element. For a detailed derivation, the interested reader is referred to Refs. 33 and 34. We want to emphasize that the length L of the waveguide super-element does not affect its computational cost, making it an ideal choice for the long, undamaged parts of the cross section.
FIG. 11.

The two types of SBFEM super-elements.

FIG. 11.

The two types of SBFEM super-elements.

Close modal
2. The polygonal super-element
The second type of super-element approximates the displacement for a star-convex polygonal sub-domain. Figure 11(b) shows a 1D finite element mesh, γ, and a special point, x c, called the scaling center. Similar to the first type of super-element, a coordinate transformation,
x = ξ ( γ ( η ) x c ) + x c ,
with ξ between 0 and 1, the product ansatz form Eq. (A1), and integration by parts are utilized to write the weak form of the partial differential Eq. (4) into another sealed boundary finite element equation in displacement:
ξ 2 E 0 ξ ξ u ̂ b ( ξ ) + ξ ( E 0 + E 1 T E 1 ) ξ u ̂ b ( ξ ) E 2 u ̂ b ( ξ ) + ξ 2 ω 2 M 0 u ̂ b ( ξ ) = 0 .
(A3)
Since there is no known closed-form solution for this type of ODE Eq. (A3), the approximation approach for polygonal super-elements is more sophisticated. A differential equation for the local stiffness matrix S ̂ b is derived and approximated by a continued-fraction iteration.32 Defining the local stiffness as a matrix expression for the inner forces, i.e.,
S ̂ b ( ξ , ω ) u ̂ b ( ξ ) = ( ξ E 0 ξ + E 1 T ) u ̂ b ( ξ ) ,
considerations regarding the scalability of the wave equation for linear elastic problems lead to the assumption that the stiffness matrix is a function of the single quantity χ = ( i ξ ω ) 2, i.e., S ̂ b ( ξ , ω ) = S ̂ b ( χ ). Substituting these assumptions into Eq. (A3) leads to
2 χ χ S ̂ b ( χ ) + ( S ̂ b ( χ ) E 1 ) E 0 1 ( S ̂ b ( χ ) E 1 T ) E 2 χ M 0 = 0 .
(A4)
Subsequently, Eq. (A4) is approximated by a continued-fraction iteration of order M,
S ̂ b ( χ ) = S ̂ ( 0 ) ( χ ) , S ̂ ( m ) ( χ ) = S 0 ( m ) + χ S 1 ( m ) χ 2 X ( m + 1 ) ( S ̂ ( m + 1 ) ( χ ) ) 1 ( X ( m + 1 ) ) T , S ̂ ( M ) ( χ ) = S 0 ( M ) + χ S 1 ( M ) ,
where S 0 ( m ) , S 1 ( m ) are the continued-fraction matrices of the mth term and X ( m ) is the mth pre-conditioner matrix. The details of the computation of the matrices S 0 ( m ) , S 1 ( m ) and X ( m ) can be found in Ref. 32. The matrices are derived either by algebraic Riccati or Lyapunov equations. Both types of algebraic equations can be solved by an eigenvalue problem.8,9 Generally, the order M of the iteration can be adjusted, but we followed the simple rule M = p, where p = 6 is the polynomial degree of the shape functions. It is worth noting that the displacement in the re-entrant corner, which leads to reduced convergence in many numerical methods, is well approximated by the discretization used here. If a re-entrant corner is located in the scaling center, optimal convergence is still observed in SBFEM.39 This property is the reason for the two polygonal super-elements around the notch (see Fig. 6).
1.
M. J.
Lowe
,
P.
Cawley
,
J.
Kao
, and
O.
Diligent
, “
The low frequency reflection characteristics of the fundamental antisymmetric Lamb wave a0 from a rectangular notch in a plate
,”
J. Acoust. Soc. Am.
112
(
6
),
2612
2622
(
2002
).
2.
M.
Lowe
and
O.
Diligent
, “
Low-frequency reflection characteristics of the s0 Lamb wave from a rectangular notch in a plate
,”
J. Acoust. Soc. Am.
111
(
1
),
64
74
(
2002
).
3.
A.
Demma
,
P.
Cawley
,
M.
Lowe
, and
A.
Roosenbrand
, “
The reflection of the fundamental torsional mode from cracks and notches in pipes
,”
J. Acoust. Soc. Am.
114
(
2
),
611
625
(
2003
).
4.
A.
Demma
,
P.
Cawley
,
M.
Lowe
,
A.
Roosenbrand
, and
B.
Pavlakovic
, “
The reflection of guided waves from notches in pipes: A guide for interpreting corrosion measurements
,”
NDT&E Int.
37
(
3
),
167
180
(
2004
).
5.
J.-B.
Ihn
and
F.-K.
Chang
, “
Pitch-catch active sensing methods in structural health monitoring for aircraft structures
,”
Struct. Health Monit.
7
(
1
),
5
19
(
2008
).
6.
X.
Zhao
,
R. L.
Royer
,
S. E.
Owens
, and
J. L.
Rose
, “
Ultrasonic Lamb wave tomography in structural health monitoring
,”
Smart Mater. Struct.
20
(
10
),
105002
(
2011
).
7.
R.
Loendersloot
and
M.
Moix-Bonet
, “
Damage identification in composite panels using guided waves
,” in
Proceedings of the 5th CEAS Air & Space Conference
, Delft, The Netherlands (
2015
), p.
14
, available at https://research.utwente.nl/en/publications/damage-identification-in-composite-panels-using-guided-waves.
8.
J.
Bulling
,
B.
Jurgelucks
,
J.
Prager
, and
A.
Walther
, “
Defect reconstruction in a two-dimensional semi-analytical waveguide model via derivative-based optimization
,”
J. Acoust. Soc. Am.
152
(
2
),
1217
1229
(
2022
).
9.
C.
Song
and
J. P.
Wolf
, “
The scaled boundary finite-element method—alias consistent infinitesimal finite-element cell method—for elastodynamics
,”
Comput. Methods Appl. Mech. Eng.
147
(
3-4
),
329
355
(
1997
).
10.
J. P.
Wolf
and
C.
Song
, “
The scaled boundary finite-element method—a primer: Derivations
,”
Comput. Struct.
78
(
1-3
),
191
210
(
2000
).
11.
C.
Song
and
J. P.
Wolf
, “
The scaled boundary finite-element method—a primer: Solution procedures
,”
Comput. Struct.
78
(
1-3
),
211
225
(
2000
).
12.
A.
Griewank
and
A.
Walther
,
Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation
(
SIAM
,
Philadelphia
,
2008
).
13.
U.
Naumann
,
The Art of Differentiating Computer Programs. An Introduction to Algorithmic Differentiation
(
SIAM
,
Philadelphia
,
2011
).
14.
X.
Wang
,
W. T.
Peter
, and
A.
Dordjevich
, “
Evaluation of pipeline defect's characteristic axial length via model-based parameter estimation in ultrasonic guided wave-based inspection
,”
Meas. Sci. Technol.
22
(
2
),
025701
(
2010
).
15.
B.
Brence
,
Y.
Lugovtsova
,
J.
Bulling
,
P.
Kraemer
, and
J.
Prager
, “
Determination of a notch depth using ultrasonic guided waves
,” in
Proceedings of the 10th European Workshop on Structural Health Monitoring—EWSHM 2022
, Palermo, Italy (
Springer
,
New York
,
2022
).
16.
S.
He
and
C.-T.
Ng
, “
Guided wave-based identification of multiple cracks in beams using a Bayesian approach
,”
Mech. Syst. Signal Process.
84
,
324
345
(
2017
).
17.
S.
Ha
and
F.-K.
Chang
, “
Adhesive interface layer effects in PZT-induced Lamb wave propagation
,”
Smart Mater. Struct.
19
(
2
),
025006
(
2010
).
18.
A. N.
Shpak
,
M. V.
Golub
,
I.
Mueller
, and
C.-P.
Fritzen
, “
Experimental and numerical aspects of lamb waves excitation and sensing by rectangular piezoelectric transducers
,” in
European Workshop on Structural Health Monitoring
(
Springer
,
New York
,
2020
), pp.
505
514
.
19.
F.
Lanza di Scalea
and
S.
Salamone
, “
Temperature effects in ultrasonic Lamb wave structural health monitoring systems
,”
J. Acoust. Soc. Am.
124
(
1
),
161
174
(
2008
).
20.
P.
Aryan
,
A.
Kotousov
,
C.-T.
Ng
, and
B.
Cazzolato
, “
A model-based method for damage detection with guided waves
,”
Struct. Control Health Monit.
24
(
3
),
e1884
(
2017
).
21.
A.
Fichtner
,
Full Seismic Waveform Modelling and Inversion
(
Springer Science & Business Media
,
Berlin
,
2010
).
22.
J.
Rao
,
M.
Ratassepp
, and
Z.
Fan
, “
Guided wave tomography based on full waveform inversion
,”
IEEE Trans. Ultrason, Ferroelect, Freq. Contr.
63
(
5
),
737
745
(
2016
).
23.
J.
Rao
,
M.
Ratassepp
, and
Z.
Fan
, “
Investigation of the reconstruction accuracy of guided wave tomography using full waveform inversion
,”
J. Sound Vib.
400
,
317
328
(
2017
).
24.
J.
Rao
,
J.
Yang
,
M.
Ratassepp
, and
Z.
Fan
, “
Multi-parameter reconstruction of velocity and density using ultrasonic tomography based on full waveform inversion
,”
Ultrasonics
101
,
106004
(
2020
).
25.
P.
Huthwaite
and
F.
Simonetti
, “
High-resolution guided wave tomography
,”
Wave Motion
50
(
5
),
979
993
(
2013
).
26.
Z.
Zeng
,
M.
Gao
,
C. T.
Ng
, and
A. H.
Sheikh
, “
Guided wave-based characterisation of cracks in pipes utilising approximate Bayesian computation
,”
Thin-Walled Struct.
192
,
111138
(
2023
).
27.
H.
Gravenkamp
,
C.
Song
, and
J.
Prager
, “
A numerical approach for the computation of dispersion relations for plate structures using the scaled boundary finite element method
,”
J. Sound Vib.
331
(
11
),
2543
2557
(
2012
).
28.
J.
Bulling
,
G.
Franosch
,
Y.
Lugovtsova
, and
J.
Prager
, “
Sensitivity of ultrasonic guided waves to elastic constants: A numerical study
,” in
European Workshop on Structural Health Monitoring
(
Springer
,
New York
,
2020
), pp.
759
768
.
29.
O.
Gauthier
,
J.
Virieux
, and
A.
Tarantola
, “
Two-dimensional nonlinear inversion of seismic waveforms: Numerical results
,”
Geophysics
51
(
7
),
1387
1403
(
1986
).
30.
L.
Métivier
,
R.
Brossier
,
Q.
Merigot
, and
E.
Oudet
, “
A graph space optimal transport distance as a generalization of Lp distances: Application to a seismic imaging inverse problem
,”
Inverse Probl.
35
(
8
),
085001
(
2019
).
31.
C.
Boehm
,
L.
Krischer
,
I.
Ulrich
,
P.
Marty
,
M.
Afanasiev
, and
A.
Fichtner
, “
Using optimal transport to mitigate cycle-skipping in ultrasound computed tomography
,” in
Medical Imaging 2022: Ultrasonic Imaging and Tomography
(
SPIE
,
Bellingham, WA
,
2022
), Vol. 12038, pp.
48
57
.
32.
D.
Chen
,
C.
Birk
,
C.
Song
, and
C.
Du
, “
A high-order approach for modelling transient wave propagation problems using the scaled boundary finite element method
,”
Int. J. Numer. Methods Eng.
97
(
13
),
937
959
(
2014
).
33.
H.
Gravenkamp
,
C.
Birk
, and
C.
Song
, “
Simulation of elastic guided waves interacting with defects in arbitrarily long structures using the scaled boundary finite element method
,”
J. Comput. Phys.
295
,
438
455
(
2015
).
34.
H.
Gravenkamp
, “
Efficient simulation of elastic guided waves interacting with notches, adhesive joints, delaminations and inclined edges in plate structures
,”
Ultrasonics
82
,
101
113
(
2018
).
35.
E.
Kausel
,
Advanced Structural Dynamics
(
Cambridge University Press
,
Cambridge, UK
,
2017
).
36.
A. B.
Bakushinsky
,
M. Y.
Kokurin
, and
A.
Smirnova
,
Iterative Methods for Ill-Posed Problems: An Introduction
(
De Gruyter
,
Berlin
,
2011
).
37.
C. H.
Bischof
,
H. M.
Bucker
,
B.
Lang
,
A.
Rasch
, and
A.
Vehreschild
, “
Combining source transformation and operator overloading techniques to compute derivatives for matlab programs
,” in
Proceedings of the Second IEEE International Workshop on Source Code Analysis and Manipulation
, Montreal, QC, Canada (
IEEE
,
New York
,
2002
), pp.
65
72
.
38.
H.
Gravenkamp
,
A. A.
Saputra
, and
S.
Duczek
, “
High-order shape functions in the scaled boundary finite element method revisited
,”
Arch. Comput. Methods Eng.
28
,
473
494
(
2021
).
39.
J.
Bulling
and
H.
Gravenkamp
, “
Comparison of different models for stress singularities in higher order finite element methods for elastic waves
,”
Proc. Appl. Math. Mech.
19
(
1
),
e201900095
(
2019
).
Published open access through an agreement with Bundesanstalt für Materialforschung und -pruefung