Dimitrios Anagnostou
email: dimitrios.anagnostou@univ-avignon.fr

George Chatzigeorgiou
email: georges.chatzigeorgiou@ensam.eu

Yves Chemisky
email: yves.chemisky@ensam.eu

Fodil Meraghni
email: fodil.meraghni@ensam.eu

Hierarchical Micromechanical Modeling of the Viscoelastic Behavior Coupled to Damage in SMC and SMC-Hybrid Composites

Keywords: Micromechanics, Viscoelasticity, Inclusion Method, Multiscale Modeling, Micro-cracks

The aim of this paper is to study, through a multiscale analysis, the viscoelastic behavior of glass reinforced sheet molding compound (SMC) composites and SMC-hybrid composites mixing two types of bundle reinforcement: glass and carbon fibers. SMC exhibit more than two distinct characteristic length scales, so that a sequence of scale transitions is required to obtain the overall behavior of the composite. An analytical procedure is used consisting of properly selected well-established micromechanical methods like the Mori-Tanaka (MTM) and the composite cylinders (CCM) accounting for each scale transition. After selecting a representative volume element (RVE) for each scale, the material response of any given length scale is described on the basis of the homogenized behavior of the next finer one. This hierarchical approach is appropriately extended to the viscoelastic domain to account for the time dependent overall response of the SMC composite material. The anisotropic damage has been introduced through a micromechanical model considering matrix penny-shape microcrack density inside bundles. The capabilities of the hierarchical modeling are illustrated with various parametric studies and simulation of experimental data for glass-based SMC composites.

Introduction

Sheet molding compounds (SMC) are fiber-reinforced thermoset polymer matrix semi-finished products obtained by a thermocompression process [START_REF] Whelan | Molding of thermosetting plastics[END_REF][START_REF] Mallick | Fiber-reinforced composites: materials, manufacturing, and design[END_REF][START_REF] Orgéas | Sheet molding compounds[END_REF][START_REF] Lamanna | Mechanical characterization of sheet moulding composites for the automotive industry[END_REF][START_REF] Schladitz | Non-destructive characterization of fiber orientation in reinforced SMC as input for simulation based design[END_REF][START_REF] Shirinbayan | Coupled effect of loading frequency and amplitude on the fatigue behavior of advanced sheet molding compound (A-SMC)[END_REF]. Since the last two decades, these composites represent an ideal choice for large structural automotive components due to their high strength-to-weight ratio. SMC composites are considered suitable materials for many engineering applications due to their enhanced properties. They have better mechanical performance compared to thermosetting compounds like bulk molding compounds (BMC) and

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

fiber-reinforced thermoplastics. The low weight and low cost together with mechanical characteristics and versatility allow to reach cost-efficient performances in many weight-critical applications. Because of advantages like the above, SMC composites are employed in the manufacturing of many functional or structural components in a wide range of industrial applications: automotive and truck industry, building constructions, electrical/electronics, transportation, aerospace, chemical engineering, medical sectors and marine industry.

SMC composites usually include short fibers or fiber bundles impregnated in the matrix phase. Traditionally, SMC are reinforced with glass fiber bundles. However, various types of fibrous reinforcement can be utilized: glass fiber bundles, carbon bundles, or hybrids. The hybrid SMC composites, containing both glass and carbon bundles, is a modern idea that very recently attracted the composite suppliers industry. The hybridization process can consider glass and other types of fibers, like carbon [START_REF] Palmer | Sheet moulding compound (SMC) from carbon fibre recyclate[END_REF][START_REF] Zaiß | Use of thermography and ultrasound for the qualuty control of SMC lightweight material reinforced by carbon fiber tapes[END_REF][START_REF] Trauth | Mechanical properties and damage evolution of a structural sheet molding compound based on a novel two step curing resin system[END_REF] or natural ones [START_REF] Sanjay | Studies on Natural/Glass Fiber Reinforced Polymer Hybrid Composites: An Evolution[END_REF]. The usual SMC glass fiber reinforced composites have been extensively studied and have practically reached their limit in terms of mechanical performance. On the other hand, the carbon fibers, which are generally stronger, are very expensive and lead to quasi-brittle response. Thus, a potential solution could be to mix glass and carbon fibers in low proportion in order to increase the mechanical performance, retain ductility and damage tolerance, while keeping also reasonable the production cost.

Modeling the overall mechanical behavior of SMC composites has received increasing attention in the past years. [START_REF] Carman | Micromechanics of short-fiber composites[END_REF] reviewed earlier investigations and developed a method where, at a first approximation, short fibers may be characterized as continuous fibers that have undergone fiber fracture. The approximation was improved by considering load transfer at the fiber ends. [START_REF] Fitoussi | A statistical micromechanical model of anisotropic damage for smc composites[END_REF] applied a statistical micro-macro relationship with the help of Mori and Tanaka's model and showed that for SMC composites the predominant damage is the debonding at the fiber/matrix interfaces. [START_REF] Derrien | Prediction of the effective damage properties and failure properties of nonlinear anisotropic discontinuous reinforced composites[END_REF] studied the tensile behavior of SMC using a micro-macro relationship taking into account the damage mechanisms. [START_REF] Morozov | Progressive damage modelling of smc composite materials[END_REF] studied the damage behavior of SMC composites using a finite element package. Next, [START_REF] Jendli | Micromechanical analysis of strain rate effect on damage evolution in sheet molding compound composites[END_REF] carried out an experimental study in conjunction with the application of well-known damage mechanics techniques in order to investigate the strain rate influence on the mechanical behavior and the damage growth of a discontinuous-random fiber composite (SMC-R26) subjected to a high-speed tensile loading. In addition, the same authors developed later [START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre smc composites[END_REF] a micromechanical model which relies upon an experimental methodology performed according to an incremental operating strategy in order to predict the constitutive law for SMC-R26. [START_REF] Larbi | Constitutive model of micromechanical damage to predict reduction in stiffness of a fatigued smc composite[END_REF] investigated the elastic behavior of SMC with a provided orientation distribution function (ODF) of the fibers acted upon by a cyclic loading. Furthermore, [START_REF] Teodorescu | An original approach of tensile behaviour and elastic properties of multiphase preimpregnated composite materials[END_REF] developed an approach based on an assumed periodicity of SMC and provided bounds of the effective elastic moduli. [START_REF] Teodorescu-Draghicescu | Homogenization and averaging methods to predict elastic properties of pre-impregnated composite materials[END_REF] applied a homogenization method and some averaging methods to compute the upper and lower limits of the homogenized elastic coefficients. [START_REF] Huang | Bridging and toughening of short fibers in smc and parametric optimum[END_REF] developed an analytical approach to predict the bridging and toughening of randomly oriented short fibers, which they use for the analysis of fiber bridging/toughening of SMC and the investigation of the effect of fiber parameters.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Most of the aforementioned works about SMC composites ignore the fact that SMC usually consist of fiber bundles reinforcements. [START_REF] Mulligan | Fibre-bundling in a short-fibre composite: 1. review of literature and development of a method for controlling the degree of bundling[END_REF] provides a comprehensive study of fiber bundling in short-fiber composites like SMC and highlight their advantages. In the present work, it is considered that a typical SMC-R composite is reinforced with randomly oriented bundles of (carbon or/and glass) fibers. In particular, we investigate the viscoelastic effective properties of SMC composites which consist of substitute (filled) matrix and planar randomly oriented bundles of randomly oriented glass fibers. In the case of SMC-hybrid composite, there are both glass and carbon fiber-bundle reinforcements.

In order to facilitate the analysis, a hierarchical approach is followed. It enables obtaining results for the effective moduli of each homogenization step (i.e. filled matrix, glass and carbon bundles) and the final time dependent effective medium as well. The modeling of such a multiscale system starts at the microscale by first considering the particulate reinforced composite of bulk matrix and fillers. In this step the classical Mori-Tanaka method (MTM) is applied [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF][START_REF] Chen | Mori-tanaka estimates of the overall elastic moduli of certain composite materials[END_REF] for obtaining the filled matrix with certain effective properties different than those of the bulk matrix. Another homogenization is also performed at the microscale, where the filled matrix is reinforced by a certain volume fraction of aligned glass (carbon) fibers. Here the composite cylinders method (CCM) is utilized [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF][START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Christensen | Mechanics of composite materials[END_REF][START_REF] Christensen | A critical evaluation for a class of micro-mechanics models[END_REF]. This second step of homogenization yields the effective properties of a glass (carbon) bundle. Next, homogenized bundles are laid up in random orientation in the filled matrix. A third step of homogenization at the mesoscale, using the generalized MTM that takes into account the random orientation of the inhomogeneities, provides the overall behavior of the final effective medium. In the third step additional nonlinear mechanisms can be integrated, like, for instance, damage at the bundles.

The paper is organized as follows: In Section 2 the problem under discussion is stated. The main objective of the present theoretical study is the micromechanical elastic/viscoelastic description of a SMC/SMC-hybrid composite exhibiting a hierarchical structure. The constituent materials, as well as their properties and geometrical aspects (i.e. aspect ratios, bundle orientation, etc.) are addressed. Further, the various associated volume fractions are defined.

In Section 3 the theoretical background for the constitutive response of a viscoelastic material is presented and two approaches are examined: the integral description with the dynamic correspondence principle and the differential description with incremental formulation. The micromechanics modeling in viscoelasticity, consisting of a sequence of scale transitions and homogenizations, is discussed in Section 4. A brief review of the classical Mori-Tanaka method with its extension in the case of random orientation and the Composite Cylinders Method is provided. The integration of damage mechanism at the scale of the bundles, based on a hybrid phenomenological-micromechanical approach is also discussed. In Section 5 the study is accompanied by various examples of parametric analyses and a simulation of real experimental data for SMC glass fiber reinforced composites. For the numerical applications, available creep data from the literature were handled properly for obtaining the Prony series representation of the viscoelastic matrix response. The conclusions of this work are summarized in the final section. The paper is accompanied by two Appendices.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Notation

Vectors will be denoted by lower-case Roman letters, second-order tensors by lower-case Greek letters, and fourthorder tensors by capital Roman letters. Whenever possible, vectors and tensors are written as boldface characters.

Indicial notation is used where necessary. The symbols I and I represent the second order and the fourth order symmetric identity tensors respectively, whose elements are given as

I i j = δ i j , I i jpq = 1 2 δ ip δ jq + δ iq δ jk , (1)
where δ i j is the Kronecker delta symbol defined by δ i j = 1 for i = j and δ i j = 0 for i j. The following notation will represent standard tensor algebraic operations.

Dot product:

σ i j n j = σ • n Double dot product:              L i jkl ε kl = L:ε L i jkl T klpq = L:T (2)
The generalized Hooke's law for a linear elastic solid is written as

σ = L:ε, ε = M:σ, (3)
where L and M are the stiffness and compliance tensors respectively and σ and ε are the Cauchy stress and infinitesimal strain tensor respectively. The inverse of a fourth order tensor A that possesses minor symmetries (A i jkl = A i jlk =

A jikl ) is defined as the tensor A -1 with the property

A: A -1 = A -1 : A = I . (4)

Statement of the Problem

The objective of the present theoretical study is to propose a micromechanical viscoelastic hierarchical modelling of a SMC-Hybrid composite, see Figure 1 a . The composite consists of randomly oriented bundles made of unidirectional aligned glass and carbon fibers. The bundle reinforcement takes place at the mesoscale of the material. It should be noted that each bundle contains one kind of fiber (that is, either glass or carbon), see Figure 1 b . The polymer matrix, itself, is reinforced with randomly distributed particles (calcium carbonate CaCO 3 and/or glass spherical fillers), see Figure 1 c . This reinforcement takes place at the microscale.

As it is described in more detail later, an hierarchical approach is followed enabling the estimation of the effective moduli at each homogenization step (i.e., filled matrix, glass and carbon bundles) and the overall viscoelstic effective medium as well (i.e., SMC-Glass, SMC-Carbon and SMC-Hybrid composite). In particular, modeling of such a system starts at the microscale, by first considering the particulate reinforced composite of bulk matrix and fillers.

This provides a filled matrix with certain effective properties different that those of the non-filled matrix. A second In what follows, a subscript A denotes axial (i.e. fiber) direction and a subscript T denotes transverse (i.e. parallel to the plane of rotational symmetry) direction. Regarding the labeling of the Poisson's ratios, the first subscript indicates the direction of the uniaxial stress, whereas the second one indicates the direction normal to the loading direction, in which the contraction takes place. In this way, ν AT is the Poisson's ratio corresponding to the contraction in the transverse direction that accompanies the uniaxial stress in the fiber direction. Further, a subscript or superscript Gf (Cf) indicates a variable associated with the glass (carbon) fibers, whereas a subscript or superscript Gb (Cb)

indicates a variable associated with the glass (carbon) bundles. A superscript m (p) indicates a variable associated with the matrix (particles; i.e. fillers). Additional notation will be introduced in Sections 3 and 4.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Involved Constituent Materials and their Properties

Polymer resin is a linear viscoelastic isotropic material. Its mechanical response can be characterized by the timedependent bulk and shear moduli, K m (t) and G m (t), respectively. As it is discussed later, it is possible to avoid time domain analysis by utilizing the so-called dynamic correspondence principle [START_REF] Fung | Foundations of solid mechanics[END_REF][START_REF] Hashin | Analysis of composite materials-a survey[END_REF][START_REF] Brinson | Polymer engineering science and viscoelasticity[END_REF]. The matrix phase may be reinforced by fillers of calcium carbonate and glass. Both of these particles are isotropic elastic materials. The glass fibers are also isotropic elastic materials whereas the carbon fibers are transversely isotropic elastic materials. The axis of rotational symmetry of the latter coincides with their geometrical axis of symmetry and their mechanical behavior can be characterized by five independent elastic moduli chosen among the axial and transverse Young's moduli, E C A and E C T , the axial and transverse shear moduli, G C A and G C T , the transverse (plane strain for lateral dilatation without axial extension) bulk modulus, K C T and the Poisson's ratios, ν C AT , ν C TA and ν C TT . Damage mechanisms in this type of composites can occur at different levels (matrix, fibers, fiber bundles, interfaces) and they can be taken into account with various ways [START_REF] Fitoussi | A statistical micromechanical model of anisotropic damage for smc composites[END_REF][START_REF] Derrien | Prediction of the effective damage properties and failure properties of nonlinear anisotropic discontinuous reinforced composites[END_REF][START_REF] Morozov | Progressive damage modelling of smc composite materials[END_REF][START_REF] Jendli | Micromechanical analysis of strain rate effect on damage evolution in sheet molding compound composites[END_REF][START_REF] Jendli | Multi-scales modelling of dynamic behaviour for discontinuous fibre smc composites[END_REF]. In this study, it is chosen to incorporate damage at the glass bundles using a recently developed hybrid phenomenological-micromechanical model [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF].

Geometrical Aspects

Here, the various geometrical aspects of the constituent materials are discussed. Starting with fillers, they are treated as spherical particles; Moreover, the length of the glass and carbon fibers is of mm order, whereas the diameter of the fiber is of µm. This results in an extremely large aspect ratio, i.e. they can be considered as long fibers inside a bundle. On the other hand, the bundles can be seen as "short reinforcements" embedded in the filled matrix.

The fiber-bundle orientations inside the filled matrix is described through an orientation distribution function (ODF), which can be implemented in the model as 2-D or 3-D function. In the examples of Section 5, a planar random bundles distribution is assumed using a uniform distribution function representing the equal probability of orientation. This choice is justified by the manufacturing process, in which all the fibers are mainly distributed in thermocompression plane. In fact, this type of process may induce very low out of plane reinforcement orientation, that can be neglected, notably, for low plate thickness.

Associated Volume Fractions

Crucial to the development of the parametric studies are the various volume fractions involved at each material scale. In the fabrication process the pure resin is reinforced by CaCO 3 and/or glass fillers. The symbol c p is reserved for this magnitude. Next, the resulting filled matrix is further reinforced either by glass and carbon bundles of fibers in the case of a SMC-hybrid composite or by glass bundles of fibers only in the case of a pure SMC composite.

At this point various volume fractions can be defined. First, the total fiber reinforcement volume fraction is denoted by c R and is defined by the total volume of fiber reinforcement (i.e. glass and carbon fibers in the case of SMC-hybrid) per total volume of the composite. This volume fraction is equal to the sum of the glass fibers volume fraction plus

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
carbon fibers volume fraction. The latter ones are defined, respectively, as the ratio of glass fiber reinforcement per total volume and carbon fiber reinforcement per total volume. The symbols c G f/t and c C f/t are reserved for these magnitudes. In addition, for the glass and carbon bundle, the volume fraction of glass and carbon fibers inside the corresponding bundle are introduced as c G f/b and c C f/b , respectively. Lastly, the volume fraction of the glass and carbon bundles in the total SMC-hybrid composite are introduced and they are denoted as c G b/t and c C b/t , respectively. The following scheme naturally explains the introduced volume fractions and their connections with the associated volumes:

c R = V G f + V C f V t = V G f V t c G f/t + V C f V t c C f/t =       V G f V G b       c G f/b       V G b V t       c G b/t +       V C f V C b       c C f /b       V C b V t       c C b/t = c G f/b c G b/t + c C f/b c C b/t (5)
Here V t is the total volume of the final composite, V G f and V C f are the volumes of the glass and carbon fibers respectively, and V G b , V C b are the volumes that the glass and carbon bundles occupy.

Theoretical Background for viscoelastic materials

Before describing the micromechanical framework, it is important to identify the constitutive laws of the various material constituents of SMC or SMC-hybrid composites. In these composites, at low or moderate stress levels the matrix material usually posseses viscoelastic behavior, while the particles and the fibers behave almost elastically.

Viscoelastic materials show hereditary behavior, i.e., their response at a given time depends on their previous load history. The constitutive equation for a linear viscoelastic material can be written in terms of Stieltjes integral equation of the form [START_REF] Christensen | Theory of Viscoelasticity, 2nd Edition[END_REF])

σ = t -∞ L(t -ξ): dε(ξ) dξ dξ, ( 6
)
where L is the time-dependent modulus of the material and σ and ε are the stress and strain tensors, respectively. In order to capture properly the material's response at a large range of frequencies, the modulus is often written in terms of a Prony series expanded form [START_REF] Brinson | Polymer engineering science and viscoelasticity[END_REF])

L(t) = L ∞ + N j=1 L j e -t/τ j , (7)
where L ∞ corresponds to the modulus tensor at very long (infinite) time, τ j are the relaxation characteristic times and L j are the Prony series moduli tensors1 . Two approaches can be considered for studying a linear viscoelastic material:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Dynamic correspondence principle

The constitutive response ( 6) can be represented in a linearized form by taking advantage of the half-sided Fourier or Laplace-Carson transform [START_REF] Fung | Foundations of solid mechanics[END_REF][START_REF] Hashin | Analysis of composite materials-a survey[END_REF][START_REF] Matzenmiller | Micromechanical modeling of viscoelastic composites with compliant fibermatrix bonding[END_REF][START_REF] Brinson | Polymer engineering science and viscoelasticity[END_REF].

The time-dependent stress and strain tensor can be generally written in terms of deviatoric (s, e) and spherical (trσ, trε) components,

σ = s + 1 3 trσ I, ε = e + 1 3 trε I, (8)
where tr denotes the trace of a second order tensor. Substituting equation ( 8) in equation ( 6), the constitutive law for a linear, isotropic viscoelastic material is expressed in terms of the time-dependent bulk modulus K(t) and shear modulus G(t) as follows:

trσ = t -∞ 3K(t -ξ) d(trε(ξ)) dξ dξ, s = t -∞ 2G(t -ξ) de(ξ) dξ dξ, K(t) = K ∞ + N j=1 K j e -t/τ j , G(t) = G ∞ + N j=1 G j e -t/τ j . ( 9
)
It is common practice in viscoelasticity to employ the dynamic correspondence principle in order to avoid time domain analysis of the problem. This tool permits the formulation of the original viscoelastic problem as a fictitious quasistatic linear elastic problem in terms of complex transformed properties in the transformed domain. The definition of the appropriate complex phase moduli are derived below for the case of isotropic phase materials. Imposing the displacement histories

u(t) = u(ω)e iωt , (10)
with ω denoting the frequency, strain fields can be determined using the strain-displacement relationship

ε = 1 2 ∇u + (∇u) T = ε(ω)e iωt . (11)
In this way, equations ( 9) are reduced to

s = 2iωG(ω) e(ω) = 2G * (ω) e(ω), trσ = 3iωK(ω) trε(ω) = 3K * (ω) trε(ω), (12)
where G and K are the half-sided Fourier transforms of the moduli, e.g.,

G (ω) = ∞ 0 G(t)e -iωt d t, K (ω) = ∞ 0 K(t)e -iωt d t. ( 13
)
The complex moduli G * and K * are given, respectively, by

G * (ω) ≡ iωG(ω) = G (ω) + iG (ω), K * (ω) ≡ iωK(ω) = K (ω) + iK (ω). (14)
A single-primed quantity indicates a storage modulus whereas a doubled-primed one indicates a loss modulus. The storage modulus is a measure of the energy stored and recovered by a viscoelastic material per cycle of sinusoidal

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT L ∞ L 1 H 1 L 2 H 2 L N H N ε ε v 1 ε v 2 ε v N Figure 2
: Rheological scheme of a spring and N Maxwell models (Maxwell branches) connected in parallel.

deformation, whereas the loss moduli is a measure of the energy dissipated as heat for a similar cycle [START_REF] Brinson | Polymer engineering science and viscoelasticity[END_REF]. Both storage and loss moduli are real quantities.

Taking the half-sided Fourier transform of equations ( 9) 3,4 and in view of equations ( 14) one obtains [START_REF] Brinson | Polymer engineering science and viscoelasticity[END_REF])

K (ω) = K ∞ + N j=1 K j ω 2 1 τ 2 j + ω 2 , K (ω) = N j=1 K j ω τ j 1 τ 2 j + ω 2 , G (ω) = K ∞ + N j=1 G j ω 2 1 τ 2 j + ω 2 , G (ω) = N j=1 G j ω τ j 1 τ 2 j + ω 2 . ( 15
)
While the linear form of the constitutive law in the space of frequencies permits fast calculations in a homogeneous viscoelastic material, in composites certain difficulties are imposed, related to the complex nature of the viscoelastic moduli. Indeed, in mean field theories like Mori-Tanaka, identifying complex Eshelby and concentration tensors is not a trivial task and most of the commercially available micromechanical softwares do not include such capability.

Incremental formalism

In linear viscoelasticity the integral formulation (6), combined with the Prony series representation of the viscoelastic moduli (7), can be represented through a rheological scheme that involves a sum of N Maxwell models (Maxwell branches) (Figure 2). Each Maxwell branch j (1 ≤ j ≤ N) consists of a linear elastic modulus L j and a 3-D dash-pot H j and is assigned with a viscoelastic strain tensor ε v j . For this rheological scheme the overall stress is provided by the equation

σ =         L ∞ + N j=1 L j         :ε - N j=1 L j :ε v j = L 0 :ε - N j=1 L j :ε v j , (16)
where

L 0 = L ∞ + N j=1
L j is the instantaneous modulus. Moreover, the evolution of each viscoelastic strain is given by

H j : εv j + L j :ε v j -L j :ε = 0, j = 1, ..., N. (17)
In isotropic linear viscoelasticity all L and H are isotropic fourth order tensors. When L j and H j are isotropic, and assuming that they have the same "Poisson's ratio" (which may be different from one Maxwell branch to another), one

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT obtains L -1 j : H j = τ j I .
Under such conditions, it can be shown that the formalisms ( 17) and ( 9) are equivalent (see for instance [START_REF] Lévesque | On the development and parameter identification of schapery-type constitutive theories[END_REF].

The evolution law (17) produces six equations with six unknowns (the viscoelastic strains) for each Maxwell branch, which need to be solved simultaneously. For calculation purposes, it is more convenient to introduce a viscoplasticity-like multiplier ṡ j for each branch which is equal to the norm of the viscoelastic strain rate εv j . Then ( 17) can be re-expressed in a plasticity-type manner

εv j = Λ j ṡ j , Λ j = H -1 j : L j :[ε -ε v j ] H -1 j : L j :[ε -ε v j ] , Φ j = H -1 j : L j :[ε -ε v j ] -ṡ j = 0, ( 18
)
where Λ j denotes the direction of the viscoelastic flow and Φ j plays the role of the viscoelastic activation criterion.

The symbol {•} denotes the usual norm of a second order tensor. Such formalism can be treated computationally using an incremental numerical scheme equivalent with those described in plasticity and viscoplasticity [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Qidwai | Numerical Implementation of a Shape Memory Alloy Thermomechanical Constitutive Model Using Return Mapping Algorithms[END_REF][START_REF] Achour | Implicit implementation and consistent tangent modulus of a viscoplastic model for polymers[END_REF][START_REF] Chatzigeorgiou | Thermomechanical Behavior of Dissipative Composite Materials[END_REF]. The advantage of these incremental methods is that they allow the identification of appropriate tangent modulus at each calculation increment, which is very useful in micromechanics schemes to estimate the overall response of a composite. Details about the numerical implementation of the linear viscoelastic model with N Prony series branches are given in Appendix A.

Multiscale framework for SMC and SMC-hybrid composites

A SMC-hybrid composite possesses more than two clearly distinct characteristic length scales and, as a result, a sequence of scale transitions is required in order to obtain its overall behavior through micromechanics methods. The whole procedure consists of various homogenizations accounting for each scale transition. In brief, after selecting a Representative Volume Element (RVE) for each scale, the material response of any given length scale is described on the basis of the homogenized behavior of the next finer one.

This hierarchical approach, which can be viewed as involving a sequence of scale transitions and homogenizations, is used in each step in order to give the "building blocks" at any level within hierarchical schemes. At each homogenization step well established methods are used, such as the Mori-Tanaka method (MTM) and the composite cylinders method (CCM). Figure 3 summarizes the procedure of the micromechanical modeling.

'Zeroth' Homogenization

In this step an effective medium from the homogenization of a composite material consisted of pure (viscoelastic) matrix and randomly distributed (elastic) spherical particles is determined, see Figure 3 a . It is assumed that no coating exists between the fillers and the matrix. Due to the existing isotropic response of both the matrix and the fillers and owing to the geometrical isotropy of the fillers (considered as spherical particles), the homogenized filled matrix exhibits an isotropic overall behavior.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

According to equations ( 12), there is a direct analogy between elasticity problems in the time domain and viscoelasticity problems in the frequency domain. Thus, it is straightforward to extend available micromechanics tools and methods of elasticity, such as the Mori-Tanaka solution, for composites with viscoelastic phases, by simply substituting complex moduli in place of elastic moduli in the final expressions [START_REF] Brinson | Comparison of micromechanics methods for effective properties of multiphase viscoelastic composites[END_REF][START_REF] Fisher | Viscoelastic interphases in polymer-matrix composites: theoretical models and finite-element analysis[END_REF].

As is well known [START_REF] Benveniste | On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media[END_REF][START_REF] Qu | Fundamentals of Micromechanics of Solidss[END_REF][START_REF] Lester | Transformation characteristics of shape memory alloy composites[END_REF][START_REF] Dvorak | Micromechanics of Composite Materials[END_REF][START_REF] Hossain | A multi-scale approach to model the curing process in magneto-sensitive polymeric materials[END_REF] the effective stiffness tensor of a composite material with M different types of inhomogeneities can be expressed in terms of the stiffness tensors of the constituent materials, their corresponding volume fraction and the global strain concentration tensor through the equation

L eff = L m + M r=1 c r L r -L m : A r , ( 19
)
where L eff denotes the effective (i.e. overall) stiffness tensor, L m is the stiffness tensor of the matrix phase, L r and c r are, respectively, the stiffness tensor and the volume fraction of the r-inhomogeneity (particle or fiber). Moreover, A r is the global strain concentration tensor which gives, essentially, the ratio between the average inhomogeneity strain and the corresponding average in the composite.

According to [START_REF] Fisher | Viscoelastic interphases in polymer-matrix composites: theoretical models and finite-element analysis[END_REF], when the constituents of a composite are viscoelastic, equation ( 19) can be formulated in the frequency domain as

L eff * = L m * + M r=1 c r L r * -L m * : A r * , (20)
where the superscript * denotes complex quantity in the form of (14).

To find the effective viscoelastic moduli of the filled matrix, the Mori-Tanaka two-phase model is utilized. From [START_REF] Qu | Fundamentals of Micromechanics of Solidss[END_REF] one has the following analytical closed-form expressions for the complex effective bulk and shear moduli (K FM * and G FM * respectively) of the filled matrix

K FM * = K m * + c p (K p * -K m * ) (3K m * + 4G m * ) 3K m * + 4G m * + 3 (1 -c p ) (K p * -K m * ) , (21)
G FM * = G m * + 5c p G m * (G p * -G m * ) (3K m * + 4G m * ) 5G m * (3K m * + 4G m * ) + 6 (1 -c p ) (G p * -G m * ) (K m * + 2G m * ) . ( 22
)
It is recalled that the superscript p denotes a quantity related with the particles. These complex moduli are transformed to the time (i.e. real) domain corresponding moduli with the help of the expressions ( 14), ( 15). In addition, the computation of the moduli of all Maxwell branches is extensively explained in the examples section.

Once the frequency to time domain transformation is performed, one can proceed to the next step that provides the effective properties of the homogenized bundles of glass fibers and/or the homogenized bundles of carbon fibers.

This homogenization, which also takes place at the microscale, is the subject of the following subsection.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

First Homogenization

In this step the overall elastic behavior of a homogenized bundle is obtained. The bundle is made of of aligned glass or carbon fibers reinforcing the filled matrix whose effective properties are provided from the 'zeroth' homogenization.

(Figure 3 b ). However, it is worth mentioning that the bundles usually have a fibers volume fraction that exceeds 70%. Due to this high fiber content, even though the matrix surrounding the fibers is viscoelastic, one could expect that the bundle exhibits an overall quasi-elastic behavior, that can be coupled to damage caused by microcracks.

Similar approximation has been considered in [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF] and the comparisons with experimental data showed excellent accuracy in the direction parallel to fibers, and quite satisfactory representation of the response normal to the fibers direction. This observation permits to simplify this step of homogenization, by assuming that, only for this case, the matrix also behaves elastically inside the bundle before damage occurs. The "elastic properties" of the filled matrix can be considered as the properties of its instantaneous response.

A. Elastic response of bundles

In order to obtain the effective elastic properties of the bundles, the Composite Cylinders Method (CCM for brevity) is utilized in the manner employed by, among others, [START_REF] Seidel | Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites[END_REF]; [START_REF] Seidel | Micromechanics modeling of the multifunctional nature of carbon nanotube-polymer nanocomposites[END_REF]; [START_REF] Chatzigeorgiou | Effective mechanical properties of fuzzy fiber composites[END_REF]. The length of the fibers is of mm order, whereas the diameter of the fibers is of µm. This results in an extremely large aspect ratio, i.e. they can be considered as long fibers.

It should be noted that the classical Mori-Tanaka scheme can provide also closed-form expressions in the case of long fibers. However, the CCM provides better estimations in the case of higher volume fractions. Also, the classical Mori-Tanaka is not applicable in the case of coated fibers and special methodologies are required [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF][START_REF] Cherkaoui | Elastic composites with coated reinforcements: a micromechanical approach for nonhomothetic topology[END_REF], in contrast to CCM where multiple coatings can be added in a straightforward manner. For these reasons, the CCM as the mean of homogenization in the current step has been selected here.

The homogenized bundles have a transversely isotropic behavior inasmuch as an isotropic (filled) matrix has been reinforced by aligned isotropic fibers (in the case of glass bundles) or by aligned transversely isotropic fibers (in the case of carbon fibers). Of course the axis of rotational symmetry coincides with the geometrical axis of symmetry of the bundle. It should be mentioned that the homogenization process provides the overall behavior of the bundles in their respective (local) coordinate system. Transformation from the local to the global (i.e. SMC or SMC-hybrid composite) coordinate system, according to the rules of tensor algebra, must take place before proceeding to the next homogenization, which gives the overall behavior of the whole composite.

The Composite Cylinders Method, originally developed by [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF], treats the problem of long fiber composites by considering an RVE of concentric cylinders, the interior being the fiber and the exterior the matrix phase. Since such a composite behaves as a transversely isotropic medium, five independent material properties are BVPs, the application of the Hill-Mandel energetic principle allows to compute accurately all the properties but the in plane shear modulus, for which the original CCM of [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] provides only bounds. To resolve this issue, [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF] have proposed for this property the use of the generalized self-consistent method.

In this technique, a third external layer is added to the RVE, which is the effective medium itself (Figure 4 e ). This permits to identify analytical solutions as a function of the unknown property, which is calculated with the help of energetic criteria. Computational details for the application of the composite cylinders method can be found in various papers [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Christensen | A critical evaluation for a class of micro-mechanics models[END_REF][START_REF] Seidel | Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites[END_REF][START_REF] Seidel | Micromechanics modeling of the multifunctional nature of carbon nanotube-polymer nanocomposites[END_REF][START_REF] Chatzigeorgiou | Effective mechanical properties of fuzzy fiber composites[END_REF].

B. Bundles containing microcracks

When dealing with polymeric composites, damage mechanisms very often occur at various scales. Damage can appear at the matrix phase, at the interfaces between matrix and fibers or matrix and bundles, or at the fibers themselves. As mentioned in the introduction, this work attempts to integrate damage that is observed at the level of the bundles by using a hybrid phenomenological-micromechanical model, developed by [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF]. It should be pointed out that the damage model in this work is a "brick" that can be easily interchanged with another damage model. Since the [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF] damage law is quite complicated and has already been published elsewhere, only the general description of the model is presented here. In fact, this constitutive law for fiber bundles has been experimentally validated in natural fibers multi-directional laminate composites.

bundle In unidirectional fiber composites damage commonly appears in the transverse direction of the fibers. Microcracks are initiated due to the debonding of the fiber-matrix interface and propagate by coalescence. The hybrid model of [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF] considers that initially the unidirectional composite is elastic, whose elasticity tensor C 0 is provided by any micromechanics approach (in this work the CCM is considered for the bundles). Then, as the loading increases, penny-shape micro-cracks start to be formed (Figure 5). The density of these cracks is denoted by γ c and they are considered as voids. Using the Mori-Tanaka homogenization scheme, the damaged fiber composite is substituted by a homogenized medium whose material properties depend on the density γ c (Figure 5). In this

+ MTM C 0 γ c micro-cracks C 0 -D(γ c ) damaged bundle

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

homogenized medium, the stress-strain constitutive law is written as

σ = C 0 -D(γ c ) :(ε -ε s ), ( 23
)
where C 0 is the initial elasticity tensor of the bundle (considered as elastic unidirectional composite), D is the fourth order reduction tensor due to the presence of the micro-cracks and ε s is the anelastic strain linked with the damage (for instance, due to micro-cracks non-closure and the resultant sliding with friction). The tensor D is a function of γ c and is computed analytically through the Mori-Tanaka scheme:

C 0 -D(γ c ) = (1 -γ c )C 0 : A 0 (γ c ). ( 24
)
In the above expression, the bundle and the voids-cracks are assumed to form a two-phase composite. The bundle's strain localization tensor A 0 and the void inclusions' strain localization tensor A c respect the consistency condition for two-phase media,

(1 -γ c ) A 0 + γ c A c = I . ( 25
)
The anelastic strain ε s follows an evolution law similar to plastic-like materials:

εs = Λ s (σ)γ c . (26)
The form of Λ s and the evolution criterion for γ c are extensively detailed in [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF]. This anisotropic damage constitutive law influences the stiffness only in the direction normal to the micro-cracks, namely the transverse and shear components. The material behavior in the direction of fibers remains elastic. To keep the manuscript in a concise form, the theoretical framework and the numerical implementation are not presented here. The reader is referred to [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF] for more details.

Second Homogenization

Next, the overall behavior of a SMC-Hybrid composite consisting of filled matrix (whose viscoelastic properties are taken from the 'zeroth' homogenization) and randomly oriented glass and carbon bundles (whose elastic properties have been computed in the first homogenization) is obtained, see Figure 3 c . Obviously, if carbon bundles are absent, the overall behavior reduces to that of a pure SMC composite. The homogenization process utilizes the Mori-Tanaka method taking into account the randomly oriented bundles, as it is described below. It should be noted that, considering the fabrication process, randomness exists in almost parallel planes. In this way, the final composite possesses a transversely isotropic elastic behavior, where the axis of rotational symmetry is normal to the aforementioned planes of randomness. In order to give more generality to the current work, a brief discussion about the Orientation Distribution Function (ODF) is included in this subsection.

As already pointed out in Section 3, working on the frequency domain implies the appearance of complex numbers in the computations. While in the case of 'zeroth' homogenization analytical expressions have been derived, in the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
case of the total SMC and SMC-hybrid composite the presence of short bundles at different orientations require cumbersome computations. Also, the presence of microcracks alter the linear representation of the global problem in the frequency domain. To overcome these difficulties, the time domain with incremental formalisms is instead adopted, which is easily applicable through available softwares for micromechanics methods.

An incremental Mori-Tanaka is used in order to obtain the overall viscoelastic response of the final composite.

The homogenized bundles are treated as short fibers, considering their low aspect ratio and they are assumed to possess ellipsoidal cross section. For ellipsoidal inclusions, analytical evaluations of the Eshelby tensor are available, in amenable form, for isotropic matrices [START_REF] Clyne | An introduction to metal matrix composites[END_REF][START_REF] Mura | Micromechanics of defects in solids[END_REF]. In the present paper the numerical algorithm of [START_REF] Gavazzi | On the numerical evaluation of eshelby's tensor and its application to elastoplastic fibrous composites[END_REF] is utilized (see also [START_REF] Desrumaux | Generalised mori-tanaka scheme to model anisotropic damage using numerical eshelby tensor[END_REF] for the calculation of the Eshelby tensor. This is a robust algorithm which has enjoyed a wide range of applications in, e.g., the prediction of the elastoplastic behavior of metal-matrix composites [START_REF] Lagoudas | Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme[END_REF] and the damaged behavior of random oriented fiber composites [START_REF] Meraghni | Implementation of a constitutive micromechanical model for damage analysis in glass mat reinforced composite structures[END_REF]. It has also been successfully implemented in numerical platforms for micromechanics tools and thermomechanical behavior studies, such as SMART+ (SMART+ development team, 2012).

Mori-Tanaka method (MTM), [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF], is a popular tool for the analysis of multi-phase materials. It has been used by a wide range of researchers to model the effective behavior of composites, and allows the average stress fields and overall effective stiffness of a composite with non-dilute concentration of inclusions to be determined. The method, its range of applicability and the connection of stiffness and compliant tensors estimations with classical variational bounds (e.g. Hashin-Shtrikman) are discussed in depth in [START_REF] Dvorak | Micromechanics of Composite Materials[END_REF].

Numerous applications have appeared in the literature with different ways of applying the Mori-Tanaka model, accounting for linear and nonlinear overall behavior. For example, [START_REF] Chen | Mori-tanaka estimates of the overall elastic moduli of certain composite materials[END_REF] derived explicit formulas for estimates of the effective moduli for many typical composite systems; [START_REF] Meraghni | Micromechanical modelling of matrix degradation in randomly oriented discontinuous-fibre composites[END_REF] applied a modified Mori-Tanaka model to study the effect of matrix degradation on the overall behavior of randomly oriented discontinuous-fiber composites; Tucker and Liang (1999) compared the predictions of MTM (as well as other micromechanical estimation schemes) with finite-element calculations for uni-directional short-fiber composites; [START_REF] Tan | The mori-tanaka method for composite materials with nonlinear interface debonding[END_REF] investigate the effect of nonlinear interface debonding on the macroscopic behavior of composite materials; [START_REF] Despringre | Analyse et modélisation des mécanismes d'endommagement et de déformation en fatigue multiaxiale de matériaux composites : polyamide renforcé par des fibres courtes[END_REF] has developed a multi-scale approach, using a modified incremental Mori-Tanaka method with 2N+2 phases, including coated reinforcements and the evolution of micro-scale damage processes in order to study the cyclic visco-damage behavior of short glass fiber reinforced polyamide (PA66/GF30). Many other references related to MTM can be found in the recent work of [START_REF] Liu | A Note on mori-tanaka's method[END_REF].

Without going into details, the Mori-Tanaka method, which takes into account inhomogeneity interactions, assumes that the average strain in the interacting inhomogeneities can be approximated by that of single inhomogeneity embedded in an infinite matrix subjected to the uniform average matrix strain; in other words, within a composite each homogeneity 'sees' a far-field strain equal to the average strain in the matrix.

In the case of M different types of unidirectional aligned inhomogeneities within the matrix, using the equivalent

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
inclusion method and the assumptions of Mori-Tanaka, it can be seen that the concentration tensor A r of the r-th inhomogeneity can be written as

A r = T r :                 1 - N q=1 c q         I + N q=1 c q T q         -1 , T r = I + S r : (L FM ) -1 : L r -L FM -1 , (27)
while for the filled matrix the concentration tensor is expressed as

A FM =                 1 - N q=1 c q         I + N q=1 c q T q         -1 . ( 28
)
In the above expressions L denotes the tangent modulus of a phase, i.e. the elastic modulus for the bundles and the algorithmic modulus for the viscoelastic filled matrix. Moreover, T r is called interaction tensor and relates the uniform strain in the r-th ellipsoidal inhomogeneity, if it was solely embedded in the filled matrix, to the average filled matrix strain. Also S r is the Eshelby tensor of the r-th inhomogeneity. Combining ( 27) and ( 19) yields

L SMC = L FM + M r=1 c r L r -L FM :T r :                 1 - N q=1 c q         I + N q=1 c q T q         -1 . ( 29
)
Due to the nonlinear nature of the problem, incremental methodology is required in which at every step the average stresses and strains of each phase are corrected until convergence is achieved. Such numerical schemes, in various forms, are available in the literature. In Appendix B the essential steps of an incremental Mori-Tanaka method are presented.

In composites such as SMC, preferential fiber orientation may be induced by the thermocompression process.

The detection of the actual distributions is not an easy task; see, e.g., [START_REF] Le | X-ray phase contrast microtomography for the analysis of the fibrous microstructure of smc composites[END_REF] for SMC and Arif (2014); Arif et al. (2014a,b) for injected reinforced polyamide. When they are available, they can be related to overall material symmetry of the composite by certain orientation distribution functions (ODFs).

Let g(θ, φ, ψ) be such an ODF with respect to the Euler angles θ, φ, ψ which connect the local (i.e., inhomogeneity) coordinate system, (x 1 , x 2 , x 3 ), and the global (i.e., composite) coordinate system, (X 1 , X 2 , X 3 ). A fourth order tensor C i jkl in the (X 1 , X 2 , X 3 ) coordinate system is related to the tensor C pqrs in the local inhomogeneity coordinate system (x 1 , x 2 , x 3 ) through the transformation rule of fourth-order tensors as

C i jkl = a ip a jq a kr a ls C pqrs , ( 30
)
where a i j are the directional cosines of angles relating the global to local coordinate systems. In terms of the Euler angles θ, φ, ψ the transformation matrix [a i j ] is given in [START_REF] Qu | Fundamentals of Micromechanics of Solidss[END_REF]. To incorporate the contribution from all possible orientations, C i jkl is multiplied by the ODF and integrated over all orientations. Lengthy computations can be avoided by utilizing the Voigt/Nye notation. Matrices in this notation are transformed using the (6 × 6) matrix X generated from the products of the directional cosines in agreement with the rules of transformations of first and second order tensors and the summation and contraction rules. This matrix can be found in [START_REF] Dvorak | Micromechanics of Composite Materials[END_REF].

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

For a multi-phase and matrix-based composite, reinforced by randomly oriented but otherwise identical ellipsoidal inhomogeneities, one may extend (29) as follows [START_REF] Entchev | Modeling porous shape memory alloys using micromechanical averaging techniques[END_REF][START_REF] Seidel | Micromechanics modeling of the multifunctional nature of carbon nanotube-polymer nanocomposites[END_REF][START_REF] Dvorak | Micromechanics of Composite Materials[END_REF])

L SMC = L FM + M r=1 c r L r -L FM :T r :                 1 - N q=1 c q         I + N q=1 c r T q         -1 , ( 31
)
where the symbol {•} denotes averaging over all possible orientations. The L FM in the above expression is provided in a global coordinate system, while c r denotes the total volume fraction of the r-th inhomogeneity. Among others, [START_REF] Norris | An examination of the mori-tanaka effective medium approximation for multiphase composites[END_REF]; [START_REF] Ferrari | On the domain of applicability of the mori-tanaka effective medium theory[END_REF]; [START_REF] Benveniste | On diagonal and elastic symmetry of the approximate effective stiffness tensor of heterogeneous media[END_REF]; [START_REF] Christensen | On the range of validity of the mori-tanaka method[END_REF]; Schjødt-Thomsen and Pyrz (2001); [START_REF] Dvorak | Micromechanics of Composite Materials[END_REF], discuss the restrictions on constituent shape alignment in the application of equation ( 31). As discussed in these references, for general ellipsoidal fiber shape the obtained macroscopic tensor L SMC may not be symmetric and proper regularization is required for symmetrize it. In this work, a regularization of the final stiffness response is considered whenever the symmetry of the macroscopic tensor is not satisfied.

In numerical methods, integrals that appear in the averaging procedure are written in a discretized form by considering a large but finite number of possible orientations. In the SMC composite, each bundle orientation is considered as a new inhomogeneity with its own volume fraction. An ODF designed for SMC composites should contain information about the in-plane angle (i.e., the plane in which the bundles are distributed). The tangent modulus L and the interaction tensor T of each bundle orientation are transformed to the global coordinate system. Then, one can apply directly the expression (29) for the extended number of "inhomogeneities".

Numerical examples

This Section presents various examples and parametric studies for SMC and SMC-hybrid composites. The elastic properties of the glass particles/fibers and the carbon fibers are provided in Table 1.

glass particles/fibers

property value

E G (MPa) 81000 ν G 0.25 carbon fibers
property value With regard to the matrix phase, the CaCO 3 -filled polyester of [START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF]. In this work the authors have performed creep tests and they have obtained an appropriate creep compliance. At reference temperature 366.3 K, the the same material, computed with the procedure described in the Appendix I of [START_REF] Park | Interconversion between Relaxation Modulus and Creep Compliance for Viscoelastic Solids[END_REF]. The figures also illustrate the numerical simulation using Prony series (material parameters are provided in Table 2, columns 3 and 4).

E C A (MPa) 241000 E C T (MPa) 14500 G C A (MPa) 22800 G C T (MPa) 4800 ν C AT 0.27
M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 1.6E-04 1.8E-04 2.0E-04 2.2E-04 2.4E-04 2.6E-04 2.8E-04 3.0E-04 -3 -2 -1 0 1 2 3 4 5 Creep compliance (1/
master curve for the creep compliance D m (t) is given by the expression

D m (t) = 1.62•10 -4 + 1.82•10 -5 t 0.18 MPa -1 ,
where the time is expressed in seconds. The Poisson's ratio is reported to be ν m = 0.35. In the sequel all the numerical examples are considered to be performed under isothermal conditions at temperature 366.3 K.

The Prony series formalism of equation ( 9) is usually calibrated with the help of relaxation tests. Assuming constant Poisson's ratio, unidirectional relaxation tests provide the relaxation modulus E m (t), from which the bulk modulus K m (t) and the shear modulus G m (t) are computed with the usual formulas

K m = E m 3(1 -2ν m ) , G m = E m 2(1 + ν m ) . (32)
Since in this case relaxation tests are not available, one needs to utilize the relation between creep compliance and relaxation modulus [START_REF] Christensen | Theory of Viscoelasticity, 2nd Edition[END_REF][START_REF] Park | Interconversion between Relaxation Modulus and Creep Compliance for Viscoelastic Solids[END_REF])

t 0 E m (t -τ)D m (τ)dτ = t, (33)
in order to identify E m (t). The convolution integral of equation ( 33) can be solved numerically in terms of the unknown relaxation modulus with the help of the numerical procedure presented in the Appendix I of [START_REF] Park | Interconversion between Relaxation Modulus and Creep Compliance for Viscoelastic Solids[END_REF].

Both the creep compliance and the relaxation modulus of the matrix phase are presented in Figure 6.

After obtaining the relaxation modulus, its Prony series representation is written in the temporal space as

E m (t) = E m ∞ + N j=1 E m j e -t/τ m j .
One can choose a sufficient discrete number N of relaxation times, spanning the time interval of the experimental curve, and then utilize the least squares method for identifying the E ∞ and the E j for j between 1 and N. The Prony

m j (s) K m j (MPa) G m j (MPa) K m j (MPa) G m j (MPa) ∞ -

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

series coefficients and relaxation times for the matrix bulk and shear moduli are shown in Table 2. The simulation of the relaxation modulus and the creep compliance using Prony series is presented in Figure 6, where it is observed an excellent match between experimental and numerical results.

'Zeroth' Homogenization

The polyester of [START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF] is already reinforced with CaCO 3 fillers, so the 'zeroth' homogenization could be omitted. Still, for demonstrating the full extent of the proposed multiscale methodology, in the following examples the matrix in the SMC composite or SMC-hybrid composite is considered to be also filled with 5% glass particles, whose properties are provided in Table 1. Figure 7 illustrates the storage and loss moduli for the bulk and shear components of the filled matrix stiffness tensor as a function of the frequency and the comparison with the corresponding moduli of the "pure" matrix. The filled matrix results have been obtained with the use of formulas (21).

In order to pass to higher scales, it is essential to identify an "analytical-type" constitutive law for the filled matrix.

This can be achieved by considering that in the frequency domain the bulk and shear moduli (both storage and loss)

for the filled matrix can be written in Prony series form as

K FM (ω) = K FM ∞ + N j=1 K FM j ω 2 1 τ 2 j + ω 2 , K FM (ω) = N j=1 K FM j ω τ j 1 τ 2 j + ω 2 , G FM (ω) = G FM ∞ + N j=1 G FM j ω 2 1 τ 2 j + ω 2 , G FM (ω) = N j=1 G FM j ω τ j 1 τ 2 j + ω 2 .
To simplify the identification procedure, the relaxation times are considered to be the same with those of the matrix phase. Then, the least squares method can be utilized in the frequency domain, where experimental data are considered the results obtained from the 'zeroth' homogenization. Both the storage and the loss moduli are utilized for identifying the K FM j coefficients for the bulk and the G FM j coefficients for the shear modulus. The obtained coefficients are shown in Table 2. Moreover, the Prony series results in the frequency domain are compared with the homogenization results of Figure 7, demonstrating that this analytical representation captures accurately the filled matrix response.

First Homogenization

As already discussed in Section 4, the glass and carbon bundles are considered to behave elastically. The "elastic" properties of the filled matrix are assumed to be those of its instantaneous response. Using Table 2 the instantaneous bulk and shear moduli are

K FM 0 = K FM ∞ + 16 j=1 K FM j = 7101.276 MPa, G FM 0 = G FM ∞ + 16 j=1 G FM j = 2445.098 MPa.
Considering fiber content of 80% volume fraction inside a glass or carbon bundle, the homogenized properties for the bundles obtained by the Composite Cylinders Method are summarized in Table 3. (1982) and the filled matrix results obtained by the 'zeroth' homogenization. The figures also illustrate the numerical simulation using Prony series for the field matrix (material parameters are provided in Table 2, columns 5 and 6).

glass bundles

property value

E Gb A (MPa) 66130 E Gb T (MPa) 35160 G Gb A (MPa) 13215 G Gb T (

Second Homogenization

In the second homogenization step the glass and/or carbon elastic bundles are embedded in a planar randomly oriented way (with a uniform ODF) inside the viscoelastic matrix to produce the final SMC or SMC-hybrid composite.

As already explained in the Section 3, transforming the problem in the Fourier space provides a linear form of all the implicated constitutive laws, but the calculations involve complex numbers that cannot be easily adopted in existing micromechanics schemes. Moreover, addition of other nonlinear mechanisms, like for instance damage at the interface between the bundles and the matrix, introduce material behaviors without analytical representation in the Fourier space. These limitations do not exist in the differential representation of the viscoelastic constitutive law, described in Section 3. Thus, for this homogenization step the incremental Mori-Tanaka presented in Appendix B is utilized, where the stress and the tangent modulus of the viscoelastic matrix are computed with the numerical procedure of Appendix A. In all the numerical examples of this subsection, the composite is subjected to a rapid uniaxial in-plane stretching (Figure 8), in which at 0.01 s the strain reaches the value of 0.001. Then the strain is kept constant for 1000 s, causing a stress relaxation response. The first parametric analysis aims at identifying the impact of the combination between glass and carbon bundles in the SMC-hybrid response. In these numerical examples, the total fiber reinforcement volume fraction is kept constant to 40%. Since the volume fraction of the fibers inside a carbon or glass bundle is 80%, the total bundle volume fraction inside the filled matrix is equal to 50% (see equation ( 5)). Varying the total glass fiber content from 0 to 40%, the total carbon fiber content is reduced analogously. The bundles are considered cylindrical and long (length/diameter equal to 1000). The stress relaxation, as a function of time, for an SMC-hybrid composite with 20% glass fibers and 20% carbon fibers is illustrated in Figure 9 and is compared with the corresponding response of the filled matrix. As it is observed, the presence of glass and carbon bundles enhance significantly the response of the composite, compared to the filled matrix, but on the same time it reduces drastically its viscous behavior: the filled matrix has a stress reduction of 20.5% between 1 and 1000 s, while the analogous reduction for the SMC-hybrid composite is only 2.3%. Figure 10 presents the obtained stress as a function of the total glass fiber volume fraction at time t =1000 s. As expected, the increase of carbon glass content increases the overall strength of the composite. From the results it is shown that, varying the ratio between carbon and glass bundles, the stress response can be between 45 and 22 MPa at 1000 s. The second parametric analysis investigates the role of the bundle length in the SMC-hybrid behavior. In this study, the SMC-hybrid composite has 50% volume fraction of filled matrix and equal volume fractions of glass and carbon bundles, which are assumed cylindrical. The bundles length/diameter ratio varies between 10 and 1000. Figure The third parametric analysis examines the influence of the bundles shape. For the needs of the micromechanics analysis, the bundles are assumed as general ellipsoids with axes a 1 , a 2 and a 3 (Figure 12). To simplify the discussion, it is considered that the principal fiber axis is the a 1 , the axis normal to a 1 and in the fibers plane is the a 2 , and the

bundles of carbon fibers bundles of glass fibers

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

axis normal to the fibers plane is the a 3 . The SMC-hybrid composite has 50% volume fraction of filled matrix and equal volume fractions of glass and carbon bundles, whose aspect ratio a 1 /a 3 is 100. The ratio a 2 /a 3 varies between 1 (cylindrical fiber) and 100. Figure 13 presents the results of this parametric analysis. As it is observed, the "flatness" of the bundle reduces the stress but not significantly, since the drop of stress in the extreme case of a 2 /a 3 ratio does not exceed 10% compared to the cylindrical bundles. Figure 13: (a) Evolution of stress with time for two bundle shapes, (b) stress as a function of the bundles ellipsoidal shape (a 2 /a 3 ratio) at time t =1000 s. The SMC-hybrid composite has 20% glass and 20% carbon fibers and the bundles have a 1 /a 3 ratio equal to 100.

Second Homogenization including damage in glass bundles

To illustrate the capability of the present multiscale method to incorporate additional inelastic mechanisms, the following examples consider that the glass bundles are experiencing damage that is taken into account with the hybrid model described in subsection 4.2.

The first example is based on real experimental data for conventional SMC glass fiber reinforced composites. The material has been provided by Owens Corning. The matrix phase is polyester, reinforced with calcium carbonate (CaCO 3 ) spherical fillers. The total reinforcement of the glass fibers in the SMC composite is 45%. The glass bundles are considered cylindrical, whose dimensions are experimentally identified as: 25 mm length and 0.189 mm diameter.

Thus, the length per diameter aspect ratio is 132.28.

Due to lack of additional information, some hypotheses about material parameters and volume fractions are necessary. With regard to the filled matrix, no experimental data were available in terms of the filler content and the viscous characteristics of the specific polyester. The analysis performed here assumes that the filled matrix has the same characteristics with the material of [START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF], whose parameters are provided in Table 2, columns 1, 2, 3 and 4. Moreover, the volume fraction of fibers inside bundles is fixed to 80%. Using equation ( 5), one obtains that the bundles inside the SMC occupy a volume fraction of 56.25%. For this conventional SMC composite, a uniform random distribution of the bundles is assumed.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

With these data, the numerical simulation can be performed as explained in the previous sections. Since the matrix of [START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF] is already filler-reinforced, the 'zeroth' homogenization step is unnecessary. The first homogenization step with the glass fibers and the matrix material give the following values for the glass bundles:

E Gb A =66.014 GPa, E Gb T = 33.591 GPa, G Gb A =12.48 GPa, G Gb T =12
.57 GPa and ν Gb AT =0.266. With regard to the damage response, the damage-related properties for the glass bundles have been identified using inverse engineering methods from the total composite response, and are included in Table 4. Explanations about the parameter identification and the physical meaning of the material constants are provided in [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF]. configuration of the void inclusions flattening ratios shape ratio 4: Damage related parameters for the glass bundles. The physical meaning of each parameter and their effect on the stiffness tensor have been presented in detail in [START_REF] Praud | Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by micro-cracks in unidirectional composites[END_REF].

δ 1 = a 1 a 2 δ 3 = a 3 a 2 δ c = δ 1 δ 3 = a 1
Three specimens have been subjected to uniaxial tensile loading (Figure 8). The experimental tests have been performed under quasi-static conditions (strain rate of 10 -3 s -1 ), which in the numerical analysis are represented with a rate of 1% strain per 10 s. The specimens response and the model simulation are shown in Figure 14. As it can be seen, the model is capable of capturing quite accurately the nonlinear overall response. A small deviation on the elastic regime is due to the high modulus of the matrix viscoelastic constitutive law given by [START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF]. A thorough experimental investigation of the pure resin could provide a better representation of the filled matrix, which will permit to capture the initial part of the curve. Such task though is out of the scope of the present manuscript.

After demonstrating the model's ability to simulate results from real experimental data, a parametric analysis can show its capacity in complicate studies of SMC composites. The hybrid composite of the second example contains 30% filled matrix, 60% glass bundles and 10% carbon bundles. Both glass and carbon bundles are long, cylindrical and randomly oriented on a plane. Since the fiber content inside the bundle is equal to 80%, the glass fibers and the carbon fibers in the overall SMC composite have volume fractions 48% and 8% respectively. The viscoelastic properties of the filled matrix are those summarized in Table 2, columns 1, 2, 5 and 6, while the elastic properties of the glass and carbon bundles are those given in Table 3. As in the previous case, the damage-related properties for the glass bundles are those of Table 4. The bundles are considered long, cylindrical with length per diameter ratio equal to 10000. In this numerical example, the composite is subjected to a uniaxial in-plane stretching (Figure 8), in which at 1000 s the strain reaches the value of 0.0032. Figure 15 demonstrates the obtained stress-strain response in the direction of stretching considering three different overall time steps for the analysis: 500 steps (time increment 2 s), 800 steps (time increment 1.25 s) and 1250 steps (time increment 0.8 s). The curves illustrate the convergence rate of the proposed method. It should be mentioned that the most demanding analysis of 1250 steps, including saving the microscopic and macroscopic information, required a computational time (actual time) of approximately 7.5 minutes in a personal Figure 16 summarizes the microscopic stress-strain response at three different material phases: the filled matrix, the glass bundles parallel to the applied stretching (0 • ) and the glass bundles normal to the applied stretching (90 • ).

The local strain and stress components illustrated in the Figure 16 are those corresponding at the global direction of the loading. The 90 • glass bundles experience significant damage, while the 0 • glass bundles remain elastic. It is worth mentioning that the current model does not take into account activation of additional damage mechanisms in the filled matrix or in the carbon bundles. Another remark is that the damage mechanism at the glass bundles is the one that mainly causes the nonlinearity in the macroscopic stress-strain curve (Figure 15). The filled matrix viscous behavior presents very reduced nonlinear response (Figure 16 b ).

Concluding Remarks

The present paper has studied the response of SMC and SMC-hybrid composites from a micromechanical point of view, by considering a hierarchical multiscale approach, in which the finer scale provides information to the higher

Figure 1 :

 1 Figure 1: Schematic representation of (a) an SMC-Hybrid composite; top view, (b) a glass or carbon fiber-bundle system and (c) the filled (i.e. particle reinforced) matrix.

  "Zeroth" homogenization: Effective properties of filled matrix. Application of Mori-First homogenization: Effective properties of a fiber-bundle system. Application of Composite Cylinders Method. fibers bundles of glass fibers (c) Second homogenization: Effective properties of SMC composite. Application of Mori-Tanaka Method.

Figure 3 :

 3 Figure 3: Homogenization of an SMC/SMC-Hybrid composite: Multiscale hierarchical modeling.

  Figure 4: Boundary value problems for the evaluation of the five overall elastic moduli, for a transversely isotropic fibrous composite.

Figure 5 :

 5 Figure5: The bundle initially is transversely isotropic elastic material. At damaged state micro-cracks appear with density γ c (void volume fraction). The stiffness reduction tensor D of the damaged bundle is computed usingMori-Tanaka.

Figure 6 :

 6 Figure 6: Evolution with time of (a) the creep compliance for the CaCO 3 -filled polyester of Jerina et al. (1982), (b) the relaxation modulus of

  time bulk modulus shear modulus bulk modulus shear modulus number j τ

Figure 7 :

 7 Figure 7: Bulk and shear moduli (both their storage and loss parts) as a function of log ω, ω in 1/s. Comparison between the matrix of Jerina et al.

Figure 8 :

 8 Figure 8: Uniaxial stretching (displacement controlled conditions) is SMC-hybrid composite.

Figure 9 :

 9 Figure 9: Stress relaxation responses of (a) the filled matrix and (b) a SMC-hybrid composite with 20% glass and 20% carbon cylindrical fibers.

Figure 10 :

 10 Figure 10: Stress versus total glass fiber volume fraction at time t =1000 s for a SMC-hybrid composite.

  11 a illustrates the evolution of the stress with time for length per diameter ratios 20 and 200, while Figure 11 b shows the stress as a function of the length per diameter ratio at time t =1000 s. As the latter Figure depicts, the fiber length/diameter ratio influences the response of the composite as well as it is below 200. Above 200 the bundles can be considered as long.

Figure 12 :

 12 Figure 11: (a) Evolution of stress with time for two length/diameter ratios. (b) Stress versus bundles length/diameter ratio at time t =1000 s. In both figures the SMC-hybrid composite contains 20% glass and 20% carbon cylindrical fibers.

Figure 14 :

 14 Figure14: Macroscopic stress-strain response of SMC glass fiber reinforced composite: comparison between experimental results on three specimens and numerical simulations with the proposed model.

Figure 15 :

 15 Figure15: Macroscopic stress-strain response of SMC-hybrid including damage in the glass bundles. The analysis was performed at three different overall time steps: 500 steps (time increment 2 s), 800 steps (time increment 1.25 s) and 1250 steps (time increment 0.8 s).

  Core i7, 2.7 GHz CPU and 8GB RAM.

  Figure 16: (a) Schematic of SMC-hybrid, showing bundles at 0 • and at 90 • with respect to the loading direction. (b)-(d) Microscopic stress-strain response of filled matrix and glass bundles at 0 • and 90 • .

Table 1 :

 1 Elastic properties of glass fibers/particles and carbon fibers. The carbon fibers are transversely isotropic.

Table 2 :

 2 Relaxation times and prony series coefficients and for a) CaCO 3 -filled polyester of[START_REF] Jerina | Viscoelastic Characterization of a Random Fiber Composite Material Employing Micromechanics[END_REF] (3rd and 4th column) and b) glass particles filled matrix (5th and 6th column). The filled matrix coefficients have been obtained from the 'zeroth' homogenization step.

			3673.976	1224.659	3925.393	1352.835
	1	1•10 -3	34.038	11.346	35.733	12.280
	2	3.1622•10 -3	46.845	15.615	49.192	16.906
	3	1•10 -2	56.010	18.670	58.821	20.215
	4	3.1622•10 -2	68.769	22.923	72.242	24.829
	5	1•10 -1	81.912	27.304	86.059	29.578
	6	3.1622•10 -1	99.461	33.154	104.544	35.933
	7	1	116.974	38.991	122.964	42.265
	8	3.1622	140.768	46.923	148.091	50.908
	9	1•10 1	162.140	54.047	170.567	58.633
	10	3.1622•10 1	194.473	64.824	204.868	70.439
	11	1•10 2	212.698	70.899	223.945	76.990
	12	3.1622•10 2	259.588	86.529	274.032	94.249
	13	1•10 3	264.648	88.216	278.840	95.870
	14	3.1622•10 3	336.314	112.105	356.246	122.591
	15	1•10 4	226.383	75.461	237.969	81.785
	16	3.1622•10 4	708.031	236.010	751.770	258.793

Table 3 :

 3 Elastic properties of glass and carbon bundles, obtained by the first homogenization.

			carbon bundles
			property	value
			E Cb A (MPa) 194120
			E Cb T (MPa) 12480
			G Cb A (MPa) 11330
		MPa) 13210	G Cb T (MPa) 4140
	ν Gb AT	0.266	ν Cb AT	0.284

Several studies(Fisher and Brinson,

2001;[START_REF] Krairi | A thermodynamically-based constitutive model, for thermoplastic polymers coupling viscoelasticity viscoplasticity and ductile damage[END_REF] allow different relaxation times and number of Prony series components in bulk and shear moduli. Without loss of generality and for simplifying the subsequent discussion, this work assumes that the whole viscoelastic modulus tensor L is described by a fixed number of Prony series components and one set of relaxation times.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

scale. The special geometrical aspects (spherical fillers, bundles consisting of long fibers, short bundles at the larger scale) and the viscoelastic nature of the matrix have enforced different homogenization techniques and different ways of accounting for the nonlinear viscoelastic response at each scale of the composite. The proposed methodology attempts to use the most efficient homogenization method according to the needs of every scale and tends to provide computational efficiency. The analytical results at the lower scales (matrix with fillers, bundles consisting of matrix and long fibers) permit to adopt a nonlinear incremental methodology at the larger scale.

The numerical examples presented in Section 5 illustrate the capability of the proposed homogenization strategy to provide complicated parametric studies, implicating various aspects such as phases volume fractions as well as form and length of the bundles, with low cost in calculations. In addition, the model takes into account damage mechanism at the scale of the bundle by integrating a proper constitutive law in the incremental Mori-Tanaka scheme.

Simulation of real experimental data in SMC glass fiber reinforced composites demonstrates that the model can be adopted easily in industrial applications. While it is true that the polyester-based composites do not show significant viscoelastic response, the proposed hierarchical modeling strategy can well capture this aspect and can be utilized in the same form for other polymer-based reinforced short bundles composite systems that exhibit a viscous behavior more prominent.

Appendix A. Numerical scheme for a linear viscoelastic material

The convex cutting plane (CCP) approach is a special form of the more general class of return mapping algorithms [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Qidwai | Numerical Implementation of a Shape Memory Alloy Thermomechanical Constitutive Model Using Return Mapping Algorithms[END_REF] and it is an implicit scheme in terms of time integration. Two important operators are introduced in this algorithm:

1. The time increment ∆, which denotes the difference between the values at two consequent time steps t n and t n+1 of a quantity x, i.e. ∆x = x(t n+1 )x(t n ).

2. The iteration increment δ, which denotes the difference between the values at two consequent iteration steps k and k + 1 of a quantity x at the time step t n+1 , i.e. δx = x k+1 (t n+1 )x k (t n+1 ).

In the iterative calculation procedure of the CCP algorithm the evolution equation ( 18) 1 is integrated explicitly, thus it considers that the viscoelastic flow directions Λ j keep the value of the previous iteration step. This technique allows to avoid complicated expressions that consider derivatives of Λ j , which introduce fourth order tensors.

The linearized forms of ( 16) and ( 18) are written as

where ∆t = t n+1t n . In the above expressions, any net quantity without an operator ∆ or δ is either constant or refers to the iteration step k. The CCP is a strain driven approach, in the sense that one provides the total strain ε at time t n and the viscoelastic strains ε v i at time t n+1 and obtains the stress σ, the viscoelastic strains and the viscoelastic tangent modulus L t at time t n . The algorithm is split in three parts:

1. Elastic prediction. At the beginning of the calculations it is assumed that no viscous effect appears. In that case all δs j are zero, the viscoelastic strains ε v j do not evolve and the stress is computed by the expression σ = L 0 :ε.

2. Viscoelastic correction. Since the viscoelastic mechanism evolves with time at every instance of the loading, the elastic prediction causes an error in the viscoelastic criteria Φ j . This error is corrected iteratively with the following procedure:

During the viscoelastic correction the term δε is considered zero, since the total strain has already been taken into account in the elastic prediction. Moreover, the viscoelastic criteria Φ j are corrected through an iterative Newton-Raphson procedure

The last expression, combined with (A.1) 3 , permits to identify the viscoelastic scalars at each Maxwell branch and at each iteration step k by the relation

Next, the stress and the viscoelastic strains are updated through the equations (A.1) 1,2 and the viscoelastic criteria Φ j are evaluated from (18) 3 . The procedure continues until all δs j are smaller than a tolerance. Usually the iterative approaches in viscoelastic materials converge very fast and most of the times only one iteration step is required.

3. Viscoelastic tangent modulus: In finite element analyses involving nonlinear materials the knowledge of the correct stress is not sufficient to proceed to the next time increment. Indeed, in implicit numerical schemes a tangent operator relating the increments of stress and total strain is required. The CCP algorithm provides a continuum tangent modulus. Once the previous two parts are finished, δε is "released" from being zero and L t is computed by utilizing equations (A.1) and assuming that δΦ j are zero. This leads to the expression

An important remark is that this tangent modulus does not preserve the major symmetries. Since the incremental Mori-Tanaka scheme (especially the calculation of the Eshelby tensor) requires symmetric tangent moduli tensors, the symmetric part of L t L t sym

is utilized here.

Appendix B. Incremental Mori-Tanaka

When a composite has nonlinear material constituents, a micromechanics framework like the Mori-Tanaka method needs to be extended in order to account for an incremental formulation that allows the computation at each step of appropriate tangent moduli.

Generally speaking, the homogenization can be seen as a two scale problem, the macroscale and the RVE (Figure macroscopic strain ε macro . The microscale (RVE) problem utilizes this information and the appropriate concentration tensors to obtain the macroscopic stress σ macro through an iterative process, as well as the necessary tangent modulus L SMC of the SMC (or SMC-hybrid) composite.

Microscale problem

Once the macroscopic strain is provided by the macroscale analysis, the microscale problem is split in the following parts:

1. For an elastic prediction, the time increment of the microscale strains ∆ε r at each phase of the composite are calculated by the linearized relations ∆ε r = A r :∆ε macro , where the concentration tensors A r are provided by equation ( 27) (for the M inhomogeneities) and (28) (for the filled matrix). The microscopic tangent moduli L r at the first iteration are considered as if all the phases are in elastic state (elastic prediction). During the iteration process their value change according to the evolution of the nonlinear mechanisms (i.e. viscoelasticity for the filled matrix).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

2. Using the microscopic strain ε r , the constitutive law of each material phase is utilized in order to calculate the microscopic stresses σ r . For nonlinear materials the stress is usually computed iteratively (see for instance Appendix A for a linear viscoelastic material).

3. Once the second step is completed, the updated tangent moduli L r , as well as the updated concentration tensors A r , for each phase are utilized in order to obtain the tangent modulus L SMC of the the SMC (or SMC-hybrid) composite from the relation (29).

Macroscale problem

Ignoring inertia and body forces, the macroscopic problem can be written in linearized incremental form as div σ macro + δσ macro = 0, (B.1)

where