Complex modulus estimation respecting causality: Application to viscoelastic bars

The identification of linear visco-elasticity models mostly focuses on the real and imaginary parts of the Young’s modulus. Many methods have been proposed in the past to identify these material model parameters from experiments. However, when these parameters are determined independently, they are likely to violate the principle of causality. The present work presents a method that accounts for the constraints of causality and positivity of dissipation rate. The proposed method is based on a finite set of n measured angular frequencies and complex moduli. It includes a noise reduction procedure and provides a rheological 2p + 1)-parameter model with p < gn which corresponds to a specific configuration of pairs of springs and dashpots. The poles of the complex modulus on the positive imaginary frequency axis are determined by p parameters which are obtained as the common positive zeros of a special class of rational functions, while the remaining parameters are obtained from a least squares fit. The level of refinement of the rheological model, expressed by p, is not an assumed value but a result of the method. The method is applied to an impact test with a Nylon bar specimen. In this case, data at the n = 29 lowest resonance frequencies resulted in a rheological model with 14 parameters (p = 6). The validity of the method is checked through supplementary experimental results at low frequencies.


Introduction
Methods for estimation of the properties of viscoelastic materials have been proposed and used since long.Such methods, devoted to estimation of complex-valued frequency-dependent material parameters such as the complex modulus, generally produce data from which discrete values of the parameters can be obtained for a set of frequencies.
For the measurement of viscoelastic properties, machines for dynamic mechanical analysis (DMA) are commercially available.A review of such machines can be found in [1].Most methods make use of the vibration of a particular structure [2][3][4][5].Wave propagation in bars has also been used to determine the viscoelastic properties of materials [6][7][8][9][10][11][12].
Often, it is not ensured that the real and imaginary parts of the complex modulus are consistent with the principle of causality.One way to achieve such consistency is to express them in terms of a rheological model constituted by linear springs and dashpots [5,13].An alternative way consists in making sure that the complex modulus is consistent with the Kramers-Kronig relations.Another fundamental constraint on the complex modulus comes from the second principle of thermodynamics: the dissipation per cycle must be positive for any harmonic load.Furthermore, the relaxation function, which is the inverse Fourier transform of the complex modulus, must be real.About these constraints, see, e.g.[14,15] In this paper we present a method for noise-corrected estimation of the complex modulus of a viscoelastic material, subject to the constraints of causality, positivity of dissipation and reality of relaxation function, given an experimental set of frequencies and corresponding values of the complex modulus.The estimation method provides the structure, the number and type of elements (springs and dashpots) and the parameters of a rheological model for the complex modulus.It will be illustrated with experimental data obtained by impacting a Nylon bar specimen.

Theory
According to Boltzmann's model of viscoelastic materials, and after Fourier transformation, the relation between stress σ (t) and strain ε (t) is given by where E (ω) is the complex modulus and ω is the angular frequency.
The complex modulus must satisfy the three constraints of causality, positivity of the dissipation rate and reality of the relaxation function.The problem to be considered is that of finding a noise-corrected complex modulus E (ω), subject to these constraints, given a finite set of angular frequencies and corresponding complex moduli obtained experimentally.
The constraint of causality implies the Kramers-Kronig relations that connect the real and imaginary parts of the complex modulus.If either the real or the imaginary part can be estimated accurately at a sufficient number of frequencies, these relations can be used to estimate the other part.With the experimental method of this paper, results for the complex modulus are obtained only for a relatively small set of frequencies.Therefore, we adopt EPJ Web of Conferences the approach of estimating simultaneously the real and imaginary parts of the complex modulus.
The constraints of causality and of positivity of the dissipation rate [14][15][16][17] are equivalent to requiring that the function E (ω) be analytic in the lower half plane and that its imaginary part be positive on the positive real axis, respectively.The constraint of reality of the relaxation function is equivalent to the requirement E (−ω) = E (ω).Here we will use the stronger constraint of complete monotonicity of the relaxation function [17,18].This constraint, which essentially means that the relaxation function and its derivatives must be monotone, is satisfied by most rheological models [17].To which extent it holds for a particular real material can be verified from experimental data.
The property of complete monotonicity of the relaxation function implies that this function is the Laplace transform of a positive function (the Bernstein-Widder Theorem, see [17]).This means that the complex modulus E (ω) is the sum of a linear function and a Stieltjes transform.More precisely, there exist two parameters α 0 > 0 and β 0 ≥ 0 a function h (s) ≥ 0, defined on the positive real axis such that Such functions (1) are analytic outside the positive imaginary axis of the complex plane and they satisfy the three constraints.Thus, they are analytic in the lower half plane and therefore satisfy the constraint of causality; their imaginary part is positive on the positive real axis and therefore the constraint of positiveness of the dissipation rate holds.Finally, the constraint of reality of the relaxation function is satisfied as E (−ω) = E (ω).
In order to check whether experimental data are of the form of Eq. ( 1), in particular for the purpose of determining the constants α 0 and β 0 , and the function h, it is convenient to apply the transformation f (z) = − E(iz) z that has the representation with where δ (s) is a delta function at the origin, and g (s) = 0 for s < 0. The function f (z) is analytic except on the positive real axis, and positive on the negative real axis.The problem of finding a function f (z) with the representation (2) given its values for finitely many z was studied in detail by Krein and Nudelman [19].Some of their results, when formulated for the complex modulus E (ω), are as follows.From the values E 1 , E 2 , . . ., E n of the function E (ω) at n angular frequencies ω 1 , ω 2 , . . ., ω n one can form two matrices with elements It can be verified that the matrix M (1) is positive semidefinite.The same properties hold for M (2) since E (0) > 0.
If one of the matrices M (1) or M (2) has a zero eigenvalue, the function E (ω) is completely determined and is a rational fraction (see [19] and below).Assume now that a is an eigenvector corresponding to a zero eigenvalue of M (1) .Then, since M (1) a = 0, we also have a, M (1) Since both terms in the left member are non-negative, they must vanish.In particular, the function h (s) must be a linear combination p l=1 γ l δ (s − s l ) with coefficients γ l ≥ 0 of delta functions δ (s − s l ) at the positive zeros s 1 , . . ., s p (p < n) of the rational function For the complex modulus given by Eq. ( 1) this implies the model where α 0 > 0 and β 0 , γ 1 , γ 2 , . . ., γ p are non-negative.The same considerations hold for the matrix M (2) and should lead to the same model.If b is an eigenvector corresponding to a zero eigenvalue of M (2) , we deduce as above that the function h (s) must also be a linear combination with non-negative coefficients of delta functions at the positive zeros of the rational function In particular, we only need to consider the common positive zeros of the rational functions ( 5) and (7).Moreover, if there are several linearly independent eigenvectors corresponding to a zero eigenvalue, we get a set of positive zeros for each of them and should keep only the ones which are common.
The model ( 6) for the complex modulus corresponds to the rheological model with 2 (p + 1) parameters shown in Fig. 1.

Implementation
From experiments, we obtain a finite set of angular frequencies ω 1 , . . ., ω n and corresponding complex moduli E 1 , . . .E n .By use of these experimental data we form the two matrices M (1) j,k and M (1)  j,k given by the first equalities of Eqs. ( 3) and (4), respectively.These self-adjoint matrices are positive semi-definite if and only if all their eigenvalues are non-negative.However, experiments commonly result in some slightly negative eigenvalues.In the experimental part of this paper, for example, the modulus of the smallest negative eigenvalue is typically a few tenths of a percent of the largest positive eigenvalue.A likely explanation for such small violations of the positive semi-definiteness is experimental noise.As advocated in [20] for a similar problem, the effects of such noise can be reduced by searching the smallest possible corrections of angular frequencies and moduli which restore the non-negativeness of all eigenvalues and establishing noise-corrected sets ω corr 1 , ..., ω corr n and E corr 1 , . . ., E corr n .If at least one of the noise-corrected matrices M (1)corr and M (2)corr has a zero eigenvector, one forms the associated rational fraction (5) or (7) and finds its positive zeros.This is repeated for the set of independent eigenvectors, if any, corresponding to zero eigenvalues of both matrices.Then one keeps the common positive zeros s 1 , . . ., s p of all these rational fractions and obtains the model given by Eq. ( 6).In this model the p positive parameters s 1 , . . ., s p are known, while α 0 > 0, β 0 ≥ 0 and the p nonnegative parameters γ 1 , γ 2 , . . ., γ p are identified by a least square fit of the model to the noise-corrected data.In this way, the method provides the structure, the number of elements (springs and dashpots) and the parameters of the rheological model.Normally, the number 2(p+1) of elements is found to be relatively low.
The experimental set of angular frequencies and corresponding values of the complex modulus may be obtained with a variety of methods, quasi-static as well as dynamic.
In what follows they will be obtained by analysing the wave propagation in a Nylon bar specimen loaded through axial impact.
3 Impact test method

Theory
Here, discrete values of the complex modulus E (ω) = E (ω) + iE (ω) will be obtained from the resonances of an impacted uniform bar specimen with density ρ and length l as shown in Fig. 2. At one end, x = 0, the bar is impacted axially by a striker that separates from the bar after the generation of a compressive primary pulse in the bar.The other end of the bar, x = l, is free.The strain ε b (t) = ε (b, t) is recorded at a distance a from the free end and b from the impacted end (a + b = l).The primary pulse should be shorter than 2a so that there is no overlap in the measured strain of this pulse and the first pulse reflected from the free end.However, it should be much longer than the diameter of the bar so that approximate 1D conditions prevail (wavelengths much longer than the diameter of the bar [21].
In the frequency domain, the strain in the bar can be expressed as where Here, ε (x, ω) is the Fourier transform of ε (x, t), and A (ω) and B (ω) are complex amplitudes of waves travelling in the directions of increasing and decreasing x, respectively, k (ω) is the wave number and α (ω) is the damping coefficient.
In order to determine an expression for the recorded strain ε1 b (ω) = ε (b, ω) associated with the primary pulse alone, the amplitudes A and B are first determined from Eq. ( 9) and the boundary conditions ε (0, ω) = ε0 (ω) and B (ω) = 0 for a semi-infinite bar x ≥ 0, where ε0 (ω) is the strain at the impacted end.With these amplitudes and x = b inserted, Eq. ( 9) give The recorded strain ε1 b (ω) = ε (b, ω) associated with the complete train of pulses is determined similarly.First, the amplitudes A and B are determined from Eq. ( 9) and the boundary conditions ε (0, ω) = ε0 (ω) and ε (l, ω) = 0 for the finite bar 0 ≤ x ≤ l.With these amplitudes, and x = b and l − b = a inserted, Eq. ( 9) gives Dividing the members of Eq. ( 12) by those of Eq. ( 11) eliminates the strain ε0 (ω) at the impacted end which normally cannot be measured.The elimination of this strain makes the difference between the method used here and that used by Othman et al. [22].Substituting ξ from the second of Eqs. ( 10 Resonance occurs at the angular frequencies ω = ω m , m =1, 2, . . .which correspond to the wave numbers, wavelengths and phase velocities respectively. It is assumed that within the m:th resonance peak α = α m ω/ω m can be taken as directly proportional to angular frequency and c (ω) = ω/k (ω) = c m can be taken as constant.By use of the third of Eqs. ( 14) we then obtain the relation k (ω) = k m ω/ω m between wave number and angular frequency within the resonance peak.Inserting these expressions for α (ω) and k (ω) into Eq.( 13) we get within the m:th resonance peak.In this relation, the spectra ε∞ b (ω) 2 and ε1 b (ω) 2 can be determined experimentally.
For each resonance m = 1, 2, . .., the resonance frequency ω m and damping coefficient α m can be estimated by minimizing the difference between the resonance peaks represented by the left and right members of Eq. (15).From ω m and m, the complex modulus E m at each resonance frequency can be finally obtained as

Experimental set-up and procedure
Impact tests were carried out with a Nylon bar specimen of length 3045 mm, diameter 10.2 mm, and density 1149 kg/m.Two pairs of strain gauges, one axial and one circumferential, were located at a distance of 1731 mm from the free end and 1314 mm from the impacted end.
In order to minimise the effects of friction, the bar was suspended horizontally by means of nine regularly spaced  strings with length 300 mm.With this arrangement, the period of oscillation of the system was about one second, which is much longer than the test duration.The striker had length 174 mm and the same diameter 10.2 mm and material as the bar, and its impact velocity was 3.7 m/s.The strain was recorded with 500 kHz sampling frequency, and the signal vanished completely before the end of the recording window.
Results were evaluated up to 8 kHz.With a phase velocity expected to be higher than 1500 m/s at this frequency, the shortest significant wave length was estimated to be greater than 0.19 m.This wave length is much larger than the diameter 0.01 m of the bar, which means that 3D effects can be neglected as assumed in Section 3.1.

Results and discussion
The estimation method presented provides a model for the complex modulus corresponding to that shown in Fig. 1, with p +1 pairs of springs and dashpots.Consequently, p represents the complexity of the material behaviour.It is not an assumed value but a result of the method.
The recorded strain is shown in Figs. 3 and 4. Figure 3 shows a long-time record, and Fig. 4 shows the primary compressive pulse followed by one tensile and one compressive pulse which have undergone one and two free-end reflections, respectively.As there is no significant overlap of these pulses, it was possible to compute the spectrum ε1 b 2 in the right member of Eq. (30) from the measured primary pulse.
The angular frequencies at resonance and the corresponding complex moduli were determined as described.Figure 5 shows that there was good agreement between 01012-p.4  the resonance peaks determined from the left and right members of Eq. (30).For the material of the Nylon bar specimen tested, a rheological model with p = 6 was obtained for n = 29.
Figures 6 and 7 show the frequency dependency of the real and imaginary parts of the complex modulus, viz., discrete experimental results, discrete noise-corrected results, and continuous results for the rheological model given by Eqs. ( 6) and ( 8) with the above sets of parameters.
From the discrete experimental results for the phase velocity and the damping coefficient, the discrete experimental results for the complex modulus (open circles) were obtained as described in Section 3. The corresponding discrete noise-corrected results (filled circles) and the continuous results for the rheological model obtained as described in Section 2.1 are shown in the same figures.
While the experimental results are scattered due to noise, there is close agreement between the noise-corrected results and the smooth curves representing the rheological model.This indicates that the noise-correction was effective and that the rheological model obtained with p = 6 provides an adequate representation of the complex modulus in the frequency range of the test.
At low frequencies, much below the range of the test, the rheological model cannot be expected to give a quantitatively accurate representation of the complex modulus.However, it is still interesting to note that at such frequencies the model predicts a rapid increase with frequency of the real part and a narrow maximum of the imaginary part of the complex modulus.Qualitatively, the observed low-frequency behaviour is consistent with the results at low frequency of servo-hydraulic tests which were carried out subsequently.

Conclusion
The estimation method of this paper has added some improvements to the identification of rheological models of the type shown in Fig. 1.First, the procedure for noise reduction is fully integrated in the method.Secondly, the method provides a rheological model with a number of elements that is in accordance with the complexity of the material.

Fig. 3 .
Fig. 3. Recorded strain in the Nylon bar specimen.(a) Long-time record of strain pulses.

Fig. 4 .
Fig. 4. Primary strain pulse followed by strain pulses which have undergone one and two free-end reflections.

DYMAT 2012 Fig. 5 .
Fig. 5. Details of the spectrum.The thin and thick curves are based on the left and the right member of Eq. (30), respectively.(upper) 1st resonance peak (lower) 25th resonance peak.

Fig. 6 .
Fig. 6.Real part of the complex modulus versus frequency.