A 2-D numerical analysis for the structure composed by viscoelastic functionally graded materials using a temporally piecewise adaptive algorithm

Journal Pre-proof

A 2-D numerical analysis for the structure composed by viscoelastic functionally graded materials using a temporally piecewise adaptive algorithm Linlin Zhang , Haitian Yang PII: DOI: Reference:

S0307-904X(20)30015-9 https://doi.org/10.1016/j.apm.2020.01.015 APM 13228

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

21 February 2019 13 December 2019 8 January 2020

Please cite this article as: Linlin Zhang , Haitian Yang , A 2-D numerical analysis for the structure composed by viscoelastic functionally graded materials using a temporally piecewise adaptive algorithm, Applied Mathematical Modelling (2020), doi: https://doi.org/10.1016/j.apm.2020.01.015

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc.

Highlights 

A new adaptive algorithm with high accuracy in the time domain is put forward for the numerical analysis of VFG problems.



A piecewise variable separation is realized, providing a flexible option of numerical skills to solve spatial problems.



A stable computing accuracy with different step sizes can be achieved in the whole time domain via an adaptive computation.



Numerical verification is provided via homogeneous and regional inhomogeneous VFG structural analysis.

1

A 2-D numerical analysis for the structure composed by viscoelastic functionally graded materials using a temporally piecewise adaptive algorithm Linlin Zhang, Haitian Yang* Dept. of Engineering Mechanics State Key Lab of Structural Analysis for Industrial Equipment Dalian (116024), P. R. China

Abstract A temporally piecewise adaptive algorithm is developed for the 2-D numerical analysis of structures composed by viscoelastic functionally graded materials (VFGMs). By expanding all variables in a discretized time interval, in addition to the more precise description of variables, a piecewise variable separation is realized and a spatially and temporally coupled problem is decoupled into a series of recursive spatial problems which are solved by FEM. A stable computing accuracy can be achieved via a recursive adaptive process when the size of time steps varies. Numerical tests concerned with homogeneous and regional inhomogeneous VFG structures are given to verify the proposed approaches, and good accordance can be observed in comparison with analytical or ABAQUS based numerical solutions. Keywords: Viscoelastic functionally graded material; Series expansion; Piecewise variable separation; Recursive adaptive process; Finite element method.

1. Introduction Functionally graded materials (FGMs) are usually tailored to optimize the structural response under different types of loading conditions [1] or formed naturally under temperature variation and aging effect [2], and have been used in various engineering aspects including shipbuilding, automotive and aerospace industries, nuclear power plants and medical devices [3,4]. The characteristic feature of FGM is that its material properties can vary through certain dimensions in prescribed models [5,6]. Such a feature has been utilized to eliminate the stress singularity, improve the thermal properties, increase the fracture toughness and prevent stratification with the conventional composites [7]. Viscoelastic functionally graded materials (VFGMs) [8] can be considered as a member of FGMs family. In comparison with the conventional FGMs, the materials properties of VFGMs vary not only with the spatial position, but also with time, therefore the response of VFG structures will vary spatially and temporally, even if the load is static [9,10]. A number of investigations concerned with VFGMs have been conducted with various engineering backgrounds [7,11–28]. For instance, Loghman et al. investigated the magnetothermoelastic creep behavior of thick-walled FGM spheres [11]; Sofiyev studied the dynamic stability of VFG and heterogeneous orthotropic viscoelastic cylindrical shells [7,12]; Mao et al. studied the creep buckling of laminated piezoelectric VFG plates [13]; Norouzi et. al. conducted a static analysis of viscoelastic cylindrical panel made up of FGM [10]; Hamed studied bending and creep buckling response of VFG beam-columns [14].

2

The mechanical analysis of problems concerned with VFGMs is temporally and spatially coupled, and the analytical solution is usually only available in the case with regular geometries and simple boundary conditions. Paulino et al. gave an analytical stress solution for an infinite VFG strip by utilizing the correspondence principle [16]. Barretta et al. gave analytical solutions of a VFG Kirchhoff plate without boundary kinematic constraints and a VFG nanobeam subjected to a torsional loading [17,18]. By combining micro and macromechanical approaches, Tsukamoto analytically described the states of inelastic thermal-stress for a rectangle ceramic-metal FG plate subjected to a through-thickness heat flow [19]. Due to the complexities of constitutive relationships and boundary geometries/conditions, most VFG problems need to be solved numerically, leading to a demand for efficient numerical methods. In the spatial domain, VFG problems can be well modeled by the finite difference method, finite element method, boundary element method and meshless method [13,15,20]. While in the time domain, step-by-step integration and finite difference are major numerical schemes, such as predictor-corrector implicit-explicit time integration algorithm [21], Runge-Kutta approach [22], incremental time-stepping method [14], lower-order difference method [23], and Newmark time-marching scheme [13]. The solution inaccuracy of these algorithms may be caused by relatively lower truncation order or improper step size which is usually hard to predict in advance. The VFG problems can also be tackled using integral transformation method and correspondence principle [24], such as in the thermodynamic analysis for a coating made of VFGMs with interface cracks [25], in the numerical simulation for an exponentially graded viscoelastic bar [2], in the mechanical analysis of the FGM beams with viscoelastic interlayer [26], and in the fracture investigation on a VFG strip [27]. Although this kind of solution is not directly carried out in the time domain, numerical inversions are usually required for most cases, and may cause a computing inaccuracy in the time domain [28] because of the difficulty of analytic integral inversions. With the above consideration, a temporally piecewise adaptive algorithm is put forward in this paper for the numerical analysis of VFG problems. By virtue of a series expansion in a discretized time interval, the variation of all variables can be described more precisely, and a piecewise variable separation is realized. Consequently, a spatially and temporally coupled problem is decoupled into a series of recursive spatial problems which can be adaptively solved by FEM or other well developed spatial skills, and a stable computing accuracy can be maintained via a recursive adaptive process when the size of time steps varies. The rest of this paper is organized as follows. The governing and recursive governing equations for the numerical analysis are provided in Section 2. The constitutive and recursive constitutive equations are given in Section 3. The recursive FE equations and recursive adaptive process are described in Section 4. Numerical examples for the homogeneous and regional inhomogeneous VFG structures are presented in Section 5. Finally, conclusions are summarized in Section 6.

2. Governing and recursive governing equations The governing equations of a 2-D problem include [29]: Equilibrium equation

[ H ]{ }  {b}  {0}

3

(1)

Relationship of strain and displacement

{ }  [ H ]T {u}

(2)

 / y   / x 0 [H ]    / y  / x   0

(3)

where {σ}={σx, σy, τxy}T and {ε}={εx, εy, γxy}T denote the vectors of stress and strain, respectively; {b}={bx, by}T and {u}={ux, uy}T stand for the vectors of body force and displacement, respectively. The boundary conditions are specified by:

{u}  {u} on u

(4)

[n]{ }={f } on 

(5)

where {𝑢̃} is the prescribed value of {u} on the boundary Γu; {𝑓̃}={𝑓̃x, 𝑓̃y}T stands for the prescribed value of traction vector on the boundary Γσ; Γ =Γu +Γσ denotes the boundary of Ω; the subscript u and σ refer to displacement and stress, respectively; and [n] stands for the matrix of outside normal along Γσ. We divide the time domain into a number of time intervals. The initial points and sizes of time intervals are defined by t0, t1, t2, …, tk …and T1, T2 …, Tk …, respectively. At a discretized time interval, in order to describe the variation of variables more precisely, all variables are expanded in the term of s.

     m  s m

(6)

     m  s m

(7)

b   bm  s m

(8)

u    um 

sm

(9)

u   um  s m

(10)

{f }   {f m }s m

(11)

m0

m0

m0

m0

m0

m 0

s

t  tk 1 Tk

(12)

where tk−1 and Tk represent the initial point and size of the kth time interval, respectively; {σm}, {εm}, {bm}, {um}, {𝑢̃m} and {𝑓̃m} denote the expanding coefficients of {σ}, {ε}, {b}, {u}, {𝑢̃} and {𝑓̃}, respectively. Substituting Eqs. (6)–(11) into Eqs. (1), (2), (4) and (5) then yields [ H ]{ m }  {bm }  {0}

(13)

{ m }  [ H ]T {um }

(14)

4

{um }  {um } on u

(15)

[n]{ m }={f m } on 

(16)

3. Constitutive and recursive constitutive equations For VFGM, an integral form of constitutive equation can be described by [8] t deij (x, ) sij (x, t )  2  (x, t   ) ]d 0 d t d  (x, )  (x, t )  3  (x, t   ) d 0 d

(17) (18)

with sij   ij   ij , eij   ij   ij

(19)

where sij and eij refer to the tensors of deviatoric stress and strain, respectively; σij and εij denote the tensors of stress and strain, respectively; δij is the Kronecker delta; x=(x, y, z) stands for the coordinate vector; t refers to the time; 𝜎̅=σkk/3 and 𝜀̅=εkk/3 refer to the tensors of spherical stress and strain, respectively; μ(x, t) and κ(x, t) are shear and bulk relaxation functions [9]. 1−2𝑣

When Poisson’s ratio v is constant, we have μ(x, t) = 1+𝑣 κ(x, t). 𝜇(𝐱,𝑡)

Let Y(x, t) = 1+𝑣 , Eq. (17) and Eq. (18) can be merged in the form [9] 

 (x, t)  [G]  0 

t

Y (x, (t   ))

d  (x, )  d  , t  0 d 

(20)

where [G] is given by

G12  G11  G22 G     symm

0  0  G33 

(21)

For the plane stress problem G11  G22  1/ (1  v 2 ); G33  1/ (2(1  v)); G12  1/ (1  v)

(22)

For the plane strain problem G11  G22  (1  v) / ((1  v)(1  2v)); G33  1/ (2(1  v)); G12  v / ((1  v)(1  2v)

(23)

For the three-parameter solid viscoelastic model, as shown in Fig. 1, Eq. (20) becomes [9] 

 (x, t)  [G]  0 

t

E ( x )+E ( x )  E1 (x)  E2 (x) E (x) 1  ( x )2 (t  )  d  (x, )  (1  2 e ) d  , t  0   E1 (x) d   E1 (x)+E2 (x) 

(24)

where E1 (x), E2 (x), η(x) are constitutive parameters concerned with spatial positions. The differential form of Eq. (24) can be written as [9]

 (x, t )  p1 (x)

d  (x, t ) d  (x, t )    [G]  q0 ( x)  ( x, t ) +q1 ( x)  , dt dt  

5

t 0

(25)

 (x, t )  [G]E2 (x)  (x, t ) ,

t 0

(26)

where q0 (x) 

E1 (x)  E2 (x) E (x)  (x)  ( x) ; q1 (x)  2 ; p1 (x)  E1 (x)+E2 (x) E1 (x)+E2 (x) E1 (x)+E2 (x)

(27)

Using Tkds =dt and substituting Eqs. (6) and (7) into Eq. (25), then yield   }s m  p1 (x) m{ m }s m1  [G ]  q0 (x)  { m }s m +q1 ( x) m{ m }s m1  m0 m 1 m0 m 1   Equating the power of the two sides of Eq. (28) then gives

{

m

 m   [ D] m   Cm  ,

(28)

m  1, 2, 3,...,

(29)

 Tk  E1 (x) 1  m1   D m1   m   ( x) p1 (x) 

(30)

where

Cm  

 D11 [ D]    symm

D12 D22

0  0  D33 

(31)

For the plane stress problem, 1 1 v v ; D22  ; D33  ; D12  s2 s3 s2 s2

(32)

q1 (x) p ( x) p ( x) ; s  1  v 2  1 ; s3  2 1  v  1 p1 (x) 2 q1 (x) q1 (x)

(33)

D11  s1 

where

s1 

For the plane strain problem, E1 (x), E2 (x) and v need to be replaced by E1 (x)/(1–v2), E2 (x)/(1–v2) and v/(1–v2), respectively. At the first time inverval,

{ 0 }  [D]{ 0 }

(34)

At the (k+1)th time interval,

{ 0 }k(  1 )  { i }k  i 0 ,  {  }  {  }   0 k(  1 ) i k i 0 

k  1, 2 , 3 , . . . ,

(35)

where subscript (k+1) and k refer to (k+1)th and kth time intervals, respectively. In [15], following equations were presented to describe material properties of a VFGM

   0  f (x)

(36)

   0(1  1/  0  et / t )  f (x) 0

(37)

where κ0, μ0, μ1, t0 are constitutive parameters, and f(x) is the graded function. Eq. (36) implies that the relationship between bulk stress and strain is linear and time

6

independent at a fixed point (x). Eqs. (36–37) can be converted into a differential form similar to Eq. (25), but with different coefficients.

 (x, t )  t0

d  (x, t ) dt

 f (x)  [Q1 ] (x, t )  f (x)  [Q2 ]

d  (x, t )

(38)

dt

where 4 2   K 0  3 G0 K 0  3 G0  2 4 [Q1 ]   K 0  G0 K 0  G0  3 3  0 0  

 0   0  ,   G0  

4 4 2 2     K 0  3 G0  3 G1 K 0  3 G0  3 G1 0   2 2 4 4   [Q2 ]  K 0  G0  G1 K 0  G0  G1 0   3 3 3 3   0 0 G0  G1    

(39) A recursive constitutive equation similar to Eq. (29) can also be derived for Eq. (38).

4. Recursive FE equations The application of principle of virtual displacement on Eq. (13) with Eq. (15–16) gives [30]





 { }T { m }d     {u}T {bm }d     {u}T { f m }d 

(40)

where δ{ε} and δ{u} refer to vectors of virtual strain and virtual displacement, respectively. By virtue of FEM, {u(t)} and δ{u} are approximated in an element. {u (t )}=[N ]{ue (t )}

(41)

 {u}=[N ] {ue }

(42)

where [N] is a matrix of shape functions; {𝑢̅e(t)} and δ{𝑢̅e} are the nodal vectors of {u(t)} and δ{u} at the element level. Substituting Eq. (41) into Eq. (2) leads to { (t )}  [H ]T [N ]{ue (t )}  [ B]{ue (t )}

(43)

Expanding {𝑢̅e(t)} in the form

ue (t )   uem  s m

(44)

m0

Substituting Eq. (7) and Eq. (44) into Eq. (2) and equating the power of the two sides then yields { m }=[H ]T [N ]uem   [ B]uem 

(45)

Substituting Eq. (29) into Eq. (40) gives  {u }T [ Le ]

T

e



e

[ B]T ([ D] m   Cm )d e   {u }T [ Le ]T  [ N ]T {bm }d e   {u }T [ Le ]

T

e

e

e



e

[ N ]T { f m}d e

(46) where {u (t )}=  {um }s m and δ{𝑢̅} stand for the global nodal vectors of displacement and m 0

virtual displacement; [Le] refers to a mapping matrix from the global nodal vector of displacement to the nodal vector of displacement of an element, i.e. {𝑢̅e(t)}=[Le]{𝑢̅(t)}. Eq. (46) can further be written as

7

[ K ]{um }  Fm ,

m  0,1, 2,...,

(47)

where [ K ]  [ Le ]T  [ B]T [ D][ B]d e [ Le ] e

e

 T T T T  [ Le ] e [ N ] {b0 }d e   [ Le ] e [ N ] { f 0 }d  e , m  0 e  e  T T {Fm }   [ Le ]  [ B] {Cm }d e   [ Le ]T  [ N ]T {bm }d e e e e  e  T T m  1, 2, 3,   [ Le ] e [ N ] { f m }d  e ,  e

(48)

(49)

At the first time inverval, {𝑢̅0} is obtained by {u0 }  [ K ]1 F0

(50)

At other time intervals, {u0 }(k 1)  {ui }k

k  1, 2， 3, ...

(51)

i 0

Since [K] in Eq. (47) refers to an elastic stiffness mtrix and is positive definite symmetric, [K] can be decomposed in the form of LDLT. On the other hand, [K] keeps constant and needs to be generated only one time in the whole computing process, thus LDLT decomposition needs to carry out only one time, with which the computational cost on solving [K]{𝑢̅m}={Fm} can be efficiently reduced in the recursive process. For static loads,{bm}=0, {𝑓̃m}=0, and {Fm }  [Le ]T



e

e

[ B]T {Cm } d e

m  1, 2， 3, ...

(52)

{Fm}can be proved to be m

{Fm }  hl [ Kl ] {um l 1

[ Kl ]  [Le ]T e



l

}

(53)

 l 1[ B]T [ D][ B] d e [ Le ]

(54)

(Tk )(m  l )! E ( x) 1 ;  2 ;  m!  ( x) p1 (x)

(55)

e

where hl  

[Kl] needs to be computed only one time in the whole computing process. let Nk stand for the number of recursive computing at kth time interval and Nkmax stand for the maximum Nk. In the whole computing process, there are

N

k

linear equations needs to solve

k

and [Kl] only needs to generate Nmax times. In comparison with other step by step algorithms, a bit more computational efforts is required to deal with {𝑢̅m}=[K]-1Fm in a discretized time interval. It is a pay for the accuracy refinement. An adaptive computing process at each time interval can be realized with the increase of m until the following criteria is satisfied, and the expansion order will adaptively vary with different step sizes.

8

um  s m s 1

2

m 1

 u  s



(56)

j

j

j 0

s 1 2

where β is a prescribed error tolerance, and || ∙ ||2 represents a L2-norm. The truncation error of expansion series is controlled by β. When β becomes smaller, the truncation order will get higher and computing accuracy will increase. Eq. (56) is checked when each {𝑢̅m} is obtained. If it is satisfied, the computation at the current time interval will stop and step into the next. When the step size Tk becomes larger, hk becomes larger and more expanding items are usually required to satisfy the criterion of Eq. (56). At some stages of recursive computation, due to the round-off error in a computer, the values of some {𝑢̅j} cannot be described accurately with the increase of power order. The inaccuracy of

{u } i 0

i k

at kth time interval caused by such a

round-off error will propagate to the (k+1)th time interval as {𝑢̅0}(k+1) in Eq. (51), and result in the failure of computing. In order to avoid such a case, a regulating skill is provided as following: In the computation, mm and mim, the upper and lower bound of m, will be prescribed in advance. If Eq. (56) is not satisfied when m=mm, a size decrement of time step is necessary to restart the recursive procedure at the current time interval; if Eq. (56) is satisfied when m
5. Numerical examples In this section, five numerical examples are presented to verify the proposed algorithm, including three examples for homogeneous VFG structures and two for regional inhomogeneous VFG structures. The computing accuracy is evaluated by

j=

n

aij  aij

i

aij

(

)2 / n ,

j  x, y, r

(57)

where δj refers to a relative error; n designates the number of feature points; subscript i refers to the ith feature point; j refers to the jth direction of the coordinates; and aij and 𝑎̅ij stand for displacement in jth direction at ith feature point, given by the proposed algorithm and reference solutions, respectively. For Examples 1, 2, 4, 5, n refers to the total number of FE nodes overall the domain; for Example 3, n refers to the total number of FE nodes on the inner and outer boundaries, respectively. 5.1. Example 1 Consider a 1-D VFG strip made of Maxwell material with E2 (x) =E0e–αx and η(x) =η0. The strip is subjected to a time dependent displacement boundary condition [24] u(0, t )  0，u( L, t )=v0t

(58)

where u(x, t) stands for the displacement along the X direction; v0 is a constant; and L is the length of the strip. The analytical solution given by [24] is taken as the reference

9

u ( x, t ) 

0v0   x x v0 xt  L  bt  e  1  e   1 1  e    LE0  L L 

(59)

where b=αLE0/[η0(e–αL–1)]. E0, η0 and α are material parameters, and given by E0=1000(N/mm2), η0=1000(N∙s/mm2), α=– 1, v0=1(mm/s) and L=1(m). As shown in Fig. 1, Maxwell model can be considered as a specific case of three-parameter solid viscoelastic model when E1 (x)= 0 FE mesh consists of 100 two-node linear elements with 101 nodes. Table 1 and Table 2 show comparisons of displacement solutions at the middle point (x= 0.5m) with different β and step sizes, respectively. δx at different time is exhibited in Table 3. 5.2. Example 2 Consider a square VFG plate subjected to a uniform tension as shown in Fig. 2. The constitutive equation is described by Eqs. (25) and (26) with E1(x) =1000 e–0.5x(N/mm2), E2(x) =2000 e–x (N/mm2), η(x) =2000 (x+1) (N∙s/mm2), v=0.3. The FE mesh consists of 400 four-node elements with 441 nodes. Numerical results given by different step sizes are exhibited in Table 4 and Table 5, and compared with a converged ABAQUS solution with 1600 four-node elements with 1681 nodes (step size=0.001s). Fig. 3 exhibits the variation of the number of recursion with step sizes, and an increase of recursion with the increase of step sizes can be observed. Fig. 4 presents a comparison of displacements at Point A (1.0, 1.0) with a series of variant step sizes (0.1s, 0.5s, 1s, 3s, 1.5s, 0.9s, 2s, 1s) and a constant step size (1s) in the whole time domain, and Fig. 5 shows the variation of the number of recursion with the step size. 5.3 Example 3 Consider the example given in [15] where a thick hollow VFG cylinder subjected to an internal pressure is treated as a plane strain problem, as shown in Fig. 6. The constitutive equation is described by Eqs. (17), (18), (36) and (37) with κ0=1280N/mm2, μ0=120 N/mm2,μ1=360 N/mm2, t0=2.5 , f(x) = eγ(r -a), γ=0.1. A finite element model containing 480 four-node elements and 527 nodes is established. Table 6 gives a comparison of radial displacement provided by the proposed algorithm and the reference solution [15]. 5.4. Example 4 This example attempts to simulate a layered inhomogeneous VFG structure similar to the asphalt concrete pavements [28,31]. As shown in Fig. 7, the structure consists of four layers. From top to bottom, the first and third layers are composed of VFGMs, while the second and fourth ones are composed of elastic materials. The constitutive equation for the VFG layers is described by Eqs. (25) and (26), and material parameters of each layer are given in Table 7. FE mesh consists of 4250 four-node elements with 4386 nodes. Numerical results with different step sizes are compared with a converged ABAQUS solution with same FE mesh (step size=0.01s), as shown in Table 8. 5.5. Example 5 Consider a 1m×1m regional inhomogeneous VFG plate consisting of four parts, as shown in Fig. 8. Part 1 is elastic material and Part 2, 3 and 4 are VFGMs with different material parameters, as listed in Table 9. The constitutive equation for the VFG parts is described by Eqs. (25) and (26). A FE model containing 10000 four-node elements with 10201 nodes is employed. A comparison between the results given by the proposed algorithm and a converged ABAQUS solution (40000

10

four-node elements with 40401 nodes, step size=0.1s) is given in Table 10. Using ABAQUS solution (step size=0.1s) as a reference, Fig. 9 gives a comparison of relative errors of displacement solutions at Point B (1.0, 1.0), provided by the proposed algorithm and ABAQUS with step sizes=10s, 20s, respectively. In addition to the numerical comparison, a comparison of 2-D contour plots of X-displacement overall the field is presented in Fig. 10, visually displaying a good accordance. Computing Remarks (1) As shown in Tables 1–6, 8 and 10, the solutions given by the proposed algorithm have good accordance with both analytical solutions and numerical solutions obtained by ABAQUS. The maximum δj is 5.1683%. (2) As shown in Tables 3–6, 8, 10, and Fig. 4, for a prescribed error tolerance β, a stable computing accuracy can be achieved via an adaptive procedure for different step sizes (variant or constant in the whole time domain). In contrast, the computing accuracy of ABAQUS is more or less affected by the step size, and the maximum deviation reaches 14% when the step size is 20s and almost 9% when the step size is 10s, as shown in Fig. 9. (3) The all numerical solutions provided by the proposed approach are mesh converged. In this sense, the source of error can mainly be attributed to the truncation error of expansion series which is controlled by β. When β becomes smaller, the truncation order will get higher and computing accuracy will increase, as exhibited in Tables 1 and 3 in comparison with the analytical solution. In order to satisfy the criteria of Eq. (57) with a prescribed β, an adaptive process is carried out via the increase of expansion order which will adaptively vary with different step sizes. Eq. (57) may be regarded as a β based convergence criterion of the expansion series, which is proved to be available in all the numerical tests. Frankly, it seems difficult to present a rigorous and general convergence criterion with a proper error tolerance because of complexity of recursive equations. The unbiased estimator, which was successfully utilized for the determination of regularization parameter in [32], is heuristic for establishing a more rigorous convergence criterion of recursive equations, although it seems a bit difficult to realize at the present stage. More efforts in this side will be made in the succeeding work.

6. Conclusions The major objective of this paper is to develop a numerical method with high accuracy for the numerical analysis of structures composed by VFGMs. The major merits of the proposed approach include: 1) Variables are described more accurately via piecewise series expansion in discretized time intervals. 2) A piecewise variable separation is realized via the above expansion, and a spatially and temporally coupled problem is decoupled into a series of spatial problems which can be numerically solved either by FEM or other well developed numerical methods in the space domain. 3) The computing accuracy is emphasized by a recursive adaptive process which is carried out via the change of expansion order with a prescribed error tolerance, and a stable solution accuracy can be achieved with constant or variant step sizes in the whole time domain. The difficulty of choosing step size can be tackled with an adjustment of step size and adaptive computation. 4) In comparison with analytical or ABQUS solutions, the proposed algorithm is verified via

11

the numerical analysis for both homogeneous and regional inhomogeneous VFG structures with good accordance.

Acknowledgements The research leading to this paper is funded by NSF [11572068, 11972109], NKBRSF [2015CB057804], NSF [11572077]. A special and deep appreciation is given to Mr. Chunjiang Ran for his suggestions and efforts in the process of manuscript revision.

References [1] M. Koizumi, FGM activities in Japan, Compos. Part B 28 (1997) 1–4. [2] E.V. Dave, G.H. Paulino, W.G. Buttlar, Viscoelastic functionally graded finite element method using correspondence principle, J. Mater. Civil Eng. 23 (2011) 39–48. [3] X. Wang, Z. Yuan, C. Jin, 3D free vibration analysis of multi-directional FGM parallelepipeds using the quadrature element method, Appl. Math. Model. 68 (2019) 383–404. [4] V.N. Van Do, C.H. Lee, Nonlinear analyses of FGM plates in bending by using a modified radial point interpolation mesh-free method, Appl. Math. Model. 57 (2018) 1–20 [5] M. Čanađija, R. Barretta, F.M. de Sciarra. On functionally graded Timoshenko nonisothermal nanobeams, Compos. Struct. 135 (2016) 286–296. [6] R. Barretta, S.A. Faghidian, R. Luciano, Free vibrations of FG elastic Timoshenko nano-beams by strain gradient and stress-driven nonlocal models, Compos. Part B 154 (2018) 20–32. [7] A.H. Sofiyev, About an approach to the determination of the critical time of viscoelastic functionally graded cylindrical shells, Compos. Part B 156 (2019) 156–165. [8] G.H. Paulino, Z.H. Jin, Viscoelastic functionally graded materials subjected to antiplane shear fracture, J. Appl. Mech. 68 (2) (2001) 284–293. [9] I. H. Shames, Elastic and inelastic stress analysis, CRC Press, (1997). [10] H. Norouzi, A. Alibeigloo, Three dimensional static analysis of viscoelastic FGM cylindrical panel using state space differential quadrature method, Euro. J. Mech. A/Solids 61 (2017) 254– 266. [11] A. Loghman, S.M.A. Aleayoub, M.H. Sadi, Time-dependent magnetothermoelastic creep modeling of FGM spheres using method of successive elastic solution, Appl. Math. Model. 36(2) (2012) 836–845. [12] A.H. Sofiyev, On the solution of the dynamic stability of heterogeneous orthotropic visco-elastic cylindrical shells, Compos. Struct. 206 (2018) 124–130. [13] Y.Q. Mao, Y.M. Fu, H.L. Dai, Creep buckling and post-buckling analysis of the laminated piezoelectric viscoelastic functionally graded plates, Euro. J. Mech. A/Solids 30 (2011) 547–558. [14] E. Hamed, Bending and creep buckling response of viscoelastic functionally graded beam-columns, Compos. Struct. 94 (10) (2012) 3043–3051. [15] S.S. Chen, C.J. Xu, G.S. Tong, A meshless local natural neighbour interpolation method to modeling of functionally graded viscoelastic materials, Eng. Anal. Bound. Elem. 52 (2015) 92–98. [16] G.H. Paulino, Z. Jin, Correspondence principle in viscoelastic functionally graded materials, J. Appl. Mech. 681 (2001) 129–132. [17] R. Barretta, L. Feo, R. Luciano, Torsion of functionally graded nonlocal viscoelastic circular nanobeams, Compos. Part B 72 (2015) 217–222. [18] R. Barretta, R. Luciano, Exact solutions of isotropic viscoelastic functionally graded

12

Kirchhoffplates, Compos. Struct. 118 (2014) 448–454. [19] H. Tsukamoto, Analytical method of inelastic thermal stresses in a functionally graded material plate by a combination of micro and macromechanical approaches. Compos. Part B 34 (2003) 561–568. [20] F. Pan, W. Li, B. Wang, X. Zhang, Viscoelastic fracture of multiple cracks in functionally graded materials, Comput. Methods Appl. Mech. Eng. 198 (33) (2009) 2643–2649. [21] M.A. Fahmy, Implicit–explicit time integration DRBEM for generalized magneto-thermoelasticity problems of rotating anisotropic viscoelastic functionally graded solids, Eng. Anal. Bound. Elem. 37 (1) (2013) 107–115. [22] H. Hilton, Tailored designer functionally graded materials for minimizing probabilistic creep buckling failures in linear viscoelastic columns with large deformations and follower loads, Aiaa/asme/asce/ahs/asc Structures, Structural Dynamics, and Materials Conference, Aiaa/asme/ahs Adaptive Structures Conference. (2013). [23] H. Ashrafi, M. Shariyat, S.M.R. Khalili, K. Asemi，A boundary element formulation for the heterogeneous functionally graded viscoelastic structures, Appl. Math. Comput. 225 (2013) 246– 262. [24] S. Mukherjee, G.H. Paulino, The Elastic-Viscoelastic Correspondence Principle for Functionally Graded Materials, Revisited, J. Appl. Mech. 70 (3) (2003) 359–363. [25] Z.Q. Cheng, S.A. Meguid, Z. Zhong, Thermo-mechanical behavior of a viscoelastic FGMs coating containing an interface crack, Int. J. Fracture 164 (2010) 15–29. [26] J. Li, B. Zheng, Q. Yang, X.J. Hu, Analysis on time-dependent behavior of laminated functionally graded beams with viscoelastic interlayer, Compos. Struct. 107 (1) (2014) 30–35. [27] Z.H. Wang, L. Zhang, L.C. Guo, A viscoelastic fracture mechanics model for a functionally graded materials strip with general mechanical properties, Euro. J. Mech. A/Solids 44 (1) (2014) 75–81. [28] E.V. Dave, Asphalt pavement aging and temperature dependent properties using functionally graded viscoelastic model, Dissertations & Theses - Gradworks, 631–632 (2009) 190. [29] R.M. Christensen, Theory of Viscoelasticity, Academic Press. New York. (1982). [30] M.P. Williamson, Finite-element analysis, Computer-Aided Engineering Journal, 1985, 2(2):66. [31] J. Liu, Y. Sun, W. Wang, J. Chen, Using the viscoelastic parameters to estimate the glass transition temperature of asphalt binders, Constr. Build. Mater. 153 (2017) 908–917. [32] S.A. Faghidian, Inverse determination of the regularized residual stress and eigenstrain fields due to surface peening, Strain. Anal. Eng. 50(2) (2015) 84–91.

13

Fig. 1. Three-parameter solid viscoelastic model

14

Fig. 2. A square VFG plate

15

16

step size= 1.0s step size= 0.1s step size= 0.01s

Number of recursion

14 12 10 8 6 4 0

2

4

6

8

10

Time (s)

Fig. 3. Variation of the number of recursion with time

16

0.20

constant step variant step

X-displacement (mm)

0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0

2

4

6

8

10

Time (s)

(a) X-displacement -0.01

constant step variant step

Y-displacement (mm)

-0.02

-0.03

-0.04

-0.05

-0.06 0

2

4

6

8

10

Time (s)

(b) Y-displacement Fig. 4. Comparisons of (a) X-displacement and (b) Y-displacement at Point A (1.0, 1.0) with variant and constant steps

17

30

constant step variable step

Number of recursion

25

20

15

10

5 0

2

4

6

8

10

Time (s)

Fig. 5. Variation of the number of recursion with variant and constant steps

18

Fig. 6. A thick hollow VFG cylinder

19

Fig. 7. A layered inhomogeneous VFG structure

20

Fig. 8. A regional inhomogeneous VFG plate

21

16 ABAQUS step size=10s ABAQUS step size=20s Proposed method step size=10s Proposed method step size=20s

14

Relative Error (%)

12 10 8 6 4 2 0 0

200

400

600

800

1000

Time (s)

(a)

12

ABAQUS step size=10s ABAQUS step size=20s Proposed method step size=10s Proposed method step size=20s

Relative Error (%)

10 8 6 4 2 0 0

200

400

600

800

1000

Time (s)

(b) Fig. 9. Variation of relative errors with time at Point B (1.0, 1.0) (a) X-displacement (b) Y-displacement

22

t=0s

t = 100 s

t=0s

VAL (mm)

VAL (mm)

0.03610

0.03610

0.03159

0.03158

0.02708

0.02705

0.02256

0.02253

0.01805

0.01800

0.01354

0.01348

0.009025

0.008950

0.004513

0.004425

0.000

-1.000E-04

VAL (mm)

0.1270

t = 100 s

VAL (mm) 0.1265

0.1111 0.1106

0.09525 0.09475

0.07938

0.07888

0.06350

t = 500 s

0.06300

0.04763

0.04713

0.03175

0.03125

0.01588

0.01538

0.000

-5.000E-04

t = 500 s

VAL (mm)

0.2800

0.2449

0.2449

0.2098

0.2098

0.1746

0.1746

0.1395

0.1395

0.1044

0.1044

0.06925

0.06925

0.03413

0.03413

-0.001000

-0.001000

t = 3000 s t = 3000 s

VAL (mm)

0.2800

VAL (mm)

VAL (mm) 0.4360

0.4360

0.3812 0.3812

0.3265

0.3265

0.2717

0.2717

0.2170

0.2170

0.1623

0.1623

0.1075

0.1075

0.05275

0.05275

-0.002000

-0.002000

(a) The proposed method

(b) ABAQUS

Fig. 10. Comparisons of 2-D contour plots of X-displacement overall field

23

Table 1 Relative error of displacement solution at the middle point with different β (Tk=1s) Time (s) 10 20 30 40 50 60 70 80 90 100

Present method (10-3m) -2

β=10

β=10

Error (%)

4.95972 9.95967 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

-3

1.07×10 1.12×10-5 1.01×10-6 6.78×10-8 1.59×10-8 4.65×10-9 3.94×10-9 3.44×10-9 3.05×10-9 2.72×10-9

-8

4.95966 9.95966 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

Error (%)

Analytical (10-3m)

4.19×10-7 3.59×10-9 1.39×10-10 2.25×10-10 2.84×10-10 3.38×10-10 3.98×10-10 4.58×10-10 5.13×10-10 5.70×10-10

4.95966 9.95966 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

Table 2 Relative error of displacement solution at middle point with different step sizes (β=10-8) Time (s) 10 20 30 40 50 60 70 80 90 100

Present method (10-3m) Tk=1s

Error (%)

Tk=5s

-7

4.19×10 3.59×10-9 1.39×10-10 2.25×10-10 2.84×10-10 3.38×10-10 3.98×10-10 4.58×10-10 5.13×10-10 5.70×10-10

4.95966 9.95966 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

4.95966 9.95966 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

Error (%)

Analytical (10-3m)

6.08×10-7 2.36×10-8 2.25×10-9 3.36×10-9 2.24×10-9 1.72×10-9 3.42×10-9 3.00×10-9 2.63×10-9 2.47×10-9

4.95966 9.95966 14.95966 19.95966 24.95966 29.95966 34.95966 39.95966 44.95966 49.95966

Table 3 δx with different step sizes and β Tk=1s

Tk=5s

Time (s)

β=10

β=10

10 20 30 40 50 60 70 80 90 100

6.61894×10-3 4.26087×10-5 5.95492×10-7 1.38762×10-7 8.34487×10-8 5.58845×10-8 3.97442×10-8 2.94924×10-8 2.25990×10-8 1.77338×10-8

9.92703×10-9 5.58258×10-11 1.62830×10-10 2.15014×10-10 2.64576×10-10 3.21938×10-10 3.74083×10-10 4.30787×10-10 4.83916×10-10 5.36727×10-10

-2

-8

24

-2

β=10

β=10-8

6.61894×10-3 4.26087×10-5 5.95492×10-7 1.38762×10-7 8.34487×10-8 5.58845×10-8 3.97442×10-8 2.94924×10-8 2.25990×10-8 1.77338×10-8

9.92703×10-9 5.58258×10-11 1.62830×10-10 2.15014×10-10 2.64576×10-10 3.21938×10-10 3.74083×10-10 4.30787×10-10 4.83916×10-10 5.36727×10-10

Table 4 δx and δy with different β (Tk=1s) δy

δx

Time (s)

β=10-2

β=10-8

β=10-2

β=10-8

1 2 3 4 5 6 7 8 9 10

0.293428% 0.198219% 0.166954% 0.194043% 0.167275% 0.149882% 0.177398% 0.100377% 0.218399% 0.312709%

0.017708% 0.006545% 0.003098% 0.002281% 0.002322% 0.002491% 0.002629% 0.002722% 0.002789% 0.002828%

0.583583% 0.375390% 0.420076% 0.336807% 0.689485% 1.253350% 1.360468% 0.665097% 1.302702% 2.235392%

0.018533% 0.008152% 0.004049% 0.002810% 0.003242% 0.003992% 0.004624% 0.005106% 0.005445% 0.005686%

Table 5 δx and δy with different β (Tk=0.1s) β=10-8 δx

Time (s)

β=10

1 2 3 4 5 6 7 8 9 10

0.020527% 0.019922% 0.055385% 0.132064% 0.120679% 0.093497% 0.068641% 0.049391% 0.035287% 0.025148%

-2

δy -8

-2

β=10

β=10

β=10-8

0.017708% 0.006545% 0.003098% 0.002281% 0.002322% 0.002491% 0.002629% 0.002722% 0.002789% 0.002828%

0.011982% 0.012334% 0.061981% 0.159806% 0.164156% 0.141168% 0.113127% 0.087371% 0.065919% 0.048869%

0.018533% 0.008152% 0.004049% 0.002810% 0.003242% 0.003992% 0.004624% 0.005106% 0.005445% 0.005686%

Table 6 δr with different step sizes (β=10-8) [15] δr Time(s) 0.01 2.64105 5.2721 7.90316 10.53421 13.16526 15.79632 18.42737 21.05842

Tk=0.01s

Tk=0.1s

0.6545% 0.7677% 0.6887% 0.7343% 0.7546% 0.8514% 0.8242% 0.8979% 0.7673%

0.6545% 0.7677% 0.6887% 0.7343% 0.7546% 0.8514% 0.8242% 0.8979% 0.7673%

25

23.68947 26.32053 28.95158 31.58263 34.21368 36.84474 39.47579 42.10684 44.7379 47.36895 50

0.8902% 0.8598% 0.8587% 0.8304% 0.9393% 0.9432% 0.8788% 0.9176% 0.9296% 0.9724% 0.8945%

0.8902% 0.8598% 0.8587% 0.8304% 0.9393% 0.9432% 0.8788% 0.9176% 0.9296% 0.9724% 0.8945%

Table 7 Material parameters Layer (from the top down)

1 2 3 4

Thickness (m) 0.2 0.3 0.4 0.8

Material

E1(x) ( N/mm2)

E2(x) ( N/mm2)

η(x) ( N·s /mm2)

v

257×y2 / 138×ey /

1355×ey 700 1219×e2y 1800

39876×eyy-2 / 21194×e-yy2 /

0.3 0.3 0.3 0.3

VFGM elastic material VFGM elastic material

Table 8 δy with different step sizes (β=10-8) Time(s) 0 30 60 90 120 150 180 210 240 270 300

δy Tk=0.1s

Tk=0.01s

0.001372% 0.00616% 0.00898% 0.00922% 0.00918% 0.00923% 0.00925% 0.00925% 0.00926% 0.00927% 0.00928%

0.001372% 0.00616% 0.00898% 0.00922% 0.00918% 0.00923% 0.00925% 0.00925% 0.00926% 0.00927% 0.00928%

26

Table 9 Material parameters

Part

Material

E1(x) ( N/mm2)

E2(x) ( N/mm2)

η(x) ( N·s /mm2)

v

1 2 3 4

VFGM VFGM elastic material VFGM

179× e-0.5y 85/x / 218× e-0.5y

1226× e-y 1355×(x+1) -1 1123 2476×e-y

2673×(y+1) 39875×ex / 27528×(y+1)

0.3 0.3 0.3 0.3

Table 10 δx and δy with different step sizes (β=10-8) Time(s) 0 300 600 900 1200 1500 1800 2100 2400 2700 3000

δx

δy

Tk=0.1s

Tk=1s

Tk=0.1s

Tk=1s

1.2472% 1.1508% 1.5734% 1.1451% 2.4670% 2.0862% 1.4283% 1.0397% 1.0194% 1.2413% 1.4729%

1.2472% 1.1508% 1.5734% 1.1451% 2.4670% 2.0862% 1.4283% 1.0397% 1.0194% 1.2413% 1.4729%

0.0719% 0.0548% 0.0527% 0.0485% 0.0471% 0.0469% 0.0468% 0.0468% 0.0468% 0.0468% 0.0468%

0.0719% 0.0548% 0.0527% 0.0485% 0.0471% 0.0469% 0.0468% 0.0468% 0.0468% 0.0468% 0.0468%

27