# Efficient Analytical Model for Time-Dependent Behavior of Layered Functionally Graded Plates with Viscoelastic Interlayers

Authors: Pengfei Liu; Peng Wu; Zhimin Wu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Layered functionally graded (FG) plate plays an important role in engineering constructions. In this work, an efficient analytical solution is proposed to investigate the time-dependent characters of layered FG plate with viscoelastic interlayers, in which each FG layer is described by the elasticity theory, and the interlayer viscoelasticity is simulated by the generalized Maxwell model. The constitutive equations in the interlayer are simplified, and then, the analytical solutions of stresses and displacements for the layered FG plate are solved by virtue of the efficient recursive matrix method. Some examples are analyzed to investigate the influences of geometric and material characteristics on the long-term behaviors for the FG plate.


---

## Body

## 1. Introduction

The layered plates are commonly used in various engineering [1–4], by right of their outstanding merits, such as high strength, lightweight, and corrosion resistance. Enough bonding stiffness between adjacent members is the guarantee of mechanical performance of layered plates. In many cases, the connection between adjacent layers is not rigid, and thus, the interfacial slip often happens [5, 6]. Besides, the bond behavior in layered plates exhibits viscoelastic property, due to the use of polymer adhesive [7, 8]. As a result, the mechanical performance of the whole layered plates is time-dependent under sustained loads [9].Although the layered plates exhibit good mechanical performance, it is necessary to reduce immediately the modulus difference between adjacent layers, which is a major cause of large interlaminar stress [10]. This requirement has led to the development of a novel material, called functionally graded (FG) materials, which can be designed to have continuous mechanical properties [11–13]. Due to its outstanding behaviors, FG materials have been applied in various areas, such as piezoelectric fiber-reinforced composite [14], heat-resistant materials [15], and nonuniform porous materials [16]. The concept of FG materials has also been introduced into the layered plates as face and/or transition layers, which can optimize the stress distribution and deformation [17, 18].A number of works have been proposed to study the mechanical responses of layered FG plates. Gunes et al. [19] presented a backpropagation artificial neural network method for investigating the three-dimensional free vibration response of an adhesively bonded wide and narrow FG plates. Demirbas and Apalak [20] performed a thermoelastic investigation of FG circular plates with adhesive bond and gave the thermal stress and strain distributions with different gradient indexes. Based on the equivalent single-layer higher-order theory, Tornabene et al. [21] analyzed the dynamic response of anisotropic doubly curved shells with arbitrary geometry and variable thickness. According to the four-variable refined plate theory, an analytical solution was deduced by Li et al. [22] for the thermomechanical bending analysis of FG sandwich plate. Wang et al. [23] derived an analytical solution for investigating the heat transfer behavior in FG Sandwich plates. Tornabene et al. [24, 25] used the higher-order theory to study the dynamic behavior of anisotropic doubly curved shells. An extended cohesive damage model was proposed by Ghimire and Chen [26] for the analysis of geometrical ratio influences on the failure mechanisms of FG sandwiches, and they found that the layered core can provide significant improvements in loading capacity. On the basis of the four-variable plate theory, Trinh el al. [27] presented a Levy solution with state space concept for static, dynamic, and buckling analyses for sandwich FG plates. Khorshidi and Karimi [28] established an analytical model based on a modified shear deformation theory for flutter investigation of sandwich plates with FG face sheets under thermal condition. A full layerwise method was employed by Nikbakht et al. [29] to obtain the yielding initiation of FG sandwich plates subjected to bi-sinusoidal distributed loading under general boundary conditions. Based on the first-order shear deformation theory, Tornabene [30] investigated the dynamic behavior of moderately thick functionally graded conical, cylindrical shells, and annular plates.In the above literature, most works in regard to layered plates were based on the assumption of perfectly bonding or static slip, while the time-dependent characters resulting from viscoelastic bonding interlayers were neglected. Furthermore, the solutions based on the simplified plate theories have considerable error for thick plates.This work proposes an efficient analytical model for investigating the time-dependent bending behavior of layered FG plates with considering viscoelastic bonding interlayer. In the analytical model, the mechanical behavior of each FG layer is expressed by the three-dimensional elasticity theory, which renounces the shear deformation assumption. Thus, it is highly accurate for structures with any thickness. The interlayer viscoelasticity is described by the generalized Maxwell model. The constitutive equations in the interlayer are simplified, and then, the analytical solutions of stresses and displacements for the layered FG plate are solved by virtue of the efficient recursive matrix method. Additionally, some examples provided by the present solution are conducted to investigate the influences of geometric and material characteristics on the time-dependent behaviors of the FG plate.


## 2. Efficient Analytical Model

In Figure1, a layered FG plate formed by p FG layers bonded by viscoelastic bonding interlayer is considered, which is simply supported and bears to a sustained load q(x,y) on its top. The three-dimensional size of the plate is a × b × H, and the thicknesses for FG layer and bonding interlayer are, respectively, hi and Δh, in which the label i means the layer index. The elastic modulus of each FG layer, denoted by Eiz, varies according to the exponential law along the thickness direction, and is expressed by(1)Eiz=E0iekiz−dib,where ki is the graded index and E0i=Eidib. The generalized Maxwell model, as shown in Figure 2, is employed to model the interlayer viscoelasticity, with time-dependent shear modulus given by(2)G∗t=G∞+∑j=1nGje−t/θG,j.Figure 1
Schematic diagram of the layered FG plate bonded by the viscoelastic interlayer.Figure 2
The configuration of the generalized Maxwell model (θi=ηi/Gj).


### 2.1. Governing Equations for a FG Layer

According to the elasticity theory, the constitutive equations for theith (i = 1, 2, …, p) FG layer are(3)σxi=λi+2Giεxi+λiεyi+λiεzi,τyzi=Giγyzi,σyi=λiεxi+λi+2Giεyi+λiεzi,τxzi=Giγxzi,σzi=λiεxi+λiεyi+λi+2Giεzi,τxyi=Giγxyi,where σxi, σyi, σzi, τxyi, τyzi, and τxzi are stresses, εxi, εyi, εzi, γxyi, γyzi, and γxzi are strains and λi and Gi are lame parameters expressed as follows:(4)λiz=μiEiz1+μi1−2μi,Giz=Eiz21+μi.The strain-displacement relations are given by(5)εxi=∂ui∂x,εyi=∂vi∂y,εzi=∂wi∂z,γxyi=∂ui∂y+∂vi∂x,γyzi=∂vi∂z+∂wi∂y,γxzi=∂ui∂z+∂wi∂x,where ui, vi, and wi represent the displacements in x, y, and z direction, respectively. The stress components of each FG layer should meet the following static equilibrium equations:(6)∂σxi∂x+∂τxyi∂y+∂τxzi∂z=0,∂σyi∂y+∂τxyi∂x+∂τyzi∂z=0,∂σzi∂z+∂τxzi∂x+∂τyzi∂y=0.The simply supported boundary conditions can be expressed by(7)vi=wi=σxi=0,atx=0,a,ui=wi=σxi=0,aty=0,b.Four steps are taken to obtain the differential equation involving displacement components. Firstly, by substituting (1) and (5) into (3) and eliminating the strain components, one has(8)σxi=E0iekiz−d0i1+μi1−2μi1−μi∂ui∂x+μi∂vi∂y+μi∂wi∂z,σyi=E0iekiz−d0i1+μi1−2μiμi∂ui∂x+1−μi∂vi∂y+μi∂wi∂z,σzi=E0iekiz−d0i1+μi1−2μiμi∂ui∂x+μi∂vi∂y+1−μi∂wi∂z,τxyi=E0iekiz−d0i21+μi∂ui∂y+∂vi∂x,τyzi=E0iekiz−d0i21+μi∂vi∂z+∂wi∂y,τxzi=E0iekiz−d0i21+μi∂ui∂z+∂wi∂x.By substituting (8) into (6) and eliminating the stress components, the governing equations only involving the displacement components are obtained:(9)21−μi1−2μi∂2ui∂x2+∂2ui∂y2+∂2ui∂z2+ki∂ui∂z+11−2μi∂2vi∂x∂y+11−2μi∂2wi∂x∂z+ki∂wi∂x=0,11−2μi∂2ui∂x∂y+∂2vi∂x2+21−μi1−2μi∂2vi∂y2+∂2vi∂z2+ki∂vi∂z+11−2μi∂2wi∂y∂z+ki∂wi∂y=0,11−2μi∂2ui∂x∂z+2kiμi1−2μi∂ui∂x+11−2μi∂2vi∂y∂z+21−2μikiμi∂vi∂y+∂2wi∂x2+∂2wi∂y2+21−μi1−2μi∂2wi∂z2+21−μi1−2μiki∂wi∂z=0.The present structural problem is displacement-based, and the stresses can be determined by (8) after the displacements are solved. The above governing equations are actually partial differential equations which can hard to be solved directly. Secondly, for simply supported boundary conditions, the displacements of each FG layer can be expanded in Fourier series as follows:(10)uix,y,z,tvix,y,z,twix,y,z,t=∑m=1∞∑n=1∞ui,mnz,tcosαmxsinβnyvi,mnz,tsinαmxcosβnywi,mnz,tsinαmxsinβny,where αm=mπ/a and βn=nπ/b. Thirdly, by applying the Fourier series expansion of (10) to (9), the partial differential equations become ordinary differential equations as follows:(11)∂2ui,mnz,t∂z2+ki∂ui,mnz,t∂z−21−μi1−2μiαm2+βn2ui,mnz,t−αmβn1−2μivi,mnz,t+αm1−2μi∂wi,mnz,t∂z+kiαmwi,mnz,t=0,−αmβn1−2μiui,mnz,t+∂2vi,mnz,t∂z2+ki∂vi,mnz,t∂z−αm2+21−μi1−2μiβn2vi,mnz,t+βn1−2μi∂wi,mnz,t∂z+kiβnwi,mnz,t=0,αm21−μi∂ui,mnz,t∂z+kiμiαm1−μiui,mnz,t+βn21−μi∂vi,mnz,t∂z+kiμiβn1−μivi,mnz,t−∂2wi,mnz,t∂z2−ki∂wi,mnz,t∂z+1−2μiαm2+βn221−μiwi,mnz,t=0.At last, a fourth-order differential equation ofwi,mn is derived out by eliminating ui,mn and vi,mn in (11):(12)∂4wi,mn∂z4+2ki∂3wi,mn∂z3+ki2−2αmn2∂2wi,mn∂z2−2kiαmn2∂wi,mn∂z+αmn4+k2αmn2μi1−μiwi,mn=0.The general solution ofwi,mn is(13)wmniz,t=∑j=14egmnijzCmnijt,where Cmn,jit are the undetermined time-dependent coefficients, which can be obtained according to the boundary and continuity conditions in the subsequent sections; the details of gmn,ji are defined in Appendix A. By substituting into (13) and (11), ui,mn and vi,mn are obtained:(14)umniz,t=∑j=14fmnijegmnijzCmnijt+∑j=56egmnijzCmnijt,vmniz,t=βnαm∑j=14fmnijegmnijzCmnijt−αmβn∑j=56egmnijzCmnijt,where the details of fmnij are listed in Appendix A. By substituting of (13), (14), and (8) into (9), the stress components involving undetermined coefficients can be written as(15)σx,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αm1−μifmnij−μiβn2fmnijαm+μigmnijegmnijzCmnijt+∑j=56αm2μi−1egmnijzCmnijt,σy,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αmμifmnij−1−μiβn2fmnijαm+μigmnijegmnijziCmnijt+∑j=56αm1−2μiegmnijziCmnijt,σz,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αmμifmnij−μiβn2fmnijαm+1−μigmnijegmnijzCmnijt,τxy,mniz,t=E0iekiz−d0i21+μi∑j=142βnfmnijegmnijzCmnijt+∑j=56βn−αm2βnegmnijzCmnijt,τyz,mniz,t=E0iekiz−d0i21+μi∑j=14βnfmnijgmnijαm+βnegmnijzCmnijt−∑j=56αmgmnijβnegmnijzCmnijt,τxz,mniz,t=E0iekiz−d0i21+μi∑j=14fmnijgmnij+αmegmnijzCmnijt+∑j=56gmnijegmnijzCmnijt.


### 2.2. Governing Equations for a Viscoelastic Interlayer

According to the Boltzmann superposition principle, the constitutive equations for theith (i = 1, 2, …, p-1) interlayer are given by(16)τxz∗ix,y,t=∫−∞tG∗t−ξ∂γxzi∗x,y,ξ∂ξdξ,τyz∗ix,y,t=∫−∞tG∗t−ξ∂γyzi∗x,y,ξ∂ξdξ.The above convolution integral equations means that the stress of viscoelastic interlayer depends on the total strain history, which leads to the heavy calculation and time-consuming. In the present work, they are simplified as(17)τxzi∗x,y,t=G∗tγxzi∗x,y,t,τyzi∗x,y,t=G∗tγyzi∗x,y,t.In comparison with the exact solution, this simplification leads to small error, and the present solution is always on the side of safety [25]. Considering that Δh is far less than hi, the shear strains can be assumed to be constant in z direction and given by(18)γxzi∗x,y,t=ui+1x,y,di+1b,t−uix,y,dit,tΔh,γyzi∗x,y,t=vi+1x,y,di+1b,t−vix,y,dit,tΔh.The shear stress continuity relationships between adjacent layers are(19)τxzix,y,dit,t=τxzi∗x,y,t=τxzi+1x,y,di+1b,t,τyzix,y,dit,t=τyzi∗x,y,t=τyzi+1x,y,di+1b,t.


### 2.3. Recursive Matrix Method

The stresses on the top and bottom surfaces of the plate are(20)σzpx,y,H,t=−qx,y,τxzpx,y,H,t=0,τyzpx,y,H,t=0,σz1x,y,0,t=0,τxz1x,y,0,t=0,τyz1x,y,0,t=0.Since the stresses and displacements are in series form,q(x,y) should also be expanded as(21)qx,y=∑m=1∞∑n=1∞qmnsinαmxsinβny,qmn=−4ab∫0a∫0bqx,ysinαmxsinβnydxdy.By substitution of (13)–(15) into (8), the general solution of stresses and displacements can be transferred into the matrix form as(22)Φmniz,t=MmnizΛmnit,where(23)Φmniz,t=ui,mnz,tvi,mnz,twi,mnz,tσzi,mnz,tτxzi,mnz,tτyzi,mnz,tT,Λmnit=Cmn,1itCmn,2itCmn,3itCmn,4itCmn,5itCmn,6itT,Mmniz=Rmn11zRmn12zRmn13zRmn14zRmn15zRmn16zRmn21zRmn22zRmn23zRmn24zRmn25zRmn26zRmn31zRmn32zRmn33zRmn34z00Rmn41zRmn42zRmn43zRmn44z00Rmn51zRmn52zRmn53zRmn54zRmn55zRmn56zRmn61zRmn62zRmn63zRmn64zRmn65zRmn66z,where the nonzero coefficients in matrix Mmniz can be found in Appendix B. By combining (17)–(19), the continuity relationships between the adjacent FG layers are rewritten in the matrix form as(24)Φmni+1di+1b,t=ΨtΦmnidit,t,where(25)Ψt=10−Δhαm0ΔhG∗t001−Δhβn00ΔhG∗t001000000100000010000001.By substitutingdib and dit into the z-coordinate in (24), respectively, we obtain(26)Φmnidib,t=MmnidibΛmnit,Φmnidit,t=MmniditΛmnit.Elimination ofΛmnit in (26) yields(27)Φmnidit,t=MmniditMmnidib−1Φmnidib,t.By reusing (24) and (27), from i = 1 to i = p, one obtains(28)Φmnpdpt,t=∏i=p2MmniditMmnidib−1ΨtMmn1d1tMmn1d1b−1Φmn1d1b,t.Four 3 × 3, submatrices are used to define the matrix multiplication in (28) as follows:(29)Cmn,11iCmn,12iCmn,21iCmn,22i=∏i=p2MmniditMmnidib−1ΨtMmn1d1tMmn1d1b−1.Thus, (28) can be written as(30)Φmnpdpt,t=Cmn,11iCmn,12iCmn,21iCmn,22iΦmnpdpt,t.By the decomposition of (30), two submatrix equations are given as(31)up,mnd1b,tvxzp,mnd1b,twyzp,mnd1b,t=Cmn,11iu1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t+Cmn,12iσz1,mnd1b,tτxz1,mnd1b,tτyz1,mnd1b,t,σzp,mnd1b,tτxzp,mnd1b,tτyzp,mnd1b,t=Cmn,21iu1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t+Cmn,22iσz1,mnd1b,tτxz1,mnd1b,tτyz1,mnd1b,t.By solving (31), the displacement boundary value on the bottom surface of the structure can be expressed as(32)u1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t=Cmn,21i−1qmn00.Similar to the relationship in (28), Φmnidit,t for any FG layer is obtained from (24) and (27) as follows:(33)Φmnidit,t=∏j=i2MmnjdjtMmnjdjt−1ΨtMmn1djtMmn1djt−1Φmn1d1b,t.The time-dependent coefficients of theith FG layer are further obtained as(34)Λmnit=Mmnid1t−1∏j=i2Mmnjd1tMmnjd1t−1ΨtMmn1d1tMmn1d1t−1Φmn1d1b,t.Finally, the substitution of (34) into (13)–(15) yields the solution of stress and displacement components for each FG layer.


## 2.1. Governing Equations for a FG Layer

According to the elasticity theory, the constitutive equations for theith (i = 1, 2, …, p) FG layer are(3)σxi=λi+2Giεxi+λiεyi+λiεzi,τyzi=Giγyzi,σyi=λiεxi+λi+2Giεyi+λiεzi,τxzi=Giγxzi,σzi=λiεxi+λiεyi+λi+2Giεzi,τxyi=Giγxyi,where σxi, σyi, σzi, τxyi, τyzi, and τxzi are stresses, εxi, εyi, εzi, γxyi, γyzi, and γxzi are strains and λi and Gi are lame parameters expressed as follows:(4)λiz=μiEiz1+μi1−2μi,Giz=Eiz21+μi.The strain-displacement relations are given by(5)εxi=∂ui∂x,εyi=∂vi∂y,εzi=∂wi∂z,γxyi=∂ui∂y+∂vi∂x,γyzi=∂vi∂z+∂wi∂y,γxzi=∂ui∂z+∂wi∂x,where ui, vi, and wi represent the displacements in x, y, and z direction, respectively. The stress components of each FG layer should meet the following static equilibrium equations:(6)∂σxi∂x+∂τxyi∂y+∂τxzi∂z=0,∂σyi∂y+∂τxyi∂x+∂τyzi∂z=0,∂σzi∂z+∂τxzi∂x+∂τyzi∂y=0.The simply supported boundary conditions can be expressed by(7)vi=wi=σxi=0,atx=0,a,ui=wi=σxi=0,aty=0,b.Four steps are taken to obtain the differential equation involving displacement components. Firstly, by substituting (1) and (5) into (3) and eliminating the strain components, one has(8)σxi=E0iekiz−d0i1+μi1−2μi1−μi∂ui∂x+μi∂vi∂y+μi∂wi∂z,σyi=E0iekiz−d0i1+μi1−2μiμi∂ui∂x+1−μi∂vi∂y+μi∂wi∂z,σzi=E0iekiz−d0i1+μi1−2μiμi∂ui∂x+μi∂vi∂y+1−μi∂wi∂z,τxyi=E0iekiz−d0i21+μi∂ui∂y+∂vi∂x,τyzi=E0iekiz−d0i21+μi∂vi∂z+∂wi∂y,τxzi=E0iekiz−d0i21+μi∂ui∂z+∂wi∂x.By substituting (8) into (6) and eliminating the stress components, the governing equations only involving the displacement components are obtained:(9)21−μi1−2μi∂2ui∂x2+∂2ui∂y2+∂2ui∂z2+ki∂ui∂z+11−2μi∂2vi∂x∂y+11−2μi∂2wi∂x∂z+ki∂wi∂x=0,11−2μi∂2ui∂x∂y+∂2vi∂x2+21−μi1−2μi∂2vi∂y2+∂2vi∂z2+ki∂vi∂z+11−2μi∂2wi∂y∂z+ki∂wi∂y=0,11−2μi∂2ui∂x∂z+2kiμi1−2μi∂ui∂x+11−2μi∂2vi∂y∂z+21−2μikiμi∂vi∂y+∂2wi∂x2+∂2wi∂y2+21−μi1−2μi∂2wi∂z2+21−μi1−2μiki∂wi∂z=0.The present structural problem is displacement-based, and the stresses can be determined by (8) after the displacements are solved. The above governing equations are actually partial differential equations which can hard to be solved directly. Secondly, for simply supported boundary conditions, the displacements of each FG layer can be expanded in Fourier series as follows:(10)uix,y,z,tvix,y,z,twix,y,z,t=∑m=1∞∑n=1∞ui,mnz,tcosαmxsinβnyvi,mnz,tsinαmxcosβnywi,mnz,tsinαmxsinβny,where αm=mπ/a and βn=nπ/b. Thirdly, by applying the Fourier series expansion of (10) to (9), the partial differential equations become ordinary differential equations as follows:(11)∂2ui,mnz,t∂z2+ki∂ui,mnz,t∂z−21−μi1−2μiαm2+βn2ui,mnz,t−αmβn1−2μivi,mnz,t+αm1−2μi∂wi,mnz,t∂z+kiαmwi,mnz,t=0,−αmβn1−2μiui,mnz,t+∂2vi,mnz,t∂z2+ki∂vi,mnz,t∂z−αm2+21−μi1−2μiβn2vi,mnz,t+βn1−2μi∂wi,mnz,t∂z+kiβnwi,mnz,t=0,αm21−μi∂ui,mnz,t∂z+kiμiαm1−μiui,mnz,t+βn21−μi∂vi,mnz,t∂z+kiμiβn1−μivi,mnz,t−∂2wi,mnz,t∂z2−ki∂wi,mnz,t∂z+1−2μiαm2+βn221−μiwi,mnz,t=0.At last, a fourth-order differential equation ofwi,mn is derived out by eliminating ui,mn and vi,mn in (11):(12)∂4wi,mn∂z4+2ki∂3wi,mn∂z3+ki2−2αmn2∂2wi,mn∂z2−2kiαmn2∂wi,mn∂z+αmn4+k2αmn2μi1−μiwi,mn=0.The general solution ofwi,mn is(13)wmniz,t=∑j=14egmnijzCmnijt,where Cmn,jit are the undetermined time-dependent coefficients, which can be obtained according to the boundary and continuity conditions in the subsequent sections; the details of gmn,ji are defined in Appendix A. By substituting into (13) and (11), ui,mn and vi,mn are obtained:(14)umniz,t=∑j=14fmnijegmnijzCmnijt+∑j=56egmnijzCmnijt,vmniz,t=βnαm∑j=14fmnijegmnijzCmnijt−αmβn∑j=56egmnijzCmnijt,where the details of fmnij are listed in Appendix A. By substituting of (13), (14), and (8) into (9), the stress components involving undetermined coefficients can be written as(15)σx,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αm1−μifmnij−μiβn2fmnijαm+μigmnijegmnijzCmnijt+∑j=56αm2μi−1egmnijzCmnijt,σy,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αmμifmnij−1−μiβn2fmnijαm+μigmnijegmnijziCmnijt+∑j=56αm1−2μiegmnijziCmnijt,σz,mniz,t=E0iekiz−d0i1+μi1−2μi∑j=14−αmμifmnij−μiβn2fmnijαm+1−μigmnijegmnijzCmnijt,τxy,mniz,t=E0iekiz−d0i21+μi∑j=142βnfmnijegmnijzCmnijt+∑j=56βn−αm2βnegmnijzCmnijt,τyz,mniz,t=E0iekiz−d0i21+μi∑j=14βnfmnijgmnijαm+βnegmnijzCmnijt−∑j=56αmgmnijβnegmnijzCmnijt,τxz,mniz,t=E0iekiz−d0i21+μi∑j=14fmnijgmnij+αmegmnijzCmnijt+∑j=56gmnijegmnijzCmnijt.


## 2.2. Governing Equations for a Viscoelastic Interlayer

According to the Boltzmann superposition principle, the constitutive equations for theith (i = 1, 2, …, p-1) interlayer are given by(16)τxz∗ix,y,t=∫−∞tG∗t−ξ∂γxzi∗x,y,ξ∂ξdξ,τyz∗ix,y,t=∫−∞tG∗t−ξ∂γyzi∗x,y,ξ∂ξdξ.The above convolution integral equations means that the stress of viscoelastic interlayer depends on the total strain history, which leads to the heavy calculation and time-consuming. In the present work, they are simplified as(17)τxzi∗x,y,t=G∗tγxzi∗x,y,t,τyzi∗x,y,t=G∗tγyzi∗x,y,t.In comparison with the exact solution, this simplification leads to small error, and the present solution is always on the side of safety [25]. Considering that Δh is far less than hi, the shear strains can be assumed to be constant in z direction and given by(18)γxzi∗x,y,t=ui+1x,y,di+1b,t−uix,y,dit,tΔh,γyzi∗x,y,t=vi+1x,y,di+1b,t−vix,y,dit,tΔh.The shear stress continuity relationships between adjacent layers are(19)τxzix,y,dit,t=τxzi∗x,y,t=τxzi+1x,y,di+1b,t,τyzix,y,dit,t=τyzi∗x,y,t=τyzi+1x,y,di+1b,t.


## 2.3. Recursive Matrix Method

The stresses on the top and bottom surfaces of the plate are(20)σzpx,y,H,t=−qx,y,τxzpx,y,H,t=0,τyzpx,y,H,t=0,σz1x,y,0,t=0,τxz1x,y,0,t=0,τyz1x,y,0,t=0.Since the stresses and displacements are in series form,q(x,y) should also be expanded as(21)qx,y=∑m=1∞∑n=1∞qmnsinαmxsinβny,qmn=−4ab∫0a∫0bqx,ysinαmxsinβnydxdy.By substitution of (13)–(15) into (8), the general solution of stresses and displacements can be transferred into the matrix form as(22)Φmniz,t=MmnizΛmnit,where(23)Φmniz,t=ui,mnz,tvi,mnz,twi,mnz,tσzi,mnz,tτxzi,mnz,tτyzi,mnz,tT,Λmnit=Cmn,1itCmn,2itCmn,3itCmn,4itCmn,5itCmn,6itT,Mmniz=Rmn11zRmn12zRmn13zRmn14zRmn15zRmn16zRmn21zRmn22zRmn23zRmn24zRmn25zRmn26zRmn31zRmn32zRmn33zRmn34z00Rmn41zRmn42zRmn43zRmn44z00Rmn51zRmn52zRmn53zRmn54zRmn55zRmn56zRmn61zRmn62zRmn63zRmn64zRmn65zRmn66z,where the nonzero coefficients in matrix Mmniz can be found in Appendix B. By combining (17)–(19), the continuity relationships between the adjacent FG layers are rewritten in the matrix form as(24)Φmni+1di+1b,t=ΨtΦmnidit,t,where(25)Ψt=10−Δhαm0ΔhG∗t001−Δhβn00ΔhG∗t001000000100000010000001.By substitutingdib and dit into the z-coordinate in (24), respectively, we obtain(26)Φmnidib,t=MmnidibΛmnit,Φmnidit,t=MmniditΛmnit.Elimination ofΛmnit in (26) yields(27)Φmnidit,t=MmniditMmnidib−1Φmnidib,t.By reusing (24) and (27), from i = 1 to i = p, one obtains(28)Φmnpdpt,t=∏i=p2MmniditMmnidib−1ΨtMmn1d1tMmn1d1b−1Φmn1d1b,t.Four 3 × 3, submatrices are used to define the matrix multiplication in (28) as follows:(29)Cmn,11iCmn,12iCmn,21iCmn,22i=∏i=p2MmniditMmnidib−1ΨtMmn1d1tMmn1d1b−1.Thus, (28) can be written as(30)Φmnpdpt,t=Cmn,11iCmn,12iCmn,21iCmn,22iΦmnpdpt,t.By the decomposition of (30), two submatrix equations are given as(31)up,mnd1b,tvxzp,mnd1b,twyzp,mnd1b,t=Cmn,11iu1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t+Cmn,12iσz1,mnd1b,tτxz1,mnd1b,tτyz1,mnd1b,t,σzp,mnd1b,tτxzp,mnd1b,tτyzp,mnd1b,t=Cmn,21iu1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t+Cmn,22iσz1,mnd1b,tτxz1,mnd1b,tτyz1,mnd1b,t.By solving (31), the displacement boundary value on the bottom surface of the structure can be expressed as(32)u1,mnd1b,tvxz1,mnd1b,twyz1,mnd1b,t=Cmn,21i−1qmn00.Similar to the relationship in (28), Φmnidit,t for any FG layer is obtained from (24) and (27) as follows:(33)Φmnidit,t=∏j=i2MmnjdjtMmnjdjt−1ΨtMmn1djtMmn1djt−1Φmn1d1b,t.The time-dependent coefficients of theith FG layer are further obtained as(34)Λmnit=Mmnid1t−1∏j=i2Mmnjd1tMmnjd1t−1ΨtMmn1d1tMmn1d1t−1Φmn1d1b,t.Finally, the substitution of (34) into (13)–(15) yields the solution of stress and displacement components for each FG layer.


## 3. Numerical Results and Discussion

In the following, the FG sandwich plate (p = 3) is taken as an example. The variables with the superscript r, e.g., σxr, are defined as the result of stress or displacement at the point of x=0⋅25a, y=0⋅25b, and z=h1. The variables with two vertical lines denote their absolute values, e.g., wi.


### 3.1. Validation of the Present Solution

The convergence property of the present solution is assessed first. The series are truncated into a finite numberM for actual calculations here. The parameters in the FG plate are taken as q(x,y) = 1 N/mm2, a = 1000 mm, b = 800 mm, ∆h = 0.5 mm, h1 = h3 = 30 mm, h2 = 40 mm, E01  =  E03  = 80 GPa, E02  = 40 GPa, k1 = k3 = 0.05, k2 = 0, and μ1  =  μ2  =  μ3  =  μ∗  = 0.3. The material of polyvinyl butyral (PVB) is chosen for the viscoelastic interlayer and its viscoelastic parameters of which are given in Table 1. Table 2 lists the present results with different series terms, respectively. It is found that the present results tend to be constant as series terms increase, and the convergence accuracy is four significant digits.Table 1
The relaxation moduli and relaxation time of the generalized Maxwell model for PVB material.
JGj∗ (MPa)θG,j (s)175.64263.256 × 10−11237.06774.949 × 10−93137.15527.243 × 10−8433.51409.864 × 10−65126.60482.806 × 10−3642.19501.644 × 10−1714.21622.265 × 10083.58223.536 × 10190.45389.368 × 103100.19126.414 × 105110.28934.135 × 107∞0.0880Table 2
Convergence analysis of the present method ast = 1 day, 1 year, and 10 years, respectively.
TMur (mm)wr (mm)σxr (MPa)σzr (MPa)τxyr (MPa)τxzr (MPa)1 day10.03731−1.068−67.96−0.289640.41−0.106550.03582−1.088−74.98−0.393338.20−0.100990.03583−1.088−74.94−0.391638.20−0.1010130.03583−1.088−74.94−0.391638.20−0.10101 year10.04072−1.152−74.16−0.305644.10−0.0418350.03866−1.180−83.43−0.437941.03−0.0394090.03867−1.180−83.30−0.430741.04−0.03946130.03867−1.180−83.30−0.431041.04−0.0394610 years10.04226−1.190−76.97−0.312945.77−0.0119650.03977−1.222−87.06−0.436642.20−0.0112190.03984−1.222−86.65−0.414342.23−0.01124130.03983−1.222−86.67−0.416442.23−0.01124Note.M means the series terms.The present solution is compared with the Kirchhoff–Love (KL) solution [26], respectively. The parameters in the FG plate are fixed at q(x,y) = sin(πx/a)sin(πy/b) N/mm2, a = b = 1000 mm, ∆h = 0.2 mm, k2 = 0, E03  =  E01ek1h1  = 80 GPa, E02  = 40 GPa, μ1  =  μ2  =  μ3  =  μ∗  = 0.3, t = 104 s, G∞∗  = 0.5 MPa, G1∗  = 500 MPa, θG,1  = 1 s, h1:h2:h3 = 1 : 2:1, and k3 = −k1. The above parameters mean that the geometric and material properties of the present structure are symmetric about the midplane in the direction of thickness. Figure 3 represents the relative error between KL solution and the present one for a different length-to-height ratio a/H and graded index k1, respectively. A good agreement is found for KL solution in thin plate case with small |k1|, while the relative error of KL results increases gradually as the plate become thick or k1 decreases. The errors for σxr, τxzr, and wr are 15.4%, 6.6%, and 15.9%, respectively, as a/H = 8 and k1 = −0.09. The imprecision of KL solutions mainly results from the neglect of the transverse shear deformation and the deviation of neutral plane of the FG layer.Figure 3
Relative errors of the KL solution compared with the present one for different length-to-height ratios ask1 = −0.03, −0.05, −0.07, and −0.09 when t = 10000 s. Note. The errors are defined by |[(KL)-Present]/Present| ×100%. (a) σxr, (b) τxzr, and (c) wr.
(a)(b)(c)


### 3.2. Parameter Research

In this section, some examples provided by the present solution are conducted to investigate the influences of geometric and material characteristics on the time-dependent behaviors of the FG plate. Some of the parameters are fixed atq(x,y) = sin(πx/a)sin(πy/b) N/mm2, a = b = 1000 mm, k1:k3 = −1, k2 = 0, E03  =  E01ek1h1, E02  = 25Gpa, and μ1  =  μ2  =  μ3  =  μ∗  = 0.3; the viscoelastic constants are given in Table 1, while the others are variable.Figure4 plots the stress and displacement distribution in z direction when t = 1 day, 1 year, and 10 years and the corresponding results in perfectly bonded (PB) case. The constants are fixed at h1 = h3 = 20 mm, h2 = 40 mm, ∆h = 0.2 mm, and k1 = −0.05. It can be obtained from Figure 4 that σxi, τxzi, τxyi, ui, and wi considerably increase with t, while σyi changes slightly with time. σxi and τxyi show zig-zag distributions, and they, in the facial layers, are obviously curve distribution, which is different from isotropic material. ui also show a zig-zag distribution but keeps straight in each FG layer. τxzi gives a multipeak distribution. Compared with the results of PB case, the maximum values of σxi, τxzi, τxyi, ui, and wi increase by 200.1%, 26.2%, 200.2%, 604.7%, and 1221%, respectively, at t = 10 years.Figure 4
The distribution of the stress and displacement in the thickness direction at different time periods. (a)σxix=0.5a,y=0.5b, (b) σyix=0.5a,y=0.5b, (c) τxzix=0,y=0.5b, (d) τxzix=0,y=0, (e) uix=0,y=0.5b, and (f) wix=0.5a,y=0.5b.
(a)(b)(c)(d)(e)(f)The effect of the graded index on the stress and displacement of the FG plate, as well as the elastic modulus distributions along the thickness direction, are illustrated in Figure5. The parameters are taken as h1 = h3 = 10 mm, h2 = 30 mm, ∆h = 0.2 mm, and t = 1 day, and the average modulus of FG facial layer, i.e., ∫0hiE0iekiydy/hi (i = 1, 3), is fixed at 80 GPa. From the results shown in Figures 5(b) and 5(c), the decline of k1 reduces the maximum value of τxzi and wi but enlarges that of σxi. As k1 goes down from 0.2 to -0.2, τxzimax and wimax fall by 10.0% and 15.7%, respectively, while σximax increases by 156.7%. In the other word, the stress and displacement distribution can be optimized by adjusting the graded index.Figure 5
The influence of graded index on the distribution of the elastic modulus, stress, and displacement along the thickness direction with fixed material consumption. (a)Eiz, (b) σxi, (c) τxzi, and (d) wi.
(a)(b)(c)(d)Figure6 shows the influences of t, k1, and Δh on the maximum value of the stress and displacement. The parameters are fixed at the same as those in Figure 5, except t, k1 and Δh are variable. It can be seen from Figure 6 that, as k1 = −0.1, σximax, τxzimax, and wimax increase monotonously and keeps invariant as t or Δh increases. Due to the increase of t and Δh resulting in a reduced shear modulus of the interlayer, and as t⟶∞, G∗t approaches to the fixed value G∞∗. Similarly, τxzimax and wimax increase monotonously and tend to definite values with the increase of k1 or t as Δh = 0.2 mm, which results from that the difference in elastic modulus between the adjacent layer surfaces grows gradually as k1 increases. For a given t, σximax decreases first and then increases as k1 grows, which is owing to the fact that the location of σximax transfers from the outer surface of the facial layer to the inner surface in the early stage, while in the medium or long term, the location of σximax transfers from the facial layer to the core layer as k1 increases. Similarly, for a small certain k1, σximax increases monotonously with t, while for a large k1, σximax firstly decreases and then increases as t goes on, which is due to the transfer of the σximax location from the facial layer to the core layer.Figure 6
The influence of time and graded index or the interlayer thickness on the maximum value of the normal stress, shear stress, and deflection. (a)σximaxΔh=0.2mm, (b) σximaxk1=−0.1, (c) τxzimaxΔh=0.2mm, (d) τxzimaxk1=−0.1, (e) wimaxΔh=0.2mm, and (f) wimaxk1=−0.1.
(a)(b)(c)(d)(e)(f)


## 3.1. Validation of the Present Solution

The convergence property of the present solution is assessed first. The series are truncated into a finite numberM for actual calculations here. The parameters in the FG plate are taken as q(x,y) = 1 N/mm2, a = 1000 mm, b = 800 mm, ∆h = 0.5 mm, h1 = h3 = 30 mm, h2 = 40 mm, E01  =  E03  = 80 GPa, E02  = 40 GPa, k1 = k3 = 0.05, k2 = 0, and μ1  =  μ2  =  μ3  =  μ∗  = 0.3. The material of polyvinyl butyral (PVB) is chosen for the viscoelastic interlayer and its viscoelastic parameters of which are given in Table 1. Table 2 lists the present results with different series terms, respectively. It is found that the present results tend to be constant as series terms increase, and the convergence accuracy is four significant digits.Table 1
The relaxation moduli and relaxation time of the generalized Maxwell model for PVB material.
JGj∗ (MPa)θG,j (s)175.64263.256 × 10−11237.06774.949 × 10−93137.15527.243 × 10−8433.51409.864 × 10−65126.60482.806 × 10−3642.19501.644 × 10−1714.21622.265 × 10083.58223.536 × 10190.45389.368 × 103100.19126.414 × 105110.28934.135 × 107∞0.0880Table 2
Convergence analysis of the present method ast = 1 day, 1 year, and 10 years, respectively.
TMur (mm)wr (mm)σxr (MPa)σzr (MPa)τxyr (MPa)τxzr (MPa)1 day10.03731−1.068−67.96−0.289640.41−0.106550.03582−1.088−74.98−0.393338.20−0.100990.03583−1.088−74.94−0.391638.20−0.1010130.03583−1.088−74.94−0.391638.20−0.10101 year10.04072−1.152−74.16−0.305644.10−0.0418350.03866−1.180−83.43−0.437941.03−0.0394090.03867−1.180−83.30−0.430741.04−0.03946130.03867−1.180−83.30−0.431041.04−0.0394610 years10.04226−1.190−76.97−0.312945.77−0.0119650.03977−1.222−87.06−0.436642.20−0.0112190.03984−1.222−86.65−0.414342.23−0.01124130.03983−1.222−86.67−0.416442.23−0.01124Note.M means the series terms.The present solution is compared with the Kirchhoff–Love (KL) solution [26], respectively. The parameters in the FG plate are fixed at q(x,y) = sin(πx/a)sin(πy/b) N/mm2, a = b = 1000 mm, ∆h = 0.2 mm, k2 = 0, E03  =  E01ek1h1  = 80 GPa, E02  = 40 GPa, μ1  =  μ2  =  μ3  =  μ∗  = 0.3, t = 104 s, G∞∗  = 0.5 MPa, G1∗  = 500 MPa, θG,1  = 1 s, h1:h2:h3 = 1 : 2:1, and k3 = −k1. The above parameters mean that the geometric and material properties of the present structure are symmetric about the midplane in the direction of thickness. Figure 3 represents the relative error between KL solution and the present one for a different length-to-height ratio a/H and graded index k1, respectively. A good agreement is found for KL solution in thin plate case with small |k1|, while the relative error of KL results increases gradually as the plate become thick or k1 decreases. The errors for σxr, τxzr, and wr are 15.4%, 6.6%, and 15.9%, respectively, as a/H = 8 and k1 = −0.09. The imprecision of KL solutions mainly results from the neglect of the transverse shear deformation and the deviation of neutral plane of the FG layer.Figure 3
Relative errors of the KL solution compared with the present one for different length-to-height ratios ask1 = −0.03, −0.05, −0.07, and −0.09 when t = 10000 s. Note. The errors are defined by |[(KL)-Present]/Present| ×100%. (a) σxr, (b) τxzr, and (c) wr.
(a)(b)(c)


## 3.2. Parameter Research

In this section, some examples provided by the present solution are conducted to investigate the influences of geometric and material characteristics on the time-dependent behaviors of the FG plate. Some of the parameters are fixed atq(x,y) = sin(πx/a)sin(πy/b) N/mm2, a = b = 1000 mm, k1:k3 = −1, k2 = 0, E03  =  E01ek1h1, E02  = 25Gpa, and μ1  =  μ2  =  μ3  =  μ∗  = 0.3; the viscoelastic constants are given in Table 1, while the others are variable.Figure4 plots the stress and displacement distribution in z direction when t = 1 day, 1 year, and 10 years and the corresponding results in perfectly bonded (PB) case. The constants are fixed at h1 = h3 = 20 mm, h2 = 40 mm, ∆h = 0.2 mm, and k1 = −0.05. It can be obtained from Figure 4 that σxi, τxzi, τxyi, ui, and wi considerably increase with t, while σyi changes slightly with time. σxi and τxyi show zig-zag distributions, and they, in the facial layers, are obviously curve distribution, which is different from isotropic material. ui also show a zig-zag distribution but keeps straight in each FG layer. τxzi gives a multipeak distribution. Compared with the results of PB case, the maximum values of σxi, τxzi, τxyi, ui, and wi increase by 200.1%, 26.2%, 200.2%, 604.7%, and 1221%, respectively, at t = 10 years.Figure 4
The distribution of the stress and displacement in the thickness direction at different time periods. (a)σxix=0.5a,y=0.5b, (b) σyix=0.5a,y=0.5b, (c) τxzix=0,y=0.5b, (d) τxzix=0,y=0, (e) uix=0,y=0.5b, and (f) wix=0.5a,y=0.5b.
(a)(b)(c)(d)(e)(f)The effect of the graded index on the stress and displacement of the FG plate, as well as the elastic modulus distributions along the thickness direction, are illustrated in Figure5. The parameters are taken as h1 = h3 = 10 mm, h2 = 30 mm, ∆h = 0.2 mm, and t = 1 day, and the average modulus of FG facial layer, i.e., ∫0hiE0iekiydy/hi (i = 1, 3), is fixed at 80 GPa. From the results shown in Figures 5(b) and 5(c), the decline of k1 reduces the maximum value of τxzi and wi but enlarges that of σxi. As k1 goes down from 0.2 to -0.2, τxzimax and wimax fall by 10.0% and 15.7%, respectively, while σximax increases by 156.7%. In the other word, the stress and displacement distribution can be optimized by adjusting the graded index.Figure 5
The influence of graded index on the distribution of the elastic modulus, stress, and displacement along the thickness direction with fixed material consumption. (a)Eiz, (b) σxi, (c) τxzi, and (d) wi.
(a)(b)(c)(d)Figure6 shows the influences of t, k1, and Δh on the maximum value of the stress and displacement. The parameters are fixed at the same as those in Figure 5, except t, k1 and Δh are variable. It can be seen from Figure 6 that, as k1 = −0.1, σximax, τxzimax, and wimax increase monotonously and keeps invariant as t or Δh increases. Due to the increase of t and Δh resulting in a reduced shear modulus of the interlayer, and as t⟶∞, G∗t approaches to the fixed value G∞∗. Similarly, τxzimax and wimax increase monotonously and tend to definite values with the increase of k1 or t as Δh = 0.2 mm, which results from that the difference in elastic modulus between the adjacent layer surfaces grows gradually as k1 increases. For a given t, σximax decreases first and then increases as k1 grows, which is owing to the fact that the location of σximax transfers from the outer surface of the facial layer to the inner surface in the early stage, while in the medium or long term, the location of σximax transfers from the facial layer to the core layer as k1 increases. Similarly, for a small certain k1, σximax increases monotonously with t, while for a large k1, σximax firstly decreases and then increases as t goes on, which is due to the transfer of the σximax location from the facial layer to the core layer.Figure 6
The influence of time and graded index or the interlayer thickness on the maximum value of the normal stress, shear stress, and deflection. (a)σximaxΔh=0.2mm, (b) σximaxk1=−0.1, (c) τxzimaxΔh=0.2mm, (d) τxzimaxk1=−0.1, (e) wimaxΔh=0.2mm, and (f) wimaxk1=−0.1.
(a)(b)(c)(d)(e)(f)


## 4. Conclusion

The elasticity theory with the recursive matrix method was proposed to analyze the time-dependent behavior of the layered functionally graded plates with viscoelastic interlayer. The following conclusions can be provided:(1)
In the thin plate case, the KL solution with small graded index approaches the present one. Since the transverse shear deformation is neglected and the neutral plane gradually deviates from the middle plane, the relative errors increase gradually as length-to-height ratio or graded index increases.(2)
As time goes on, the influence of the interlayer shear modulus degeneration for the adhesive bonding case is obvious, which leads to the reduced interface shear stress and the increased interfacial slip as well as the greater deflection.(3)
By adjusting the graded index, the stress and displacement distribution of the case with fixed material consumption can be optimized and the location of the maximum value of the normal stress transfers between the facial layer and the core layer. The maximum value of deflection and shear stress decreases with the reduction of the graded index and the interlayer thickness.



---

*Source: 1832525-2022-05-28.xml*