Device and method for determining the elasticity of soft-solids

Abstract

The invention comprises a device and method to estimate the elasticity of soft elastic solids from surface wave measurements. The method is non-destructive, reliable and repeatable. The final device is low-cost and portable. It is based in audio-frequency shear wave propagation in elastic soft solids. Within this frequency range, shear wavelength is centimeter sized. Thus, the experimental data is usually collected in the near-field of the source. Therefore, an inversion algorithm taking into account near-field effects was developed for use with the device. Example applications are shown in beef samples, tissue mimicking materials and in vivo skeletal muscle of healthy volunteers.

Claims

1. A method for determining the elasticity of a soft-solid, the method comprising: using a wave source for exciting low-amplitude audible frequency waves in a selected location of a free surface of the soft-solid whose elasticity is to be determined; recording the time-traces of the surface displacement with a plurality of contact vibration sensors that are linearly arranged and equally spaced to each other, wherein the distance between the point at which the waves are being excited and the nearest sensor is the same as the distance between sensors; computing the phase velocity of the surface wave by estimating the phase-shift between sensors and a reference signal sent to the source; and converting the phase velocity to an elasticity value by using an inversion algorithm, wherein the soft-solid is an isotropic solid, and wherein the inversion algorithm takes into account near-field effects and employs: a dispersion curve comprising a range of frequencies, or a single frequency value.

2. A method for determining the elasticity of a soft-solid, the method comprising: using a wave source for exciting low-amplitude audible frequency waves in a selected location of a free surface of the soft-solid whose elasticity is to be determined; recording the time-traces of the surface displacement with a plurality of contact vibration sensors that are linearly arranged and equally spaced to each other, wherein the distance between the point at which the waves are being excited and the nearest sensor is the same as the distance between sensors; computing the phase velocity of the surface wave by estimating the phase-shift between sensors and a reference signal sent to the source; and converting the phase velocity to an elasticity value by using an inversion algorithm, wherein the soft-solid is a transversely isotropic solid, and wherein the inversion algorithm takes into account near-field effects and the propagation direction of surface waves and employs: a dispersion curve comprising a range of frequencies, or a single frequency value.

Description

DESCRIPTION OF THE FIGURES

(1) FIG. 1: shows a linear arrangement of vibration sensors (2) and the head piece (17). This is a non-limiting example of piezoelectric vibration sensor, constituted by a PVDF flexible piezoelectric film, it has a mass on the end with a small extension attached to it. Dimensions in millimeters are shown below to the left. The distance d between sensors as well as the distance d′ between the source and the first sensor are also displayed.

(2) FIG. 2: is a flowchart for obtaining data with the present invention. Dotted lines indicate that pressure sensor, temperature sensor and distance sensor are optional elements.

(3) FIG. 3A-3D: shows two alternative assemblies for arranging the vibration sensors on the medium of interest. 3A-3C corresponds to inferior, superior and lateral views, respectively. Scale bar=12 mm.

(4) FIG. 3D is a strap for attaching the sensors to live muscles (14).

(5) FIG. 4: shows the coordinate system orientation for the transversely isotropic solid. The symmetry axis is oriented along the x.sub.1 direction parallel to the fibers (either, muscle fibers or long molecule chains). The rectangular head piece oriented with an angle ϕ with respect to the x.sub.1 direction is also displayed.

(6) FIG. 5: Displays the directivity diagram of shear waves created by a rectangular head piece in a semi-infinite isotropic elastic solid. The maximum amplitude is attained at 34°. Below this angle, the amplitude decreases rapidly to 0.

(7) FIG. 6: is a schematic representation of bulk shear wave propagation within the solid and mode conversion to SP surface wave. The coordinate system is also displayed.

(8) FIG. 7: shows the surface wave phase velocity measured in an isotropic solid as function of frequency for a position x.sub.1<x.sub.s. Dots, experimental values. Full line: Inversion method #1.

(9) FIG. 8: shows the surface wave phase velocity measured in an isotropic solid as function of frequency for a position x.sub.1≈x.sub.s. Dots: experimental values. Full line: Inversion method #3.

(10) FIG. 9: shows the surface wave phase velocity measured in an isotropic solid as function of frequency for a position x.sub.1>2x.sub.s. Dots: Experimental values. Full line: Inversion method #4.

(11) FIG. 10: is a block diagram representing the inversion algorithm for an isotropic solid.

(12) FIG. 11A-11B: displays an example of experimental data showing the dependence of c.sub.55 with temperature in beef samples.

(13) FIG. 12A-12E shows the relative elastic changes during age maturation process in different beef samples. All curves were normalized with respect to the value of c.sub.55 at the first measurement. The errorbars represent one standard deviation over 10 measurements.

(14) FIG. 13A-13C show a comparison of the performance of the present invention with results obtained by a standard Warner-Bratzler shear force (WBSF) test.

(15) FIG. 14 displays an example of using the present invention to rate different beef cuts according to their elasticity value.

(16) FIG. 15A-15C show examples of using the present invention to monitor the elasticity of skeletal muscle in vivo in dynamic situations.

(17) FIG. 16 illustrates an example of the time trace of a recorded signal.

THEORY OF OPERATION

(18) In linear elasticity theory, the stress and strain within an elastic solid are related by the generalized Hooke's law:
τ.sub.mj=C.sub.mjklϵ.sub.kl  1
where the Einstein's convention over repeated indices apply. In this equation τ.sub.ij is the stress tensor, ϵ.sub.kl is the strain tensor defined as

(19) ${ϵ}_{\mathrm{kl}}=\frac{1}{2}\left(\frac{\partial u_k}{\partial u_l}+\frac{\partial u_l}{\partial u_k}\right),$
u.sub.k is the k=1, 2, 3 component of the displacement field and C.sub.mjkl is the stiffness tensor. It has 81 (3.sup.4) elements representing elastic constants. However, since τ.sub.mj=τ.sub.jm and ϵ.sub.kl=ϵ.sub.lk, this number reduces to 36 independent coefficients. Moreover, since the symmetry of the derivative of the strain energy with respect to the strain tensor, this number reduces further to 21 independent coefficients. Thus, the 21 independent elastic coefficients define the general anisotropic elastic solid. In order to avoid working with a fourth rank tensor, it is usual to represent the independent constants of the stiffness tensor by two indices α and β with values 1 to 6 corresponding to a 6×6 array with the following convention:

(20) (11).Math.1 (22).Math.2 (33).Math.3

(21) (23)=(32).Math.4 (31)=(13).Math.5 (12)=(21).Math.6

(22) Thus, C.sub.mjkl=C.sub.αβ with α related to (mj) and β related to (kl).

(23) The fundamental relation of dynamics applied to this system gives:

(24) $\begin{array}{cc}\rho \frac{{\partial }^2{u}_m}{\partial t^2}=\frac{\partial {\tau }_{mj}}{\partial x_j}& 2\end{array}$
where ρ is the material density. When using the Hooke's law (1) to express the stress tensor, this equation reads:

(25) $\begin{array}{cc}\rho \frac{{\partial }^2{u}_m}{\partial t^2}=C_{mjkl}\frac{{\partial }^2{u}_l}{\partial x_j\partial x_k}& 3\end{array}$
which is a system of three second order differential equation accounting for wave propagation in three dimensional anisotropic bodies. Plane wave solutions to this equation are expressed as:

(26) $\begin{array}{cc}{u}_m={\hat{e}}_{m}F\left(t-\frac{{\hat{n}}_j{x}_j}{V}\right)& 4\end{array}$
where ê.sub.m is a unit vector in the direction of particle displacement and is referred as the wave polarization and {circumflex over (n)}=(n.sub.1, n.sub.2, n.sub.3) is a unit vector in the direction of wave propagation. Inserting equation (4) into equation (3) gives:
ρV.sup.2ê.sub.m=C.sub.mjkln.sub.jn.sub.kê.sub.l;   5
Introducing the second rank tensor Γ.sub.ml=C.sub.mjkln.sub.jn.sub.k this equation becomes the Christoffel equation:
ρV.sup.2ê.sub.m=Γ.sub.mlê.sub.l  6
Thus, the polarization and phase velocity of a plane wave propagating in direction {circumflex over (n)} are the eigenvector and eigenvalue of the Christoffel tensor Γ.sub.ml respectively. Since this tensor is symmetric, its eigenvalues are real and its eigenvectors are orthogonal with each other.

(27) Isotropic Solid

(28) For an isotropic solid, the stiffness matrix takes the form:

(29) $\begin{array}{cc}{C}_{\alpha \beta }=\left(\begin{array}{cccccc}{c}_{11}& c_{112}& c_{12}& 0& 0& 0\\ c_{12}& c_{11}& c_{12}& 0& 0& 0\\ c_{12}& c_{12}& c_{11}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{44}& 0\\ 0& 0& 0& 0& 0& c_{44}\end{array}\right)& 7\end{array}$
where c.sub.11=c.sub.12+2c.sub.44. Thus, an isotropic solid has only two independent elastic constants c.sub.12=λ and c.sub.44=μ referred as Lame constants. In mechanical engineering literature two other constants are employed, the Young's modulus Y and the Poisson's ratio σ. They are related to Lamé constants by:

(30) $\begin{array}{cc}Y=\mu \frac{3\lambda +2\mu }{\lambda +\mu }& 8\\ \sigma =\frac{\lambda }{2\lambda +\mu }& 9\end{array}$
An incompressible elastic solid is defined as a solid with Poisson's ratio σ≈½. In terms of Lame constants this means λ>>μ and thus, c.sub.11>>c.sub.44. The objective of the present invention when applied to isotropic soft-solids is to estimate the value of μ=c.sub.44 from surface waves.

(31) For an isotropic solid, the Christoffel tensor is independent of the propagation direction. It reads:

(32) $\begin{array}{cc}{\Gamma }_{il}=\left(\begin{array}{ccc}{c}_{44}& 0& 0\\ 0& c_{44}& 0\\ 0& 0& c_{11}\end{array}\right)& 10\end{array}$
Therefore, the eigenvalues are V.sub.T=√{square root over (c.sub.44/ρ)} (degenerate) and V.sub.L=√{square root over (c.sub.11/ρ)} corresponding to polarizations perpendicular and parallel to the propagation direction respectively. Thus, V.sub.T is the velocity of transverse or shear waves and V.sub.L is the velocity of longitudinal or compressional waves. Due to the relation c.sub.11>>c.sub.44, V.sub.L>>V.sub.T in a soft-solid. If the material density of the sample is known, the value of c.sub.44 can be estimated by measuring the velocity V.sub.T of shear waves and inverting the first eigenvalue relation written above: c.sub.44=ρV.sub.T.sup.2. Thus, the task has changed to estimate V.sub.T from measuring the surface displacement field.

(33) When inserting the stiffness matrix given in equation (7) into the wave equation (3), the wave equation for isotropic solids is obtained. It can be written in vectorial form as:

(34) $\begin{array}{cc}\rho \frac{{\partial }^2\stackrel{.fwdarw.}{u}}{\partial t^2}=\left(c_{11}-c_{44}\right)\nabla \left(\nabla .Math.\stackrel{.fwdarw.}{u}\right)+c_{44}{\nabla }^2\stackrel{.fwdarw.}{u}& 11\end{array}$
Taking the divergence in the equation above gives:

(35) $\begin{array}{cc}\rho \frac{{\partial }^2\left(\nabla .Math.\stackrel{.fwdarw.}{u}\right)}{\partial t^2}=c_{11}{\nabla }^2\left(\nabla .Math.\stackrel{.fwdarw.}{u}\right)& 12\end{array}$
Thus, ∇.Math.{right arrow over (u)} it is an irrotational wave that propagates with velocity V.sub.L, i.e., the compressional wave velocity. Let D be a characteristic dimension of the sample being tested (e.g. its length or width). If the excitation frequency f.sub.0 of the wave is chosen such that f.sub.0<<V.sub.L/D, then the wavelength of the irrotational wave ∇.Math.{right arrow over (u)} is much larger than D. Therefore, ∇.Math.{right arrow over (u)} it is approximately constant within the sample and it is possible to write ∇(∇.Math.{right arrow over (u)})≅0. Under this condition, the wave equation (11) becomes:

(36) 0$\begin{array}{cc}\rho \frac{{\partial }^2\stackrel{.fwdarw.}{u}}{\partial t^2}\cong c_{44}{\nabla }^2\stackrel{.fwdarw.}{u}& 13\end{array}$
That is, low-frequency waves in soft-solids propagate almost as shear waves. This last assertion is true for bulk wave propagation in an infinite solid. If the sample is limited by a free surface, surface waves also propagate.

(37) Without loss of generality, consider a surface wave propagating along the x.sub.1 direction. The displacement field is given by:
u.sub.m=ψ.sub.m(kx.sub.3)e.sup.i(ωt−kx.sup.1.sup.)  14
where m=1,3 and u.sub.2 ≡0. Inserting (14) into the wave equation (3) and using the stiffness matrix for an isotropic solid (7), gives:
(ρV.sup.2−c.sub.11)ψ.sub.1(kx.sub.3)+c.sub.44ψ.sub.1″(kx.sub.3)−(c.sub.12+c.sub.44)ψ.sub.3′(kx.sub.3)=0  15
(ρV.sub.2−c.sub.44)ψ.sub.3(kx.sub.3)+c.sub.11ψ.sub.3″(kx.sub.3)−i(c.sub.12+c.sub.44)ψ.sub.1′(kx.sub.3)=0  16
where V=ω/k is the velocity of the surface wave and the primes over the functions indicates derivative with respect to their argument. The solution to this system is given by:
ψ.sub.1(kx.sub.3)=A.sub.1e.sup.−χ.sup.1.sup.kx.sup.3+iχ.sub.3A.sub.3e.sup.−χ.sup.3.sup.kx.sup.3  17
ψ.sub.3(kx.sub.3)=−iχ.sub.1A.sub.1e.sup.−χ.sup.1.sup.kx.sup.3+A.sub.3e.sup.−χ.sup.3.sup.kx.sup.3  18
in which the coefficients A.sub.l and A.sub.3 are determined by the mechanical boundary conditions and

(38) $\begin{array}{cc}{\chi }_1=\sqrt{1-{\left(\frac{V}{V_L}\right)}^2};{\chi }_3=\sqrt{1-{\left(\frac{V}{V_T}\right)}^2}& 19\end{array}$

(39) For a free surface, the relevant boundary conditions at x.sub.3=0, are: τ.sub.m3=0, where m=1,3. Using (1) and (7), the boundary conditions read:
2χ.sub.1A.sub.1+i(1+χ.sub.3.sup.2)A.sub.3=0  20
i(−c.sub.12+c.sub.11χ.sub.1.sup.2)A.sub.1+χ.sub.3(c.sub.12−c.sub.11)A.sub.3=0  21
In order to avoid the trivial solution, the determinant of the coefficient matrix must be zero. This gives rise to the secular equation for the surface waves:
2χ.sub.1χ.sub.3(c.sub.12−c.sub.11)+(1+χ.sub.3.sup.2)(c.sub.11χ.sub.1.sup.2−c.sub.12)=0  22
Replacing in this equation the values of χ.sub.1, χ.sub.3, c.sub.11 and c.sub.12 in terms of the velocities, the secular equation becomes the well-known Rayleigh equation:

(40) $\begin{array}{cc}\frac{V}{V_T}\left[{\left(\frac{V}{V_T}\right)}^6-8{\left(\frac{V}{V_T}\right)}^4+\left(24-\frac{16V_T}{V_L}\right){\left(\frac{V}{V_T}\right)}^2-16\left(1-{\left(\frac{V_T}{V_L}\right)}^2\right)\right]=0& 23\end{array}$
It is known that if V.sub.L/V.sub.T>1.8, this equation has one real root and two complex conjugate roots. The real root corresponds to the velocity V.sub.R of a propagative surface wave (Rayleigh wave). The complex roots correspond also to a physical surface wave termed as the leaky surface wave with propagation velocity V.sub.LS. By using the relation V.sub.L>>V.sub.T in equation (23), it is found that:
V.sub.R=0.96V.sub.T;V.sub.LS=(1.97±i0.57)V.sub.T  24
Therefore, the Rayleigh velocity has a simple relation with the shear velocity: V.sub.R≅0.96V.sub.T. Thus, many authors use, as inversion algorithm to estimate c.sub.44 from surface wave measurements, the following expression:
c.sub.44=ρ(V.sub.R/0.96).sup.2  25
However, this relationship only holds in a semi-infinite medium and in the far-filed where the leaky wave is negligible due to its complex velocity. Real samples are, of course, finite. Thus, in order to meet the conditions for equation (25) to be valid, many authors employ high frequency wave propagation. In this way, the wavelength of the Rayleigh wave is much lesser than the sample's height and the medium is considered semi-infinite Nevertheless, the use of high frequencies is not desirable for many reasons. Firstly, because the surface wave only senses a small portion of the whole sample (its penetration depth is limited to one wavelength). In addition, real samples are attenuating. The attenuation coefficient grows as a power of frequency. Thus, the higher the frequency, the smaller the propagation distance of the wave. Finally, for high frequencies, the compressional wave is not negligible and guided wave propagation is not avoidable.

(41) If the wavelength of the shear wave is comparable to the height of the sample being tested, the inversion algorithm used in equation (25) is no longer valid. In such case, many authors appeal to the Rayleigh-Lamb model of guided wave propagation in plates. Since low frequencies are employed, inversion algorithms based on the zeroth order modes are used because these are the only modes without a cutoff frequency. However, data is usually collected in the near-field of the source and for short times. Within these conditions, the Rayleigh-Lamb modes are not yet developed. This fact leads to biases in the estimation of V.sub.T. Therefore, in the present invention, a different inversion algorithm is used to retrieve the shear wave velocity V.sub.T from surface displacements.

(42) Consider a homogeneous soft-solid isotropic elastic plate of height 2h which is excited by a source located at the free surface in x.sub.3=0 as displayed in FIG. 6. The source vibrates normally to the free surface of the plate, i.e. in the x.sub.3 direction. If the vibration frequency is chosen appropriately, only shear waves propagate into the bulk of the plate as expressed in equation (13). The directivity diagram H(θ) of the shear wave for a source acting normally to the free surface is given by:

(43) $\begin{array}{cc}H\left(\theta \right)=\frac{\cos \left(\theta \right)\sin \left(2\theta \right)\left[{\left(\frac{V_L}{V_T}\right)}^2{\sin }^2\left(\theta \right)-1\right]}{{\left(2{\sin }^2\left(\theta \right)-1\right)}^2-4{\sin }^2\left(\theta \right)\sqrt{{\sin }^2\left(\theta \right)-1}\sqrt{{\sin }^2\left(\theta \right)-{\left(\frac{V_T}{V_L}\right)}^2}}& 26\end{array}$
and displayed in FIG. 5. The angle for the maximum value in the lobe is 34° with respect to the source's axis. The displacement amplitude of the wave decreases rapidly for angles below 34°. Thus, only waves propagating above this direction will produce appreciable reflected waves. This wave reaches the opposite surface at x.sub.3=2h and reflects back into the medium with the same angle. It reaches the free surface x.sub.3=0 at a distance x.sub.s from the source given by: x.sub.s=4h tan 34°≅2.7h. Since the compressional wave is negligible, the surface displacement field at x.sub.3=0 and 0<x.sub.1<x.sub.s has not contributions from the opposite surface. Thus, within this region only surface waves propagate. This include the Rayleigh wave but also the leaky surface wave. The leaky surface wave attenuates as it propagates along the x.sub.1 direction. Using the values of V.sub.LS from equation (24), it is found that its phase velocity is twice the shear wave velocity and its characteristic attenuation distance is given by:

(44) $\begin{array}{cc}\zeta =\frac{1}{\mathrm{Im}\left[\frac{\omega }{V_{LS}}\right]}\cong 1.1{\lambda }_T& 27\end{array}$
where λ.sub.T is the shear wavelength. Thus, the attenuation distance of the leaky surface wave along the x direction is comparable to the shear wavelength in soft-solids. Therefore, the inversion algorithm expressed in equation (25) is valid only if ζ<x.sub.1<x.sub.s. However, the value of is not negligible at low frequencies. For example, typical values of the involved velocities in soft-tissues are V.sub.L˜1500 m/s and V.sub.T˜10 m/s. Therefore, the shear wavelength may vary from 10 to 1 cm for excitation frequencies in the range 100-1000 Hz. Then, the leaky wave is not negligible in the near-field. Therefore, interference between the Rayleigh and the leaky surface wave is possible. This fact has consequences over the phase velocity of the surface field as shown below.

(45) Note that, depending on frequency, the value of can be greater than x.sub.s. A critical frequency f.sub.c can be defined such as ζ=x.sub.s. The value of f.sub.c is given by:

(46) $\begin{array}{cc}1.1\frac{V_T}{f_c}=2.7h.Math.f_c=0.4\frac{V_T}{h}& 28\end{array}$

(47) From the considerations above, the x.sub.3-component of the surface displacement u.sub.3 (x.sub.1<x.sub.s, x.sub.3=0) for frequencies below f.sub.c can be written as:
u.sub.3(x.sub.1<x.sub.s)=exp(−ik.sub.Rx.sub.1)−exp(−x.sub.1/ζ)exp(−ik.sub.LSx.sub.1)  29
corresponding to the sum of the Rayleigh and leaky surface wave, where k.sub.Ls is the real part of the wavenumber of the leaky surface wave. For low frequencies, k.sub.Rx.sub.1<<1 and therefore:

(48) $\begin{array}{cc}{u}_3\left(x_1<x_s\right)\cong \left[1-\exp \left(-x_1\text{/}\zeta \right)+\frac{k_R^2}{2}\left({\delta }^2\exp \left(-x_1\text{/}\zeta \right)-1\right)x_1^2\right]+{\mathrm{ik}}_R\left[\left(1-\delta \exp \left(-x_1\text{/}\zeta \right)\right)x_1+\left(\frac{{\delta }^3\exp \left(-x_1\text{/}\zeta \right)-1}{6}\right)k_R^2{x}_1^3\right]& 30\end{array}$
where δ is defined as δ=k.sub.LS/k.sub.R≅½. Defining

(49) $N\left(x_1,\omega ,V_T\right)=k_R\left[\left(1-\delta \exp \left(-x_1\text{/}\zeta \right)\right)x_1+\left(\frac{{\delta }^3\exp \left(-x_1\text{/}\zeta \right)-1}{6}\right)k_R^2{x}_1^3\right],\mathrm{and}$$D\left(x_1,\omega ,V_T\right)=1-\exp \left(-x_1\text{/}\zeta \right)+\frac{k_R^2}{2}\left({\delta }^2\exp \left(-x_1\text{/}\zeta \right)-1\right)x_1^2,$
the phase ϕ(x.sub.1, ω, V.sub.T) of this wave can be expressed as:

(50) $\begin{array}{cc}\varphi \left(x_1,\omega ,V_T\right)=\mathrm{atan}\left(\frac{N\left(x_1,\omega ,V_T\right)}{D\left(x_1,\omega ,V_T\right)}\right)& 31\end{array}$
Thus, the phase velocity V of the surface wave for a given frequency ω.sub.0 is expressed as:

(51) $\begin{array}{cc}V\left(x_1,{\omega }_0,V_T\right)={{\omega }_0\left(\frac{\partial \varphi }{\partial x_1}\right)}^{-1}={\omega }_0\left(\frac{D^2\left(x_1,{\omega }_0,V_T\right)+N^2\left(x_1,{\omega }_0,V_T\right)}{N^{\prime }\left(x_1,{\omega }_0,V_T\right)D\left(x_1,{\omega }_0,V_T\right)-N\left(x_1,{\omega }_0,V_T\right)D^{\prime }\left(x_1,{\omega }_0,V_T\right)}\right)& 32\end{array}$
where the primes over the function indicates derivative with respect to x.sub.1. Note that this expression gives a different dispersion curve for each value of x.sub.1. Thus, it is not the dispersion curve for Rayleigh-Lamb modes. At this stage, two inversion methods are envisaged to retrieve V.sub.T from the experimental values of V. First inversion method. If, for a given position x.sub.1<x.sub.s, an experimental dispersion curve V(ω) for frequencies ω>ω.sub.c=2πf.sub.c is available, then equation (32) can be fitted in a least-squares sense to the experimental data. The value of V.sub.T that minimizes the sum of quadratic differences is the best estimation for the shear wave velocity. This was the inversion method employed in FIG. 7. Second inversion method. If, for a given position x.sub.1<x.sub.s, a single value of V is available, corresponding to a single frequency ω.sub.0>ω.sub.c, then equation (32) can be used to build the curve V(V.sub.T). The abscissa corresponding to the intersection between the experimental and theoretical values gives the estimation for the shear wave velocity V.sub.T.

(52) So far, the inversion methods proposed above only consider the surface displacement field for x.sub.1<x.sub.s. If x.sub.1>x.sub.s, the reflected shear waves are not negligible and must be taken into account in the inversion method. To this end, the effects of an impinging shear wave on the free surface must be considered first.

(53) Consider an impinging shear wave at the free surface x.sub.3=0. This wave produces a reflected shear wave and, due to mode conversion, a reflected compressional wave. If A.sub.s.sup.i and θ.sub.s.sup.i are the amplitude and angle of the incident shear wave, the amplitude A.sub.p.sup.r and angle θ.sub.p.sup.r of the reflected compressional wave are given by:

(54) 0$\begin{array}{cc}\sin \left({\theta }_p^r\right)=\frac{V_L}{V_T}\sin \left({\theta }_s^i\right)& 33\\ A_p^r=-\frac{2{\left(V_L\text{/}{V}_T\right)}^2\sin \left(2{\theta }_s^i\right)\cos \left(2{\theta }_s^i\right)}{\sin \left(2{\theta }_p^r\right)\sin \left(2{\theta }_s^i\right)+{\left(V_L\text{/}{V}_T\right)}^2{\cos }^2\left(2{\theta }_s^i\right)}{A}_s^i& 34\end{array}$
Since V.sub.L>>V.sub.T, the angle of is complex even if θ.sub.s.sup.i<<1. Let θ.sub.r.sup.p=π/2−iγ where γ is real and positive. Consider now the components of the wavevector {right arrow over (k)}.sup.p of the reflected compressional wave along x.sub.1 and x.sub.3:
k.sub.1.sup.p=k.sup.p sin(θ.sub.p.sup.r)=k.sup.p sin(π/2−iγ)=k.sup.p cos h(γ)
k.sub.3.sup.p=k.sup.p cos(θ.sub.p.sup.r)=k.sup.p cos(π/2−iγ)=ik.sup.p sin h(γ)
Thus, whatever the incident angle, the impinging shear wave produces an evanescent compressional wave confined to the free surface of the soft-solid. Thus, it is a surface wave which we call the SP wave (FIG. 6). The phase velocity V.sub.SP of this wave along the x.sub.1-direction can be computed as:

(55) $V_{SP}=\frac{\omega }{k_1^p}=\frac{V_T}{\sin \left({\theta }_s^i\right)}$
Therefore, the phase velocity of the SP wave varies from nearly infinite for θ.sub.s.sup.i≅0 to V.sub.T for θ.sub.s.sup.i≅π/2. After some computation, the amplitude of the SP wave is given by:

(56) $A_p^r\cong -\frac{2\sin \left(2{\theta }_s^i\right)\cos \left(2{\theta }_s^i\right)}{2i{\sin }^2\left({\theta }_s^i\right)\sin \left(2{\theta }_s^i\right)+{\cos }^2\left({\theta }_s^i\right)}{A}_s^i$
where the approximation is valid if γ>>1. This is always the case since θ.sub.s.sup.i≥34°. Thus, due to the complex denominator, the SP wave is out of phase with the reflected shear wave. The phase difference depends upon the incident angle θ.sub.s.sup.i of the shear wave. We expect to observe the SP wave whenever the amplitude A.sub.s.sup.i of the shear wave is maximum. According to the directivity pattern of the shear wave, it has a maximum for θ.sub.s.sup.i≅34°. At this angle, the phase velocity of the SP wave is V.sub.Sp≅1.8V.sub.T. Thus, the surface displacement field in the vicinity of x.sub.s is a complicated superposition of different wave types, each with its own phase velocity. If ω>ω.sub.c the x.sub.3 component of the surface field can be expressed as the superposition of the Rayleigh wave, the reflected shear wave and the SP wave as:
u.sub.3(x.sub.1≈x.sub.s)=e.sup.−ik.sup.R.sup.x.sup.1+A.sub.s(x.sub.1)[e.sup.−ik.sup.s.sup.1.sup.x.sup.1+|A.sub.r.sup.p(θ.sub.s)|e.sup.−i(k.sup.sp.sup.x.sup.1.sup.+η)]  35
where k.sub.s.sup.1 is the component of the shear wave vector along x.sub.1, η is the phase angle of the SP wave with respect to the shear wave and A.sub.s(x.sub.1) is the amplitude of the reflected shear wave at θ≅θ.sub.s=34°. Taking into account the directivity and the geometrical attenuation, it is given by:

(57) $\mathrm{As}\left(x_1\right)=\frac{H\left({\theta }_s\right)}{{\left(x_1^2+16h^2\right)}^{1\text{/}2}}$
Now, for low frequencies, equation (35) can be expressed as:
u.sub.3(x.sub.1≈x.sub.s)≅1+A.sub.s(x.sub.1)|A.sub.r.sup.p(θ.sub.s)|(cos(η)−k.sub.sp sin(η))−½((k.sub.R.sup.2+A.sub.s(x)(k.sub.s.sup.1).sup.2+A.sub.s(x.sub.1)|A.sub.r.sup.p(θ.sub.s)|k.sub.sp.sup.2 cos(η))x.sub.1.sup.2+⅙A.sub.s(x.sub.1)|A.sub.r.sup.p)(θ.sub.s)|k.sub.sp.sup.2 sin(η)x.sub.1.sup.3−i[A.sub.s(x.sub.1)|A.sub.r.sup.p(θ.sub.s)|sin(η)+(k.sub.R+A.sub.s(x)|A.sub.r.sup.p(θ.sub.s)|k.sub.sp cos(η)+A.sub.s(x)k.sub.s.sup.1)x.sub.1−½A.sub.s(x.sub.1)|A.sub.r.sup.p(θ.sub.s)|k.sub.sp.sup.2 sin(η)x.sub.1.sup.2−−⅙(k.sub.R.sup.3+A.sub.s(x.sub.1)(k.sub.s.sup.1).sup.3+A.sub.s(x.sub.1)|A.sub.r.sup.p(θ.sub.s)|k.sub.sp.sup.3 cos(η))x.sub.1.sup.3]  36
Defining N.sub.s(x)=Im[u.sub.z] and D.sub.s(x)=Re[u.sub.2], the phase velocity is expressed as:

(58) $\begin{array}{cc}V\left(x_1\approx x_s,{\omega }_0,V_T\right)={\omega }_0\left(\left(\frac{D_s^2\left(x_1,{\omega }_0,V_T\right)+N_s^2\left(x_1,{\omega }_0,V_T\right)}{\begin{array}{c}{N}_s^{\prime }\left(x_1,{\omega }_0,V_T\right)D_s\left(x_1,{\omega }_0,V_T\right)-\\ N_s\left(x_1,{\omega }_0,V_T\right)D_s^{\prime }\left(x_1,{\omega }_0,V_T\right)\end{array}}\right)\right)& 37\end{array}$
Thus, a third inversion method is possible:

(59) Third inversion method. If, for a given position x.sub.1≈x.sub.s, an experimental dispersion curve V(ω) for frequencies ω>ω.sub.c is available, then equation (37) can be fitted in a least-squares sense to the experimental data. The value of V.sub.T that minimizes the sum of quadratic differences is the best estimation for the shear wave velocity. This was the inversion method employed in FIG. 8.

(60) A third possibility regarding the position x.sub.1 at which the surface displacement is measured must be considered. If x.sub.1>>x.sub.s a few reflections back and forth in the interfaces of the sample have taken place. Thus, in this zone, the Rayleigh-Lamb modes have developed. Due to the source type and polarization used in this invention, it favors the antisymmetric modes. In addition, since low-frequencies are employed, only the zeroth order mode propagate. Thus, other two inversion methods are possible, in the same spirit of methods #1 and #2: Fourth inversion method. If, for a given position x.sub.1>>x.sub.s, an experimental dispersion curve V(ω) is available, then the zeroth order antisymmetric Rayleigh-Lamb mode can be fitted in a least-squares sense to the experimental data. The value of V.sub.T that minimizes the sum of quadratic differences is the best estimation for the shear wave velocity. This was the inversion method employed in FIG. 9. Fifth inversion method. If, for a given position x.sub.1>>x.sub.s, a single value of V is available, corresponding to a single frequency ω.sub.0, then the zeroth order antisymmetric Rayleigh-Lamb mode can be used to build the curve V(V.sub.T). The abscissa corresponding to the intersection between the experimental and theoretical values gives the estimation for the shear wave velocity V.sub.T.

(61) As a final observation for isotropic solids, we note here that, due to the interference of the different surface wave types, the amplitude of the surface displacement is a complicated function of position x and frequency ω. Thus, an estimation of the attenuation coefficient by fitting the amplitude of the displacement field to the usual exponential decay does not make sense, at least in the near field. If far-field (x.sub.1>>x.sub.s) measurements are available, then, the amplitude should be corrected by diffraction before trying to estimate the attenuation coefficient.

(62) Transversely Isotropic Solid

(63) Potential applications of the present invention include estimating the elastic properties of skeletal semi-tendinous muscles, e.g., meat samples or application in vivo to external human muscles such as biceps, triceps, quadriceps, etc. Due to the parallel fiber orientation in these kind of muscles, they can be modelled as transversely isotropic materials. In addition, some long-chain polymers with chains aligned within a preferred direction are also modelled as transversely isotropic materials.

(64) The stiffness matrix has five independent elastic modulus for this kind of materials. Let x.sub.1 be the axis of symmetry, i.e., the orientation axis of muscular fibers or long molecule-chains. Thus, the axes x.sub.2 and x.sub.3 are perpendicular to the fibers as shown in FIG. 4. For this choice of axes orientation, the stiffness matrix reads:

(65) $\begin{array}{cc}{C}_{\alpha \beta }=\left(\begin{array}{cccccc}{c}_{11}& c_{12}& c_{12}& 0& 0& 0\\ c_{12}& c_{22}& c_{22}-2c_{44}& 0& 0& 0\\ c_{12}& c_{22}-2c_{44}& c_{22}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& 0\\ 0& 0& 0& 0& 0& c_{55}\end{array}\right)& 38\end{array}$
The constants c.sub.11 and c.sub.22 are related to the compressional wave propagating along and perpendicular to the symmetry axis respectively: V.sub.L.sup.∥=√{square root over (c.sub.11/ρ)} and V.sub.L.sup.⊥=√{square root over (c.sub.22/ρ)}. In many incompressible transversely isotropic solids (such as skeletal muscle for example), these two values of compressional wave velocities are almost equal each other. Thus, c.sub.22≅c.sub.11, i.e. the solid is isotropic regarding the compressional waves. The constant c.sub.44 is related to a shear wave propagating perpendicular to the symmetry axis with perpendicular polarization V.sub.T.sup.⊥=√{square root over (c.sub.44/ρ)}. The constant c.sub.55 is related to a shear wave propagating parallel to the symmetry axis with perpendicular polarization V.sub.T.sup.⊥=√{square root over (c.sub.55/ρ)}. As in the isotropic solid case, V.sub.L>>V.sub.T, whatever the polarization. Thus, c.sub.22−2c.sub.44≅c.sub.22≅c.sub.11. Finally, the constant c.sub.12 is related to wave propagation (either compressional or shear waves) in directions out of the principal axes.

(66) Consider now an anisotropic soft-solid where V.sub.L>>V.sub.T whatever the propagation direction of the waves. The wave equation for this case cannot be written in a simple manner as for the isotropic case. However, it is still valid that the contribution of the compressional waves are negligible at low-frequencies. Thus, low-frequency waves propagate almost as shear waves. The objective of the present invention, when applied to transversely isotropic solids, is to estimate either c.sub.44, c.sub.55 or both of them.

(67) A. Estimation of c.sub.44

(68) If wave propagation takes place in a plane perpendicular to the symmetry axis (plane (x.sub.2, x.sub.3) in FIG. 4), the shear wave velocity is V.sub.T.sup.⊥ and is independent of the propagation direction within this plane. Using the aforementioned relations between elastic constants for the soft solid case, i.e., c.sub.22=c.sub.11 and c.sub.22−2c.sub.44=c.sub.11, the Rayleigh wave equation takes exactly the same form as in the isotropic case with V.sub.L=√{square root over (c.sub.11/ρ)} and V.sub.T.sup.⊥=√{square root over (c.sub.44/ρ)}. Therefore, if the source does not excite the u.sub.1 component of the field, the inversion method to estimate c.sub.44 proceeds exactly as for the isotropic soft-solid. To achieve this condition, the invention is equipped with a line source (item(1), FIG. 3). The length of the line source is much longer than the shear wavelength (i.e. the line source is “infinite”). If the source is aligned parallel to the x.sub.1 direction, it excites the u.sub.2 and u.sub.3 components of the field, but not the u.sub.1 component.

(69) B. Estimation of c.sub.55

(70) For estimating c.sub.55, consider the line source oriented parallel to the x.sub.2 axis. For this case, the wave propagation takes place in a plane that includes the symmetry axis (plane (x.sub.1, x.sub.3) in FIG. 4). The secular equation for the Rayleigh wave under this condition is given by:

(71) $\begin{array}{cc}\left(\vartheta -\rho V^2\right)\sqrt{1-\frac{\rho V^2}{c_{55}}}-\rho V^2=0& 39\end{array}$
where ∂=c.sub.11+c.sub.22−2c.sub.12. Since the medium is considered isotropic for compressional waves, the constant c.sub.12 is no longer an independent constant. It is related to other constants by:
c.sub.12=c.sub.11−c.sub.55
Since c.sub.22=c.sub.11, the secular equation (39) is expressed as:

(72) $\begin{array}{cc}\left(2-{\left(\frac{V}{V_T^{\parallel }}\right)}^2\right)\sqrt{1-{\left(\frac{V}{V_T^{\parallel }}\right)}^2}-{\left(\frac{V}{V_T^{\parallel }}\right)}^2=0& 40\end{array}$
As for the isotropic case, this equation has one real root (corresponding to the Rayleigh wave) and two complex conjugate roots corresponding to the leaky surface wave:
V.sub.R=0.84V.sub.T.sup.∥;V.sub.LS=(1.42±0.6i)V.sub.T.sup.∥  41

(73) Note that the attenuation distance for the leaky wave ζ is larger than for the isotropic case, ζ=1.6λ.sub.T. Another difference with respect the isotropic solid concerns the directivity pattern of the shear wave. There is no simple expression for the directivity pattern in the (x.sub.1, x.sub.3) plane. However, some authors have computed it numerically. It is shown that the main lobe is oriented towards 60° from the source. Therefore the value of x.sub.s is now given by x.sub.s≅7h. Thus, the value of the critical frequency at which ζ=x.sub.s is given by f.sub.c=0.23V.sub.T/h, which is lower than the isotropic case. Therefore, if x.sub.1<x.sub.s and ω>ω.sub.c=2πf.sub.c, the inversion method to obtain c.sub.55 proceeds as given by equation (32) for the isotropic solid but changing the values of V.sub.R and V.sub.LS by the ones given in equation (41).

(74) Inversion procedures for x.sub.1≥x.sub.s are complicated since no analytic expressions are available for the directivity of shear waves. In addition, Rayleigh-Lamb modes in the (x.sub.1, x.sub.3) plane are difficult to compute even numerically. Therefore, the present invention does not include inversion methods for these cases.

(75) The lack of inversion algorithm for x.sub.1>x.sub.s, does not forbids the estimation of c.sub.55 in common applications of the invention. Consider for example skeletal muscles such as biceps branchii. The mean value of its height in adults is 2h=28±7 mm. Thus, x.sub.s≅7h≅100 mm. The mean height for other skeletal muscles like vastus lateralis or vastus medialis is even larger than 28 mm. Therefore, if the data is collected at positions x<100 mm, the value of c.sub.55 can be estimated by the inversion methods proposed in the present invention.

Application Examples

(76) Temperature Dependence of c.sub.55

(77) The elasticity of soft tissues depends on the temperature. In order to compare the elasticity of different beef samples it is important that all are taken at the same temperature. However, this is not always the case. The temperature in slaughter houses can vary within a range between 3 and 10° C. In this example, the temperature dependence of c.sub.55 in beef samples is shown in FIGS. 11A-11B. A consisted decrease in elasticity with temperature rise is observed in the samples tested.

(78) Age Maturation Process

(79) Age maturation is a fundamental process in meat industry. It is mediated by the action of many enzymatic systems. After rigor mortis, these enzymatic systems produce a progressive softening of beef allowing to reach the desired tenderness requested by consumers. There is a lack of nondestructive monitoring method of enzymatic maturation. In this example, the maturation process of 5 different cuts, consisting in 5 samples of each, is monitored with the present invention. The samples were kept vacuum sealed inside a cold room between 0 and 4° C. during 21 days. All samples were measured once a day. FIGS. 12A to 12E display the normalized elasticity with respect to the first day for each cut. The errorbars represent one standard deviation of the 5 samples for each cut. After initial fluctuations, a decrease in elasticity is observed for all samples. Roughly, an asymptotic value is reached around 14 days of maturation for all samples.

(80) Comparison with Warner-Bratzler Shear Force Test (WBSF)

(81) The WBSF test is a destructive method. It measures the force that a sharp inverted “V” shaped knife needs to cut the beef sample while moving at constant speed. It is the standard test in meat industry to quantify the tenderness. In this example, the results of WBSF tests in two different beef samples are shown (FIG. 13A). The results of the present invention applied to the same samples are shown in FIG. 13B. FIG. 13C displays the ratio between the elastic modulus of the two samples. This ratio is in agreement with the ratio expected from the WBSF tests (FIG. 13A)

(82) Sorting Beef Cuts According to Tenderness

(83) The tenderness of beef samples represents a variable that defines, in a high percentage, its commercial value. The standard procedure to evaluate tenderness is the WBSF test. In the previous example, it was shown that the shear force test and the elasticity are correlated. Thus, the present invention can be employed to discriminate beef samples according to their tenderness. FIG. 14 displays the elasticity value of different beef cuts and their classification according to 4 grades of tenderness: tough, medium tough, tender and extremely tender.

(84) Elasticity Estimation in Skeletal Muscle in Vivo

(85) Estimation of elastic properties of skeletal muscle in dynamic conditions is important in sport science to evaluate the health state of muscles. Thus, the present invention would be beneficial in this and related areas. In this example, the present invention is used to monitor changes in c.sub.55 in biceps brachial of three healthy volunteers during isometric contraction. The volunteers were asked to contract the muscle to 40% of their maximum voluntary contraction (MVC, measured previously with an isokinetic dynamometer) in 20 seconds. Then, keep this contraction for other 5 seconds and finally reduce the contraction to 0% in 20 seconds. FIGS. 15A-15C show the obtained results for each volunteer where it can be observed that the present invention is capable of following elasticity changes in dynamic situations.

(86) Determination of V.sub.T in Tissue Mimicking Phantoms

(87) In this section we present experimental results for the application of the present invention in agar-gelatin based phantoms (isotropic solid). These types of phantoms are widely used to simulate the mechanical properties of soft tissues. They were made from a mixture of agar (1% w/w) and gelatin (3% w/w) in hot water (>80° C.). Alcohol and antibiotics were also added for conservation purposes. Thus, we elaborated a sample with parallelepiped shape of height 2h=20 mm, long and width 100 and 120, respectively. FIGS. 7, 8 and 9 display the dispersion curve for the phase velocity of the surface wave, obtained for the configurations x.sub.1<x.sub.s, x.sub.1≈x.sub.s y x.sub.1>2x.sub.s, respectively. Each point of these curves was obtained by exciting the surface actuator with 6 cycles of sinusoid at the corresponding frequency. Then the phase velocity V for each receiving channel was computed. In this way, by adjusting the experimental data to the theoretical curves provided by the appropriate inversion method according to the experimental configuration, the corresponding values of V.sub.T were obtained.