# An optimisation method for separating and rebuilding one-dimensional dispersive waves from multi-point measurements

Marie-Noël Bussac, Pierre Collet, Gérard Gary, Ramzi Othman

► To cite this version:

Marie-Noël Bussac, Pierre Collet, Gérard Gary, Ramzi Othman. An optimisation method for separating and rebuilding one-dimensional dispersive waves from multi-point measurements: Application to elastic or viscoelastic bars. Journal of the Mechanics and Physics of Solids, 2002, 50 (2), pp.321-349. 10.1016/S0022-5096(01)00057-6. hal-01005238

HAL Id: hal-01005238

https://hal.science/hal-01005238

Submitted on 3 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An optimisation method for separating and rebuilding one-dimensional dispersive waves from multi-point measurements. Application to elastic or viscoelastic bars

Marie-Noëlle Bussac<sup>a</sup>, Pierre Collet<sup>a</sup>, Gérard Gary<sup>b</sup>, Ramzi Othman<sup>b</sup>

<sup>a</sup>Laboratoire de Physique Théorique, UMR 7644, Ecole Polytechnique, 91128 Palaiseau, France

<sup>b</sup>Laboratoire de Mécanique des Solides, UMR 7649, Ecole Polytechnique, 91128 Palaiseau, France

When using a classical SHPB (split Hopkinson pressure bar) set-up, the useful measuring time is limited by the length of the bars, so that the maximum strain which can be measured in material testing applications is also limited. In this paper, a new method with no time limits is presented for measuring the force and displacement at any station on a bar from strain or velocity measurements performed at various places on the bar. The method takes the wave dispersion into account, as must inevitably be done when making long time measurements. It can be applied to one-dimensional and single-mode waves of all kinds propagating through a medium (flexural waves in beams, acoustic waves in wave guides, etc.). With bars of usual sizes, the measuring time can be up to 50 times longer than the time available with classical methods. An analysis of the sensitivity of the results to the accuracy of the experimental data and to the quality of the wave propagation modelling was also carried out. Experimental results are given which show the efficiency of the method.

Keywords: Kolsky bar; C. Impact testing; B. Stress wave; B. Viscoelastic material; A. Dynamics

## 1. Introduction

The SHPB (split Hopkinson pressure bar) has become a standard experimental technique for performing tests under dynamic loading conditions. Its success is mainly due to the accuracy of the measurements it yields. This technique is based on the work of
Hopkinson (1914), who recorded a pressure-pulse profile using a long thin bar, and this approach has been widely adopted since the critical study subsequently published by Davies (1948). The practical configuration consisting of two long bars with a short specimen between them, which is widely used nowadays, was introduced by Kolsky (1949). After being initially developed for compression tests, the technique was applied to tensile loading (Harding et al., 1960) and torsion loading situations (Duffy et al., 1971). To improve the accuracy of the basic force and displacement measurements, wave dispersion effects in elastic bars have been studied (Davies, 1948; Yew and Chen, 1978; Follansbee and Franz, 1983; Gorham, 1983), and interactive simulation methods have been proposed to obtain an exact wave shift (Zhao and Gary, 1996). Other aspects involving the analysis of the specimen response, such as three-dimensional effects (Davies and Hunter, 1963; Klepaczko, 1969; Dharan and Hauser, 1970; Bertholf and Karnes, 1975; Malinowski and Klepaczko, 1986) and transient effects (Lindholm, 1964; Conn, 1965; Bell, 1996; Jahsman, 1971) have also been studied in recent decades. Here we deal with the measurement made at the ends of the bars.

Any measuring technique involving the use of bars requires a knowledge of the characteristics of the two elementary waves which propagate in opposite directions. Once they have been characterised, they can be time shifted to the appropriate cross-sections (the bar specimen interface, for example) and all the mechanical values required can be calculated. The SHPB technique involves the use of long bars and a short loading pulse, so that there exists a cross-section where the total incident pulse and the first part of the reflected waves (having the same duration) can be recorded separately. Therefore there exists a maximum observation time which depends on the length of the bar. The measuring time available when a classical SHPB set-up is used is thus limited (Kolsky, 1963) to $\Delta T \le L/c$, where c is the wave speed and L the length of the bar. Consequently, when testing material behaviour at a given average strain rate, the maximum measurable strain is also limited ($e_{max} \le \dot{c}\Delta T$).

To increase the useful measuring time when working with a SHPB, some authors have analysed the multiple reflections occurring in bars. Campbell and Duby (1956) described a method based on a one-dimensional elastic wave theory. Lundberg and Henchoz (1977) have also proposed a simple explicit formula (based on a one-dimensional wave propagation assumption) separating the two elementary waves to measure the particle velocity, and using the signals recorded at two different cross-sections in a bar. This method has also been applied (Lundberg and Blanc, 1988) in a study on the viscoelastic properties of materials and (Lundberg et al., 1990) to the prediction of the wave propagation in a bar with a non-uniform impedance (due to a temperature gradient, for instance) and has been used successfully in high temperature SHPB testing (Bacon et al., 1991, 1994; Bacon and Brun 2000; Lataillade et al., 1994). A “one-point measurement” method has even been proposed, where a free end is used as the second measurement point (Park and Zhou, 1999). However, as several authors have pointed out (Campbell and Duby, 1956; Lundberg and Henchoz, 1977), these methods are valid only if the wave dispersion can be neglected. This condition is satisfied when the bar is thin and the measuring time is short (about 1 ms in these studies with a bar having a diameter of 10 mm).
Zhao and Gary (1994,1997) presented a 2-strain gauge measurement method in which wave dispersion and attenuation are taken into account for use with elastic or viscoelastic bars. With this method, the calculations are carried out in the time domain by means of an iterative formula. The two strain measurements are decomposed into finite time intervals. The first interval at the first strain gauge contains the single wave travelling in the positive direction. This wave is time shifted to the location of the second strain gauge using a Fourier transform, since the frequency components do not all propagate with the same phase velocity. Upon subtracting the results from the strain measured at the second strain gauge, the wave propagating in the negative direction is calculated at the position of this gauge. Likewise, waves propagating in the opposite directions can be evaluated at all possible time intervals by applying an iterative process. Although this method gives satisfactory results, it is an approximate method: the length of the time intervals used corresponds to the time taken by the wave to travel between the two strain gauges, and because of wave dispersion, this time is not exactly the same for all the frequency components of the signal.

Bacon (1999) subsequently proposed a mixed approach. The “two-wave frequency formula” proposed by Zhao and Gary (1997: see also below Eq. (5) in Section 2.2) can be used with most frequency components. In the case of frequency components for which the denominator of this formula equals zero, an iterative calculation is used in the time domain.

In the present paper, an exact mathematical method of solving the same equation in the whole frequency domain is presented. This method was used to analyse the sensitivity of the results to the accuracy of the experimental data and to the quality of the propagation modelling. In addition, a general *n*-strain gauge formula is proposed to stabilise the effects of the noise.

## 2. Two point measurements. A new method for solving the “frequency equations”

### 2.1. Introduction of the wave propagation model

Even with a classical SHPB set-up (where the waves measured propagate once at most along the bar), the one-dimensional wave propagation theory has turned out to be not accurate enough (Davies, 1948; Follansbee and Franz, 1983; Gorham, 1983; Gong et al., 1990; Safford, 1992; Lifshitz and Leber, 1994). With the methods involving longer propagation distances, an accurate propagation theory taking the wave dispersion into account is even more necessary, as pointed out by Lundberg and Henchoz (1977).

The analytical solution to the problem of longitudinal waves propagating along an infinite elastic bar obtained by Pochhammer (1876) and Chree (1889) has been used in several cases (Davies, 1948; Yew and Chen, 1978; Follansbee and Franz, 1983; Gorham, 1983; Gong et al., 1990; Lifshitz and Leber, 1994; Zhao and Gary, 1996). This solution has been extended to bars consisting of linear viscoelastic materials (Zhao and Gary, 1995). Bars of this kind (made of Nylon or PMMA) are needed to improve the impedance ratio when soft materials such as polymer foam are studied. As mentioned
by Davies (1948), it has been confirmed that the infinite bar theory can be applied with sufficiently accurate results to the classical SHPB set-up.

It is therefore natural to use this wave propagation theory to solve problems of the kind we are dealing with here. This approach relies on the assumption that the constitutive law governing the bars is linear, so that all the mechanical variables can be described in terms of their harmonic components. The pulse propagating along the bar can be expressed in space and time by a Fourier integral involving the initial pulse, the phase velocity and the damping (Hunter, 1960). We will use the Fourier transform defined by

$$ \tilde{g}(\omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} g(t)e^{-i\omega t} dt. $$

The inverse Fourier transform is then given by

$$ g(t) = \int_{-\infty}^{\infty} \tilde{g}(\omega)e^{i\omega t} dt. $$

For instance, if we take an ascending (descending) wave propagating along the axis of the bar (x-axis), the associated axial component $\epsilon_{asc}(x, t)$ of the strain tensor $\underline{\epsilon}_{asc}(x, t)$, at any point x and time t can be expressed as follows, where $\tilde{\epsilon}$ denote the frequency component of the Fourier transform

$$ \epsilon_{asc}(x,t) = \int_{-\infty}^{+\infty} \tilde{\epsilon}_{asc}(\omega)e^{i[\omega t - \tilde{\xi}(\omega)x]} d\omega. \quad (1) $$

Likewise, for the descending wave:

$$ \epsilon_{des}(x,t) = \int_{-\infty}^{+\infty} \tilde{\epsilon}_{des}(\omega)e^{i[\omega t + \tilde{\xi}(\omega)x]} d\omega. \quad (2) $$

As strain and velocity measurements are always made at the surface of the bars, their values will be taken at this location.

The dispersion equation $\tilde{\xi} = \xi(\omega)$ between the wave number $\tilde{\xi}$ and the angular frequency $\omega$, which describes the propagation of each frequency component, is defined by the so-called frequency equation based on the solution to the problem of three-dimensional waves propagating along an infinite bar. When dealing with viscoelastic bars, the wave number $\tilde{\xi}$ is a complex function.

Using Eqs. (1) and (2) to describe two points A and B, where a and b are the abscissas of sections A and B, we obtain

$$ \tilde{\epsilon}_{ascB}(\omega) = \tilde{\epsilon}_{ascA}(\omega)e^{-i\tilde{\xi}(\omega)(b-a)}. \quad (3) $$

This shows that ascending waves at two points are simply correlated, and the same applies to descending waves. The wave is shifted between A and B by multiplying the frequency components by a term depending on the dispersion relation.

## 2.2. Presentation of the frequency formulae

The frequency components of the strain can be expressed as the sum of the components of “ascending” and “descending” waves:
$$ 
\begin{aligned}
\tilde{e}_A(\omega) &= \tilde{e}_{\text{ascA}}(\omega) + \tilde{e}_{\text{desA}}(\omega), \\
\tilde{e}_B(\omega) &= \tilde{e}_{\text{ascB}}(\omega) + \tilde{e}_{\text{desB}}(\omega).
\end{aligned}
\qquad (4) $$

Using Eq. (3), we can express the “descending” wave at A by that at B and also the “ascending” wave at B by that at A, so that with $\Delta = (b - a)$, Eq. (4) leads to

$$ \tilde{e}_{\text{ascA}}(\omega) = \frac{\tilde{e}_B(\omega) - \tilde{e}_A(\omega)e^{i\xi(\omega)\Delta}}{e^{-i\xi(\omega)\Delta} - e^{i\xi(\omega)\Delta}}, \qquad (5a) $$

$$ \tilde{e}_{\text{desB}}(\omega) = \frac{\tilde{e}_A(\omega) - \tilde{e}_B(\omega)e^{i\xi(\omega)\Delta}}{e^{-i\xi(\omega)\Delta} - e^{i\xi(\omega)\Delta}}. \qquad (5b) $$

Both ascending and descending waves can therefore be determined at any point on the bar. In particular, at point $x = 0$:

$$ \tilde{e}_{\text{ascO}}(\omega) = A(\omega) = \tilde{e}_{\text{ascA}}(\omega)e^{i\xi(\omega)a} = \frac{\tilde{e}_B(\omega)e^{i\xi(\omega)a} - \tilde{e}_A(\omega)e^{i\xi(\omega)b}}{e^{-i\xi(\omega)\Delta} - e^{i\xi(\omega)\Delta}}, \qquad (6a) $$

$$ \tilde{e}_{\text{desO}}(\omega) = B(\omega) = \tilde{e}_{\text{desB}}(\omega)e^{-i\xi(\omega)b} = \frac{\tilde{e}_A(\omega)e^{-i\xi(\omega)b} - \tilde{e}_B(\omega)e^{-i\xi(\omega)a}}{e^{-i\xi(\omega)\Delta} - e^{i\xi(\omega)\Delta}}. \qquad (6b) $$

At any point, the strain is now simply written in the form

$$ \tilde{e}_X(\omega) = A(\omega)e^{-i\xi(\omega)x} + B(\omega)e^{i\xi(\omega)x}. \qquad (7) $$

In the case of a linear elastic or viscoelastic material and a single (ascending) wave, the stress and the particle velocity are related to the strain by the linear relations

$$ \tilde{\sigma}_{\text{asc}}(\omega) = E(\omega)\tilde{e}_{\text{asc}}(\omega), \qquad (8) $$

$$ \tilde{v}_{\text{asc}}(\omega) = \frac{-\omega}{\xi(\omega)} \tilde{e}_{\text{asc}}(\omega), \qquad (9) $$

where, in the most general case, $E(\omega)$ is the complex Young’s modulus.

As explained above, these parameters are measured at the bar surface. According to the dispersion considered in the following, the stress and the particle velocity are not uniform along the radius of the bar, at any section, except at a free end where the stress is zero. This has been discussed by Davies (1948). This aspect of the question is not directly involved in any of the subsequent calculations. It is important to mention here that the final measurements made with a SHPB, in this version as well as in the classical one, generally consider that average forces and displacements can be deduced to a very good approximation (Davies, 1948) from formulae (8) and (9). In that case the stress and the velocity are assumed to have negligible variations across the section of the bar.

The general formulas for the stress and the particle velocity at a point $x$ are:

$$ \tilde{\sigma}_X(\omega) = E(\omega)(A(\omega)e^{-i\xi(\omega)x} + B(\omega)e^{i\xi(\omega)x}), \qquad (10) $$

$$ \tilde{v}_X(\omega) = \frac{-\omega}{\xi(\omega)} (A(\omega)e^{-i\xi(\omega)x} - B(\omega)e^{i\xi(\omega)x}). \qquad (11) $$
## 2.3. Equivalence between the frequency formula and the “Lundberg–Henchoz” formulae

As mentioned above, Lundberg and Henchoz (1977) have proposed the following explicit recursive formulas (based on a one-dimensional wave propagation assumption) giving the strain and the particle velocity at one point (say the origin), using two strain signals recorded at two different cross-sections in a bar. As mentioned by these authors, it is easy with these formulas to separate the two elementary waves and to measure the particle velocity, using two signals recorded at two different cross-sections in a bar.

$$ \epsilon_0(t) = \epsilon_0(t - 2T) + \epsilon_A(t + T_A) - \epsilon_A(t + T_A - 2T_B) \\ + \epsilon_B(t - T_B) - \epsilon_B(t - T_B + 2T_A) \quad (12) $$

$$ v_0(t)/c = v_0(t - 2T)/c - v_A(t + T_A) - v_A(t + T_A - 2T_B) \\ + v_B(t - T_B) + v_B(t - T_B + 2T_A) \quad (13) $$

where $T_A = a/c$, $T_B = b/c$, $T = T_B - T_A$, and c the elastic wave speed.

These formulas (12) and (13) are equivalent to formulas (7) and (11). They can easily be derived from the latter for the non-dispersive situation.

From Eqs. (6a) (6b) and (7), assuming that $0 < a < b < L$, (where L is the length of the bar), it follows that

$$ \tilde{\epsilon}_X(\omega) = \tilde{\epsilon}_A(\omega) \frac{e^{-i\xi(\omega)(x-b)} - e^{i\xi(\omega)(x-b)}}{e^{-i\xi(\omega)(a-b)} - e^{i\xi(\omega)(a-b)}} + \tilde{\epsilon}_B(\omega) \frac{e^{-i\xi(\omega)(x-a)} - e^{i\xi(\omega)(x-a)}}{e^{-i\xi(\omega)(b-a)} - e^{i\xi(\omega)(b-a)}} $$

or

$$ \tilde{\epsilon}_X(\omega) = \tilde{\epsilon}_A(\omega) \frac{\sin(\xi(\omega)(x-b))}{\sin(\xi(\omega)(a-b))} + \tilde{\epsilon}_B(\omega) \frac{\sin(\xi(\omega)(x-a))}{\sin(\xi(\omega)(b-a))}, \quad (14) $$

$\tilde{\epsilon}_A(\omega)$ and $\tilde{\epsilon}_B(\omega)$ are analytical bounded functions in the lower half-plane of complex numbers $\text{Im}(\omega) < 0$. This property expresses causality in time, since the stress has to be 0 for $t < 0$. It is also assumed that $\xi(\omega)$ is an analytical function in the lower half-plane and maps this half plane into itself. Strictly speaking, these conditions cannot be satisfied in the case of a cylindrical bar if only one of the modes involved is considered: causality can only be ascertained from the complete solution of the equation of propagation including all the modes. We assume, however, in the present paper that considering only the first mode is a legitimate approximation.

We can now write the same formula in the lower half-plane in the following form:

$$ \tilde{\epsilon}_X(\omega) = \tilde{\epsilon}_A(\omega) \frac{e^{-i\xi(\omega)(b-a)} - e^{i\xi(\omega)(x-b)}}{1 - e^{-2i\xi(\omega)(b-a)}} \\ - \tilde{\epsilon}_B(\omega) \frac{e^{-i\xi(\omega)(b-a)} - e^{i\xi(\omega)(x-a)}}{1 - e^{-2i\xi(\omega)(b-a)}}. $$
Since $b-a>0$, if $\operatorname{Im}(\omega) < 0$ then $\operatorname{Im}(\xi(\omega)(b-a)) < 0$, and hence

$$|e^{-2i\xi(\omega)(b-a)}| < 1.$$

It is then possible to develop the fractions in convergent series. With $\Delta = b-a$, it follows that

$$ \begin{aligned} \tilde{e}_X(\omega) = & \sum_{n=0}^{n=\infty} \tilde{e}_A(\omega) e^{-i\xi(\omega)(x-a) - 2in\xi(\omega)\Delta} - \sum_{n=0}^{n=\infty} \tilde{e}_A(\omega) e^{-i\xi(\omega)(2b-a-x) - 2in\xi(\omega)\Delta} \\ & + \sum_{n=0}^{n=\infty} \tilde{e}_B(\omega) e^{-i\xi(\omega)(b-x) - 2in\xi(\omega)\Delta} - \sum_{n=0}^{n=\infty} \tilde{e}_B(\omega) e^{-i\xi(\omega)(x+b-2a) - 2in\xi(\omega)\Delta}. \end{aligned} \quad (15) $$

It is now possible to obtain the time evolution using the inverse Fourier transform. In the non-dispersive case, this yields

$$ \begin{aligned} \varepsilon_X(t) = & \sum_{n=0}^{\infty} \varepsilon_A(t - (x - a + 2n\Delta)/c) - \sum_{n=0}^{\infty} \varepsilon_A(t - (2b - a - x + 2n\Delta)/c) \\ & + \sum_{n=0}^{\infty} \varepsilon_B(t - (b - x + 2n\Delta)/c) - \sum_{n=0}^{\infty} \varepsilon_B(t - (x + b - 2a + 2n\Delta)/c). \end{aligned} \quad (16) $$

We can write this formula replacing $t$ by $t - 2\Delta/c$:

$$ \begin{aligned} \varepsilon_X(t - 2\Delta/c) = & \sum_{n=0}^{\infty} \varepsilon_A(t - (x - a + 2(n+1)\Delta)/c) - \dots \\ \varepsilon_X(t - 2\Delta/c) = & \sum_{n=1}^{\infty} \varepsilon_A(t - (x - a + 2n\Delta)/c) - \dots \end{aligned} \quad (17) $$

For a finite value of $t$, the terms with larger $N$ vanish, since $\varepsilon_A$ and $\varepsilon_B$ vanish for $t < 0$. Only a finite number $N$ of terms needs to be considered ($N < (ct/2\Delta) + 1$). (17) can then be subtracted from (16). All the terms cancel except the first term of (16). Therefore

$$ \begin{aligned} \varepsilon_X(t) = & \varepsilon_X(t - 2\Delta/c) + \varepsilon_A(t - x + a/c) - \varepsilon_A(t + x + a/c - 2b/c) \\ & + \varepsilon_B(t + x - b/c) - \varepsilon_B(t - x - b/c + 2a/c). \end{aligned} \quad (18) $$

For $x=0$, $T_A = a/c$, $T_B = b/c$ and $T = T_B - T_A = \Delta/c$, this formula is exactly the same as formula (12) proposed by Lundberg and Henchoz (1977): in their paper, $a$ and $b$ are replaced by $x_1$ and $x_2$ and $A$ and $B$ by 1 and 2.

Formula (13) can be derived in a similar way.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic diagram of the integration method.</figcaption></figure>

## 2.4. A new deconvolution method for solving the "frequency formula"

As observed by Bacon (1999), the two functions in formulas (5a) and (5b) are not defined if the denominator equals zero, as occurs when $\xi(\omega)=n\pi/\Delta$. The technique used by Bacon is based on a frequency component discretisation approach where unsuitable frequencies can be processed in a special way.

Since the poles of the denominator are on the real axis, we propose an integration method that operates in the complex domain. If causality is ensured by choosing an integration path in the lower half-plane, and since $\xi(\omega)$ is an analytical function, the results of this calculation will not depend on the integration path (Morse and Feshbach, 1953).

With formula (6a), the method involves going through the following steps (Fig. 1). $\tilde{f}$ and $\tilde{f}^{-1}$ denote the Fourier transform operation and the inverse Fourier transform operation, respectively. The frequency $\omega'$ is in the complex domain: $\omega' = \omega - i\eta$, $\omega$ is real and $\eta$ is a small positive real number.

The dispersion relation $\xi(z)$ is then calculated in the complex domain for the complex values of $z = \omega - i\eta$. When the dispersion relation is known only on the real axis, its value near the real axis can be approximated by the first term in its Taylor expansion.
The same technique is used on both formulas (6a) and (6b). Knowing only the "ascending" strain at point A, the "descending" strain at the same point is also obtained by determining the difference with the measured strain.

This technique can also be used in the case of non-dispersive media, where the dispersion relation is simply written as

$$\xi(\omega) = \omega/c_0,$$

where $c_0$ is the wave speed.

## 2.5. Numerical check of the validity

The method will first be illustrated by describing a purely numerical simulation in order to prevent any experimental artefacts from interfering with the checking process.

This procedure consists in building the waves that are produced by a given known loading speed. From the characteristics of these waves, simulated strain records can be built at given stations on the simulated bar.

These strains are then used as data for the deconvolution method presented above. To check the validity of the method, it is necessary to ascertain whether the same limit conditions are recovered at the ends of the bar as those used to build the waves. In particular, the stress at the free end of the bar has to be equal to zero.

### 2.5.1. Simulation of a simple test

The loading that we chose to simulate here corresponds to an idealised situation where the impact of a striker occurs at one end of a bar. The bar is then subjected to a given constant step-speed with a finite duration at one (say the input) end. Before and after the loading, the input end is stress free and the other (say the output) end is always stress free.

When the loading time is shorter than the time taken by the waves to make a round trip along the bar and back (say $\Delta T$), this corresponds to the idealised impact of a striker shorter than the bar, which resembles the example proposed by Bacon (1999). It leads to waves the shape of which is familiar to people dealing with Hopkinson bars.

In addition, this case allows for a step by step construction of the wave: using relation (9), the speed step is converted into a quasi-strain step loading. Using relation (3), the wave is then shifted along the bar. At the free ends (after the loading, both ends are stress free), the wave reflects and the sign of the strain changes.

The relation between the strain and the particle velocity (9) is frequency dependent. For a speed-step longer than $\Delta T$, the exact construction of the reflected wave at the input end would need to be analysed in the frequency domain. As in the case dealt with here to illustrate the method, the result does not depend on the loading, and we shall apply a speed step shorter than $\Delta T$ for the sake of simplicity.

The case of a 3-m long bar is investigated, as shown in Fig. 2.

Strain gauge stations 1, 2 and 3 are 0.7, 1.41 and 2.3 m, respectively, from the impact end. Station 2 is not placed in the middle of the bar to avoid any resonance, as will be discussed below in Section 3.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Sketch of the bar.</figcaption></figure>

The simulated loading is produced by the impact of a striker made of the same material and having the same diameter as the bars. The loading speed is 12 m/s. (the constructed speed at the striker end of the bar therefore has a value of 6 m/s).

In real test situations, the signals are recorded digitally with a 1 µs sampling step. All the following numerical simulations will therefore be carried out with this sampling step.

Three different situations will be investigated here, depending on the constitutive model used to describe the bars.

In the first case (case 1), the bar and the striker (1 m long) are assumed to be elastic and ideally thin so that the waves can propagate without any dispersion. The wave speed is $c_0 = 5000$ m/s. The loading time is therefore 400 µs.

In the second case (case 2), the bar and the striker (2.5 m long) are assumed to be elastic, with a diameter of 40 mm. To perform the wave shifting process used to build the simulated waves, it is necessary to determine the dispersion relation $\xi(\omega)$. As described in Zhao and Gary (1995), this dispersion relation is calculated as the first mode of the Pochhammer-Three equation recalled below:

$$ (2\alpha/a)(\beta^2 + \xi^2)J_1(\alpha.a)J_1(\beta.a) - (\beta^2 - \xi^2)^2J_0(\alpha.a)J_1(\beta.a) - 4\xi^2\alpha.\beta.J_1(\alpha.a)J_0(\beta.a) = 0, \quad (19) $$

where

$$ \alpha^2 = \frac{\rho\omega^2}{\tilde{\lambda}(\omega) + 2\tilde{\mu}(\omega)} - \xi^2, \quad \beta^2 = \frac{\rho\omega^2}{\tilde{\mu}(\omega)} - \xi^2, $$

$J_0, J_1$ are zero and first order Bessel's functions. $a$ is the radius of the bar. $\tilde{\lambda}$ and $\tilde{\mu}$ are the complex Lamé coefficients of the material of the bar. When the bar is purely elastic, the formula is the same and the Lamé coefficients are constants.

The values for the material are similar to those of aluminum (velocity of low frequency waves: $c_0 = 5000$ m/s, Poisson's ratio $\nu = 0.34$, density = 2800 kg/m³). The duration of the loading pulse is 400 µs.

The subsequent dispersion relation $\xi(\omega)$ is very similar to the relation found in Bancroft (1941) and Davies (1948) for the phase velocity $c(\omega)$ as a function of the angular frequency.

In the third case (case 3), the striker (1 m long) and the bars are both assumed to be viscoelastic and to have a diameter of 40 mm. The dispersion relation for the bar is calculated with values similar to those of Nylon. The complex Young's modulus is the complex modulus of the standard rheological model presented in Fig. 3a with the following data: $E = 4.395 \times 10^9$ Pa, $E_v = 2 \times 10^{10}$ Pa, $\eta = 9 \times 10^5$ Pa s. The Poisson's ratio is not frequency dependent: $\nu = 0.4$.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) A standard rheological model for Nylon. (b) Dispersion relation and damping factor for a simulated viscoelastic bar.</figcaption></figure>

The corresponding velocity of the low frequency waves is 1770 m/s. The duration of the loading pulse is therefore 1130 µs.

The subsequent dispersion relation $\xi(\omega)$ is shown in Fig. 3b, where the phase velocity and the damping are plotted as functions of the angular frequency. The damping factor is the ratio of the amplitude of a single frequency signal (after propagating 1 m) to its initial value.

This model is not very different from that used to process the experimental data for Nylon (see Section 5). It shows the existence of a fairly strong damping factor. Consequently, angular frequencies greater than 30 000 rad/s. (frequencies of more than 5 kHz) will almost disappear after the wave has propagated a few meters along the bar.

Basic deformation waves are then constructed at strain gauges 1 and 3, allowing, by superposition, the computation of the strains at these gauges. The corresponding, computed strains at gauges 1 and 3, for cases 1, 2 and 3 are shown in Figs. 4, 5 and 6, respectively.

### 2.5.2. Reconstruction of basic waves and given loading conditions at ends of the bar

Using the data shown above and formulas (6a), (6b) and (7) with the algorithm presented in Fig. 1, the stresses and velocities are reconstructed at both ends of the bar.

Figs. 4, 5 and 6 show computed strains at gauges 1 and 3 and the reconstructed stresses associated with cases 1, 2 and 3, respectively.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Ideally *thin* elastic bar: Simulated strains at stations 1 and 3 and reconstructed stresses at the ends of the bar with the “2-strain gauge method”. (The stress at the free end has been shifted down by 10 MPa to make the figure more readable).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Thick elastic bar: Simulated strains at stations 1 and 3 and reconstructed stresses at the ends of the bar with the “2-strain gauge method”. (The stress at the free end has been shifted down by 10 MPa to make the figure more readable).</figcaption></figure>

Case 1: thin elastic bar.

Case 2: thick elastic bar.

The above illustration is based on the use of a short striker to make the figures, showing the strains, easier to understand. As the purpose of the method is to overcome the limitations of the classical Hopkinson bar method, the following example will be given with a longer striker, for the case of the elastic bar. In this case, there is overlap between the ascending and descending waves (as seen in Fig. 5). With the same 3-m long bar, the striker is 2.5 m long. The duration of the loading pulse is 1000 µs.

Case 3: thick viscoelastic bar.

For the three cases investigated, it turns out, as expected, that the reconstructed force at the output end is still equal to zero, and the reconstructed force at the input end is also still equal to zero after the loading step.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Viscoelastic bar: Simulated strains at stations 1 and 3 and reconstructed stresses at the ends of the bar with the "2-strain gauge method". (The stress at the free end has been shifted down by 1 MPa to make the figure more readable).</figcaption></figure>

Using formula (11), the velocities at the ends of the bar (or at any other point) were also calculated. The results are again in perfect agreement with those expected.

### 3. Sensitivity of the two point method to the accuracy of measurements

Strain gauge measurements are intrinsically subject to electronic noise and suffer from the fact that the gain in the amplifying system cannot be accurately determined. We have also observed during experimental tests that the expected zero strain (that occurs after the onset of the loading when the striker is short enough) is sometimes not associated with an exact zero voltage measurement. Since we have no clear explanation for this, we will also look at the possible effects of this factor.

These effects will be investigated using case 2 (elastic bar) and a smaller striker (1 m long). It is indeed to be expected that the high frequency damping observed in case 3 will decrease the effects of the noise so that the worst possible situation as far as the noise is concerned is the case of a purely elastic bar.

#### 3.1. *Initial signal processing*

During this investigation, we must keep in mind the experimental situation, and the following simulations have to approximate the real case. In addition, the algorithm to be tested (and checked) here should be the same as that used in a real situation.

Fig. 7 shows a real signal recorded by the middle *strain gauge* on a real aluminium bar and the frequency spectrum of this signal.

The peaks observed in the frequency spectrum correspond to the resonance of the bar (and reflect the time taken by the waves to travel along the bar). It should be stressed here that almost no energy is detected at frequencies higher than 30 kHz. Before being further processed, all the signals measured can then be safely low-pass
<figure><img src="image_9.png" /><figcaption>Fig. 7. Impact of a steel ball (diameter 20 mm) on an aluminium bar (diameter 40 mm, length 3 m).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. “2-strain gauge method” reconstructed stress at the impact end of the bar. Ideally elastic bar. (a) Influence of noise; (b) Influence of the imperfect gain in one amplifier (the gain in the second amplifier was increased by 2%).</figcaption></figure>

filtered at frequencies of less than 100 kHz. The same procedure will be adopted in the subsequent simulations.

## 3.2. Noise perturbations

White noise, the mean amplitude of which is 1% of the maximum strain, is added to the simulated signals. The strain at any station j becomes

$$ \hat{e}_j(t) = \varepsilon_j(t) + w_j(t), $$

where $\varepsilon_j(t)$ is the exact constructed strain and $w_j(t)$ is the white noise. The white noises at different stations are not correlated but have the same amplitude.

The reconstructed force at the impact end of the bar is shown in Fig. 8a.
This numerical simulation shows that the noise is amplified, as often occurs when using a mathematical deconvolution technique. The effects of the noise can be investigated using formula (16). This expression is linear, and hence to obtain the noise on the reconstructed function $e_x(t)$, it suffices to replace $e_A(t)$ and $e_B(t)$ by white noise. We recall here that by causality, Eq. (16) contains only a finite number of terms of the order $ct/\Delta$ (see also Eq. (12)). Therefore, since white noises are independent at different times, the variance of the noise in $e_x(t)$ is proportional to $t$ (for large $t$), and hence the average amplitude is proportional to $\sqrt{t}$.

### 3.3. Imprecise knowledge of the amplifier gain

The value of the gain applied to the second strain gauge signal (at station 3) is increased by 2%, and the gain at the first strain gauge station is still equal to 1.

Reconstructed force at the impacted end of the bar is shown in Fig. 8b.

Although the effects on the reconstructed force were significant, a rather low frequency response was induced. This frequency is directly correlated with the distance between the two gauges. The strain (proportional to the stress) at position $x$ is given by formula (14).

If the amplification of the strain at point $B$ has been imprecisely set, the strain $\hat{e}_B(\omega)$ can be written as

$$ \hat{e}_B(\omega) = (1 + \gamma)\tilde{e}_B(\omega). $$

The difference between the exact value and that calculated with an imprecise amplification is

$$ \Delta e_x(\omega) = \gamma \tilde{e}_B(\omega) \frac{\sin(\xi(\omega)(x-a))}{\sin(\xi(\omega)(b-a))}. $$

A resonance will occur at values of the denominator around zero:

$$ \xi(\omega)(b-a) = n\pi. $$

Neglecting the damping, this gives

$$ \omega = c(\omega)n\pi/(b-a). $$

The first resonance in our example is obtained for $\omega = 9800$ rad/s, which corresponds to a period of 0.64 ms, and is similar to the period observed in Fig. 8b.

This effect of an imprecise gain adjustment has a very typical signature which should make it easy to detect if this occurs. Another point is that the gains can be accurately adjusted in a real test situation, before running the deconvolution method: as real signals start from zero, part of the first ascending wave can always be observed at each gauge station (as can be seen in Figs. 4a–c and 7). Since we know the dispersion relation,
<figure><img src="image_11.png" /><figcaption>Fig. 9. “2-strain gauge method” reconstructed velocity at impact end of the bar. Ideally elastic bar; linear shift of the second amplifier.</figcaption></figure>

the first strain gauge signal can then be shifted to exactly match the other strain gauge signals, which makes it possible to adjust the gains accordingly. The effects of poorly determined gains will be left aside from now on.

### 3.4. Imprecise null strain measurement

The zero strain level of the strain gauge at station 3 is shifted by a value increasing with time, starting at zero and reaching up to 2% of the maximum strain at the end of the test.

We have indeed observed this kind of distortion in real tests. It is a very low frequency phenomenon that can be simulated in this way.

Reconstructed velocity at the impacted end of the bar is shown in Fig. 9, where it is compared with the reconstructed velocity from perfect signals.

A non-physical variation of the average reconstructed velocity is observed. It has been checked that the effect of a shifted measured strain therefore turns out to be much stronger in the case of reconstructed velocities than in the case of the reconstructed stresses.

As in the case of the noise, we again consider Eqs. (12) and (13). Since these equations are linear, a variation $\delta\epsilon_B$ alone leads to variations $\delta\epsilon_0$ and $\delta v_0$ satisfying

$$\delta\epsilon_0(t) = \delta\epsilon_0(t - 2T) + \delta v_B(t - T_B) - \delta v_B(t - T_B + 2T_A)$$

and

$$\delta v_0(t) = \delta v_0(t - 2T) + c(\delta v_B(t - T_B) + \delta v_B(t - T_B + 2T_A)).$$

We now assume that $\delta\epsilon_B(t)$ is proportional to $t\theta(t)$, where $\theta$ is the Heaviside function (equal to 0 for $t < 0$ and to 1 otherwise). Since $\delta\epsilon_0(t) = \delta v_0(t) = 0$ for $t < 0$, we can solve the above equations by iteration. One can easily check that $\delta\epsilon_0(t)$ increases linearly with time and $\delta v_0(t)$ quadratically.
<figure><img src="image_12.png" /><figcaption>Fig. 10. "2-strain gauge method" reconstructed stresses at both ends of the bar. Ideally elastic bar. (a) Influence of imperfect gauge positions; (b) Influence of the imperfect dispersion relation.</figcaption></figure>

### 3.5. Imprecise knowledge of the positions of the gauges

In a real situation, it should be possible to measure the positions of the strain gauge stations with a satisfactory level of accuracy. With a 3-m long bar, an accuracy of 1 mm can be expected. To investigate the influence of this parameter, a simulation was carried out with an error of this magnitude in the gauge position. The waves were then reconstructed with the position of station 1 at 0.701 m (overestimated by 1 mm) and station 3 was placed by 2.299 m (underestimated by 1 mm).

Reconstructed forces at the end of the bar are shown in Fig. 10a.

Even a small error (1‰) in the position of the strain gauges can therefore be seen to have had clearly visible effects on the reconstructed forces: the strain gauge positions have to be accurately known. Nevertheless it does not induce perturbations increasing with time.

### 3.6. Imprecise knowledge of the dispersion relation

It is very difficult to exactly determine the dispersion relation, as previously established by other authors (Gorham, 1983; Zhao and Gary, 1997). In preliminary attempts to use the two gauge methods in a real test situation, we have noticed a significant influence upon the results of material parameters (in particular Poisson's ratio) used to compute the dispersion relation. It was therefore necessary to investigate the influence of an imprecisely determined dispersion relation.

In the following simulation, the dispersion was only modified by 2‰. Namely, the $\xi(\omega)$ data have been transformed into $0.998\xi(\omega)$.

The reconstructed forces at the ends of the bar obtained using this modified relation are shown in Fig. 10b.

This figure shows that an exact knowledge of the dispersion relation is required when performing long time measurements. This is in line with previous reports (Zhao and Gary, 1997) on studies using an approximate method.
<figure><img src="image_13.png" /><figcaption>Fig. 11. "N(3)-strain gauge method" reconstructed stress at the impact end of the bar. Ideally elastic bar: (a) 1% noise; (b) All the imprecise measurements have been combined.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. "N(3)-strain gauge method" reconstructed velocity at impact end of the bar. Ideally elastic bar; linear shift of the second amplifier.</figcaption></figure>

# 4. Multi-point measurement method

## 4.1. N-strain gauge wave separation method

Previous simulations (Figs. 9a–15b) have shown that various inaccuracies in the measurements or a lack of knowledge of the bar's material properties can affect the results of the deconvolution process. In particular, Fig. 8a shows that the noise is amplified during the process: any signals recorded at two gauge stations will yield a solution for waves propagating in the bar, and in particular, a solution for the loading conditions at the ends of the bar.

It is therefore worth introducing more relevant information by adding extra gauges. This idea was also recently used by Hillström and co-authors in the problem of the experimental determination of the complex modulus of a material (Hillström et al., 2000).
<figure><img src="image_15.png" /><figcaption>Fig. 13. "N(3)-strain/P(1)-velocity method" reconstructed stresses and velocities at the impacted end of the bar. All the imprecise measurements have been combined (including the zero shift).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Scheme of the experimental set-up.</figcaption></figure>

By considering only the effect of white noise, one can reconstruct propagating waves using the maximum likelihood method.

We assume that the measured strain at station j is given by

$$ \hat{\varepsilon}_j(t) = \varepsilon_j(t) + w_j(t), $$

where $w_j(t)$ are statistically independent white noises with the same variance. The maximum number N will be assumed below to be greater than 2. The Maximum Likelihood functional (Van Trees, 1992) is then proportional to F given (for N stations) by

$$ F = \int \sum_{j=1}^{N} (\hat{\varepsilon}_j(t) - \varepsilon_j(t))^2 dt. $$
<figure><img src="image_17.png" /><figcaption>Fig. 15. Short loading pulse, recorded data: strains, displacement at station 4 and corresponding filtered calculated velocity (3 kHz).</figcaption></figure>

Using the Parseval's relation and formula (7), we obtain

$$F = \frac{1}{2\pi} \int \sum_{j=1}^{N} |\hat{e}_j(\omega) - A(\omega)e^{-i\xi(\omega)x_j} - B(\omega)e^{i\xi(\omega)x_j}|^2 d\omega.$$

The maximum likelihood estimator is given by the functions $A(\omega)$ and $B(\omega)$, which minimise $F$. In other words, in a given set of functions (in this case, the solutions of the evolution equation with possibly some additional causality constraints) we look for a minimiser of the maximum likelihood functional. Defining the functions $h_1$, $h_2$ and $g$ (on the real $\omega$ line) by

$$h_1(\omega) = \sum_{j=1}^{N} e^{-i(\xi(\omega) - \bar{\xi}(\omega))x_j}, \quad h_2(\omega) = \sum_{j=1}^{N} e^{i(\xi(\omega) - \bar{\xi}(\omega))x_j},$$

$$g(\omega) = \sum_{j=1}^{N} e^{i(\xi(\omega) + \bar{\xi}(\omega))x_j},$$

we obtain

$$F = \frac{1}{2\pi} \int d\omega \left[ h_1(\omega) \overline{A(\omega)A(\omega)} + h_2(\omega) \overline{B(\omega)B(\omega)} + g(\omega) \overline{A(\omega)B(\omega)} \right. \\ 
+ \overline{g(\omega)B(\omega)A(\omega)} - A(\omega) \sum_{j=1}^{N} \hat{e}_j(\omega) e^{-i\xi(\omega)x_j} - B(\omega) \sum_{j=1}^{N} \hat{e}_j(\omega) e^{i\xi(\omega)x_j} \\ 
\left. - \overline{A(\omega)} \sum_{j=1}^{N} \hat{e}_j(\omega) e^{i\xi(\omega)x_j} - \overline{B(\omega)} \sum_{j=1}^{N} \hat{e}_j(\omega) e^{-i\xi(\omega)x_j} + \sum_{j=1}^{N} \hat{e}_j(\omega) \hat{e}_j(\omega) \right].$$
It is now convenient to use a matrix notation. Let $\langle | \rangle$ denote the (complex) scalar product in dimension two. We define two-dimensional vectors as follows

$$ \mathbf{X}(\omega) = \begin{pmatrix} A(\omega) \\ B(\omega) \end{pmatrix} \quad \text{and} \quad \mathbf{E}(\omega) = \begin{pmatrix} \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{i\overline{\xi}(\omega)x_j} \\ \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{-i\overline{\xi}(\omega)x_j} \end{pmatrix} $$

and a two by two Hermitian matrix:

$$ T(\omega) = \begin{pmatrix} h_1(\omega) & g(\omega) \\ \overline{g}(\omega) & h_2(\omega) \end{pmatrix}. $$

With these notations, we have

$$ F = \frac{1}{2\pi} \int d\omega \left[ \langle \mathbf{X}(\omega) | T(\omega) \mathbf{X}(\omega) \rangle - \langle \mathbf{E}(\omega) | \mathbf{X}(\omega) \rangle - \langle \mathbf{X}(\omega) | \mathbf{E}(\omega) \rangle + \sum_{j=1}^{N} \hat{e}_j(\omega) \overline{\hat{e}_j(\omega)} \right]. $$

After some algebraic manipulations, one obtains

$$ F = K + \frac{1}{2\pi} \int d\omega [\langle \mathbf{X}(\omega) - T(\omega)^{-1} \mathbf{E}(\omega) | T(\omega)(\mathbf{X}(\omega) - T(\omega)^{-1} \mathbf{E}(\omega)) \rangle], $$

where K is a quantity which does not depend on $\mathbf{X}$. Since $T(\omega)$ will be a positive definite matrix (except for $\omega = 0$, see the discussion below), the minimum value of $F$ is reached only for the vector

$$ \mathbf{X}(\omega) = T(\omega)^{-1} \mathbf{E}(\omega) \qquad (20) $$

or more explicitly

$$ A(\omega) = \frac{h_2(\omega) \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{-i\overline{\xi}(\omega)x_j} - g(\omega) \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{-i\overline{\xi}(\omega)x_j}}{h_1(\omega)h_2(\omega) - g(\omega)\overline{g(\omega)}} \quad (20a) $$

and

$$ B(\omega) = \frac{h_1(\omega) \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{-i\overline{\xi}(\omega)x_j} - \overline{g(\omega)} \sum_{j=1}^{N} \hat{e}_j(\omega) \mathbf{e}^{i\overline{\xi}(\omega)x_j}}{h_1(\omega)h_2(\omega) - g(\omega)\overline{g(\omega)}} \quad (20b) $$

Here we have some comments to make. First of all, one can easily check that if, for $j = 1, \dots, N$

$$ \hat{e}_j(\omega) = A'(\omega)\mathbf{e}^{-i\overline{\xi}(\omega)x_j} + B'(\omega)\mathbf{e}^{i\overline{\xi}(\omega)x_j}, $$

one obtains $A = A'$ and $B = B'$.
One notices immediately that (since $\xi(0) = 0$) the denominator is equal to zero for $\omega = 0$, and this is a double root. However the numerators of A and B are also equal to zero for $\omega = 0$ and therefore the pole is either single or absent. We now show that unless resonance occurs, $(\xi(\omega)x_j = \alpha + 2p_j\pi$ for $j = 1, \dots, N$ where the numbers $p_j$ are integers), there are no other zeros in the denominator on the real line.

Note that if a and b are (non zero) complex numbers one has

$$ (\tilde{a}\tilde{b}) \begin{pmatrix} h_1(\omega) & g(\omega) \\ \bar{g}(\omega) & h_2(\omega) \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \sum_{l}^{N} |ae^{-i\xi(\omega)x_l} + be^{i\xi(\omega)x_l}|^2. $$

This last sum can be zero if and only if each term is equal to zero, namely if simultaneously for $j = 1, \dots, N$ we have

$$ e^{2i\xi(\omega)x_j} = -\frac{a}{b}. $$

Since the numbers $x_j$ are real and different, it follows that $a/b$ must have a modulus equal to one, $\xi(\omega)$ must be real and the above resonance conditions must be satisfied. If this is not the case, we conclude that the matrix $T(\omega)$ is invertible for any non zero real $\omega$.

We are aware that although, as demonstrated below, the above estimator gives good control over the noise, it has several drawbacks. First we have had to perform the integration over the whole time axis in the definition of $F$ in order to be able to use the Fourier transform. This requires, however, that the signal should be observed for infinite time and decay to zero fast enough at infinity for convergence of the integral to occur. Although, it is possible to use only a finite time interval in the definition of the maximum likelihood functional, the minimisation procedure then becomes much more involved.

Although we have argued above that there is no singularity of the denominator except at the origin, the presence of complex conjugate terms will generally break the analyticity in the lower half-plane and hence disrupt the causality. We have nevertheless used a line parallel to the real axis and slightly below it as the integration contour in order to circumvent the singularity at the origin (here we assume that the function $\xi(\omega)$ is analytical in the lower half-plane). A more rigorous causal approach would be to minimise $F$ under the constraint that A and B must be analytical and bounded in the lower half plane. This would require mathematical developments which extend far beyond the scope of this paper.

We also checked experimentally that the white noise assumption was a good approximation in the frequency range where the level of the signal was non negligible.

To check the efficiency of this formula, we consider case 2, where an additional strain gauge signal was simulated at station 2 (Fig. 2).

### 4.1.1. White noise

As in Sections 2.1 and 3.1, a white noise, the amplitude of which was 1% of the maximum strain, was added to the simulated signals at the three stations.

The force at the impacted end of the bar reconstructed using formula (20) is shown in Fig. 11.
In comparison with Fig. 8, this figure shows, as was expected, that the noise is no longer amplified.

### 4.1.2. *Imprecise measurements other than of the null strain*

As in a real test, we assume that all the effects occur together, including a white noise. They are then simultaneously introduced into the simulated waves at the three strain gauge stations as follows:

* white noise at the three strain gauge stations: amplitude equal to 1% of the maximum strain;

* gauge position: 0.701 m (instead of 0.700 m) at station 1, 1.4105 m (instead of 1.410 m) at station 2, and 2.299 m (instead of 3.000 m) at station 3;

* imprecise knowledge of the dispersion relation as in Section 3.5.

Using formula (20), the result shown in Fig. 11 (above) is obtained.

The maximum error observed in the force at the impacted end of the bar is less than 4% of the maximum force. The error level can therefore be said to be acceptable. Using formula (20) yields significant improvements in comparison with the 2-strain gauge formula presented above, which has been used by previous authors (Lundberg and Henchoz, 1977; Zhao and Gary, 1997; Bacon, 1999).

### 4.1.3. *Imprecise measurement of the null strain*

The zero strain recorded by the gauge at station 3 was *increased* linearly in time, from zero up to 2% of the maximum strain at the end of the test (as in the 2-strain gauge simulation). At the next simulated measurement station 2, the zero strain level was *decreased* linearly, from zero down to -2% of the maximum strain recorded, at the end of the test.

Despite the fact that formula (20) was not designed for this purpose (i.e., for the imprecise measurement of the null strain), it was used to reconstruct the velocity at the free ends of the bar. This is shown in Fig. 12.

Fig. 12 shows that a significant difference was observed in the case of the velocity. It amounted to around 20% of the average velocity at the end of the bar.

Similar results were observed as with 2-strain gauges.

## 4.2. "N-strain/P-velocity gauge" wave separation method

Using more than 2-strain gauges obviously does not satisfactorily solve the errors introduced by imprecise null strain measurements. The most significant effect of this factor was that observed on the velocities at the bar ends. In a previous paper, Zhao and Gary (1997) have shown that performing strain measurements at one station and velocity measurements at another station significantly improves the stability of the reconstructed velocities. In line with this idea, we investigated the use of both strain and velocity measurements.

Using the maximum likelihood method and formulas (7) and (11), a "N-strain/P-velocity gauge" formula can be obtained in a similar way to the "N-strain gauge
formula". One gets

$$
\begin{align}
A(\omega) &= \frac{h_2(\omega)E_1(\omega) - g(\omega)E_2(\omega)}{h_1(\omega)h_2(\omega) - g(\omega)g(\omega)}, \\
B(\omega) &= \frac{h_1(\omega)E_2(\omega) - \overline{g(\omega)}E_1(\omega)}{h_1(\omega)h_2(\omega) - g(\omega)g(\omega)} \tag{21}
\end{align}
$$

with

$$
h_1(\omega) = \sum_{j=1}^{N} e^{-i(\xi(\omega) - \bar{\xi}(\omega))x_j} + \left|\frac{c(\omega)}{c_0}\right|^2 \sum_{k=1}^{p} e^{-i(\xi(\omega) - \bar{\xi}(\omega))x_k},
$$

$$
h_2(\omega) = \sum_{j=1}^{N} e^{i(\xi(\omega) - \bar{\xi}(\omega))x_j} + \left|\frac{c(\omega)}{c_0}\right|^2 \sum_{k=1}^{p} e^{i(\xi(\omega) - \bar{\xi}(\omega))x_k},
$$

$$
g(\omega) = \sum_{j=1}^{N} e^{i(\xi(\omega) + \bar{\xi}(\omega))x_j} - \left|\frac{c(\omega)}{c_0}\right|^2 \sum_{k=1}^{p} e^{i(\xi(\omega) + \bar{\xi}(\omega))x_k},
$$

$$
E_1(\omega) = \sum_{j=1}^{N} e^{i \bar{\xi}(\omega) x_j} \hat{e}_j(\omega) - \frac{1}{c_0^2} \sum_{k=1}^{p} e^{i \bar{\xi}(\omega) x_k} c(\omega) \hat{v}_k(\omega),
$$

$$
E_2(\omega) = \sum_{j=1}^{N} e^{-i\bar{\xi}(\omega)x_j} \hat{e}_j(\omega) + \frac{1}{c_0^2} \sum_{k=1}^{p} e^{-i\bar{\xi}(\omega)x_k} c(\omega) \hat{v}_k(\omega).
$$

This formula was tested with simulated measurements at the same previous three
stations. Velocity measurements were carried out at station 3 (instead of strain mea-
surements). At stations 1 and 2, the same simulated signals as those mentioned in
Section 3 are used, where the imprecise measurement of the null strain (used in Sec-
tion 3.3) is added to other imprecise measurements (including the imprecise knowledge
of the measurement station positions).

A white noise (2% of the striker speed) is also added to the simulated velocity.

Formula (21) is used with these simulated signals to reconstruct the stresses and
velocities at the impacted end of the bar, shown in Fig. 13.

It was observed that using two strain measurements and one velocity measurement
reduces the disturbance due to the imperfections in much the same way as using
three strain measurements. Using a velocity measurement improves the correction of
the errors due to the imprecise knowledge of the zero strain, as might be expected.
Directly measuring the velocity in at least one point prevents an excessively large shift
of the reconstructed velocities from occurring at the ends of the bar.

It can be concluded that with an “*N*-strain, *P*-velocity” method, the imprecision of
the initial measurements is not amplified by the deconvolution procedure. It is therefore
possible to use this method successfully in real experimental environments.
**5. Experimental illustration of the multi-point method**

Two experimental tests were performed. The first one was of the kind simulated above, with an initial short loading pulse. The second test was performed using a long time loading pulse in order to demonstrate the validity of the method when the elementary waves overlap at the gauge stations right from the onset of the loading.

The set up used for this purpose is described in Fig. 14.

A 3.022-m long Nylon bar 40 mm in diameter was used. Strain gauge stations 1, 2 and 3 were 0.404, 1.508 and 2.414 m from the impact end of the bar, respectively.

Two special optical extensometers ("TSI Zimmer PM100H") were used to measure the displacement of an optical contrast painted on the bar 2.025 m (station 4) from the impact end and the displacement of the free end (B).

In the first test, the striker, which was 0.40 m long, was made of the same material and had the same diameter as the bar. The initial speed of the striker was 3.03 m/s.

In the second one, the striker was made of steel. It was 1.8 m long and had a diameter of 45 mm, which was greater than the diameter of the bars. Given that its impedance was much greater than that of the bar, the duration of the loading pulse was longer than a round trip along the bar and back (say $\Delta T$). The initial speed of the striker was 3.27 m/s.

The velocity at point 4 was obtained by performing numerical differentiation on the recorded displacement. The strains recorded at stations 1, 2 and 3 and the velocity at point A were used as the input data in formula (21) ($N(3)$-strains/$P(1)$ velocities).

The corresponding data are shown in Fig. 15 for the first test, and in Fig. 17 for the second one.

The dispersion relation used in both cases was previously established for the same Nylon bars in (Zhao and Gary, 1995). It is used here as the first mode solution of relation (19). The constitutive linear viscoelastic model was derived from a one-dimensional rheological model consisting of four Voigt elements and a spring connected in series (same as in Zhao and Gary, 1995). Assuming that we were dealing with a homogeneous isotropic material and a constant Poisson's ratio, the function $\xi(\omega)$ was constructed.

With these signals (Figs. 15 and 17), formula (21) was used to reconstruct the stresses and velocities at the free bar end. The technique was improved as follows.

Because of the numerical derivation used to measure the velocity at station 4, and also because the displacement measurement is very sensitive to the quality of the lighting on the optical target, the velocity showed a high frequency noise. From the simulations presented above, we knew that the velocity measurements served mainly to correct the effects of the imprecise zero strain measurements. These effects were known to result in a low frequency shift of the reconstructed velocities (see Fig. 12). The velocity data used in the calculation were therefore low pass filtered at 3 kHz. The corresponding velocity signal is shown in Figs. 15 and 17. In the same way, formula (21) was used only for the low frequency velocity and strain signals (below 3 kHz). The higher frequency components (above 3 kHz) of the basic waves were reconstructed using only strain signals and formula (20). All the frequency components were then combined before performing the inverse Fourier transform, which yields elementary waves in the time domain, and then strains and velocities at any point.
<figure><img src="image_18.png" /><figcaption>Fig. 16. Short loading pulse: N(3)-strain P(1)-velocity gauge method reconstructed stresses at the ends of the bar and displacement at the free end, based on experimental data. (The stress at the free end has been shifted down by 0.2 MPa to make the figure more readable.)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Long loading pulse, recorded data: strains, displacement at station 4 and corresponding filtered calculated velocity (3 kHz).</figcaption></figure>

The displacement of the free bar end was deduced from the corresponding calculated velocity and compared with the direct measurement (point B, see Fig. 14).

Both results are shown in Fig. 16 for the shorter loading pulse and in Fig. 18 for the longer loading pulse. The zero force expected at the bar end can be seen here to have actually occurred. In comparison with the force recorded at the impact end of the bar, the error in the reconstructed force measurement turns out to be less than 3.5%.

The reconstructed displacement was very similar to the directly measured value: the relative error in the reconstructed displacement was less than 2.5%.

The band pass of these measurements was greater than 20 kHz (based on a spectral analysis of the signals). With the short striker, the total measurement time was 50
<figure><img src="image_20.png" /><figcaption>Fig. 18. Long loading pulse: N(3)-strain P(1)-velocity gauge method reconstructed stresses at the ends of the bar and displacement at the free end, based on experimental data. (The stress at the free end has been shifted down by 1 MPa to make the figure more readable.)</figcaption></figure>

times longer than that available with the corresponding Hopkinson bar. With the long striker, the total displacement of the free end of the bar was 25 times longer than that of the corresponding Hopkinson bar with the same initial speed of loading.

These results can therefore be said to be highly satisfactory.

## 6. Conclusion

In the present paper, a multi-point method (multi-strain and/or multi-velocity measurements) is presented for reconstructing one-dimensional waves in bars. It is proposed to perform the wave separation in the frequency domain for all frequencies simultaneously, including those with which the denominator equals zero. This method is exact when used with the single-mode dispersive propagation model commonly applied to Hopkinson bars. The sensitivity of this method to imprecise measurements affecting the accuracy of the results increasingly with time is assessed.

This method yields consistent results (the inaccuracy due to imprecise measurements does not increase with time). It is illustrated here by applying it successfully to the analysis of a real test on a Nylon bar.

As mentioned by previous authors (Zhao and Gary, 1997; Bacon, 1999), the significant increase in the observation time available when using measuring techniques based on the use of bars such as SHPB set-ups makes it possible to obtain precise measurements at medium strain rates in a test range in between that of mechanical testing machines and that of Hopkinson bars.

## Acknowledgements

The authors benefited from many discussions with Dr. R.H. Blanc at Laboratoire de Mécanique et d'Acoustique du CNRS, BP71, 13402 Marseille cedex, F., Dr. H. Zhao
at Laboratoire de Mécanique et Technologie de Cachan, (University Pierre & Marie Curie - Paris 6), 61, Avenue du Président Wilson, 94235 Cachan cedex, F. and Dr. S.M. Walley at Cambridge Cavendish Laboratory, Cambridge, GB). They also want to thank the following former students at Ecole Polytechnique, who performed helpful simulations to test the new two-point method: Cedric Bourillet, Christelle Callamand, Max-André Delannoy and Alexandre Schaer.

This work has been partially supported by the "Conseil Régional de l'île de France"

References

Bacon, C., 1999. Separating waves propagating in an elastic or viscoelastic Hopkinson pressure bar with three-dimensional effects. Int. J. Impact Eng. 22, 55-69.

Bacon, C., Carlsson, J., Lataillade, J.L., 1991. Evaluation of force and particle velocity at the heated end of rod subjected to impact loading. J. Phys. IV France 1 Colloq. C3, 395-402.

Bacon, C., Färm, J., Lataillade, J.L., 1994. Dynamic fracture toughness determined from load-point displacement. Exp. Mech. 34 (9), 217-223.

Bacon, C., Brun, A., 2000. Methodology for a Hopkinson bar test with a non-uniform viscoelastic bar. Int. J. Impact Eng. 24, 219-230.

Bancroft, D., 1941. The velocity of longitudinal waves in cylindrical bars. Phys. Rev. 59, 588-593.

Bell, J.F., 1996. An experimental diffraction grating study of quasi-static hypothesis of SHPB experiments. J. Mech. Phys. Solids 14, 309-327.

Bertholf, L.D., Karnes, J., 1975. Two-dimensional analysis of the split Hopkinson pressure bar system. J. Mech. Phys. Solids 23, 1-19.

Campbell, J.D., Duby, J., 1956. The yield behaviour of mild steel in dynamic compression. Proc. Roy. Soc. London A 236, 24-40.

Chree, C., 1889. The equations of an isotropic elastic solid in polar and cylindrical co-ords, their solutions and applications. Cambridge Philos. Soc. Trans. 14, 250-369.

Conn, A.F., 1965. On the use of thin wafers to study dynamic properties of metals. J. Mech. Phys. Solids 13, 311-327.

Davies, R.M., 1948. A critical study of Hopkinson pressure bar. Philos. Trans. Roy. Soc. A 240, 375-457.

Davies, E.D.H., Hunter, S.C., 1963. The dynamic compression testing of solids by the method of the split Hopkinson pressure bar. J. Mech. Phys. Solids 11, 155-179.

Dharan, C.K.H., Hauser, F.E., 1970. Determination of stress-strain characteristics at very high strain rates. Exp. Mech. 10, 370-376.

Duffy, J., Campbell, J.D., Hawley, R.H., 1971. On the use of a torsional split Hopkinson bar to study rate effects in 1100-0 aluminium. J. Appl. Mech. 38, 83-91.

Follansbee, P.S., Franz, C., 1983. Wave propagation in the split Hopkinson pressure bar. J. Eng. Mater. Technol. 105, 61-66.

Gong, J.C., Malvern, L.E., Jenkins, D.A., 1990. Dispersion investigation in the split Hopkinson pressure bar. J. Eng. Mater. Technol. 112, 309-314.

Gorham, D.A., 1983. A numerical method for the correction of dispersion in pressure bar signals. J. Phys. E 16, 477-479.

Harding, J., Wood, E.D., Campbell, J.D., 1960. Tensile testing of materials at impact rate of strain. J. Mech. Eng. Sci. 2, 88-96.

Hillström, L., Mossberg, M., Lundberg, B., 2000. Identification of complex modulus from measured strains on an axially impacted bar using least squares. J. Sound Vib. 230 (3), 689-707.

Hopkinson, B., 1914. A method of measuring the pressure in the detonation of high explosives or by the impact of bullets. Philos. Trans. Roy. Soc. A 213, 437-452.

Hunter, S.C., 1960. Viscoelastic waves. In: Sneddon, I.N., Hill, R. (Eds.), Progress in Solid Mechanics, Vol. 1. North-Holland, Amsterdam, pp. 1-57.
Jahsman, W.E., 1971. Reexamination of the Kolsky technique for measuring dynamic material behaviour. J. Appl. Mech. 38, 77-82.

Klepaczko, J.R., 1969. Lateral inertia effects in the compression impact experiments. Reports of Inst. Fund. Technical Research, No. 17, Warsaw.

Kolsky, H., 1949. An investigation of the mechanical properties of materials at very high rates of loading. Proc. Phys. Soc. B 62, 676-700.

Kolsky, H., 1963. Stress Waves in Solids. Clarendon Press, Oxford.

Lataillade, J.L., Bacon, C., Collombet, F., Delaet, M., 1994. The benefit of Hopkinson bar techniques for the investigation of composite and ceramic materials. In: Kinra, W.K., Clifton, R.J., Johnson, G.C. (Eds.), Wave Propagation and Emerging Technologies, AMD-Vol. 188. ASME, New York, pp. 85-93.

Lifshitz, J.M., Leber, H., 1994. Data processing in the split Hopkinson pressure bar tests. Int. J. Impact Eng. 15, 723-733.

Lindholm, U.S., 1964. Some experiments with the split Hopkinson pressure bar. J. Mech. Phys. Solids 12, 317-335.

Lundberg, B., Blanc, R.H., 1988. Determination of mechanical material properties from the two-point response of an impacted linearly viscoelastic rod specimen. J. Sound Vib. 126, 97-108.

Lundberg, B., Carlsson, J., Sundin, K.G., 1990. Analysis of elastic waves in non-uniform rods from two-point strain measurement. J. Sound Vib. 137, 483-493.

Lundberg, B., Henchoz, A., 1977. Analysis of elastic waves from two-point strain measurement. Exp. Mech. 17, 213-218.

Morse, P.M., Feshbach, H., 1953. Methods of Theoretical Physics. McGraw-Hill, New York.

Park, S.W., Zhou, M., 1999. Separation of elastic waves in split Hopkinson bars using one-point strain measurements. Exp. Mech. 39, 287-294.

Pochhammer, L., 1876. Über die Fortpflanzungsgeschwindigkeiten kleiner Schwingungen in einem unbergrenzten isotropen Kreiszylinder. J. Reine Angew. Math. 81, 324-336.

Safford, N.A., 1992. Materials testing up to $10^5 s^{-1}$ using a miniaturised Hopkinson bar with dispersion corrections. In: Guanren, Z., Shihui, H. (Eds.), Proceedings of the Second International Symposium on Intense Dynamic Loading and its Effects. Sichuan University Press, China, pp. 378-383.

Van Trees, H.L., 1992. Detection, Estimation and Modulation Theory. Vols. I, III. Wiley, New York, 1968. Krieger, Malabar, 1992.

Yew, E.H., Chen, C.S., 1978. Experimental study of dispersive waves in beam and rod using FFT. J. Appl. Mech. 45, 940-942.

Zhao, H., Gary, G., 1994. Une nouvelle méthode de séparation des ondes pour l'analyse des essais dynamiques. C. R. Acad. Sci. Paris Sér II 319, 987-992.

Zhao, H., Gary, G., 1995. A three dimensional analytical solution of longitudinal wave propagation in an infinite linear viscoelastic cylindrical bar. Application to experimental techniques. J. Mech. Phys. Solids 43 (8), 1335-1348.

Zhao, H., Gary, G., 1996. On the use of SHPB techniques to determine the dynamic behaviour of materials in the range of small strains. Int. J. Solids Struct. 33 (23), 3363-3375.

Zhao, H., Gary, G., 1997. A new method for the separation of waves. Application to the SHPB technique for an unlimited measuring duration. J. Mech. Phys. Solids 45, 1185-1202.