Closed-form expressions for effective viscoelastic properties of fiber-reinforced composites considering fractional matrix behavior

Mechanics of Materials 127 (2018) 14–25

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Research paper

Closed-form expressions for effective viscoelastic properties of fiber-reinforced composites considering fractional matrix behavior

T

⁎

U. Hofer , M. Luger, R. Traxl, R. Lackner Material Technology Innsbruck (MTI), University of Innsbruck, Technikerstraße 13, Innsbruck A-6020, Austria

A R T I C LE I N FO

A B S T R A C T

Keywords: Fractional viscoelasticity Fiber-reinforced-polymers Effective properties

For determination of effective elastic properties of fiber-reinforced materials the so-called Chamis-equations are standardly used (Chamis, 1989). In this paper, these equations are extended towards viscoelasticity, making use of the correspondence principle (Lee, 1955). The determination of the effective properties takes place in the Laplace–Carson domain, with the back-transformation to the time domain being performed analytically, providing closed-form expressions for the effective viscoelastic behavior of the composite. Fractional viscoelastic models are used to describe the time-dependent behavior of the matrix, well-capturing the experimentallyobserved viscoelastic behavior with a comparatively low number of parameters. The performance of the proposed determination of effective viscoelastic properties is assessed by means of alternative methods such as Eigenstrain-based methods and numerical simulations employing unit-cell (UC) models.

1. Introduction Fractional viscoelastic models (see Fig. 1) are standardly employed to describe creep and relaxation of a variety of materials (Podlubny, 1999; Mainardi, 2010): Yancey and Pindera (1990) used a fractional Maxwell model to represent creep data of epoxy resin; in (Dinzart and Lipiński, 2010) a fractional Zener model is used to model creep of PMMA (polymethyl methacrylate) and SAN (styrene acrylonytrile copolymer); Pichler et al. (2012) identified a fractional Maxwell model for the viscoelastic behavior of bitumen, with Oeser et al. (2008) and Celauro et al. (2012) using fractional Burgers models to fit creep data of asphalt mixtures; moreover, soft biological tissues are found to obey a fractional Kelvin-Voigt model as proposed in (Meral et al., 2010). Representing the power-law model, i.e. the fractional Maxwell model, by a mechanical spring-dashpot combination, Deseri et al. (2014) found a unique expression for the corresponding free energy, making the fractional Maxwell model even more attractive for being used in viscoelastic material modeling. For materials consisting of more than one material phase, its effective (homogenized) elastic as well as viscoelastic behavior may be determined from the behavior of its constituents by means of homogenization methods: Hashin and Rosen (1964) derived expressions for the effective moduli of fiber-reinforced materials by using the composite cylinder assemblage model, consisting of cylindrical fibers embedded in a matrix. The elastic strain energy is formulated for a

⁎

kinematically admissible strain field and a statically admissible stress field, leading to lower and upper bounds for the strain energy. Expressions for the effective elastic moduli are finally found by comparing the strain energy of the cylinder assemblage model with the strain energy of the effective material. The same author (Hashin, 1966) extended the aforementioned method to viscoelastic behavior by applying the correspondence principle (Lee, 1955), where the viscoelastic matrix material behavior was represented by a Maxwell model. The backtransformation to time domain was performed analytically assuming rigid fiber behavior. In (Upadhyaya and Upadhyay, 1991), Hashin’s composite cylinder model was employed to obtain the effective viscoelastic properties of unidirectionally-reinforced composites, with the viscoelastic material behavior of the matrix extended towards a Prony series, i.e. a generalized Maxwell model. Using the equivalent Eigenstrain method (Nemat-Nasser and Hori, 1993), the effective elastic properties of a medium with inclusions are determined. An Eigenstrain field is applied onto an effective medium in such a way, that it represents the strain field of the respective heterogenous material. By comparing the stress and strain fields, respectively, in the heterogenous and homogenous medium, the elastic stiffness tensor is obtained. For ellipsoidal inclusions, the solution given in (Eshelby, 1957) leads to explicit expressions for the elastic properties. For cylindrical inclusions instead, i.e. for fiber-reinforced composites, explicit expressions for the elastic properties were presented in (Nemat-Nasser and Hori, 1993), which have been adapted

Corresponding author. E-mail address: [email protected] (U. Hofer).

https://doi.org/10.1016/j.mechmat.2018.08.005 Received 5 February 2018; Received in revised form 24 July 2018; Accepted 16 August 2018 Available online 18 August 2018 0167-6636/ © 2018 Elsevier Ltd. All rights reserved.

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

(Luciano and Barbero, 1995; Barbero and Luciano, 1995) in order to obtain analytical expressions for the effective relaxation moduli of unidirectionally-reinforced composites considering a Burgers model for the viscoelastic behavior of the matrix material. By application of the correspondence principle, the viscoelastic problem is solved in the Laplace–Carson domain analogously to the elastic problem. In an adaption of the equivalent Eigenstrain method, i.e. the selfconsistent scheme, the effective medium is used for the medium surrounding the inclusions instead of the pure matrix (Budiansky, 1965; Hill, 1965). This scheme was extended towards viscoelasticity in (Laws and McLaughlin, 1978) for spheroidal inclusions and cylindrical fibers by applying the correspondence principle. Finally, Dinzart and Lipiński (2010) extended the self-consistent method towards fractional viscoelastic behavior, representing the polymeric constituents by a fractional Zener model. Experimental validation was exclusively performed in the Laplace–Carson domain employing cyclic-test results; hence, back-transformation of the effective relaxation moduli to the time domain became obsolete. Wang and Pindera (2016) extended an elastic homogenization theory introduced by Drago and Pindera (2008) towards viscoelasticity by using the correspondence principle. The viscoelastic matrix behavior was represented by a fractional Maxwell model, the effective viscoelastic properties were back-transformed to the time domain numerically. In this paper, a closed-form solution for the effective viscoelastic behavior of fiber-reinforced composites considering fractional viscoelastic matrix behavior is presented. For this purpose, the Chamisequations – standardly employed for determination of the effective elastic properties – are extended towards viscoelasticity. The homogenization is performed in the Laplace–Carson domain. By means of back-transformation to the time domain, closed-form expressions for the effective viscoelastic moduli are found. This paper is organized as follows: In Section 2, fractional viscoelasticity is briefly reviewed, dealing with the fractional Burgers model, covering the fractional Zener, Maxwell, and Kelvin-Voigt models as special cases. The extension of the Chamis-equations towards viscoelasticity is treated in Section 3, finally providing closed-form expressions for the effective viscoelastic properties. In Section 4, the performance of the extended Chamis-equations is assessed by alternative homogenization methods for determining the effective viscoelastic behavior, such as the equivalent Eigenstrain method and numerical UCbased simulations. The paper closes with concluding remarks in Section 5.

p = K ɛvol,

τ2 , ν

μ1

μ1

τ2 , ν

(c)

τ1 ,ν

(d)

Fig. 1. Mechanical representations of fractional models: (a) fractional Burgers model, (b) fractional Kelvin-Voigt model, (c) fractional Zener model, and (d) fractional Maxwell model.

s (t ) = 2

∫0

t

R (t − t ′)

de (t ′) dt ′ dt ′

(2)

assuming s (t ) = 0 for t < 0. In Eq. (2), R(t) refers to the relaxation modulus, with R = G in case of elastic behavior. The deviatoric relaxation modulus R(t) depends on the underlying viscoelastic model, reading for the fractional Burgers model considered in this paper (Mainardi and Spada, 2011): ν

2

⎛ t ⎞⎤ R (t ) = 2μ1 ∑ Ri Eν ⎡ ⎢−⎜ τ0 ρ ⎟ ⎥ i=1 ⎣ ⎝ i⎠ ⎦

(3)

where the following abbreviations are used:

μ1 , μ2

rμ =

rτ =

τ1 , τ2

τ0ν =

(τ1 τ2) ν τ1ν + τ2ν

rμ rτν rτν 1 1 , δ1 = , + 1 + rμ (1 + rτν )2 1 + rμ 1 + rμ 1 + rτν 1 rτ , γ0 = δ2 = 1 + rμ 1 + rτ

δ0 =

δ12 1 δ1 δ ± − 0 , ρ1ν = − , 2δ2 δ2 z1 4δ22 z2 + γ0 z1 + γ0 , R2 = , R1 = z2 − z1 z1 − z2

z1,2 = −

ρ2ν = −

1 , z2 (4)

with μ1, μ2 [Pa], τ1, τ2 [s], and ν [-] referring to the parameters of the springs and fractional dashpots as illustrated in Fig. 1. The MittagLeffler function employed in Eq. (3)is defined as (Mainardi and Spada, 2011): ∞

s = 2 Ge,

(b)

μ2

Eα (z ) =

For many materials, the viscoelastic behavior is mainly driven by the deviatoric part of the stress state (Lakes and Wineman, 2006). Therefore, constitutive equations are split into bulk and shear behavior, linking volumetric and deviatoric stresses to the respective strains, reading in case of isotropic elastic behavior:

⇔

μ2

τ2 , ν

(a)

2. Fractional viscoelasticity

σ = : ɛ

μ2

τ1 ,ν

μ1

∑ k=0

zk , α > 0, z ∈  Γ(αk + 1)

(5)

Considering 2μ1 being equal to the elastic shear modulus G, the dimensionless relaxation modulus r(t) is introduced as

r (t ) =

(1)

with K and G as bulk and shear modulus, respectively. In Eq. (1), 1 p = 3 (σ11 + σ22 + σ33), ɛvol = ɛ11 + ɛ22 + ɛ33, and s and e represent the stress and strain deviator, respectively. In case of viscoelastic deviatoric behavior, the respective stress history is expressed by the convolution integral, reading (Christensen, 1982)

R (t ) = G

2

ν

t ⎞⎤ ⎟ ⎥ ⎣ ⎝ 0 ρi ⎠ ⎦

∑ Ri Eν ⎡⎢−⎛ τ ⎜

i=1

(6)

3. Viscoelastic extension of Chamis-equations The Chamis-equations provide the effective elastic properties (Young’s moduli, shear moduli, Poisson’s ratios) for a unidirectionallyreinforced composite (Chamis, 1989):

15

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al. eff E11 = EF 11 VF + (1 − VF ) EM EM eff eff E22 = E33 = 1 − VF + VF EM / EF 22 GM eff eff G12 = G13 = 1 − VF + VF GM / GF 12 GM eff = G23 1 − VF + VF GM / GF 23

with kF = KF / GM and gF = GF / GM . Inserting Eqs. (10) and (11) into the Chamis-equations (7) gives the effective material properties in dimensionless form as: eff kF gF E11 kM r (t ) = 9VF + 9(1 − VF ) 3kM + r (t ) GM 3kF + gF

eff eff ν12 = ν13 = VF νF 12 + (1 − VF ) νM eff ν23 =

eff E22 eff 2G23

(7)

eff eff eff 0 0 0 ⎤ ⎡ C11 C12 C12 ⎥ E11 ⎫ ⎢ eff eff C C 0 0 0 ⎥⎧ 22 23 ⎢ ⎪ E22 ⎪ eff ⎢ C22 0 0 0 ⎥ ⎪ E33 ⎪ ⎥ ⎢ = eff ⎨ Σ12 ⎬ ⎢ C44 0 0 ⎥ ⎨ 2E12 ⎬ ⎪ Σ13 ⎪ ⎢ ⎥ ⎪ 2E ⎪ eff C44 0 ⎥ ⎪ 13 ⎪ ⎪ ⎪ ⎢ 2E32 ⎭ Σ ⎩ 32 ⎭ ⎢ eff ⎥ ⎩ C66 ⎦ ⎣ sym

eff D23

= −

eff D12 =−

,

eff ν23 eff E22

,

eff D44

=

eff ν12 eff E11

1 eff G12

,

eff D22 =

,

eff D66

=

1 eff E22

⎟

VF ) gF +

−1

⎛⎜ 1 − VF + VF ⎞⎟ ⎤ kF ⎠ ⎥ VF r (t )) ⎝ kM ⎦

In order to extend the Chamis-equations towards viscoelasticity, the socalled correspondence principle (Lee, 1955) is exploited, where the viscoelastic problem is solved by considering its corresponding elastic solution in the Laplace–Carson domain. Accordingly, all parameters appearing in Eq. (12) are replaced by the respective transformed expressions. While constant parameters are not affected by the Laplace–Carson transformation, i.e. •* =  [•] = •, the transformation of the relaxation modulus of the fractional Burgers model (see Eq. (6)) becomes: 2

r * (p) =  [r (t )] =

(pρ τ ) ν

∑ Ri (pρ τ i) ν0 + 1

(13)

i 0

i=1

Inserting Eq. (13) into Eq. (12) gives – after algebraic manipulations – the effective dimensionless viscoelastic properties in the Laplace–Carson domain in the form:

,

•eff * (p) = c0 +

1 (9)

c1 + c2 pν + c3 p2ν c4 + c5 pν + c6 p2ν c1 c3 c4 c6

c2 ν p + p2ν c3 c5 ν + c p + p2ν 6 c1 c2 ν + p + p2ν c3 c3 (pν − p1ν )(pν − p2ν ) c1 c + c2 pν + p2ν c3 3 ⎛ ⎜ ν p2ν − p1ν ⎝p

+

= c0 +

c3 c6

= c0 +

c3 c6

EM =

= c0 +

c3 c6

νM

= c0 +

1 1 1 c1 ⎛ ⎞ − ν ⎜ ⎟ c6 p2ν − p1ν ⎝ pν − p2ν p − p1ν ⎠

Considering deviatoric viscoelastic behavior of the matrix material, the parameters EM and νM in Eq. (7) are expressed in terms of the constant bulk modulus KM and the respective relaxation modulus RM(t) replacing the shear modulus GM, reading

EM (t ) 9KM RM (t ) 9kM r (t ) 9KM GM → EM (t ) = = , GM 3KM + GM 3KM + RM (t ) 3kM + r (t ) 3KM − 2RM (t ) 3KM − 2GM = → νM (t ) = , 2(3KM + RM (t )) 2(3KM + GM ) 3kM − 2r (t ) νM (t ) = , 2(3kM + r (t ))

ν pν ⎞ 1 c2 ⎛ p − ν ⎜ ⎟ c6 p2ν − p1ν ⎝ pν − p2ν p − p1ν ⎠

+

2ν p2ν ⎞ 1 c3 ⎛ p − ν ⎜ ⎟ c6 p2ν − p1ν ⎝ pν − p2ν p − p1ν ⎠

EF νF

where the roots

9KF GF = , 3KF + GF

p1ν

and

c4 c + 5 pν + p2ν = 0 c6 c6 (11)

1 1 ⎞ − ν ⎟ − p2ν p − p1ν ⎠

+

(10) where the dimensionless quantities kM = KM / GM and r (t ) = R (t )/ GM were used. For the fibers, considering isotropic elastic behavior, Young’s modulus and Poisson’s ratio become

9kF gF EF = GM 3kF + gF 3kF − 2gF 3KF − 2GF = , νF = 2(3KF + GF ) 2(3kF + gF )

−1 (12)

(8)

eff G23

2gF r (t )

⎡2 eff ν23 = ⎢ + 3 9((1 − ⎣

where, due to the transversal isotropy (perpendicular to the fiber dieff eff eff eff eff eff = C12 , C33 = C22 , and C55 = C44 rection), C13 were considered. The components of the material compliance matrix  , with  =  −1, are linked to the engineering parameters given in Eq. (7) as

eff E11

⎜

−1

⎧ Σ11 ⎫ ⎪ Σ22 ⎪ ⎪ Σ33 ⎪

1

1 1 ⎞⎤ VF ) ⎛ + 3r (t ) ⎠ ⎥ ⎝ 9kM ⎦

eff VF 1 − VF ⎞ G eff G12 = 23 = ⎛⎜ + ⎟ GM GM g r (t ) ⎠ ⎝ F 3kF − 2gF 3kM − 2r (t ) eff ν12 = VF + (1 − VF ) 2(3kF + gF ) 2(3kM + r (t ))

− 1,

where the indices M and F refer to quantities of the matrix and fibers, respectively, considering possible fiber-anisotropy as observed e.g. for carbon fibers (Kriz and Stinchcomb, 1979); VF denotes the volume fraction of the fibers, E, G, and ν the Young’s modulus, shear modulus, and Poisson’s ratio, respectively. The effective elastic properties given in Eq. (7) link the homogenized strain tensor E with the homogenized stress tensor Σ, reading in matrix notation:

eff D11 =

−1

eff E22 1 1 ⎞ ⎡ = ⎢ VF ⎜⎛ + ⎟ + (1 − GM 3gF ⎠ ⎝ 9kF ⎣

p2ν

→

of the quadratic equation are given by ν p1,2 =−

1 c5 ± 2 c6

and the constants c0, c1…c6 are given for

16

(14)

1 c52 c − 4, 4 c62 c6

(15)

U. Hofer et al.

Mechanics of Materials 127 (2018) 14–25

eff E11 : GM

after simplifications – to the effective dimensionless viscoelastic properties in the time domain1:

•eff =

c0 = 9

kF gF VF 3kF + gF

•eff (t ) =  −1[•eff * (p)]

c1 = 0 c2 = 9(1 − VF ) δ2 γ0 kM τ0ν

− t νEν,1 + ν (p1ν t ν )]

c3 = 9(1 − VF ) δ2 kM τ02ν c4 = 3δ0 kM c5 = (δ2 γ0 + 3δ1 kM ) τ0ν

•eff =

+

c2/ c6 [Eν (p2ν t ν ) − Eν (p1ν t ν )] p2ν − p1ν

+

c3/ c6 [t −νEν,1 − ν (p2ν t ν ) − t −νEν,1 − ν (p1ν t ν )] p2ν − p1ν

3kM ) τ02ν

c6 = δ2 (1 + eff E22

c1/ c6 [t νEν,1 + ν (p2ν t ν ) p2ν − p1ν

= c0 +

1 ⎞ c1/ c6 ⎡ 1 ⎛ Eν (p2ν t ν ) − ν Γ(1) p2ν − p1ν ⎢ p ⎠ ⎣ 2⎝

= c0 +

:

GM c0 = 0 c1 = 0 c2 = 9δ2 γ0 gF kF kM τ0ν VF ) δ0 gF kF kM

c5 = {3(1 −

⎟

1 ⎛ 1 ⎞⎤ Eν (p1ν t ν ) − Γ(1) ⎠ ⎥ p1ν ⎝ ⎦ c /c + ν 2 6 ν [Eν (p2ν t ν ) − Eν (p1ν t ν )] p2 − p1 −

c3 = 9δ2 gF kF kM τ02ν c4 = 3(1 −

⎜

VF ) δ1 gF kF kM + δ2 γ0 [gF kF +

+

VF (gF kM + 3kF kM

− gF kF )]} τ0ν

⎜

⎟

1 c3/ c6 ⎡ ⎛ ν ⎞ p2 Eν (p2ν t ν ) + Γ(1 − ν ) ⎠ p2ν − p1ν ⎢ ⎣⎝ ⎜

⎟

1 ⎞⎤ − ⎛p1ν Eν (p1ν t ν ) + Γ(1 − ν) ⎠ ⎥ ⎝ ⎦ ⎜

c6 = δ2 {gF kF (1 + 3kM ) + [3kF kM − gF (kF − kM

c1/ c6 + c2/ c6 p2ν + c3/ c6 p22ν ⎞ = c0 + Eν (p2ν t ν ) ⎜⎛ ⎟ p2ν (p2ν − p1ν ) ⎠ ⎝ c / c + c2/ c6 p1ν + c3/ c6 p12ν ⎞ ν ν ⎛ 1 6 − Eν (p1 t ) ⎜ ⎟ p1ν (p2ν − p1ν ) ⎝ ⎠ c /c + 1ν 6ν p1 p2

+ 3kF kM )] VF } τ02ν •eff =

eff G12

: GM c0 = 0 c1 = 0 c2 = δ2 γ0 gF τ0ν c3 = δ2 gF τ02ν c4 = (1 −

(17)

VF ) δ0 gF

c5 = [(1 −

VF ) δ1 gF +

c6 = δ2 [(1 −

VF ) gF +

using the following property (Haubold et al., 2011)

VF δ2 γ0] τ0ν VF ] τ02ν

Eα, β (z ) = zEα, α + β (z ) +

eff •eff = ν12 :

∞

Eα, β (z ) =

−

c0 = − 1

c3 = 9δ2 kF kM [(1 − c4 = 6(1 −

VF ) δ1 gF + VF ) gF −

VF δ2 γ0] τ0ν VF ] τ02ν

VF ) δ1 gF kF kM + δ2 γ0 [(1 −

VF ) gF kF

VF (gF kM + 3kF kM )]} τ0ν

c6 = 2δ2 [(1 −

(19)

VF ) gF kF (1 + 3kM ) +

c1/ c6 + c2/ c6 p2ν + c3/ c6 p22ν p2ν (p2ν − p1ν ) c1/ c6 + c2/ c6 p1ν + c3/ c6 p12ν p1ν (p2ν − p1ν ) c3 c6

+

c1/ c6 p1ν p2ν (20)

where Eν (t = 0) = 1 was considered, retrieving Eq. (12) for r (t ) = 1. Eq. (17) covers as well the more simple models given in Fig. 1: setting τ1 = ∞ gives the fractional Zener model, τ1 = ∞ and μ1 = ∞ the fractional Kelvin-Voigt model, whereas the fractional Maxwell model is

VF ) δ0 gF kF kM

c5 = 2{3(1 − +

= c0 +

VF ) δ0 gF kF kM

c2 = 9kF kM [(1 −

zk , α > 0, β > 0, z ∈  , Γ(αk + β )

•eff (t = 0) = c0 +

eff : ν23

c1 = 9(1 −

(18)

From Eq. (17) the elastic case is recovered for t = 0, giving

c6 = 2δ2 (1 + 3kM ) τ02ν =

∑ k=0

c3 = δ2 (1 − VF )(3kM − 2) τ02ν c4 = 6δ0 kM c5 = 2(δ2 γ0 + 3δ1 kM ) τ0ν •eff

1 Γ(β )

of the generalized Mittag-Leffler function (Mainardi and Spada, 2011):

3kF − 2gF 2(3kF + gF ) = 3δ0 kM (1 − VF ) = (1 − VF )(3δ1 kM − 2δ2 γ0) τ0ν

c0 = VF c1 c2

⎟

1

For back-transformation of Eq. (14) into the time domain, the following transformation rules were employed (Mainardi and Spada, 2011; Rogosin, 2015):

VF (3kF kM + gF kM )] τ02ν

Eq. (14) is transformed term-wise into the time domain, leading –

0

p ⎤ 1 ⎤ =  −1⎡ α = t αEα,1 + α (λt α )  −1⎡ α ⎢ ⎢ ⎣ p − λ⎥ ⎦ ⎣ p − λ⎥ ⎦ pα ⎤ = Eα (λt α )  −1⎡ α ⎢ − λ⎥ p ⎣ ⎦ p2α ⎤ = t −αEα,1 − α (λt α )  −1⎡ α ⎢ ⎣ p − λ⎥ ⎦

17

(16)

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

As regards the latter, UCs representing different types of packing, i.e. diamond, square, and hexagonal, are considered, as well as different values of the fiber volume fraction (see Fig. 2). Periodic boundary conditions are applied to ensure deformation, stress, and strain periodicity (Anthoine, 1995) by coupling surface nodes of the UC to respective reference nodes by constraining the displacements as:

obtained by τ2 = ∞ and μ 2 = ∞; the respective relaxation moduli in the Laplace–Carson and time domain as well as the parameters c0, c1…c6 for these models are given in Appendix A. 4. Model assessment The closed-form expressions for the extended Chamis-equations presented in the previous section are assessed by

• •

ui(j+) = ui(j−) + ujiREF ,

the equivalent Eigenstrain method, specialized for cylindrical inclusions and fractional viscoelasticity (see Appendix B) and numerical simulations employing a unit-cell (UC) model.

(21)

ui(j+)

is the i-th component of the displacement vector at the j-th where positive surface, ui(j−) at the opposite side, and ujiREF is the i-th component of the displacement vector at the reference node for the j-th

Fig. 2. Employed UC models: (a-c) diamond, (d-f) square, and (g-k) hexagonal packing, illustrated for the fiber volume fractions under investigation, (l): definition of geometric properties of the UC.

18

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

square and diamond UCs bound the results of the hexagonal UC. F = 0.75, close to the theoretical maximum VF ≈ 0.785, the square and diamond UCs show locking behavior and thus overestimate the stiffness: the square UC for the eff eff , the diamond UC for the shear modulus G23 , transverse modulus E22 eff and both UCs for the shear modulus G12 . The hexagonal cell, on the other hand, although being close to the theoretical maximum of ≈ 0.907, seems to produce realistic results also for a fiber volume fraction of VF = 0.90, presumably due to the fine discretization. Therefore, the hexagonal UC is assumed to be the most realistic and is, in the following, used as reference model.

Table 1 Coefficients of Prony series. i

ri

τi

1

[·10−3 ] 2.37

1.59·10−3

2

3.78

2.82·10−2

3

6.52

4 5 6 7 8

13.04 25.78 49.12 86.95 131.38

5.00·10−1 8.83 1.55 · 102 2.68 · 103 4.55 · 104 7.49 · 105

• For a fiber volume fraction of V

[s]

In Figs. 4 to 6, the effective relaxation functions obtained from the extended Chamis-equations, Eigenstrain method, and the hexagonal UC – related to the respective instantaneous value at t = 0 – are shown. The following observations are made:

REF REF REF = u31 = u23 = 0, surface. The UC is additionally constraint by u21 in order to avoid rigid body motions. The homogenized strain tensor E, depending solely on the displacement vectors of the three reference nodes and the UC geometry, is given by

Eij =

uijREF Li

,

• For the longitudinal Young’s modulus E

eff 11

(22)

•

where Li denotes the i-th dimension of the UC (see Fig. 2(l)). The homogenized stress tensor Σ is determined from the reaction forces at the reference nodes, reading

Σij =

RjiREF Aj

•

(23)

with RjiREF denoting the i-th component of the reaction force at the reference node of the j-th surface and Aj as the area of the j-th surface (see Fig. 2(l)). Displacements are prescribed at the reference nodes in such a way, that only one component of the homogenized strain tensor E is non-zero. Accordingly, six load cases are required in order to determine eff , the six non-zero components of the effective stiffness matrix (C11 eff C12 ... , see Eq. (8)) and, thereafter, the effective Young’s moduli, shear moduli, and Poisson’s ratios from  =  −1 and Eq. (9). For the assessment of the closed-form expressions derived in the previous section, a matrix material following a fractional Maxwell model is considered, with τ1 = 5·107 s and ν = 0.25 (see Fig. 1(d)). Within the numerical simulations, the fractional matrix behavior is approximated by a Prony series, reading n

r (t ) = 1 −

i=1

eff To summarize, for the longitudinal Young’s modulus E11 and eff , the results obtained from the analytical models and Poisson’s ratio ν12 eff , the numerical model almost coincide. For the transversal modulus E22 eff eff eff the shear moduli G12 and G23 and Poisson’s ratio ν23 , on the other hand, deviations between the extended Chamis-equations, the Eigenstrain method and the numerical simulation can be observed, being attributed to the model representation of the material’s morphology perpendicular to the fiber direction, i.e. the arrangement of the fibers:

t

∑ ri (1 − e− τi ),

(24)

• In case of the Chamis-equation, the fibers and the interjacent matrix are assumed to act as serially-connected springs. • Deriving the equivalent Eigenstrain method, a unit-cell representing

where n = 8 is chosen in this paper with the respective model parameters given in Table 1. The dimensionless material parameters are set equal to kM = 3, kF = 40, and gF = 30, representing realistic values e.g. for glass fibers embedded in an epoxy matrix. As regards the effective elastic behavior, the following observations are made (see Fig. 3):

•

• The elastic properties obtained from both the Chamis-equations and • •

and the Poisson’s ratio the analytically-obtained results coincide with the numericallyeff , reobtained UC-based results2. Regarding the shear modulus G23 sults obtained from both analytical methods slightly overestimate the relaxation behavior obtained from the UC-model. eff , the extended Chamis-equations predict more proAs regards E22 nounced relaxation behavior than the Eigenstrain method, with the results obtained by the Eigenstrain method being closer to the numerically-obtained results. eff eff , the extended Chamis-equations and the Eigenstrain For G12 and ν23 eff , the extended method predict different results. While for G12 Chamis-equations predict a more pronounced relaxation behavior, eff the contrarious behavior is observed for ν23 . For low values of the fiber volume fraction, the Eigenstrain method is closer to the UCmodel response, whereas for high fiber volume fractions, the relaxation curves obtained from the extended Chamis-equations are closer to the results obtained from the hexagonal UC. eff ν12 ,

Eigenstrain method agree well with the elastic properties obtained from the UCs, being closest for the hexagonal UC. For Poisson’s ratio eff ν23 , the Eigenstrain method performs better for lower, the extended Chamis-equations for higher values of the fiber volume fraction. Regarding the extended Chamis-equations, the fiber volume fraction eff eff has a linear influence on E11 and ν12 and a nonlinear influence on the other moduli, corresponding to the appearance of VF in Eq. (7). eff eff eff , and ν12 , G12 , square and diamond UCs lead to the same For E11 eff eff eff , on the other hand, the results of results. For E22 , G23 , and ν23

a circular cross-section of the fiber within a rectangular matrix block is assumed, resulting in a square fiber arrangement (see Fig. 2 d–f) with the fiber and the matrix acting as serially-connected springs. Finally, within the numerical simulations – employing a hexagonal UC (see Fig. 2 g–k) – the fibers and the interjacent matrix are acting as mixed parallel/serially-connected springs.

The different representations of fiber arrangement in the transverse direction significantly influence the homogenized behavior. Fig. 3 shows the large scatter within the UC-based elastic results, illustrating the strong influence of fiber arrangement (diamond, square, or hexagonal, see Fig. 2 a–k) onto the results of the transverse effective properties.

2 Observable deviations between analytical and UC-based results are explained by approximating the fractional Maxwell model by a Prony series in the FE-simulations (see Eq. (24) and Table 1).

19

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

Fig. 3. Chamis-equations vs. Eigenstrain method and UC-model: effective elastic material properties for varying values of VF, related to the respective property of the matrix.

Fig. 4. Extended Chamis-equations vs. Eigenstrain method and UC-model (hexagonal packing) for VF = 0.25: effective relaxations functions, related to their respective instantaneous value at t = 0 .

20

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

Fig. 5. Extended Chamis-equations vs. Eigenstrain method and UC-model (hexagonal packing) for VF = 0.50 : effective relaxations functions, related to their respective instantaneous value at t = 0 .

Fig. 6. Extended Chamis-equations vs. Eigenstrain method and UC-model (hexagonal packing) for VF = 0.75: effective relaxations functions, related to their respective instantaneous value at t = 0 .

21

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

5. Concluding remarks

similar quality in comparison to numerically-obtained curves. Due to the closed-form character of the extended Chamis-equations, application is straight-forward and computationally cheap, not requiring a – in general numerically-performed – inverse Laplace–Carson transformation. The consideration of fractional viscoelasticity offers the opportunity to describe a large variety of viscoelastic materials.

Closed-form expressions for the effective viscoelastic relaxation moduli of unidirectionally-reinforced composites were presented by extending the Chamis-equations, considering fractional viscoelastic behavior of the matrix material. The performance of the so-obtained closed-form expressions was assessed by alternative methods such as the Eigenstrain method and numerical simulations, assuming isotropic elastic behavior for the fibers and deviatoric viscoelastic matrix behavior following a fractional Maxwell model. In general, the extended Chamis-equations predict the effective elastic and viscoelastic properties well, compared to numerical (UCbased) simulations. For lower values of the fiber volume fraction, the effective relaxation curves obtained by the Eigenstrain method are closer to the numerically-obtained results. For higher values of the fiber volume fraction, on the other hand, the predicted effective relaxation curves from both Chamis-equations and Eigenstrain method are of

Acknowledgments This research was funded by the K-Regio project “Innovative Tube Design - Entwicklung und Optimierung von neuen High-Tech Faserverbundstrukturen für den industriellen Einsatz auf Basis modellund simulationsbasierter Methoden” in cooperation with Thöni Industriebetriebe GmbH, superTEX composites GmbH, Intales GmbH. Financial support by the Tyrolean Government, the European Regional Development Fund (ERDF) and the University of Innsbruck is gratefully acknowledged.

Appendix A. Extended Chamis-equations for fractional Zener, Maxwell, and Kelvin-Voigt model The relaxation moduli in the Laplace–Carson domain for the fractional Zener, Maxwell, and Kelvin-Voigt model, respectively, are obtained from r*(p) in Eq. (13) as

* (p) = lim r * (p) = rZE

μ2 ⎡ μ (pτ2) ν ⎤ 1+ 1 ν + (1 + μ / μ ) ⎥ ( ) μ1 + μ 2 ⎢ μ pτ 2 1 2 ⎦ 2 ⎣

* (p) = rMX

(pτ1) ν 1 + (pτ1) ν

τ 1→∞

* (p) = rKV

lim r * (p) =

τ2→∞ μ2 →∞

lim r * (p) = 1 + (pτ2) ν

τ 1→∞ μ1 →∞

(A.1)

By back-transformation, the dimensionless relaxation moduli are obtained in the time domain, reading:

rZE (t ) =

ν μ2 ⎡ μ1 ⎡ ⎛ t ⎞ ⎤⎤ + E − + μ μ 1 (1 / ) ν 1 2 ⎥ μ1 + μ 2 ⎢ μ2 ⎢ ⎝ τ2 ⎠ ⎥ ⎣ ⎦⎦ ⎣ ⎜

⎟

ν

t rMX (t ) = Eν ⎡−⎛ ⎞ ⎤ ⎢ ⎝ τ1 ⎠ ⎥ ⎣ ⎦ (t / τ2)−ν rKV (t ) = 1 + Γ(1 − ν ) ⎜

⎟

(A.2)

The coefficients c0, c1…c6, introduced in the extended Chamis-equations (14) for the effective viscoelastic properties, are specialized in Tables A1 to A3 for the fractional Zener, Maxwell, and Kelvin-Voigt model. Since the coefficients c3 and c6 are always zero, Eq. (14) reduces to

c1 + c2 pν c4 + c5 pν pν c1 1 c = c0 + + 2 c5 c4/ c5 + pν c5 c4/ c5 + pν

•eff * (p) = c0 +

(A.3)

which, back-transformed to time domain, becomes

•eff (t ) =  −1[•eff * (p)] = c0 + = c0 +

c1 ν c c c t E ν, ν + 1 ⎛ − 4 t ν ⎞ + 2 E ν ⎛ − 4 t ν ⎞ c5 c5 ⎝ c5 ⎠ ⎝ c5 ⎠ ⎜

⎟

⎜

⎟

c1 c c c + ⎛ 2 − 1 ⎞ Eν ⎛− 4 t ν ⎞ c4 c4 ⎠ ⎝ c5 ⎠ ⎝ c5 ⎜

⎟

⎜

⎟

(A.4)

22

9kM (1 − VF )

3kM (1 + rμ) + 1

1 + 3kM

c4

c5 τ2ν

gF kF (1 −

gF kF (1 −

9gFkFkM

9gFkFkM

0

eff E22

VF )·[1 + 3kM (1 + rμ)] + (gF + 3kF ) kM VF

23

0

9kM (1 − VF )

3kM

1 + 3kM

c4

c5 τ1ν

gF kF gF + 3kF

c2 τ1ν

9VF

c1

c0

eff E11

(gF kF + 3gF kF kM )·(1 −

3gF kF kM (1 −

9gFkFkM

0

0

eff E22

VF ) VF ) + (gF + 3kF ) kM VF

VF )·(1 + 3kM ) + (gF + 3kF ) kM VF

Table A2 Chamis-coefficients for fractional Maxwell model .

9kM (1 − VF )

gF kF gF + 3kF

c2 τ2ν

9VF

c1

c0

eff E11

Table A1 Chamis-coefficients for fractional Zener model.

gF (1 −

gF (1 −

gF

gF

0

eff G12

gF (1 −

gF (1 −

gF

0

0

eff G12

VF ) +

VF )

VF

VF )(1 + rμ) + VF ) +

VF

VF

VF (2gF − 3kF ) 2 (gF + 3kF )

VF (2gF − 3kF ) 2 (gF + 3kF )

2 + 6kM

6kM

(3kM − 2)(1 − VF )

3kM (1 − VF )

−

eff ν12

2(1 + 3kM )

2[1 + 3kM ·(rμ + 1)]

(3kM − 2)·(1 − VF )

[3kM − 2(rμ + 1)]·(1 − VF )

−

eff ν12

VF ) +

VF ]

VF ) +

VF ) + 2(gF + 3kF ) kM VF

VF )]

VF )

VF ) +

VF )

VF ) + 2(gF + 3kF ) kM VF

VF ]

VF ) + 2(gF + 3kF ) kM VF

(2gF kF + 6gF kF kM )·(1 −

6gF kF kM (1 −

9kF kM [gF (1 −

9gF kF kM (1 −

−1

eff ν23

[2gF kF (1 + 3kM )]·(1 −

[6gF kF kM (1 + rμ)+2gF kF ](1 −

9kF kM ·[gF (1 −

9kF kM [gF (rμ + 1)·(1 −

−1

eff ν23

U. Hofer et al.

Mechanics of Materials 127 (2018) 14–25

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

Table A3 Chamis-coefficients for fractional Kelvin-Voigt model . eff E11

c0

gF kF gF + 3kF

9VF

eff E22

eff G12

0

0

eff ν12

−

VF (2gF − 3kF ) 2 (gF + 3kF )

eff ν23

−1

c1

9kM (1 − VF )

9gFkFkM

gF

(3kM − 2)·(1 − VF )

9kF kM [ VF +gF (1 −

c2 τ2ν

9kM (1 − VF )

9gFkFkM

gF

− 2(1 − VF )

9kF kM VF

c4

1 + 3kM

gF kF (1 + 3kM )·(1 −

(1 − gF ) VF + gF

2 + 6kM

(2gF kF + 6gF kF kM )·(1 −

c5 τ2ν

1

gF kF (1 −

2

2gF kF (1 −

VF ) + kM (gF + 3kF ) VF

VF ) + kM (gF + 3kF ) VF

VF

VF )]

VF ) + 2(gF + 3kF ) kM VF

VF ) + 2(gF + 3kF ) kM VF

Appendix B. Equivalent Eigenstrain method The equivalent Eigenstrain method was originally formulated for elliptical inclusions (Eshelby, 1957), later on adapted for consideration of cylindrical inclusions (Nemat-Nasser et al., 1982), in order to simulate fibers embedded in a matrix, and extended towards determination of the effective viscoelastic behavior by applying the correspondence principle (Luciano and Barbero, 1995). Considering deviatoric viscoelastic behavior only (setting the bulk modulus constant) and employing a fractional Maxwell model to represent the viscoelastic behavior of the matrix, the expressions for the effective properties – originally formulated in terms of the Lamé parameters λ and μ (see (Luciano and Barbero, 1995)) – are expressed in terms of bulk and shear modulus, reading in the Laplace–Carson domain: eff * 2 (p) 2L12 E11 = L11 − L22 + L23 GM eff * 2 (p) E22 [2L11 (L22 + L23) − 4L12 ](L22 − L23 + 2L44 ) = 2 GM L11 (3L22 + L23 + 2L44 ) − 4L12 eff * (p) G12 1 (L22 − L23 + 2L44 ) = 4 GM eff * (p) G23 = L66 GM L12 eff * (p) = ν12 L22 + L23 eff * (p) = ν23

2 L11 (L22 + 3L23 − 2L44 ) − 4L12 2 L11 (3L22 + L23 + 2L44 ) − 4L12

(B.1)

with

* (p) gS 2 − 2S3 S6 S 2 − S72 ⎞ 1 bS7 + aS6 − agS3 4rM a2 − b2 + 3 + + 26 − VF ⎛⎜ ⎟ 2 2 * * * (p)2 ⎠ D grM (p) cgrM (p) g rM 3 ⎝ 4c b (gS3 − S6 + S7) 2r * (p) ab + b2 ⎤ 1 + VF ⎡ − L12 (p) = kM − M ⎢ * (p) 3 2cgrM 4c 2 ⎥ ⎣ ⎦D L11 (p) = kM +

* (p) 2rM aS7 ab + b2 ⎞ 1 + VF ⎛⎜ − ⎟ * 3 2 cgr ( p ) 4c 2 ⎠ D M ⎝ 4r * (p) aS6 − agS3 ⎞ 1 a2 − b2 + − VF ⎛⎜ L22 (p) = kM + M ⎟ 2 * (p) ⎠ D 3 2cgrM ⎝ 4c L23 (p) = kM −

−1

* (p) + rM * (p) VF ⎡ L44 (p) = rM ⎢ ⎣ * (p) + rM * (p) VF L66 (p) = rM

* (p) S3 * (p) + 2gF S3 − 2rM * (p)) S7 ⎤ rM 4(3kM + rM − * (p) * (p) ⎥ 3kM + 4rM gF − rM ⎦

* (p) gF − rM * (p) gF S3 + (1 − S3) rM

a (S62 − S72 − 2gS3 S6) (a2 − b2) S6 + (ab + b2S7) (b2 − a2) S3 a3 − 3ab2 − 2b3 aS32 + + + + * (p) * (p) * (p)2 * (p)2 8c 3 2c 2grM 2c 2rM 2cg 2rM 2crM * (p)) a (p) = 6(−gF − 3kF + 3kM + rM * (p) b (p) = 6gF − 9kF + 9kM − 6rM * c (p) = 27(kF − kM )(gF − rM (p))

D (p) =

g (p) =

* (p) 3kM + 4rM * (p) 3kM + rM

S3 = 0.49247 − 0.47603VF − 0.02748VF2 S6 = 0.36844 − 0.14944VF − 0.27152VF2 S7 = 0.12346 − 0.32035VF + 0.23517VF2

(B.2)

24

Mechanics of Materials 127 (2018) 14–25

U. Hofer et al.

* (p) is the Laplace-Carson transformed relaxation modulus of the matrix material. The relaxation modulus of the underlying fractional where rM * (p) in Eq. (A.1)) is inserted into Eq. (B.2) and, finally, into Eq. (B.1). Back-transformation of the so-obtained effective Maxwell model (see rMX properties to the time domain is performed numerically, using Talbot’s algorithm (Abate and Whitt, 2006). The respective instantaneous value of the * (p) = 1 in Eq. (B.2). relaxation function is obtained by setting rM

Lakes, R.S., Wineman, A., 2006. On Poisson’s ratio in linearly viscoelastic solids. J. Elast. 85 (1), 45–63. Laws, N., McLaughlin, R., 1978. Self-consistent estimates for the viscoelastic creep compliances of composite materials. Proc. Royal Soc. Lond. A Math. Phys. Eng. Sci. 359 (1697), 251–273. Lee, E.H., 1955. Stress analysis in visco-elastic bodies. Q Top Q. Appl. Math. 13 (2), 183–190. Luciano, R., Barbero, E.J., 1995. Analytical expressions for the relaxation moduli of linear viscoelastic composites with periodic microstructure. J. Appl. Mech. 62 (3), 786–793. Mainardi, F., 2010. Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models. Imperial College Press, London. Mainardi, F., Spada, G., 2011. Creep, relaxation and viscosity properties for basic fractional models in rheology. Eur. Phys. J. Spec. Top. 193 (1), 133–160. Meral, F.C., Royston, T.J., Magin, R., 2010. Fractional calculus in viscoelasticity: an experimental study. Commun. Nonlinear Sci. Numer. Simul. 15 (4), 939–945. Nemat-Nasser, S., Hori, M., 1993. Micromechanics: overall properties of heterogeneous materials. Elsevier, Amsterdam. Nemat-Nasser, S., Iwakuma, T., Hejazi, M., 1982. On composites with periodic structure. Mech. Mater. 1 (3), 239–267. Oeser, M., Pellinen, T., Scarpas, T., Kasbergen, C., 2008. Studies on creep and recovery of rheological bodies based upon conventional and fractional formulations and their application on asphalt mixture. Int. J. Pavement Eng. 9 (5), 373–386. Pichler, C., Lackner, R., Aigner, E., 2012. Generalized self-consistent scheme for upscaling of viscoelastic properties of highly-filled matrix-inclusion composites - application in the context of multiscale modeling of bituminous mixtures. Compos. Part B Eng. 43 (2), 457–464. Podlubny, I., 1999. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, San Diego. Rogosin, S., 2015. The role of the mittag-Leffler function in fractional modeling. Mathematics 3 (2), 368–381. Upadhyaya, P., Upadhyay, C.S., 1991. A three-dimensional micromechanical model to predict the viscoelastic behavior of woven composites. Compos. Struct. 93 (11), 2733–2739. Wang, G., Pindera, M.-J., 2016. Locally-exact homogenization of viscoelastic unidirectional composites. Mech. Mater. 103, 95–109. Yancey, R.N., Pindera, M.-J., 1990. Micromechanical analysis of the creep response of unidirectional composites. J. Eng. Mater Technol. 112 (2), 157–163.

References Abate, J., Whitt, W., 2006. A unified framework for numerically inverting laplace transforms. INFORMS J. Comput. 18 (4), 408–421. Anthoine, A., 1995. Derivation of the in-plane elastic characteristics of masonry through homogenization theory. Int. J. Solids Struct. 32 (2), 137–163. Barbero, E.J., Luciano, R., 1995. Micromechanical formulas for the relaxation tensor of linear viscoelastic composites with transversely isotropic fibers. Int. J. Solids Struct. 32 (13), 1859–1872. Budiansky, B., 1965. On the elastic moduli of some heterogeneous materials. J. Mech. Phys. Solids 13 (4), 223–227. Celauro, C., Fecarotti, C., Pirrotta, A., Collop, A.C., 2012. Experimental validation of a fractional model for creep/recovery testing of asphalt mixtures. Constr. Build. Mater. 36, 458–466. Chamis, C., 1989. Mechanics of composite materials: past, present, and future. J. Compos. Technol. Res. 11 (1), 3–14. Christensen, R.M., 1982. Theory of viscoelasticity: an introduction, 2nd ed. Academic Press, New York. Deseri, L., Di Paola, M., Zingales, M., 2014. Free energy and states of fractional-order hereditariness. Int. J. Solids Struct. 51 (18), 3156–3167. Dinzart, F., Lipiński, P., 2010. Self-consistent approach of the constitutive law of a twophase visco-elastic material described by fractional derivative models. Arch. Mech. 62 (2), 135–156. Drago, A.S., Pindera, M.-J., 2008. A locally exact homogenization theory for periodic microstructures with isotropic phases. J. Appl. Mech. 75 (5). 051010–01–14 Eshelby, J., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. Royal Soc. Lond. A Math. Phys. Eng. Sci. 241 (1226), 376–396. Hashin, Z., 1966. Viscoelastic fiber reinforced materials. Am. Inst. Aeronaut. Astronaut. (AIAA) J. 4 (8), 1411–1417. Hashin, Z., Rosen, B.W., 1964. The elastic moduli of fiber-reinforced materials. J. Appl. Mech. 31 (2), 223–232. Haubold, H.J., Mathai, A.M., Saxena, R.K., 2011. Mittag-Leffler functions and their applications. J. Appl. Math. 2011, 1–51. Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech. Phys. Solids 13 (4), 213–222. Kriz, R.D., Stinchcomb, W.W., 1979. Elastic moduli of transversely isotropic graphite fibers and their composites. Exp. Mech. 19 (2), 41–49.

25