1 Introduction

The rapid assessment of the elastic properties of materials has never been more prevalent than today, due to the acceleration in material development. Advancement of material technology has expanded the number of materials and applications for which these materials can be used. In particular, carbon fibre reinforced polymers (CFRP) have given way to a wide variety of complex materials with material characteristics that vary according to use. For example, CFRP applications range from a very specific highly technology environments (aerospace, Formula 1 and nuclear) to mass market applications (sport equipment, automotive and energy applications). At the same time the advent of commercial 3D printing has added a wide variety of materials that can be customised for a specific use, with highly varying material properties. A further area of interest is the characterisation of biological material where structures can be both anisotropic and heterogeneous making it difficult to assess these materials with traditional engineering methods. The speed and rate of material development has also been enhanced due to ever increasing automation. Thus it is vital to evaluate the material properties (Elasticity, Yield, Hardness, etc.) of these components before they are used in the field.

The focus of this paper is to provide an alternative method for the measurement of the out-of-plane Young’s Modulus (Elastic Modulus - E). The Young’s Modulus of a material defines the relationship between stress and strain, and provides a direct measure of material stiffness. A number of techniques have been used for the determination of the Young’s Modulus, using both destructive and non-destructive testing procedures. These methods include: resonant frequency (in-plane) [1], instrumental indentation [2, 3], traditional engineering methods (tensile tests) and ultrasound techniques [4].

Resonant frequency testing has focused on the evaluation of the Young’s Modulus with regard to thin film materials (in recent times), which generally focus on the calculation of the Young’s modulus of a micro-lever using the value of resonant frequencies obtained in prior experiments, although requiring both a numerical modal analysis as well as evaluation of experimental results. While instrumental indentation is a destructive method that infers from an indentations size and depth the Young’s Modulus.

Traditional engineering methods are very common, highly controlled testing techniques, which generally focus on: tensile, compression and three point bending tests. These tests are highly suited for isotropic materials, as a tensile test can only provide one Young’s modulus (in the direction of tension) but determination of the out-of-plane Young’s Modulus becomes very difficult when investigating anisotropic materials. For example composite materials; as these structures are relatively thin in the out-of-plane direction, further compounding the issue is the manufacturing process which requires individual plies (± 0.25 mm per ply) to be placed on top of one another in the out-of-plane direction, thus requiring an enormous number of plies to be placed to meet the requirements for a tensile test in the out-of-plane direction.

Ultrasound methods provide various advantages over the previously mentioned techniques as smaller specimens can be used with less complicated geometries than mechanical testing methods [5]. Furthermore, the tests are non-destructive, thus more than one elastic constant can be evaluated per specimen. These methods often rely on the determination of the speed of sound through the specimen from which the elastic properties can be obtained [6].

This paper proposes a novel impedance based method to determine the out-of-plane Young’s Modulus of anisotropic materials. The method is based on the evaluation of the acoustic impedance (Z) of a material with the use of a two microphone impedance tube. Impedance tubes are generally used for characterisation of sound absorbing materials, such as wool or polymers, used within building and transport applications.

An impedance tube is a simple apparatus that allows for a plane wave to be generated and propagate down a pipe, this acoustic wave then reflects off test materials. The phase interference between the waves in the pipe, which are incident upon and reflected from the test sample result in the formation of a standing wave pattern in the pipe. Under conditions of no absorption, the incident and reflected waves have the same amplitude resulting in nodes of zero pressure and antinodes with double the pressure. Whereas, if there is absorption by the test sample, the incident and reflected waves have different amplitudes and the nodes do not have zero pressure. The pressure amplitudes at nodes and antinodes can be measured with microphone probes. These measurements can be used to calculate a materials reflection coefficient (R), its absorption coefficient (α), and its impedance (Z) [7, 8].

One of the most common ways to determine the acoustic impedance of materials is with ultrasound, where the impedance is directly related to the longitudinal speed of sound through the material. Recently various authors have developed alternative techniques based on ultrasound in order to determine: the acoustic impedance [9, 10], non-acoustic properties of materials [11] and the Young’s Modulus [12, 13]. One of the issues with ultrasound testing of impedance and Young’s Modulus testing is its ability to evaluate materials with high structural damping or porosity. Acoustic impedance and Young’s Modulus analysis testing requires the accurate determination of the speed of sound through a known thickness. Porous materials result in the scattering of the longitudinal wave which makes it difficult to identify the back wall echo and thus accurately determine the impedance or Young’s Modulus of these types of materials (foams, rubber etc.).

The proposed method provides a solution that is applicable to a wide range of materials from soft foams to hard composite materials, and thus provides an advantage over traditional techniques such as ultrasound for the determination of the out-of-plane Young’s Modulus.

The paper is organized in four sections. The first section of the paper describes the principles and basic assumptions of the proposed method. The method is then numerically validated in the second section. A coupled structural-acoustic FE model is developed including a RVE model of composite material as the tested sample for which the elastic properties are well-know from the literature. The proposed method is applied to the simulation data in a further post-processing step, where the elastic modulus is estimated and compared with the expected one. An experimental measurement campaign on different materials is then reported in the last section of the paper. The experimental setup and the samples used in the experimental campaign are presented. Results evaluated with the proposed method for the tested materials are finally discussed and compared with those given by standard method (tensile testing and ultrasound time of flight).

2 Elastic Modulus Estimation Method

The proposed method is based on the two microphone impedance tube because a plane standing wave condition is required. Starting from the acoustic pressure measurements in two different locations through the pipe, the Transfer Function Method is applied to estimate the reflection coefficient and acoustic impedance. It is demonstrated how the acoustic impedance can be related to elastic properties of the material including the propagation speed of sound.

2.1 Impedance Tube

According to the ASTM E 1050 [14], a two microphone impedance tube is expressed as a one-dimensional wave guide, with a sound source at one end and the acoustic load (typically the test sample) placed at the other end which guarantees a plane standing wave into the duct. In order to maintain plane wave propagation into the duct, the characteristic dimension of the tubes cross section satisfy the following condition:

$$ \mathrm{d}<\frac{\mathrm{K}\ \mathrm{c}}{{\mathrm{f}}_{\mathrm{u}}} $$
(1)

where K is a geometrical parameter associated to the cross section shape (0.5 for rectangular or square cross section), c the speed of sound and fu the upper frequency limit. In order for a plane wave to fully develop before reaching the pressure acquisition points within the tube a minimum of 3d must be allowed between sound source and the nearest pressure acquisition point. The acquisition points spacing was prescribed to be less than 80% of shortest half wave length of interest

$$ \mathrm{s}<0.4\frac{\ c}{{\mathrm{f}}_{\mathrm{u}}} $$
(2)

2.2 Transfer Function Method

The Transfer Function Method for acoustic properties measurements was proposed in 1980 by Chung and Blaser [15, 16]. The theory behind this method involve the decomposition of a standing wave into its incident and reflected components using a simple transfer-function relation between the acoustic pressures measured in two different locations as shown in Fig. 1.

Fig. 1
figure 1

Impedance Tube configuration with two microphones and a rigid back wall

Full size image

According to Fig. 1, the complex sound pressure at the microphone positions can be calculated as

$$ {p}_1={\widehat{P}}_i{e}^{\left(-j{k}_0{x}_1\right)}+{\widehat{P}}_r{e}^{\left(j{k}_0{x}_1\right)} $$
(3)
$$ {p}_2={\widehat{P}}_i{e}^{\left(j{k}_0{x}_2\right)}+{\widehat{P}}_r{e}^{\left(-j{k}_0{x}_2\right)} $$
(4)

where \( {\widehat{P}}_i \)\( {\widehat{P}}_r \) are the amplitude of the incident and reflecting wave respectively, k0 the wave number in air and x1 and x2 are the microphones coordinates according to the considered reference system.

The reflection coefficient (R) can be expressed through the transfer function (H12) between the complex pressures at two microphone positions and the transfer function of the incident (Hi) and reflecting waves (Hr):

$$ R=\frac{H_{12}-{H}_i}{H_r-{H}_{12}}{e}^{\left(2j{k}_0\left(L+s\right)\right)} $$
(5)

Where s is the microphone spacing and L is the distance between the sample and the closed microphones. Moreover, the transfer function of the incident and reflecting wave at the two microphone positions can be expressed as:

$$ {H}_i=\frac{P_i\left({x}_2\right)}{P_i\left({x}_1\right)}=\frac{{\widehat{P}}_i{e}^{\left(-j{k}_0{x}_2\right)}}{{\widehat{P}}_i{e}^{\left(-j{k}_0{x}_1\right)}}={e}^{\left(j{k}_0s\right)} $$
(6)
$$ {H}_r=\frac{P_r\left({x}_2\right)}{P_r\left({x}_1\right)}=\frac{{\widehat{P}}_r{e}^{\left(j{k}_0{x}_2\right)}}{{\widehat{P}}_r{e}^{\left(j{k}_0{x}_1\right)}}={e}^{-\left(j{k}_0s\right)} $$
(7)

where Pi and Pr are the complex pressures of the incident and reflecting waves.

2.3 Impedance Calibration

Considering that, the Transfer Function Method for the post-processing of data is a complex ratio between the acoustic pressure responses measured by the two microphones system, any mismatch on the amplitude and phase will affect the accuracy of the measurements. Therefore, a correction factor was be defined. So at the first step of the calibration procedure the mics are placed in a rigid plate located at the end of the tube. In this case they were exposed to the same sound pressure and phase, so the ratio of the fast Fourier transform of the two time signals for the two microphones represent a correction factor (Hc) for the amplitude and phase mismatches

$$ {H}_c=\frac{fft\left({p}_2\right)}{fft\left({p}_1\right)} $$
(8)

which will be applied to the complex measured transfer function between the sound pressures measured by the microphones placed in the standard position (Fig. 1).

$$ {H}_{12}=\frac{{\widehat{H}}_{12}}{H_c} $$
(9)

where \( {\widehat{H}}_{12}= fft\left({\widehat{p}}_{mic1}\right)/ fft\left({\widehat{p}}_{mic2}\right) \) is the uncorrected (measured) transfer function.

Moreover, the thermal viscous losses, which induce attenuation on the propagating incident and reflected waves within the tube are not negligible [14]. In order to take into account tube attenuation, a complex wave number must be introduced where the imaginary part represent the attenuation constant.

$$ {\mathrm{k}}_0=\mathrm{k}-\mathrm{j}{\mathrm{k}}^{{\prime\prime} } $$
(10)

If the real part is linearly dependent on the frequency and inversely proportional to the speed of sound (k = 2πf/c), the imaginary part is not a priori known, except to empirical relationship [14, 17].

Therefore the attenuation constant can be measured [17].

Considering an acoustically rigid and high reflection material sample and measuring the transfer function (corrected for the amplitude and phase mismatches), the Eq. (5) can be solved assuming a reflection coefficient (R) of 1 and including a complex wave number

$$ {\mathrm{H}}_{12}=\frac{{\mathrm{e}}^{-\mathrm{j}{\mathrm{k}}_0\mathrm{L}}+{\mathrm{e}}^{\mathrm{j}{\mathrm{k}}_0\mathrm{L}}}{{\mathrm{e}}^{-\mathrm{j}{\mathrm{k}}_0\left(\mathrm{L}+\mathrm{s}\right)}+{\mathrm{e}}^{\mathrm{j}{\mathrm{k}}_0\left(\mathrm{L}+\mathrm{s}\right)}}=\frac{\mathit{\cos}\left[\left(\mathrm{k}-\mathrm{j}{\mathrm{k}}^{{\prime\prime}}\right)\mathrm{L}\right]}{\mathit{\cos}\left[\left(\mathrm{k}-\mathrm{j}{\mathrm{k}}^{{\prime\prime}}\right)\left(\mathrm{L}+\mathrm{s}\right)\right]} $$
(11)

Equation (11) can be solved numerically for k” by means of a Newton-Raphson iteration scheme and assuming an initial guess the empirical value

$$ {\mathrm{k}}_0^{\hbox{'}\hbox{'}}=\frac{A\sqrt{\mathrm{f}}}{\mathrm{c}\ \mathrm{d}} $$
(12)

where d is the tube diameter and A is a constant (A = 0.0194 [18]).

When the imaginary part is estimated, the complex wavenumber can be included into the Eq. (5) for the experimental data post-processing.

2.4 Elastic Modulus Estimation

In this section the proposed method is presented starting from the sound wave propagation at the interface of two different media.

Consider the air-solid interaction shown in Fig. 2, with a plane incident sound wave. Due to the difference in acoustic properties between two media, the energy associated to the incident wave will be in part reflected and in part transmitted at the interface.

Fig. 2
figure 2

System of reflection-transmission at air-solid interaction

Full size image

At the interface three boundary conditions must be satisfied. First, the particle velocity normal to the bounding surface and the pressure variation must be continuous. Continuity of normal stress must be guaranteed and finally zero tangential stress is required since the fluid cannot support viscous stress (assuming an ideal non-viscous liquid). These boundary conditions can be written in term of velocity potentials, which can be expressed as plane wave solutions to the wave equation [19, 20] and solved in terms of reflection and transmission coefficients. In particular, the reflection coefficient can be written as

$$ R\left(\theta \right)=\frac{Z_L co{s}^22{\theta}_S+{Z}_S si{n}^22{\theta}_S-{Z}_1}{Z_L co{s}^22{\theta}_S+{Z}_S si{n}^22{\theta}_S+{Z}_1} $$
(13)

where θS is the angle of the propagating shear waves.

Since in a solid the wave propagation is associated to longitudinal and shear waves, ZL and ZS represent the longitudinal and shear characteristic impedance respectively, while Z1 is the characteristic impedance of the fluid

$$ {Z}_1=\frac{\rho_1{V}_1}{\mathit{\cos}{\theta}_i},{Z}_L=\frac{\rho_2{V}_L}{\mathit{\cos}{\theta}_L},{Z}_S=\frac{\rho_2{V}_S}{\mathit{\cos}{\theta}_S} $$
(14)

So the characteristic impedance is a function of the density (ρ1 density of fluid - ρ2 density of solid media) and the relative propagation velocity (V1 speed of sound in the fluid, VLspeed of propagating longitudinal waves, VS speed of propagating shear waves).

Because of longitudinal and shear wave are involved, an effective impedance (Zeff) can be defined

$$ {Z}_{eff}={Z}_L co{s}^2\left(2{\theta}_s\right)+{Z}_S si{n}^2\left(2{\theta}_s\right) $$
(15)

so that the reflection coefficient becomes

$$ R\left(\theta \right)=\frac{Z_{eff}-{Z}_1}{Z_{eff}+{Z}_1} $$
(16)

The geometrical laws of reflection of sound waves are the same as those that apply to light waves unless the length of the sound waves is not comparable with the linear dimensions of the reflecting objects [21]. As the wavelength is small relative to the dimension of the reflectors, the Snell’s law can be applied

$$ \frac{{\mathrm{sin}\uptheta}_{\mathrm{i}}}{{\mathrm{V}}_{\mathrm{i}}}=\frac{{\mathrm{sin}\uptheta}_{\mathrm{S}}}{{\mathrm{V}}_{\mathrm{S}}} $$
(17)

In the particular case of plane and normal incident waves, the incidence angle is zero (θi=0), so according to Eq. (17) results in reflection and transmission angles of zero and only longitudinal waves are transmitted into the media, so from Eq. (16) the reflection coefficient can be simplified as follows:

$$ \mathrm{R}=\frac{{\mathrm{Z}}_2-{\mathrm{Z}}_1}{{\mathrm{Z}}_2+{\mathrm{Z}}_1} $$
(18)

where Z1 is the characteristic impedance of fluid (Z1 = ρf Cf, with Cf speed of sound in the fluid) and Z2 is characteristic impedance of solid (Z2 = ρsCLs) which is a function of the longitudinal sound speed in the solid.

The speed of sound describes the speed of sound waves passing through an elastic medium and varies with the medium employed as well as with the properties of the medium, in particular it is closely related to the stiffness of the material. When considering solid materials, there is non-zero stiffness both for volumetric and shear deformations. Hence, in a solid it is possible to generate sound waves with different velocities dependent on the deformation mode. When the lateral dimension (thickness) of the medium is much smaller than the incident wavelength the speed of sound can be expressed by the relation

$$ {\mathrm{C}}_{\mathrm{Ls}}=\sqrt{\frac{\mathrm{E}}{\uprho}} $$
(19)

Considering a normal incident plane wave, as we demonstrated only longitudinal waves will pass through the solid, so the speed of sound in Eq. (19) represents the longitudinal speed of sound. Moreover the Elastic Modulus included in Eq. (19) represents the out of plane component of elastic properties of the material (E3).

Now combining the Eqs. (18) and (19) and assuming the acoustic impedance under plane incident wave condition (Z2 = ρsCLs), the Elastic Modulus of material can be expressed as:

$$ E=\frac{1}{\rho_s}\frac{{\left(1+R\right)}^2}{{\left(1-R\right)}^2}{Z}_1 $$
(20)

In conclusion the elastic modulus of the material can be estimated using acoustic measurement. Starting from the sound pressure measurements in two different position through an impedance tube and applying the transfer function method the acoustic impedance and the reflection coefficient can be measured and the elastic modulus can be estimated then through the Eq. (20).

3 Numerical Simulation

In order to validate the proposed method, a numerical model, based on the Representative Volume Element (RVE), was developed. An RVE model of composite material for which the elastic properties were known from literature, was developed and included in the FE simulation.

First a coupled acoustic-structural harmonic analysis was considered in order to simulate the impedance tube test rig. The simulation data in terms of sound pressures acquired during the simulation was used as input data and the elastic modulus was estimated by the proposed method and compared with the expected one.

Moreover a static structural analysis was implemented using the same RVE model in order to simulate a standard tensile testing and estimate the elastic modulus according to the Hooke’s law. The results was then compared with the expected one and with the elastic modulus extracted by the proposed method.

3.1 Representative Volume Element (RVE)

The Representative Volume Element (RVE) [22,23,24,25] is a powerful and commonly used approach for numerical modelling and elastic behaviour analysis of composite materials and general periodic structures [26,27,28,29]. Different RVE models have been developed, but in this work the Homogenization Theory was taken into account which was shown to give more accurate estimates of effective stiffness for periodic composite material [26].

The homogenization theory is based on two assumptions. The fields vary on multiple scales due to the existence of microstructure, and the microstructure is assumed spatially periodic. Starting from such assumption the field variables are approximated by an asymptotic expansion

$$ {u}_i^{\eta}\left({x}_i,{y}_i\right)={u}_{0i}\left({x}_i,{y}_i\right)+\eta {u}_{1i}\left({x}_i,{y}_i\right)+{\eta}^2{u}_{2i}\left({x}_i,{y}_i\right)+\dots $$
(21)

where \( {\mathrm{u}}_{\mathrm{i}}^{\upeta} \)is the exact value of the field variable, u0i is the macroscopic or average value of the field variable, u1i, u2i, etc. are perturbations in the field variables due to the microstructure, xi are the global level coordinates, yi are the micro level coordinates, and η the ratio of the microstructure size to the total size of the analysis region. In the elasticity theory u0i would be the continuum level displacements while u1i would be the microstructural displacements. The relation between the macro and micro coordinate can be expressed

$$ {y}_i=\frac{x_i}{\eta } $$
(22)

Applying the derivatives with respect to xi to the expansion of the displacement (Eq. (21)) and neglecting the higher order terms of the weak form of equilibrium equations governing the mechanical behaviour of composite material, the small deformation strain tensor can be written as

$$ {\varepsilon}_{ij}(u)={\overline{\varepsilon}}_{ij}+{\varepsilon}_{ij}^{\ast } $$
(23)

where

$$ {\overline{\varepsilon}}_{ij}(u)=\frac{1}{2}\left[\frac{\partial {u}_{0i}}{\partial {x}_j}+\frac{\partial {u}_{0j}}{\partial {x}_i}\right]\kern0.5em and\kern0.5em {\varepsilon}_{ij}^{\ast }(u)=\frac{1}{2}\left[\frac{\partial {u}_{1i}}{\partial {y}_j}+\frac{\partial {u}_{1j}}{\partial {y}_i}\right] $$
(24)

εij is the microstructural strain tensor, \( {\overline{\varepsilon}}_{ij} \) is the macroscopic strain tensor and \( {\varepsilon}_{ij}^{\ast } \) is the fluctuating microstructural strain tensor which is assumed to vary periodically.

In the same way, the virtual strain tensor can be expressed

$$ {\varepsilon}_{ij}(v)={\varepsilon}_{ij}^0+{\varepsilon}_{ij}^1 $$
(25)

where

$$ {\varepsilon}_{ij}^0(v)=\frac{1}{2}\left[\frac{\partial {v}_{0i}}{\partial {x}_j}+\frac{\partial {v}_{0j}}{\partial {x}_i}\right]\kern0.5em \mathrm{and}\kern0.5em {\varepsilon}_{ij}^1(v)=\frac{1}{2}\left[\frac{\partial {v}_{1i}}{\partial {y}_j}+\frac{\partial {v}_{1j}}{\partial {y}_i}\right] $$
(26)

Assuming the virtual displacement varying only on microscopic level and is constant on the macroscopic level (\( {\varepsilon}_{ij}^0(v) \)=0), and assuming η goes to zero in the limit, then the standard weak form of equilibrium equation can be written as

$$ \underset{\varOmega }{\int}\frac{1}{V_{RVE}}\underset{V_{RVE}}{\int }{D}_{ij kl}\left({\varepsilon}_{ij}^1(v)\right)\left({\overline{\varepsilon}}_{kl}(u)+{\varepsilon}_{kl}^{\ast }(u)\right)d{V}_{RVE} d\varOmega =0 $$
(27)

where Ω is the total of macroscopic and microscopic domain of composite material, VRVE is the volume of RVE and Dijkl is the stiffness tensor. Equation (27) is satisfied if the integral over RVE is zero, so the equilibrium equations may be written

$$ \underset{V_{RVE}}{\int }{D}_{ij kl}{\varepsilon}_{ij}^1(v){\varepsilon}_{kl}^{\ast }(u)d{V}_{RVE}=-\underset{V_{RVE}}{\int }{D}_{ij kl}{\varepsilon}_{ij}^1(v){\overline{\varepsilon}}_{kl}(u)d{V}_{RVE} $$
(28)

For linear problem in 3-dimensions, \( {\overline{\varepsilon}}_{kl}(u) \), which is unknown apriori, can be written as a linear combination of unit strains as following

$$ {\displaystyle \begin{array}{c}{\overline{\varepsilon}}_{pm}^{11}=-\left[\begin{array}{ccc}1& 0& 0\\ {}0& 0& 0\\ {}0& 0& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{12}=-\left[\begin{array}{ccc}0& 1& 0\\ {}0& 0& 0\\ {}0& 0& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{13}=-\left[\begin{array}{ccc}0& 0& 1\\ {}0& 0& 0\\ {}0& 0& 0\end{array}\right]\\ {}{\overline{\varepsilon}}_{pm}^{21}=-\left[\begin{array}{ccc}0& 0& 0\\ {}1& 0& 0\\ {}0& 0& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{22}=-\left[\begin{array}{ccc}0& 0& 0\\ {}0& 1& 0\\ {}0& 0& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{23}=-\left[\begin{array}{ccc}0& 0& 0\\ {}0& 0& 1\\ {}0& 0& 0\end{array}\right]\\ {}{\overline{\varepsilon}}_{pm}^{31}=-\left[\begin{array}{ccc}0& 0& 0\\ {}0& 0& 0\\ {}1& 0& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{32}=-\left[\begin{array}{ccc}0& 0& 0\\ {}0& 0& 0\\ {}0& 1& 0\end{array}\right],{\overline{\varepsilon}}_{pm}^{33}=-\left[\begin{array}{ccc}0& 0& 0\\ {}0& 0& 0\\ {}0& 0& 1\end{array}\right]\end{array}} $$
(29)

Applying the symmetry of strain states and substituting the unit strains on the right hand side of Eq. (28), the stress tensor can be derived as

$$ {\sigma}_{ij}^{\ast kl}={D}_{ij pm}{\overline{\varepsilon}}_{pm}^{kl} $$
(30)

A set of six auxiliary problems have to be solved for \( {\varepsilon}_{kl}^{\ast kl} \) in Eq. (28) using Eq. (30)

$$ \underset{V_{RVE}}{\int }{D}_{ij kl}{\varepsilon}_{ij}^1(v){\varepsilon}_{pm}^{\ast kl}(u)d{V}_{RVE}=\underset{V_{RVE}}{\int }{\varepsilon}_{ij}^1(v){\sigma}_{ij}^{\ast kl}d{V}_{RVE} $$
(31)

To ensure periodicity of the strain field \( {\varepsilon}_{pm}^{\ast kl} \) in Eq. (28) the value of \( {\varepsilon}_{ij}^{\ast } \) for any arbitrary \( {\overline{\varepsilon}}_{kl} \) can be written as

$$ {\varepsilon}_{ij}=-{\varepsilon}_{ij}^{\ast kl}{\overline{\varepsilon}}_{kl} $$
(32)

So the relation between local RVE strain and the average strain can be written

$$ {\varepsilon}_{ij}={M}_{ij kl}{\overline{\varepsilon}}_{kl}\kern0.5em with\kern0.5em {M}_{ij kl}=\frac{1}{2}\left({\delta}_{ik}{\delta}_{il}+{\delta}_{il}{\delta}_{jk}\right)-{\varepsilon}_{ij}^{\ast kl} $$
(33)

Finally the effective stiffness tensor which relates average stress and average strain can be estimated from Mijkl by writing the Hooke’s law at microscopic level and then integrating over the volume of RVE and divide by the volume of RVE

$$ {\overline{D}}_{ijkl}=\frac{1}{\mid {V}_{RVE}\mid}\underset{V_{RVE}}{\int }{D}_{ijpm}{M}_{pmkl}d{V}_{RVE} $$
(34)

3.2 Static Structural Analysis

A static structural analysis was used to extract a preliminary indication of the elastic modulus of the sample. Considering an RVE approach, and assuming the material to behave as a linear elastic solid, a constant axial loading and fixed support boundary conditions was applied on the sample (Fig. 3). The two ends of the RVE along the transverse direction are marked as a fixed area and load area respectively. A displacement boundary condition was prescribed at the fixed area (ux = uy = uz = 0), while a force was applied as the loading condition to the load area. The applied load is confined in the elastic range, so measuring the stress induced by the loading conditions and the relative strain status, the Hooke’s law can provide an estimation about the elastic modulus of the sample:

Fig. 3
figure 3

Loading and boundary condition applied on Representative Volume Element for the static structural analysis

Full size image
$$ {\mathrm{E}}_3=\frac{\upsigma_3}{\upvarepsilon_3} $$
(35)

Where σ3 and ε3 are the tensile stress and strain respectively along the loading direction.

In order to obtain the averaged stress-strain results for a certain prescribed load and boundary conditions, periodic boundary condition [30,31,32,33,34] need to be applied on the RVE. The periodic boundary conditions were implemented by creating identical mesh on opposite faces of RVE to ensure the same coordinates are tangential to the face.

3.3 Coupled Acoustic-Structural Analysis

Then a coupled acoustic structural analysis was carried to simulate in-silico the impedance tube measurements and verify the proposed approach numerically. In acoustic fluid-structural interaction (FSI) problems, the structural dynamic equation (Eq. (36)) must be considered along with the Navier-Stokes equations of fluid momentum (Eq. (37)) and the flow continuity equation (Eq. (38).

$$ \left[M\right]\left\{\overset{..}{u}(t)\right\}+\left[C\right]\left\{\overset{.}{u}(t)\right\}+\left[K\right]\left\{u(t)\right\}=\left\{f(t)\right\} $$
(36)
$$ \rho \frac{dv}{dt}=-\nabla p+\nabla \cdot S+\rho b $$
(37)
$$ \frac{\partial \rho }{\partial t}=-\nabla \cdot \left(\rho v\right)+Q $$
(38)

Where ρ is the density, v is the velocity vector, Q is mass source, S is the viscous stress tensor, p the pressure, b is the body force and [M], [C], [K] are the mass, damping, stiffness matrices respectively. Combining the Navier-Stokes and the continuity equations, the acoustic wave equation can be written as

$$ \nabla \cdot \left(\frac{1}{\rho_0}\nabla p\right)-\frac{1}{\rho_0{c}^2}\frac{\partial^2p}{\partial {t}^2}+\nabla \cdot \left[\frac{4\mu }{3{\rho}_0}\nabla \left(\frac{1}{\rho_0{c}^2}\frac{\partial p}{\partial t}\right)\right]=-\frac{\partial }{\partial t}\left(\frac{Q}{\rho_0}\right)+\nabla \cdot \left[\frac{4\mu }{3{\rho}_0}\nabla \left(\frac{Q}{\rho_0}\right)\right] $$
(39)

where μ is the dynamic viscosity and ρ0 is the mean fluid density.

The discretized structural equation (Eq. (36)) and the acoustic wave equation (Eq. (40)) must be considered simultaneously at the FSI interface. So it is convenient writing the equation in a matrix form. The finite element formulation is obtained by testing wave equation (Eq. (39)) using the Galerkin procedure [35]. Equation (39) is multiplied by testing function w and integrated over the volume of the domain [36] which yields the following:

$$ {\displaystyle \begin{array}{c}{\iiint}_{\varOmega_f}\frac{1}{\rho_0{c}^2}w\frac{\partial^2p}{\partial {t}^2} dv+{\iiint}_{\varOmega_f}\nabla \left(\frac{4\mu }{3{\rho}_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right) dv+{\iiint}_{\varOmega_f}\nabla w\cdot \left(\frac{1}{\rho_0}\nabla p\right) dv+\\ {}-{\iint}_{\varGamma_f}w\left(\frac{1}{\rho_0}+\frac{4\mu }{3{\rho}_0^2{c}^2}\frac{\partial }{\partial t}\right)\hat{n}\cdot \nabla pds+{\iint}_{\varGamma_f}w\frac{4\mu }{3{\rho}_0^2}\hat{n}\cdot \nabla Qds=\\ {}={\iiint}_{\varOmega_f}w\frac{1}{\rho_0}\frac{\partial Q}{\partial t} dv+{\iiint}_{\varOmega_f}\nabla w\cdot \left(\frac{4\mu }{3{\rho}_0^2}\nabla Q\right) dv\end{array}} $$
(40)

where dv is the volume differential of acoustic domain (Ωf), ds is the surface differential of acoustic domain boundary (Γf) and \( \widehat{n} \) is the normal unit vector to the boundary.

From the equation of momentum conservation, the normal acceleration of the fluid, expressed in terms of the normal displacement, is given

$$ \frac{\partial {v}_{n,f}}{\partial t}=\widehat{n}\cdot \frac{\partial v}{\partial t}=\widehat{n}\cdot \frac{\partial^2{u}_f}{\partial {t}^2}=-\left(\frac{1}{\rho_0}+\frac{4\mu }{3{\rho}_0^2{c}^2}\frac{\partial }{\partial t}\right)\widehat{n}\cdot \nabla p+\frac{4\mu }{3{\rho}_0^2}\widehat{n}\cdot \nabla Q $$
(41)

Substituting Eq. (41) into Eq. (40) yields the acoustic wave Eq. (39) and can be expressed as follows:

$$ {\displaystyle \begin{array}{c}{\iiint}_{\varOmega_f}\frac{1}{\rho_0{c}^2}w\frac{\partial^2p}{\partial {t}^2} dv+{\iiint}_{\varOmega_f}\nabla w\cdot \left(\frac{4\mu }{3{\rho}_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right) dv+{\iiint}_{\varOmega_f}\nabla w\cdot \left(\frac{1}{\rho_0}\nabla p\right) dv+\\ {}+{\iint}_{\varGamma_f}w\hat{n}\cdot \frac{\partial^2{u}_f}{\partial {t}^2} ds={\iiint}_{\varOmega_f}w\frac{1}{\rho_0}\frac{\partial Q}{\partial t} dv+{\iiint}_{\varOmega_f}\nabla w\cdot \left(\frac{4\mu }{3{\rho}_0^2}\nabla Q\right) dv\end{array}} $$
(42)

The finite element approximating shape functions for the spatial variation of the pressure and displacement components are given by:

$$ P={\left\{N\right\}}^T\left\{p\right\} $$
(43)
$$ u={\left\{{N}^{\hbox{'}}\right\}}^T\left\{u\right\} $$
(44)

with {N} is the element shape function for pressure, {N'} is the element shape function for displacement and {p}, {u} are the nodal pressure and displacements vectors.

Using such formulation into the Eq. (42), the finite element statement of the wave (Eq. (39)) is expressed in matrix notation as following:

$$ \left[{M}_f\right]\left\{\overset{..}{p}\right\}+\left[{C}_f\right]\left\{\overset{.}{p}\right\}+\left[{K}_f\right]\left\{p\right\}+{\overline{\rho}}_0{\left[R\right]}^T\left\{\overset{..}{u}\right\}=\left\{{f}_f\right\} $$
(45)

Where

[Mf]:

\( {\overline{\rho}}_0{\iiint}_{\varOmega_f}\frac{1}{\rho_0{c}^2}\left\{N\right\}{\left\{N\right\}}^T dv \) is the acoustic fluid mass matrix

[Cf]:

\( {\overline{\rho}}_0{\iiint}_{\varOmega_f}\frac{4\mu }{\rho_0^2{c}^2}{\left[\nabla N\right]}^T\left[\nabla N\right] dv \) is the acoustic fluid damping matrix

[Kf]:

\( {\overline{\rho}}_0{\iiint}_{\varOmega_f}\frac{1}{\rho_0}{\left[\nabla N\right]}^T\left[\nabla N\right] dv \) is the fluid stiffness matrix

[R]T:

\( {\iint}_{\varOmega_f}\left\{N\right\}{\left\{n\right\}}^T{\left\{{N}^{\hbox{'}}\right\}}^T ds \) is the acoustic fluid boundary matrix

[ff]:

\( {\overline{\rho}}_0{\iiint}_{\varOmega_f}\frac{1}{\rho_0{c}^2}\left\{N\right\}{\left\{N\right\}}^T dv\left\{\dot{q}\right\}+{\overline{\rho}}_0{\iiint}_{\varOmega_f}\frac{4\mu }{\rho_0^2{c}^2}{\left[\nabla N\right]}^T\left[\nabla N\right] dv\left\{q\right\} \) is the acoustic fluid load with \( {\overline{\rho}}_0 \) the acoustic fluid mass density constant.

So the acoustic and structural matrices are coupled for fluid-structure interaction

$$ \left[\begin{array}{cc}{M}_s& 0\\ {}{\overline{\rho}}_0{R}^T& {M}_f\end{array}\right]\left\{\begin{array}{c}\overset{..}{u}\\ {}\overset{..}{p}\end{array}\right\}+\left[\begin{array}{cc}{C}_s& 0\\ {}0& {C}_f\end{array}\right]\left\{\begin{array}{c}\overset{.}{u}\\ {}\overset{.}{p}\end{array}\right\}+\left[\begin{array}{cc}{K}_s& -R\\ {}0& {K}_f\end{array}\right]\left\{\begin{array}{c}u\\ {}p\end{array}\right\}=\left\{\begin{array}{c}{f}_s\\ {}{f}_f\end{array}\right\} $$
(46)

where the subscript “s” and “f” represent the structural and fluid parts respectively.

3.4 Numerical Simulation Results

In order to validate the numerical models two different materials were considered of which the mechanical properties are well known [37]: a thermoplastic polymer (Polypropylene) and a unidirectional carbon-epoxy composite. For the polypropylene isotropic mechanical properties were considered. Regarding the composite, the matrix was considered to be isotropic, while the anisotropic material parameters for the fiber been applied [22]. The material properties are summarized in Tables 1 and 2.

Table 1 Material Properties for composite T300/914
Full size table
Table 2 Material Properties for Polypropylene
Full size table

A RVE with a cross section of 0.014–0.009 mm was built. For the epoxy-carbon composite, the carbon fibers were included. In particular a RVE with a hexagonal distribution of fiber was modelled. The fibers presented a diameter of 7.0 μm and a fiber volume fraction of 0.6 (Fig. 4). Symmetric and periodic boundary conditions in y and x coordinates, according to the local reference system of Fig. 4 must be applied in order to guaranty the periodicity of the microstructural strain tensor (\( {\varepsilon}_{ij}^{\ast } \)).

Fig. 4
figure 4

RVE model of unidirectional carbon-epoxy composite

Full size image

Regarding the static-structural analysis, a unidirectional tension loading was applied in z direction to estimate the out of plane (E3) elastic modulus. The mesh was fine enough to accurately represent the geometry of the reinforced fibers. Under the prescribed loading condition, the maximum tensile stress on the RVE and the relative strain status in the z direction were used to estimate the out of plane elastic modulus according to Hooke’s law (Eq. (35)). The simulation result in terms of tensional stress (σ3) and strain (ε3) for both materials are shown in Figs. 5 and 6.

Fig. 5
figure 5

Simulation results in terms of equivalent tensile stress σ3 and equivalent strain ε3 for the carbon-epoxy RVE

Full size image
Fig. 6
figure 6

Simulation results in terms of equivalent tensile stress σ3 and equivalent strain ε3 for the Polypropylene RVE

Full size image

Regarding the coupled acoustic-structural analysis the implemented FE model is shown in Fig. 7. The fluid volume into the tube was modelled as an acoustic body. Radiation boundary conditions were used at the input plane of the fluid domain, therefore acoustic waves normal to the boundary will be absorbed and not reflected into the acoustic domain. The acoustic load was applied in the same input plane and modelled using an acoustic mass source. On the other end of the acoustic body, the previous described RVEs were placed as samples and fluid-structure boundary conditions were applied. In order to avoid any translation of the sample a fixed support boundary conditions was applied on the external surface of the test sample. Finally, a coupled acoustic-structural harmonic analysis was set up to measure the pressure profile in two different planes into the fluid volume, according to constrains of Eq. (1) and Eq. (2).

Fig. 7
figure 7

Couple acoustic-structural FE model

Full size image

Starting from the acquired pressure time histories during the simulation, the reflection coefficient of the RVE can be extracted following the transfer function method, and then the acoustic impedance of the RVE can be estimated (Figs. 8 and 9) and the proposed method was applied in the post-processing of data to estimate the elastic modulus.

Fig. 8
figure 8

Acoustic-structural simulation results for carbon-epoxy RVE in terms of fft of pressure time histories at two different acquisition points, Absorption and Reflection Coefficients and Acoustic impedance

Full size image
Fig. 9
figure 9

Acoustic-structural simulation results for polypropylene RVE in terms of fft of pressure time histories at two different acquisition points, Absorption and Reflection Coefficients and Acoustic impedance

Full size image

The results in terms of elastic modulus estimated by the proposed method and the tensile testing simulation compared with the elastic properties declared by the materials datasheet are reported in Table 3.

Table 3 Simulation results for the carbon-epoxy and polypropylene RVE, compared with the elastic properties declared by the data sheet
Full size table

The results obtained from the numerical simulation both static-structural and acoustic-structural analysis are in good agreement with the elastic properties declared by the data-sheet of the relative materials. The proposed approach produced good results for isotropic (polypropylene) and orthotropic (epoxy-carbon composite) materials. The results from the acoustic-structural analysis match the results obtained with the static-structural model which simulate the standard tensile tests, commonly used for the elastic characterization of materials. Therefore, it was shown that using a standard impedance tube test rig both material acoustic characterization and elastic property evaluation of materials can be achieved.

4 Experimental Measurements Campaign

Here the results of an experimental measurement campaign are reported.

Different materials were tested and the proposed method was used to estimate the elastic modulus. In particular, two types of materials (foam and composite) were used to demonstrate the ability of the method to evaluate the Young’s Modulus of soft porous isotropic materials and out of plane Young’s Modulus of hard complex anisotropic materials.

The estimated elastic modulus were compared with those measured by standard testing methods. Regarding the composite samples the ultrasound time of flight (according to ASTM E 494 [38]) was used while for the foam sample, where such method produces inaccurate results because of the scattering of the longitudinal wave through the holes, a tensile test was used (according to ASTM D 3574 [39]).

4.1 Experimental Setup

Four samples was tested in the present measurement campaign.

The tested materials include two acoustic soft foams with different densities (Foam 1 (ρ = 20 kg/m3) – Foam 2 (ρ = 25 kg/m3)), and two composites. The composite materials consist of epoxy resin (HexFlow RTM 6) reinforced by carbon fibre (HexForce Carbon Fabrics G1157) (sample 1) and Vynilester reinforced by glass fibre (sample 2).

The test rig used in the present work was the standard two microphones impedance tube described in the ASTM standard [14]. Such equipment is normally used for the measurement of the acoustic properties of porous and fibre materials.

When testing elastic properties of thermoplastic, polymer or composite materials within the test rig, there is potential for critical errors to occur. These errors are generally related to the phase and amplitude mismatch between the microphone measurements and attenuation of the tube walls associated to viscous thermal losses. In this case, these materials behave as acoustically rigid and the reflection coefficient is normally greater than 0.98–0.99, thus limiting the accuracy of the measurement. As the elastic properties are a function of the longitudinal propagation speed in the material which is strictly related to the reflection coefficient, the measurement of such coefficient must as accurate as possible.

A calibration procedure was applied to take into account the phase mismatch and the viscous thermal losses at the walls of the tube in order to increase the accuracy of the experimental measurements.

Cylindrical samples with a diameter of 50.8 mm were tested in an aluminium impedance tube with an internal diameter of 50.8 mm, so the experimental results are confined to the working frequency range (250–3800 Hz). A loudspeaker at one end of the tube was used to generate a broadband signal; in particular the samples were excited by a 3 s sine-sweep in the considered working frequency range. In addition, the amplitude-phase mismatch of the microphones measurements and the thermal-viscous losses were compensated according to the following proposed calibration procedure. Figure 10 shows the calibration effect on the measured reflection coefficient. The effect of using no correction data is to underestimate the reflection coefficient which lead to a under estimation of the acoustic impedance. At low frequencies the uncalibrated reflection coefficient estimation error was over 3%, a critical error for high reflective materials such as composites and polymeric materials.

Fig. 10
figure 10

Effect of the calibration effect on the measured Reflection coefficient of composite material

Full size image

The sample placing into the tube were done carefully in order to fit the sample itself with the internal wall of the tube. Any relative moving of the sample with the tube and any leakage must be avoid. The downstream section of the tube consisted of a hard-back termination attached to the downstream surface of the sample.

The two upstream measurement position were spaced of 30 mm, the second microphone position were separated by a 140 mm from the incident surface of the sample and the first microphone position were placed 330 mm from the speaker in order to guarantee the completely development of the standing plane wave into the tube.

Regarding the tensile testing carried out on the foam samples (Foam 1 and Foam 2) a dog bone shaped sample was used (refer to Fig. 11) with the material elastic modulus determined from measured stress-strain curves.

Fig. 11
figure 11

Dog bone shape samples for tensile testing for the Foam 1 sample and Foam 2 sample

Full size image

For the ultrasound time of flight measurements an ultrasonic pulse was sent which propagates through the sample with a certain longitudinal sound speed. Measuring the time between the front wall and back wall echo of the samples, the time of flight was estimated. Moreover, since the sample thickness is known, the longitudinal sound wave into the material can be calculated and used to evaluate the elastic modulus of the tested sample.

4.2 Experimental Results

Since the considered test rig is designed to guarantee a plane standing wave, the hypothesis of normal incident at the fluid-sample interface is fulfilled and longitudinal sound wave will propagated into the sample, thus estimation of the out of plane elastic properties is possible.

The calibrated sound pressures measured at two microphone positions were used to estimate the reflection coefficient and the acoustic impedance of the samples using the Transfer Function Method.

To reduce the effect of errors related to mounting the samples, material variability and the uncertainty of experimental measurements, five different measurements were performed for each sample. The results were consistent in terms of the reflection coefficient and elastic modulus. Moreover, since the tests were performed under plane wave conditions and only longitudinal waves propagate through the sample, the out of plane elastic properties are a function of the fibre volume fraction and not a function of fibre orientation.

The measured reflection coefficient and the relative acoustic impedance for the tested materials are shown in Figs. 12, 13, 14 and 15.

Fig. 12
figure 12

Measured reflection coefficient and acoustic impedance for Sample 1 (Epoxy-Carbon fiber composite)

Full size image
Fig. 13
figure 13

Measured reflection coefficient and acoustic impedance for Sample 2 (Vynilester-Glass fiber composite)

Full size image
Fig. 14
figure 14

Measured reflection coefficient and acoustic impedance soft acoustic Foam 1 (ρ = 20 kg/m3)

Full size image
Fig. 15
figure 15

Measured reflection coefficient and acoustic impedance soft acoustic Foam 2 (ρ = 25 kg/m3)

Full size image

Moreover in Fig. 16 are shown the stress-strain curves measured during the tensile testing on the foam sample (Foam 1 and Foam 2).

Fig. 16
figure 16

Stress-Strain curves for the tested acoustic soft foam

Full size image

Finally, the initial pulse and the back wall echo measured for both composite samples (Sample 1 and Sample 2) are shown in Fig. 17.

Fig. 17
figure 17

Initial pulse and the back wall echo for the Sample 1 (Epoxy-carbon fibre and Sample 2 (vynilester-glass sample)

Full size image

In Table 4 the elastic properties of the tested samples were estimated with the proposed method and compared with the relative elastic properties estimated with the commonly used methods.

Table 4 Elastic properties of tested materials evaluated using the proposed method and compared with the results obtained through the commonly used methods
Full size table

Considering the acoustic soft foam the results in terms of elastic modulus estimated with the proposed method match with the elastic properties estimated by tensile tests with errors of under 8% for the Foam 1 and under 2% for Foam 2. It is clear that, the random distribution of porosity within the foam, allows for it to be considered isotropic.

When composite materials were tested, the results are always in good agreement with the results obtained with the time of flight approach both in terms of longitudinal sound speed and elastic modulus, but in this case the measurement uncertainties lead to error of 12% and 6% on the elastic modulus and 4% and 3% on the longitudinal speed respectively for the epoxy-carbon composite and vynilester-glass composite. In this case the mismatch of the results are associated with thermal-viscous losses through the impedance tube and to the amplitude-phase mismatch of the microphones. Reduced reflection coefficient accuracy occurs with highly reflective materials leading to a reduction in the accuracy of the estimated elastic properties (as explained in the previous section).

5 Conclusions

This paper proposes a method for the estimation of the elastic properties of porous, polymeric and composite materials using sound wave by measuring the reflection coefficient and the acoustic impedance in an impedance tube. The procedure described is based on the Transfer Function Method that derives information about the materials reflection coefficient. Under such hypothesis, the reflection coefficient is a function of the longitudinal speed of sound, of the acoustic impedance and as a consequence and the elastic properties of the material. A procedure was proposed to reduce the errors and the uncertainties on the measurements for high acoustic reflective materials such as composite or polymeric materials. The method was first tested in-silico using representative volume element and then experimentally verified.

The proposed approach would allow for the estimation of the elastic modulus for materials where ultrasound time of flight cannot be applied, for example porous materials. An error of only 2% is achieved when comparing the elastic modulus results with those obtained by destructive test (tensile testing).

For the first time, it was shown that using sound waves the elastic properties of composite materials can be evaluated where no standard test is available. In this case, the measured elastic modulus represents the out of plane elastic modulus (E3) and measurement errors were around 3% for the longitudinal speed of sound and 10% for the elastic modulus with respect to ultrasonic testing.

The proposed methods is a reliable, effective and quick non-destructive method for the Elastic Modulus evaluation providing clear advantage over traditional techniques such as standard destructive methods or ultrasound time of flight measurements.