Toward Analytical Homogenized Relaxation Modulus for Fibrous Composite Material with Reduced Order Homogenization Method

1  Introduction

Fiber-reinforced composites offer advantages such as lightweight properties, high strength, extended maintenance intervals, and corrosion resistance, making them widely used in the aerospace and automotive industries [1–5]. Thermosetting or thermoplastic matrices with cylindrical microstructures exhibit transversely isotropic viscoelastic behavior, as documented in the literature [6,7]. The overall mechanical behavior of fibrous composites is analyzed using multiscale methods. These methods provide detailed insights into the underlying mechanisms, enabling material performance optimization and guiding innovative material design.

Initially, homogenization methods for viscoelastic composites were based on the correspondence principle. Micromechanical theoretical schemes developed for elastic composites were extended to viscoelastic composites to estimate their effective elastic moduli in the Laplace domain. Various approaches, including mean-field homogenization methods [8–10], the asymptotic homogenization method [11–13], the variational principle [14,15], and the sequential linearization approach [16], have been employed to simulate viscoelastic composites. Results obtained in the Laplace domain are then transformed back to the time domain to determine the effective response of the viscoelastic materials. However, performing an analytical inverse transform is usually too complex, so the macroscopic response in the time domain can only be obtained through numerical inversion [12] or approximate inversion methods [17], which can be costly and less accurate.

To overcome the limitations of the Laplace transformation, various homogenization methods have been developed to operate directly in the time domain [18–20]. Mean-field homogenization methods, initially formulated for linear elastic composite materials, have since been extended to address nonlinear behaviors, such as elastoplastic [21,22] and viscoelastic [23] responses, among others. These extensions are typically implemented by linearizing local constitutive laws incrementally, an approach necessitated by the Eshelby problem’s assumption that mechanical interactions within phases are purely elastic. Additionally, stress and strain fields within each phase are approximated as phase averages, treating them as homogeneous. However, in nonlinear problems, where stress and strain distributions are inherently non-uniform across phases, this approximation can lead to an over-stiffness issue. To address this challenge, techniques such as the isotropization method [24] and higher-order theories [25] have been proposed. Bleiler et al. [26] introduced a tangent second-order homogenization method to estimate the effective mechanical response of fibrous composites with hyperelastic, anisotropic, and incompressible phases. Furthermore, Barral et al. [27] advanced the Mori-Tanaka method by incorporating a modified transformation field analysis (TFA) approach first proposed by Dvorak and his coworkers [28] and introducing a specialized coating between fibers and the matrix, effectively mitigating over-stiffness in the analysis of fiber-reinforced composites. In the TFA method, the inelastic strains are considered as given eigenstrains. This method assumes that eigenstrains in each phase of the composite are piecewise constant, this lower-order approximation can lead to inaccurate solutions for overall properties. To overcome the problem, Michel et al. [29,30] have introduced the nonuniform transformation field analysis, where full field calculations take place and the plastic behavior is traced with the help of several plastic modes. Strategies such as subdividing each phase into numerous partitions, using asymptotically consistent eigenstrain fields have also been employed [31].

The asymptotic homogenization method is a mathematically rigorous methodology developed by Sanchez-Palencia [32], Bensoussan and Alain [33]. Compare to the Mean-field homogenization method, it effectively captures the interactions between different phases within composite materials, enabling more precise predictions of their overall performance. It has been proved that the considered fields converge towards the homogeneous macroscopic solution as the micro-structural parameter ξ=x/y tends to 0. The classic asymptotic homogenization method can obtain the linear homogenized modulus by solving a linear unit cell problem with given phase material elastic properties. Thus, one can use the homogenized elastic properties in a macroscopic simulation without downscaling or upscaling between the macro- and meso-scales. In this manner, the linear homogenization process is named as “off-line” process. However, if material nonlinearities at the mesoscale are considered, a unit cell problem must be solved at every integration point of all the macroscopic elements during each macroscopic Newton-Raphson iteration. Additionally, the unit cell problem itself requires a Newton-Raphson iterative process due to material nonlinearities. This repeating process is known as the “on-line” process since it depends on the “on-the-fly” kinematic measurements at the macro-scale, significantly increasing the computational cost. Yu et al. [34] used the asymptotic homogenization method to solve thermo-viscoelastic problems, considering both multiple spatial and temporal scales. Zhai et al. [18] developed a time-domain asymptotic homogenization method based on the integral form of the Kelvin–Voigt viscoelastic model and used it for woven fabric composites.

Multiscale computational techniques based on the finite element method such as FE2, offer the possibility of computing the macro-structural response of heterogeneous materials with an arbitrary microscopic geometry and constitutive behavior. A scheme based on the discrete homogenization method are developed to predict the effective mechanical properties of 3D dry textiles [35]. Recently, the isogeometric analysis method (IGA), proposed by Hughes and his coworkers [36,37], has been incorporated into computational homogenization [38–40]. The method aims to unify the representation of geometric models and mesh models by directly employing the spline functions as the shape functions. The errors produced by geometric approximations inherent in conventional finite-element and finite-volume techniques could be greatly reduced because the same model is used for both modeling and analysis. Cylindrically orthotropic [41] and viscoelasticity [42] are considered for fiber-reinforced composites using IGH methods.

To decrease computational costs at the meso-scale without significantly compromising solution accuracy, the reduced-order-homogenization (ROH) multiscale method was developed, combining the asymptotic homogenization method with transformation field analysis (TFA) [43,44]. The ROH method was first used in composite plasticity and failure analysis [43,45], and higher-order ROH methods have been developed for three-scale nonlinear problems [46–49]. In this paper, a two-phase model with the assumption of one partition per phase is employed for fiber-reinforced composites. Due to the unique structure of the fiber-reinforced unit cell, analytical transformation influence functions can be derived. Based on ROH theory, an analytical homogenization constitutive model for unidirectional fibrous composites is proposed, eliminating the need for computational homogenization. As demonstrated in the following section, the proposed model does not exhibit significant over-stiffness and helps reduce computational time.

This paper is organized as follows. In Section 2, we introduce the analytical viscoelastic fibrous composites constitutive model based on the ROH homogenization method. Special cases where the fiber volume fraction approaches zero are considered to confirm the model’s consistency under extreme conditions. We investigate the impact of material parameters of each component on the relaxation times of the homogenized viscoelastic material. Additionally, numerical implementation specifics for the proposed model are detailed in this section. Section 3 presents numerical verifications of the proposed model through three sets of numerical examples: unit cell tests, three-point beam bending tests, and torsion tests. Experimental results from the literature are used to validate the proposed constitutive model. Conclusions are provided in Section 4.

2  Methodology

In this section, we introduce an analytical framework to determine homogenized viscoelastic material properties for the fibrous composite material based on the so-called reduced order homogenization (ROH) method. Here we assume the fiber phase follows transversely isotropic symmetry and keeps elastic that

σijf=Lijklfεklf(1)

On the other hand, the matrix phase is viscoelastic and described by the standard solid model [50], and the stress response is written as the following form:

σijm=Lijklm,∞εklm+∫0te−(t−η)/τ1Lijklm,isoε˙klm(η)dη(2)

where

Lijklm,∞=2μm,∞Iijkldev+3κm,∞Iijklvol,μm,∞=Em,∞2(1+νm),κm,∞=Em,∞3(1−2νm)(3a)

Lijklm,iso=2μisomIijkldev+3κisomIijklvol,μisom=Eisom2(1+νm),κisom=Eisom3(1−2νm)(3b)

Here Em,∞ is the equilibrium Young’s modulus of the matrix phase, and νm is the corresponding Poisson’s ratio. Eisom is the isotropic relaxation modulus. Finally, the fourth-order tensors Iijklvol and Iijkldev are defined as

Iijklvol=13δijδkl,Iijkldev=12(δikδjl+δilδjk)−13δijδkl(4)

respectively. Here the quantities associated with the fiber and the matrix phase are marked with superscription Ef and Em, respectively.

Our goal is to obtain an analytical homogeneous relaxation tensor Rijklc(t) that

σijc=Lijklc,∞εklc+∫0tRijklc(t−η)ε˙klc(η)dη=Lijklc,∞εklc+Rijklc(t)∗ε˙klc(t)(5)

based on given material properties of the fiber and the matrix phase without computational homogenization process, where σijc, εijc are the macroscopic stress and strain tensor, respectively, Lijklc,∞ is homogenized linear elastic stiffness tensor, and f(t)∗g(t) denotes the convolution of function f(t) and g(t). Here the macroscopic properties are marked with superscription εc.

2.1 Reduced Order Homogenization (ROH) for Fibrous Composite Material

The reduced-order homogenization (ROH) method is based on two key assumptions: 1) It is possible to distinguish between two length scales associated with macroscopic and microscopic phenomena; 2) The microstructure is sufficiently regular to be considered periodic; 3) The material within the same partition maintains a consistent status throughout, which reduces computational costs at the meso-scale. A complete derivation of the ROH method would be too extensive for this section, so we will present only the essential formulas. For a detailed explanation of the methodology, please refer to reference [44]. The fibrous composite material consists of the fiber and matrix phases. We introduce eigenstrains μij defined as [44]

μij=εij−Mijklσkl(6)

where Mijkl is the elastic compliance tensor. The ROH aims to solve the following system of nonlinear equations [44]:

εijf−Pijkl,fmμklm−Pijkl,ffμklf=Eijkl,fεklc(7a)

εijm−Pijkl,mmμklm−Pijkl,mfμklf=Eijkl,mεklc(7b)

where Eijkl,f and Eijkl,m are the elastic strain influence functions, Pijkl,fm, Pijkl,mm, Pijkl,ff and Pijkl,mf are the eigenstrain influence functions, and εijc is the macroscopic strain tensor.

Since the fiber phase keeps elastic, we can obtain μijf=0 directly. The eigenstrains of the matrix phase are given as

μijm=εijm−Mijklσijm=εijm−(Lijklm)−1σklm=εijm−(Lijklm,∞)−1σklm=−ψm∫0te−(t−η)/τ1ε˙ijm(η)dη(8)

by the standard solid model defined in Eq. (2), with consideration of material isotropy Eqs. (3a) and (3b), and define ψm≡Eisom/Em,∞. Here we define

Lijklm=Lijklm,∞(9)

By linear elastic stiffness tensor of the fiber phase and the matrix phase given in Eqs. (1) and (9), respectively, we can first obtain the homogenized elastic stiffness tensor Lijklc that

Lijklc=cfLijmnfEmnkl,f⏟≡Lijklcf+cmLijmnmEmnkl,m⏟≡Lijklcm(10)

where cf and cm=1−cf are the volume fraction for the fiber and the matrix phase, respectively. The macro-scopic averaged stress then is given by volumetric averaging of phase stresses that

σijc=cfLijklfεklf+cmLijklmεklm=Lijklcεklc+Aijklmμklm(11)

Aijklm=cfLijstfPstkl,fm+cmLijstm(Pstkl,mm−Istkl)(12)

Usually, the fibrous composite material is assumed to satisfy the transversely isotropic symmetry so that the linear elastic stiffness tensor can be written with Hill’s constants nc, lc, mc, kc, and pc [51], with the fibers aligned along the longitudinal x1−axis:

{Lijklc}=[nclclclckc+mckc−mclckc−mckc+mcmcpcpc](13)

The fiber phase also is transversely isotropic, where its Hill’s constants are written as nf, lf, mf, kf, and pf. For fibrous composite material, the following universal connections are satisfied [51]:

kc−kflc−lf=kc−(λm,∞+μm,∞)lc−λm,∞=lc−cflf−cmlmnc−cfnf−cmnm=kf−(λm,∞+μm,∞)lf−λm,∞(14)

where λm,∞ and μm,∞ are the Lamé constants of the matrix phase.

2.2 Analytical Influence Functions for Fibrous Composite Material

The elastic strain and eigenstrain influence functions are usually evaluated by solving a number of finite element problems over the given fibrous unit cell [44]. In this manuscript, we aim to provide an analytical method to evaluate those functions with given linear homogenized stiffness tensor Lijklc, which can be estimated by Mori-Tanaka or self-consistent method [51].

First of all, the elastic strain influence functions satisfy the following constraint that [44]

cfEijkl,f+cmEijkl,m=Iijkl(15)

One can obtain the elastic strain influence functions by Eqs. (10) and (15) [44]:

Eijkl,f=1cf(Lijmnf−Lijmnm)−1(Lmnklc−Lmnklm)(16a)

Eijkl,m=1cm(Lijmnm−Lijmnf)−1(Lmnklc−Lmnklf)(16b)

In other words, for a two-phase unit cell, as long as Lijklc is given, the elastic strain influence functions can be obtained directly.

Next, we continue to evaluate the eigenstrain influence functions. Define the following supportive tensors:

{E~ijpq,f}≡[100−νatf00−νatf00000](17)

and

{E~ijpq,m}≡[100cfνatf/cmf220cfνatf/cm0f33f23f13f12](18)

fij=Lijijc/Lijijcm,ij=22,33,23,13,12,no summation over repeated indices(19)

Pijkl,mm=Iijkl−Eijmn,m(E~mnkl,m)−1,Pijkl,f=Iijkl−Eijkl,m−Pijkl,mm(20a)

Pijkl,fm=(Eijmn,f−E~ijmn,f)(E~mnkl,m)−1,Pijkl,ff=Iijkl−Eijkl,f−Pijkl,fm(20b)

However, Pijkl,mf and Pijkl,ff are not used since μijf=0.

Due to the symmetry of fibrous unit cell, one can use only 5 non-zero parameters to define one influence tensor. For the linear strain influence tensors, we have (E=f,m)

[100E11,22E22,22E33,22E11,33E22,33E33,33E23,23E13,13E12,12]⇒[100E1E2E3E1E3E2E4E5E5](21)

and for the eigenstrain influence tensors, we also have (E=f,m)

[100P2211,mmP2222,mmP2233,mmP3311,mmP3322,mmP3333,mmP2323,mmP1313,mmP1212,mm]⇒[100P1mmP2mmP3mmP1mmP3mmP2mmP4mmP5mmP5mm](22)

From Eqs. (21) and (22), we can quickly obtain the following results:

ε11m=ε11f=ε11c(23)

which means the axial strains of the matrix and the fiber phases are the same and equals to the macroscopic axial strain ε11c.

2.3 Viscoelastic Model for the Fibrous Composite Material

In this section, we propose a homogenized viscoelastic model with the analytical elastic strain and eigenstrain influence functions. Substituting Eqs. (8) into (7b) with μijf=0 yields

εijm+Pijkl,mmψmm∫0te−(t−η)/τ1ε˙klm(η)dη=Eijkl,mεklc,ij=22,33,23,13,12,summation convention over k,l(24)

where ij=11 is not considered due to Eq. (23). In the rest of this section, we derive system of ordinary differential equations (ODEs) for ε˙ijm with ij=22,33,23,13,12. First, taking the derivative of Eq. (24) with respect to time t yields

ε˙ijm(t)+Pijkl,mmψmm[ε˙klm(t)−1τ1∫0te−(t−η)/τ1ε˙klm(η)dη]=Eijkl,mε˙klc(t)(25)

Thus we obtain

Pijkl,mmψm∫0te−(t−η)/τ1ε˙klm(η)dη=τ1[Iijkl+ψmPijkl,mm]ε˙klm(t)−τ1Eijkl,mε˙klc(t)(26)

Next, substituting Eqs. (26) into (24) yields

εijm+τ1[Iijkl+ψmPijkl,mm]ε˙klm(t)=Eijkl,m[εklc+τ1ε˙klc(t)],ij=22,33,23,13,12(27)

By solving the system of ordinary differential equations in Eq. (27), we obtain the matrix phase strain εijm, where ij=22,33,23,13,12. The detailed derivation is provided in Appendix A. The solutions for the normal strains ε22m(t) and ε33m(t) can be expressed as

ε22m(t)=E1mε11c(t)+E2mε22c(t)+E3mε33c(t)+[−E1mτ2−τ1τ2−τ1τ2ψmP1mm]e−t/τ2∗ε˙11c(t)+[−12(E2m+E3m)τ2−τ1τ2e−t/τ2−12(E2m−E3m)τ3−τ1τ3e−t/τ3]∗ε˙22c(t)+[−12(E2m+E3m)τ2−τ1τ2e−t/τ2+12(E2m−E3m)τ3−τ1τ3e−t/τ3]∗ε˙33c(t)(28a)

ε33m(t)=E1mε11c(t)+E3mε22c(t)+E2mε33c(t)+[−E1mτ2−τ1τ2−τ1τ2ψmP1mm]e−t/τ2∗ε˙11c(t)+[−12(E2m+E3m)τ2−τ1τ2e−t/τ2+12(E2m−E3m)τ3−τ1τ3e−t/τ3]∗ε˙22c(t)+[−12(E2m+E3m)τ2−τ1τ2e−t/τ2−12(E2m−E3m)τ3−τ1τ3e−t/τ3]∗ε˙33c(t)(28b)

The solutions for the shear strains

ε23m(t)=E4mε23c(t)−τ4−τ1τ4E4me−t/τ4∗ε˙23c(t)(29a)

ε13m(t)=E5mε13c(t)−τ5−τ1τ5E5me−t/τ5∗ε˙13c(t)(29b)

ε12m(t)=E5mε12c(t)−τ5−τ1τ5E5me−t/τ5∗ε˙12c(t)(29c)

respectively. With

τ2=τ1(1+ψmP2mm+ψmP3mm),τ3=τ1(1+ψmP2mm−ψmP3mm),τ4=τ1(1+ψmP4mm),τ5=τ1(1+ψmP5mm)(30)

here τk,k=2,3,4,5 are characteristic relaxation times for the fibrous composite materials.

So far, we are able to find the matrix phase strain εijm by the macroscopic strain εijc and the corresponding strain rate ε˙ijc. In order to obtain an analytical result for macroscopic stress defined in Eq. (11), we continue to derive the eigenstrain of the matrix phase. For μ11m, we can directly obtain

μ11m=−ψm∫0te−(t−η)/τ1ε˙11m(η)dη=−ψm∫0te−(t−η)/τ1ε˙11c(η)dη=−ψme−t/τ1∗ε˙11c(t)(31)

by considering Eq. (23). For the rest components, expand Eq. (7b) into

[P2mmP3mmP3mmP2mm]{μ22m(t)μ33m(t)}=[ε22m(t)−E1mε11c(t)−E2mε22c(t)−E3mε33c(t)−P1mmμ11mε33m(t)−E1mε11c(t)−E3mε22c(t)−E2mε33c(t)−P1mmμ11m](32a)

P4mmμ23m(t)=ε23m(t)−E4mε23c(t)(32b)

P5mmμ13m(t)=ε13m(t)−E5mε13c(t)(32c)

P5mmμ12m(t)=ε12m(t)−E5mε12c(t)(32d)

Finally, one can introduce a fourth-order tensor Yijkl,m(t) to represent the eigenstrain μijm(t) as the following:

μijm(t)=Yijkl,m(t)∗ε˙klc(t)(33)

First, from Eq. (31), we can obtain the components corresponding to μ11m that

Y1111,m(t)=−ψme−t/τ1(34)

From Eqs. (28a) and (28b), we obtain

Y2211,m(t)=Y3311,m(t)=P1mmψmP2mm+P3mme−t/τ1+[−τ2−τ1τ2E1mP2mm+P3mzm−τ1τ2P1mmψmP2mm+P3mm]e−t/τ2(35a)

Y2222,m(t)=−τ2−τ1τ2E2m+E3m2(P2mm+P3mm)e−t/τ2−τ3−τ1τ3E2m−E3m2(P2mm−P3mm)e−t/τ3(35b)

Y3322,m(t)=−τ2−τ1τ2E2m+E3m2(P2mm+P3mm)e−t/τ2+τ3−τ1τ3E2m−E3m2(P2mm−P3mm)e−t/τ3(35c)

Y2233,m(t)=Y3322,m(t),Y3333,m(t)=Y2222,m(t)(35d)

In addition, from Eq. (29b) and (29c), we obtain

Y2323,m(t)=−τ4−τ1τ4E4mP4mme−t/τ4,Y1313,m(t)=Y1212,m(t)=−τ5−τ1τ5E5mP5mme−t/τ5(35e)

and all other components of Yijkl are zero. In summary, the components of the eigenstrain μijm(t) are written as the following:

μ11m(t)=Y1111,m(t)∗ε˙11c(t)μ22m(t)=Y2211,m(t)∗ε˙11c(t)+Y2222,m(t)∗ε˙22c(t)+Y2233,m(t)∗ε˙33c(t)μ33m(t)=Y3311,m(t)∗ε˙11c(t)+Y3322,m(t)∗ε˙22c(t)+Y3333,m(t)∗ε˙33c(t)μ23m(t)=Y2323,m(t)∗ε˙23c(t)μ13m(t)=Y1313,m(t)∗ε˙13c(t)μ12m(t)=Y1212,m(t)∗ε˙12c(t)(36)

Substituting Eqs. (33) into (11) yields

σijc=Lijklcεklc+AijmnmYmnkl,m(t)⏟Rijklc(t)∗ε˙klc(t)(37)

2.4 Model Consistency in an Extreme Case

In this section, we exam the homogenized viscoelastic properties of the fibrous composite material by letting cf=0 and cm=1. Under this situation, the homogenized viscoelastic properties should approach the viscoelastic response of pure matrix phase as well.

By setting cm=1, one obtains Lijklc=Lijklcm=Lijklm through Eq. (10) and Eijkl,m=Iijkl. In addition, we have Aijklm=−Lijklm by Eq. (12). The value of Eijkl,f does not have any meaning. Similarly, by Eq. (18), we obtain E~ijpq,m=Iijkl and thus leads to Pijkl,mm=0 by Eq. (20a). This result guarantees that τk=τ1, k=2,3,4,5 by Eq. (30).

Considering the aforementioned conditions, we can obtain εijm=εijc from Eqs. (23), (28a) and (29a), by applying ε~22m(t)=ε~33m(t)=0, E2m=E4m=E5m=1, E1m=E3m=0. Finally, one can obtain Yijkl,m=−ψmIijkl and

σijc=σijm=Lijklm,∞εklm+ψmLijklm∗ε˙klm(t)=Lijklm,∞εklm+Lijklm,iso∗ε˙klm(t)(38)

which traces back to the pure matrix case.

2.5 Influence of Component Parameters on Homogenization Relaxation Times

In this section, we analytically calculate the homogenization relaxation times τi,i=1,2,…,5 using known component parameters, investigating how fiber and matrix properties influence the overall viscoelastic response. It is important to note that relaxation time is closely related to material viscosity. A smaller viscosity results in a shorter relaxation time. When viscosity (η) approaches zero, the material behaves like an inviscid fluid, whereas when it approaches infinity, the material behaves like an elastic solid. We observe that changes in component parameters lead to variations in the homogenization relaxation times, with certain parameters exerting more significant influence compared to others.

We first investigate the influence of fiber parameters on homogenized relaxation times. Parameters are chosen as detailed in Table 1, with one parameter varied within specified ranges. Fig. 1a illustrates the impact of fiber volume fraction on viscoelastic behavior: as the fiber fraction increases, the relaxation times also increase, indicating that, in our research, fibrous composites with higher fiber content tend to exhibit more solid-like characteristics in the homogenized state. It is noteworthy that relaxation time associated with the fiber direction is always equal to the matrix relaxation time, as indicated by Eq. (36), implying that the viscoelastic behavior in this direction is solely influenced by matrix components, while viscoelastic behavior in other directions is influenced by both fiber and matrix components. Fig. 1b demonstrates how homogenized relaxation time varies with changes in νtf=0.20,0.25,0.35,0.40, revealing that νtf has a limited impact on viscosity behavior. A slight decrease in homogenized relaxation time is observed with increasing νtf.

Figure 1: Influence of fiber parameters on homogenized relaxation times: (a) Fiber volume; (b) Poisson ratio in the transverse plane νtf

Next, we explore the influence of matrix parameters on homogenized relaxation times. As shown in Fig. 2a, an increase in Eisom leads to higher homogenized relaxation times due to the corresponding increase in viscosity. Similarly, Fig. 2b demonstrates that higher E∞m results in increased homogenization relaxation times, paralleling the trend observed in Fig. 1a. Fig. 2c indicates that νm has minimal impact on viscosity behavior, albeit showing a slight increase with higher νm.

Figure 2: Influence of matrix parameters on homogenized relaxation times: (a) Relaxation modulus Eisom; (b) Elastic modulus E∞m; (c) Poisson ratio νm

2.6 Numerical Algorithm

In this section, we briefly introduce the numerical algorithm for the homogenized viscoelastic model of fibrous composite material. In finite element analysis, we update material stress as well as consistent tangent modulus over one increment [tn,tn+1]. The stress at tn+1, namely, σijc,n+1 is given as

σijc,n+1=Lijklcεklc,n+1+Aijklmμklm,n+1=Lijklcεklc,n+1+AijklmYklst,m(t)∗ε˙stc,n+1(t)=Lijklcεklc,n+1+Rijklc(t)∗ε˙klc,n+1(t)(39)

again Rijklc(t)≡AijmnmYmnkl,m(t) is the homogeneous relaxation tensor. The second term in the right-hand-side of Eq. (39) is given by the classic integration form:

AijklmYklst,m(t)∗ε˙stc,n+1(t)=Aijklm∫0tn+1Yklst,m(tn+1−η)ε˙stc(η)dη=Aijklm∫0tnYklst,m(tn+1−η)ε˙stc(η)dη+Aijklm∫tntn+1Yklst,m(tn+1−η)ε˙stc(η)dη(40)