Deterministic and stochastic damage detection via dynamic response analysis

TL;DR

This paper introduces a method for damage detection in elastic materials using dynamic response analysis, combining deterministic and Bayesian stochastic calibration of wave parameters to identify damage.

Contribution

It presents a novel approach that integrates deterministic and Bayesian stochastic methods for calibrating acoustic wave parameters in damage detection.

Findings

Successful calibration of wave speed and damping coefficient

Bayesian posterior distributions enable damage testing

Experimental validation supports method effectiveness

Abstract

The paper proposes a method of damage detection in elastic materials, which is based on analyzing the time-dependent (dynamic) response of the material excited by an acoustic signal. A case study is presented consisting of experimental measurements and their mathematical analysis. The decisive parameters (wave speed and damping coefficient) of a mathematical model of the acoustic wave are calibrated by comparing the measurement data with the numerically evaluated exact solution predicted by the mathematical model. The calibration is done both deterministically by minimizing the square error over time and stochastically by a Bayesian approach, implemented through the Metropolis-Hastings algorithm. The resulting posterior distribution of the parameters can be used to construct a Bayesian test for damage.

Equations

$${\partial }_{tt}u(z,t)-c^2{\partial }_{zz}u(z,t)=0$$

$$c^2=(\lambda +2\mu )/\rho ,$$

$$u(z,t)=\{\begin{array}{ll}v(z,t)& \text{ if }t\in [0,T_{\text{ex}}]\\ w(z,t)& \text{ if }t\in (T_{\text{ex}},T_{\text{end}}],\end{array}$$

$$\{\begin{array}{rlrl}{\partial }_{tt}v(z,t)+b{\partial }_{t}v(z,t)-c^2{\partial }_{zz}v(z,t)& =0& & \\ v(z,0)& =0& {\partial }_{t}v(z,0)& =0\\ v(0,t)& =f(t)& {\partial }_{z}v(L,t)& =0\end{array}$$

$$\{\begin{array}{rlrl}{\partial }_{tt}w(z,t)+b{\partial }_{t}w(z,t)-c^2{\partial }_{zz}w(z,t)& =0& & \\ w(z,T_{\text{ex}})-v(z,T_{\text{ex}})& =0& {\partial }_{t}w(z,T_{\text{ex}})-{\partial }_{t}v(z,T_{\text{ex}})& =0\\ {\partial }_{z}w(0,t)& =0& {\partial }_{z}w(L,t)& =0.\end{array}$$

$$F(t)=\{\begin{array}{ll}f(t)& t\in (0,T_{\text{ex}})\\ 0& \text{else}.\end{array}$$

$$V(z,t)=\{\begin{array}{ll}v(z,t)& t\ge 0\\ 0& t<0.\end{array}$$

$$\displaystyle\left\{\begin{aligned} \partial_{tt}V(z,t)+b\,\partial_{t}V(z,t)-c^{2}\,\partial_{zz}V(z,t)&=0\\ V(0,t)&=F(t)&\partial_{z}V(L,t)&=0\\ V(z,\cdot)\big{|}_{t<0}&\equiv 0.\end{aligned}\right.$$

$$\{\begin{array}{rl}-{\tau }^2\tilde{V}(z,\tau )+i\tau b\tilde{V}(z,\tau )-c^2{\partial }_{zz}\tilde{V}(z,\tau )& =0\\ \tilde{V}(0,\tau )=\tilde{f}(\tau )({\partial }_z\tilde{V})(L,\tau )& =0.\end{array}$$

$$\tilde{V}(z,\tau )=C_1(\tau ){\mathrm{e}}^{-B(\tau )z}+C_2(\tau ){\mathrm{e}}^{B(\tau )z},$$

$$C_1(\tau )=\tilde{F}(\tau )\frac{{\mathrm{e}}^{B(\tau )L}}{{\mathrm{e}}^{B(\tau )L}+{\mathrm{e}}^{-B(\tau )L}}\text{and}{C}_2(\tau )=\tilde{F}(\tau )\frac{{\mathrm{e}}^{-B(\tau )L}}{{\mathrm{e}}^{B(\tau )L}+{\mathrm{e}}^{-B(\tau )L}}$$

$$\tilde{V}(z,\tau )=\tilde{F}(\tau )\frac{{\mathrm{e}}^{(L-z)B(\tau )}+{\mathrm{e}}^{-(L-z)B(\tau )}}{{\mathrm{e}}^{LB(\tau )}+{\mathrm{e}}^{-LB(\tau )}}.$$

$$V(z,t)={\int }_{-\infty }^{\infty }{\mathrm{e}}^{i\tau t}\tilde{V}(z,\tau )\mathrm{d}to0.0pt-\tau$$

$$\frac{{\mathrm{e}}^{(L-z)B(\tau )}+{\mathrm{e}}^{-(L-z)B(\tau )}}{{\mathrm{e}}^{LB(\tau )}+{\mathrm{e}}^{-LB(\tau )}}$$

$$w(z,T_{\text{ex}})={\int }_{-\infty }^{\infty }{\mathrm{e}}^{i\tau T_{\text{ex}}}\tilde{V}(z,\tau )\mathrm{d}to0.0pt-\tau \text{and}({\partial }_{t}w)(z,T_{\text{ex}})={\int }_{-\infty }^{\infty }{\mathrm{e}}^{i\tau T_{\text{ex}}}i\tau \tilde{V}(z,\tau )\mathrm{d}to0.0pt-\tau .$$

$$W(z,t)=\{\begin{array}{ll}w(z,t)& z\ge 0\\ w(-z,t)& z<0\end{array}$$

$$V_e(z,T_{\text{ex}})=\{\begin{array}{ll}V(z,T_{\text{ex}})& z\ge 0\\ V(-z,T_{\text{ex}})& z<0.\end{array}$$

$$\{\begin{array}{rl}{\partial }_{tt}W(z,t)+b{\partial }_{t}W(z,t)-c^2{\partial }_{zz}W(z,t)& =0\\ W(z,T_{\text{ex}})-V_e(z,T_{\text{ex}})& =0{\partial }_{t}W(z,T_{\text{ex}})-{\partial }_t{V}_e(z,T_{\text{ex}})=0\\ {\partial }_{z}W(-L,t)& =0{\partial }_{z}W(L,t)=0.\end{array}$$

$$W_e(z,t)=\{\begin{array}{ll}⋮& ⋮\\ W(z-2L,t)& \text{if }z\in [-3L,-L)\\ W(z,t)& \text{if }z\in [-L,L)\\ W(z-2L,t)& \text{if }z\in [L,3L).\\ ⋮& ⋮\end{array}$$

$$\{\begin{array}{rl}{\partial }_{tt}W(z,t)+b{\partial }_{t}W(z,t)-c^2{\partial }_{zz}W(z,t)& =0\\ W(z,T_{\text{ex}})-V_e(z,T_{\text{ex}})& =0{\partial }_{t}W(z,T_{\text{ex}})-{\partial }_t{V}_e(z,T_{\text{ex}})=0\\ W(-L,t)-W(L,t)& =0{\partial }_{z}W(-L,t)-{\partial }_{z}W(L,t)=0.\end{array}$$

$$a_{k}(t)=A_{k}\ \exp\bigg{(}-\frac{b}{2}t+i\sqrt{\frac{c^{2}k^{2}\pi^{2}}{L^{2}}-\frac{b^{2}}{4}}t\bigg{)}+B_{k}\ \exp\bigg{(}-\frac{b}{2}t-i\sqrt{\frac{c^{2}k^{2}\pi^{2}}{L^{2}}-\frac{b^{2}}{4}}t\bigg{)}.$$

$$\arg \underset{b,c}{\min }{\Vert g_{\text{meas}}-g_{\text{comp}}\Vert }_{L_2([0,T_{\text{end}}])}.$$

$$\mathbb{P}(\mathbf{\theta }\in B|\mathbf{x}=\mathbf{\alpha })=\frac{{\int }_B{f}_{\mathbf{x}|\mathbf{\theta }=\mathbf{\beta }}(\mathbf{\alpha })f_{\mathbf{\theta }}(\mathbf{\beta })\mathrm{d}\mathbf{\beta }}{{\int }_{\Theta }{f}_{\mathbf{x}|\mathbf{\theta }=\mathbf{\beta }}(\mathbf{\alpha })f_{\mathbf{\theta }}(\mathbf{\beta })\mathrm{d}\mathbf{\beta }},$$

$$f_{\mathbf{\theta }|\mathbf{x}=\mathbf{\alpha }}(\mathbf{\beta })=\frac{f_{\mathbf{x}|\mathbf{\theta }=\mathbf{\beta }}(\mathbf{\alpha })f_{\mathbf{\theta }}(\mathbf{\beta })}{f_{\mathbf{x}}(\mathbf{\alpha })}.$$

$$f_{\mathbf{\theta }|\mathbf{x}=\mathbf{\alpha }}(\mathbf{\beta })=\kappa f_{\mathbf{x}|\mathbf{\theta }=\mathbf{\beta }}(\mathbf{\alpha })f_{\mathbf{\theta }}(\mathbf{\beta }),$$

$$\pi ({\mathbf{\eta }}^k,{\mathbf{\beta }}^{k-1})=\{\begin{array}{ll}\min \left\{1,\frac{p({\mathbf{\eta }}^k)q({\mathbf{\beta }}^{k-1},{\mathbf{\eta }}^k)}{p({\mathbf{\beta }}^{k-1})q({\mathbf{\eta }}^k,{\mathbf{\beta }}^{k-1})}\right\}& \text{if }p({\mathbf{\beta }}^{k-1})q({\mathbf{\beta }}^{k-1},{\mathbf{\eta }}^k)>0\\ 1& \text{otherwise}\end{array}$$

$$\sup\left\{\Big{|}{\boldsymbol{P}}^{m}(A|{\boldsymbol{\beta}}^{0})-\int_{A}{p({\boldsymbol{\beta}})}\,\mathrm{d}{\boldsymbol{\beta}}\Big{|},{A\subset\Theta\text{ measurable}},{\boldsymbol{\beta}}^{0}\in\Theta\right\}\leq Mr^{m},$$

$${\int }_{\Theta }(1+\left|f(\mathbf{\beta })\right|{)}^{2}p(\mathbf{\beta })\mathrm{d}\mathbf{\beta }<\infty$$

$$\sqrt{N}\left(\frac{1}{N}\sum _{m=1}^{N}f({\mathbf{\beta }}^m)-{\int }_{\Theta }f(\mathbf{\beta })p(\mathbf{\beta })\mathrm{d}\mathbf{\beta }\right)\stackrel{d}{\to }\mathcal{N}(0,{\sigma }_f^2),$$

$$b_{\text{min}}$$

$$c_{\text{min}}$$

$$f_{\mathbf{x}|\mathbf{\theta }=\mathbf{\beta }}(\mathbf{\alpha })=C\exp \left(-\frac{(\mathcal{M}(\mathbf{\beta })-\mathbf{\alpha }{)}^T{\Sigma }^{-1}(\mathcal{M}(\mathbf{\beta })-\mathbf{\alpha })}{2}\right),$$

Full text

1 Introduction

Using elastic waves in solids for identifying material properties has a long history going back to the 19th century with the pioneering work of Rayleigh, Lamb and Love [1]. In the second half of the 20th century, wave propagation analysis became a central field in seismology [2] and in the material sciences [3]. Especially the analysis of ultrasonic waves, as a part of vibration based condition monitoring [4], has become a prominent field in structural health monitoring.

Specifically with respect to structural health monitoring, changes in material properties, such as damage or fatigue of a structure, are sought to being detected by means of the response of the structure to an excitation by an acoustic wave. The focus of this paper will be on acoustic waves induced by a piezoelectric transducer, and measured at the same position or other sensor locations.

In the literature, a wealth of methods have been developed, among them: (a) modal analysis [3]; (b) frequency based methods [5, 6]; (c) general spectral methods [7, 8]; (d) comparison with finite element solutions [9, 10, 11, 12]; (d) Rayleigh, Lamb and Love waves [12]; (e) nonlinear analyses [13]; (f) guided waves [14]; (g) wavelet analysis [15]; (h) time series methods [16, 17]; (k); statistical indicators [18]; artificial intelligence [4, 19], just to cite a few representative papers. For further information, we refer to the monographs [7, 20, 21, 22], the classical survey [23] and the more recent survey papers [4, 8, 24, 25, 26, 27].

Ultrasonic imaging can be handled in one-, two- and three-dimensional elasticity. Specific lower dimensional applications are to damaged rods [11, 28, 29] as well as two-dimensional composites [14, 30], including damage location using sensor configurations [31]. In industrial applications, it is common practice to measure the transit time of the reflected pressure wave, induced by a piezoelectric transducer, and to move the transducer from point to point along the structure. This amounts to measuring the propagation speed of an approximately one-dimensional signal.

Especially in seismology and geotechnics [32, 33], the random structure of the soil has been incorporated in the study of wave propagation in more recent years. The same applies to material science, see [34] and references therein.

The present paper is part of a research project addressing linear wave propagation in random media. The main features of our approach are: (a) models are set up in a way so that exact solutions can be used and (b) stochastic analyses are employed to obtain indicators for changes in material parameters and for damage detection. The program encompasses unidirectional propagation (transport equations), one-dimensional acoustic waves and wave propagation in three-dimensional linearly elastic solids. In the present paper, we present a theoretical and experimental analysis of ultrasonic inspection of composite plates. The acoustic wave produced by a piezoelectric transducer is modelled as a one-dimensional wave issuing from the (moving) location of the transducer, which also acts as recording sensor for the reflected wave. The corresponding initial-boundary value problem for the damped wave equation is solved explicitly, using Fourier transform techniques. The propagation speed and damping coefficient are directly related to elastic parameters of the material and can be calibrated by comparison with the time-dependent signal. In addition, a stochastic, Bayesian parameter estimation allows one to design hypothesis tests for critical thresholds of the model parameters. A three-dimensional analysis by means of Fourier integral operators will be presented in a forthcoming paper.

In the specific case of ultrasonic pulse echo measurements, the presented work includes an experimental part of performing the measurements, and a theoretical part of setting up the wave propagation model, calibrating the parameters, finding the posterior distributions of the calibrated parameters, and setting up Bayesian hypothesis tests for damage.

The experimental set-up included four carbon fiber composite plates, three of which were damaged by a localized impact. The plates were measured with a piezoelectric transducer in impulse echo mode. As mentioned, the wave propagation through the material at a single measuring location is modeled by means of a one-dimensional wave equation with a damping term included. This plane wave assumption is justified through the respective dimensions of the transducer head and the thickness of the plates. Having set up the wave propagation model, including the incoming and reflecting boundary conditions, the two parameters (wave speed and damping coefficient) are calibrated by fitting the numerically evaluated exact solution to the measurements.

The calibration is first performed deterministically at each location by minimizing the mean square error over the observation time. Repeating this procedure on a grid of points on the plate allows one to identify locations in which the parameters deviate from the overall mean values.

However, a more stringent approach to determining the statistical properties of the parameter values at single locations is stochastic parameter calibration by means of Bayesian methods. Indeed, based on a priori bounds on the parameters and the likelihood function (given by the probability distribution of the error between certain measured and computed features of the response), the posterior (joint) distribution of the wave speed and the damping coefficient can be obtained at each location. What is more, the posterior density admits to determine credible regions, by means of which Bayesian hypothesis tests can be designed. One may rephrase the null hypothesis of undamaged material as a non-critical domain for the model parameters; the posterior probability of the null hypothesis is nothing but the posterior probability of the non-critical domain. The so designed Bayesian test, performed at a 1% rejection threshold, enables the location of damaged points on the plate. The results are in accordance with the deterministic approach, but the Bayesian test contains more statistical information.

The outline of this paper is as follows. In the second section the measurement set-up will be described in more detail.

The third section is devoted to setting up the mathematical model of the waves traveling through the plate. We construct a Fourier transform based solution operator and we show that the model is capable of reproducing the measured signals.

In the fourth section we describe methods for parameter estimation and testing. In the first part we deterministically calibrate the wave speed and the damping coefficient to get a best fit with the measured signal at each location of the plate, using minimization of the mean square error between measured signal and computed solution over time. In this way, damage can be localized by observing deviations of the parameters from their nominal values. In the second part we apply Bayesian methods to compute the posterior distribution of the parameters. To this end, we implement the Metropolis-Hastings algorithm, generating a Markov chain the stationary distribution of which is the posterior distribution. Thereby, the uncertainty of the parameter estimation for each single location on the plate can be quantified. Finally, in the third part, we construct a Bayesian test for damage, which allows one to identify damaged areas on the plate. The paper concludes with a discussion and a summary.

The paper is based on the work [35] of the second author. Partial results on transport equations and on acoustic waves in three-dimensional elasticity have been reported in [36, 37, 38, 39].

2 Ultrasonic impulse echo

As mentioned in the introduction, the measurement data were obtained by recording the response of four carbon fiber composite plates to an ultrasonic signal. Three of the specimen had been damaged by a high speed impact and the fourth plate was undamaged. The dimensions of the specimen were 105 mm times 95 mm with a thickness of about 13 mm. (The exact thickness is not relevant for the analysis.) All four plates were scanned with a horizontal resolution of $5\times 5$ millimeters. In addition, the damaged area was scanned with a horizontal resolution of $1\times 1$ mm. The damaged specimen analyzed further in the paper is shown in Figure 1(a).

The principle of an ultrasonic impulse echo measurement is the following: A piezoelectric transducer produces an ultrasonic pulse at the surface of the plate which then goes through the plate as an elastic wave. The wave is reflected at the bottom and comes back to the top. The transducer measures the amplitude at the top over time, which is recorded as voltage by means of an oscilloscope. In order to have good contact between transducer and plate a water film was placed on top of the plate. As pulse generating and recording device an Olympus${}^{\mbox{\footnotesize\textregistered}}$ Epoch 600 Ultrasonic Flaw Detector was used. The pulser voltage was set to 400 V; a transducer of diameter 22 mm and resonance frequency of 1 MHz was employed. The oscilloscope had a resolution of $400$ MHz in time. For each scan, the recorded period was $35$ $\mathrm{\mu}$s and thus $14000$ data points per scan were acquired. The amplitude resolution of the oscilloscope was $512$ points, where number $256$ represented the zero line. For the subsequent mathematical analysis, the signal was normalized by its maximal value of $256$ points. Accordingly, the analyzed signal was dimensionless with a range between $-1$ and $1$. In accordance with widespread usage (e.g., [12, Chapter 10], [14]), [22, Chapter 5]), we shall refer to it as normalized amplitude.

The transmission from the oscilloscope to the PC sometimes produced scrambled signals. To detect such measurement errors, the transmitted signal was analyzed by an automated script. A faulty signal typically had a jump at a certain time point. Thus, the successive differences in time were considered. During the excitation period larger differences were accepted, whereas in the arriving echoes the differences had to be small, since the signal should be continuous. Measurements of signals classified as erroneous were repeated.

Furthermore, the signal of the transducer not contacting the plate was subtracted from the signal touching the plate, since the transducer head measured the vibrations within itself, too. The resulting adjusted signal depicts the actual vibrations in the plate (see Figure 2).

As a final remark we would like to point out that the excitation is done by sending an electrical rectangular pulse to the transducer head. This pulse causes the transducer head to generate a sinusoidal vibration. However, this pulse is stronger than the oscilloscope can measure, which produced an overflow in amplitude direction for the first 4 $\mu$s. This means that the initial pulse cannot be measured exactly and by subtracting the signals it is set to zero (although it is nonzero). As a consequence, the first period of the initial pulse is measured as zero. Therefore, the echoes have one more period than the initial pulse in the measurements, as can be seen in Figure 3 (located on the time axis between $12$ $\mu s$ and $13$ $\mu s$).

3 Mathematical modeling: 1D telegraph equation

To set up a simple mathematical model that nevertheless captures the features of interest, we assume that the undamaged plate is a linearly elastic, homogeneous, isotropic continuum of constant density $\rho$. The in-plane coordinates will be denoted by $x$ and $y$, the vertical direction by $z$ with $z=0$ corresponding to the upper surface of the plate. Further, we assume that the transducer induces a plane wave moving in $z$-direction of the plate. This can be justified, since the diameter of the transducer is larger than the thickness of the plate (approximately $2:1$). Plane waves can be reduced to the one-dimensional acoustic wave equation

[TABLE]

for the displacement $u(z,t)$ in $z$-direction at time $t$. The wave speed is related to the Lamé constants $\lambda$, $\mu$ of the material through

[TABLE]

see [1, Sections 1.2 and 1.3]. The acoustic wave equation neglects attenuation which is due to absorption and to scattering of waves at inhomogeneities of the material [21, Chapter 6]. The attenuation can also be observed in the experimental data, see Figure 3(a). Thus it is justified to add a damping term to the wave equation, resulting in the telegraph equation. The bottom side of the plate is assumed to be stress free, which implies that the $z$-derivative of $u$ vanishes.

During the time interval $[0,T_{\text{ex}}]$ the excitation due to the piezoelectric transducer produces a displacement $f(t)$ at the top side of the plate [20]. We assume $f$ to be smooth and to be of compact support in $(0,T_{\text{ex}})$. During the observation period $(T_{\text{ex}},T_{\text{end}}]$, a stress free boundary on top side is assumed, and the time-dependent displacement induces the signal recorded by the transducer.

As noted in Section 2, the exact value of the plate thickness is irrelevant to the analysis. (One-sided ultrasonic measurements do not allow one to measure the absolute values of the plate thickness and the wave speed simultaneously.) We assume from now on that the thickness of the plate has the constant value $L$, which also serves as our (virtual) length unit.

This results in the following equations, valid for $0\leq z\leq L$:

[TABLE]

In this set-up, $f(t)$ for $t<T_{\text{ex}}$ and $w(0,t)$ for $t>T_{\text{ex}}$ correspond to the readings of the oscilloscope and couple the measurements with the solution to (2a)–(2c). The output of the oscilloscope is voltage. As mentioned in the introduction, the signal is normalized by its maximally measurable value. We refer to this dimensionless quantity as normalized amplitude.

The graphs in Figure 3 confirm the good coherence between the measured signal and the solution of (2a)–(2c) at $z=0$, the top of the plate. The parameters $b$ and $c$ were calibrated such that the mean square error of the computed and the measured signal was minimal. The calibration procedure will be explained in detail in the next section. As mentioned before, the measured echoes have one more period in the time between $12$ $\mu s$ and $13$ $\mu s$ than the computed ones. This is due to the fact that the force term cannot be measured in the first 4 $\mu s$ and is set to zero.

Various analytical methods for solving equations (2b) and (2c) are available (see e.g. [40]). It turned out that for numerical reasons, problem (2b) is advantageously solved by applying Fourier transform in time direction, which can be evaluated very quickly. Problem (2c) is more easily solved by Fourier series expansion.

We start with solving (2b). In order to apply the Fourier transform in the time direction, we first have to extend the time domain to the full time axis. We define

[TABLE]

At this stage, we let $v$ be the solution of (2b) for $t\in[0,\infty)$ (not only $t\in[0,T_{\text{ex}}]$). We further extend $v$ to negative times by setting

[TABLE]

If $V$ is at least twice continuously differentiable, it satisfies

[TABLE]

We formally deduce the solution: If we apply the Fourier transform in time direction (denoted as ${\mathcal{F}}_{t\to\tau}[V]=\widetilde{V}$), we get the ordinary differential equation

[TABLE]

This is solved by

[TABLE]

where $B(\tau)=\frac{1}{c}\sqrt{-\tau^{2}+ib\tau}$, and $\sqrt{\cdot}$ is the principle branch of the complex root, and $\sqrt{0}=0$. The constants $C_{1}$ and $C_{2}$ are given by

[TABLE]

and consequently

[TABLE]

Applying the inverse Fourier transform yields the formal solution

[TABLE]

where \,\mathrm{d}\hbox to0.0pt{\hss\raisebox{3.65973pt}{\scriptstyle-!}}\tau is shorthand for $\,\mathrm{d}\tau/2\pi$.

Proposition. Let $V$ be as above. Then (4) is a convergent integral and $V$ is infinitely differentiable and satisfies (3).

Proof.

We will only sketch how to prove this. The steps are as follows:

• •

Since $F$ is a Schwartz function (infinitely differentiable, faster decay than any negative power of $\left|t\right|$), the time Fourier transform of $F$ is a Schwartz function, too. Thus, it suffices to show that

[TABLE]

is bounded in $\tau$. Then (4) is a convergent integral.

• •

To show that $V$ is smooth one needs to show that differentiation under the integral sign is justified.

• •

It is trivial to show that (4) satisfies the boundary condition. Finally, we use the Paley-Wiener theorem, which gives a sufficient condition on $\widetilde{V}$ such that V(z,\cdot)\big{|}_{t<0}\equiv 0.

A detailed proof of all steps can be found in [35]. ∎

We can compute initial values of (2c) by

[TABLE]

Both the Fourier transform and the inverse Fourier transform can be numerically approximated by the fft respective ifft algorithm. Since the force term was smooth enough, the transformed signal decayed fast enough to get a good numerical approximation.

For the numerical calculation we set $L=1$. (As noted at the beginning of Section 3, this is no restriction.) As spatial discretization we chose $\Delta z=0.001$. The time discretization was the same as from the oscilloscope, i.e. $T_{\text{max}}=35\;\mu$s and $\Delta t=0.0025\;\mu$s. Furthermore, we set $T_{\text{ex}}=11.8\;\mu$s. This is the time after which the force term $f$ was zero in all measurements.

In order to have a smooth input signal $f$, the measured signal of the transducer was regularized by smoothing out the high frequencies. For this purpose we multiplied $\widetilde{f}$ with a Tukey window (see e.g. [41]), removing frequencies beyond 6.5 MHz.

We use the reflection principle (see e.g. [40, Chapter 2, Section 2]) to extend (2c) to the strip $[-L,L]\times[T_{\text{ex}},\infty)$ as

[TABLE]

and

[TABLE]

Then $W$ satisfies

[TABLE]

Since $W$ is even and continuously differentiable, the homogeneous Neumann boundary condition is equivalent to periodic boundary conditions. This can be shown by iteratively applying the reflection principle and extending $W$ to the whole space by

[TABLE]

Then ${W}_{e}$ is $2L$-periodic. On the other hand, if $W$ is even and periodic, then $W(z+L)=W(z-L)=W(-z+L)=W(-z-L)$, and thus $W_{z}(L)=W_{z}(-L)=0$.

Thus, we can equivalently solve the periodic boundary value problem

[TABLE]

We make the ansatz $W(z,t)=\sum_{k\in\mathbb{Z}}a_{k}(t){\mathrm{e}}^{i{k\pi}{}z/L}$, which is $2L$-periodic in space. Plugging it into the equation yields

[TABLE]

The constants $A_{k}$ and $B_{k}$ are determined from a Fourier series expansion of ${V}_{e}(z,T_{\text{ex}})$ and $\partial_{t}{V}_{e}(z,T_{\text{ex}})$.

For the numerical calculation we used the same grid size: $L=1,\Delta z=0.001.$ The time discretization was the same as from the oscilloscope, i.e. $T_{\text{max}}=35$ $\mu$s and $\Delta t=0.0025$ $\mu$s. The computed signals at an undamaged and a damaged location are shown in Figure 3. The evaluation of the numerical solution takes approximately 2.8 seconds on a PC.

4 Parameter estimation and tests

The mathematical model described above will serve as an input-output model. As input we have the measured excitation $f$ and the (unknown) material parameters $b$ and $c$. As output we get $g_{\text{comp}}(t)=u(0,t)$. Furthermore, we will call $g_{\text{meas}}(t)$ the measured signal.

In the following subsection we will calibrate the material parameters $b$ and $c$ at each measurement location. In the second subsection we will estimate the posterior distribution of $b$ and $c$ for every measurement location. According to our length convention, the dimension of $c$ is L/$\mu$s. The dimension of $b$ is 1/$\mu$s.

4.1 Parameter estimation via optimization

In order to estimate $b$ and $c$ in a given location, we solve the non-convex optimization problem

[TABLE]

The optimal parameters were computed using the Nelder-Mead algorithm. Because the problem is nonconvex we had to be careful to choose a good initial value. As initial value, we chose the pair $(b,c)=(0.2,0.22)$, which seemed to work well for most locations. Figure 3 shows the comparison between the measured data and the computed signal with optimized parameters at a single location. The spatial resolution of the optimized parameters $b$ and $c$ of the plate specimen under investigation is displayed in Figure 4.

We observe that the wave speed decreases and the damping coefficient increases in the damaged region. The damage affects the stiffness of the material, which makes waves travel more slowly. Indeed, in an elastic medium the Lamé constants and hence the square of the wave speed are proportional to the modulus of elasticity $E$, compare Equation (1) and [1, Table 2.1], and damage generally goes together with a decrease of $E$ [42, Section 2.1.2], [43]. On the other hand, one expects delamination within the damaged region, which leads to higher damping, that is, attenuation of the signal [44, 45].

4.2 Bayesian parameter estimation

The least square error optimization described above produces a point estimator for the parameters $b$ and $c$. The Bayesian inference described in the sequel provides more information about the estimated parameters. Using prior knowledge of the parameters this method gives a probability distribution of the parameters. This can be used to get an estimator for $b$ and $c$ but also to assess the variation of the parameters.

For this reason consider a state space $X\subset{\mathbb{R}}^{n}$ and a parameter space $\Theta\subset{\mathbb{R}}^{m}$. Let ${\boldsymbol{x}}$ resp. ${\boldsymbol{\theta}}$ be random variables with values in $X$ resp. $\Theta$ and marginal density $f_{{\boldsymbol{x}}}$ resp. $f_{{\boldsymbol{\theta}}}$. Then, by Bayes’ theorem

[TABLE]

where $B$ is a Borel measurable set and $f_{{{\boldsymbol{x}}}|{{\boldsymbol{\theta}}={\boldsymbol{\beta}}}}$ is the conditional density of ${\boldsymbol{x}}$ given that ${\boldsymbol{\theta}}={\boldsymbol{\beta}}$. In terms of densities one has

[TABLE]

We use this rule for Bayesian inference [46] with a fixed model: If ${\boldsymbol{\alpha}}\in X$ is the measured data obtained in an experiment, then the posterior probability density function of the parameters ${{\boldsymbol{\theta}}}\in\Theta$ is

[TABLE]

where $\kappa$ is a constant depending only on ${\boldsymbol{\alpha}}$. Here, $f_{{\boldsymbol{\theta}}}$ is the prior probability distribution of the parameters ${{\boldsymbol{\theta}}}$ and $f_{{{\boldsymbol{x}}}|{{\boldsymbol{\theta}}}}$ is the likelihood function.

We assume that $f_{{\boldsymbol{x}}|{\boldsymbol{\theta}}}$ is known. In principle, one could compute the posterior distribution directly. However, usually – and also in our case – the model is computationally expensive to evaluate, so a direct evaluation of the integrals in (5) is problematic. We will employ the Metropolis-Hastings algorithm (MHA) instead.

This algorithm allows one to generate a Markov chain with a (up to a normalization constant) given stationary distribution [47]. Some mild assumptions guarantee the convergence of the Markov chain to the stationary distribution. The generated sample can then be used to estimate the conditional probability $f_{{\boldsymbol{\theta}}|{\boldsymbol{x}}}$ by various methods, e.g. by applying a smooth kernel density estimator.

For the MHA one needs an initial probability density $p_{0}$ and a target probability density $p$. Furthermore, one needs a proposal probability density $q(\cdot,{\boldsymbol{y}})$, which may depend on parameter ${\boldsymbol{y}}\in\Theta$. Then, the Markov chain is generated as follows.

• •

Sample the first Markov chain link ${\boldsymbol{\beta}}_{0}$ according to the given initial probability distribution $p_{0}$.

• •

The $k$th chain link ${\boldsymbol{\beta}}^{k}$ is generated as follows: Generate a candidate ${\boldsymbol{\eta}}^{k}$ according to the proposal distribution $q(\cdot,{\boldsymbol{\beta}}^{k-1})$. Then compute

[TABLE]

and set ${\boldsymbol{\beta}}^{k}={\boldsymbol{\eta}}^{k}$ with probability $\pi(\boldsymbol{\eta}^{k},{\boldsymbol{\beta}}^{k-1})$. Otherwise set ${\boldsymbol{\beta}}^{k}={\boldsymbol{\beta}}^{k-1}$.

Finally, one arrives at the Markov chain $({\boldsymbol{\beta}}^{0},{\boldsymbol{\beta}}^{1},\ldots,{\boldsymbol{\beta}}^{N})$ having asymptotically the distribution $p$.

Naturally, the question arises, what criteria guarantee the convergence to the stationary distribution and, secondly, how fast is the convergence.

For the convergence theory we refer to [35, 48, 49, 50, 51]. A possible answer can be summarized as follows:

• •

If the initial density $p_{0}$ equals the target density $p$, then ${\boldsymbol{\beta}}^{k}$ is $p$-distributed for $k\in\left\{0,\ldots,N\right\}$.

• •

Let $\Theta^{+}=\left\{{\boldsymbol{\beta}}\in\Theta:p({\boldsymbol{\beta}})>0\right\}$. If $\Theta^{+}$ has finite Lebesgue measure and $p({\boldsymbol{\beta}})$ and $q({\boldsymbol{\beta}},{\boldsymbol{y}})$ are bounded away from zero on $\Theta^{+}$, then there exists a constant $M$ and $r<1$ such that

[TABLE]

where ${\boldsymbol{P}}^{m}(\cdot|{\boldsymbol{\beta}}^{0})$ is the probability measure of the $m$th output of the Metropolis-Hastings algorithm with given ${\boldsymbol{\beta}}^{0}$.

Furthermore, for every function $f$ with

[TABLE]

there exists a constant $\sigma_{f}$, such that for $N\to\infty$

[TABLE]

where ${\mathcal{N}}(0,\sigma_{f}^{2})$ is the normal distribution with variance $\sigma_{f}^{2}$ and the arrow indicates convergence in distribution.

Concerning the second issue, one should point out that the efficiency strongly depends on the proposal distribution $q(\cdot,{\boldsymbol{y}})$. For instance if $q(\cdot,{\boldsymbol{y}})\approx p$ the MHA will work most efficiently. Of course, usually this is not the case. However, it is noted in [52] that the efficiency does not depend on the type of the proposal distribution, but much more on the spread, which usually can be controlled by the variance of the proposal distribution. If the variance is too small, the sample gets too correlated. But if the variance is too large, the acceptance rate is too low and one has only very few accepted samples. So the question may as well be, “what is the optimal acceptance rate for the MHA”.

One of the first papers which tries to give an answer to this question, is [53]. The authors show that for a standard normally distributed target distribution the optimal acceptance rate is approximately $44\%$ for dimension $d=1$ and $23\%$ for dimension $d\rightarrow\infty$. The concept has been generalized to a multidimensional random variable with independent components in [54]. In this case the number $23\%$ for very high dimensions is not always attained, depending on the target distribution. The authors state that it remains still unclear what is an optimal acceptance rate for a general target distribution.

There is an additional issue concerning the efficiency of the MHA: the burn-in phase. Like stated above, if the initial distribution is the target distribution, we have perfect sampling. However, if the initial distribution is far away from the target distribution, the first few output samples of the MHA are not distributed according to the target distribution, and thus one deletes the first few outputs of the algorithm. For the interested reader we refer to [55], where several methods of output analysis (such as variance ratio method, spectral method or cumulative sum method) are described. These methods can help to determine if the chain has converged to the stationary distribution.

In the case of the 1D telegraph equation, we are interested in the posterior distribution of the parameters ${\boldsymbol{\theta}}=(b,c)$. Since we do not have any prior knowledge about the parameters except reasonable bounds on their range, we assume the uniform distribution with bounds

[TABLE]

which were chosen based on numerical experiments.

The choice of the features ${\boldsymbol{\alpha}}$ is non-trivial. Unlike in the subsection before, one cannot take the square difference between the measured signal and some “mean” signal, since there is no such thing as a “mean” signal. If one takes the arithmetic mean of all measured signals, by the phase shift, this is likely to be close to zero. Thus, we chose a set of characteristic features of the signal. It turned out that the phase angle and the amplitude of the three most dominant (in terms of amplitude) frequencies of the first echo in the discrete Fourier spectrum are feasible as model features. These were computed as follows: We set the signal before and after the first echo to zero (i.e., between 11.8 $\mu$s and 22 $\mu$s) and compute the discrete Fourier transform. The frequencies are between $1.12$ and $1.19$ MHz (i.e. the $33$rd to $35$th entry). So, the features were ${\boldsymbol{\alpha}}_{\text{}}=(\varphi_{33},\varphi_{34},\varphi_{35},r_{33},r_{34},r_{35})$ which are the amplitudes and the phase angles. We assume to have a normally distributed error with zero mean and covariance matrix $\Sigma$. In [56, Chapter 2] it is suggested to choose $\Sigma$ from measured data. Since we had an undamaged plate at hand, we used the scans of that plate to compute $\Sigma$. In fact, the features $(\widetilde{{\boldsymbol{\alpha}}}_{1,\text{meas}},\ldots,\widetilde{{\boldsymbol{\alpha}}}_{n,\text{meas}})$ obtained at the $n$ grid points of the undamaged plate can be seen as realizations of the feature values at any single location (assuming homogeneity of the underlying random field) and hence can serve to estimate the covariance matrix $\Sigma$.

Although the support of the probability density of the normal distribution is the whole space (contrary to the features) the likelihood function is chosen as

[TABLE]

where $C$ is a constant and the function ${\mathcal{M}}$ is the numerical solution operator from the previous section, computing the features $(\varphi_{33},\varphi_{34},\varphi_{35},r_{33},r_{34},r_{35})$ for given parameters $b$ and $c$. The choice can be justified since the variances were very small and thus $f_{{\boldsymbol{x}}|{\boldsymbol{\theta}}={\boldsymbol{\beta}}}$ decreased fast enough and no further cut-offs were necessary.

So for the $i$th measurement location we have

[TABLE]

where $D$ is a constant.

The Metropolis-Hastings algorithm was performed with the following proposal distribution:

[TABLE]

where

[TABLE]

and

[TABLE]

Experiments with the model showed that the posterior distribution of $b$ and $c$ is concentrated in the center of the prior domain. So we set the initial guess at ${\boldsymbol{\beta}}^{0}=\frac{1}{2}({b_{\text{min}}+b_{\text{max}}},{c_{\text{min}}+c_{\text{max}}})$. For this setup the condition of an exponential convergence of the Markov chain to the target distribution is satisfied. Since ${\boldsymbol{\beta}}^{0}$ is not $f_{{\boldsymbol{\theta}}|{\boldsymbol{x}}}$-distributed, we neglect the first $100$ samples of the algorithm as we interpret this as burn-in phase. After the burn-in phase, the length of the Markov chain was $N=1000$. In Figure 5(a) one can see a typical example of the Markov chain at one certain measurement location. After the burn-in phase it settles in a certain region and stays there. Figure 5(b) shows the joint posterior distribution of $b$ and $c$ at this measurement location.

Figure 6 depicts the spatial resolution of the mean value of the posterior distribution in each scan location. One can observe that the image is not as clear as in Figure 4. This comes from the fact that the chosen features were only the phase angle and amplitude of the most dominant frequencies of the first echo. One may get different results with other features.

4.3 Bayesian damage test

We now formulate the following Bayesian hypothesis test (see e.g. [57, Chapter 8, §1]). As null hypothesis we assume that the material is undamaged. This means that the material parameters $(b,c)$ are in the domain

[TABLE]

defined by critical thresholds $b_{\text{crit}}$, $c_{\text{crit}}$. Then the posterior probability of the null hypothesis is just

[TABLE]

We chose a confidence level of $1\%$ and, thus, the null hypothesis is rejected if ${\mathbb{P}}({\boldsymbol{\theta}}\in\Theta_{0}|{\boldsymbol{x}}={\boldsymbol{\alpha}})<0.01$.

In real life applications, the thresholds $b_{\text{crit}}$, $c_{\text{crit}}$ will be given by engineering requirements on the material properties of the plates. Alternatively, they could be estimated from the response of an undamaged plate. In our case, the threshold $b_{\text{crit}}$ resp. $c_{\text{crit}}$ was chosen from the optimized parameters (cf. Subsection 4.1) such that $99\%$ parameters of the undamaged region of the same plate were below $b_{\text{crit}}$ resp. above $c_{\text{crit}}$. The pragmatic reason for this choice was that we could not use the estimated parameters from another, undamaged plate, since all plates were subjected to a grinding procedure which resulted in a thickness difference of $2-3\%$ between different plates. Due to our modelling assumption of constant plate thickness $L$, taking data from different plates would produce an error of the wave speed in the same range, and thus would affect the choice of $c_{\text{crit}}$. For that reason we used the undamaged region of the damaged plate as reference.

Figure 7 shows the spatial resolution of the posterior probability of the null hypothesis. It is seen that the null hypothesis is rejected in the area enclosed by the black square, which coincides with the actually damaged area of the plate under investigation (compare with Figure 1), and the null hypothesis is accepted in the undamaged parts. Thus the damage is correctly identified, at least within the accuracy admitted by the grid size.

5 Discussion

As outlined in the introduction, a wealth of methods for damage detection based on ultrasonic measurements have been developed and are in use. The approach presented in this paper is a refined analysis of data acquired by so-called A-scans [12, 21]. An A-scan is a unidirectional measurement, usually in the direction perpendicular to the surface of the material. A change in the transit time of the reflected wave indicates a damaged area. It is common industrial practice to check the production quality, e.g., of carbon fiber composite plates, by rapidly performed A-scans. The intention of the present paper is to extract additional information by analyzing the time dependent, dynamic signal as recorded by the oscilloscope. This is done by setting up a simple model for one-dimensional wave propagation and calibrating wave speed and damping coefficient to the recorded data. This way two independent indicators are obtained, as compared to the usual practice of monitoring the transit time only. An additional benefit of our approach is that our model can be solved analytically, thereby enabling a rapid computation of the exact solution. The low computational cost is a big advantage in the required optimization step or the evaluation of the posterior distribution, respectively. Using exact analytical solutions, when possible, has also been recommended in the literature [12], [25, Section 4].

Our analysis is solely based on the data recorded by the oscilloscope. It does not require knowledge of the absolute value of the plate thickness nor the value of the wave speed, so no further (mechanical or geometrical) measurements are needed (in the case of plates of constant thickness). On the other hand, it does not allow one to infer these quantities either. Another limitation is that our method does not admit the determination of the elastic parameters (the Lamé constants) of the material.

A possible application could be to automated scanning, where the decision about the presence of damage could be based on the hypothesis test presented in Subsection 4.3. Finally, there is also some theoretical interest in our approach: the mathematical model is extremely simple, yet it proved to be capable of describing A-scans and processing information needed for damage detection.

6 Summary

The paper addressed the possibility of parameter calibration and damage detection based on the time dependent response of a structure under acoustic excitation. The work extended from actual experimental impulse echo measurements to establishing a mathematical model of wave propagation and reflection, solving the model equations by an efficient numerical procedure based on methods from Fourier analysis, and finally calibrating the model parameters by comparing the computed response to (features) of the measured data. The calibration was done by deterministic optimization and by a Bayesian approach using the Metropolis-Hastings algorithm. The posterior distribution of the parameters can be used to design a hypothesis test detecting damage and its location on the structure. An analogous method for three-dimensional elastic solids has been worked out in [35], announced in [36] and will be the topic of a future publication.

Acknowledgements

Foremost we wish to thank our industrial partner INTALES GmbH Engineering Solutions, Natters, Austria, for continuous support and advice and for providing the plate specimen. We thank Jonathan Halmen for carrying out the ultrasonic measurements and also for installing and programming the control of the positioning table. We thank Peter Paulini for providing the equipment and for valuable advice on ultrasonic flaw detection. The project was supported by grant No. 4602529, Bridge Program, The Austrian Research Promotion Agency (FFG), together with INTALES GmbH Engineering Solutions. The second author acknowledges support through grant P-27570-N26 of The Austrian Science Fund (FWF). The computational results presented have been achieved (in part) using the HPC infrastructure LEO of the University of Innsbruck.