1 Introduction

With the development of the material science and technology, composite materials have been widely used in engineering owing to their excellent mechanical properties. Recently, the homogenization method to predict the effective mechanical properties of composites and the relative multiscale analysis method for the elastic problem have attracted a lot of attention from scientists and engineers. However, many composite materials, such as polymer and concrete, exhibit time and rate dependence which can be modelled in the framework of viscoelasticity at a normal temperature. As the temperature rises, the viscoelasticity of composite materials is more significant.

There have been some reports about estimating the non ageing linear viscoelastic mechanical prop- erties of composite materials^{[1, 2, 3]}. In these references, the effective properties of non ageing viscoelastic composite materials are predicted by the homogenization method routinely through the correspondence principle. Actually, the Laplace-Carson transform changes the non ageing viscoelastic problem into the elastic case, and the homogenization schemes of the elastic composites can be used directly in the Laplace domain. Then, through the inversion of the Laplace-Carson transform, the effective behavior in the time domain is derived. However, when the phases of composites exhibit ageing, the correspon- dence principle is no longer valid. This is the major difficulty in extending the existing homogenization and multiscale method of the elastic problem to the composites with ageing. To avoid this shortcoming of the Laplace-Carson transform, some new techniques in the time domain have been designed to deal with the non-aging viscoelastic problem^{[4, 5, 6, 7]}. Some of these techniques have been extended to ageing linear viscoelastic problem successfully^{[8, 9]}. These methods are a homogenization approach to the ef- fective behavior of the composites actually, and most of them only contain the information of volume fractions for composites.

The asymptotic homogenization method^[10] is a mathematical theory originating from the study of partial differential equations with the rapidly varying coefficients and is an alternative approach to evaluate the effective properties of composite materials. Compared with the other micromechanics methods, such as Mori-Tanaka method^[11] and self-consistent method^[12], the asymptotic homogeniza- tion method has a rigorous mathematical background, and it can contain almost all the microcosmic structure information of the composites. It should be pointed out that the above methods cannot effectively capture the microcosmic oscillation information in the composite materials. In other words, the local stress and strain field cannot be more precisely calculated by most of these methods. To solve this problem, the researchers^{[13, 14, 15, 16]} have established a second-order two-scale (SOTS) method based on the asymptotic homogenization method, to predict the thermal and mechanical properties of the composite materials. The SOTS method has been proved to be able to effectively capture the local strain and stress information in the microscopic structure through the two-order displacement solution in the elastic case. This paper is devoted to designing an SOTS computational method to predict the ageing linear viscoelastic performance of composite materials with the small periodic structure operat- ing in the time domain directly. To deal with the ageing linear viscoelastic problem, a new high-order asymptotic expansion of the displacement solution is constructed at first. Further, the expression of homogenization coefficient, the homogenization equation in the macroscopic domain, and the governing equations in the microscopic scale are derived analytically. Through solving the problems mentioned above, the displacement approximate solution can be obtained.

The remainder of this paper is organized as follows. In Section 2, the SOTS formulations for calculating the effective relaxation modulus and displacement approximate solutions of the ageing viscoelastic problem in the periodic structure are derived formally. Then, the error analysis of the displacement approximate solutions is presented, and the algorithm flow of SOTS method is listed in Section 3. Finally, some numerical results are demonstrated to show that the SOTS method is valid and effective.

For convenience, we use the Einstein summation convention on repeated indices in this paper.

2 SOTS formulation

In this section, the SOTS formulation is derived for calculating the ageing viscoelastic properties and displacement fields for composite materials with the periodic structure. The ageing viscoelastic initial boundary value problem can be expressed as follows:

where G_ijhk^ε(x, t, τ ) denotes the relaxation tensor, providing the stress response along i, j at the time t resulting from a unit strain step along h, k occurring at the time τ . For τ > t, it can be extended as G_ijhk^ε(x, t, τ ) = 0. u^ε(x, t) denotes the displacement vector, f_i(x, t) denotes the body force, and Γ1 and Γ2 denote the stress and displacement boundaries, respectively. The structure is periodic and can be denoted as , where Z is the set of the three-dimensional integer vectors, and Y is the unit cell, as shown in Fig. 1.

At first, we make the following suppositions:

(A) G_ijhk^ε(x, t, τ ) is symmetric, G_ijhk^ε(x, t, τ ) = G_ijkh^ε(x, t, τ ) = G_jihk^ε(x, t, τ ),

(B) G_ijhk^ε(x, t, τ ) is bounded, and there exist two positive constants λ₁ and λ₂ such that

where ηij is an arbitrary symmetric matrix,

(C) G_ijhk^ε(x, t, τ ) and u^ε(x, t) are smooth functions in the domain Ω × (0, T),

(D)

Under the above suppositions, set y = x/ε ∈ Y , and y is the local coordinate on the unit cell, which contains the microcosmic information. Then, the relaxation modulus G_ijhk^ε(x, t, τ ) can be expressed as G_ijhk^ε(x, t, τ ) = G_ijhk^ε(y, t, τ ). As the elastic case of the asymptotic homogenization method, the displacement solution u^ε(x, t) is supposed to be expressed as the following form:

where u^ε(x, t) is the homogenization solution in the whole domain and denotes the macroscopic behavior of composites.

Due to y = x/ε , according to the chain rule , substituting Eq. (3) into the first equation in Eq. (1) and matching the terms of the same order of " , we can obtain a series of equations. From the coefficients of ε⁻¹, the following equation is obtained:

Then, from the coefficients of ε⁰, we obtain the following equation:

In the elastic case, u₁ is defined as u₁(x, y) = N_α1 (y) . For the aging viscoelastic problem in this paper, we define u₁ as the following expression constructively:

where N_α1 is a continuous function and defined in the unit cell Y . Introducing Eq. (6) into Eq. (4), we exchange the order of differential and integral operators by using the continuity of functions. Then, applying the subscript transform, the following equation is obtained:

According to the above suppositions, the integrand is a continuous function, and Eq. (7) holds for the arbitrary upper limit of the integral t. It should be pointed out that u₀ (x, t) is the homogenization solution, and is not identically zero. Then, the following equation is obtained:

Attach the following periodic boundary condition, and N_α1 (y, t, τ ) is the solution of the following partial differential equation: This is an elliptic partial differential equation with the parameter t. According to the Lax-Milgram lemma, Korn’s inequalities, and the assumptions (A)-(D), it is easy to prove that the problem (9) has a unique solution for arbitrary fixed t, , α₁, and m^[17]. Similarly, u₂ is defined as the following expression: Introducing Eqs. (6) and (10) into Eq. (5), we obtain Then, making integral on both sides of Eq. (11) with the local variable y in the unit cell Y , according to the continuity of functions and symmetry of the unit cell, we obtain where $\tilde G$ is called the homogenization relaxation modulus and has the following expression: Equation (12) is named as the homogenization equation, and the homogenization displacement solution u₀(x, t) is the solution of this equation.

Introducing Eq. (12) into Eq. (11), similar to the previous derivation for Eq. (8), we obtain

Attach the following periodic boundary condition, and N_α1α2 (y, t,τ) is the solution of the following partial differential equation: It also can be proved that the problem (15) has a unique solution for the arbitrary fixed t, , α1, α2, and m.

Now we can define the two-scale approximate solutions for the ageing linear viscoelastic problem in the periodic structure as

where u₁^ε is named as the first-order two-scale approximate solution, and u₂^ε is named as the SOTS approximate solution.

3 Error analysis and algorithm procedure 3.1 Error analysis of two-scale approximate solutions

In this subsection, we give a brief error analysis of the two-scale approximate solutions to show the necessity of calculating the second-order solution. To compare the first-order two-scale approximate solution u₁^ε with the solution of the original problem u^ε, introducing u^ε(x, t) − u₁^ε(x, t) into the original ageing viscoelastic problem, we obtain

Notice that the order of the residual on the right side of the above equation is O(1), and in the engineering, the parameter ε is a constant. The first-order approximate solution cannot effectively catch the local oscillatory behavior, and therefore it is necessary to develop the high order approximate solution.

Then, to compare the SOTS solution u₂^ε with the solution of the original problem u^ε, introducing u^ε(x, t) − u₂^ε(x, t) into the original ageing viscoelastic problem, we obtain

where The order of the residual on the right side of Eq. (19) is O(ε), so the second-order approximate solution has a better approximation to the original solution compared with the first-order solution. In addition, we give the main convergence theorem for the second-order approximate solution.

Theorem 1 Suppose that ⊂ Rⁿ is the union of entire periodic cells. Let u^ε(x, t) be the solution of Eq. (1), u₀(x, t) be the solution of homogenized Eq. (12), and u₂^ε (x, t) be the approximate solution which is defined by Eq. (17). Under the assumptions (A)-(D), we obtain the following error estimate:

where C is a constant independent of ε, but dependent of T.

Remark 1 In a similar way of proving the elastic problem^[18], we can easily obtain the error estimate at the arbitrary time point. Then, taking the infinite norm in the interval (0, T), Theorem 1 can be obtained easily.

3.2 Algorithm procedure for SOTS method

In this subsection, we describe the algorithm procedure of the SOTS method for the ageing linear viscoelastic problem. In the space domain, we use the finite element method to solve the equations, and the trapezoidal numerical integration is used to deal with the time integral in the equations. The linear tetrahedral element is used in the finite element method. The finite element algorithms for solving problems (9) and (15) have been introduced in detail in Ref. ^[19]. For the homogenization problem (12), we can solve the problem on each subinterval separately and get a recurrent formula containing the solution on all previous subintervals on the right-hand side, i.e.,

The whole procedure of SOTS method is as follows:

(1) Divide the time domain into several subintervals

and form the geometry of the structure, then, partition the unit cell Y into the finite element mesh.

(2) For every ti (i = 0, 1, 2, · · · , n), compute N_α1 according to the problem (9) by the finite element method in the unit cell Y . Then, the homogenization coefficient eG is computed by the formula (13).

(3) Through the recursion formula in the time domain, the homogenization displacement solution u₀ is obtained by solving the problem (12) in the whole structure .

(4) Applying the same mesh in the step (1), compute N_α1α2 by solving the problem (15) in the unit cell.

(5) From Eqs. (16) and (17), the two-scale approximate solutions are evaluated, respectively.

4 Numerical examples

In order to demonstrate the validity and feasibility of the SOTS method for the ageing viscoelastic problem, some numerical examples are given in this section. The macrostructure and the unit cell Y = [0, 1]³ are shown in Fig. 1. We consider that ε = 1/5.

4.1 Example 1

We consider that the matrix is an ageing viscoelastic material and has the relaxation modulus and Poisson’s ratio γ = 0.2. The particles have elastic modulus E = 120, Poisson’s ratio γ = 0.3. Poisson’s ratio is considered as a constant in the whole computational process. We set the body force f_i = (0, 0, −1 000), and the surface force p_i = (0, 0, 200) is applied on the top surface of in the z-direction. At first, according to the formula (13), the homogenization relaxation modulus is evaluated in the unit cell. The numerical result for the load time τ = 0 is shown in Fig. 2. The homogenization result is compared with the Voigt-Reuss (V-R) upper and lower bounds, and it shows that the method in this paper is effective to evaluate the homogenization relaxation modulus of composite materials.

Then, the homogenization displacement solution u₀(x, t), first-order approximate solution u₁^ε (x, t), and second-order approximate solution u₂^ε (x, t) can be calculated, respectively, based on the previous formulations. Since it is difficult to find the analytical solution of Eq. (1), we have to take u^ε(x, t) as their finite element solution ue in the very fine mesh, and implement the tetrahedron partition for . It should be pointed out that u₀(x, t) is evaluated by solving the homogenization equation in a coarse mesh. The numbers of elements and nodes are listed in Table 1. Set the error e₀ = u_e−u₀, e₁ = u_e−u₁^ε , and e₂ = ue−u₂^ε . For convenience, we introduce the norm . The relative numerical errors of the homogenization, first-order, and SOTS displacement solutions in the norm H1 at the time t = 30 are listed in Table 2. Figure 3 shows the numerical results of displacement solutions u₀₃, u₁₃^ε, u₂₃^ε, and u_e3 at the intersection z = 0.5, the calculating time and the load time are t = 30 andτ= 0, respectively. Figure 4 illustrates the numerical results of the displacement solution u₂^ε at the different calculating time at the intersection z = 0.5. Table 2 and Fig. 3 show that the SOTS approximate solution is in good agreement with the finite element solution in a very fine mesh, and the second-order solution u₂₃^ε is much better than the homogenization solution u₀, though it is just a little better than u₁₃^ε in this example. The reason is that the exact solution for this example has little oscillating behaviour actually. Another example is shown in the next subsection to show the necessity of second-order approximate solution. Figure 4 shows that the displacement is increasing as the time increases, due to the aging viscoelastic properties of the matrix. But the rate of change of the displacement is gradually reducing, and this is actually determined by the quality of the relaxation modulus in this example. Further, we choose the different values of ε and compute the relative errors of the second-order approximate solutions. Figure 5 shows the results for ε = 1/2, 1/3, 1/5, and the relative errors are decreasing as ε decreases.

4.2 Example 2

In this example, the mechanical parameters of matrix are. The particles have elastic modulus E = 5, Poisson’s ratio γ = 0.3. We set f_i = (0, 0, −10 000) and p_i = (0, 0, 2 000). The relative errors of displacement field in the norm are listed in Table 3. Figure 6 shows the numerical results of displacement solutions u₀3, u₁₃^ε, u₂₃^ε and ue3 at the intersection z = 0.5, the calculating time and the load time are t = 30 and τ = 0, respectively. From Table 3 and Fig. 6, it is evident that the first-order displacement solution does not have a good agreement with the finite element solution based on the fine mesh, when the exact solution of the original problem has a large oscillation behavior, but the second-order solution can catch the oscillation behavior effectively. All the results obtained from the above numerical examples demonstrate that the SOTS method presented in this paper is effective to evaluate the ageing linear viscoelastic problem for the composite materials with the periodic structure.

4.3 Application of SOTS method

In this subsection, we use the SOTS method to predict the creeping properties of concrete. Concrete can be considered as a composite material which consists of the cement matrix and aggregate, and the aggregate particles are distributed in the cement randomly. Yu et al.^[20] have given an effective computational method to generate the unit cell with random particles, which is used in this example to generate the geometric model of concrete. The aggregate material in this example is granite, and it is considered as an elastic material. The mechanical properties of granite aggregate are given by E = 43.5 GPa and Poisson’s ratio γ = 0.2. The cement matrix usually has the viscoelastic properties, and the relaxation modulus can be expressed as G(t,τ) = E( )K(t,τ). In this example, E(τ) and K(t,τ) have the following expressions:

and Poisson’s ratio γ = 0.2 for the cement matrix. The numerical results predicted by the SOTS method for the effective modulus of concrete are given in Fig. 7, and the volume fraction of granite aggregate is 20%. Figure 8 shows the influence of different load time τ on the effective modulus. Further, we set a 12 MPa stress load on the upper surface at τ= 3 d. The strain and stress distributions in the concrete are calculated, and the result in a local cell for t = 60 d is shown in Fig. 9. These numerical results show that the SOTS method in this paper is effective to predict the creeping properties of concrete.

5 Conclusions

In this paper, the SOTS computational method is presented for the ageing linear viscoelastic prob- lem in composite materials with a periodic structure. The formulation of SOTS method is given, including the SOTS solutions of displacement field and the homogenization relaxation modulus. The error analysis of the two-scale approximate solutions is briefly given, and the second-order solution has a better approximation to the solution of the original problem compared with the first-order solution theoretically. By evaluating the second-order cell function N_α1α2 , we can obtain more detailed information of the displacement, strain, and stress, which indicates that the SOTS method can provide a better result. The numerical examples also show that the SOTS method is effective to solve the ageing linear viscoelastic problem in composite materials.