- Email: [email protected]
Identification of stiffness and damping properties of plates by using the local equation of motion
Journal of Sound and Vibration 333 (2014) 2454–2468
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Identification of stiffness and damping properties of plates by using the local equation of motion Frédéric Ablitzer a,n, Charles Pézerat a, Jean-Michel Génevaux a, Jérôme Bégué b a b
LUNAM Université, Université du Maine, CNRS UMR 6613, LAUM, Avenue Olivier Messiaen, 72085 Le Mans Cedex 9, France CETIM (Centre Technique des Industries Mécaniques), Technocampus EMC2, Chemin du Chaffault, 44340 Bouguenais, France
a r t i c l e in f o
abstract
Article history: Received 15 January 2013 Received in revised form 12 December 2013 Accepted 13 December 2013 Handling Editor: I. Trendafilova Available online 23 January 2014
This paper deals with the identification of stiffness and damping properties of vibrating structures by an inverse method inspired from the Force Analysis Technique (FAT). The proposed approach uses a local equation of motion assumed a priori, which provides a relative straightforward relationship between the displacement field and material properties. The spatial derivatives of the displacement in the equation are calculated using finite differences. As this operation amplifies measurement noise, a regularization step is applied before solving the inverse problem. A procedure is proposed to automatically adjust the level of regularization. The method also allows one to identify local stiffness and damping on a heterogeneous structure. Illustrations for both homogeneous and heterogeneous cases are shown using simulated and measured displacement fields. & 2013 Elsevier Ltd. All rights reserved.
1. Introduction Composite materials, which provide high stiffness and light weight, are increasingly used in the aerospace and automotive industry. During the conception of structures including such materials, an extensive use of finite element models is made to predict the dynamic and vibroacoustic behavior. To achieve predictions that are as accurate as possible, one of the key challenges is to enter the right material characteristics into the model. As a wide variety of composite materials exist, their elastic and damping characteristics are rarely available in tables and must be obtained experimentally. Moreover, if the finite element method uses plate elements, the determination of the characteristics of the equivalent plate (stiffness and damping for all types of loads) needs the choice of kinematic assumptions. This choice influences the results of the model and has to be confirmed by measurements. Measurement techniques based on modal analysis theory, which is well established [1], are widely used. The measurement of modal parameters (natural frequencies, mode shapes and damping ratio) on a test structure, generally combined to a numerical model, allows one to identify elastic and damping properties of the constitutive material [2]. However, there are several limitations to such methods. First, as it is a global approach, it is not possible to characterize the material directly on the target structure itself: experiments are generally carried out on specific test specimens, a precise knowledge of boundary conditions being crucial. Secondly, modal parameters are difficult to measure in the mid-frequency domain, characterized by a strong modal overlap. Because the stiffness and damping of composites may strongly vary with frequency, it is necessary to obtain material properties in a wide frequency range instead of extrapolating results obtained at low frequencies. In recent years, works have been done on the estimation of wavenumbers in a structure, from which material properties such as loss factor or dynamic stiffness may be deduced [3–6]. Such wave-based approaches overcome most limitations
n
Corresponding author. E-mail address: [email protected] (F. Ablitzer).
0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.12.013
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2455
of the modal approach. In particular, they are applicable at any frequency and do not require any knowledge about boundary conditions. Their principle is to find the wavenumber giving the best fit between the displacement field measured at a given driving frequency and a theoretical field consisting of a sum of waves. This theoretical field can be viewed as the solution of an equation of motion involving material properties, which may be related to the identified wavenumber. This paper presents a new method to identify material properties from the displacement field measured locally on a structure. The underlying idea is to use the local equation of motion itself, which provides a relative straightforward relationship between the material properties and the displacement field. The method is inspired from the Force Analysis Technique (FAT, also known by the french acronym RIFF), initially aimed at identifying dynamic loads acting on structures [7–11]. In the Force Analysis Technique, the measured displacement field is injected into a discrete form of the local equation of motion, in which the spatial derivatives are approximated by a finite-difference scheme, leading to a straightforward calculation of the force acting in the area considered. This simple principle, combined with an appropriate regularization procedure (windowing and filtering), has been successfully applied to beams [7], plates [8,9], shells [12,13], and thereafter extended to more complex structures by using a finite element operator instead of an analytical equation of motion [14]. To localize and quantify external forces, the Force Analysis Technique obviously requires the a priori knowledge of material properties involved in the equation of motion. The originality of the present method is to identify unknown material properties from the equation of motion by considering an area of the structure which is free of any external force. The key ingredients for the identification of the material properties (approximation of partial derivatives, regularization) are the same as those used in FAT, complemented with specific developments to compensate the lack of a priori knowledge about the structure. In particular, a novel automatic regularization technique is proposed. The paper focuses on the application to plates. However, the method may be extended to other types of structures for which a differential equation exists. The continuation of this paper is structured as follows. First, the general principle of the method is exposed more in detail. Then, the feasibility of the proposed approach is assessed using simulated displacement fields. The effect of noise is illustrated and a procedure to automatically adjust the regularization is described. The applicability of the method is assessed for both homogeneous and heterogeneous structures. The last part presents experimental results. 2. Principle of the identification technique The equation of motion of an isotropic thin plate in the harmonic regime is [15] 4 ∂ w ∂4 w ∂4 w D ρhω2 wðx; yÞ ¼ f ðx; yÞ; þ 2 þ ∂x4 ∂x2 ∂y2 ∂y4
(1)
where D is the flexural stiffness, ρ the density, h the thickness, ω the angular frequency, wðx; yÞ the transverse displacement and f ðx; yÞ the force per unit area. The material is assumed to be linear viscoelastic. Therefore, the flexural stiffness 3
D¼
Eð1 þ jηÞh 12ð1 ν2 Þ
(2)
involves a complex Young0 s modulus Eð1 þ jηÞ, where j is the unit imaginary number and η denotes the loss factor, which characterizes material damping. The equation of motion (1) is local and therefore valid everywhere in the structure. More particularly, considering an area on which no external force is applied, i.e., f ðx; yÞ ¼ 0, Eq. (1) becomes 4 D ∂ w ∂4 w ∂4 w þ2 þ ¼ wðx; yÞ: (3) ρhω2 ∂x4 ∂x2 ∂y2 ∂y4 Eq. (3) simply expresses the fact that the bilaplacian ∇4 wðx; yÞ equals the displacement wðx; yÞ up to a constant multiplier. The constant multiplier D=ρhω2 depends on the characteristics of the plate, and on the angular frequency ω. It should be noted that the quantity D=ρh is complex and may vary with frequency. It contains information about the stiffness and damping of the material. Thus, assuming that the displacement wðx; yÞ is known, it seems possible to determine the quantity D=ρh from Eq. (3), at any frequency, without any knowledge about the exact location and amplitude of the external force, or about the boundary conditions. The only requirement is that no external force is applied at the considered location (x,y). In practice, however, there are several obstacles to this simple idea. The first one is due to the presence of fourth-order partial derivatives in Eq. (3). Whereas the transverse displacement at a given location may be easily obtained by using an accelerometer or a laser vibrometer, its partial derivatives are less straightforward to obtain. To overcome this difficulty, one may measure the displacement at discrete abscissas ðxi ; yi Þ over a regular meshgrid, so that the partial derivatives may be approximated by a finite difference scheme. A discretized form of Eq. (3) is D 4x ¼ wi;j : δi;j þ 2δ2x2y þ δ4y (4) i;j i;j 2 ρhω 4y 2x2y The expressions of δ4x may be found in Appendix A. The approximation of the partial derivatives introduces i;j , δi;j and 2δi;j a bias which is likely to affect the value of D=ρh calculated from Eq. (4). Note that the corrected finite differences scheme proposed in [11] can also be applied.
2456
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
A second difficulty is due to the noise which corrupts the measured displacements. Therefore, the quantity D=ρh identified by the method is likely to vary depending on the point considered. To overcome this difficulty, one may use a set of N data points and find a unique value of D=ρh that satisfy Eq. (4) in a least-square sense, i.e., 2 3 2 3 w1 δ4x þ 2δ2x2y þ δ4y 1 1 6 1 7 D 6 7 6 7 ⋮ (5) 4 5 ρhω2 ¼ 4 ⋮ 5: 2x2y 4y 4x w δN þ 2δN þ δN N However, identifying material properties through the least squares solution of Eq. (5) is not sufficient to reduce the effect of measurement noise. As will be shown in the next section, the effect of noise which is present in the measured displacements wi;j is considerably amplified by the fourth-order partial derivatives, making this simple identification technique not applicable. With an appropriate denoising treatment, it is yet possible to overcome this obstacle, as shown in the following section. 3. Proof-of-concept using simulations 3.1. Homogeneous structures 3.1.1. Direct problem We consider a simply supported isotropic thin plate (Fig. 1), whose characteristics are given in Table 1. The plate is driven at angular frequency ω by a harmonic point force at location (x0 ¼ 0:05 m, y0 ¼ 0:05 m) with amplitude F 0 ¼ 1 N. The displacement field is calculated by modal expansion, i.e., M
wexact ðx; yÞ ¼ ∑
N
∑
4F 0
m ¼ 1 n ¼ 1 Lx Ly ρh
sin ðmπx=Lx Þ sin ðnπy=Ly Þ sin ðmπx0 =Lx Þ sin ðnπy0 =Ly Þ ; ω2mn ω2
(6)
pffiffiffiffiffiffiffiffiffiffiffi where ωmn ¼ D=ρhððmπ=Lx Þ2 þ ðnπ=Ly Þ2 Þ. The number of modes M and N is chosen so as to satisfy a spatial Shannon0 s criterion. The displacements are calculated at discrete abscissas over a Mx M y ¼ 51 51 meshgrid, delimited by x A ½0:2 m; 0:45 m and y A ½0:1 m; 0:35 m. This area does not contain the excitation point, which is a requirement to apply the proposed identification technique. The corresponding spatial step is Δx ¼ Δy ¼ 5 mm. In order to simulate measurement uncertainties, noise can be added to the displacements. The noisy displacements at every frequency are obtained as follows: wnoisy
i;j
¼ wexact
jΔφ i;j þ ΔA i;j e
i;j
;
(7)
where ΔA i;j is a Gaussian random value with mean zero and standard deviation schosen in order to reach a desired signalto-noise ratio (SNR), and Δφ i;j is a uniform random value in the range ½0; 2π. The signal-to-noise ratio is calculated as follows: SNR ¼ 20 log10
∑i;j jwexact i;j j2 ∑i;j jwnoisy
2 i;j wexact i;j j
:
(8)
wi;j without subscript refers to noisy displacements, whereas subscript “exact” is kept when referring to exact displacements. 3.1.2. Inverse problem The aim of the inverse problem is to recover some of the plate properties from the available displacement field. First, the feasibility of the proposed approach in the case of non-noisy data is shown. Then, we illustrate the problem arising from noisy measurements and propose a denoising treatment allowing the identification of properties. The moduli of the discrete field wexact i;j calculated by solving the direct problem, as well as the discrete approximation of 2x2y the bilaplacian operator ∇4i;j ¼ ðδ4x þ δ4y i;j þ 2δi;j i;j Þ, are represented in Fig. 2(a) and (b). A very good fit between both fields can be visually observed, yet with different colorscales, as expected from Eq. (4) which implies that jD=ρhω2 jjδ4x i;j þ 2δ2x2y þδ4y i;j i;j j ¼ jwi;j j. This linear relationship between the moduli may also be observed in Fig. 2(c), in which the exact displacements and the bilaplacian are plotted against each other for all points in the considered domain. At this particular
Table 1 Geometrical and material properties set for the direct problem. Length Width Thickness Density Poisson0 s ratio Young0 s modulus Loss factor
Lx ¼ 0.5 m Ly ¼ 0.4 m h ¼ 5 mm ρ ¼1500 kg m 3 ν ¼ 0:3 E ¼5 GPa η ¼ 0:05
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2457
y force
x
observation area
Fig. 1. Geometry of plate used for numerical simulation.
2x2y
×10−6
4y
0.35
0.35
6
0.1 0.2
0.45
4
yi (m)
yi (m)
1
0.5
2 0.1 0.2
0
xi (m)
wi, j (m)
×10−6 1.5
−3 δ4x i, j + 2δi, j + δi, j (m )
wexact i, j (m)
wexact i, j (m)
1 0.5 0
0.45 xi (m) 2x2y
×10−6
4y
×10−6 1.5
−3 δ4x i, j + 2δi, j + δi, j (m )
0.35
0.35
yi (m)
yi (m)
600 400
wi, j (m)
800
1
0.5
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
1 0.5
200 0.1 0.2
0.45 xi (m)
0
0.1 0.2
0.45 xi (m)
0
0
0 500 1000 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
Fig. 2. Colormaps of the moduli of the two terms in Eq. (5), which are also plotted against each other, at f¼ 1000 Hz. (a)–(c) without noise; (d)–(f) with noise on wi;j (SNR ¼40 dB). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
frequency, knowing ρ, h and ν and solving Eq. (5) we can calculate Young0 s modulus E ¼5.03 GPa and the loss factor η ¼ 4:99 10 2 . The identified values are very close to the actual ones (Table 1). The moduli of the same discrete fields are then plotted in the case of noisy data, with an SNR of 40 dB. Whereas the noisy displacement field wi;j (Fig. 2(d)) is not distinguishable from the exact field wexact i;j (Fig. 2(a)), the discrete bilaplacian ∇4i;j calculated from noisy displacements exhibits a high randomness. Not only the amplitude is considerably amplified compared to the non-noisy case (the range in values is scaled by about 130), but the distribution of energy in the wavenumber domain is also affected. The bilaplacian obtained from non-noisy displacements (Fig. 2(b)) exhibits a relative high spatial wavelength compared to the spatial step. This wavelength corresponds to the natural flexural wavelength of the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi structure λf ¼ 2π=kf at the considered frequency (kf ¼ 4 ðρh=DÞω2 is the natural wavenumber). The bilaplacian calculated now from noisy displacements (Fig. 2(e)) contains wavelengths in the order of the spatial step, as inferred from the huge variations in amplitude among neighboring points. As observed in Fig. 2(f), the linear relationship between the moduli of the displacements and the bilaplacian is not verified. The identification of Young0 s modulus E and the loss factor η in a large frequency range, f A ½12:5; 5000 Hz, is then performed. At each frequency, the displacement field wi;j is obtained by solving the direct problem. Then, the discrete bilaplacian ∇4i;j is calculated and Eq. (5) is solved in a least squares sense. Results obtained from the displacements without
2458
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
10
0.1 from exact displacements from noisy displacements
0.08
6
0.06 η
E (GPa)
8
4
0.04
2
0.02
0
0
1000
2000 3000 f (Hz)
4000
5000
0
0
1000
2000 3000 f (Hz)
4000
5000
Fig. 3. (a) Young0 s modulus E and (b) loss factor η identified. The arrows on both sides indicate the actual values of E and η.
and with noise (SNR ¼40 dB) are shown in Fig. 3. When applying the method to non-noisy data, both elastic and damping properties can be identified with good accuracy. The small deviation from the actual values, observed especially on Young0 s modulus, is due to the bias error introduced by the approximation of partial derivatives in the equation of motion. It should be reminded that no approximation regarding the derivatives in the operator is used for the direct problem, which is solved by modal decomposition. When using noisy displacements, the identification of Young0 s modulus is no longer reliable. A very low value is found 2x2y over the considered frequency range (Fig. 3(a)). This is due to the fact that the term δ4x þ δ4y i;j þ 2δi;j i;j in Eq. (4) is very high, because of the amplification of noise, which has to be counterbalanced by a small value of stiffness D. The identification of the loss factor using noisy data seems to be compromised only at low frequencies (Fig. 3(b)). When frequency increases, the values identified are closer to the actual ones. This result is the consequence of the finite difference scheme, which acts as a low-pass filter in the wavenumber domain when the wavelength becomes small compared to the spatial step [11]. However, this effect is not enough to ensure the robustness of the method, as it does not allow a good identification of stiffness. Moreover, the bias error due to the finite difference scheme increases as the wavelength decreases compared to the spatial step. The denoising treatment proposed in the next section is intended to keep the accuracy allowed by a fine spatial discretization, without suffering the consequences of noise amplification. 3.1.3. Regularization: windowing and filtering As seen in the previous section, the noise in the displacement field wi;j is considerably amplified by the finite difference scheme. The calculated discrete bilaplacian ∇4i;j is therefore corrupted by noise at high wavenumbers, impairing the identification of the elastic and damping properties. In order to eliminate high wavenumbers components, a low-pass filter in the wavenumber domain is applied to the displacement field wi;j and the bilaplacian ∇4i;j . For this purpose, a similar procedure as the one applied in the Force Analysis Technique is followed. First, each field is windowed by a bidimensional Tukey window ψ 2D i;j . This preliminary step replaces the truncation of the field at the edges by a smooth variation from zero amplitude. Then, the windowed field is convolved by the finite spatial response hi;j of a low-pass filter with cut-off wavenumber kc. The expression of the window and the filter response are given in Appendix B. In the Force Analysis Technique, the cut-off wavenumber is generally chosen proportional to the natural wavenumber kf, i.e., kc ¼ akf . The coefficient a is usually chosen between 1 and 4, depending on the signal-to-noise ratio: the higher the SNR, the higher the cut-off wavenumber. This way of choosing the level of filtering is allowed by the fact that the properties of the plate and therefore its natural wavenumber are a priori known. We propose a method able to determine some of these properties which are unknown. Therefore, a new criterion to select the cut-off wavenumber has to be determined. The effect of cut-off wavenumber on the windowed filtered fields is illustrated in Fig. 4. The four values of the considered cut-off wavenumber kc correspond to a ¼4, 3, 2 and 1. As kc diminishes, each field contains less and less components with small wavelengths. Whereas the effect is hardly distinguishable on the displacement field, it appears clearly on the bilaplacian field. At kc ¼ 4kf (Fig. 4(a)), the filter keeps too much information at high wavenumbers to make any deterministic pattern observable in j∇4i;j j. At kc ¼ 3kf (Fig. 4(b)), the pattern looks like the one of jwi;j j but is still corrupted by small wavelengths. A good agreement between the displacement field and the bilaplacian is obtained at kc ¼ 2kf (Fig. 4(c)). The filter is effective enough. At kc ¼kf (Fig. 4(d)), the filter brings no more significant noise reduction, but continues to remove energy to both fields, as pointed out by the change in colorscale of wi;j . The effect of the filter is also well observable when the displacements and the bilaplacian are plotted against each other. As kc decreases, less and less values with high amplitude are kept in j∇4i;j j. Whereas the scale for jwi;j j is the same as in Fig. 2(f) (no filter), the scale for j∇4i;j j is diminished by a factor of 100. At the same time the dispersion of data points is reduced, especially for j∇4i;j j. At kc ¼ 2kf (Fig. 4(c)), the same linear relationship as in Fig. 2(c) appears, yet with a small dispersion and a reduction in amplitude along both axes due to windowing and filtering. Diminishing further kc removes energy to both fields (Fig. 4(d)). At this point, it appears that the cut-off wavenumber kc should be chosen low enough to diminish the noise level in ∇4i;j , but not too low in order to keep as much information as possible. A value satisfying this qualitative criterion may be simply
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
8
4 2
2
wi, j (m)
0.45
0.1 0.2
×10−7
wi, j (m)
xi (m)
wi, j (m)
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
0.1 0.2
×10−6
0.35
xi (m) wi, j (m)
0.45
0
×10−7
0.35
1
yi (m)
0.5
4
0.5
2 1
0.1 0.2
xi (m)
0.45
0
0.1 0.2
×10−6 1.5 wi, j (m)
0.5
0.45
0
3
0.45
1
xi (m)
1
×10−6 1.5 wi, j (m)
0.1 0.2
×10−6
yi (m)
yi (m)
10 8 6 4 2
0
0.45
1 0.5
0
0.35
0.35
0.1 0.2
xi (m)
2
2 1.5
yi (m)
xi (m)
0
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
4
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
xi (m)
0.45
0
0.1 0.2
×10−6 1.5 wi, j (m)
0.1 0.2
4y
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
xi (m)
0.45
0
×10−6 1.5 wi, j (m)
4
yi (m)
yi (m)
6
2x2y
−3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
6
yi (m)
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
yi (m)
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
2459
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
Fig. 4. Effect of cut-off wavenumber kc on displacement wi;j and discrete bilaplacian ∇4i;j represented as colormaps and plotted against each other, in modulus. (a) kc ¼ 4kf , (b) kc ¼ 3kf , (c) kc ¼ 2kf , and (d) kc ¼kf. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
found by varying kc and observing plots like those represented in Fig. 4. For example, for the frequency and noise level considered in this figure, a value close to kc ¼ 2kf would be appropriate. The corresponding Young0 s modulus and loss factor identified are E¼5.03 GPa and η ¼ 5:06 10 2 , respectively. These values are very close to those found in the absence of noise. This shows the efficiency of regularization for this simulated noisy displacement field. 3.1.4. Automatic determination of optimum cut-off wavenumber The range of appropriate values of kc depends on the driving frequency and noise level, according to the Force Analysis Technique. The task of manually searching a value of kc for each frequency in a wideband measurement would be tedious. Therefore, it appears necessary to define a calculable indicator that expresses the qualitative criterion described above. To this aim, the solution of Eq. (5) is now focused on. Young0 s modulus (Fig. 5(a)) and the loss factor (Fig. 5(b)) are identified from the least squares solution of Eq. (5) for different values of cut-off wavenumber. The axis corresponding to kc is shown in a logarithmic scale, as the value for which the regularized solution converges to the non-regularized one (kc =kf ≳20) is very high compared to the range of interest. The identified Young0 s modulus is very small for high values of cut-off wavenumber (kc =kf ≳6). On the contrary, it is very close to the actual Young0 s modulus for low values of cut-off wavenumber (kc =kf ≲3), even when kc is significantly lower than the flexural wavenumber of the plate kf. This result may seem strange, as the energy of the displacement field in the wavenumber domain is essentially contained in the vicinity of kf. However, it can be explained by the fact that Eq. (3) remains valid when applying a linear filter to its terms, as the quantity D=ρhω2 does not depend on the spatial location. Thus, as long as the residual noise (i.e., with components lower that kc in the wavenumber domain) is not high enough to impair the identification, there is no reason that Young0 s modulus found should significantly differ from the actual one. The identified loss factor is also far from the actual value when the cut-off wavenumber is high. However, as the cut-off wavenumber decreases, the identification becomes satisfactory slightly earlier (kc =kf ≲5) than in the case of Young0 s modulus. This could be due to the fact that the imaginary part of the displacement field is much smaller than the real part, giving less weight to the amplification by the finite differences scheme. Thus, the filter needs not to be as restrictive as the one for Young0 s modulus identification. However, as the identification of both Young0 s modulus and structural loss factor is generally desired, the acceptable range in cut-off wavenumber has to be common to both quantities. 2 My 4 2 x The sum of squared residuals associated to the least squares solution, R ¼ ∑M i ¼ 0 ∑j ¼ 0 9∇i;j ðD=ρhω Þ wi;j 9 , is represented in Fig. 5(c). As expected, the residual is maximum when the cut-off wavenumber kc is very high, due to the fact that much noise is kept, especially in ∇4i;j . On the contrary, it is minimum at low values of kc, indicating a very good fit. It can be seen
2460
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
0.1
8
0.08
6
0.06 η
E (GPa)
10
4
0.04
2
0.02
0 10−1
100
101
0 10−1
102
100
kc / k f
102
101
102
kc / k f
8 ∂R/ ∂kc (normalized)
1
6 R (m2 )
101
4 2 0 10−1
100
101
0.5
threshold 0 10−1
102
100
kc / k f
kc / k f
Fig. 5. Effect of the level of regularization (cut-off wavenumber kc) on the results, (a) Young0 s modulus E, (b) loss factor η, (c) residual term R and (d) its derivative. The circles in (a) and (b) point out the values of E and η at the optimum cut-off wavenumber, for which the slope of residual reaches 5% of its maximum value.
10
0.1 SNR = 40 dB SNR = 30 dB SNR = 20 dB
0.08
6
0.06 η
E (GPa)
8
4
0.04
2
0.02
0
0 0
1000
2000 3000 f (Hz)
4000
5000
0
1000
2000 3000 f (Hz)
4000
5000
Fig. 6. (a) Young0 s modulus E and (b) loss factor η identified from noisy displacements, with a denoising treatment (windowing and filtering). Three levels of noise are considered. The arrows on both sides indicate the actual values of E and η.
that the value of cut-off wavenumber, at which the residual starts to rise significantly, kc =kf 3, corresponds to the value at which the identified Young0 s modulus starts to deviate from the actual Young0 s modulus (see Fig. 5(a)). The sudden rise in residual thus appears as a good detection indicator of the value of cut-off wavenumber at which the filter becomes efficient enough. This sudden rise may be easily found by considering the slope of the residual curve, ∂R=∂kc , which is obtained by discrete differentiation. The evolution of the slope with the cut-off wavenumber (Fig. 5(d)) presents the shape of a bell in the opt region where the residual increases. A so-called “optimum” cut-off wavenumber kc may be chosen as the one for which the slope becomes higher than a certain threshold. In the following, this threshold is chosen as 5 percent of the maximum slope (see Fig. 5(d)).
3.1.5. Results on noisy data with regularization The denoising procedure previously described is now applied on the same noisy displacement field as in Fig. 3 (SNR ¼40 dB), and in the case of higher noise levels (SNR¼30 and 20 dB). For the three levels of noise considered (Fig. 6), a very good identification of Young0 s modulus is achieved, except at low frequencies (Fig. 6(a)). As a good identification was achieved in the case of non-noisy data (see Fig. 3), the low-frequency fail is attributed to the filtering operation. Indeed it has been verified that the identification of Young0 s modulus becomes accurate when the size of the spatial response of the filter, i.e., twice the cut-off wavelength λc ¼ 2π=kc , is smaller than the size of the observation area. The denoising procedure is also
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2461
efficient for the identification of loss factor (Fig. 6(b)). However, the relative error on loss factor highly depends on the level of noise, whereas the relative error on Young0 s modulus is almost the same for the three noise levels. This is due to the fact that the loss factor is very small compared to unity. Therefore, the imaginary part of the flexural stiffness 3 D ¼ Eð1 þ jηÞh ð12ð1 ν2 ÞÞ 1 has a smaller contribution in the equation of motion than the real part. Fig. 7 shows the optimum cut-off wavenumber determined from the procedure explained previously, as well as the real pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi part of the natural flexural wavenumber kf ¼ 4 ðρh=DÞω2 . The evolution of optimum cut-off wavenumber with frequency follows the same trend as the flexural wavenumber, yet with a scale factor depending on the noise level. The noisier are the displacements, the closer is the optimum cut-off wavenumber to the flexural wavenumber. It should be noted that the same signal-to-noise ratio has been applied to all frequencies at which the displacement field was calculated. In an experimental context, the signal-to-noise ratio may vary with frequency and cannot be known precisely. Thus, the distance between the optimum cut-off wavenumber and the flexural wavenumber identified may serve as an indicator of the quality of the measurement at a given frequency. 3.2. Heterogeneous structures As the proposed identification technique is based on the local equation of motion, it may be applied to heterogeneous structures. The heterogeneity considered here corresponds to the presence of several subdomains in the structure, which have different material properties. However, it is assumed that the material is isotropic and homogeneous throughout each subdomain. As a consequence, the vibration within each subdomain is governed by an equation of the same form as Eq. (1). Such a structure may arise, for instance, from applying a local damping treatment to a plate (e.g. constrained-layer patch in viscoelastic material). This kind of heterogeneity has to be distinguished from a small-scale heterogeneity (e.g. honeycomb panel). In such a case, however, the method remains applicable as long as the behavior of the structure in every subdomain may be represented by equivalent parameters Deq and ðρhÞeq in an equation of motion like Eq. (1). In the following, the application of the method to a heterogeneous structure is illustrated using a simulated displacement field. 3.2.1. Direct problem using FEM We consider a plate with dimensions Lx ¼0.9 m, Ly ¼0.3 m and thickness h¼5 mm. In order to introduce a heterogeneity, the plate is partitioned into two areas, as shown in Fig. 8. The first area, covering most of the surface of the plate, has a relatively low damping (loss factor η ¼ 0:05). The second area, delimited by a circle with diameter d ¼0.14 m located in the middle of the plate, has a higher damping (η ¼ 0:3). Other material properties, which are common to both areas, are given in Table 2. With such a set of parameters, the real part of the flexural wavenumber kf, and therefore the flexural wavelength, do not differ significantly from one area to the other. The plate is excited by a harmonic point force at location (xF ¼ 0:07 m, yF ¼ 0:07 m). 250
opt
kc (SNR = 40 dB)
k (rad/m)
200
opt
kc (SNR = 30 dB)
150
kc (SNR = 20 dB)
100
Re[k f ]
opt
50 0
0
1000
2000
3000
4000
5000
f (Hz) opt
Fig. 7. Evolution of the optimum cut-off wavenumber kc
for several SNR, versus the frequency. – – – Real part of the natural flexural wavenumber kf.
low damping (η = 0.05)
high damping (η = 0.3)
y force Ød x
observation area
Fig. 8. Geometry of plate composed of two materials using the same characteristics except for the damping, which is used for numerical simulation.
2462
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
Table 2 Geometrical and material properties set for the direct problem. Lx ¼ 0.9 m Ly ¼ 0.3 m h ¼ 5 mm ρ ¼ 1500 kg m 3 ν ¼ 0:3 E ¼5 GPa
Length Width Thickness Density Poisson0 s ratio Young0 s modulus
regularize wi, j
solve
wi, j reg Ei, j
calculate derivatives 4 i, j
ηi, j
4 reg i, j
windowing
filtering
Fig. 9. Procedure to identify local material properties. The displacement field at f ¼1500 Hz with noise (SNR ¼40 dB) and a cut-off wavenumber kc ¼150 rad m 1 are used.
The resulting displacement field in the range ½10; 5000 Hz is obtained through a finite element simulation by using COMSOL. The mesh consists of 28,898 triangular elements. The size of each element if of the order of 5 mm. When considering the highest frequency, the mesh presents approximatively 12 elements per wavelength (λf 0:06 m). The numerical results are then processed with Matlab. The considered observation area is delimited by x A ½0:31; 0:89 m and y A ½0:01; 0:29 m. In order to apply the method, the displacement field has to be known over a uniform grid. An interpolation from nodal values is therefore required. This interpolation is carried out from the nodal transverse displacements into a 117 57 grid (Δx ¼ Δy ¼ 5 mm). It should be noted that this operation introduces a small error, which is yet considered negligible compared to the noise added afterwards. 3.2.2. Inverse problem on noisy data The purpose of the inverse problem is to recover the local material properties of the structure from the noisy displacement field (SNR ¼40 dB). The procedure, represented schematically in Fig. 9, is similar to the one described above. First, the spatial derivatives ∇4i;j are calculated from the discrete displacement field wi;j . Then, both terms are windowed and filtered to regularize the problem (the resulting fields are denoted by superscript “reg”). Finally, the quantity D=ρh in the equation of motion is obtained in a local manner, i.e., wreg D i;j ; (9) ¼ ω2 4 ρh i;j ∇i;j reg from which the local Young0 s modulus Ei;j and the local loss factor ηi;j may be then identified. Contrary to the case of homogeneous structures, where the local equation of motion is solved in a least-square sense, there is no residual associated to the solution of Eq. (9). Thus, the procedure previously described to automatically adjust the cut-off wavenumber is not applicable here. However, rather than characterizing stiffness and damping in a large frequency range, the proposed approach is intended to obtain a map of these properties at a single frequency. Therefore, the cut-off wavenumber is adjusted manually. Eq. (9) does not impose any condition on frequency to achieve a desired spatial accuracy when calculating a map of properties. However, it should be kept in mind that the regularization step involves the convolution of each field by the spatial response of the low-pass filter. This operation is likely to affect the sharpness of the map at the intersection between
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
|wi, j |
λc
ηi, j
0.3
y (m) λc
0 0.3
0.9
ηi, j
0.3
x (m)
0 0.3
x (m)
ηi, j
λc
0 0.3
0.9
0 0.3
0.9
0.3
y (m)
y (m)
0.3
x (m)
ηi, j
x (m)
λc
0.9
ηi, j
y (m)
0 0.3
λc
y (m)
0.3
y (m)
0.3
2463
x (m)
0
0.25
0 0.3
0.9
x (m)
0.9
0.5
Fig. 10. For a frequency f ¼1500 Hz with a simulated heterogeneous plate (SNR ¼40 dB), (a) the displacement field, (b)–(f) the evolutions of the identified local loss factor with the cut-off wavenumber, kc ¼ 50, 75, 100, 125 and 150 rad m 1 (the corresponding cut-off wavelength λc is indicated above each colormap). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
|wi, j |
λc
ηi, j
0.3
y (m) λc
0 0.3
0.9
ηi, j
0.3
x (m)
λc
0 0.3
x (m)
0.9
0 0.3
0.9
ηi, j
0.3
y (m)
y (m)
0.3
x (m)
ηi, j
x (m)
λc
0.9
ηi, j
y (m)
0 0.3
λc
y (m)
0.3
y (m)
0.3
0 0.3
0
x (m)
0.25
0.9
0 0.3
x (m)
0.9
0.5
Fig. 11. For a frequency f ¼5000 Hz with a simulated heterogeneous plate (SNR ¼40 dB), (a) the displacement field, (b)–(f) the evolutions of the identified local loss factor with the cut-off wavenumber, kc ¼ 50, 75, 100, 125 and 150 rad m 1 (the corresponding cut-off wavelength λc is indicated above each colormap). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
two regions with different properties. To illustrate the phenomenon, the inverse problem is applied for different values of cut-off wavenumber at two frequencies. The local loss factor is obtained from the noisy displacement field at the frequencies f ¼1500 Hz (Fig. 10) and f¼5000 Hz (Fig. 11) with five different levels of regularization, from kc ¼ 50 to 150 rad m 1 (the flexural wavenumbers at the considered frequencies are respectively kf ¼58 rad m 1 for f¼1500 Hz and kf ¼107 rad m 1 for f¼5000 Hz). Although the region with higher damping may be clearly localized on each figure, the value of cut-off wavenumber strongly influences the quality of the spatial reconstruction: at kf ¼50 rad m 1 (Figs. 10(b) and 11(b)), the contours of the inclusion are poorly reconstructed. At kf ¼75–100 rad m 1 (Fig. 10(c) and (d)) or kf ¼75–150 rad m 1 (Fig. 11(c) to (f)), the map reproduces the shape of the inclusion in a better way and the value of loss factor in each region may be inferred from the colorscale with satisfying accuracy. For f ¼1500 Hz and above kf ¼125 rad m 1 (Fig. 10(e) and (f)), the quantification of loss factor from the map is less accurate because of the noise, which is not removed efficiently.
2464
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
measurement area plate shaker impedance head
w
y
x F
scanning laser vibrometer
Fig. 12. Measurement set-up.
Two kinds of spurious patterns can also be distinguished on such maps. The first one, which does not depend on kc, consists in spurious lines in the map at the location of nodal lines in the displacement field (see e.g. Fig. 10(c) and (d) and Fig. 11(c)–(f)). This is explained by the fact that nodal lines of the displacement field are also nodal lines of the bilaplacian field (both fields are equal up to the constant multiplier D=ρhω2 , even once windowed and filtered). Thus, the estimation of ðD=ρhÞi;j from Eq. (9) consists in the ratio between two very small values. As these values are systematically below or comparable to the noise remaining after regularization, their ratio differs strongly from the actual value of ðD=ρhÞ at the considered points. The second kind of spurious pattern consists in “grain” in the map, at a size comparable to the cut-off wavelength λf. This is caused by the presence of spurious components at small wavelengths in the bilaplacian field, which are not contained in the displacement field and not removed by the filter (this phenomenon is well observable in Fig. 4(b)). 4. Experiments In the previous section, the identification technique has been validated on data from simulations. The purpose of the present section is to apply the procedure on real measurement data. The plate to test is suspended from a frame to approximate free boundary conditions (Fig. 12). A Polytec Scanning Vibrometer PSV 300 is used to measure the displacement field. The excitation is provided by a Brüel & Kjær 4810 shaker, supplied in power by a B&K 2718 amplifier. Although no knowledge of excitation level is required by the method, a B&K 8001 impedance head is used to provide both phase reference and input force measurement, in conjunction with a B&K NEXUS conditioning amplifier. 4.1. Homogeneous PMMA plate A first experiment was carried out on a polymethyl methacrylate (PMMA) plate. The plate dimensions are 0.5 m 0.4 m and thickness h¼3 mm. The material density is ρ¼1190 kg/m3. Because the material is transparent, one face of the plate was covered with silver paint to allow reflection of the laser beam. The displacements were measured over a meshgrid consisting of 71 71 ¼ 5041 points. The dimensions of the measurement area were 0.36 m 0.36 m, corresponding to a spatial step Δx ¼ Δy 5 mm. The plate was excited by the shaker at its bottom edge, outside the measurement area. The excitation signal was a periodic chirp in the frequency range [16, 6400] Hz. The displacement at each point was acquired using 10 averages to improve the signal-to-noise ratio. Young0 s modulus and loss factor identified by the method are shown in Fig. 13. For both properties, the dispersion around the mean curve is low. Regardless of small local variations inherent to the method, the evolution of these properties with frequency follows the same trend as the one observed by other authors for PMMA [16]. opt The optimum cut-off wavenumber kc does not appear to be proportional to the identified flexural wavenumber kf (Fig. 14), contrary to the results get by the simulations (see Fig. 7). This is due to the fact that the signal-to-noise ratio varies with frequency in the experiment. However, the optimum cut-off wavenumber is relatively high compared to the optimum opt wavenumber (kc =kN 1:7 in the range [2000, 6400] Hz), which confirms the good quality of the measurement. 4.2. Heterogeneous composite plates To assess the possibility to obtain a map of local properties from experimental data, two plates made by CETIM were used. One plate was made of PP thermoplastic matrix reinforced by glass fibers. It contains 10 patches of Santoprene™ elastomer. The other was made of a thermosetting epoxy matrix reinforced by glass fibers. It contained eight patches of an
10
0.1
8
0.08
6
0.06
2465
η
E (GPa)
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
4
0.04
2
0.02
0
0
1000 2000 3000 4000 5000 6000 f (Hz)
0
0
1000 2000 3000 4000 5000 6000 f (Hz)
Fig. 13. (a) Young0 s modulus E and (b) loss factor η identified on the PMMA plate by the method.
250
kf opt kc
k (rad/m)
200 150 100 50 0
0
1000
opt
Fig. 14. Optimum cut-off wavenumber kc
2000
3000 4000 f (Hz)
5000
6000
automatically found by the procedure, compared with the identified flexural wavenumber kf.
|wi, j |
y (m)
0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
(D/ ρh)i, j (m4 .s−2 )
ηi, j 1
0.3
40
0.3
0.2
30
0.2
20 0.1 0
10 0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
0
y (m)
y (m)
50
0.5
0.1 0
0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
0
Fig. 15. (a) Displacement field, (b) identified specific flexural rigidity, and (c) identified local loss factor. Plate made of glass fiber and PP matrix with 10 patches of Santoprene™ elastomer, f ¼6700 Hz, kc ¼150 rad m 1.
2466
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
Fig. 16. (a) Displacement field, (b) identified specific flexural rigidity, and (c) identified local loss factor. Plate made of glass fiber and epoxy matrix with eight patches of damping elastomer, f¼ 2000 Hz, kc ¼ 100 rad m 1.
elastomer specifically formulated so as to provide high damping. Both plates were made by co-curing the composite and elastomer. Contrary to classical damping treatments (e.g. constrained viscoelastic layer) which may be applied to a structure, the damping function is here directly embedded in the material. The plates were driven by a periodic chirp in the frequency range [500, 10,000] Hz. The excitation point was at the bottom edge, outside the scanning area. The displacements were measured on a fine mesh (around 20,000 points) with five averages. The measurement duration was in the order of several hours. On the first plate, the inverse problem was carried out with the displacement field obtained at frequency f¼6700 Hz on a mesh containing 211 101 ¼ 21; 311 points (Fig. 15(a)). The problem was regularized with cut-off wavenumber kc ¼150 rad m 1. The identified local specific flexural rigidity (Fig. 15(b)) clearly varies from the regions containing only composite material to those containing an elastomer core. Yet, no significant rise in the local loss factor is observed in the regions containing the elastomer (Fig. 15(c)). For the second plate, the displacements obtained at frequency f¼2000 Hz on a mesh containing 239 81 ¼ 19; 359 points were used in the inverse problem (Fig. 16(a)). A cut-off wavenumber kc ¼100 rad m 1 was chosen. Whereas the identified local specific flexural rigidity (Fig. 16(b)) does not vary significantly in the area considered, the local loss factor is much higher at the locations containing damping elastomer (Fig. 16(c)). It should be kept in mind that the test plates were made of composite materials. The inverse method should then be used with an a priori equation of motion of an orthotropic plate and the test should be done with several sources in order to generate waves in all directions. It is then particularly interesting to see that the fact to consider an equivalent isotropic plate here is sufficient for the understanding of the behavior of the inserted patches. The two different behaviors, modification of the local flexural stiffness (Fig. 15(b)) or the local damping (Fig. 16(c)), can be attributed to the quality of the adhesion between the patches and the composite layers, and to the material properties of each elastomer.
5. Conclusion In this paper, an adaptation of the Force Analysis Technique has been developed to identify stiffness and damping properties of plates. If the general form of the equation of motion is known, it is possible to identify structural parameters by measuring the vibratory displacement field over an area in which no external force is applied. While keeping the overall regularization procedure as in the Force Analysis Technique, i.e., windowing and filtering, the method has been improved by a novel automatic determination of the optimum cut-off wavenumber. Contrary to modal methods, the proposed approach is independent of boundary conditions and may be applied at any frequency, not necessarily a resonance. It is also valid in the mid-frequency domain, where the modal overlap is high. Moreover, it can be applied in the case of heterogeneous structures to obtain a map of local properties. The applicability of the method has been verified on simulated plates and experimental analysis on real plates. In this paper, the equation of motion of an isotropic plate has been considered. Following the same principle, the method can be applied to beams. A natural continuation of this work is to consider orthotropic plates, which are commonly encountered with composite materials.
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2467
Acknowledgments This study was funded by Région Pays de la Loire within the framework of the AMORTI project. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the manuscript, and Alain Lemascon for his careful reading of the manuscript. Appendix A. Finite differences scheme The fourth-order partial derivatives of the displacement field involved in the local equation of motion are calculated as follows: ∂4 w 1 δ4x i;j ¼ 4 wi þ 2;j 4wi þ 1;j þ 6wi;j 4wi 1;j þwi 2;j ∂x4 Δx ∂4 w 1 δ4y i;j ¼ 4 wi;j þ 2 4wi;j þ 1 þ 6wi;j 4wi;j 1 þwi;j 2 ∂y4 Δy ∂4 w 1 δ2x2y ¼ 2 2 wi þ 1;j þ 1 2wi þ 1;j þ wi þ 1;j 1 2wi;j þ 1 þ 4wi;j 2wi;j 1 þwi 1;j þ 1 2wi 1;j þwi 1;j 1 ; i;j ∂x2 ∂y2 Δx Δy
(A.1)
where Δx and Δy are the spatial steps in x and y directions. It clearly appears from these expressions that the spatial derivatives cannot be obtained at points close to the edges of the observation area. Thus the effective area in which the inverse problem can be applied only contains ðM x 4Þ ðM y 4Þ points. Appendix B. Expression of the spatial window and filter response In the regularization step of the inverse method, each field wi;j and ∇4i;j is first multiplied by a spatial Tukey window ψ 2D i;j . As the exploitable area for the inverse problem is delimited by i A ½3; M x 2 and j A ½3; M y 2, the window is implemented as follows: 1D ψ 2D ðxi;j x1;1 2Δx Þ ψ 1D ðyi;j y1;1 2Δy Þ; i;j ¼ ψ
where
8 πx > > 0:5 1 cos > > α > > > <1 1D ψ ðxÞ ¼ πðx L þ 2αÞ > > 0:5 1 cos > > α > > > : 0
(B.1)
if 0 r x o α; if α rx rL α; if L α ox r L;
(B.2)
else
is the expression of the unidimensional Tukey window, L ¼ Lx 4Δx and 8 < λc if L Z 2λc ; α¼ L else: : 2
(B.3)
Eq. (B.3) expresses the fact that the Tukey window is replaced by a Hanning window when the cut-off wavelength becomes too large compared to the size of the observation area. Once windowed, each field is convolved with the discrete spatial response hi;j ¼ hðxi;j ; yi;j Þ of a low-pass filter calculated using the following expression: 8
2 > kc kc x kc y 2π 2π < 1 þ cos sin k 1 þ cos ð x Þ sin ð k y Þ if x and y A ; ; c c 2 2 kc kc (B.4) hðx; yÞ ¼ 4π 2 xyN f > : 0 else; where Nf is a normalization parameter, chosen such that ∬ hðx; yÞ dx dy ¼ 1: References [1] D. Ewins, Modal Testing: Theory and Practice, Research Studies Press, Baldock, United Kingdom, 1984. [2] M. Matter, T. Gmuer, J. Cugnoni, A. Schorderet, Numerical-experimental identification of the elastic and damping properties in composite plates, Composite Structures 90 (2009) 180–187. [3] J. McDaniel, W. Shepard, Estimation of structural wave numbers from spatially sparse response measurements, Journal of the Acoustical Society of America 108 (2000) 1674–1682. [4] J. Berthaut, M. Ichchou, L. Jezequel, K-space identification of apparent structural behaviour, Journal of Sound and Vibration 280 (2005) 1125–1131. [5] M. Rak, M. Ichchou, J. Holnicki-Szulc, Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer, Journal of Sound and Vibration 310 (2008) 801–811.
2468
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
[6] M.N. Ichchou, O. Bareille, J. Berthaut, Identification of effective sandwich structural properties via an inverse wave approach, Engineering Structures 30 (2008) 2591–2604. [7] C. Pezerat, J.L. Guyader, 2 Inverse methods for localization of external sources exciting a beam, Acta Acustica 3 (1995) 1–10. [8] C. Pezerat, J.L. Guyader, Force analysis technique: reconstruction of force distribution on plates, Acustica 86 (2000) 322–332. [9] C. Pezerat, J.L. Guyader, Identification of vibration sources, Applied Acoustics 61 (2000) 309–324. [10] C. Pezerat, Q. Leclere, N. Totaro, M. Pachebat, Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique, Journal of Sound and Vibration 326 (2009) 540–556. [11] Q. Leclere, C. Pezerat, Vibration source identification using corrected finite difference schemes, Journal of Sound and Vibration 331 (2012) 1366–1377. [12] M. Djamaa, N. Ouelaa, C. Pezerat, J.L. Guyader, Mechanical radial force identification of a finite cylindrical shell by an inverse method, Acta Acustica united with Acustica 92 (2006) 398–405. [13] M.C. Djamaa, N. Ouelaa, C. Pezerat, J.L. Guyader, Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering, Journal of Sound and Vibration 301 (2007) 560–575. [14] C. Renzi, C. Pezerat, J.L. Guyader, Vibratory source identification by using the Finite Element Model of a subdomain of a flexural beam, Journal of Sound and Vibration 332 (2013) 545–562. [15] J.L. Guyader, Vibration in Continuous Media, ISTE Ltd, London, United Kingdom, 2006. [16] Y. Liao, V. Wells, Estimation of complex modulus using wave coefficients, Journal of Sound and Vibration 295 (2006) 165–193.
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
Identification of stiffness and damping properties of plates by using the local equation of motion Frédéric Ablitzer a,n, Charles Pézerat a, Jean-Michel Génevaux a, Jérôme Bégué b a b
LUNAM Université, Université du Maine, CNRS UMR 6613, LAUM, Avenue Olivier Messiaen, 72085 Le Mans Cedex 9, France CETIM (Centre Technique des Industries Mécaniques), Technocampus EMC2, Chemin du Chaffault, 44340 Bouguenais, France
a r t i c l e in f o
abstract
Article history: Received 15 January 2013 Received in revised form 12 December 2013 Accepted 13 December 2013 Handling Editor: I. Trendafilova Available online 23 January 2014
This paper deals with the identification of stiffness and damping properties of vibrating structures by an inverse method inspired from the Force Analysis Technique (FAT). The proposed approach uses a local equation of motion assumed a priori, which provides a relative straightforward relationship between the displacement field and material properties. The spatial derivatives of the displacement in the equation are calculated using finite differences. As this operation amplifies measurement noise, a regularization step is applied before solving the inverse problem. A procedure is proposed to automatically adjust the level of regularization. The method also allows one to identify local stiffness and damping on a heterogeneous structure. Illustrations for both homogeneous and heterogeneous cases are shown using simulated and measured displacement fields. & 2013 Elsevier Ltd. All rights reserved.
1. Introduction Composite materials, which provide high stiffness and light weight, are increasingly used in the aerospace and automotive industry. During the conception of structures including such materials, an extensive use of finite element models is made to predict the dynamic and vibroacoustic behavior. To achieve predictions that are as accurate as possible, one of the key challenges is to enter the right material characteristics into the model. As a wide variety of composite materials exist, their elastic and damping characteristics are rarely available in tables and must be obtained experimentally. Moreover, if the finite element method uses plate elements, the determination of the characteristics of the equivalent plate (stiffness and damping for all types of loads) needs the choice of kinematic assumptions. This choice influences the results of the model and has to be confirmed by measurements. Measurement techniques based on modal analysis theory, which is well established [1], are widely used. The measurement of modal parameters (natural frequencies, mode shapes and damping ratio) on a test structure, generally combined to a numerical model, allows one to identify elastic and damping properties of the constitutive material [2]. However, there are several limitations to such methods. First, as it is a global approach, it is not possible to characterize the material directly on the target structure itself: experiments are generally carried out on specific test specimens, a precise knowledge of boundary conditions being crucial. Secondly, modal parameters are difficult to measure in the mid-frequency domain, characterized by a strong modal overlap. Because the stiffness and damping of composites may strongly vary with frequency, it is necessary to obtain material properties in a wide frequency range instead of extrapolating results obtained at low frequencies. In recent years, works have been done on the estimation of wavenumbers in a structure, from which material properties such as loss factor or dynamic stiffness may be deduced [3–6]. Such wave-based approaches overcome most limitations
n
Corresponding author. E-mail address: [email protected] (F. Ablitzer).
0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.12.013
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2455
of the modal approach. In particular, they are applicable at any frequency and do not require any knowledge about boundary conditions. Their principle is to find the wavenumber giving the best fit between the displacement field measured at a given driving frequency and a theoretical field consisting of a sum of waves. This theoretical field can be viewed as the solution of an equation of motion involving material properties, which may be related to the identified wavenumber. This paper presents a new method to identify material properties from the displacement field measured locally on a structure. The underlying idea is to use the local equation of motion itself, which provides a relative straightforward relationship between the material properties and the displacement field. The method is inspired from the Force Analysis Technique (FAT, also known by the french acronym RIFF), initially aimed at identifying dynamic loads acting on structures [7–11]. In the Force Analysis Technique, the measured displacement field is injected into a discrete form of the local equation of motion, in which the spatial derivatives are approximated by a finite-difference scheme, leading to a straightforward calculation of the force acting in the area considered. This simple principle, combined with an appropriate regularization procedure (windowing and filtering), has been successfully applied to beams [7], plates [8,9], shells [12,13], and thereafter extended to more complex structures by using a finite element operator instead of an analytical equation of motion [14]. To localize and quantify external forces, the Force Analysis Technique obviously requires the a priori knowledge of material properties involved in the equation of motion. The originality of the present method is to identify unknown material properties from the equation of motion by considering an area of the structure which is free of any external force. The key ingredients for the identification of the material properties (approximation of partial derivatives, regularization) are the same as those used in FAT, complemented with specific developments to compensate the lack of a priori knowledge about the structure. In particular, a novel automatic regularization technique is proposed. The paper focuses on the application to plates. However, the method may be extended to other types of structures for which a differential equation exists. The continuation of this paper is structured as follows. First, the general principle of the method is exposed more in detail. Then, the feasibility of the proposed approach is assessed using simulated displacement fields. The effect of noise is illustrated and a procedure to automatically adjust the regularization is described. The applicability of the method is assessed for both homogeneous and heterogeneous structures. The last part presents experimental results. 2. Principle of the identification technique The equation of motion of an isotropic thin plate in the harmonic regime is [15] 4 ∂ w ∂4 w ∂4 w D ρhω2 wðx; yÞ ¼ f ðx; yÞ; þ 2 þ ∂x4 ∂x2 ∂y2 ∂y4
(1)
where D is the flexural stiffness, ρ the density, h the thickness, ω the angular frequency, wðx; yÞ the transverse displacement and f ðx; yÞ the force per unit area. The material is assumed to be linear viscoelastic. Therefore, the flexural stiffness 3
D¼
Eð1 þ jηÞh 12ð1 ν2 Þ
(2)
involves a complex Young0 s modulus Eð1 þ jηÞ, where j is the unit imaginary number and η denotes the loss factor, which characterizes material damping. The equation of motion (1) is local and therefore valid everywhere in the structure. More particularly, considering an area on which no external force is applied, i.e., f ðx; yÞ ¼ 0, Eq. (1) becomes 4 D ∂ w ∂4 w ∂4 w þ2 þ ¼ wðx; yÞ: (3) ρhω2 ∂x4 ∂x2 ∂y2 ∂y4 Eq. (3) simply expresses the fact that the bilaplacian ∇4 wðx; yÞ equals the displacement wðx; yÞ up to a constant multiplier. The constant multiplier D=ρhω2 depends on the characteristics of the plate, and on the angular frequency ω. It should be noted that the quantity D=ρh is complex and may vary with frequency. It contains information about the stiffness and damping of the material. Thus, assuming that the displacement wðx; yÞ is known, it seems possible to determine the quantity D=ρh from Eq. (3), at any frequency, without any knowledge about the exact location and amplitude of the external force, or about the boundary conditions. The only requirement is that no external force is applied at the considered location (x,y). In practice, however, there are several obstacles to this simple idea. The first one is due to the presence of fourth-order partial derivatives in Eq. (3). Whereas the transverse displacement at a given location may be easily obtained by using an accelerometer or a laser vibrometer, its partial derivatives are less straightforward to obtain. To overcome this difficulty, one may measure the displacement at discrete abscissas ðxi ; yi Þ over a regular meshgrid, so that the partial derivatives may be approximated by a finite difference scheme. A discretized form of Eq. (3) is D 4x ¼ wi;j : δi;j þ 2δ2x2y þ δ4y (4) i;j i;j 2 ρhω 4y 2x2y The expressions of δ4x may be found in Appendix A. The approximation of the partial derivatives introduces i;j , δi;j and 2δi;j a bias which is likely to affect the value of D=ρh calculated from Eq. (4). Note that the corrected finite differences scheme proposed in [11] can also be applied.
2456
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
A second difficulty is due to the noise which corrupts the measured displacements. Therefore, the quantity D=ρh identified by the method is likely to vary depending on the point considered. To overcome this difficulty, one may use a set of N data points and find a unique value of D=ρh that satisfy Eq. (4) in a least-square sense, i.e., 2 3 2 3 w1 δ4x þ 2δ2x2y þ δ4y 1 1 6 1 7 D 6 7 6 7 ⋮ (5) 4 5 ρhω2 ¼ 4 ⋮ 5: 2x2y 4y 4x w δN þ 2δN þ δN N However, identifying material properties through the least squares solution of Eq. (5) is not sufficient to reduce the effect of measurement noise. As will be shown in the next section, the effect of noise which is present in the measured displacements wi;j is considerably amplified by the fourth-order partial derivatives, making this simple identification technique not applicable. With an appropriate denoising treatment, it is yet possible to overcome this obstacle, as shown in the following section. 3. Proof-of-concept using simulations 3.1. Homogeneous structures 3.1.1. Direct problem We consider a simply supported isotropic thin plate (Fig. 1), whose characteristics are given in Table 1. The plate is driven at angular frequency ω by a harmonic point force at location (x0 ¼ 0:05 m, y0 ¼ 0:05 m) with amplitude F 0 ¼ 1 N. The displacement field is calculated by modal expansion, i.e., M
wexact ðx; yÞ ¼ ∑
N
∑
4F 0
m ¼ 1 n ¼ 1 Lx Ly ρh
sin ðmπx=Lx Þ sin ðnπy=Ly Þ sin ðmπx0 =Lx Þ sin ðnπy0 =Ly Þ ; ω2mn ω2
(6)
pffiffiffiffiffiffiffiffiffiffiffi where ωmn ¼ D=ρhððmπ=Lx Þ2 þ ðnπ=Ly Þ2 Þ. The number of modes M and N is chosen so as to satisfy a spatial Shannon0 s criterion. The displacements are calculated at discrete abscissas over a Mx M y ¼ 51 51 meshgrid, delimited by x A ½0:2 m; 0:45 m and y A ½0:1 m; 0:35 m. This area does not contain the excitation point, which is a requirement to apply the proposed identification technique. The corresponding spatial step is Δx ¼ Δy ¼ 5 mm. In order to simulate measurement uncertainties, noise can be added to the displacements. The noisy displacements at every frequency are obtained as follows: wnoisy
i;j
¼ wexact
jΔφ i;j þ ΔA i;j e
i;j
;
(7)
where ΔA i;j is a Gaussian random value with mean zero and standard deviation schosen in order to reach a desired signalto-noise ratio (SNR), and Δφ i;j is a uniform random value in the range ½0; 2π. The signal-to-noise ratio is calculated as follows: SNR ¼ 20 log10
∑i;j jwexact i;j j2 ∑i;j jwnoisy
2 i;j wexact i;j j
:
(8)
wi;j without subscript refers to noisy displacements, whereas subscript “exact” is kept when referring to exact displacements. 3.1.2. Inverse problem The aim of the inverse problem is to recover some of the plate properties from the available displacement field. First, the feasibility of the proposed approach in the case of non-noisy data is shown. Then, we illustrate the problem arising from noisy measurements and propose a denoising treatment allowing the identification of properties. The moduli of the discrete field wexact i;j calculated by solving the direct problem, as well as the discrete approximation of 2x2y the bilaplacian operator ∇4i;j ¼ ðδ4x þ δ4y i;j þ 2δi;j i;j Þ, are represented in Fig. 2(a) and (b). A very good fit between both fields can be visually observed, yet with different colorscales, as expected from Eq. (4) which implies that jD=ρhω2 jjδ4x i;j þ 2δ2x2y þδ4y i;j i;j j ¼ jwi;j j. This linear relationship between the moduli may also be observed in Fig. 2(c), in which the exact displacements and the bilaplacian are plotted against each other for all points in the considered domain. At this particular
Table 1 Geometrical and material properties set for the direct problem. Length Width Thickness Density Poisson0 s ratio Young0 s modulus Loss factor
Lx ¼ 0.5 m Ly ¼ 0.4 m h ¼ 5 mm ρ ¼1500 kg m 3 ν ¼ 0:3 E ¼5 GPa η ¼ 0:05
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2457
y force
x
observation area
Fig. 1. Geometry of plate used for numerical simulation.
2x2y
×10−6
4y
0.35
0.35
6
0.1 0.2
0.45
4
yi (m)
yi (m)
1
0.5
2 0.1 0.2
0
xi (m)
wi, j (m)
×10−6 1.5
−3 δ4x i, j + 2δi, j + δi, j (m )
wexact i, j (m)
wexact i, j (m)
1 0.5 0
0.45 xi (m) 2x2y
×10−6
4y
×10−6 1.5
−3 δ4x i, j + 2δi, j + δi, j (m )
0.35
0.35
yi (m)
yi (m)
600 400
wi, j (m)
800
1
0.5
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
1 0.5
200 0.1 0.2
0.45 xi (m)
0
0.1 0.2
0.45 xi (m)
0
0
0 500 1000 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
Fig. 2. Colormaps of the moduli of the two terms in Eq. (5), which are also plotted against each other, at f¼ 1000 Hz. (a)–(c) without noise; (d)–(f) with noise on wi;j (SNR ¼40 dB). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
frequency, knowing ρ, h and ν and solving Eq. (5) we can calculate Young0 s modulus E ¼5.03 GPa and the loss factor η ¼ 4:99 10 2 . The identified values are very close to the actual ones (Table 1). The moduli of the same discrete fields are then plotted in the case of noisy data, with an SNR of 40 dB. Whereas the noisy displacement field wi;j (Fig. 2(d)) is not distinguishable from the exact field wexact i;j (Fig. 2(a)), the discrete bilaplacian ∇4i;j calculated from noisy displacements exhibits a high randomness. Not only the amplitude is considerably amplified compared to the non-noisy case (the range in values is scaled by about 130), but the distribution of energy in the wavenumber domain is also affected. The bilaplacian obtained from non-noisy displacements (Fig. 2(b)) exhibits a relative high spatial wavelength compared to the spatial step. This wavelength corresponds to the natural flexural wavelength of the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi structure λf ¼ 2π=kf at the considered frequency (kf ¼ 4 ðρh=DÞω2 is the natural wavenumber). The bilaplacian calculated now from noisy displacements (Fig. 2(e)) contains wavelengths in the order of the spatial step, as inferred from the huge variations in amplitude among neighboring points. As observed in Fig. 2(f), the linear relationship between the moduli of the displacements and the bilaplacian is not verified. The identification of Young0 s modulus E and the loss factor η in a large frequency range, f A ½12:5; 5000 Hz, is then performed. At each frequency, the displacement field wi;j is obtained by solving the direct problem. Then, the discrete bilaplacian ∇4i;j is calculated and Eq. (5) is solved in a least squares sense. Results obtained from the displacements without
2458
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
10
0.1 from exact displacements from noisy displacements
0.08
6
0.06 η
E (GPa)
8
4
0.04
2
0.02
0
0
1000
2000 3000 f (Hz)
4000
5000
0
0
1000
2000 3000 f (Hz)
4000
5000
Fig. 3. (a) Young0 s modulus E and (b) loss factor η identified. The arrows on both sides indicate the actual values of E and η.
and with noise (SNR ¼40 dB) are shown in Fig. 3. When applying the method to non-noisy data, both elastic and damping properties can be identified with good accuracy. The small deviation from the actual values, observed especially on Young0 s modulus, is due to the bias error introduced by the approximation of partial derivatives in the equation of motion. It should be reminded that no approximation regarding the derivatives in the operator is used for the direct problem, which is solved by modal decomposition. When using noisy displacements, the identification of Young0 s modulus is no longer reliable. A very low value is found 2x2y over the considered frequency range (Fig. 3(a)). This is due to the fact that the term δ4x þ δ4y i;j þ 2δi;j i;j in Eq. (4) is very high, because of the amplification of noise, which has to be counterbalanced by a small value of stiffness D. The identification of the loss factor using noisy data seems to be compromised only at low frequencies (Fig. 3(b)). When frequency increases, the values identified are closer to the actual ones. This result is the consequence of the finite difference scheme, which acts as a low-pass filter in the wavenumber domain when the wavelength becomes small compared to the spatial step [11]. However, this effect is not enough to ensure the robustness of the method, as it does not allow a good identification of stiffness. Moreover, the bias error due to the finite difference scheme increases as the wavelength decreases compared to the spatial step. The denoising treatment proposed in the next section is intended to keep the accuracy allowed by a fine spatial discretization, without suffering the consequences of noise amplification. 3.1.3. Regularization: windowing and filtering As seen in the previous section, the noise in the displacement field wi;j is considerably amplified by the finite difference scheme. The calculated discrete bilaplacian ∇4i;j is therefore corrupted by noise at high wavenumbers, impairing the identification of the elastic and damping properties. In order to eliminate high wavenumbers components, a low-pass filter in the wavenumber domain is applied to the displacement field wi;j and the bilaplacian ∇4i;j . For this purpose, a similar procedure as the one applied in the Force Analysis Technique is followed. First, each field is windowed by a bidimensional Tukey window ψ 2D i;j . This preliminary step replaces the truncation of the field at the edges by a smooth variation from zero amplitude. Then, the windowed field is convolved by the finite spatial response hi;j of a low-pass filter with cut-off wavenumber kc. The expression of the window and the filter response are given in Appendix B. In the Force Analysis Technique, the cut-off wavenumber is generally chosen proportional to the natural wavenumber kf, i.e., kc ¼ akf . The coefficient a is usually chosen between 1 and 4, depending on the signal-to-noise ratio: the higher the SNR, the higher the cut-off wavenumber. This way of choosing the level of filtering is allowed by the fact that the properties of the plate and therefore its natural wavenumber are a priori known. We propose a method able to determine some of these properties which are unknown. Therefore, a new criterion to select the cut-off wavenumber has to be determined. The effect of cut-off wavenumber on the windowed filtered fields is illustrated in Fig. 4. The four values of the considered cut-off wavenumber kc correspond to a ¼4, 3, 2 and 1. As kc diminishes, each field contains less and less components with small wavelengths. Whereas the effect is hardly distinguishable on the displacement field, it appears clearly on the bilaplacian field. At kc ¼ 4kf (Fig. 4(a)), the filter keeps too much information at high wavenumbers to make any deterministic pattern observable in j∇4i;j j. At kc ¼ 3kf (Fig. 4(b)), the pattern looks like the one of jwi;j j but is still corrupted by small wavelengths. A good agreement between the displacement field and the bilaplacian is obtained at kc ¼ 2kf (Fig. 4(c)). The filter is effective enough. At kc ¼kf (Fig. 4(d)), the filter brings no more significant noise reduction, but continues to remove energy to both fields, as pointed out by the change in colorscale of wi;j . The effect of the filter is also well observable when the displacements and the bilaplacian are plotted against each other. As kc decreases, less and less values with high amplitude are kept in j∇4i;j j. Whereas the scale for jwi;j j is the same as in Fig. 2(f) (no filter), the scale for j∇4i;j j is diminished by a factor of 100. At the same time the dispersion of data points is reduced, especially for j∇4i;j j. At kc ¼ 2kf (Fig. 4(c)), the same linear relationship as in Fig. 2(c) appears, yet with a small dispersion and a reduction in amplitude along both axes due to windowing and filtering. Diminishing further kc removes energy to both fields (Fig. 4(d)). At this point, it appears that the cut-off wavenumber kc should be chosen low enough to diminish the noise level in ∇4i;j , but not too low in order to keep as much information as possible. A value satisfying this qualitative criterion may be simply
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
8
4 2
2
wi, j (m)
0.45
0.1 0.2
×10−7
wi, j (m)
xi (m)
wi, j (m)
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
0.1 0.2
×10−6
0.35
xi (m) wi, j (m)
0.45
0
×10−7
0.35
1
yi (m)
0.5
4
0.5
2 1
0.1 0.2
xi (m)
0.45
0
0.1 0.2
×10−6 1.5 wi, j (m)
0.5
0.45
0
3
0.45
1
xi (m)
1
×10−6 1.5 wi, j (m)
0.1 0.2
×10−6
yi (m)
yi (m)
10 8 6 4 2
0
0.45
1 0.5
0
0.35
0.35
0.1 0.2
xi (m)
2
2 1.5
yi (m)
xi (m)
0
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
4
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
xi (m)
0.45
0
0.1 0.2
×10−6 1.5 wi, j (m)
0.1 0.2
4y
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
xi (m)
0.45
0
×10−6 1.5 wi, j (m)
4
yi (m)
yi (m)
6
2x2y
−3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
6
yi (m)
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
yi (m)
2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m ) 0.35
2459
1 0.5 0
0 5 10 2x2y 4y −3 δ4x i, j + 2δi, j + δi, j (m )
Fig. 4. Effect of cut-off wavenumber kc on displacement wi;j and discrete bilaplacian ∇4i;j represented as colormaps and plotted against each other, in modulus. (a) kc ¼ 4kf , (b) kc ¼ 3kf , (c) kc ¼ 2kf , and (d) kc ¼kf. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
found by varying kc and observing plots like those represented in Fig. 4. For example, for the frequency and noise level considered in this figure, a value close to kc ¼ 2kf would be appropriate. The corresponding Young0 s modulus and loss factor identified are E¼5.03 GPa and η ¼ 5:06 10 2 , respectively. These values are very close to those found in the absence of noise. This shows the efficiency of regularization for this simulated noisy displacement field. 3.1.4. Automatic determination of optimum cut-off wavenumber The range of appropriate values of kc depends on the driving frequency and noise level, according to the Force Analysis Technique. The task of manually searching a value of kc for each frequency in a wideband measurement would be tedious. Therefore, it appears necessary to define a calculable indicator that expresses the qualitative criterion described above. To this aim, the solution of Eq. (5) is now focused on. Young0 s modulus (Fig. 5(a)) and the loss factor (Fig. 5(b)) are identified from the least squares solution of Eq. (5) for different values of cut-off wavenumber. The axis corresponding to kc is shown in a logarithmic scale, as the value for which the regularized solution converges to the non-regularized one (kc =kf ≳20) is very high compared to the range of interest. The identified Young0 s modulus is very small for high values of cut-off wavenumber (kc =kf ≳6). On the contrary, it is very close to the actual Young0 s modulus for low values of cut-off wavenumber (kc =kf ≲3), even when kc is significantly lower than the flexural wavenumber of the plate kf. This result may seem strange, as the energy of the displacement field in the wavenumber domain is essentially contained in the vicinity of kf. However, it can be explained by the fact that Eq. (3) remains valid when applying a linear filter to its terms, as the quantity D=ρhω2 does not depend on the spatial location. Thus, as long as the residual noise (i.e., with components lower that kc in the wavenumber domain) is not high enough to impair the identification, there is no reason that Young0 s modulus found should significantly differ from the actual one. The identified loss factor is also far from the actual value when the cut-off wavenumber is high. However, as the cut-off wavenumber decreases, the identification becomes satisfactory slightly earlier (kc =kf ≲5) than in the case of Young0 s modulus. This could be due to the fact that the imaginary part of the displacement field is much smaller than the real part, giving less weight to the amplification by the finite differences scheme. Thus, the filter needs not to be as restrictive as the one for Young0 s modulus identification. However, as the identification of both Young0 s modulus and structural loss factor is generally desired, the acceptable range in cut-off wavenumber has to be common to both quantities. 2 My 4 2 x The sum of squared residuals associated to the least squares solution, R ¼ ∑M i ¼ 0 ∑j ¼ 0 9∇i;j ðD=ρhω Þ wi;j 9 , is represented in Fig. 5(c). As expected, the residual is maximum when the cut-off wavenumber kc is very high, due to the fact that much noise is kept, especially in ∇4i;j . On the contrary, it is minimum at low values of kc, indicating a very good fit. It can be seen
2460
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
0.1
8
0.08
6
0.06 η
E (GPa)
10
4
0.04
2
0.02
0 10−1
100
101
0 10−1
102
100
kc / k f
102
101
102
kc / k f
8 ∂R/ ∂kc (normalized)
1
6 R (m2 )
101
4 2 0 10−1
100
101
0.5
threshold 0 10−1
102
100
kc / k f
kc / k f
Fig. 5. Effect of the level of regularization (cut-off wavenumber kc) on the results, (a) Young0 s modulus E, (b) loss factor η, (c) residual term R and (d) its derivative. The circles in (a) and (b) point out the values of E and η at the optimum cut-off wavenumber, for which the slope of residual reaches 5% of its maximum value.
10
0.1 SNR = 40 dB SNR = 30 dB SNR = 20 dB
0.08
6
0.06 η
E (GPa)
8
4
0.04
2
0.02
0
0 0
1000
2000 3000 f (Hz)
4000
5000
0
1000
2000 3000 f (Hz)
4000
5000
Fig. 6. (a) Young0 s modulus E and (b) loss factor η identified from noisy displacements, with a denoising treatment (windowing and filtering). Three levels of noise are considered. The arrows on both sides indicate the actual values of E and η.
that the value of cut-off wavenumber, at which the residual starts to rise significantly, kc =kf 3, corresponds to the value at which the identified Young0 s modulus starts to deviate from the actual Young0 s modulus (see Fig. 5(a)). The sudden rise in residual thus appears as a good detection indicator of the value of cut-off wavenumber at which the filter becomes efficient enough. This sudden rise may be easily found by considering the slope of the residual curve, ∂R=∂kc , which is obtained by discrete differentiation. The evolution of the slope with the cut-off wavenumber (Fig. 5(d)) presents the shape of a bell in the opt region where the residual increases. A so-called “optimum” cut-off wavenumber kc may be chosen as the one for which the slope becomes higher than a certain threshold. In the following, this threshold is chosen as 5 percent of the maximum slope (see Fig. 5(d)).
3.1.5. Results on noisy data with regularization The denoising procedure previously described is now applied on the same noisy displacement field as in Fig. 3 (SNR ¼40 dB), and in the case of higher noise levels (SNR¼30 and 20 dB). For the three levels of noise considered (Fig. 6), a very good identification of Young0 s modulus is achieved, except at low frequencies (Fig. 6(a)). As a good identification was achieved in the case of non-noisy data (see Fig. 3), the low-frequency fail is attributed to the filtering operation. Indeed it has been verified that the identification of Young0 s modulus becomes accurate when the size of the spatial response of the filter, i.e., twice the cut-off wavelength λc ¼ 2π=kc , is smaller than the size of the observation area. The denoising procedure is also
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2461
efficient for the identification of loss factor (Fig. 6(b)). However, the relative error on loss factor highly depends on the level of noise, whereas the relative error on Young0 s modulus is almost the same for the three noise levels. This is due to the fact that the loss factor is very small compared to unity. Therefore, the imaginary part of the flexural stiffness 3 D ¼ Eð1 þ jηÞh ð12ð1 ν2 ÞÞ 1 has a smaller contribution in the equation of motion than the real part. Fig. 7 shows the optimum cut-off wavenumber determined from the procedure explained previously, as well as the real pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi part of the natural flexural wavenumber kf ¼ 4 ðρh=DÞω2 . The evolution of optimum cut-off wavenumber with frequency follows the same trend as the flexural wavenumber, yet with a scale factor depending on the noise level. The noisier are the displacements, the closer is the optimum cut-off wavenumber to the flexural wavenumber. It should be noted that the same signal-to-noise ratio has been applied to all frequencies at which the displacement field was calculated. In an experimental context, the signal-to-noise ratio may vary with frequency and cannot be known precisely. Thus, the distance between the optimum cut-off wavenumber and the flexural wavenumber identified may serve as an indicator of the quality of the measurement at a given frequency. 3.2. Heterogeneous structures As the proposed identification technique is based on the local equation of motion, it may be applied to heterogeneous structures. The heterogeneity considered here corresponds to the presence of several subdomains in the structure, which have different material properties. However, it is assumed that the material is isotropic and homogeneous throughout each subdomain. As a consequence, the vibration within each subdomain is governed by an equation of the same form as Eq. (1). Such a structure may arise, for instance, from applying a local damping treatment to a plate (e.g. constrained-layer patch in viscoelastic material). This kind of heterogeneity has to be distinguished from a small-scale heterogeneity (e.g. honeycomb panel). In such a case, however, the method remains applicable as long as the behavior of the structure in every subdomain may be represented by equivalent parameters Deq and ðρhÞeq in an equation of motion like Eq. (1). In the following, the application of the method to a heterogeneous structure is illustrated using a simulated displacement field. 3.2.1. Direct problem using FEM We consider a plate with dimensions Lx ¼0.9 m, Ly ¼0.3 m and thickness h¼5 mm. In order to introduce a heterogeneity, the plate is partitioned into two areas, as shown in Fig. 8. The first area, covering most of the surface of the plate, has a relatively low damping (loss factor η ¼ 0:05). The second area, delimited by a circle with diameter d ¼0.14 m located in the middle of the plate, has a higher damping (η ¼ 0:3). Other material properties, which are common to both areas, are given in Table 2. With such a set of parameters, the real part of the flexural wavenumber kf, and therefore the flexural wavelength, do not differ significantly from one area to the other. The plate is excited by a harmonic point force at location (xF ¼ 0:07 m, yF ¼ 0:07 m). 250
opt
kc (SNR = 40 dB)
k (rad/m)
200
opt
kc (SNR = 30 dB)
150
kc (SNR = 20 dB)
100
Re[k f ]
opt
50 0
0
1000
2000
3000
4000
5000
f (Hz) opt
Fig. 7. Evolution of the optimum cut-off wavenumber kc
for several SNR, versus the frequency. – – – Real part of the natural flexural wavenumber kf.
low damping (η = 0.05)
high damping (η = 0.3)
y force Ød x
observation area
Fig. 8. Geometry of plate composed of two materials using the same characteristics except for the damping, which is used for numerical simulation.
2462
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
Table 2 Geometrical and material properties set for the direct problem. Lx ¼ 0.9 m Ly ¼ 0.3 m h ¼ 5 mm ρ ¼ 1500 kg m 3 ν ¼ 0:3 E ¼5 GPa
Length Width Thickness Density Poisson0 s ratio Young0 s modulus
regularize wi, j
solve
wi, j reg Ei, j
calculate derivatives 4 i, j
ηi, j
4 reg i, j
windowing
filtering
Fig. 9. Procedure to identify local material properties. The displacement field at f ¼1500 Hz with noise (SNR ¼40 dB) and a cut-off wavenumber kc ¼150 rad m 1 are used.
The resulting displacement field in the range ½10; 5000 Hz is obtained through a finite element simulation by using COMSOL. The mesh consists of 28,898 triangular elements. The size of each element if of the order of 5 mm. When considering the highest frequency, the mesh presents approximatively 12 elements per wavelength (λf 0:06 m). The numerical results are then processed with Matlab. The considered observation area is delimited by x A ½0:31; 0:89 m and y A ½0:01; 0:29 m. In order to apply the method, the displacement field has to be known over a uniform grid. An interpolation from nodal values is therefore required. This interpolation is carried out from the nodal transverse displacements into a 117 57 grid (Δx ¼ Δy ¼ 5 mm). It should be noted that this operation introduces a small error, which is yet considered negligible compared to the noise added afterwards. 3.2.2. Inverse problem on noisy data The purpose of the inverse problem is to recover the local material properties of the structure from the noisy displacement field (SNR ¼40 dB). The procedure, represented schematically in Fig. 9, is similar to the one described above. First, the spatial derivatives ∇4i;j are calculated from the discrete displacement field wi;j . Then, both terms are windowed and filtered to regularize the problem (the resulting fields are denoted by superscript “reg”). Finally, the quantity D=ρh in the equation of motion is obtained in a local manner, i.e., wreg D i;j ; (9) ¼ ω2 4 ρh i;j ∇i;j reg from which the local Young0 s modulus Ei;j and the local loss factor ηi;j may be then identified. Contrary to the case of homogeneous structures, where the local equation of motion is solved in a least-square sense, there is no residual associated to the solution of Eq. (9). Thus, the procedure previously described to automatically adjust the cut-off wavenumber is not applicable here. However, rather than characterizing stiffness and damping in a large frequency range, the proposed approach is intended to obtain a map of these properties at a single frequency. Therefore, the cut-off wavenumber is adjusted manually. Eq. (9) does not impose any condition on frequency to achieve a desired spatial accuracy when calculating a map of properties. However, it should be kept in mind that the regularization step involves the convolution of each field by the spatial response of the low-pass filter. This operation is likely to affect the sharpness of the map at the intersection between
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
|wi, j |
λc
ηi, j
0.3
y (m) λc
0 0.3
0.9
ηi, j
0.3
x (m)
0 0.3
x (m)
ηi, j
λc
0 0.3
0.9
0 0.3
0.9
0.3
y (m)
y (m)
0.3
x (m)
ηi, j
x (m)
λc
0.9
ηi, j
y (m)
0 0.3
λc
y (m)
0.3
y (m)
0.3
2463
x (m)
0
0.25
0 0.3
0.9
x (m)
0.9
0.5
Fig. 10. For a frequency f ¼1500 Hz with a simulated heterogeneous plate (SNR ¼40 dB), (a) the displacement field, (b)–(f) the evolutions of the identified local loss factor with the cut-off wavenumber, kc ¼ 50, 75, 100, 125 and 150 rad m 1 (the corresponding cut-off wavelength λc is indicated above each colormap). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
|wi, j |
λc
ηi, j
0.3
y (m) λc
0 0.3
0.9
ηi, j
0.3
x (m)
λc
0 0.3
x (m)
0.9
0 0.3
0.9
ηi, j
0.3
y (m)
y (m)
0.3
x (m)
ηi, j
x (m)
λc
0.9
ηi, j
y (m)
0 0.3
λc
y (m)
0.3
y (m)
0.3
0 0.3
0
x (m)
0.25
0.9
0 0.3
x (m)
0.9
0.5
Fig. 11. For a frequency f ¼5000 Hz with a simulated heterogeneous plate (SNR ¼40 dB), (a) the displacement field, (b)–(f) the evolutions of the identified local loss factor with the cut-off wavenumber, kc ¼ 50, 75, 100, 125 and 150 rad m 1 (the corresponding cut-off wavelength λc is indicated above each colormap). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
two regions with different properties. To illustrate the phenomenon, the inverse problem is applied for different values of cut-off wavenumber at two frequencies. The local loss factor is obtained from the noisy displacement field at the frequencies f ¼1500 Hz (Fig. 10) and f¼5000 Hz (Fig. 11) with five different levels of regularization, from kc ¼ 50 to 150 rad m 1 (the flexural wavenumbers at the considered frequencies are respectively kf ¼58 rad m 1 for f¼1500 Hz and kf ¼107 rad m 1 for f¼5000 Hz). Although the region with higher damping may be clearly localized on each figure, the value of cut-off wavenumber strongly influences the quality of the spatial reconstruction: at kf ¼50 rad m 1 (Figs. 10(b) and 11(b)), the contours of the inclusion are poorly reconstructed. At kf ¼75–100 rad m 1 (Fig. 10(c) and (d)) or kf ¼75–150 rad m 1 (Fig. 11(c) to (f)), the map reproduces the shape of the inclusion in a better way and the value of loss factor in each region may be inferred from the colorscale with satisfying accuracy. For f ¼1500 Hz and above kf ¼125 rad m 1 (Fig. 10(e) and (f)), the quantification of loss factor from the map is less accurate because of the noise, which is not removed efficiently.
2464
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
measurement area plate shaker impedance head
w
y
x F
scanning laser vibrometer
Fig. 12. Measurement set-up.
Two kinds of spurious patterns can also be distinguished on such maps. The first one, which does not depend on kc, consists in spurious lines in the map at the location of nodal lines in the displacement field (see e.g. Fig. 10(c) and (d) and Fig. 11(c)–(f)). This is explained by the fact that nodal lines of the displacement field are also nodal lines of the bilaplacian field (both fields are equal up to the constant multiplier D=ρhω2 , even once windowed and filtered). Thus, the estimation of ðD=ρhÞi;j from Eq. (9) consists in the ratio between two very small values. As these values are systematically below or comparable to the noise remaining after regularization, their ratio differs strongly from the actual value of ðD=ρhÞ at the considered points. The second kind of spurious pattern consists in “grain” in the map, at a size comparable to the cut-off wavelength λf. This is caused by the presence of spurious components at small wavelengths in the bilaplacian field, which are not contained in the displacement field and not removed by the filter (this phenomenon is well observable in Fig. 4(b)). 4. Experiments In the previous section, the identification technique has been validated on data from simulations. The purpose of the present section is to apply the procedure on real measurement data. The plate to test is suspended from a frame to approximate free boundary conditions (Fig. 12). A Polytec Scanning Vibrometer PSV 300 is used to measure the displacement field. The excitation is provided by a Brüel & Kjær 4810 shaker, supplied in power by a B&K 2718 amplifier. Although no knowledge of excitation level is required by the method, a B&K 8001 impedance head is used to provide both phase reference and input force measurement, in conjunction with a B&K NEXUS conditioning amplifier. 4.1. Homogeneous PMMA plate A first experiment was carried out on a polymethyl methacrylate (PMMA) plate. The plate dimensions are 0.5 m 0.4 m and thickness h¼3 mm. The material density is ρ¼1190 kg/m3. Because the material is transparent, one face of the plate was covered with silver paint to allow reflection of the laser beam. The displacements were measured over a meshgrid consisting of 71 71 ¼ 5041 points. The dimensions of the measurement area were 0.36 m 0.36 m, corresponding to a spatial step Δx ¼ Δy 5 mm. The plate was excited by the shaker at its bottom edge, outside the measurement area. The excitation signal was a periodic chirp in the frequency range [16, 6400] Hz. The displacement at each point was acquired using 10 averages to improve the signal-to-noise ratio. Young0 s modulus and loss factor identified by the method are shown in Fig. 13. For both properties, the dispersion around the mean curve is low. Regardless of small local variations inherent to the method, the evolution of these properties with frequency follows the same trend as the one observed by other authors for PMMA [16]. opt The optimum cut-off wavenumber kc does not appear to be proportional to the identified flexural wavenumber kf (Fig. 14), contrary to the results get by the simulations (see Fig. 7). This is due to the fact that the signal-to-noise ratio varies with frequency in the experiment. However, the optimum cut-off wavenumber is relatively high compared to the optimum opt wavenumber (kc =kN 1:7 in the range [2000, 6400] Hz), which confirms the good quality of the measurement. 4.2. Heterogeneous composite plates To assess the possibility to obtain a map of local properties from experimental data, two plates made by CETIM were used. One plate was made of PP thermoplastic matrix reinforced by glass fibers. It contains 10 patches of Santoprene™ elastomer. The other was made of a thermosetting epoxy matrix reinforced by glass fibers. It contained eight patches of an
10
0.1
8
0.08
6
0.06
2465
η
E (GPa)
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
4
0.04
2
0.02
0
0
1000 2000 3000 4000 5000 6000 f (Hz)
0
0
1000 2000 3000 4000 5000 6000 f (Hz)
Fig. 13. (a) Young0 s modulus E and (b) loss factor η identified on the PMMA plate by the method.
250
kf opt kc
k (rad/m)
200 150 100 50 0
0
1000
opt
Fig. 14. Optimum cut-off wavenumber kc
2000
3000 4000 f (Hz)
5000
6000
automatically found by the procedure, compared with the identified flexural wavenumber kf.
|wi, j |
y (m)
0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
(D/ ρh)i, j (m4 .s−2 )
ηi, j 1
0.3
40
0.3
0.2
30
0.2
20 0.1 0
10 0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
0
y (m)
y (m)
50
0.5
0.1 0
0
0.1
0.2
0.3
0.4 x (m)
0.5
0.6
0.7
0
Fig. 15. (a) Displacement field, (b) identified specific flexural rigidity, and (c) identified local loss factor. Plate made of glass fiber and PP matrix with 10 patches of Santoprene™ elastomer, f ¼6700 Hz, kc ¼150 rad m 1.
2466
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
Fig. 16. (a) Displacement field, (b) identified specific flexural rigidity, and (c) identified local loss factor. Plate made of glass fiber and epoxy matrix with eight patches of damping elastomer, f¼ 2000 Hz, kc ¼ 100 rad m 1.
elastomer specifically formulated so as to provide high damping. Both plates were made by co-curing the composite and elastomer. Contrary to classical damping treatments (e.g. constrained viscoelastic layer) which may be applied to a structure, the damping function is here directly embedded in the material. The plates were driven by a periodic chirp in the frequency range [500, 10,000] Hz. The excitation point was at the bottom edge, outside the scanning area. The displacements were measured on a fine mesh (around 20,000 points) with five averages. The measurement duration was in the order of several hours. On the first plate, the inverse problem was carried out with the displacement field obtained at frequency f¼6700 Hz on a mesh containing 211 101 ¼ 21; 311 points (Fig. 15(a)). The problem was regularized with cut-off wavenumber kc ¼150 rad m 1. The identified local specific flexural rigidity (Fig. 15(b)) clearly varies from the regions containing only composite material to those containing an elastomer core. Yet, no significant rise in the local loss factor is observed in the regions containing the elastomer (Fig. 15(c)). For the second plate, the displacements obtained at frequency f¼2000 Hz on a mesh containing 239 81 ¼ 19; 359 points were used in the inverse problem (Fig. 16(a)). A cut-off wavenumber kc ¼100 rad m 1 was chosen. Whereas the identified local specific flexural rigidity (Fig. 16(b)) does not vary significantly in the area considered, the local loss factor is much higher at the locations containing damping elastomer (Fig. 16(c)). It should be kept in mind that the test plates were made of composite materials. The inverse method should then be used with an a priori equation of motion of an orthotropic plate and the test should be done with several sources in order to generate waves in all directions. It is then particularly interesting to see that the fact to consider an equivalent isotropic plate here is sufficient for the understanding of the behavior of the inserted patches. The two different behaviors, modification of the local flexural stiffness (Fig. 15(b)) or the local damping (Fig. 16(c)), can be attributed to the quality of the adhesion between the patches and the composite layers, and to the material properties of each elastomer.
5. Conclusion In this paper, an adaptation of the Force Analysis Technique has been developed to identify stiffness and damping properties of plates. If the general form of the equation of motion is known, it is possible to identify structural parameters by measuring the vibratory displacement field over an area in which no external force is applied. While keeping the overall regularization procedure as in the Force Analysis Technique, i.e., windowing and filtering, the method has been improved by a novel automatic determination of the optimum cut-off wavenumber. Contrary to modal methods, the proposed approach is independent of boundary conditions and may be applied at any frequency, not necessarily a resonance. It is also valid in the mid-frequency domain, where the modal overlap is high. Moreover, it can be applied in the case of heterogeneous structures to obtain a map of local properties. The applicability of the method has been verified on simulated plates and experimental analysis on real plates. In this paper, the equation of motion of an isotropic plate has been considered. Following the same principle, the method can be applied to beams. A natural continuation of this work is to consider orthotropic plates, which are commonly encountered with composite materials.
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
2467
Acknowledgments This study was funded by Région Pays de la Loire within the framework of the AMORTI project. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the manuscript, and Alain Lemascon for his careful reading of the manuscript. Appendix A. Finite differences scheme The fourth-order partial derivatives of the displacement field involved in the local equation of motion are calculated as follows: ∂4 w 1 δ4x i;j ¼ 4 wi þ 2;j 4wi þ 1;j þ 6wi;j 4wi 1;j þwi 2;j ∂x4 Δx ∂4 w 1 δ4y i;j ¼ 4 wi;j þ 2 4wi;j þ 1 þ 6wi;j 4wi;j 1 þwi;j 2 ∂y4 Δy ∂4 w 1 δ2x2y ¼ 2 2 wi þ 1;j þ 1 2wi þ 1;j þ wi þ 1;j 1 2wi;j þ 1 þ 4wi;j 2wi;j 1 þwi 1;j þ 1 2wi 1;j þwi 1;j 1 ; i;j ∂x2 ∂y2 Δx Δy
(A.1)
where Δx and Δy are the spatial steps in x and y directions. It clearly appears from these expressions that the spatial derivatives cannot be obtained at points close to the edges of the observation area. Thus the effective area in which the inverse problem can be applied only contains ðM x 4Þ ðM y 4Þ points. Appendix B. Expression of the spatial window and filter response In the regularization step of the inverse method, each field wi;j and ∇4i;j is first multiplied by a spatial Tukey window ψ 2D i;j . As the exploitable area for the inverse problem is delimited by i A ½3; M x 2 and j A ½3; M y 2, the window is implemented as follows: 1D ψ 2D ðxi;j x1;1 2Δx Þ ψ 1D ðyi;j y1;1 2Δy Þ; i;j ¼ ψ
where
8 πx > > 0:5 1 cos > > α > > > <1 1D ψ ðxÞ ¼ πðx L þ 2αÞ > > 0:5 1 cos > > α > > > : 0
(B.1)
if 0 r x o α; if α rx rL α; if L α ox r L;
(B.2)
else
is the expression of the unidimensional Tukey window, L ¼ Lx 4Δx and 8 < λc if L Z 2λc ; α¼ L else: : 2
(B.3)
Eq. (B.3) expresses the fact that the Tukey window is replaced by a Hanning window when the cut-off wavelength becomes too large compared to the size of the observation area. Once windowed, each field is convolved with the discrete spatial response hi;j ¼ hðxi;j ; yi;j Þ of a low-pass filter calculated using the following expression: 8
2 > kc kc x kc y 2π 2π < 1 þ cos sin k 1 þ cos ð x Þ sin ð k y Þ if x and y A ; ; c c 2 2 kc kc (B.4) hðx; yÞ ¼ 4π 2 xyN f > : 0 else; where Nf is a normalization parameter, chosen such that ∬ hðx; yÞ dx dy ¼ 1: References [1] D. Ewins, Modal Testing: Theory and Practice, Research Studies Press, Baldock, United Kingdom, 1984. [2] M. Matter, T. Gmuer, J. Cugnoni, A. Schorderet, Numerical-experimental identification of the elastic and damping properties in composite plates, Composite Structures 90 (2009) 180–187. [3] J. McDaniel, W. Shepard, Estimation of structural wave numbers from spatially sparse response measurements, Journal of the Acoustical Society of America 108 (2000) 1674–1682. [4] J. Berthaut, M. Ichchou, L. Jezequel, K-space identification of apparent structural behaviour, Journal of Sound and Vibration 280 (2005) 1125–1131. [5] M. Rak, M. Ichchou, J. Holnicki-Szulc, Identification of structural loss factor from spatially distributed measurements on beams with viscoelastic layer, Journal of Sound and Vibration 310 (2008) 801–811.
2468
F. Ablitzer et al. / Journal of Sound and Vibration 333 (2014) 2454–2468
[6] M.N. Ichchou, O. Bareille, J. Berthaut, Identification of effective sandwich structural properties via an inverse wave approach, Engineering Structures 30 (2008) 2591–2604. [7] C. Pezerat, J.L. Guyader, 2 Inverse methods for localization of external sources exciting a beam, Acta Acustica 3 (1995) 1–10. [8] C. Pezerat, J.L. Guyader, Force analysis technique: reconstruction of force distribution on plates, Acustica 86 (2000) 322–332. [9] C. Pezerat, J.L. Guyader, Identification of vibration sources, Applied Acoustics 61 (2000) 309–324. [10] C. Pezerat, Q. Leclere, N. Totaro, M. Pachebat, Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique, Journal of Sound and Vibration 326 (2009) 540–556. [11] Q. Leclere, C. Pezerat, Vibration source identification using corrected finite difference schemes, Journal of Sound and Vibration 331 (2012) 1366–1377. [12] M. Djamaa, N. Ouelaa, C. Pezerat, J.L. Guyader, Mechanical radial force identification of a finite cylindrical shell by an inverse method, Acta Acustica united with Acustica 92 (2006) 398–405. [13] M.C. Djamaa, N. Ouelaa, C. Pezerat, J.L. Guyader, Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering, Journal of Sound and Vibration 301 (2007) 560–575. [14] C. Renzi, C. Pezerat, J.L. Guyader, Vibratory source identification by using the Finite Element Model of a subdomain of a flexural beam, Journal of Sound and Vibration 332 (2013) 545–562. [15] J.L. Guyader, Vibration in Continuous Media, ISTE Ltd, London, United Kingdom, 2006. [16] Y. Liao, V. Wells, Estimation of complex modulus using wave coefficients, Journal of Sound and Vibration 295 (2006) 165–193.