Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 3, pp. 741-750, Warsaw 2013 STRESS DISTRIBUTION IN TWO-LAYERED HALF-SPACE WITH PERIODICAL STRUCTURE CAUSED BY HERTZ PRESSURE Waldemar Kołodziejczyk, Roman Kulchytsky-Zhyhailo Bialystok University of Technology, Faculty of Mechanical Engineering, Białystok, Poland e-mail: waldekk@pb.edu.pl; r.kulczycki@pb.edu.pl The paper presents analytical solutions for displacements and stresses in an elastic layered half-space with periodic structure in the case of an axi-symmetrical contact problem. The solutions are arrived at employing the Hankel transform. Finally, a comparison is made between the results obtained using the classical theory of elasticity and those obtained within the framework of the homogenizedmodel with microlocal parameters. Key words: layered structures, contact, stress distribution 1. Introduction Extensive applications of layered composite materials in both civil andmechanical engineering, geophysics as well as their widespread presence in nature (varved clays, sandstone-slates, thin layered limestones) havemade compositematerials an object of intensive research. To be able to create or verify strength hypotheses of the materials, adequate knowledge of stress distribution in the composites caused bymechanical factors, eg. pressure distribution, is essential. The formulation of contact problems for periodically layered composites based on the classi- cal theory of elasticity leads to boundary value problems for partial differential equations with discontinuous, strongly oscillating coefficients. As this approach is rather complicated, the clas- sical formulation of the problem is often replaced with approximated models (Sun et al., 1967; Achenbach, 1975; Bensoussan et al., 1978; Christensen, 1980; Sanchez-Palenzia, 1980) in which mechanical properties of themediumare determined on the basis ofmechanical and geometrical properties of the analyzed composite components. In some problems, a periodically layered structure is replaced with a homogeneous trans- versely isotropic medium (Postma, 1955). In this model, calculated strains and stresses are averaged in the fundamental layer. On the other hand, the homogenized model with microlocal parameters (Matysiak and Woźniak, 1987; Woźniak, 1987) makes it possible to evaluate local values of strains and stresses in each sublayer. Owing to the simplification of the mathematical solution, the homogenized model with microlocal parameters has been applied to resolve some contact problems of mechanics and geomechanics (Kaczyński and Matysiak, 1988, 1993, 2001; Kulchytsky-Zhyhailo andMatysiak, 1995; Matysiak and Pauk, 1995). However, the application of thehomogenizedmodel to a periodically layered structure results in errors.Toestimate these errors, it is necessary to compare solutions for theperiodically layered structure based on two approaches: (i) homogenized model with microlocal parameters (ii) non-homogeneous mediumwith separate elastic layers. In the case of heat transfer and theory of elasticity, such investigations were carried out in Kulchytsky-Zhyhailo and Kołodziejczyk (2004, 2007), Kulchytsky-Zhyhailo and Matysiak (2005a,b). It was shown that the solutions obtained using these models were consistent if the 742 W. Kołodziejczyk, R. Kulchytsky-Zhyhailo ratio δ of thickness of the fundamental layer to the size of heating area or loading areawas small (δ < 0.2). One of the fundamental issues of contact mechanics is the problem of pressing a rigid sphere into an elastic half-space. In the case of a homogeneous half-space, the problem is reduced to the loading of the surface by Hertz pressure p(r)= 3 2 p0 √ 1−r2 (1.1) where: p0 –mean contact pressure, r –dimensionless cylindrical coordinate (r < 1)with respect to the radius a of the contact area. Similar results were obtained for a periodically layered half-space using the homogeni- zed model with microlocal parameters (Kulchytsky-Zhyhailo and Matysiak, 1995). For a non- -homogeneous half-space, expression (1.1) is a good description of pressure distribution for the case when δ < 0.2 (Kulchytsky-Zhyhailo andKołodziejczyk, 2007). It should be noted here that the authors focused their attention on two parameters only, i.e. pressure distribution and size of the contact area. In the presentpaper, the authors focus on stress analysis in a layeredmediumwithperiodical structure caused by pressure (1.1) and compare the results obtained using the two models. 2. Formulation of the problem Letus consider a two-layered elastic half-spacewithaperiodic structure loadedbyHertz pressure (Fig. 1). A repeated fundamental layer is composed of two homogeneous elastic sublayers with thickness l1 and l2.Mechanical properties of the sublayers are characterized byYoung’smoduli E1 and E2 and Poisson’s ratios ν1 and ν2. Fig. 1. Diagram of the analyzed problem 3. Homogenized model Initially, the basis for consideration is the homogenized model withmicrolocal parameters. The governing equations of the homogenized model given in terms of macro-displacements can be written in the following form (Kaczyński andMatysiak, 2001): — equations A1 (∂2u ∂r2 + 1 r ∂u ∂r − u r2 ) +(A3+A5) ∂2w ∂r∂z +A5 ∂2u ∂z2 =0 (A3+A5) ( ∂2u ∂r∂z + 1 r ∂u ∂z ) +A4 ∂2w ∂z2 +A5 (∂2w ∂r2 + 1 r ∂w ∂r ) =0 (3.1) Stress distribution in two-layered half-space ... 743 —boundary conditions on the surface of the half-space σzz =−p(r)H(1−r) σrz =0 z=0 (3.2) — regularity conditions in infinity u,w→ 0 r2+z2 →∞ (3.3) where A1 = ( λ̃+2µ̃− [λ]2 λ̂+2µ̂ ) A3 = λ̃− [λ]([λ]+2[µ]) λ̂+2µ̂ A5 = µ̃− [µ]2 µ̂ A4 =(λ̃+2µ̃)− ([λ]+2[µ])2 λ̂+2µ̂ λ̃= ηλ1+(1−η)λ2 [λ] = η(λ1−λ2) λ̂= ηλ1+ η2 1−η λ2 µ̃= ηµ1+(1−η)µ2 [µ] = η(µ1−µ2) µ̂= ηµ1+ η2 1−η µ2 λi = Eiνi (1+νi)(1−2νi) µi = Ei 2(1+νi) and η = l1/l, u, w are dimensionless displacements, σ – stress tensor, H(·) – Heviside’s unit step function. Relations between the stresses andmacro-displacements are as follows (Kaczyński, 1994) σzz =A3 (∂u ∂r + u r ) +A4 ∂w ∂z σ (i) rr =Ki ∂u ∂r +Li u r +Mi ∂w ∂z σ (i) ψψ =Li ∂u ∂r +Ki u r +Mi ∂w ∂z σrz =A5 (∂u ∂z + ∂w ∂r ) (3.4) and Ki =λi +2µi −hiλi [λ] λ̂+2µ̂ Li =λi −hiλi [λ] λ̂+2µ̂ Mi =λi −hiλi [λ]+2[µ] λ̂+2µ̂ i=1,2 h1 =1 h2 =− η 1−η (3.5) where i is the sublayer number in the fundamental layer. From relations (3.4), it follows that the presented model yields separate equations for computation of stress tensor components in the first and second sublayer of the fundamental layer. However, for σzz and σrz they are identical (index i can be omitted). For the stress components σrr and σψψ having jump discontinuities on boundaries between the layers we have two formulas for each. Solutions of the set of equations (3.1) can be expressed as follows (Elliot, 1949a,b) u= ∂Ψ1 ∂r + ∂Ψ2 ∂r w=κ1 ∂Ψ1 ∂z +κ2 ∂Ψ2 ∂z (3.6) where κk = γ2kA1−A5 A3+A5 (3.7) Ψk are elastic potentials satisfying equations ∂2Ψi ∂r2 + 1 r ∂Ψi ∂r + ∂2Ψi ∂z2 γ2i =0 i=1,2 (3.8) 744 W. Kołodziejczyk, R. Kulchytsky-Zhyhailo γ2i are roots of the characteristic equation A1A5γ 4+(A23+2A3A5−A1A4)γ 2+A4A5 =0 (3.9) In order to find a solution to equations (3.8), Hankel’s transform is used f̃(s,z)= ∞∫ 0 rf(r,z)J0(sr) dr (3.10) where J0 is the Bessel function. A general solution to the set of equations (3.8) satisfying conditions (3.3) takes the form Ψ̃i(s,z)=Ci(s)exp ( − sz γi ) i=1,2 (3.11) Expressing (3.4)-(3.6) in Hankel’s transform domain and satisfying boundary conditions (3.2), we obtain C1 =C (p) −1p̃(s) 1 s2 C2 =C (p) 0 p̃(s) 1 s2 (3.12) where p̃(s)= 3 2 p0 √ π 2 J3/2(s) √ s3 C (p) −1 = γ1 A5(γ2−γ1) A3+A5 A3+A1γ 2 1 C (p) 0 =− γ2 A5(γ2−γ1) A3+A5 A3+A1γ 2 2 Taking into account (3.4), (3.6), (3.11), we obtain u(r,z) = ∞∫ 0 Ur(s,z)p̃(s)J1(sr) ds w(r,z) = ∞∫ 0 Uz(s,z)p̃(s)J0(sr) ds σ(i)rr = ∞∫ 0 S (i) r1 (s,z)p̃(s)J0(sr)s ds− 1 r ∞∫ 0 S (i) r2 (s,z)p̃(s)J1(sr) ds i=1,2 σ (i) ψψ = ∞∫ 0 (S (i) r1 −S (i) r2 )p̃(s)J0(sr)s ds+ 1 r ∞∫ 0 S (i) r2 (s,z)p̃(s)J1(sr) ds i=1,2 σzz = ∞∫ 0 Sz1(s,z)p̃(s)J0(sr)s ds σrz = ∞∫ 0 Sz2(s,z)p̃(s)J1(sr) ds (3.13) where Ur(s,z) =− 2∑ k=1 C (p) k−2exp ( − sz γk ) Uz(s,z)=− 2∑ k=1 κkγ −1 k C (p) k−2exp ( − sz γk ) S (i) r1 (s,z) = 2∑ k=1 ( Mi κk γ2k −Ki ) C (p) k−2exp ( − sz γk ) S (i) r2 (s,z) = 2µi 2∑ k=1 C (p) k−2exp ( − sz γk ) Sz1(s,z)= 2∑ k=1 [ A4 κk γ2k −A3 ] C (p) k−2exp ( − sz γk ) Sz2(s,z)=A5 2∑ k=1 (1+κk)γ −1 k C (p) k−2exp ( − sz γk ) Stress distribution in two-layered half-space ... 745 4. Classical model Let us consider a layered half-space as n fundamental layers (2n sublayers) joinedwith a homo- geneous part of the half-space (Fig. 2). It was shown (Kulchytsky-Zhyhailo and Kołodziejczyk, 2007) that over a fixed number of separately considered layers, the mechanical properties of the homogeneous part of the half-space can be quite arbitrarily chosen. They may include the properties of either the first or the second layer. However, considering the fact that the homoge- neous part of the half-space possesses averaged properties of the homogenizedmediumdescribed above, it is possible to obtain a sufficiently exact approximation solution of the problem for a significantly smaller number of separately considered layers. Fig. 2. Diagram of the analyzed problem The equations for separately considered layers can be written as follows ∆uk − uk r2 +dk ∂θk ∂r =0 ∆wk +dk ∂θk ∂z′ =0 k=1,2, . . . ,2n (4.1) where θk = ∂uk/∂r+uk/r+∂wk/∂z ′, k – layer number, d2j =1/(1−2ν1), d2j−1 =1/(1−2ν2), j=1,2, . . . ,n. Differential equations (3.1) for the homogenizedmedium“0” shouldbe attached to the above system of equations. The solution to the equations should satisfy the following conditions: — boundary conditions on the surface z′ =nδ σ(n)zz =−p(r)H(1−r) σ (n) rz =0 z ′ =nδ (4.2) —boundary conditions between the layers andbetween the layers and the homogenizedmedium uk =uk+1 wk =wk+1 σ (k) rz =σ (k+1) rz σ (k) zz =σ (k+1) zz z ′ = zk k=0,1, . . . ,2n−1 (4.3) — regularity conditions in infinity uk,wk → 0 for r 2+z′2 →∞ (4.4) where zk is the coordinate z ′ of the point of intersection of the upper plane limiting the k-th layer with z′ axis. The solution to the classical elasticity theory problem formulated in this way was obtained using the Hankel transform. Applying the algorithm described in Kulchytsky-Zhyhailo andKo- łodziejczyk (2007), we obtain expressions to calculate displacements and stresses in form (3.13), 746 W. Kołodziejczyk, R. Kulchytsky-Zhyhailo where functions Ur,Uz, Sz1, Sz2, Sr1, Sr2 should be replaced with the following 2U (k) z′ (s,z ′)=DkskC (p) 4k−3(s)+DkckC (p) 4k−2(s)+2sskC (p) 4k−1(s)+2sckC (p) 4k (s) 2U(k)r (s,z ′)= µ1 µk S (k) r2 (s,z ′)= ( (2+dk)sk +Dkck ) C (p) 4k−3(s) + ( (2+dk)ck +Dksk ) C (p) 4k−2(s)+2sckC (p) 4k−1(s)+2sskC (p) 4k (s) µ1 µk S (k) r1 (s,z ′)= ( (1+2dk)sk +Dkck ) C (p) 4k−3(s)+ ( (1+2dk)ck +Dksk ) C (p) 4k−2(s) +2sckC (p) 4k−1(s)+2sskC (p) 4k (s) µ1 µk S (k) z′1(s,z ′)= (sk +Dkck)C (p) 4k−3(s)+(ck +Dksk)C (p) 4k−2(s) +2sckC (p) 4k−1(s)+2sskC (p) 4k (s) µ1 µk S (k) z′2(s,z ′)= ( (1+dk)ck +Dksk ) C (p) 4k−3(s)+ ( (1+dk)sk +Dkck ) C (p) 4k−2(s) +2sskC (p) 4k−1(s)+2sckC (p) 4k (s) (4.5) where ck =cosh(s(zk −z ′)) sk =sinh(s(zk −z ′)) Dk = dks(zk −z ′) µk = Ek 2(1+νk) Expressions (4.5) andthecorrespondingequations for thehomogenizedhalf-space“0” include (8n+2) unknown functions Ck(s), k = −1,0,1, . . . ,8n. Satisfying boundary conditions (4.2) and (4.3), we obtain a system of algebraic equations allowing us to determine the unknown functions (Kulchytsky-Zhyhailo and Kołodziejczyk, 2007). 5. Numerical results and discussion Thestress distributionarising in thenon-homogeneoushalf-space obtainedwithin the framework of the homogenized model depends on four dimensionless parameters: ratio E1/E2 of Young’s moduli of the materials of individual layers, Poisson’s ratios ν1 and ν2, the thickness ratio of the upper layer to the thickness of the periodicity cell η = l1/l = δ1/δ. When considering the layered medium, the stress distribution is additionally dependent on the parameter δ. In order to confine the range of the considered parameters, it was assumed that ν1 = ν2 = 0.3. To emphasize the differences that arise between the solutions obtained from bothmodels the ratio ofYoung’smoduluswas assumed E1/E2 =8. Since the results change only insignificantlywithin the considered range of δ parameter when altering the sequence of layers in the periodicity cell, we restrict our analysis to the case when E1/E2 > 1. The calculations show that the stress distributions σzz and σrz for the homogenized model are very close to the corresponding stress distribution for the layered medium. All the distri- butions are similar to those obtained for a homogeneous medium in the classical Hertz contact problem. However, significantly different results can be observed for the stress distribution of σrr and σψψ. These stress distributions along z-axis are shown in Fig. 3 (at z-axis σrr =σψψ). For a non-homogeneous medium on each boundary between layers these stress components un- dergo a jump. This discontinuity is illustrated by black vertical lines. The two gray lines are obtained on the basis of the homogenizedmodel, and they describe the stress distribution in the layers with smaller and greater Young’s modulus respectively. As seen in Fig. 3, the stress di- stribution obtained from the homogenizedmodel describes the stress distribution in the layered mediumwith sufficient accuracy. The difference between the results decreases with reduction of δ parameter. Thus, a suggestion can bemade that a detailed analysis of stress distribution in a Stress distribution in two-layered half-space ... 747 non-homogeneousmedium can be based on the homogenizedmodel, which will radically reduce the time of computation. Fig. 3. Stress distribution σ rr /p max along z-axis for different values of δ; black line – solution for the layeredmedium, gray lines – solution for the homogenizedmodel, dashed lines – boundaries between layers; E1/E2 =8, ν1 = ν2 =0.3, η=0.5 When analyzing the stress level, special attention should be paid to the existence of regions of tensile stresses that are determined on the basis of principal stress σ1 as well as to the distribution of the second invariant of the deviatoric stress tensor J2 = 1 √ 6 √ (σ1−σ2)2+(σ1−σ3)2+(σ2−σ3)2 (5.1) Typical distributions of dimensionless parameters σ1 and J2 with respect to parameter pmax =(3/2)p0 are shown in Figs. 4 and 5. In layers with greater Young’s modulus two regions of tensile stresses occur. The first is found in the vicinity of the points of the unloaded surface of the half-space (Fig. 4a). Maximum of the σ1 stress, similarly to the contact problem for the homogeneousmedium,occurs on theboundaryof the loaded circle.The stress canbedetermined from the expression σ (1) 1 (1,0)= p0µ1 √ A1A4+A3 (5.2) The other region arises at a certain depth. Tensile stresses in this region are generally much smaller than the stress at points z=0, r=1. Fig. 4. Tensile stress distribution σ1 in layers with: (a) greater Young’s modulus, (b) smaller Young’s modulus In layers with smaller Young’s modulus, the region of tensile stresses occurs beneath the unloaded surface of the half-space (Fig. 4b). In this case, there exist two local maxima of σ1 748 W. Kołodziejczyk, R. Kulchytsky-Zhyhailo stresses. The first at the boundary of the loading area, the other at a certain small depth in the vicinity of the boundary of the contact surface. It should be noted that for points z=0, r=1 we have a simple relation σ (1) 1 /σ (2) 1 =µ1/µ2. The distribution of J2 parameter in the layers with greater Young’s modulus has three local maxima (Fig. 5a). Thefirst one occurs at somedepth on z-axis, the second one – at point z=0, r =0, the third one – at points z =0, r =1. The distribution of J2 parameter in layers with smaller Young’s modulus (Fig. 5b) is similar the one that occurs for the homogeneous isotropic medium in the classical Hertz contact problem. The only maximum arises on z-axis at some depth beneath the surface. Fig. 5. Distribution of J2 in layers with: (a) greater Young’s modulus, (b) smaller Young’s modulus Figures 6 and 7 show maximum values of relative parameters σ1/σ (i) 1,H and J2/J (i) 2,H as a function of η for three values E1/E2 = 2,4,8 described by lines 1,2,3 respectively. σ (i) 1,H, J (i) 2,H are corresponding parameters given by the classical Hertz contact problem for the homogeneous mediumwithmechanical properties of the i-th layer (i=1,2) (Timoshenko andGoodier, 1951; Johnson, 1985) σ (i) 1,H = 1 3 p (i) max,H(1−2νi) J (i) 2,H(νi =0.3)= 0.356p (i) max,H (5.3) where pmax,H is the maximum pressure in the classical Hertz contact problem of pressing of a sphere into a homogeneous isotropic half-space and p (i) max,H pmax = [ A1(γ1+γ2)µi (A1A4−A 2 3)(1−νi) ]2/3 (5.4) Black lines in Figs. 6a and 6b show values of the parameter on the boundary of load circle. Gray lines showmaximum values on z-axis. It can be seen in Fig. 6a that tensile stresses on z-axis and on the boundary of the circle of contact achieve comparable values for small values of η and relatively large difference in values of Young’s moduli. For η close to 1, the tensile stress does not appear at all. It follows from Fig. 6b that the local maximum J2 existing on z-axis dominates the local maximumat theboundaryof the loading circle.One exception is the small range 0.37<η< 0.47 (E1/E2 = 8). It should be noted that gray line 3 consists of two parts. The left hand part corresponds to the case when themaximum J2 at the center of the loading circle prevails in the maximum J2 on z-axis at some depth. The part on the right refers to the opposite situation. When η → 1, the maximum on z-axis is at the same depth as in the classical Hertz contact problem, i.e. z=0.48 (ν1 =0.3). With decreasing η, the depth of maximum J2 increases (e.g. E1/E2 =4 from 0.48 (η→ 1) to 0.66 (η=0.01)). Stress distribution in two-layered half-space ... 749 Fig. 6. Distribution of (a) σ1, (b) J2 in layers with greater Young’s modulus as function of η Fig. 7. Distribution of (a) σ1, (b) J2 along z-axis in layers with smaller Young’s modulus as function of η Lines in Figs. 7a and 7b show maximum σ1 and J2 in the layers with smaller Young’s modulus. Curves 2 and 3 in Fig. 7a consist of two parts. The one on the left corresponds to the case when maximum σ1 is on the boundary of the loaded circle, the one on the right – when maximum σ1 is at some depth under the unloaded surface. When η→ 0, the maximum J2 in layers with smaller Young’s modulus occurs at the same depth as in the Hertz contact problem.With increasing η, the depth of maximum J2 decreases (i.e E1/E2 =4 from 0.48 (η→ 0) to 0.25 (η=0.99). 6. Conclusions The paper presents the stress state in a layered half-space with periodical structure caused by Hertz pressure. Two different models of the layered medium have been considered: one – a homogenized model in which mechanical properties of the medium are determined using both themechanical and geometrical properties of the analyzed composite components, and the other – where a certain number of layers are considered as separate homogeneous and interacting entities satisfying the conditions of ideal mechanical contact. In conclusion, it is demonstrated that stresses in layers of different mechanical properties are quite different. Additionally, the location of tensile stress regions σ1 as well as the distribution of the second invariant of the deviatoric stress tensor J2 are presented. Acknowledgements The investigations described in this paper are a part of the research projects S/WM/2/08 and S/WM/1/08 realized at Bialystok University of Technology. 750 W. Kołodziejczyk, R. Kulchytsky-Zhyhailo References 1. Achenbach J.D., 1975, A Theory of Elasticity with Microstructure for Directionally Reinforced Composites, CISMCourses and Lectures, Springer, NewYork 2. Bensoussan A., Lions J.L., Papanicolaou G., 1978,Asymptotic Analysis for Periodic Struc- tures, North-Holland, Amsterdam 3. Christensen R.M., 1980,Mechanics of Composite Materials, J.Wiley and Sons 4. Elliot D.A., 1949a, Three-dimensional stress distributions in aeolotropic hexagonal crystals, Cambridge Phil. Society, 44, 522-533 5. ElliotD.A., 1949b,Axial symmetric stress distributions in aeolotropic hexagonal crystals,Cam- bridge Phil. Society, 45, 621-630 6. Johnson K.L., 1985,Contact Mechanics, Cambridge Univesity Press, Cambridge 7. Kaczyński A., 1994, Three-dimensional thermoelastic problems of interface crack in periodic two-layered composites,Engineering Fracture Mechanics, 48, 783-800 8. Kaczyński A., Matysiak S.J., 1988, Plane contact problems for periodic two-layered elastic composite, Ingenieur-Archiv, 58, 137-147 9. Kaczyński A., Matysiak S.J., 1993, Rigid sliding punch on a periodic two-layered elastic half- -space, Journal of Theoretical and Applied Mechanics, 31, 2, 295-306 10. Kaczyński A., Matysiak S.J., 2001, On 3D punch problems for a periodic two-layered elastic half-space, Journal of Theoretical and Applied Mechanics, 39, 3, 523-538 11. Kulchytsky-ZhyhailoR., KołodziejczykW., 2004,A bi-dimensional contact problemof the theory of elasticity relating to indentation of the die into an elastic laminated semispace,Friction and Wear, 25, 2, 125-133 12. Kulchytsky-Zhyhailo R., Kołodziejczyk W., 2007, On axisymmetrical contact problem of pressure of a rigid sphere into a periodically two-layered semi-space, International Journal of Mechanical Sciences, 49, 704-711 13. Kulchytsky-Zhyhailo R., Matysiak S.J., 1995, On three-dimensional problems of multilay- ered periodic elastic composites, Journal of Theoretical and Applied Mechanics, 33, 4, 771-781 14. Kulchytsky-Zhyhailo R., Matysiak S.J., 2005a, On heat conduction problem in a semi- -infinite periodically laminated laminated layer, International Communications in Heat and Mass Transfer, 32, 1/2, 123-132 15. Kulchytsky-Zhyhailo R., Matysiak S.J., 2005b, On some heat conduction problem in a periodically two-layered body, International Communications in Heat and Mass Transfer, 32, 332-340 16. Matysiak S.J., Pauk V.J., 1995, Plane contact problem for periodic laminated composite invo- lving frictional heating,Archive of Applied Mechanics, 66, 1/2, 82-89 17. Matysiak S.J., Woźniak C., 1987, Micromorphic effects in modelling of periodic multilayered elastic composites, International Journal of Engineering Science, 25, 549-559 18. Postma G.W., 1995,Wave propagation in a stratifiedmedium, Geophysics, 20, 780-806 19. Sanchez-Palenzia E., 1980, Non-homogeneous media and vibration theory, Lecture Notes in Physics, 27, Springer Verlag 20. Sun C.T., Achenbach J.D., Herrmann G., 1967, Continuum theory for a laminatedmedium, Trans. ASME, Journal of Applied Mechanics, 467-475 21. Timoshenko S., Goodier J.N., 1951,Theory of Elasticity, McGraw-Hill, NewYork 22. Woźniak C., 1987, A nonstandardmethod of modelling thermoelastic composites, International Journal of Engineering Science, 25, 483-499 Manuscript received September 27, 2012; accepted for print December 28, 2012