1 Introduction

Functionally graded materials (FGMs), as advanced composite materials, have extensively been used in recent years. The concept of the FGMs originated in Japan during the space plane project in 1984 (Miyamoto et al. 1999). An FGM is a mixture of two or more constituent phases with a variable composition. Microstructures of the FGMs are often fabricated in three major types: continuously graded microstructures, discretely graded microstructures, and multi-phase graded microstructures (Suresh and Mortensen 1998). The simplest FGMs are composed of two constituent materials whose volume fractions change gradually so that sum of their volume fractions equals unity in each arbitrary point of the component. Unlike the traditional composites, which are piecewise-homogeneous mixtures or layered systems, in majority of the FGMs the mechanical properties vary continuously. An FGM can be customized to an application by specifying the form of the gradation of the material properties to fit the design requirements.

Owing to their special features and application potentials, the FGMs have drawn attention of numerous scientists and engineers in broad areas of research. The FGMs may be ideal in applications where the operating conditions are severe. Some typical examples of the FGM components are the heat-engine components and the rocket heat shields. Under these conditions, the FGMs may exhibit time-dependent behaviors (Mukherjee and Paulino 2003; Khazanovich 2008). Utilizing the FGMs’ potentials requires development of more accurate numerical procedures and appropriate modeling techniques. A realistic phenomenological model for describing the time-dependent behavior of the materials is the linear viscoelasticity. Owing to complexities of the viscoelastic models which include the time as an independent variable, the available exact solutions have been obtained for only a few simplified problems (Mukherjee and Paulino 2003; Khazanovich 2008). Rigorous predictions of the viscoelastic behaviors usually rely on numerical approaches such as the finite difference, finite element, boundary element, and meshless methods.

During the recent years, the boundary element method (BEM) has found considerable applications in various engineering fields, such as the contact mechanics, viscoelasticity, thermoelasticity, elastoplasticity, elastodynamics, fracture mechanics, and geomechanics (Paris and Canas 1997; Sutradhar et al. 2008). Many important problems involve structures subjected to both mechanical and thermal loadings. Some researchers have studied thermoelasticity of the solids using the BEM (Becker 1992; Gao 2003; Hematiyan 2007). Moreover, numerous researchers have developed BEMs to model behaviors of the homogeneous viscoelastic media (Schanz 1999; Perez-Gavilan and Aliabadi 2001; Huang et al. 2006; Mesquita and Coda 2007; Ashrafi and Farid 2009, 2010; Cezario et al. 2011; Ashrafi et al. 2012). They have presented BEM solutions for Laplace transforms of the primary variables of the elastic problems and then returned the solutions to the time domain through proper numerical inversions. In the boundary integral equation methods, only the boundary has to be discretized and there is no need to discretize the domain under consideration into elements. Furthermore, application of the BEM is mainly convenient for analysis of majority of the engineering problems, including those with infinite domains. Because only the boundary surfaces have to be discretized and the boundary conditions at the infinity are automatically satisfied by adopting a proper fundamental solution.

The first computational works on constitutive modeling of the thermoviscoelastic behaviors of the FGMs were accomplished using the finite element method (FEM), due to its capability to trace variations of the material properties in a simpler manner (Mahmoud et al. 2011). Dave et al. (2011) proposed a functionally graded viscoelastic finite element formulation based on the correspondence principle, using a special graded element. Singh et al. (2011) focused on numerical simulation of the inhomogeneities/discontinuities (cracks, holes and inclusions) in the FGMs using an extended FEM. Moreover, Sladek et al. (2006) developed a meshless local Petrov–Galerkin method for analysis of the continuous heterogeneous and linear viscoelastic materials. They derived local integral equations for the subdomains in the Laplace-transformed solution domain. The time-dependent behavior of the homogeneous thermoviscoelastic materials can accurately be predicted by the BEM. Some researchers introduced BEM solutions in incremental forms in the time domain. Although several works have been developed to effectively model the constitutive behavior of the FGMs, most of these researches have been limited to the elastic behavior. In contrast, few works can be found in literature on analysis of the continuous heterogeneous and linear viscoelastic media. Application of the BEM to the heterogeneous media is a formidable task, since fundamental solutions correspond to concentrated loads imposed on such media are difficult to obtain. The fundamental solution of the heat transfer problem in the heterogeneous media has been presented by Gray et al. (2003) and Kuo and Chen (2005). The fundamental solutions for 2D and 3D elastostatic problems have recently been developed for layered media and for exponentially graded materials (Chan et al. 2004; Criado et al. 2007, 2008; Ashrafi et al. 2013). Green functions of the exponentially graded media have been derived using Fourier transform techniques. Furthermore, crack analysis of the 2D continuous and linear elastic FGMs has been accomplished by Yue et al. (2003) and Gao et al. (2008).

Present research introduces a direct boundary element formulation and a relevant Green’s function (fundamental solution) for quasistatic problems of thermoviscoelastic FGMs whose material properties vary continuously along a given direction. In this regard, the shear modulus of the FGM is assumed to obey an exponential function along a specified direction and Poisson’s ratio is assumed to be constant. The formulation is developed based on the differential constitutive equations of the heterogeneous Kelvin–Voigt viscoelastic solid model. The proposed solution algorithm is capable of solving quasistatic problems of components made of functionally graded thermoviscoelastic materials with arbitrary boundary conditions. Finally, three applied examples are included to evaluate accuracy and efficiency of the proposed boundary element formulations for the thermoviscoelastic FGMs.

2 Constitutive relations of thermoviscoelastic FGMs

Some important engineering materials simultaneously store and dissipate the mechanical energy when they are subjected to externally applied loads. This type of behavior is called viscoelasticity. Intrinsically, constitutive equations of the viscoelasticity incorporate the elastic deformations, equivalent viscous flows, and the patterns that characterize the combined behavior. The viscoelastic behavior of the materials in the uniaxial stress field may be modeled using elastic and viscous elements. Moreover, flow of the thermal energy within a solid leads to conduction heat transfer which in turn establishes thermal distributions and thermal gradients.

Employing adequate combinations of the two fundamental elements, i.e., the elastic spring and viscous dashpot, it is easy to construct proper viscoelastic models. One simple combination is Kelvin–Voigt viscoelastic model. This model constitutes a simple representation of the viscoelastic behavior, and at the same time, it is suitable for modeling of the heterogeneous materials. The Kelvin–Voigt model for the solids consists of a spring and a dashpot in parallel, as shown in Fig. 1. For this model, we have

$$ \sigma_{ij} = \left( {E + \eta \partial_{t} } \right)\varepsilon_{ij} $$
(1)

where E and η are material constants, as shown in Fig. 1. According to the Kelvin–Voigt solid model, the total stress and strain components may be related to the elastic, viscoelastic, mechanical or thermal components as follows

$$ \sigma_{ij} = \sigma_{ij}^{e} + \sigma_{ij}^{v} $$
(2)
$$ \varepsilon_{ij} = \varepsilon_{ij}^{M} + \varepsilon_{ij}^{T} $$
(3)
$$ \varepsilon_{ij}^{M} = \varepsilon_{ij}^{e} = \varepsilon_{ij}^{v} $$
(4)

where

$$ \varepsilon_{ij}^{T} = \alpha (x)\,\delta_{ij} (T - T_{0} ) $$
(5)

in which the superscripts v, e, M, and T denote the viscous, elastic, mechanical, and thermal components, respectively. α is the coefficient of thermal expansion and T 0 is the reference temperature (Boley and Weiner 1985).

Fig. 1
figure 1

A sketch of the heterogeneous Kelvin–Voigt viscoelastic Model

Full size image

Generally, for infinitesimal deformations, the stress–strain relations of a linear FGM may be written as

$$ \sigma_{ij}^{e} = C_{ijkm} (x)\varepsilon_{km} - \kappa_{ij} (x)(T - T_{0} ) $$
(6)

where C ijkm and κ ij are the heterogeneous elastic and thermoelastic coefficients tensors of the FGM, respectively and are functions of the spatial variable x. It may be easily verified that the most general fourth–order isotropic C ijkm and κ ij tensors have the following form

$$ C_{ijkm} \left( x \right) = \lambda \left( x \right)\delta_{ij} \delta_{km} + \mu \left( x \right)\left( {\delta_{ik} \delta_{jm} + \delta_{im} \delta_{jk} } \right) $$
(7)
$$ \kappa_{ij} (x) = [3\lambda (x) + 2\mu (x)]\alpha (x)\delta_{ij} $$
(8)

where λ and μ are Lame’s parameters (Mase and Mase 1999) and δ is Kronecker’s delta symbol. If exponential gradations are chosen for Lame’s parameters in terms of the x-direction, one may write

$$ \lambda \left( x \right) = \lambda_{0} \;\exp \left( {2\gamma x} \right) $$
(9)
$$ \mu \left( x \right) = \mu_{0} \;\exp \left( {2\gamma x} \right) $$
(10)

where γ is the material gradation exponent. Similarly, for the viscous stress components, the following constitutive equation may be used

$$ \sigma_{ij}^{v} = K_{ijmn} \left( x \right)\dot{\varepsilon }_{mn} $$
(11)

in which K ijmn represents the heterogeneous isotropic viscosity coefficients tensor of the material which is a function of the spatial variable x, and can be defined as

$$ K_{ijmn} \left( x \right) = \left[ {\beta_{\lambda } \lambda \left( x \right)} \right]\delta_{ij} \delta_{mn} + \left[ {\beta_{\mu } \mu \left( x \right)} \right]\left( {\delta_{im} \delta_{jn} + \delta_{in} \delta_{jm} } \right) $$
(12)

where β λ and β μ are the hydrostatic and deviatoric viscosity coefficients, respectively (Mase and Mase 1999).

Substituting Eqs. 6 and 11 into Eq. 2, leads to

$$ \sigma_{ij} = C_{ijkl} (x)\varepsilon_{kl} - \kappa_{ij} (x)(T - T_{0} ) + K_{ijkl} (x)\dot{\varepsilon }_{kl} $$
(13)

In the numerical solution section, a further simplification is employed for the viscosity coefficients, i.e. β λ  = β μ  = β, in order to obtain the boundary integral equations. In this case, Eq. 13 reduces to

$$ \sigma_{ij} = C_{ijkl} (x)\;(\varepsilon_{kl} + \beta \,\dot{\varepsilon }_{kl} ) - \kappa_{ij} (x)(T - T_{0} ) $$
(14)

It may be reminded that the viscosity coefficients can be identified from results of the creep shear and uniaxial tensile tests (Flügge 1975).

The governing equation of the transient heat conduction, in the absence of the heat generation, can be expressed as

$$ \nabla^{2} T(x,t) - \frac{\rho c}{\vartheta }\frac{\partial T(x,t)}{\partial t} = 0 $$
(15)

where \( \vartheta \), c, and ρ are the thermal conductivity, specific heat capacity at a constant volume, and mass density, respectively.

3 BEM form of the governing equations

The BEM is based on the boundary integral equations. There are several methods for deriving the BEM formulations: the reciprocal theorem, the weighted residual concept, and the variational approaches (Sutradhar et al. 2008). In the present research, the boundary integral equations are derived using the weighted residual concept and extended to the heterogeneous Kelvin–Voigt solid model.

The quasistatic viscoelastic integral equations for the boundary and interior points may be obtained imposing the weighted residual technique on the equilibrium equations in terms of the stress components

$$ \sigma_{ij,i} = B_{j} = 0 $$
(16)

where B j is the body force acting in j direction. Following Galerkin’s technique, the fundamental solution of an elastic infinite body may be adopted as a weight function (Paris and Canas 1997). Therefore, Eq. 16 can be weighted over the considered domain D as

$$ 0 = \int\limits_{D} {\psi_{i} } \left\langle {\sigma_{ij,j} + B_{i} } \right\rangle {\text{dv}} $$
(17)

where ψ i is the heterogeneous displacement fundamental solution. Determination of the fundamental solution of the original partial differential equation is a vital issue for reducing a boundary value problem to a boundary integral equation. After the fundamental solution is used in conjunction with the corresponding Green’s function, one may reduce the problem to a BEM one. This solution represents a response to a unit concentrated load applied at a point of the infinite domain. The fundamental solution of the 2D elastostatic components made of exponentially graded materials has recently been developed as

$$ \psi = e^{ - 2\gamma x} \left\{ {C_{0} K_{0} \left( {\left| \gamma \right|\left| r \right|} \right) + C_{1} K_{1} \left( {\left| \gamma \right|\left| r \right|} \right) + \psi^{NS} } \right\} $$
(18)

where

$$ C_{0} = \frac{3 - 4\upsilon }{{8\pi \mu_{0} (1 - \upsilon )}}\left[ {\begin{array}{*{20}c} 1 & 0 \\ 0 & 1 \\ \end{array} } \right]\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,C_{1} = \frac{{\frac{{r^{2} }}{\left| r \right|}\left| \gamma \right| - \frac{{\gamma^{2} }}{\left| \gamma \right|}\left| r \right|}}{{16\pi \mu_{0} (1 - \upsilon )}}\left[ {\begin{array}{*{20}c} 1 & 0 \\ 0 & 1 \\ \end{array} } \right] $$
(19)
$$ \psi^{NS} = \frac{1}{\pi }\int\limits_{ - \infty }^{\infty } {} \int\limits_{ - \infty }^{\infty } {} U\left( {\xi_{1} ,\xi_{2} } \right)e^{ - i\;\xi \cdot r} \;d\xi_{1} \;d\xi_{2} $$
(20)

and K 0 and K 1 are the modified Bessel functions (Chan et al. 2004).

Integrating Eq. 17 by parts and applying the divergence theorem yields

$$ 0 = \int\limits_{\partial D}^{{}} {} \psi_{i} \sigma_{ij} n_{j} ds - \int\limits_{D}^{{}} {} \psi_{i, j} \sigma_{ij} dv + \int\limits_{D}^{{}} {} \psi_{i} B_{i} dv $$
(21)

where ∂D is the boundary of the FGM component, and n j is component of the outward normal vector of the boundary. Knowing that,

$$ \varepsilon_{ij}^{\psi } \sigma_{ij} = 0.5\left( {\psi_{i,j} + \psi_{j,i} } \right)\sigma_{ij} = \psi_{i,j} \sigma_{ij} $$
(22)

where \( \varepsilon^{\psi } \) is the fundamental strain term, Eq. 21 reduces to

$$ 0 = \int\limits_{\partial D}^{{}} {} \psi_{i} t_{i} ds - \int\limits_{D}^{{}} {} \varepsilon_{ij}^{\psi } \sigma_{ij} dv + \int\limits_{D}^{{}} {} \psi_{i} B_{i} dv $$
(23)

By substituting the thermoviscoelastic constitutive equation of the FGM (Eq. 14) into the integral Eq. 23, one has

$$ 0 = \int\limits_{\partial D} {\psi_{i} \;t_{i} \;ds} - \int\limits_{D} {\varepsilon_{ij}^{\psi } C_{ijkl} (x)\varepsilon_{kl} \;dv - \beta \int\limits_{D} {\varepsilon_{ij}^{\psi } C_{ijkl} (x)\varepsilon_{kl} \;dv} } \; + \int\limits_{D} {\varepsilon_{ij}^{\psi } \kappa_{ij} (x)(T - T_{0} )} \;dv + \int\limits_{D} {\psi_{i} \;B_{i} } \;dv $$
(24)

On the other hand, noting that,

$$ \varepsilon_{ij}^{\psi } \;C_{ijkl} (x)\;\varepsilon_{kl} = \sigma_{kl}^{\psi } \,\varepsilon_{kl} = \sigma_{kl}^{\psi } \,u_{k,l}\, {\text {and}} \,\varepsilon_{ij}^{\psi } \kappa_{ij} (x) = \sigma_{ij}^{T\psi } $$
(25)

Eq. 24 leads to

$$ 0 = \int\limits_{\partial D} {\psi_{i} \;t_{i} \;ds} - \int\limits_{D} {\sigma_{ij}^{\psi } u_{i,j} \;dv - \beta \int\limits_{D} {\sigma_{ij}^{\psi } \dot{u}_{i,j} \;dv} } + \int\limits_{D} {\sigma_{ij}^{T\psi } (T - T_{0} )} \;dv + \int\limits_{D} {\psi_{i} \;B_{i} } \;dv $$
(26)

By integrating the second and third terms of Eq. 26 by parts, one may obtain

$$ \begin{aligned} 0 & = \int\limits_{\partial D} {\psi_{i} \;t_{i} \;ds} - \int\limits_{\partial D} {\sigma_{ij}^{\psi } u_{i} n_{j} \;ds - \beta \int\limits_{\partial D} {\sigma_{ij}^{\psi } \dot{u}_{i} n_{j} \;ds} } \; + \int\limits_{D} {\sigma_{ij,j}^{\psi } u_{i} \;dv + \beta \int\limits_{D} {\sigma_{ij,j}^{\psi } \dot{u}_{i} \,dv} } \\ & \quad + \int\limits_{\partial D} {P_{j}^{\psi } \,T\,ds} - \int\limits_{\partial D} {Q_{j}^{\psi } } \;q\,ds + \int\limits_{D} {\psi_{i} \;B_{i} } \;dv \\ \end{aligned} $$
(27)

The displacement and traction fields corresponding to the point force solution can be written as

$$ \psi_{i} = \psi_{ij} \left( {x,z} \right)e_{j} t_{i}^{\psi } = t_{ij}^{\psi } \left( {x,z} \right)e_{j} $$
(28)

where e i is the component of the unit vector associated with a unit positive force in the i direction applied at the specified z point. Eq. 27 can be rewritten by using the fundamental equilibrium equation, i.e.

$$ \sigma_{ij,j}^{\psi } = - B_{i}^{\psi } \left( z \right) = - \Updelta \left( {x,z} \right)e_{i} $$
(29)

in which Δ (x, z) is Dirac delta function whose z and x are the field and the source arguments, respectively (Paris and Canas 1997). Also the kernels P ψ and Q ψ are introduced for the temperature and heat flux, respectively (Aliabadi and Wen, 2011). Consequently, Eq. 27 may be rewritten as

$$ \subset_{ij} \,u_{i} (x) + \beta \subset_{ij} \,\dot{u}_{i} (x) = \int\limits_{\partial D} {\psi_{ij} t_{i} ds} - \int\limits_{\partial D} {t_{ij}^{\psi } u_{i} ds} - \beta \int\limits_{\partial D} {t_{ij}^{\psi } \dot{u}_{i} \;ds} + \int\limits_{\partial D} {P_{j}^{\psi } \,T\,ds} - \int\limits_{\partial D} {Q_{j}^{\psi } } \;q\,ds + \int\limits_{D} {\psi_{ij} \;B_{i} } \;dv $$
(30)

where the free term (\( \subset_{ij} \)) may be found in the elastostatic formulations presented in the BEM handbooks (Paris and Canas 1997; Sutradhar et al. 2008). Equation 30 may be considered as an integral constitutive equation for the viscoelastic FGMs and may be employed to model the heterogeneous Kelvin–Voigt solids. The domain integral of the body forces may be easily related to its equivalent boundary integral that may be determined based on the boundary values. The simple and robust method, which is called the radial integration method, is used for transforming the domain integrals into their equivalent boundary integrals (Gao 2002). Therefore, any 2D domain integrals may be evaluated without the need to discretize the domain into internal elements. By neglecting the body force B i , the relevant domain integral may be omitted completely. The boundary integral form of Eq. 15 may similarly be transformed into the boundary through employing the weighted residual method. Thus one has

$$ \subset_{ij} T_{i} (x,t) + \frac{1}{\rho \,c}\int\limits_{0}^{\tau } {\int\limits_{\partial D} {T_{ij}^{\psi } q_{i} \,ds\,dt} } = \frac{1}{\rho \,c}\int\limits_{0}^{\tau } {\int\limits_{\partial D} {q_{ij}^{\psi } T_{i} \,ds\,dt} } + \int\limits_{D} {T_{ij}^{\psi } T_{i} \,dv} $$
(31)

where q is the heat flux. Moreover, the temperature and heat flux fundamental solutions T Ψ and q Ψ may be derived as

$$ T_{ij}^{\psi } = \frac{1}{4\pi \;\alpha (x,t - \tau )}e^{{\left( {\frac{{ - r^{2} }}{4\,\alpha (x,t - \tau )}} \right)}} $$
(32)
$$ q_{ij}^{\psi } = - k\;\frac{{\partial T_{ij}^{\psi } }}{\partial n} $$
(33)

For the integral Eq. 31 to become a purely boundary–type, the point p has to be moved to the boundary, therefore \( \subset_{ij} \) should be 0.5.

To derive the stress integral equations for the interior points, the strain integral equations have to be derived first. At the interior points, the displacement integral equation is given by

$$ u_{i} (x) + \beta \dot{u}_{i} (x) = \int\limits_{\partial D} {\psi_{ij} t_{i} ds} - \int\limits_{\partial D} {t_{ij}^{\psi } u_{i} ds - \beta \int\limits_{\partial D} {t_{ij}^{\psi } \dot{u}_{i} \;ds} } + \int\limits_{\partial D} {P_{j}^{\psi } \,T\,ds} - \int\limits_{\partial D} {Q_{j}^{\psi } } \;q\,ds $$
(34)

Using tensor notations, both the small strains and strain rates may be related to the displacements and displacement rates as (Simo and Hughes 1997)

$$ \varepsilon_{kl} = \frac{1}{2}\left\{ {u_{k,l} + u_{l,k} } \right\} $$
(35)
$$ \dot{\varepsilon }_{kl} = \frac{1}{2}\left\{ {\dot{u}_{k,l} + \dot{u}_{l,k} } \right\} $$
(36)

The stress equation for an internal point may be obtained by using derivatives of the displacement Eq. 34 with respect to the source point location. Therefore, the stress field of the functionally graded viscoelastic component may be obtained combining the constitutive Eq. 14 and the boundary integral Eq. 34, which results in

$$ \sigma_{ij} (x) = \int\limits_{\partial D} {d_{kij}^{\psi } \,t_{k} ds} - \int\limits_{\partial D} {s_{kij}^{\psi } \,u_{k} ds - \beta \int\limits_{\partial D} {s_{kij}^{\psi } \,\dot{u}_{k} \;ds} } - \int\limits_{D} {\varepsilon_{ij}^{\psi } \kappa_{ij} (x)(T - T_{0} )} \;dv+ \int\limits_{\partial D} {P_{ij}^{\psi } \,T\,ds} - \int\limits_{\partial D} {Q_{ij}^{\psi } } \;q\,ds - \frac{E\alpha }{1 - 2\upsilon }\;T\delta_{ij}$$
(37)

where \( P_{ij}^{\psi } \), \( Q_{ij}^{\psi } \), \( d_{ijk}^{\psi } \), and \( s_{ijk}^{\psi } \) are the fundamental kernels which were developed based on the heterogeneous displacement and traction fundamental solutions.

It is only possible to solve the boundary element formulations for very simple problems, analytically. The first step in the numerical discretization is dividing the boundary ∂D into N e elements. After choosing identical number of the source and nodal points, all integrals may be calculated numerically and consequently, the discretized BEM equations may be written in a matrix form.

4 The numerical solution algorithm

To solve the time-dependent differential governing equations of the functionally graded domain (Eqs. 30, 31 and 37), it is necessary to approximate the velocity in the time domain by a time marching treatment (Simo and Hughes 1997). This is carried out as follows

$$ \dot{\mathbf{u}}_{n + 1} = \frac{{\mathbf{u}_{n + 1} - \mathbf{u}_{n} }}{\Updelta t} $$
(38)

In this regard, the time domain is discretized into time steps smaller than the fundamental vibration period of the structure. The solution begins by incorporating the initial conditions. The quantities (u n+1, T n+1, q n+1) obtained by solving the matrix system of equations at the end of each time step i.e. t n+1, may be used as initial conditions for the next time step. The presented algorithm has been implemented in the MATLAB commercial software.

5 Results and discussions

In the present section, three numerical problems are used to investigate quasistatic thermoviscoelastic behavior of FGM structures subjected to mechanical or thermal loads. Verification of the results as well as evaluating the accuracy and efficiency of the proposed formulations is accomplished through the first two examples.

5.1 Example 1—A viscoelastic strip subjected to a uniform loading

As a first numerical example, a heterogeneous viscoelastic strip shown in Fig. 2 is considered. The viscoelastic strip is subjected to a uniform tensile load P H(t) at the end (x = L), wherein H(t) is the Heaviside step function. The opposite end of the strip is fixed in the x–direction, while its other boundaries are free of tractions. The length and width of the strip are L = 3 and b = 1, respectively. In discretization of the boundaries, 40 nodes are used as shown in Fig. 2. The numerical results of the presented BEM formulation for time history of longitudinal displacement u of node 16 are illustrated and compared with analytical results of Flügge (1975) in Fig. 3. As it may be readily deduced from Fig. 3, there is a good agreement between the results.

Fig. 2
figure 2

A functionally graded viscoelastic strip subjected to a uniform tension

Full size image
Fig. 3
figure 3

A comparison between results of the present viscoelastic boundary element method (VBEM) and analytical results of Flügge (1975) for time history of the longitudinal displacement of node 16 of the homogeneous viscoelastic strip

Full size image

In the next stage, the strip is assumed to be constructed from a heterogeneous material with an exponential gradation in the x-direction. The material constants are adopted as \( \lambda_{0} = 1/6 \) \( \mu_{0} = 0.25 \) and β = 10 and the material heterogeneity exponent γ, introduced in Eqs. 12 and 13, is chosen as 0.25 and 0.5. By increasing the heterogeneity exponent γ, Lame’s elastic constants increase, and in turn, the longitudinal displacements of the functionally graded strip decrease. The longitudinal displacement time histories of the heterogeneous viscoelastic strip are illustrated in Fig. 4, for different gradation exponents. To evaluate behavior of the strip after an abrupt unloading, a case wherein the tensile load is suddenly removed at a prescribed time, e.g., at t = 5, is considered. The relevant time history of the longitudinal displacement u of node 16 is depicted in Fig. 5, for the heterogeneity exponent γ = 0.2.

Fig. 4
figure 4

Time histories of the longitudinal displacement of node 16 of the heterogeneous viscoelastic strip, for different gradation exponents

Full size image
Fig. 5
figure 5

Time history of longitudinal displacement of the node 16 of the heterogeneous viscoelastic strip, under a tensile load that is suddenly removed at t = 5 s

Full size image

5.2 Example 2—A viscoelastic cylindrical vessel under internal pressure

In the second numerical example, a thick-walled cylindrical vessel with exponential material gradation in the r-direction under radial internal pressure PH(t) is considered, as shown in Fig. 6. Since all materials exhibit viscoelastic behaviors, investigation of the pressurized viscoelastic cylindrical vessels is of technical importance. The considered vessel is assumed to be sufficiently long such that a plane strain condition assumption holds. In this case, it is necessary to prevent moving of both ends of the cylinder in the axial direction (z-direction) so that the condition (ε z = 0) holds. Let us consider a thick-walled cylinder of inner radius R 1  = 300 (mm) and outer radius R 2  = 500 (mm) that is subjected to an internal pressure P = 2 MPa. Due to symmetry of the structure, only a quarter of the cylinder has to be discretized, as shown in Fig. 6.

Fig. 6
figure 6

A functionally graded cylindrical vessel subjected to a radial internal pressure (P)

Full size image

The numerical solution is obtained using 30 boundary elements for approximation of the displacements, and 20 time steps for the time-marching process to arrive at t = 120 (days). Radial displacements of the inner and outer surfaces of the homogeneous viscoelastic cylinder obtained by the present formulations are compared with analytical results of Flügge (1975) in Figs. 7 and 8, respectively. It may be observed that a good agreement exists between present results and results of Flügge (1975).

Fig. 7
figure 7

A comparison between present results and results of Flügge (1975) for time history of the radial displacement of inner surface of the homogeneous viscoelastic cylinder

Full size image
Fig. 8
figure 8

A comparison between present results and results of Flügge (1975) for time history of the radial displacement of the outer surface of the homogeneous viscoelastic cylinder

Full size image

Next, influences of the material properties heterogeneity on the radial displacements of the inner and outer surfaces of the functionally graded viscoelastic cylinder are investigated. In this regard, an exponential material properties gradation in the radial direction with γ = 0.2, λ 0 = 28 (MPa), μ 0 = 7 (MPa), υ = 0.4, and β = 14 is considered.

The radial displacements of the inner and outer surfaces of the exponentially graded viscoelastic cylindrical vessel predicted by the present approach are shown in Figs. 9 and 10, respectively. In comparison with the homogeneous cylindrical vessel, the radial displacements of the inner and outer surfaces are significantly smaller in the heterogeneous vessel, especially when comparing the responses for large enough times. It is worth reminding that the radial material gradation is practical in engineering manufacturing.

Fig. 9
figure 9

Time history of the radial displacement of the inner surface of the heterogeneous viscoelastic vessel, for an exponential radial gradation of the material properties

Full size image
Fig. 10
figure 10

Time history of the radial displacement of the outer surface of the heterogeneous viscoelastic vessel, for an exponential radial gradation of the material properties

Full size image

5.3 Example 3—Thermal loading of a viscoelastic cylindrical vessel

Finally, a thick hollow cylinder with exponential radial material properties gradation subjected to a radial thermal gradient is considered (Fig. 6). Geometric parameters of the present cylinder are identical with those of the foregoing example. Parameters of the exponential radial gradation of the mechanical and thermal properties of the materials are considered as: γ = 0.25 and 0.5, λ 0 = 1,000 (MPa), μ 0 = 90 (MPa), υ = 0.4, α = 10−5 (1/ °C), ρ = 3,500 (kg/m3), k = 85 (W/m °C) and β = 14. Temperature of the outer boundary is assumed to be 23 °C while temperature of the inner surface is maintained at 400 °C.

Due to the symmetry, only one quarter of the hollow cylinder was discretized. The displacement and stress fields are determined through using 42 spatial boundary elements and 2,000 time steps to arrive at t = 8,000 (days) through a time-marching process.

Time variations of the temperature of the outer surface of the cylinder are shown in Fig. 11. As shown in Fig. 11, as time elapses, the temperature tends to approach a steady-state quantity asymptotically. Time variations of the radial displacements of the inner and outer surfaces of the exponentially graded viscoelastic cylindrical vessel, obtained based on the present formulation, are shown for different gradation exponents in Figs. 12 and 13, respectively.

Fig. 11
figure 11

Time variations of the temperature of the outer surface of the heterogeneous viscoelastic vessel, for different heterogeneity exponents

Full size image
Fig. 12
figure 12

Time variations of the radial displacement of the inner surface of the heterogeneous viscoelastic vessel, for different heterogeneity exponents

Full size image
Fig. 13
figure 13

Time variations of the radial displacement of the outer surface of the heterogeneous viscoelastic vessel, for different heterogeneity exponents

Full size image

Moreover, time variations of the tangential stress of the inner and outer surfaces of the exponentially graded viscoelastic vessel are shown in Figs. 14 and 15, respectively. In comparison with a homogeneous vessel, trends of time variations of the stresses of the inner and outer surfaces are quite different in an exponentially graded cylinder. It may be inferred from results illustrated in Figs. 14 and 15 that generally, magnitude of the tangential stress increases with an increase in the heterogeneity exponent.

Fig. 14
figure 14

Time variations of the tangential stress of the inner surface of the heterogeneous viscoelastic vessel

Full size image
Fig. 15
figure 15

Time variations of the tangential stress of the outer surface of the heterogeneous viscoelastic vessel

Full size image

6 Concluding remarks

A boundary element formulation is presented for quasistatic thermoelastic analysis of components fabricated from heterogeneous exponentially graded viscoelastic materials. A new numerical solution scheme is introduced that is more appropriate for the present BEM-based thermoviscoelastic analysis of the FGMs. Influence of time variations of the temperatures which affects the viscoelastic properties of FGMs remarkably, is taken into account adequately by the presented BEM formulation. For validating the proposed formulation and presenting new practical results, three numerical examples are presented. Verification of the results is accomplished by comparing present results with analytical results of the homogeneous structures. Through these examples, differences between time variations of the displacements and stresses of the homogeneous and heterogeneous components are discussed. Present approach does not require employing relaxation functions and mathematical transformations, and it is capable of investigation of quasistatic thermoviscoelastic problems with arbitrary time-dependency of the loads and arbitrary boundary and initial conditions. One of the main advantages of the presented approach is its integral appearance that depends on the boundary values only. The resulting governing equations may be easily solved by adopting a linear time approximation for the velocity. In this regard, no domain approximations were assumed.