The computations of reflection coefficients of multilayer structure based on the reformulation of Thomson-Haskell method

Ultrasonics 52 (2012) 1019–1023

Contents lists available at SciVerse ScienceDirect

Ultrasonics journal homepage: www.elsevier.com/locate/ultras

The computations of reflection coefficients of multilayer structure based on the reformulation of Thomson-Haskell method Jian Chen, Xiaolong Bai, Keji Yang, Bing-Feng Ju ⇑ The State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou 310027, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 10 April 2012 Received in revised form 21 July 2012 Accepted 2 August 2012 Available online 25 August 2012 Keywords: Multilayer structure Ultrasonic NDT technique Reformulation of Thomson-Haskell method Reflection coefficients

a b s t r a c t A reformulation of the Thomson-Haskell method is presented for calculating the reflection coefficients of multilayer structure immersing in the coupling fluid. Instead of directly multiplying the layer propagator matrix, the new method splits the layer propagator matrix and excursively determines the interface stiffness matrix starting from the bottom half-space with known stiffness. A formulation for the reflection coefficients is derived based on the obtained interface stiffness matrix of the top layer. This scheme can be applied to a single solid layers or layered structures containing both fluid and solid layers. It keeps the simplicity but naturally excludes the exponential growth term and thus can be applied at any frequency range. Its validity and feasibility were experimentally proved by the measurement of the reflection coefficients of a three layered structure of aluminum–glass–aluminum and a sandwiched layer structure of two 250 lm stainless plates filled with 100 lm deionized water based on the inversion of V(z, t) technique. The result of experiments is consistent with the theoretical calculation. The reformulation of the Thomson-Haskell method offers an efficient and effective solution for calculating the acoustic reflection coefficients of multilayer structures of any configurations. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Multilayer structures have been widely used in safety–critical applications because of their efficiency and advantage in weight reduction and load distribution, for example, the primary and secondary structures of aircraft, submarine, ultra-thin foils and even some kinds of new biomaterials [1]. Ultrasonic non-destructive evaluation the acoustical material properties of multilayer structure have varieties applications in judging the adhesive quality, characterizing the protective coatings, determining the uniformity of ultra-thin foils and so on. The Thomson-Haskell transfer matrix method was firstly developed by Thomson and later refined by Haskell, which provided a systematic treatment of initial-boundary value problems. The study of wave propagation in layered structure has received much attention and found wide applications [2]. However, for the cases of high frequency or thick layer, as the matrix components become very large and the secular function loses significant figures, it is unable to obtain accurate values of its roots and will lead instability in high frequencies computation [3]. To eliminate this limitation, many customary remedies, including the fast delta matrix methods, the fast Schwab–Knopoff method, the Abo-Zena method, the reflection and transmission matrix method and generalized RT

⇑ Corresponding author. Tel./fax: +86 571 8795 1730. E-mail address: [email protected] (B.-F. Ju). 0041-624X/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ultras.2012.08.004

coefficient method, were proposed [4,5]. Another alternative method is the global matrix method, originally proposed for layered isotropic media and if properly implemented, it is unconditionally stable. Wang et al. further reformulated as a recursive stiffness matrix and demonstrated that the computation efficiency is the same as that of the original Thomson-Haskell transfer matrix method [6,7]. For problems where not all submatrices of the stiffness matrix are needed, e.g., for a layered half space, an impedance matrix method can be utilized [8]. In this paper, a reformulation of Thomson-Haskell method was presented to compute the acoustic reflection coefficients of layered structure. It splits the layer propagator matrix and determines the interface stiffness, and has many advantages over other available complex arithmetic methods. Firstly, it keeps the original simplicity and achieves the same stability as those of the global matrix method and the fast-generalized RT method for it naturally excludes the exponential growth terms. Secondly, it requires only the elementary matrices constituting the propagator matrix, not necessary for the layer stiffness matrix. Finally, due to its concise sub-matrix form, the reformulation of Thomson-Haskell method is more efficient than other available complex arithmetic methods, and the computational speeds rise by 20%, comparing to the original Thomson-Haskell method and fast-generalized RT method [2]. To verify its validity and feasibility, an experimental measurement of the reflection coefficients for the layered structures is carried out, which are consistent with the corresponding theoretical calculations. Since the acoustic reflection coefficients involve a wealth of

1020

J. Chen et al. / Ultrasonics 52 (2012) 1019–1023

readily available information of material acoustic properties, the computation of reflection and transmission coefficients of the layered structure as a help for analysis and interpretation is of great importance [9], the reformulation of Thomson-Haskell method offers a new NDT way to evaluate the multilayer structures.

2

1

6 cam 6 Tm ¼ 6 4 qm cm cam

cbm

1

1

cam

7 7 7: qm ð1  cm Þ 5 1

ca

qm ð1  cm Þ qm cm m qm cm cbm qm ð1  cm Þ

qm ð1  cm Þ

3

cbm

qm cm cbm ð4Þ

2. The reformulation of Thomson-Haskell method

And

A typical multilayer structure is schematically shown in Fig. 1, in which the layers could be fluid, solid, or even be mixed solid–liquid layers. Each layer is assumed to be homogeneous and with uniform elasticity. The arbitrary mth layer is the one between planes zm1 and zm, which has the density qm, the elastic constants km and lm, and the thickness hm. It has longitudinal and transverse pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi wave speed of cpm ¼ ðkm þ 2lm Þ=qm and csm ¼ lm =qm , respectively. The dilatational and shear potentials £m and /m in the mth layer which is related to the state vector is given as:

Kþm ðzÞ ¼ diagðeiam z ; eibm z Þ;

ð5aÞ

Km ðzÞ ¼ diagðeiam z ; eibm z Þ;

ð5bÞ

;m ¼ am e

iðkx xþam ðzzm1 ÞÞ

um ¼ cm e

iðkx xþbm ðzzm1 ÞÞ

iðkx xam ðzzm1 ÞÞ

þ bm e

;

ð1Þ

iðkx xbm ðzzm1 ÞÞ

ð2Þ

þ dm e

;

2

2

where am ¼ ððx=cpm Þ2  kx Þ1=2 ; bm ¼ ððx=csm Þ2  kx Þ1=2 ; kx is the x component of wave number k0 of water. For a plain-strain problem, the state vector of the mth layer at any point in the layer is given as:



um ðzÞ

rm ðzÞ where

"

 ¼

T 1m

T 2m

T 3m

T 4m

#"

Kþm ðzÞ 0 0 Km ðzÞ

#"

þ

;m 

;m

# ¼ T m Km ; m ;

um ¼ ðikx Þ ½uxm uzm T ; m ¼ ð Þ2 ½ xm zm T , ½um m T is þ  and ½;m ;m T ¼ ½am cm bm dm T is down-going and 1

r

x

r r

r

ð3Þ the state

up-going vector wave amplitude vector. For the solid layers, the terms of Tm can be expressed as [3]:

where 1

am ¼ ððx=cpm Þ2  k2x Þ2 ; bm ¼ ððx=csm Þ2  k2x Þ1=2 ;

ð6Þ

and

cm ¼ 2ðkx =ks Þ2 ; ks ¼ x=csm ; cam ¼ am =kx ; cbm ¼ bm =kx

ð7Þ

For the fluid layer, the terms of Tm can be expressed as follows:

2

0 6 ca 6 Tm ¼ 6 m 4 0 qm

0

0

0 0

ca

0

3

07 7 7; 05

m

0

0 qm

ð8Þ

0

and

Kþm ðzÞ ¼ diagðeiam z ; 0Þ; Km ðzÞ ¼ diagðeiam z ; 0Þ

ð9Þ

Through Eq. (3), the transfer relationship of the state vectors at two opposite interfaces of the mth layer can be established as follows:



uum

"

 ¼

rum

T 1m

T 2m

T 3m

T 4m

#

"

Em

0

0

Eþm

T 1m

T 2m

T 3m

T 4m

#1 "

udm

rdm

#

" ¼ Bm

udm

#

rdm ð10Þ

where the superscripts l and d represent the upper and lower interfaces of an arbitrary layer, Eþ= ¼ Kþ= m m ðhm Þ and Bm is a propagator matrix. For the solid layer, the analytical solution for the inverse matrix T 1 m is as follows: 3 2 cm ðcm  1Þ=cam 1=ðqm cam Þ 1=qm b cm 1=qm 1=ðqm cbm Þ 7 16 7 6 ðcm  1Þ=cm T 1 7; 6 m ¼ a a 4 2 cm ðcm  1Þ=cm 1=ðqm cm Þ 1=qm 5

ðcm  1Þ=cbm

cm

1=qm

1=ðqm cbm Þ ð11Þ

and

3 0 1=qm 0 0 7 7 7: 0 1=qm 5

2

T 1 m

0 1=cam 6 0 160 ¼ 6 2 4 0 1=cam 0

0

0

ð12Þ

0 T 1 m

For the case of fluid layer, is the Moore–Penrose pseudoinverse matrix of Tm rather than the inverse of Tm. The continuity conditions between the mth and (m + 1)th layers imply ½udm rdm T ¼ ½uumþ1 rumþ1 T . This immediately leads to the Thomson-Haskell method.



uu1

ru1

Fig. 1. Geometry and coordinate system of the multilayer structure.



¼ B1 B2 ;    ; Bn



uun

run

 ð13Þ

Starting from the bottom half-space, the recursive algorithm can be applied to progressively derive the boundary value at the top layer. þ Since the propagator matrix Bm depends on both terms of E m and Em , the recursive algorithm would have numerical difficulties with the

1021

J. Chen et al. / Ultrasonics 52 (2012) 1019–1023

Fig. 2. The schematic view of a scanning acoustic microscope system for acoustic reflection coefficients measurement.

existing evanescent waves, it is an inherent shortcoming of the original Thomson-Haskell method. For the purpose of overcoming the inherent computation instability of original Thomson-Haskell at high frequency, an efficient and effective reformulation of the Thomson-Haskell method is presented. It naturally excludes the exponential growth terms and is unconditionally stable at any frequency range and can be applied to both the solid layer and fluid layer. By introducing the interface stiffness matrix Sm and let rum ¼ Sm uum , it has:



uum



rum



 I ¼ uu : Sm m

ð14Þ

According to Eqs. (10) and (14) and interface continuity conditions, there is:



I



Sm

" uum ¼

T 1m T 3m

T 2m T 4m

#

Em

0

0

Eþm

"

T 1m T 3m

T 2m T 4m

#1 

I Smþ1



uumþ1 ;

ð15Þ

Two intermediate matrices P 1m ; P2m are introduced:

"

P1m

#

" ¼

P2m

T 1m

T 2m

T 3m

T 4m

#1 

I Smþ1

 ;

ð18Þ

and Eq. (17) can be expressed as follows:





I Sm

" uum

¼

Em P1m

#

Eþm P2m

uumþ1

ð19Þ

And then

Sm ¼ Eþm P2m ðP1m Þ1 Eþm

ð20Þ

Combining Eqs. (16) and (20), it can be derived:

Sm ¼ ðT 3m þ T 4m Sm ÞðT 1m þ T 2m Sm Þ1 ;

ð21Þ

where m ¼ n  1; n  2;    ; 2; 1. For the case of fluid layer, the inversion in Eqs. (20) and (21) are Moore–Penrose pseudoinverse matrix.

where m ¼ n  1; n  2;    ; 2; 1 . Again by introducing auxiliary stiffness matrix Sm and auxiliary  um , and splitting Eq. (15) leads to the followdisplacement vector u ing equations:



I



Sm

" uum

¼

T 1m

T 2m

T 3m

T 4m

Em

0

0

Eþm

#

I



Sm

uum ;

ð16Þ

and



I Sm



uum ¼



"

T 1m

T 2m

T 3m

T 4m

#1 



I Smþ1

uumþ1

ð17Þ

Table 1 Acoustic properties of the concerned materials. Materials

Vl (m/s)

Vt (m/s)

q (kg/m3)

Stainless steel Aluminum Glass Water

5650 6350 5640 1480

3060 3230 3020 –

7900 2700 2300 1000

Vl: Longitudinal velocity; Vt: Transverse velocity; q: Density.

Fig. 3. (a) Measured and (b) calculated reflection coefficients R(h, x) for a three layered structure of aluminum–glass–aluminum.

1022

J. Chen et al. / Ultrasonics 52 (2012) 1019–1023

Starting from the bottom half-space with known stiffness Sn, the recursive algorithm deals with interface stiffness matrix of one layer using Eqs. 18, 20, and 21 each time, until the top layer is reached. The proposed method requires five 2  2 matrix multiplications and two 2  2 matrix inversions at each recursion. As for the original Thomson-Haskell method (The same for recursive stiffness matrix method), each recursion requires eight 2  2 matrix multiplications. In the case of the fast-generalized RT method, it requires seven 2  2 matrix multiplications and one 2  2 matrix inversion. Therefore, the proposed algorithm is the more efficient. 3. The computations of reflection coefficients for multilayer structure A plane acoustic wave of unity amplitude is incident on the fluid–solid interface at an angle h as shown in Fig. 1 and the wave potential in the fluid is given by [10]:

;f ¼ eiðkx xþkz zÞ þ Reiðkx xkz zÞ ;

ð22Þ

where kx = kf sin h, kz = kf cos h, and kf = x/Vf. Vf is the longitudinal wave speed of the fluid. For the solid substrate, the stiffness of bottom half-space Sn can be obtained from Eq. (3) as follows:

Sn ¼ T 3n ðT 1n Þ1

ð23Þ

For the case that the solid substrate not exists, i.e. the fluid is the substrate for the case of the layered structure is immersed in the fluid. The stiffness of half-space Sn can be expressed explicitly from the boundary-condition of solid–fluid interface at the bottom halfspace of [11].

" Sn ¼

0 0 0 qf kx =kz

# ð24Þ

;

where qf is the density of the fluid. With the known Sn, the recursive algorithm is applied to recursively transfer the interface stiffness by using of Eqs. 18, 20, and 21 until the uppermost layer (i.e. the fluid–solid interface) is reached. Finally, the relationship between the displacement and stress at the surface of the uppermost layer is:

ru1 ¼ S1 uu1 :

ð25Þ

Given the wave potential in the fluid as described in Eq. (22), the state vector at Z = 0 can be derived as:

3 u0 6 uz 7 6 kz ð1  RÞ 7 7 ikx x 6 1 7 6 kx 7e : 6 x 7¼6 5 4 r1 5 4 0 z qf ð1 þ RÞ r1 2

ux1

3

2

ð26Þ

Combining Eqs. (25) and (26) and eliminating u0, the acoustic reflection coefficients R can be obtained:

R¼

22 12 21 11 kz ðS11 1 S1  S1 S1 Þ  qS1 kx 22 12 21 11 kz ðS11 1 S1  S1 S1 Þ þ qS1 kx

;

ð27Þ

where Sij1 is the element of the matrix S1, i, j = 1, 2. 4. Experimental validation To verify the feasibility of the reformulation of the ThomsonHaskell method, two representative multilayer structures were selected. The first specimen is made by depositing 30 lm aluminum layers on both surfaces of a 1.0 mm glass disc. It is a typical pure aluminum–glass–aluminum solid multilayered structure. Another specimen is also a three layered structure, that is two 250 lm

Fig. 4. (a) Measured and (b) calculated reflection coefficients R(h, x)for a sandwiched layer structure of two 250 lm stainless plates filled with 100 lm deionized water.

stainless plates filled with 100 lm water. It is a typical solid–liquid multilayered structure. The measurements of the reflection coefficients of the layered structures based on the V(z, t) inversion technique with self-developed scanning acoustic microscopy are performed [12]. The experimental set-up for measurement of acoustic reflection coefficients is shown in Fig. 2. The first specimen is measured with a point focusing transducer of 50 MHz nominal center frequency, and another is measured with 25 MHz one. The half aperture angles of the two transducers are approximately 8°. Due to the band limit of the transducer, the results for very low frequencies and high frequencies are not evident and the proper frequency range is essential. So the frequency band ranges from 15 to 32 MHz was assigned for the 25 MHz transducer, and 40–55 MHz for the 50 MHz transducer. The measured reflection coefficients are then compared with the theoretical ones calculated with the proposed algorithm. The acoustic properties of the materials for the calculation and measurement are listed in Table 1. The measured reflection coefficients with respect to incident angles h and frequency x for the solid layered structure of aluminum–glass–aluminum and the theoretical calculation using the reformulation of Thomson-Haskell method is given in Fig. 3. Fig. 4 shows the reconstructed reflection coefficients R(h, x) for the sandwiched layer structure of two 250 lm stainless plates filled with 100 lm deionized water. From Figs. 3 and 4, the mode traces of the layered structure show virtually identical, while the contrast in the two-dimensional plane is not as sharp as the theoretical calculation. There are several reasonable explanations responsible for these differences. The first is that the attenuation has not been taken into account when executing theoretical calculations, which leads the sharpness of theoretical mode traces. The second would be the finite V(z, t)V(z) acquisition that degrades the inversion resolution. Thirdly, the practical experimental conditions such as the parallelism, the orthogonality of the specimen surface and beam axis are not as ideal as that in theoretical calculation.

5. Summary A reformulation of the Thomson-Haskell method is presented for calculating the reflection coefficients of layered structure

J. Chen et al. / Ultrasonics 52 (2012) 1019–1023

immersing in the coupling fluid. Instead of directly multiplying the layer propagator matrix, the new method splits the layer propagator matrix and excursively determines the interface stiffness matrix starting from the bottom half-space with known stiffness Sn. A formulation for the reflection coefficients is derived based on the obtained interface stiffness matrix of the top layer. This scheme can be applied to a pure solid layers or layered structures containing both fluid and solid layers. It keeps the simplicity but naturally excludes the exponential growth term and thus can be applied at any frequency range. Its validity and feasibility were experimentally proved by the measurement of the reflection coefficients for layered structures based on the inversion of V(z, t) technique. The result of experiments is consistent with the theoretical simulation. The reformulation of the Thomson-Haskell method offers an efficient and effective solution for calculating the acoustic reflection coefficients of layered structures of arbitrary configurations.

Acknowledgments This work is supported by the National Natural Science Foundation of China project 51175465 and the Zhejiang Provincial Natural Science Foundation of China under Grants No. Z1110393. It was also supported by the Fundamental Research Funds for the Central Universities.

1023

References [1] G.J. Yao, H.R. Zhang, K.X. Wang, A calculation approach for acoustic impedance and time-of-flight of an embedded elastic thin layer by ultrasound, J. Phys. D: Appl. Phys. 39 (2006) 4643–4650. [2] T.Y. Liu, Efficient reformulation of Thomson-Haskell method for computation of surface waves in layered half-space, Bull. Seismol. Soc. Am. 100 (2010) 2310–2316. [3] J.W. Dunkin, Computation of modal solutions in layered, elastic media at high frequencies, Bull. Seismol. Soc. Am. 55 (1965) 335–358. [4] D.H. Pei, J.N. Louie, S.K. Pullammanappallil, Improvements on computation of phase velocities of Rayleigh waves based on the generalized R/T coefficient method, Bull. Seismol. Soc. Am. 98 (2008) 280–287. [5] E.L. Tan, A concise and efficient scattering matrix formalism for stable analysis of elastic wave propagation in multilayered anisotropic solids, Ultrasonics 41 (2003) 229–236. [6] L. Wang, S.I. Rokhlin, Stable reformulation of transfer matrix method for wave propagation in layered anisotropic media, Ultrasonics 39 (2001) 413–424. [7] S.I. Rokhlin, L. Wang, Stable recursive algorithm for elastic wave propagation in layered anisotropic media: stiffness matrix method, J. Acoust. Soc. Am. 112 (2002) 822–834. [8] B. Hosten, M. Castaings, Surface impedance matrices to model the propagation in multilayered media, Ultrasonics 41 (2003) 501–507. [9] R. Croce, P. Calmon, L. Paradis, Modeling of propagation and echo formation in a multilayered structure, Ultrasonics 38 (2000) 537–541. [10] Z.Q. Gu, J.D. Achenbach, A. Madan, K. Martin, M.E. Graham, Modeling and acoustic microscopy measurements for evaluation of the adhesion between a film and a substrate, Thin Solid Films 394 (2001) 188–200. [11] T. Kundu, A.K. Mal, R.D. Weglein, Calculation of the acoustic material signature of a layered solid, J. Acoust. Soc. Am. 77 (1985) 353–361. [12] J. Chen, X.L. Bai, K.J. Yang, B.-F. Ju, Angular measurement of acoustic reflection coefficients by the inversion of V(z, t) data with high frequency time-resolved acoustic microscopy, Rev. Sci. Instrum. 83 (2012) 014901.