A three-scale asymptotic analysis for ageing linear viscoelastic problems of composites with multiple configurations

Applied Mathematical Modelling 71 (2019) 223–242

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A three-scale asymptotic analysis for ageing linear viscoelastic problems of composites with multiple configurations Zhiqiang Yang a,∗, Yi Sun a, Yizhi Liu a, Tianyu Guan a, Hao Dong b a b

Department of Astronautic Science and Mechanics, Harbin Institute of Technology, Harbin 150001 China School of Mathematics and Statistics, Xidian University, Xi’an 710071 China

a r t i c l e

i n f o

Article history: Received 10 May 2018 Revised 28 January 2019 Accepted 13 February 2019 Available online 23 February 2019 Keywords: Three-scale asymptotic expansion Homogenization Ageing linear viscoelastic Multiple configurations

a b s t r a c t A novel three-scale asymptotic expansion is proposed in this work to investigate ageing linear viscoelastic problems of composites with multiple periodic configurations. Here, the composite structures are established by periodical distribution of local cells on microscale and mesoscale. The new three-scale asymptotic expansion formulas based on classical homogenized methods in time domain are constructed at first, and the microscale and mesoscale functions are also derived in detail. Further, two distinct homogenized parameters defined on microscale and mesoscale domains are obtained by upscaling methods, and the newly developed homogenized equations are given on whole structure in time domain. Then, the three-scale strain and stress solutions are constructed by assembling the unit cell solutions and homogenized solutions. Also, the efficient finite element algorithms based on the three-scale asymptotic expansion and homogenized method are brought forward. Finally, some representative numerical examples are evaluated to verify the proposed methods. They show that the three-scale asymptotic expansions developed in this work are effective and valid for predicting the ageing linear viscoelastic properties of the composite materials with multiple configurations. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Composite materials with viscoelastic properties have wide utilization in actual engineering and industry. In polymer materials, damping properties can be analyzed in connection with strength properties. In concrete composites, investigating creep magnitude can permit reduction of the related damage [1]. These types of composites are usual heterogeneous at certain scale [2], and the macroscopic properties of the composites are controlled by a fine-scale structure. In addition, several biological and synthetic viscoelastic materials often possess multiscale character and exhibit more than two length scales. The associated experimental relaxation needs extremely considerable expense and can last for months or years. Inevitably, in order to design a high-performance concrete, it is very necessary to find an effective prediction model and numerical algorithm to accurately describe the thermal and mechanical behavior for the composite structures. Generally speaking, traditional homogenization methods analyze the macroscopic effective behavior of composites by making the original equations with rapidly varying coefficients to the homogenized equations with effective parameters, which can not only save the computation time but also maintain the calculation precision [3–11]. Based on the classical

∗

Corresponding author. E-mail address: [email protected] (Z. Yang).

https://doi.org/10.1016/j.apm.2019.02.021 0307-904X/© 2019 Elsevier Inc. All rights reserved.

224

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

homogenization methods [4–6], different kinds of multiscale methods for composites with periodic distribution were developed, referring to the heterogeneous multiscale method [12,13], two-scale convergence method [14], multiscale finite element method [15,16] and variational multiscale method [17,18]. However, these multiscale methods just expanded the first-order asymptotic expansions for the periodical problems. In some cases [9,19–22], the first-order multiscale asymptotic solutions can not accurately obtain the local oscillations behaviour of composites for some physical and mechanical fields. To solve the difficulties, Allaire and Habibi [23] and Cuiand co-workers [19–22] constructed a higher-order multiscale asymptotic solution by introducing some suitable high order unit cell functions, and the local thermal and mechanical information inside the composites can be captured effectively. In addition, for the viscoelastic composites with periodic microstructures, the homogenized methods and the associated multiscale methods were extended to investigate the problems in the past years. Yi et al. [24] utilized the classical homogenized methods to analyze the effective macroscopic properties for the viscoelastic composites. Liu et al. [25], Cai and Sun [26] developed the homogenization methods to predict the viscoelastic properties of layered materials and three-dimensional braided composites, respectively. Tran et al. [27] employed a computational homogenization method for computing the linearly non-aging viscoelastic and highly heterogeneous materials. Zhang et al. [28] built a statistical higher-order multiscale model to simulate the viscoelastic problems of composites with random distribution of particles. Later, Nguyen et al. [29] established a micromechanical approach based on asymptotic expansion methods to study the viscoelastic model of porous materials. Also, Sanahuja [30] and Maghous and Creus [31] developed the homogenization methods and multiscale methods in the study of ageing problems and thermosviscoelasticity problems in heterogeneous materials. Lavergnea et al. [32] proposed an extension of the Mori-Tanaka and Ponte Castañeda-Willis homogenization technique for the aging linear viscoelastic composites with ellipsoidal inclusions. Zhang et al. [33] presented a second-order two-scale computation in time domain to effectively predict the aging linear viscoelastic properties of composite structures with a periodical distribution. Obviously, from the work mentioned previously, the homogenized methods and the relevant higher-order multiscale algorithms can evaluate accurate prediction of the viscoelastic and thermos-viscoelasticity properties of composites with complex microstructures. Besides, such numerical techniques and algorithms can perform effective analysis of the strain and stress fields on the macroscopic scale according to the effective parameters obtained at the mesoscopic scale. In addition, some composite structures often have different length scales in practical engineering computations, such as hierarchical layered materials including lotus leaves, kidney’s glomeruli, bones and concrete materials, the classical homogenization methods and the higher-order two-scale methods are not efficient to obtain the macroscopic thermal and mechanical behavior for the composites. To resolve these shortcomings, a reiterated homogenization method was constructed by Bensoussan et al. [4] to effectively consider the micro–meso–macro scales of heterogeneities. In theory, Allaire and Briane [34], Holmbom et al. [35], Trucu et al. [36] and Abdulle and Bai [37] analyze the three-scale convergence for the distinct multiscale problems. Fabricius and co-workers [38] extended the reiterated homogenization method to study the surface effect in hydrodynamic lubrication. Liu and co-workers [39,40] and Yan and co-workers [41] established a statistic threescale model on the basis of traditional homogenized methods to simulate the macroscopic effective coefficients for fiber reinforced concrete. Cruz and co-workers [42,43] applied the reiterated homogenization methods to predict the macroscopic performance of heat conduction problems in inhomogeneous materials with multiple configurations. Cao [44] considered a linear elastic problem with various length scales by an iterated two-scale homogenization method, and proposed a detailed mathematical verification for some convergence analysis. from Cao’ work [44], the iterated two-scale homogenization method is very efficient to acquire the macroscopic effective parameters from microscale domain to macroscale domain step by step. Mahnken and Dammann [45,46] developed a three-scale homogenization model for temperature-dependent viscoelastic effect accompanied by curing in a resin transfer molding process, and effective macroscopic coefficients are given by the homogenization model for a representative unit cell on the heterogeneous microscale and mesoscale. RamírezTorres et al. [47] used multiscale asymptotic methods to investigate the effective properties of hierarchical composites with periodic microstructure at distinct length scales. Unfortunately, for these types of structures with multiple configurations, the two-scale methods and reiterated homogenization can be utilized to analyze the macroscopic effective coefficients from smallest scale to largest scale step by step, but they cannot efficiently investigate their thermal and mechanical behaviors. In other words, the local stress and strain distributions of the composites cannot be more accurately captured by these methods. Thus, motivated by the multiscale homogenization method [4–6] and the work in [39], Cui and co-workers [48,49] proposed a higher-order three-scale asymptotic approach for the elastic problems of heterogeneous materials with multiple spatial scales. Later, Yang et al. [50] developed the three-scale method to study the coupled conduction-radiation problems of porous structures with periodic configuration. Obviously, by higher-order correctors and upscaling techniques, the microscale and mesoscale information of thermal and mechanical behaviors inside the composites can be caught more effectively. However, the previous higher-order multiscale analysis and algorithms cannot be directly used to analyze the viscoelastic composites owing to the complicated time integration for the problems and the specific microstructure of the composites. Thus, this work will develop a new high-order asymptotic expansion that is completely different from the multiscale algorithms reported in Yang et al. [50] to investigate the aging viscoelastic problems including the new homogenized equations and new high-order unit cell functions. Further, to handle the ageing linear viscoelastic problems in this work, an efficient three-scale numerical algorithm and technique considering the microscale and mesoscale of the displacement solutions is needed to be established in time domain directly. In addition, to our knowledge, the correspondence principle [51] is an important mathematical technique to evaluate the non-ageing linear viscoelastic problems since it permits to use the numerical algorithms originally developed for the

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

225

Fig. 1. Illustration of the three-scale composite structure.

elastic case. Unfortunately, the correspondence principle will be invalid as the composites emerge ageing effect. Inevitably, the previous three-scale analysis and numerical technique cannot be effectively applied to the aging linear viscoelastic problems due to the complicated integral form in time domain for the original equations. The reason is that we could not directly establish the homogenized equations and the auxiliary cell equations by the multiscale methods in Laplace domain. Moreover, the microscale and mesoscale information of the viscoelastic problems cannot be accurately obtained by the classical three-scale methods. But all in all, there is a lack of enough researches on ageing linear viscoelastic problems of heterogeneous materials with multiple configurations. Also, we have not found effective high-order multiscale methods to discuss the aging linear viscoelastic problems of composite structures with multiple spatial scales. Therfore, in this work, we will mainly investigate the ageing linear viscoelastic properties of composite materials with multiple small-scale configurations operating in time domain directly, as shown in Fig. 1. The aim of this study is devoted to designing a newly micro–meso–macro multiscale algorithm to give a better approximation. Moreover, the CPU time and computer memory can be saved drastically. Also, it should be pointed out that the higher-order methods given in this work include the efficient high-order unit cell functions compared to traditional three-scale methods in [4,34–47], which can make the proposed approaches accurately catch the microscopic oscillation behaviors and remarkably change the numerical precision for solving these multiscale problems. This paper is organized as follows. In the next section, the detailed governing equations of aging linear viscoelastic problems are given. In Section 3, the three-scale asymptotic expansion formulas for the aging viscoelastic problems are derived. Section 4 illustrates the three-scale finite element procedures in detail. In Section 5, some typical examples are proposed to confirm the accuracy and effectiveness of the three-scale algorithms. Section 6 is devoted to some conclusions. 2. Aging linear viscoelastic models in composites with multiple configurations In this study, composite structures are characterized by three different length scales or two different periodical unit cell, as illustrated in Fig. 1. The macroscopic effective properties of the aging viscoelastic composites at microscale and mesoscale are investigated by two small local cells: the microscale cell ε 2 Z with size ε 2 and the mesoscale cell ε 1 Y with size ε 1 , and ε 2 < <ε 1 . The microscopic domain ε 2 Z contains two or more material phases and establishes the matrix on the mesoscopic scale ε 1 Y. Then, the macroscopic domain  is constructed by mesoscopic domain ε 1 Y and microscopic domain ε 2 Z. According to the detailed microstructure description, the aging linear viscoelastic problems for a composite structure with multiple spatial scales are given as follows:

−

 ∂  ε1 ε2 σ (uε1 ε2 (x, t )) = fi (x, t ), in  × (0, T ), i = 1, 2, 3, ∂xj ij

σiεj1 ε2 (uε1 ε2 (x, t )) =



t

0

1 ε2 Gεi jkl (x, t, τ )

∂ eεkl1 ε2 (uε1 ε2 (x, τ )) dτ , ∂τ

(1) (2) ε ε

1 2 where uε1 ε2 (x, t ) = (uε1 ε2 1 (x, t ), uε1 ε2 2 (x, t ), uε1 ε2 3 (x, t ))T are the displacement fields; Gi jkl (x, t, τ )(i, j, k, l = 1, 2, 3) denote ε ε T the coefficients of viscoelastic relaxation tensor; f(x, t) = (f1 (x, t), f2 (x, t), f3 (x, t)) is the body force; ekl1 2 (uε1 ε2 (x, τ )) (k, l = 1, 2, 3 ) denotes the strain computed from displacement uε1 ε2 (x, τ ) given by

ε1 ε2

ekl

1 (uε1 ε2 (x, τ )) = 2



 ∂ u ε1 ε2 k ( x , τ ) ∂ u ε1 ε2 l ( x , τ ) + . ∂ xl ∂ xk

226

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

The boundary conditions are specified as follows:

v j σiεj1 ε2 (x, t ) = pi (x, t ), on 1 × (0, T ), uε1 ε2 (x, t ) = u¯ (x, t ),

on 2 × (0, T ),

u ε1 ε2 ( x , 0 ) = 0,

x ∈ ,

(3)

where v = (vj ), j = 1, 2, 3 denotes the unit outward normal on  1 , pi (x, t) is the tractions on  1 .  1 and  2 define the parts of the boundary domain, and satisfy  1 ∪ 2 = ∂ ,  1 ∩ 2 = ∅, as illustrated in Fig. 1. u¯ (x, t ) denotes the value of a displacement ε y component on  2 . Let y = εx and z = εx = ε1 , x denotes macroscale coordinate of the structure, y is mesoscale coordinate 1

2

2

ε ε

1 2 in ε 1 Y and z is microscale coordinate in ε 2 Z. Then, relaxation moduli Gi jkl (x, t, τ ) = Gi jkl ( εx1 , εx2 , t, τ ) = Gi jkl (y, z, t, τ ) are Y × Z-periodic piecewise differential functions, and satisfy the conditions of symmetry 1 ε2 1 ε2 1 ε2 Gεi jkl (x, t, τ ) = Gεi jlk (x, t, τ ) = Gεjikl (x, t, τ ),

and ellipticity 1 ε2 λ1 ηi j ηi j ≤ Gεi jkl (x, t, τ )ηi j ηkl ≤ λ2 ηkl ηkl ,

where (ηij ) is an arbitrary real symmetric matrix in R3×3 , λ1 and λ2 are the positive constants independent of ε 1 and ε 2 . Also, the existence and uniqueness of (1)–(3) were investigated by Orlik [52,53]. 3. Three-scale asymptotic expansion formulas In this section, the three-scale asymptotic expansion formulas are given in detail for computing the aging linear viscoelastic problems of heterogeneous materials with multiple configuration. Then, enlightened by the idea proposed in [4–6,48,49], uε1 ε2 (x, t ) can be obtained into the series in the following form:

uε1 ε2 (x, t ) = u0 (x, t ) + ε1 u1 (x, y, t ) + ε12 u2 (x, y, t ) + ε2 u3 (x, y, z, t ) + ε1 ε2 u4 (x, y, z, t ) + ε22 u5 (x, y, z, t ) + ε12 ε2 u6 (x, y, z, t ) + ε1 ε22 u7 (x, y, z, t ) + ε12 ε22 u8 (x, y, z, t ) + ε1−2 ε23 u9 (x, y, z, t ) + ε1−1 ε23 u10 (x, y, z, t ) + ε1−3 ε24 u11 (x, y, z, t ) + · · · Let

 ∂2 Gi jkl (y, z, t, τ ) dτ , ∂ τ ∂ xl 0    t  t ∂ ∂2 ∂ ∂2 =− Gi jkl (y, z, t, τ ) dτ − Gi jkl (y, z, t, τ ) dτ , ∂yj 0 ∂ τ ∂ xl ∂xj 0 ∂ τ ∂ yl  t  ∂ ∂2 =− Gi jkl (y, z, t, τ ) dτ , ∂yj 0 ∂ τ ∂ yl    t  t ∂ ∂2 ∂ ∂2 =− Gi jkl (y, z, t, τ ) dτ − Gi jkl (y, z, t, τ ) dτ , ∂xj 0 ∂ τ ∂ zl ∂zj 0 ∂ τ ∂ xl    t  t 2 ∂ ∂ ∂ ∂2 =− Gi jkl (y, z, t, τ ) dτ − Gi jkl (y, z, t, τ ) dτ , ∂yj 0 ∂ τ ∂ zl ∂zj 0 ∂ τ ∂ yl   t 2 ∂ ∂ =− G (y, z, t, τ ) dτ , ∂ z j 0 i jkl ∂ τ ∂ zl

∂ A1 = − ∂xj A2 A3 A4 A5 A6



(4)

t

(5)

Also, the chain rule can be given as:

∂ ∂ 1 ∂ 1 ∂ → + + , ∂xj ∂ x j ε1 ∂ y j ε2 ∂ z j Inserting (4) and (6) into (1), one obtains

−

 ∂ 2 u ε1 ε2 ( x , τ ) dτ ∂ τ ∂ xl 0   t   2  1 ∂ 1 ∂ ∂ ∂ u ε1 ε2 ( x , τ ) 1 ∂ 2 u ε1 ε2 ( x , τ ) 1 ∂ 2 u ε1 ε2 ( x , τ ) 1 ε2 =− + + Gεi jkl (x, t, τ ) + + dτ ∂ x j ε1 ∂ y j ε2 ∂ z j 0 ∂ τ ∂ xl ε1 ∂ τ ∂ yl ε2 ∂ τ ∂ zl  t 1 1 1 =− A1 u ε1 ε2 ( x , τ ) + A u ( x , τ ) + 2 A3 u ε1 ε2 ( x , τ ) + A4 u ε1 ε2 ( x , τ ) ε1 2 ε1 ε2 ε2 ε1 0

∂ ∂xj



t



1 ε2 Gεi jkl (x, t, τ )

(6)

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

1 1

+

ε1 ε2

A5 u ε1 ε2 ( x , τ ) +

and



227



1

A u ( x, τ ) d τ = f i ( x, t ) 2 6 ε1 ε2

(7)

ε2



t ∂ 2 u0 ∂ 3 u0 d τ , u2 ( x, y, t ) = N pq (y, t, τ ) dτ , ∂τ ∂ x ∂ τ ∂ x p ∂ xq p 0 0    t ∂ N p ( y ) ∂ 2 u0 ∂ 2 u0 u3 ( x, y, z, t ) = Mh (z, t, τ ) + Mh (z, t, τ ) dτ , ∂ yh ∂ τ ∂ x p ∂τ ∂ xh 0    t ∂ N pq (y, t, τ ) ∂ 3 u0 ∂ 3 u0 h h p u4 ( x, y, z, t ) = M (z, t, τ ) + M (z, t, τ )N (y, t, τ ) dτ , ∂ yh ∂ τ ∂ x p ∂ xq ∂τ ∂ xh ∂ x p 0   t ∂ 2 N pq (y ) ∂ 3 u0 ∂ 3 u0 u5 ( x, y, z, t ) = Mhr (y, z, t, τ ) + Mhr (y, z, t, τ ) ∂ y ∂ y ∂ τ ∂ x ∂ x ∂τ ∂ xh ∂ xr r p q 0 h

u1 ( x, y, t ) =

t

N p (y, t, τ )

∂ N p (y, t, τ ) ∂ 3 u0 ∂ N p (y, t, τ ) ∂ 3 u0 + Mhr (y, z, t, τ ) ∂ yh ∂ τ ∂ xr ∂ x p ∂ yr ∂ τ ∂ xh ∂ x p  3 pq ∂ u0 ∂ N (y, t, τ ) ∂ 3 u0 + Dh (y, z, t, τ )N r (y, t, τ ) + Dh (y, z, t, τ ) dτ , ∂τ ∂ xh ∂ xr ∂ yh ∂ τ ∂ x p ∂ xq + Mhr (y, z, τ )

(8)

where u0 = u0 (x, t) denotes the homogenization solutions based on macroscopic domain . Np (y, t, τ ), Npq (y, t, τ ), Mh (z, t, τ ), Mhr (y, z, t, τ ) and Dh (y, z, t, τ )(p, q, h, r=1, 2, 3) denote the unit cell solutions defined on the mesoscopic and microscopic domain, respectively. Also, they all can be derived in the following. Substituting the asymptotic expansion formulas (4) into (1), and considering the terms ε m ε n (m, n = −2, −1, 0, 1, 2) yield to the following equalities:



∂ − ∂xj =

t 0

1

ε22

A6 u0 + 1

+ ε1

ε2 1

+ ε12 +

1 Gi jkl (y, z, t, τ ) 2

ε2

1

ε12

1 1

ε1 ε2



A5 u0 +

  ∂ u ε1 ε2 k ( x , τ ) ∂ u ε1 ε2 l ( x , τ ) + dτ ∂ xl ∂ xk 1

ε12

A3 u0 + ε1

( A4 u1 + A5 u2 + A6 u4 ) +

1

ε1

1

ε22

A6 u1 + ε12

1

ε22

A6 u2 +

1

ε2

( A4 u0 + A5 u1 + A6 u3 )

( A2 u0 + A3 u1 + A5 u3 ) + ( A1 u0 + A2 u1 + A3 u2 + A4 u3 + A5 u4 + A6 u5 )

(A4 u2 + A6 u6 ) + ε1 (A1 u1 + A2 u2 + A4 u4 + A5 u6 + A6 u7 ) + ε12 (A1 u2 + A4 u6 + A6 u8 )

ε2 ( A3 u3 + A6 u9 ) +

1

ε13

ε22 (A5 u9 + A6 u11 ) +

1

ε1

ε2 (A2 u3 + A3 u4 + A5 u5 + A6 u10 ) + O(ε2 ) = fi (x, t )

(9)

Obviously, it follows from A6 u0 = 0, A5 u0 = 0, A3 u0 = 0, A6 u1 = 0 and A6 u2 = 0 that u0 (x, t) are independent of y and z, u1 (x, y, t) and u2 (x, y, t) are independent of z, respectively. Also, from (9), the following expressions can be obtained:

O

1 ε2 

O

O O

1

ε2 1 ε1 

ε12

1

ε2

O ( ε1 ) : O

 :

A4 u1 + A5 u2 + A6 u4 = 0, A2 u0 + A3 u1 + A5 u3 = 0,

:

 ε10 ε20 : A1 u0 + A2 u1 + A3 u2 + A4 u3 + A5 u4 + A6 u5 = fi (x, t ),

 O

ε1

A4 u0 + A5 u1 + A6 u3 = 0,

:

 2 ε1 :

 :

A4 u2 + A6 u6 = 0,

(10)

(11)

(12) (13)

(14)

A1 u1 + A2 u2 + A4 u4 + A5 u6 + A6 u7 = 0,

(15)

A1 u2 + A4 u6 + A6 u8 = 0,

(16)

228

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

 O O

1

ε12

1

ε1 

 ε2 :

A3 u3 + A6 u9 = 0,

 ε2 :

A2 u3 + A3 u4 + A5 u5 + A6 u10 = 0,

 2 ε : 3 2

1

O

(17) (18)

A5 u9 + A6 u11 = 0,

ε1

(19)

Then, the equalities (10)–(13) will be investigated in the following to obtain the microscopic and mesoscopic local functions, homogenized coefficients and homogenized equations, successively. First of all, by virtue of (8), Eq. (10) can be rewritten as:



t



0

    p h ∂ Mln ∂ Gi jnh (y, z, t, τ ) ∂ ∂ u0m ∂ ∂ u0n ∂ Nnm + G (y, z, t, τ ) + d τ = 0. ∂zj ∂ z j i jkl ∂ zk ∂τ ∂ xk ∂ yk ∂ x p

(20)

Since equality (20) holds for an arbitrary upper limit of time t, then the following partial differential equations can be obtained [33]:

⎧ ⎨

  h ∂ Mln (z, t, τ ) ∂ Gi jnh (y, z, t, τ ) ∂ Gi jkl (y, z, t, τ ) = , ∂zj ∂ zk ∂zj ⎩ h Mln (z, t, τ ) = 0, −

in Z,

(21)

z ∈ ∂ Z,

h (z, t, τ ) (h = 1, 2, 3 ) [44]. From Lax–Milgram theorem, (21) has one unique solution Mln

Next, taking (8) into (11), the following equation which represents the vanishing ε1 ε2−1 -term is formulated as follows:



t



0

   pq h ∂ Mln ∂ Gi jnh (y, z, t, τ ) ∂ ∂ 3 u0m ∂ Nnm ∂ 3 u0m p + Gi jkl (y, z, t, τ ) Nnm + dτ = 0 ∂zj ∂zj ∂ zk ∂ τ ∂ xk ∂ x p ∂ yk ∂ τ ∂ x p ∂ xq

(22)

Owing to (21), the equality (22) always holds. In analogy, and from (8) and (12), following partial differential equations are derived for the ε 1 −1 -term:

 p  2 ∂ Mlm ∂ ∂ u0m Gi jmp (y, z, t, τ )+Gi jkl (y, z, t, τ ) dτ ∂ zk ∂ τ ∂ x p 0 ∂yj   p 2  t h ∂ Mln ∂ ∂ Nnm ∂ u0m + Gi jnh (y, z, t, τ ) + Gi jkl (y, z, t, τ ) dτ ∂ zk ∂ yh ∂ τ ∂ x p 0 ∂yj   2  t p ∂ ∂ ∂ u 0m h ∂ Nnm + Gi jkl (y, z, t, τ ) Mln dτ = 0 ∂ yk ∂ yh ∂ τ ∂ x p 0 ∂zj 

t

(23)

Further, g(z) denotes an integrable function on microscopic domain Z ∈ Rn (n = 1, 2, 3 ) and let “−” represent the homogenized operator given by

g¯ =

1 |Z |



Z

g(z )dz

(24)

Applying the homogenized operator (24) to both sides of (23), and the homogenized equations can be obtained on the mesoscale Y in the following:

G¯ 1i jmp (y, t, τ )

∂yj

+

  ∂Np ∂ ¯1 Gi jkl (y, t, τ ) lm = 0. ∂yj ∂ yk

where

1 G¯ 1i jkl (y, t, τ ) = |Z |



 Z

(25)

∂ Mkpl Gi jkl (y, z, t, τ ) + Gi jpq (y, z, t, τ ) dz, ∂ zq

(26)

and G¯ 1i jkl (y, t, τ )(i, j, k, l = 1, 2, 3) is the mesoscopic homogenized parameters. From (25) and the periodic boundary conditions on Y, Nlm (y, t, τ ) ( p = 1, 2, 3 ) satisfies the partial differential equations given by: p

⎧   p 1 ⎨− ∂ G¯ 1 (y, t, τ ) ∂ Nlm (y, t, τ ) = G¯ i jmp (y, t, τ ) , ∂ y j i jkl ∂ yk ∂yj ⎩ p Nlm (y, t, τ )

in Y, is Y − periodic,

Similarly, one can obtain the following equality from (8) and (13):

(27)

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

 t













229





p pq ∂ Nlm ∂ Nlm ∂ 2 u0m ∂ ∂ 2 u0m ∂ ∂ 3 u0m ∂ ∂ 3 u0m ∂ p Gi jkl + Gi jkl + Gi jkl + Gi jkl Nlm ∂ x ∂ τ ∂ x ∂ x ∂ y ∂ τ ∂ x ∂ y ∂ y ∂ τ ∂ x ∂ x ∂ y ∂ τ ∂ x p ∂ xk p p q 0 j j j j k k k       p p h h ∂ Mln ∂ Mlm ∂ ∂ Nnm ∂ 2 u0m ∂ ∂ 3 u0m ∂ ∂ 2 u0m h ∂ Nnm + Gi jkl + Gi jkl Mln + Gi jkl ∂xj ∂ zk ∂ yh ∂ τ ∂ x p ∂zj ∂ yh ∂ τ ∂ xk ∂ x p ∂xj ∂ zk ∂ τ ∂ xh        3  pq pq h 3 3 ∂ Mln ∂ Nnm ∂ u0m ∂ ∂ u0m ∂ ∂ ∂ ∂ u0m h h ∂ Nnm + Gi jkl Mlm + Gi jkl + Gi jkl Mln ∂zj ∂ τ ∂ xk ∂ xh ∂yj ∂ zk ∂ yh ∂ τ ∂ x p ∂ xq ∂zj ∂ yk ∂ yh ∂ τ ∂ x p ∂ xq       pq h hr 3 3   ∂ Mln p ∂ u0m ∂ Mln ∂ u0m ∂ ∂ ∂ ∂ ∂ 2 Nnm ∂ 3 u0m p h + Gi jkl Nnm + Gi jkl Mln Nnm + Gi jkl ∂yj ∂ zk ∂ τ ∂ xh ∂ x p ∂zj ∂ yk ∂ τ ∂ xh ∂ x p ∂zj ∂ zk ∂ yh ∂ yr ∂ τ ∂ x p ∂ xq       p p hr 3 hr 3 hr ∂ Mlm ∂ u0m ∂ ∂ ∂ M ∂ Nnm ∂ u0m ∂ ∂ M ∂ Nnm ∂ 3 u0m + Gi jkl + Gi jkl + Gi jkl ∂zj ∂ zk ∂ τ ∂ xh ∂ xr ∂zj ∂ zk ∂ yh ∂ τ ∂ xr ∂ x p ∂zj ∂ zk ∂ yr ∂ τ ∂ xh ∂ x p     pq h h 3 3 ∂ Dln p ∂ u0m ∂ Dln ∂ Nnm ∂ u0m ∂ ∂ + G N + G d τ = f i ( x, t ). (28) ∂ z j i jkl ∂ zk nm ∂ τ ∂ xh ∂ x p ∂ z j i jkl ∂ zk ∂ yh ∂ τ ∂ x p ∂ xq

−

Also, imposing the homogenized operator (24) to both sides of equality (28) yields to the following equation:





 ∂ 3 u0m ∂ 3 u0m ∂  ¯1 ∂ N p ∂ 3 u0m p + + G¯ 1i jnh (y, t, τ ) nm Gi jkl (y, t, τ )Nlm ∂τ ∂ x j ∂ x p ∂ y j ∂ τ ∂ xk ∂ x p ∂ yh ∂τ ∂ x j ∂ x p 0   3  pq ∂ ¯1 ∂N ∂ u0m + d τ = − f i ( x, t ), G (y, t, τ ) nm ∂ y j jinh ∂ yh ∂ τ ∂ x p ∂ xq t

G¯ 1i jmp (y, t, τ )

(29)

Similarly to (24), one can define the integrable function d(y) on mesoscopic domain Y ∈ Rn (n = 1, 2, 3 ) given by

1 d¯ = |Y |



Y

d (y )dy.

(30)

Thus, utilizing the homogenized operator (30) to the equality (29), and the homogenized equations on macroscopic  are formulated given by:

    ⎧ ∂  t ¯2 1 ∂ ∂ u0k ( x, τ ) ∂ u0l ( x, τ ) ⎪ ⎪ − + d τ = f i ( x, t ), ⎪ 0 Gi jkl (t, τ ) ⎪ 2 ∂τ ∂ xl ∂ xk ⎨ ∂xj 0 v j σi j ( x , t ) = p i ( x , t ) , ⎪ ⎪ ⎪u0 (x, t ) = u¯ (x, t ), ⎪ ⎩ u0 ( x, 0 ) = 0,

where

G¯ 2i jkl

1 (t, τ ) = |Y |

 Y

G¯ 1i jkl

(y, t, τ ) +

G¯ 1i jpq

( x, t ) ∈  × ( 0, T ), ( x , t ) ∈ 1 × ( 0 , T ) , ( x , t ) ∈ 2 × ( 0 , T ) , x ∈ ,

∂ Nplk (y, t, τ ) dy ∂ yq

(31)

(32)

where G¯ 2i jkl (t, τ )(i, j, k, l = 1, 2, 3) is defined as the macroscopic homogenized coefficients. Furthermore, fi (x, t) (i = 1, 2, 3) can be taken by (31), and Eq. (29) can be rewritten in the following:





 ∂ 3 u0m ∂ 3 u0m ∂  ¯1 ∂ N p ∂ 3 u0m p + + G¯ 1i jnh (y, t, τ ) nm Gi jkl (y, t, τ )Nlm ∂τ ∂ x j ∂ x p ∂ y j ∂ τ ∂ xk ∂ x p ∂ yh ∂τ ∂ x j ∂ x p 0   3   t pq 3 ∂ ¯1 ∂N ∂ u0l ∂ u0m + dτ = dτ . G (y, t, τ ) nm G¯ 2i jkl (t, τ ) ∂ y j i jnh ∂ yh ∂ τ ∂ x p ∂ xq ∂τ ∂ x j ∂ xk 0 t

G¯ 1i jmp (y, t, τ )

(33)

Then, from (33), the mesoscale functions Nlm (y, t, τ ) ( p, q = 1, 2, 3 ) satisfies the following problems: pq

⎧   pq ∂ Nlm (y, t, τ ) ∂ ¯1 ⎪ ⎪ = G¯ 2iqmp (t, τ ) − G¯ 1iqmp (y, t, τ ) G (y, t, τ ) ⎪ ⎪ ∂ yk ⎨ ∂ y j i jkl  ∂Np ∂  ¯1 p ¯ 1 (y, t, τ ) lm , in Y, − ( y , t, τ ) N − G G ⎪ iqlk ⎪ lm ∂ yk ⎪ ∂ y j i jlq ⎪ ⎩ pq Nhm (y, t, τ ) is Y − periodic.

(34)

Motivated by the ideas in Refs. [6,20,21], the problems (27) and (34) have the unique solutions. Finally, substituting (29) into (28) and taking some simple derivation and combination, the equality is obtained as follows:

230

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242



    h ∂ G¯ 1i jnh (y, t, τ ) ∂ Gi jnh (y, z, t, τ ) ∂ Dhln ∂ Mln ∂ ∂ − G (y, z, t, τ ) + − − G (y, z, t, τ ) ∂ z j i jkl ∂ zk ∂yj ∂yj ∂ y j i jkl ∂ zk 0        t pq h hr 2 ∂ Mln ∂ Mlm ∂ ∂ 2 u0m ∂ ∂ Nnm p ∂ u0m 1 ¯ − G (y, z, t, τ ) + Nnm dτ + G (y, z, t, τ ) Gihmr (y, t, τ ) − ∂ z j i jkl ∂ yk ∂ yh ∂ x p ∂ xq ∂ xh ∂ x p ∂ z j i jkl ∂ zk 0  pq p r 2 2 2 2  ∂ u0m ∂ Mlm ∂  ∂ Nmn ∂ u0n ∂ Nmn ∂ u0n r − Gihmr (y, z, t, τ ) − Gihlk (y, z, t, τ ) − G (y, z, t, τ )Mlm + + ∂ zk ∂ z j i jlh ∂ xh ∂ xr ∂ yh ∂ yr ∂ x p ∂ xq ∂ yh ∂ xr ∂ x p  p 2 ∂ N ∂ u0n + mn dτ = 0 (35) ∂ yr ∂ xh ∂ x p 

t

To make sure (35) always hold, Dh (y, z, t, τ ) and Mhr (y, z, t, τ ) (h, r = 1, 2, 3) should satisfy the following equations, respectively.

⎧   ∂ G¯ 1i jnh (y, t, τ ) ⎪ ∂ Dhln (y, z, t, τ ) ∂ ⎪ ⎪ G ( y , z , t, τ ) = i jkl ⎪ ⎪ ∂zj ∂z ∂yj ⎪ ⎪  k  ⎪ h ⎪ ∂ Mln ∂ ⎨ ∂ Gi jnh (y, z, t, τ ) − − G (y, z, t, τ ) ∂yj ∂ y j i jkl ∂ zk ⎪   ⎪ h ⎪ ∂ ∂ Mln ⎪ ⎪ ⎪ ⎪− ∂ z j Gi jkl (y, z, t, τ ) ∂ yk , in Z, ⎪ ⎪ ⎩ h Dln (y, z, t, τ ) = 0, z ∈ ∂ Z,

(36)

⎧   hr ∂ Mlm (y, z, t, τ ) ∂ ⎪ ⎪ G (y, z, t, τ ) = G¯ 1ihmr (y, t, τ ) ⎪ ⎪ ∂ z j i jkl ∂ zk ⎪ ⎪ ⎪ r ⎨ ∂ Mlm −Gihmr (y, z, t, τ ) − Gihlk (y, z, t, τ ) ∂ zk ⎪   ⎪ ∂ ⎪ r ⎪ − G (y, z, t, τ )Mlm , in Z, ⎪ ⎪ ∂ z j i jlh ⎪ ⎩ hr Mlm (y, z, t, τ ) = 0, z ∈ ∂ Z.

(37)

Now, the three-scale asymptotic expansion solutions for the aging viscoelastic problems (1) can be determined by the following forms:

uε11 ε2 (x, t ) = u0 + ε1



t 0

N p (y, t, τ ) 

∂ 2 u0 dτ , ∂τ ∂ x p

∂ 3 u0 dτ , ∂τ ∂ x p ∂ xq 0    t ∂ N p (y, t, τ ) ∂ 2 u0 ∂ 2 u0 ε1 ε2 ε1 ε2 h h u3 ( x, t ) = u2 ( x, t ) + ε2 M (z, t, τ ) + M (z, t, τ ) dτ ∂ yh ∂ τ ∂ xp ∂τ ∂ xh 0    t ∂ 3 u0 ∂ N pq (y, t, τ ) ∂ 3 u0 + ε1 ε2 Mh (z, t, τ ) + N p (y, t, τ ) dτ ∂ y ∂ τ ∂ x ∂ x ∂τ ∂ xh ∂ x p p q 0 h  t  ∂ 2 N pq (y, t, τ ) ∂ 3 u0 ∂ 3 u0 + ε22 Mhr (y, z, t, τ ) + Mhr (y, z, t, τ ) ∂ y ∂ y ∂ τ ∂ x ∂ x ∂τ ∂ xh ∂ xr r p q 0 h uε21 ε2 (x, t ) = uε11 ε2 (x, t ) + ε12

t

N pq (y, t, τ )

∂ N p (y, t, τ ) ∂ 3 u0 ∂ N p (y, t, τ ) ∂ 3 u0 + Mhr (y, z, t, τ ) ∂ yh ∂ τ ∂ xr ∂ x p ∂ yr ∂ τ ∂ xh ∂ x p   3 pq ∂ u0 ∂ N (y, t, τ ) ∂ 3 u0 + Dh (y, z, t, τ )N r (y, t, τ ) + Dh (y, z, t, τ ) dτ , ∂τ ∂ xh ∂ xr ∂ yh ∂ τ ∂ x p ∂ xq + Mhr (y, z, t, τ )

ε ε

ε ε

(38)

ε ε

where u11 2 (x, t ), u21 2 (x, t ) and u31 2 (x, t ) are the first-order, the second-order two-scale and three-scale homogenization approximate solutions, respectively. Remark 3.1. The three-scale asymptotic solutions are equivalent to the solutions of original problem (1) in O(ε1−2 ε2 )-order pointwise sense, as shown in Appendix A. This is the main motivation that it is essential to establish the three-scale asymptotic expansions in this work. Also, the examples presented in Section 5 clearly demonstrate that it is important to contain the three-scale corrector terms to effectively predict the macroscopic properties of the composites.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Further,

ε1 ε2

σi j

ε1 ε2

( u3

from

(38),

the

three-scale

asymptotic

ε ε

(x, t )) (i, j = 1, 2, 3 ) can be formulated as follows:



1 2  t

eεkl1 ε2 (uε31 ε2 (x, t )) =

σiεj1 ε2 (uε31 ε2 (x, t )) =

0

ε ε

ekl1 2 (u31 2 (x, t )) (k, l = 1, 2, 3 )

strain

231

and

stress

fields

 ∂ u 3,ε1 ε2 k ( x , t ) ∂ u 3,ε1 ε2 l ( x , t ) + , ∂ xl ∂ xk

Gi jkl (y, z, t, τ )

∂ ε1 ε2 ε1 ε2 e (u3 (x, τ ))dτ . ∂τ kl

(39)

4. Numerical implementation for three-scale asymptotic expansion In this section, the numerical procedures are proposed to obtain the three-scale finite element formulas for the aging viscoelastic problems (1), and investigate the multiscale effective behavior of the problems. 4.1. Three-scale finite element formulas (1) Finite element (FE) formulas of unit cell functions and homogenized coefficients per scale hh FE solutions of microscopic functions Mln 0 (z, t, τ ) (h = 1, 2, 3 ) can be evaluated by computing the FE virtual work equation on microscale cell Z given by

 Z

Gi jkl (y, z, t, τ )

 hh0 ∂ Mln (z, t, τ ) ∂ vi ∂ vi dz = −Gi jnh (y, z, t, τ ) dz, ∀vi ∈ Sh0 (Z ). ∂ zk ∂zj ∂zj Z

(40)

where Sh0 (Z ) ⊆ H01 (Z ) represents the FE space of Z, Z is partitioned into FE set Sh0 , h0 is the mesh size, as illustrated in 1h Fig. 2(a). Thus, the finite element solutions G¯ 0 (y, t, τ ) (i, j, k, l = 1, 2, 3 ) of the mesoscale homogenized coefficients can be i jkl

obtained as follows: h0 (y, t, τ ) = G¯ 1i jkl

1 |Z |

 Z

(Gi jkl (y, z, t, τ ) + Gi jpq (y, z, t, τ )

∂ Mkplh0 )dz. ∂ zq

(41)

ph h

Also, the FE solutions of unit cell functions Nlm 0 1 (y, t, τ ) ( p = 1, 2, 3 ) can be obtained by solving the following FE virtual work equation on mesoscale cell Y

 Y

h0 (y, t, τ ) G¯ 1i jkl

 ph0 h1 ∂ Nlm (y, t, τ ) ∂ vi ∂v h0 dy = −G¯ 1i jkl (y, t, τ ) i dy, ∀vi ∈ Sh1 (Y ). ∂ yk ∂yj ∂yj Y

(42)

where Sh1 (Y ) ⊆ H01 (Y ) denotes the FE space of unit cell Y, Y is partitioned into FE set Sh1 , h1 is the mesh size of Y and 2h h is displayed in Fig. 2(b). Further, the approximate solutions G¯ 0 1 (t, τ ) (i, j, k, l = 1, 2, 3 ) of the macroscopic homogenized i jkl

coefficients are given by h0 h1 G¯ 2i jkl

1 (t, τ ) = |Y |



 Y

h0 G¯ 1i jkl

(y, t, τ ) +

h0 G¯ 1i jpq

∂ Nplkh0 h1 (y, t, τ ) dy. ∂ yq

(43)

Similarly to (40) and (41), the additional meso- and microscale Npq (y, t, τ ), Mhr (y, z, t, τ ) and Dh (y, z, t, τ )(p, q, h, r=1, 2, 3) can be derived by solving the (34), (36) and (37), successively. (2) FE formulas of homogenized solutions FE solutions of the homogenized problems (31) are obtained by following FE virtual work equations related with the homogenized coefficients (32) on 

 

t

 0

h0 h1 (t, τ ) G¯ 2i jkl

  ∂ ∂ uh02l (x, τ ) ∂ vi dτ dx = pi (x, t )vi dx + fi (x, t )vi dx, ∀vi ∈ Sh2 (), ∂τ ∂ xk ∂xj 1 

(44)

where Sh2 () indicates the FE space with mesh size h2 based on the domain , and is depicted in Fig. 2(c). As for the space domain, the FE methods are utilized to solve the homogenized Eq. (44), and the trapezoidal integration is used to handle the time integral for the same equations. Further, the homogenized problems (44) are evaluated on each subinterval separately, and they will be changed a recurrent formula including the solution on all previous subintervals on the right-hand side, i.e.,

 

t

 tn0 −1

h0 h1 (t, τ ) G¯ 2i jkl

  ∂ ∂ uh02l (x, τ ) ∂v dτ dx = pi (x, t )vdx + fi (x, t )vdx ∂τ ∂ xk ∂xj 1  −

n−1    n0 =1 

tn0

tn0 −1

h0 h1 (t, τ ) G¯ 2i jkl

∂ ∂ uh02l (x, τ ) ∂v dτ dx. ∂τ ∂ xk ∂xj

(45)

232

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 2. (a) Meshes of microscopic cell; (b) meshes of mesoscopic cell; (c) meshes of the homogenization domain .

Fig. 3. (a) Domain ; (b) mesoscopic cell Y; (c) microscopic cell Z.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 4. The first component of the displacement functions at time t= 30.0: (a)uh02 (x, t ); (b)uε21,hε2,h 0

1 ,h2

233

(x, t ); (c)uε31,hε02,h1 ,h2 (x, t ); (d) uFE .

4.2. FE procedures for three-scale asymptotic expansions In this subsection, the FE procedures of three-scale asymptotic expansions for predicting the aging viscoelastic properties of the problems (1) are described in the following: (1) The time domain is divided into several subintervals

tn = n t, 0 = t0 < t1 < t2 · · · < tn = T

(46)

(2) Verify the distribution of grains and the geometry of the structure in a microscale cell. Then, Z can be partitioned into FE meshes set. (3) For every ti (i = 0, 1, 2, …, n), evaluate the FE solutions of Mh (z, ti , τ ) (h = 1, 2, 3) according to the unit cell problems (40) based on the microscale domain Z. Thus, the mesoscale homogenized coefficients G¯ 1i jkl (y, ti , τ ) (i, j, k, l = 1, 2, 3 ) can be obtained by the formula (41). (4) Analyze the distribution of grains in mesoscopic domain Y and the geometric structure. Further, partition Y into finite element meshes set. (5) Evaluate the Np (y, ti , τ ) (p = 1, 2, 3) by the local cell problems (42) on the mesoscale cell Y. Meanwhile, the macroscopic homogenized coefficients G¯ 2i jkl (ti , τ ) (i, j, k, l = 1, 2, 3 ) are solved by the formula (43). (6) Calculate the FE solutions u0 (x, ti ) by homogenized aging viscoelastic problems (44) on the basis of G¯ 2i jkl (ti , τ ) (i, j, k, l = 1, 2, 3 ) computed from step 5. (7) Evaluate Npq (y, ti , τ )(p, q = 1, 2, 3) by solving the cell problems (34) using the same meshes as step 4. (8) Appling the same FE meshes proposed in step 2, Dh (y, z, ti , τ ) and Mhr (y, z, ti , τ )(h, r=1, 2, 3) are solved by the microscale cell problems (36) and (37) on Z, respectively.

234

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242 Table 1 FE information for the aging viscoelastic problems. Elements

Original equation 86,400

Mesoscale cell 782

microscale cell 824

Homogenized equation 20,0 0 0

Nodes

43,633

432

453

10,201

Table 2 Computational results for distinct methods in norm L2 .

ε 1 =1/6, ε 2 =1/36

error0 L2 /uF E L2

error2 L2 /uF E L2

error3 L2 /uF E L2

0.0665982

0.0434386

0.0374053

Table 3 Computational results for distinct methods in semi-norm H1 .

ε 1 =1/6, ε 2 =1/36

|error0 |H1 /|uF E |H1

|error2 |H1 /|uF E |H1

|error3 |H1 /|uF E |H1

0.73972

0.598839

0.288798

Fig. 5. Strain fields distribution at time t= 30.0: (a)eεkl1 ε2 (uh02 (x, t )); (b)eεkl1 ε2 (uε21,hε2,h 0

∂2u

1 ,h2

(x, t )); (c)eεkl1 ε2 (uε31,hε02,h1 ,h2 (x, t )); (d)eεkl1 ε2 (uF E ).

(9) Calculate the high-order partial derivatives ∂ x ∂ 0x (p, q = 1, 2, 3) by the average technique on relative elements (see p q Ref. [54]). (10) Evaluate the three-scale strain and stress fields by the formulas (39) according to assembling the unit cell solutions and homogenized solutions determined from step 2 to step 7.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 6. ε1 = 1/10, ε2 = 1/60, the first component of the displacement function at time t= 30.0: (a)uh02 (x, t ); (b)uε21,hε2,h 0

235

1 ,h2

(x, t ); (c)uε31,hε02,h1 ,h2 (x, t ); (d)uFE .

5. Numerical examples and discussions In this section, some typical numerical examples are proposed to verify the three-scale asymptotic homogenization presented in this paper. Example 1. Considering the ageing linear viscoelastic problems (1), where  is depicted in Fig. 3(a), and the mesoscopic domain Y and microscopic domain Z are illustrated in Fig. 3(b) and (c), respectively. The length and width of the rectangular inclusions in Y are 1/3, and the radii of the circular inclusions in Z are 0.3. Besides, ε1 = 1/6 and ε2 = 1/36. Set f(x, t) = (0, t−τ −10 0 0 0)T , pi (x, t) = 0 and u¯ (x, t ) = (0, 0 )T . In the example, the mechanical parameters of matrix are G(t, τ ) = 150(1 + e− 10 ) and ν = 0.2. The elastic moduli of inclusions are E2 =1, E3 =2 and Poisson’s ratio is ν = 0.3. Obviously, it is quite difficult and even impossible to get the theoretical solutions of the aging viscoelastic problems (1) in the periodic composites, we have to take uε1 ε2 (x, t ) by the FE solutions uFE in a fine mesh. Here, the linear triangular elements are employed for the viscoelastic problems (1) by fine meshes and the corresponding homogenized problems (31) by a coarse mesh. Table 1 gives the FE mesh information for the multiscale problems. h ε ε ε ε Also, it should be noted that u02 (x, t ) is the FE solutions of the homogenized Eq. (31), u21,h 2,h ,h (x, t ) and u31,h 2,h ,h (x, t ) 0

1

2

0

1

2

are the finite element solutions of second-order two-scale expansion and the three-scale asymptotic expansion. Further, the following notations are defined in this work for simplicity

error0 = uFE − uh02 (x, t ), error2 = uFE − uε21,hε2,h 0

1 ,h 2

(x, t ), error3 = uFE − uε31,hε02,h1 ,h2 (x, t ).

236

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 7. ε1 = 1/10, ε2 = 1/60, strain fields distribution at time t= 30.0: (a)eεkl1 ε2 (uh02 (x, t )); (b)eεkl1 ε2 (uε21,hε2,h 0

1 ,h2

(x, t )); (c)eεkl1 ε2 (uε31,hε02,h1 ,h2 (x, t )); (d)eεkl1 ε2 (uF E ).

And, the expressions are introduced as follows:

vL2 () =

 

2

|v| dx

1/2 , |v|H 1 () =

  

 1 / 2 . |∇v|2 dx

Therefore, the numerical errors of the three methods in L2 -norm and H1 -norm are depicted in Tables 2 and 3, respectively. h ε ε ε ε Fig. 4 displays the computational results for u02 (x, t ), u21,h 2,h ,h (x, t ), u31,h 2,h ,h (x, t ) and uFE in this example and at time 0

1

2

0

1

2

t= 30. ε ε Fig. 5 displays the computational results for the strain fields ekl1 2 (uε1 ε2 (x, t )) at time t= 30.0. h

ε ε

Fig. 6 illustrates the computational results for the u02 (x, t ), u21,h 2,h 0

1 ,h2

(x, t ), uε31,hε2,h 0

1 ,h2

(x, t ) and uFE in this example at

time t= 30.0, and ε1 = 1/10, ε2 = 1/60. ε ε Fig. 7 illustrates the computational results for the strain fields ekl1 2 (uε1 ε2 (x, t )) at time t= 30.0, and ε1 = 1/10, ε2 = 1/60. By analyzed the numerical results in Figs. 4–7 and Table 2–3, it can be found that the homogenized methods and secondorder two-scale methods have not good agreement with the FE solution, but the three-scale asymptotic expansion performs excellently and clearly demonstrates better approximations of the FE results than the results obtained by traditional homogenized methods and second-order two-scale methods. In addition, it should be pointed out that the second-order two-scale methods can effectively catch the mesoscale behaviors for the aging viscoelastic problems, however, the microscale information of the composite materials cannot be captured accurately.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

237

Fig. 8. (a) Domain ; (b) mesoscopic unit cell Y; (c) microscopic unit cell Z.

Table 4 FE information for the aging viscoelastic problems. Elements

Mesoscale cell 84,389

microscale cell 82,457

Nodes

14,368

14,029

Example 2. The 3D composites with multiple spatial scales are shown in Fig. 8(a), and the mesoscopic and microscopic domains are illustrated in Fig. 8(b) and (c), respectively. The sizes of the long axes of the ellipsoidal inclusions are chosen as 1/3, middle axes 1/5 and short axes 1/4. The length, width and height of the rectangular inclusions in Y are 1/3, and the radii of the spherical inclusions in Z are 0.3. The mechanical parameters are the same with that proposed in Example 1. Also, the linear tetrahedral elements are adopted for the mesoscale problems and microscale problems (26) and (32). Besides, the information of the 3-D FE meshes is given in Table 4.

In order to study the aging viscoelastic performance of the heterogeneous materials with spatial scales, the effective relaxation moduli are evaluated and the representative results in mesoscale and macroscale are shown in Fig. 9(a) and (b), respectively. In addition, the macroscopic effective relaxation moduli obtained for the ellipsoidal inclusions at the mesoscale are illustrated in Fig. 10. As a result of the numerical analysis in Figs. 8–10 and Table 4, the three-scale asymptotic approximate solutions are adequately valid to predict the macroscopic effective behavior of the aging viscoelastic problems in the heterogeneous materials. In particular, the distinct inclusions in the mesoscopic domain can also be accurately obtained by the three-scale asymptotic expansion developed in this work.

238

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 9. Effective relaxation moduli for (a) mesoscopic domain and (b) macroscopic domain for the composites.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

Fig. 10. Macroscopic effective relaxation moduli for the ellipsoidal inclusion at the mesoscale (a) τ = 0 (b) τ = 15.

239

240

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

6. Conclusions In this paper, the efficient three-scale asymptotic expansion is given for predicting the aging viscoelastic performance of composite materials with multiple small-scale configuration. By the reiterated homogenization method, the three-scale asymptotic expansions formulas are derived in detail, including the three-scale asymptotic solutions of displacement fields, the homogenization relaxation modulus and the three-scale FE algorithms for strain and stress fields in time domain. Finally, the reasonableness of the three-scale model and the availability of the developed three-scale asymptotic expansions are verified by comparing with the FE methods with fine meshes. According to the FE verifications and analysis, it is well illustrated that the three-scale model is available to simulate the problems of the composites by introducing the microscale and mesoscale higher-order corrections and the two kinds of homogenized parameters. Numerical examples clearly show that the computational results evaluated by the three-scale asymptotic methods agree fairly well with those obtained by the direct FE methods. Through adding higher-order correction terms and homogenized solutions, the detailed oscillation information into the microscopic and mesoscopic domain of composites can be obtained effectively, especially, the local strain and stress distributions. In addition, the three-scale asymptotic expansion and related FE techniques presented in this study can be extended to handle the other nonlinear FE analysis such as elastoplasticity and damage analysis with complex microstructures. Also, this approach can be developed to analyze the heterogeneous structure with a random configuration. Acknowledgments This work is supported by the Fundamental Research Funds for the Central Universities and also supported by the National Natural Science Foundation of China (11701123). Appendix A. Theoretical analysis of the three-scale asymptotic solutions ε ε

ε ε

ε ε

To compare u11 2 (x, t )and u21 2 (x, t ) with the original solutions, we take uε1 ε2 (x, t ) − u11 2 (x, t ) and uε1 ε2 (x, t ) − u2 (x, t ) into the aging viscoelastic problems (1), and have the following equalities: ε1 ε2





Lε1 ε2 uε1 ε2 (x, t ) − uε11 ε2 (x, t ) = fi (x, t ) + +

1

ε1

1

ε2

( A4 u0 + A5 u1 ) + ε1

1

ε2

A4 u1

( A2 u0 + A3 u1 + ( A1 u0 + A2 u1 ) + ε1 A1 u1 ,

(A.1)

and





Lε1 ε2 uε1 ε2 (x, t ) − uε21 ε2 (x, t ) = fi (x, t ) +

1

ε2

( A4 u0 + A5 u1 ) + ε1

+ (A1 u0 + A2 u1 + A3 u2 ) + ε12 where Lε1 ε2 = − ∂∂x ( j

t

1

ε2

1

ε2

( A4 u1 + A5 u2 ) +

1

ε1

( A2 u0 + A3 u1 )

(A4 u2 + ε1 (A1 u1 + A2 u2 ) + ε12 A1 u2 ,

(A.2)

Gi jkl (y, z, t, τ ) 12 ( ∂ τ∂∂ x + ∂ τ∂∂ x ) ). It is not difficult to see that the residual (A.1) and (A.2) are the order l k −1 O(ε2 ). But in the actual engineering computations, ε 1 and ε 2 are the fixed smaller constant rather than tending to zero. The error O(ε2−1 ) is not accepted for the engineers who want to capture the microscopic behavior of the composites. So, the first-order and second-order two-scale solutions are not adopted for some practical engineering computations. ε ε Substituting uε1 ε2 (x, t ) − u31 2 (x, t ) into (1), we can obtain that



2

2

0



Lε1 ε2 uε1 ε2 (x, t ) − uε31 ε2 (x, t ) =

ε12 ε2 ε2 A u + ε (A u + A u + A u ) + ε12 (A1 u2 ) + 2 A3 u3 + (A2 u3 + A3 u4 + A5 u5 ) ε2 4 2 1 1 1 2 2 4 4 ε1 ε1 (A.3)

Note that the residual of (A.3) is the order O(ε1−2 ε2 ). The three-scale asymptotic expansion solutions are equivalent to the solutions of the problems (1) in O(ε1−2 ε2 ) -order pointwise sense, also, referring to numerical examples proposed in Section 5. It is the reason that we consider the three-scale asymptotic expansions in this paper. References [1] [2] [3] [4] [5] [6] [7]

P.S. Mangat, M.M. Azari, A theory of the creep of steel fiber reinforced cement matrices under compression, J. Mater. Sci. 20 (1985) 1119–1133. S. Kuznetsov, J. Fish, Mathematical homogenization theory for electroactive continuum, Int. J. Numer. Methods Eng. 91 (2012) 1199–1226. R.M. Christensen, A critical evaluation for a class of micro-mechanics models, J. Mech. Phys. Solids 38 (3) (1990) 379–404. A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, American Mathematical Society, Rhode Island, 2011. D. Cioranescu, P. Donato, An Introduction to Homogenization, Oxford University Press, 1999. O.A. Oleinik, A.S. Shamaev, G.A. Yosifian, Mathematical problems in elasticity and homogenization. North-Holland, Amsterdam, 1992. V.V. Jikov, S.M. Kozlov, O.A. Oleinik, Homogenization of Differential Operators and Integral Functions, Springer, Berlin, 1994.

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

241

[8] I. Temizer, P. Wriggers, Homogenization in finite thermoelasticity, J. Mech. Phys. Solids 59 (2011) 344–372. [9] Z.H. Li, Q. Ma, J.Z. Cui, Second-order two-scale finite element algorithm for dynamic thermo-mechanical coupling problem in symmetric structure, J. Comput. Phys. 314 (2016) 712–748. [10] R.U. PatilB, K. MishraI, V. Singh, A new multiscale XFEM for the elastic properties evaluation of heterogeneous materials, Int. J. Mech. Sci. 122 (2017) 277–287. [11] J. Fish, Q. Yu, K. Shek, Computational damage mechanics for composite materials based on mathematical homogenization, Int. J. Numer. Methods Eng. 45 (1999) 1657–1679. [12] W.N. E, B. Engquist, The heterogenous multiscale methods, Commun. Math. Sci. 1 (2003) 87–132. [13] A. Abdulle, A. Nonnenmacher, Adaptive finite element heterogeneous multiscale method for homogenization problems, Comput. Method Appl. Math. 200 (2011) 2710–2726. [14] G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23 (6) (2003) 1482–1518. [15] T.Y. Hou, X.H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1997) 169–189. [16] H.W. Zhang, J.K. Wu, Z.D. Fu, Extended multiscale finite element method for elasto-plastic analysis of 2D periodic lattice truss materials, Comput. Mech. 45 (2010) 623–635. [17] T.J.R. Hughes, Multiscale phenomena: green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Method Appl. Math. 127 (1995) 387–401. [18] N. Zabaras, B. Ganapathysubramanian, A stochastic multiscale framework for modeling flow through random heterogeneous porous media, J. Comput. Phys. 228 (2009) 591–618. [19] X.G. Yu, J.Z. Cui, The prediction on mechanical properties of 4-step braided composites via two-scale method, Compos. Sci. Technol. 67 (2007) 471–480. [20] Z.Q. Yang, J.Z. Cui, Y.F. Nie, Q. Ma, The second-order two-scale method for heat transfer performances of periodic porous materials with interior surface radiation, CMES: Comput. Model Eng. 88 (5) (2012) 419–442. [21] Z.Q. Yang, J.Z. Cui, Y. Sun, J.R. Ge, Multiscale computation for transient heat conduction problem with radiation boundary condition in porous materials, Finite Elem. Anal. Des. 102–103 (2015) 7–18. [22] X.F. Guan, X. Liu, X. Jia, Y. Yuan, J.Z. Cui, H.A. Mang, A stochastic multiscale model for predicting mechanical properties of fiber reinforced concrete, Int. J. Solids Struct. 56–57 (2015) 280–289. [23] G. Allaire, Z. Habibi, Second order corrector in the homogenization of a conductive-radiative heat transfer problem, Discret. Contin. Dyn. – B 18 (1) (2013) 1–36. [24] Y.M. Yi, S.H. Park, S.K. Youn, Asymptotic homogenization of viscoelastic composites with periodic microstructures, Int. J. Solids Struct. 35 (17) (1998) 2039–2055. [25] S.T. Liu, K.Z. Chen, X.A. Feng, Prediction of viscoelastic property of layered materials, Int. J. Solids Struct. 41 (13) (2004) 3675–3688. [26] Y.M. Cai, H.Y. Sun, Prediction on viscoelastic properties of three-dimensionally braided composites by multi-scale model, J. Mater. Sci. 48 (19) (2013) 6499–6508. [27] A.B. Tran, J. Yvonnet, Q.C. He, C. Toulemonde, J. Sanahuja, A simple computational homogenization method for structures made of linear heterogeneous viscoelastic materials, Comput. Methods Appl. Mech. Eng. 200 (2011) 2956–2970. [28] Y. Zhang, Y.F. Nie, Y.T. Wu, The statistical two-scale method for predicting viscoelastic properties of composites with consistent random distribution of particles, Appl. Mech. Mater. 697 (2015) 3–6. [29] S.T. Nguyena, M.Q. Thai, M.N. Vu, Q.D. To, A homogenization approach for effective viscoelastic properties of porous media, Mech. Mater. 100 (2016) 175–185. [30] J. Sanahuja, Effective behaviour of ageing linear viscoelastic composites: homogenization approach, Int. J. Solids Struct. 50 (2013) 2846–2856. [31] S. Maghous, G.J. Creus, Periodic homogenization in thermoviscoelasticity: case of multilayered media with ageing, Int. J. Solids Struct. 40 (2003) 851–870. [32] F. Lavergnea, K. Saba, J. Sanahuja, M. Bornert, C. Toulemonde, Homogenization schemes for aging linear viscoelastic matrix-inclusion composite materials with elongated inclusions, Int. J. Solids Struct. 80 (2016) 545–560. [33] Y. Zhang, J.Z. Cui, Y.F. Nie, Second-order two-scale computational method for ageing linear viscoelastic problem in composite materials with periodic structure, Appl. Math. Mech. – Engl. Ed. 37 (2) (2016) 253–264. [34] G. Allaire, M. Briane, Multiscale convergence and reiterated homogenization, Proc. R. Soc. Edinb A 126 (1996) 297–342. [35] A. Holmbom, N. Svanstedt, N. Wellander, Multiscale convergence and reiterated homogenization of parabolic problems, Appl. Math. – CZECH 50 (2005) 131–151. [36] D. Trucu, M.A.J. Chaplain, A. Marciniak-Czochra, Three-scale convergence for processes in heterogeneous media, Appl. Anal. 91 (2012) 1351–1373. [37] A. Abdulle, Y. Bai, Fully discrete analysis of the heterogeneous multiscale method for elliptic problems with multiple scales, IMA J. Numer. Anal. 35 (1) (2015) 133–160. [38] A. Almqvist, E.K. Essel, J. Fabricius, P. Wall, Reiterated homogenization applied in hydrodynamic lubrication, Proc. Inst. Mech. Eng. J – J. Eng. 222 (2008) 827–841. [39] X.F. Guan, X. Liu, X. Jia, Y. Yuan, J.Z. Cui, H.A. Mang, A stochastic multiscale model for predicting mechanical properties of fiber reinforced concrete, Int. J. Solids Struct. 56–57 (2015) 280–289. [40] J.L. Zhang, X. Liu, Y. Yuan, H.A. Mang, Multiscale modeling of the effect of the interfacial transition zone on the modulus of elasticity of fiber-reinforced fine concrete, Comput. Mech. 55 (2015) 37–55. [41] Q. Chen, H.H. Zhu, Z.G. Yan, J. Woody Ju, Z.W. Jiang, Y.Q. Wang, A multiphase micromechanical model for hybrid fiber reinforced concrete considering the aggregate and ITZ effects, Constr. Build. Mater. 114 (2016) 839–850. [42] E.I. Rodríguez, M.E. Cruz, J. Bravo-Castillero, Reiterated homogenization applied to heat conduction in heterogeneous media with multiple spatial scales and perfect thermal contact between the phases, J. Braz. Soc. Mech. Sci. 38 (2016) 1333–1343. [43] E.S. Nascimento, M.E. Cruz, J. Bravo-Castillero, Calculation of the effective thermal conductivity of multiscale ordered arrays based on reiterated homogenization theory and analytical formulae, Int. J. Eng. Sci. 119 (2017) 205–216. [44] L.Q. Cao, Iterated two-scale asymptotic method and numerical algorithm for the elastic structures of composite materials, Comput. Method Appl. Math. 194 (2005) 2899–2926. [45] R. Mahnken, C. Dammann, A three-scale framework for fibre-reinforced-polymer curing Part I: microscopic modeling and mesoscopic effective properties, Int. J. Solids Struct. 100–101 (2016) 341–355. [46] R. Mahnken, C. Dammann, A three-scale framework for fìbre-reinforced-polymer curing part II: mesoscopic modeling and macroscopic effective properties, Int. J. Solids Struct. 100–101 (2016) 356–375. [47] A. Ramírez-Torres, R. Penta, R. Rodríguez-Ramos, J. Merodio, F.J. Sabina, J. Bravo-Castilleroc, R. Guinovart-Díazc, L. Preziosi, A. Grilloa, Three scales asymptotic homogenization and its application to layered hierarchical hard tissues, Int. J. Solids Struct. 130–131 (2018) 190–198. [48] Z.H. Yang, Y. Zhang, H. Dong, J.Z. Cui, X.F. Guan, Z.Q. Yang, High-order three-scale method for mechanical behavior analysis of composite structures with multiple periodic configurations, Compos. Sci. Technol. 152 (2017) 198–210. [49] Y. Zhang, J.Z. Cui, Y.F. Nie, H. Dong, Z.H. Yang, High-order triple-scale method for composite structures of the configurations with small periodicities of two-levels. 2016, WCCM XII & APCOM VI, Seoul. [50] Z.Q. Yang, Y. Sun, J.Z. Cui, Z.H. Yang, T.Y. Guan, A three-scale homogenization algorithm for coupled conduction-radiation problems in porous materials with multiple configurations, Int. J. Heat Mass Transf. 125 (2018) 1196–1211.

242

Z. Yang, Y. Sun and Y. Liu et al. / Applied Mathematical Modelling 71 (2019) 223–242

[51] H.V. Nguyen, J. Pastor, D. Muller, A method for predicting linear viscoelastic mechanical behavior of composites a comparison with other methods and experimental validation, Eur. J. Mech. A – Solid 14 (1995) 939–960. [52] J. Orlik, Homogenization for viscoelasticity of the integral type with aging and shrinkage, Tech. Rep. (1998). [53] J. Orlik, Existence and stability estimate for the solution of the ageing hereditary linear viscoelasticity problem, Abstr. Appl. Anal. 2009 (2009) 1–19. [54] J.Z. Cui, The two-scale expression of the solution for the structure with several sub-domains of small periodic configurations. Invited Presentation “Workshop on Scientific Computing 99” June 27-30, 1996, Hong Kong.