Abstract:
A method to measure the complex frequency-dependent dilatational and shear wavenumbers of a material under a static compressional force. The material is first vibrated in a vertical and horizontal directions while obtaining transfer functions in each direction. The two transfer functions are combined with a theoretical model to estimate a dilatational wavenumber and a shear wavenumber. The wavenumbers can be utilized to give the complex dilatational wavespeed, complex shear wavespeed, complex Lamé constants, complex Young&#39;s modulus, complex shear modulus, and complex Poisson&#39;s ratio.

Description:
STATEMENT OF GOVERNMENT INTEREST 
   The invention described herein may be manufactured and used by or for the Government of the United States of America for governmental purposes without the payment of any royalties thereon or therefore. 

   BACKGROUND OF THE INVENTION 
   (1) Field of the Invention 
   The present invention relates to a method to measure (or estimate) the complex frequency-dependent dilatational and shear wavenumbers of a single slab of material subjected to large static compressional forces. More particularly, this invention provides a method to determine complex dilatational wavespeed, complex shear wavespeed, complex Lamé constants, complex Young&#39;s modulus, complex shear modulus, and complex Poisson&#39;s ratio. 
   (2) Description of the Prior Art 
   Measuring the mechanical properties of slab-shaped materials are important because these parameters significantly contribute to the static and dynamic response of structures built with such materials. One characteristic that most elastomeric solids possess is that, when they are subjected to large static forces (or pressure), their rigidity changes. Materials that have one set of mechanical properties at a pressure of one atmosphere can have very different properties when subjected to increased pressure. The ability to determine the pressure dependence of material properties is extremely important for modeling the behavior of systems comprised of these materials. 
   Resonant techniques have been used to identify and measure longitudinal and shear properties for many years. These methods are based on comparing measured eigenvalues to modeled eigenvalues and calculating the resulting material properties. These methods do not account for static pressure or large compressive forces. Additionally, they typically require long, slender materials to perform the measurement process. Comparison of analytical models to measured frequency response functions are also used to estimate stiffness and loss parameters of a structure. When the analytical model agrees with one or more frequency response functions, the parameters used to calculate the analytical model are considered accurate. If the analytical model is formulated using a numerical method, a comparison of the model to the data can be difficult due to dispersion properties of the materials. These methods do not take into account large compressive forces. 
   In the prior art, some efforts have been made to measure material properties under large pressures. These methods consist of placing materials in pressurized settings, insonifying them, and then measuring their response. These methods are difficult because they have to be conducted under great atmospheric pressure that can adversely affect the instrumentation. Safety issues can also arise in connection with laboratory testing at extreme pressures. Finally, a mass loaded long thin rod has been studied with respect to the bar wavespeed and corresponding Young&#39;s modulus; however, this work does not investigate shear motion. 
   Accordingly, there is a need for a method of measuring mechanical properties of slab-shaped materials placed under pressure. 
   SUMMARY OF THE INVENTION 
   Accordingly, in this invention, a method to measure the complex frequency-dependent dilatational and shear wavenumbers of a material under a static compressional force is provided. The material is first vibrated in both vertical and horizontal directions while obtaining transfer functions in each direction. The two transfer functions are combined with a theoretical model to estimate a dilatational wavenumber and a shear wavenumber. The wavenumbers can be combined to give the complex dilatational wavespeed, complex shear wavespeed, complex Lamé constants, complex Young&#39;s modulus, complex shear modulus, and complex Poisson&#39;s ratio. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     A more complete understanding of the invention and many of the attendant advantages thereto will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in conjunction with the accompanying drawings wherein: 
       FIG. 1  shows apparatus for measurement of transfer functions in a vertical direction according to the current invention; 
       FIG. 2  shows apparatus for measurement of transfer functions in a horizontal direction according to the current invention; 
       FIG. 3  is a diagram of the coordinate system of used with a test specimen in the model; 
       FIG. 4A  is a plot of the transfer function magnitude versus input frequency for the vertical direction test; 
       FIG. 4B  is a plot of the transfer function phase angle versus input frequency for the vertical direction test; 
       FIG. 5A  is a plot of the transfer function magnitude versus input frequency for the horizontal direction test; 
       FIG. 5B  is a plot of the transfer function phase angle versus input frequency for the horizontal direction test; 
       FIG. 6  is a contour plot of the absolute value of the dilatational wavenumber on an real-imaginary coordinate system of the dilatational wavenumbers at 2000 Hz; 
       FIG. 7  is a contour plot of the absolute value of the dilatational wavenumber on an real-imaginary coordinate system of the dilatational wavenumbers at 5000 Hz; 
       FIG. 8A  is a plot of the real dilatational wavenumber versus frequency; 
       FIG. 8B  is a plot of the imaginary dilatational wavenumber versus frequency; and 
       FIG. 9A  is a plot of the real shear wavenumber versus frequency; 
       FIG. 9B  is a plot of the imaginary shear wavenumber versus frequency; 
       FIG. 10  is a plot of the real and imaginary Young&#39;s modulus versus frequency; 
       FIG. 11  is a plot of the real and imaginary shear modulus versus frequency; and 
       FIG. 12  is a plot of the Poisson&#39;s ratio versus frequency. 
   

   DESCRIPTION OF THE PREFERRED EMBODIMENT 
   The test procedure consists of vibrating a mass-loaded, slab-shaped test specimen  10  with a shaker  12  in two different directions, vertical  14 A and horizontal  14 B, as shown in  FIGS. 1 and 2 , respectively. It is noted that the load mass  16  attached to the top of the test specimen  10  must be sufficiently stiffer than the specimen  10  that it can be modeled as lumped parameter expression rather than a continuous media system. A typical example would be a steel load mass  16  attached above a rubber-like material test specimen  10 . This example results in a ratio between the two stiffnesses of greater than 100. Lower ratios result in less accurate estimations. Vibrating the combined specimen  10  and load mass  16  causes different waveforms to propagate in the specimen  10 . The inverse method developed here allows for the data from the experiments to be manipulated so that the complex dilatational and shear wavenumbers can be measured for the specimen  10 . This test is usually done at multiple frequencies (swept sine) so any frequency dependencies can be identified and measured. Input vibration data is collected from the shaker  12 . A sensor  18  is mounted on load mass  16  and another sensor  20  is mounted on shaker  12  for collecting transfer function data. In  FIG. 1 , the test is set up for monitoring the vertical transfer function.  FIG. 2  shows the test as set up for monitoring the horizontal transfer function. Sensors  18  and  20  should be oriented properly to capture the motion being measured. Other test configurations using directions other than vertical and horizontal are possible; however, the test setups shown are preferred for ease of set up and calculation. These sensors  18  can be either accelerometers that record accelerations, or laser velocimeters that record velocities. In the swept sine mode, transfer functions of acceleration divided by acceleration or velocity divided by velocity are both equal to displacement divided by displacement. The time domain data collected from the sensors  18  and  20  are Fourier transformed into the frequency domain and then recorded as complex transfer functions, typically using a spectrum analyzer  22 . 
   The motion of the test specimen shown in  FIGS. 1 and 2  is governed by the equation: 
                   μ   ⁢           ⁢       ∇   2     ⁢   u       +       (     λ   +   μ     )     ⁢       ∇   ∇     ·   u         =     ρ   ⁢         ∂   2     ⁢   u       ∂     t   2             ,           (   1   )             
 
where λ and μ are the complex Lamé constants (N/m 2 );
         ρ is the density (kg/m 3 );   t is time (s);   · denotes a vector dot product; and   u is the Cartesian coordinate displacement vector of the material.       

   The coordinate system of the test configuration is shown in FIG.  3 . Note that using this orientation results in b=0 and a having a value less than zero. The thickness of the specimen, h, is a positive value. Equation (1) is manipulated by writing the displacement vector u as 
               u   =     {             u   x     ⁡     (     x   ,   y   ,   z   ,   t     )                   u   y     ⁡     (     x   ,   y   ,   z   ,   t     )                   u   z     ⁡     (     x   ,   y   ,   z   ,   t     )             }       ,           (   2   )             
 
where x is the location along the plate (m), y is the location into the plate (m), and z is the location normal to the plate (m), as shown in FIG.  3 . The symbol ∇ is the gradient vector differential operator written in three-dimensional Cartesian coordinates as 
               ∇     =         ∂     ∂   x       ⁢     i   x       +       ∂     ∂   y       ⁢     i   y       +       ∂     ∂   z       ⁢     i   z             ,           (   3   )             
 
with i x  denoting the unit vector in the x-direction, i y  denoting the unit vector in the y-direction, and i z  denoting the unit vector in the z-direction; ∇ 2  is the three-dimensional Laplace operator operating on vector u as
 
∇ 2 =∇ 2   u   x   i   x =∇ 2   u   y   i   y =∇ 2   u   2   i   z ,  (4) 
 
and operating on scalar u as 
                   ∇   2     ⁢     u     x   ,   y   ,   z         =       ∇     ·     ∇     u     x   ,   y   ,   z             =           ∂   2     ⁢     u     x   ,   y   ,   z           ∂     x   2         +         ∂   2     ⁢     u     x   ,   y   ,   z           ∂     y   2         +         ∂   2     ⁢     u     x   ,   y   ,   z           ∂     z   2               ;           (   5   )             
 
and the term ∇·u is called the divergence and is equal to 
               ∇     ·   u       =         ∂     u   x         ∂   x       +       ∂     u   y         ∂   y       +         ∂     u   z         ∂   z       .               (   6   )             
 
The displacement vector u is written as
 
u=∇φ+∇×{overscore (ψ)},  (7) 
 
where φ is a dilatational scalar potential, ×denotes a vector cross product, and {overscore (ψ)} is an equivoluminal vector potential expressed as 
               ψ   -&gt;     =       {             ψ   x     ⁡     (     x   ,   y   ,   z   ,   t     )                   ψ   y     ⁡     (     x   ,   y   ,   z   ,   t     )                   ψ   z     ⁡     (     x   ,   y   ,   z   ,   t     )             }     .             (   8   )             
 
The problem is formulated as a two-dimensional system, thus y=0, u y (x,y,z,t)=0, and ∂(·)/∂y≡0. Expanding equation (7) and breaking the displacement vector into its individual nonzero terms yields 
                 u   x     ⁡     (     x   ,   z   ,   t     )       =         ∂     ϕ   ⁡     (     x   ,   z   ,   t     )           ∂   x       -       ∂       ψ   y     ⁡     (     x   ,   z   ,   t     )           ∂   z                 (   9   )             
 
and 
                 u   z     ⁡     (     x   ,   z   ,   t     )       =         ∂     ϕ   ⁡     (     x   ,   z   ,   t     )           ∂   z       +         ∂       ψ   y     ⁡     (     x   ,   z   ,   t     )           ∂   x       .               (   10   )             
 
   Equations (9) and (10) are next inserted into equation (1), which results in 
                 c   d   2     ⁢       ∇   2     ⁢     ϕ   ⁡     (     x   ,   z   ,   t     )           =         ∂   2     ⁢     ϕ   ⁡     (     x   ,   z   ,   t     )           ∂     t   2                 (   11   )             
 
and 
                 c   s   2     ⁢       ∇   2     ⁢       ψ   y     ⁡     (     x   ,   z   ,   t     )           =         ∂   2     ⁢       ψ   y     ⁡     (     x   ,   z   ,   t     )           ∂     t   2                 (   12   )             
 
where equation (11) corresponds to the dilatational component and equation (12) corresponds to the shear component of the displacement field. Correspondingly, the constants c d  and c s  are the complex dilatational and shear wave speeds, respectively, and are determined by 
               c   d     =         λ   +     2   ⁢           ⁢   μ       ρ               (   13   )             
 
and 
               c   s     =         μ   ρ       .             (   14   )             
 
The relationship of the Lamé constants to the Young&#39;s and shear moduli is shown as 
             λ   =       E   ⁢           ⁢   υ         (     1   +   υ     )     ⁢     (     1   -     2   ⁢           ⁢   υ       )                 (   15   )             
 
and 
               μ   =     G   =     E     2   ⁢     (     1   +   υ     )             ,           (   16   )             
 
where E is the complex Young&#39;s modulus (N/m 2 ), G is the complex shear modulus (N/m 2 ), and ν is the Poisson&#39;s ratio of the material (dimensionless).
 
   The conditions of infinite length and steady-state response are now imposed, allowing the scalar and vector potential to be written as
 
φ( x,z,t )=Φ( z )exp( ikx )exp( iωt )  (17) 
 
and
 
ψ y ( x,z,t )=ψ( z )exp( ikx )exp( iωt )  (18) 
 
where i is the square root of −l, ω is frequency (rad/s), and k is wavenumber with respect to the x axis (rad/m). Inserting equation (17) into equation (11) yields 
                       ⅆ   2     ⁢     Φ   ⁡     (   z   )           ⅆ     z   2         +       α   2     ⁢     Φ   ⁡     (   z   )           =   0     ,           (   19   )             
 
where 
               α   =         k   d   2     -     k   2           ,           (   20   )             
 
and 
               k   d     =       ω     c   d       .             (   21   )             
 
Inserting equation (18) into equation (12) yields. 
                       ⅆ   2     ⁢     Ψ   ⁡     (   z   )           ⅆ     z   2         +       β   2     ⁢     Ψ   ⁡     (   z   )           =   0     ,           (   22   )             
 
where 
               β   =         k   s   2     -     k   2           ,           (   23   )             
 
and 
               k   s     =       ω     c   s       .             (   24   )             
 
   The solution to equation (19) is
 
Φ( z )= A ( k,ω )exp( iβz )+ D ( k ,ω)exp(− iαz )+ B ( k,ω )exp(− iαz ),  (25) 
 
and the solution to equation (22) is
 
ψ( z )= C ( k ,ω)exp( iβz )+ D ( k ,ω)exp(− iβz ),  (26) 
 
where A(k,ω), B(k,ω), C(k,ω), and D(k,ω) are wave response coefficients that are determined below. The displacements can now be written as functions of the unknown constants using the expressions in equations (9) and (10). They are 
                       u   z     ⁡     (     x   ,   z   ,   t     )       =       ⁢         U   z     ⁡     (     k   ,   z   ,   ω     )       ⁢     exp   ⁡     (   ikx   )       ⁢     exp   ⁡     (     i   ⁢           ⁢   ω   ⁢           ⁢   t     )                     =       ⁢     {       i   ⁢           ⁢     α   ⁡     [         A   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (     i   ⁢           ⁢   α   ⁢           ⁢   z     )         -       B   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (       -   i     ⁢           ⁢   α   ⁢           ⁢   z     )           ]         +                       ⁢     ik   [         C   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (     i   ⁢           ⁢   β   ⁢           ⁢   z     )         +                             ⁢     D   ⁢     (     k   ,   ω     )     ⁢     exp   ⁡     (       -   i     ⁢           ⁢   β   ⁢           ⁢   z     )         ]     }     ⁢     exp   ⁡     (   ikx   )       ⁢     exp   ⁡     (     i   ⁢           ⁢   ω   ⁢           ⁢   t     )         ,                 (   27   )             
 
and 
                       u   x     ⁡     (     x   ,   z   ,   t     )       =       ⁢         U   x     ⁡     (     k   ,   z   ,   ω     )       ⁢     exp   ⁡     (   ikx   )       ⁢     exp   ⁡     (     i   ⁢           ⁢   ω   ⁢           ⁢   t     )                     =       ⁢     {       i   ⁢           ⁢     k   ⁡     [         A   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (     i   ⁢           ⁢   α   ⁢           ⁢   z     )         +       B   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (       -   i     ⁢           ⁢   α   ⁢           ⁢   z     )           ]         -                       ⁢     i   ⁢           ⁢     β   [         C   ⁡     (     k   ,   ω     )       ⁢     exp   ⁡     (     i   ⁢           ⁢   β   ⁢           ⁢   z     )         -                             ⁢     D   ⁢     (     k   ,   ω     )     ⁢     exp   ⁡     (       -   i     ⁢           ⁢   β   ⁢           ⁢   z     )         ]     }     ⁢     exp   ⁡     (   ikx   )       ⁢       exp   ⁡     (     i   ⁢           ⁢   ω   ⁢           ⁢   t     )       .                   (   28   )             
 
Specific boundary conditions are now needed to individually solve for the case of vertical and horizontal motion. These are formulated separately below.
 
   For the case of vertical motion, the base at z=a is vibrated vertically using a shaker, as shown in FIG.  1 . Four boundary conditions are necessary to formulate this problem. Because the mass is attached to the material, the tangential (horizontal) motion at the top of the plate (z=b) is zero and this equation is written as
 
 u   x ( x,b,t )=0.  (29) 
 
The normal stress at the top of the specimen is equal to and opposite the load created by the mass in the z direction. This expression is 
                   τ   zz     ⁡     (     x   ,   b   ,   t     )       =           (     λ   +     2   ⁢   μ       )     ⁢       ∂       u   z     ⁡     (     x   ,   b   ,   t     )           ∂   z         +     λ   ⁢       ∂       u   x     ⁡     (     x   ,   b   ,   t     )           ∂   x           =       -   M     ⁢         ∂   2     ⁢       u   z     ⁡     (     x   ,   b   ,   t     )           ∂     t   2               ,           (   30   )             
 
where M is mass per unit area (kg/m 2 ) of the attached mass. The tangential motion at the bottom of the plate (z=a) is zero and this equation is written as
 
 u   x ( x,a,t )=0.  (31) 
 
and the normal motion at the bottom of the plate is prescribed as a system input. This expression is
 
 u   z ( x,a,t )= U   0  exp(iωt).  (32) 
 
Assembling equations (1)-(32) and letting b=0 yields the four-by-four system of linear equations that model the system. They are
 
 Ax=b ,  (33) 
 
where the entries of equation (33) are
 
 A   11   =ik,   (34) 
 
 A   12   =A   11 ,  (35) 
 
 A   13   =−iβ,   (36) 
 
  A   14   =A   13 ,  (37)
 
 A   21 =−α 2 λ− 2 α 2   −λk   2   −iMω   2 α,  (38) 
 
 A   22 =−α 2 λ− 2 α 2   −λk   2   −iMω   2 α,  (39) 
 
 A   23 =−2 kβμ−iMω   2   k,   (40) 
 
 A   24 =−2 kβμ−iMω   2   k,   (41) 
 
 A   31   =A   11  exp( iαa )  (42) 
 
 A   32   =A   11  exp(− iαa )  (43) 
 
 A   33   =A   13  exp( iβa ),  (44) 
 
 A   34   =A   13  exp(− iβa ),  (45) 
 
 A   41   =i α exp( iαa ), (46) 
 
 A   42   =iα  exp(− iβa ), (47) 
 
 A   43   =ik  exp( iβa ),  (48) 
 
  A   44   =ik  exp(− iβa ),  (49)
 
 x   11   =A ( k ,ω),  (50) 
 
 x   21   =B ( k ,ω),  (51) 
 
 x   31   =C ( k ,ω),  (52) 
 
 x   41   =D ( k ,ω),  (53) 
 
 b   11 =0,  (54) 
 
 b   21 =0,  (55) 
 
 b   31 =0,  (56) 
 
and
 
 b   41   =U   0 ,  (57) 
 
Using equations (34)-(57) the solution to the constants A(k,ω), B(k,ω), C(k,ω), and D(k,ω) can be calculated at each specific wavenumber and frequency using
 
x=A −1 b.  (58) 
 
Noting that for vertical motion, k=0, and using the coefficients from equation (58), the transfer function between the vertical base displacement and the vertical mass displacement can be written as 
                   T   1     ⁡     (   ω   )       =       1       R   1     ⁡     (   ω   )         =           U   z     ⁡     (     0   ,   b   ,   ω     )         U   0       =     1       cos   ⁡     (       k   d     ⁢   h     )       -       (     M   ρ     )     ⁢     k   d     ⁢     sin   ⁡     (       k   d     ⁢   h     )                   ,           (   59   )             
 
where T 1 (ω) or R 1 (ω) correspond to the data from the vertical motion experiment.
 
   The next step is to solve the inverse problem for vertical motion. This involves using the experimental data and equation (59) to estimate the dilatational wavenumber. Equation (59) can be rewritten as 
               f   ⁡     (     k   d     )       =     0   =       cos   ⁡     (       k   d     ⁢   h     )       -       (     M   ρ     )     ⁢     k   d     ⁢     sin   ⁡     (       k   d     ⁢   h     )         -     R   1                 (   60   )             
 
where the problem now becomes finding the zeros of the right-hand side of equation (60), or, in the presence of actual data that contains noise, finding the relative minima of the right-hand side of equation (60) and determining which relative minimum corresponds to dilatational wave propagation and which relative minima are extraneous. Because equation (60) has a number of relative minima, zero finding algorithms are not applied to this function, as they typically do not find all of the minima locations and are highly dependent on initial starting locations. The best method to find all of the minima locations is by plotting the absolute value of the right-hand side of equation (60) as a surface with the real part of dilatational wavenumber k d  on one axis and the imaginary part of k d  on the other axis. In order to do this, the user should start at a low frequency where the aliasing minimum has not yet appeared. In the specific example shown herein, this is below 3850 Hz for the dilatational wave and below 1550 Hz for the shear wave. At these lower frequencies, the minimum furthest to the left will correspond to dilatational wave propagation. As the frequency increases, extraneous minima will appear to the left of the minimum that corresponds to dilatational wave propagation, however, the wave propagation minimum will always be close to the previous test frequency wave propagation minimum provided that the frequency increments are relatively small. At a resolution of 0.5 rad/m for the materials in the example herein, this requires a frequency increment of 37.3 Hz for the dilatational measurement and 14.4 Hz for the shear measurement. Different test specimens and top masses require different increments. Additionally, the real part of the wavenumber is a monotonically increasing function with respect to frequency, so at each increase in frequency, the new wavenumber to be estimated has to be greater than the old wavenumber that was previously estimated. This process is further illustrated as related to the discussion concerning FIG.  6  and  FIG. 7  below.
 
   For the case of horizontal motion, the base at z=a is vibrated horizontally using a shaker, as shown in FIG.  2 . Four boundary conditions are necessary to formulate this problem. Because the mass is attached to the material, the shear (tangential) stress at the top of the plate is equal to opposite the load created by the mass in the x direction. This expression is 
                   τ   zx     ⁡     (     x   ,   b   ,   t     )       =       μ   ⁡     [         ∂       u   x     ⁡     (     x   ,   b   ,   t     )           ∂   z       +       ∂       u   z     ⁡     (     x   ,   b   ,   t     )           ∂   x         ]       =       -   M     ⁢         ∂   2     ⁢       u   x     ⁡     (     x   ,   b   ,   t     )           ∂     t   2               ,           (   61   )             
 
where M is mass per unit area (kg/m 2 ) of the attached mass. The normal motion at the top of the plate (z=b) is zero and this equation is written as
 
 u   z ( x,b,t )=0.  (62) 
 
The tangential motion at the bottom of the plate (z=a) is prescribed as a system input and this equation is written as
 
 u   x ( x,a,t )= V   0  exp( iωt ).  (63) 
 
and the normal motion at the bottom of the plate is zero. This expression is
 
 u   z ( x,a,t )=0.  (64) 
 
   Assembling equations (1)-(28) and (62)-(64) and letting b=0 yields the four-by-four system of linear equations that model the system. They are
 
 Ax=b,   (65) 
 
where the entries of equation (61) are
 
 A   11 =−2μ kα−iω   2   Mk,   (66) 
 
 A   12 =−2μ kα−iω   2   Mk,   (67) 
 
 A   13 =μβ 2   −μk   2   +iω   2   Mβ,   (68) 
 
 A   14 =μβ 2   −μk   2   −iω   2   Mβ,   (69) 
 
 A   21   =iα,   (70) 
 
 A   22   =−A   21 ,  (71) 
 
  A   23   =ik,   (72)
 
 A   24   =A   23 ,  (73) 
 
 A   31   =A   23  exp( iαa ),  (74) 
 
 A   32   =A   23  exp(− iαa ),  (75) 
 
 A   33   =−iβ  exp( iβa ),  (76) 
 
 A   34   =−i β exp(− iβa ),  (77) 
 
 A   41   =A   21  exp( iαa ),  (78) 
 
 A   41   =−A   21  exp(− iαa ),  (79) 
 
 A   43   =A   23  exp( iβa ),  (80) 
 
 A   44   =A   23  exp(− iβa ),  (81) 
 
 x   11   =A ( k ,ω), (82) 
 
 x   21   =B ( k ,ω), (83) 
 
  x   31   =C ( k ,ω), (84)
 
 x   41   =D ( k ,ω), (85) 
 
 b   11 =0,  (86) 
 
 b   21 =0,  (87) 
 
 b   31   =V   0 ,  (88) 
 
and
 
 b   41 =0.  (89) 
 
Using equations (67)-(89) the solution to the constants A(k,ω), B(k,ω), C(k,ω), and D(k,ω) can be calculated at each specific wavenumber and frequency using
 
 x=A   −1   b.   (90) 
 
Noting that for horizontal motion, k=0, and using the coefficients from equation (90), the transfer function between the horizontal base displacement and the horizontal mass displacement can be written as 
                 T   2     ⁡     (   ω   )       =       1       R   2     ⁡     (   ω   )         =           U   x     ⁡     (     0   ,   b   ,   ω     )         V   0       =       1       cos   ⁡     (       k   s     ⁢   h     )       -       (     M   ρ     )     ⁢     k   s     ⁢     sin   ⁡     (       k   s     ⁢   h     )             .                 (   91   )             
 
where T 2 (ω) or R 2 (ω) correspond to the data from the horizontal motion experiment.
 
   The next step is to solve the inverse problem for horizontal motion. This involves using the data and equation (91) to estimate the shear wavenumber. Equation (91) can be rewritten as 
               f   ⁡     (     k   s     )       =     0   =       cos   ⁡     (       k   s     ⁢   h     )       -       (     M   ρ     )     ⁢     k   s     ⁢     sin   ⁡     (       k   s     ⁢   h     )         -       R   2     .                 (   92   )             
 
It is noted that this equation is identical, except for the subscripts, to equation (60). The shear wavenumber is estimated using the same procedure that was used to estimate the dilatational wavenumber above.
 
   The material properties can be determined from the wavenumbers. First, the dilatational and shear wavespeeds are determined using 
               c   d     =     ω     k   d               (   93   )             
 
and 
                 c   s     =     ω     k   s         ,           (   94   )             
 
respectively. The Lamé constants are calculated using equations (13) and (14) written as 
             μ   =     ρ   ⁢           ⁢     c   s   2               (   95   )             
 
and 
             λ   =       ρ   ⁢           ⁢     c   d   2       -     2   ⁢   ρ   ⁢           ⁢       c   s   2     .                 (   96   )             
 
Poisson&#39;s ratio is then calculated using 
             υ   =       λ     2   ⁢     (     μ   +   λ     )         .             (   97   )             
 
Young&#39;s modulus can be calculated with 
             E   =       2   ⁢     μ   ⁡     (       2   ⁢   μ     +     3   ⁢   λ       )           2   ⁢     (     μ   +   λ     )                 (   98   )             
 
and the shear modulus can be determined using
 
 G=μ.   (99) 
 
The above measurement method can be simulated by means of a numerical example. Soft rubber-like material properties of the test specimen are used in this simulation. The material has a Young&#39;s modulus E of [(1e8−i2e7)+(5e3∫−i3e2∫)]N/m 2  where f is frequency in Hz, Poisson&#39;s ratio ν is equal to 0.40 (dimensionless), density ρ is equal to 1200 kg/m 3 , and a thickness h of 0.1 m. The top mass is a 0.0254 m (1 inch) steel plate that has a mass per unit area value M of 199 kg/m 2 .  FIG. 4  is a plot of the transfer function of the system for vertical motion and corresponds to equation (59).  FIG. 5  is a plot of the transfer function of the system for horizontal motion and corresponds to equation (91). In  FIGS. 4 and 5 , the top plot is the magnitude in decibels and the bottom plot is the phase angle in degrees.
 
     FIG. 6  is a contour plot of the absolute value of equation (60) expressed in decibels versus real dilatational wavenumber on the x axis and imaginary dilatational wavenumber on the y axis at 2000 Hz. The estimated dilatational wavenumber, read directly from the plot at the location the minimum value appears and marked with an arrow, is 27.89+2.61i rad/m. The actual value of the dilatational wavenumber is 27.99+2.60i rad/m, which is slightly different from the estimated value due to the surface discretization of equation (60).  FIG. 7  is a contour plot of equation (60) at 5000 Hz. At this frequency, an extraneous minimum appears on the left-hand side of the plot. However, because the real part of the wavenumber must be increasing with increasing frequency, the minimum corresponding to dilatational wave propagation is located at the arrow marked spot and is equal to 65.86+5.60i rad/m, as compared to an actual value of 65.77+5.62 i rad/m. Again, the difference between the two values can be attributed to the discretization of the surface. 
     FIG. 8  is plot of actual (solid line) and estimated (x symbols) dilatational wavenumber versus frequency.  FIG. 9  is plot of actual (solid line) and estimated (+symbols) shear wavenumber versus frequency. In  FIGS. 8 and 9 , the top plot is the real part of the wavenumber and the bottom part is the imaginary part of the wavenumber.  FIG. 10  is a plot of actual (solid line) and estimated (real part—× symbols, imaginary part—o symbols) Young&#39;s modulus versus frequency.  FIG. 11  is a plot of actual (solid line) and estimated (real part—× symbols, imaginary part—o symbols) shear modulus versus frequency. In  FIGS. 10 and 11 , the imaginary part of the modulus all have a negative sign but are depicted with positive signs for plotting purposes.  FIG. 12  is a plot of actual (solid line) and estimated (square symbols) of the real part of Poisson&#39;s ratio versus frequency. Because the numerical example is formulated using a Poisson&#39;s ratio that is strictly real, no imaginary component is shown in this plot. Imaginary values of Poisson&#39;s ratio are possible and have been shown to theoretically exist. 
   This invention gives the ability to estimate complex dilatational and shear wavespeeds of a material that is slab-shaped and subjected to compressive forces. It also allows estimation of complex Lamé constants of a material that is slab-shaped and subjected to compressive forces. Complex Young&#39;s and shear moduli of a material that is slab-shaped and subjected to compressive forces can be estimated using this invention. The invention also allows estimation of the complex Poisson&#39;s ratio of a material that is slab-shaped and subjected to compressive forces. The advantage of this patent is that it does not require a testing configuration that has to be pressurized. 
   Obviously many modifications and variations of the present invention may become apparent in light of the above teachings. In light of the above, it is therefore understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described.