Frequency-dependent scaling from mesoscale to macroscale in viscoelastic random composites

(a) Random microstructure

Mathematically, the random material is taken as a set of all the realizations B(ω) parametrized by sample events ω of the Ω space B={B(ω);ω∈Ω}.2.1Any realization B(ω) of the composite B={B(ω);ω∈Ω}, while spatially disordered (i.e. heterogeneous), follows deterministic laws of mechanics. Given our interest in scaling effects, we consider finite-size mesoscale domains of the random medium Bδ={Bδ(ω);ω∈Ω},δ=Ld,2.2where L is the domain of size, d is the typical microstructural length (e.g. grain size), in a two- or three-dimensional Euclidean space ED, D={2,3}, where D represents the dimensionality of the problem (two dimensional or three dimensional). One extreme case δ=1 signifies the description at the level of one grain, while the second extreme δ→∞ is the RVE limit, and, in general, one wants to know with what precision is the RVE attained at any finite δ for specific physical and geometric parameters.

In the world of real and man-made materials, viscoelastic microstructures are very diverse and complex, just one example being shown in figure 1a. The description of such materials is developed in terms of TRFs of material properties such as the relaxation modulus C_klmn (or, equivalently, the creep compliance tensor S_klmn) Cklmn:Ω×ED×T→Vor{Cklmn(ω,x,t); ω∈Ω,x∈B⊂ED,t∈T},2.3where V=S2(S2(ED)) (D=2 or 3) is the space of all symmetric fourth-rank tensors over ED; T is the time domain. The constitutive equations, ∀ω∈Ω, x∈B, are expressed in temporal Stieltjes convolution σkl(t)=∫−∞tCklmn(t−τ)ddτϵmn dτandϵkl(t)=∫−∞tSklmn(t−τ)ddτσmn dτ,2.4where t is the time and τ is the dummy variable. For simplicity of notation, we have written σ_kl(t) and ϵ_kl(t) for the TRF σkl(ω,x,t) and TRF ϵkl(ω,x,t).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

(a) Microstructure of a real viscoelastic composite (asphalt concrete mixture). (b) Realization of the random, two-phase chessboard on a 40×40 lattice. (c) Numerically generated micrograph of Gaussian correlated microstructure [33].

It is seen that we are interchangeably using the symbolic (C) and the subscript (C_i…) notations for tensors, as the need arises.

The random field C is defined over a probability space (Ω,F,P) and a material domain subset B_δ of ED, taking values in (a subset of) a finite-dimensional real Hilbert space V . Also, F is the σ-algebra assumed to be rich enough to support any random field encountered in applications, while P is the probability measure (equivalently, a probability distribution). That is, with C(ω,x,t):Ω×ED×T→V denoting the TRF, this function, for any fixed x0∈B, is measurable and, as commonly done in probability theory, C(ω,t) denotes a single realization of TRF, while C(x,t) stands for a random tensor at a given spatial location x and time t. Henceforth, in this subsection the explicit dependence of C on time is dropped as it is specified by time-independent parameters (g¯n, k¯m, τ_i) in the Prony series model (2.20). In other words, the viscoelastic properties of every material point are fixed and TRF can be characterized as the TRF of the stiffness tensor. Thus, we assume C to be second-order and mean-square continuous. Next, we let ⟨C(x)⟩ be the ensemble (or statistical) average (i.e. a mathematical expectation with respect to a given probability distribution) of the field at a given x and R(x,y):=⟨[C(x)−⟨C(x)⟩]⊗[C(y)−⟨C(y)⟩]⟩2.5be the two-point covariance function of C. Assuming these two functions to be invariant with respect to arbitrary translations implies that ⟨C(x)⟩∈V is constant in space, while R(x,y)∈V⊗V depends only on the difference h=x−y, i.e. C is wide-sense stationary (WSS).

Let K=O(3) be the group of rotations and reflections in ED, and let (V,g) be an orthogonal representation of K, g being the group action. Assume that, for all k∈K and for all x∈RD, we have ⟨C(kx)⟩=g(k)⟨C(x)⟩andR(kx)=g(k)R(x)g−1(k).2.6Such a field is called wide-sense isotropic. In what follows, we collect the key consequences of this assumption for the expectation and two-point correlation functions of several TRFs.

In the special case of anti-plane viscoelasticity, one has to take V=S2(E2), i.e. the space of all symmetric tensors of the second rank over E2, with the representation g(k)C=kCk−1. Then ⟨r(x)⟩=rI with I being the unit second-rank tensor. Using the theory of invariants, Lomakin [34] proved that the two-point correlation Rklmn(x,y):=⟨[Ckl(x)−⟨Ckl(x)⟩][Cmn(y)−⟨Cmn(x)⟩]⟩ is represented in terms of five continuous scalar functions K₁, …, K5:[0,∞)→R with K₃(0)=K₄(0)=K₅(0)=0, along with five corresponding tensor functions Lklmnn, as Rklmn(h)=∑n=15Lklmnn(h)Kn(∥h∥).2.7Malyarenko & Ostoja-Starzewski [35] found an equivalent representation using different functions Mklmnn(x) (n=1,…,5) which, while not as simple as those of Lomakin, lead to spectral expansions of tensor-valued homogeneous and isotropic random fields similar to those in [36].

In the case of in-plane or three-dimensional viscoelasticity, one has to work with V=S2(S2(ED)) (D=2 or 3), whose elements are fourth-rank relaxation tensors C_klmn. Its two-point covariance function is an eighth-rank tensor Rkℓmnprst(x,y):=⟨[Ckℓmn(x)−⟨Ckℓmn(x)⟩][Cprst(y)−⟨Cprst(x)⟩]⟩.2.8Taking the TRF C_klmn to be WSS and wide-sense isotropic, we have Rkℓmnprst(h). Using group representation theory, Malyarenko & Ostoja-Starzewski [37] have recently proved that this covariance is given, most generally, in terms of 29 continuous scalar functions K_n, along with 29 corresponding tensor functions Lkℓmnprstn, Rkℓmnprst(h)=∑n=129Lkℓmnprstn(h)Kn(∥h∥).2.9With such an explicit representation of a TRF of second or fourth rank established, one can turn, say, to simulations. Substantial work on specific cases of TRFs has been reported in [17,18] and references therein. In the following, we work with random viscoelastic media such as delineated here.

Given that random fields are sets of realizations over space domains, it is also important to be able to determine the probabilistic characteristics (mean, correlation function, moments, distributions, etc.) of a random field in terms of just one realization C(ω) of the relaxation modulus. Roughly speaking, a WSS TRF is said to be ergodic in the mean if the spatial average of C(ω,x,t) over B_δ converges to the expected value when the volume |B_δ| tends to infinity, lim|Bδ|→∞1|Bδ|∫BδCkℓmn(ω,x,t) dV≡Ckℓmn(ω,t)¯=⟨Ckℓmn(x,t)⟩≡∫ΩCkℓmn(ω,x,t) dP(ω).2.10In applications, C(ω,t)¯ is computed from a finite number of sampling points N (taken over one realization ω), while ⟨C(x,t)⟩ is computed from a finite number M of available realizations ω (taken at a chosen sampling point x).

(b) Random chessboard model

In keeping with the previous work on the viscoelasticity of random materials in the time domain [9], the computational studies reported further below will focus on a so-called random chessboard (or checkerboard) in two dimensions; see figure 1b for a sample realization at mesoscale δ=40 and a nominal volume fraction of 50% of either phase. This microstructure is easy to simulate with square-shaped finite elements, yet it presents an interesting material system at that volume fraction where both inclusions and quite large clusters of grains are encountered (the percolation occurs at approx. 59%). Technically speaking, the random chessboard comes from a Bernoulli process generated on a Cartesian lattice, with each square cell occupied, independently of realizations at all other cells, with probability p and 1−p by phases a and b, respectively. Henceforth, we adopt p=12. For a square L×L lattice, the number of different realizations (or the size of the entire sample space Ω) is |Ω|=2^L×L. Each elementary event ω gives rise to one realization B_δ(ω), which occurs with probability 1/2^L×L. The Bernoulli process is strict white noise in two dimensions, i.e. without any spatial correlation between adjacent points. However, our methodology also applies to microstructures with spatial correlations such as that shown in figure 1c [33].

The mechanical properties of either phase are now specified through its relaxation modulus tensor C_klmn(t) (or, again equivalently, its creep compliance tensor S_klmn(t)). The random medium is now given in terms of a two-state TRF Cklmn:Ω×E2×T→{Cklmn(a),Cklmn(b)}or{Cklmn(ω,x,t); ω∈Ω,x∈B⊂E2,t∈T},2.11where the superscripts (a) and (b) denote the two possible component phases. This is a special case of (2.3). The correlation function of this Bernoulli-type random chessboard is essentially a strict white noise type, a special case of the general model introduced above.

(c) Complex moduli of linear viscoelastic solids

Under steady-state harmonic oscillation, the constitutive equations (2.4) can be simplified to a quasi-elastic type, e.g. [38]. Thus, if a linear viscoelastic material is subjected to a time-harmonic strain function with circular frequency γ, the steady-state response in stress is also harmonic with the same γ, and vice versa, ϵkl(t)=ϵklA eiγtandσkl(t)=σklA eiγt,2.12where ϵklA, σklA are the strain and stress amplitudes; i=−1. Upon substituting (2.12) into the convolution integral and eliminating e^iγt on both sides, the quasi-elastic relationships for viscoelasticity under steady-state oscillation are derived as: σklA=Cklmn∗(γ)ϵmnAandϵklA=Sklmn∗(γ)σmnA,2.13the tensors Cklmn∗ and Sklmn∗ are now called the complex modulus tensor and complex compliance tensor, respectively. In (2.13), all the terms σklA, ϵklA, Cklmn∗ and Sklmn∗ are time independent, but functions of the loading frequency (γ), and they can all be complex numbers.

The physical meaning of complex moduli is better illustrated by expressing equation (2.13) in polar coordinates, |σklA| eiϕσ=|Cklmn∗(γ)| eiϕC|ϵmnA| eiϕϵand|ϵklA| eiϕϵ=|Sklmn∗(γ)| eiϕS|σmnA| eiϕσ,2.14where |.| denotes the magnitude of the complex number and ϕ represents the rotation angle with respect to the polar axis. Thus, one can get two groups of separate equations: |Cklmn∗(γ)|=|σklA||ϵmnA|,|Sklmn∗(γ)|=|ϵklA||σmnA|2.15and ϕC=ϕσ−ϕϵ,ϕS=ϕϵ−ϕσ.2.16

Clearly, the absolute value of complex moduli characterizes the ratio between magnitudes of stress and strain, while its associated argument ϕ, also called a phase angle in viscoelasticity, represents the time lag between stress and strain. In most viscoelastic solids, the stress and strain signals are not in-phase, which makes ϕ≠0. The above equations also justify the experimental approaches to measure complex moduli in dynamic testing. A comparison of the equations in (2.15) and (2.16) leads to relationships for interchange of viscoelastic properties in the frequency domain, which are much more straightforward than time domain properties, |Sklmn∗(γ)|=1|Cklmn∗(γ)|,ϕC=−ϕS.2.17Technically speaking, the complex modulus and complex compliance are just the Fourier transforms of the relaxation modulus and creep compliance in the time domain. Hence, once either one of these is determined, the other should automatically follow.

Assuming the viscoelastic response on the microscale (i.e. of each phase) to be isotropic, the time-dependent stiffness tensor C(t) is characterized by the relaxation shear modulus μ(t) and relaxation bulk modulus k(t) C(t)=2μ(t)K+3k(t)J,2.18where J and K are, respectively, the spherical and deviatoric parts of the unit fourth-order tensor J=13δklδmn;K=Iklmn−13δklδmn;Iklmn=12(δkmδln+δknδlm),2.19and μ(t) and k(t) are each modelled via Prony series μ(t)=μ0[1−∑n=1Ng¯n(1−e−t/τn)]andk(t)=k0[1−∑m=1Mk¯m(1−e−t/τm)],2.20with g¯n,k¯m, τ_i being the parameters fitted through experimental tests. In general, the time dependence of shear and bulk is not necessarily the same. On the other hand, once all the coefficients in the Prony series are determined, the frequency domain properties are also set. Taking, for example, the shear modulus μ(t), its storage and loss moduli are Re[μ∗(γ)]=1−∑n=1Ng¯n+∑n=1Ng¯nτn2γ21+τn2γ2andIm[μ∗(γ)]=∑n=1Ng¯nτnγ1+τn2γ2,2.21with the magnitude and phase angle of complex properties being |μ∗(γ)|=Im[μ∗]2+Re[μ∗]2andtan⁡(ϕμ∗)=Im[μ∗]Re[μ∗].2.22The same type of formula for the bulk modulus k* holds as well.

(a) Apparent properties of a heterogeneous body

For a heterogeneous viscoelastic body at finite scales, B_δ(ω), the constitutive equations (2.13) in the homogenized material systems are generally defined in terms of volume averages (with an overbar), σ¯klA=Cklmn,δ∗(γ,ω)ϵ¯mnAandϵ¯klA=Sklmn,δ∗(γ,ω)σ¯mnA.3.1Note that the expression σ¯klA denotes the amplitude of the volume-averaged stress, which is evaluated by taking the volume average first and then deriving the amplitude. The subscript δ in Cklmn∗ and Sklmn∗ demonstrates the scale dependence of established properties, to be further discussed below. Generally, σ¯klA and ϵ¯klA are functions of the given mesoscale, loading frequency γ and realization ω∈Ω. Here, according to Huet [6], the tensors Cklmn,δ∗(γ,ω) and Sklmn,δ∗(γ,ω) may also be called apparent complex modulusand compliance, respectively. We can determine such properties from the standpoint of the Hill–Mandel macrohomogeneity condition, which guarantees the equivalence between energetic and mechanical approaches in setting up constitutive equations. In the time-dependent material-like viscoelasticity, this condition [6] is σ:ϵ˙¯=σ¯:ϵ˙¯,3.2where the overdot represents the time derivative. At finite scales, (3.2) holds under either of three types of boundary loading applied to a specific realization Bδ(ω)∈Bδ through:

• — kinematic uniform boundary condition (KUBC) uk(x,t)=ϵkl0(t)xl,∀x∈∂Bδ(ω);3.3

• — static uniform boundary condition (SUBC) tk(x,t)=σkl0(t)nl,∀x∈∂Bδ(ω);3.4

• — kinematic-static (or mixed-orthogonal) boundary condition (uk(x,t)−ϵkl0(t)xl)(tk(x,t)−σkl0(t)nl)=0,∀x∈∂Bδ(ω).3.5

Here σkl0(t) and ϵkl0(t) are the imposed constant tensor functions of time (but not space). (Actually, in viscoelasticity considered here, the loading (3.4) is quasi-static, but we call it ‘static’ consistent with the earlier works of Huet.)

Equations (3.3) and (3.4) also justify the mean (volume average) strain and stress theorems: ϵkl(t)¯=ϵkl0(t)andσkl(t)¯=σkl0(t).3.6

Given that Bδ is a statistical ensemble of realizations of B_δ(ω) for any mesoscale δ, the boundary conditions (3.3)–(3.5) specify, respectively, three types of stochastic initial-boundary value problems (BVPs). Henceforth, we focus on the first two of these as they can deliver bounds (indeed, scale-dependent bounds) on the macroscopic response. More to the point, the mesoscale (apparent) properties of a heterogeneous body, i.e. B_δ(ω), are generally scale and boundary condition dependent. In continuum elasticity, bounds of elastic responses of random composites have been proved using variational principles in combination with assumptions of spatial ergodicity (or somewhat more restrictive mixing) and stationarity of the microstructure [3–5]. For viscoelastic random materials, similar bounds for the relaxation modulus and creep compliance in the time domain have been derived using extended viscoelastic minimum principles [6–8] and also demonstrated numerically for a random chessboard microstructure [9].

In the time-harmonic case σkl(x,t)=σ^kl(x) eiγt and ϵkl(x,t)=ϵ^kl(x) eiγt; the hat being used to indicate the spatially dependent quantities. Thus, working in the frequency domain, the linear viscoelastic problem is equivalent to a spatial problem (neglecting the inertial effects and body forces) σ^kl,l=0,in B,u^k=g^k,on ∂Bδuandσ^klnl=h^k,on ∂Bδt}3.7with the constitutive equations and strain–displacement relations also satisfying at all locations: σ^kl=Cklmn∗ϵ^mn3.8and ϵ^kl=12(u^k,l+u^k,l).3.9In (3.7), ∂Bδu and ∂Bδt are, respectively, the displacement-controlled and traction-controlled parts of the entire boundary ∂B_δ, such that ∂Bδu∪∂Bδt=∂B and ∂Bδu∩∂Bδt=∅. The boundary conditions are also assumed to have harmonic dependence: uk(x,t)=gk∗exp⁡(iγt) on ∂Bδu and tk(x,t)=hk∗exp⁡(iγt) on ∂Bδt. From (3.7) to (3.9), it is clear that the equations in viscoelastic steady-state oscillation retain the same form as those of ordinary elasticity except that the stiffness tensor is now complex.

(b) Hierarchies of mesoscale bounds

Focusing on the steady-state response, the Hill–Mandel condition (3.2) becomes σ^:ϵ^¯=σ^¯:ϵ^¯,3.10with (3.6) replaced by ϵ^¯kl=ϵkl0andσ^¯kl=σkl0.3.11The apparent properties are defined from a modification of (3.3) and (3.4):

• — uk(x,t)=ϵkl0sin⁡(γt)xl ∀x∈∂Bδ(ω) and σ¯klA=Cklmn,δ∗(γ,ω)ϵmn0;3.12

• — tk(x,t)=σkl0sin⁡(γt)nl ∀x∈∂Bδ(ω) and ϵ¯klA=Sklmn,δ∗(γ,ω)σmn0.3.13

In the above ϵkl0 and σkl0 are real-valued constant tensors, while σ¯klA and ϵ¯klA are complex-valued amplitudes of the volume average stress and strain. To be more specific, if the volume average stress is fitted as σ¯kl(t)=σkl0sin⁡(γt+ϕ), then σ¯klA=σkl0 eiϕ.

The derivation of scale-dependent bounds follows a similar procedure to that for elastic random materials and is illustrated in two dimensions with the help of figure 2; the same approach then holds for a three-dimensional situation. To this end, consider a square-shaped mesodomain of a composite body B_δ(ω) and evenly partition it into four square-shaped bodies B_δs(ω), s=1,…,4, with δ_s=δ/2. Next, introduce two types of BVPs in terms of KUBC applied over the mesodomain B_δ(ω): unrestricted uk(x)=ϵ^kl0xl,∀x∈∂Bδ(ω),3.14and restricted (^r) ukr(x)=ϵ^kl0xl,∀x∈∂Bδs(ω), s=1,…,4.3.15Given that the stress and strain fields of the BVP under the restricted condition (3.15) are admissible under the unrestricted condition (3.14), whereas the stress and strain fields of the BVP under (3.14) are inadmissible under (3.15), we find the inequality U(ω,ϵ^kl0)≤Ur(ω,ϵ^kl0),3.16where U(ω,ϵ^kl0) and Ur(ω,ϵ^kl0) stand for energies stored in B_δ(ω) under (3.14) and (3.15), respectively. Upon ensemble averaging, we obtain ⟨U(ϵ^kl0)⟩≤⟨Ur(ϵ^kl0)⟩.3.17

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Partition of the domain into four subdomains, each of mesoscale δ^′=δ/2.

Being statistically averaged, the two media compared in (3.17) are homogeneous, so that now the third minimum theorem of Christensen [38] applies. As a result, and in view of the statistical isotropy of the microstructure, these statistically averaged media are isotropic so that we obtain inequalities on apparent relaxation shear (μδ∗) and bulk (kδ∗) moduli Re(⟨μδ∗(γ)⟩)≤Re (⟨μδs∗(γ)⟩),Im(⟨μδ∗(γ)⟩)≤Im(⟨μδs∗(γ)⟩), ∀δs=δ23.18and Re(⟨kδ∗(γ)⟩)≤Re(⟨kδs∗(γ)⟩),Im(⟨kδ∗(γ)⟩)≤Im(⟨kδs∗(γ)⟩), ∀δs=δ2.3.19Note that it is preferable to work with these two moduli rather than with the shear modulus and Poisson ratio, because the Poisson ratio can be regarded as frequency independent only in very specific situations [39,40].

Just as for the elastic heterogeneous media, the conclusions in (3.18) and (3.19) can also be applied to non-commensurate partitions, leading to scale-dependent hierarchies for complex shear and bulk responses, ∀δ^′<δ, Re(μ∞∗(γ))≤⋯≤Re(⟨μδ∗(γ)⟩)≤Re(⟨μδ′∗(γ)⟩)≤⋯≤Re(⟨μ1∗(γ)⟩)andIm(μ∞∗(γ))≤⋯≤Im(⟨μδ∗(γ)⟩)≤Im(⟨μδ′∗(γ)⟩)≤⋯≤Im(⟨μ1∗(γ)⟩)}3.20Re(k∞∗(γ))≤⋯≤Re(⟨kδ∗(γ)⟩)≤Re(⟨kδ′∗(γ)⟩)≤⋯≤Re(⟨k1∗(γ)⟩)andIm(k∞∗(γ))≤⋯≤Im(⟨kδ∗(γ)⟩)≤Im(⟨kδ′∗(γ)⟩)≤⋯≤Im(⟨k1∗(γ)⟩).}3.21The tensors μ∞∗ and k∞∗ denote the responses for an infinitely large domain (i.e. macroscale or RVE).

The foregoing approach may be employed to assess scaling trends of compliances under static boundary conditions. One needs to introduce unrestricted SUBC tk(x)=σkl0nl(x),∀x∈∂Bδ,3.22and restricted (^r) SUBC tkr(x)=σkl0nl(x),∀x∈∂Bδs(ω), s=1,…,4,3.23so as, using the principle of minimum complementary energy, followed by statistical averaging, to obtain an inequality ⟨U∗(σ^kl0)⟩≤⟨U∗r(σ^kl0)⟩.3.24Here the terms on the left- and right-hand sides stand, respectively, for energies stored in B_δ under (3.22) and (3.23). The inequality hinges on the fact that the stress and strain fields of the BVP under the restricted condition (3.23) are admissible under the unrestricted condition (3.22), but not vice versa.

Again, since we now deal with two homogeneous and isotropic media, the fourth minimum theorem of Christensen [38] applies. As a result, we obtain inequalities on apparent shear (Jδ∗) and bulk (Lδ∗) compliances Re(J∞∗(γ))≤⋯≤Re(⟨Jδ∗(γ)⟩)≤Re(⟨Jδ′∗(γ)⟩)≤⋯≤Re(⟨J1∗(γ)⟩)andIm(J∞∗(γ))≥⋯≥Im(⟨Jδ∗(γ)⟩)≥Im(⟨Jδ′∗(γ)⟩)≥⋯≥Im(⟨J1∗(γ)⟩)}3.25Re(L∞∗(γ))≤⋯≤Re(⟨Lδ∗(γ)⟩)≤Re(⟨Lδ′∗(γ)⟩)≤⋯≤Re(⟨L1∗(γ)⟩)andIm(L∞∗(γ))≥⋯≥Im(⟨Lδ∗(γ)⟩)≥Im(⟨Lδ′∗(γ)⟩)≥⋯≥Im(⟨L1∗(γ)⟩).}3.26

Since all the terms in (3.25)₂ and (3.26)₂ are negative, the inequality in its absolute values is reversed in order. On account of (3.20) and (3.25), there results a hierarchy of absolute values of complex shear modulus ∣⟨J1∗⟩∣−1≤⋯≤∣⟨Jδ′∗⟩∣−1≤∣⟨Jδ∗⟩∣−1≤⋯≤∣⟨J∞∗⟩∣−1=∣⟨μ∞∗⟩∣≤⋯≤∣⟨μδ∗⟩∣≤∣⟨μδ′∗⟩∣≤⋯≤∣⟨μ1∗⟩∣, ∀δ′<δ,3.27while (3.21) and (3.26) lead to a hierarchy of absolute values on complex bulk modulus ∣⟨L1∗⟩∣−1≤⋯≤∣⟨Lδ′∗⟩∣−1≤∣⟨Lδ∗⟩∣−1≤⋯≤∣⟨L∞∗⟩∣−1=∣⟨k∞∗⟩∣≤⋯≤∣⟨kδ∗⟩∣≤∣⟨kδ′∗⟩∣≤⋯≤∣⟨k1∗⟩∣, ∀δ′<δ.3.28

As mentioned earlier, the mesoscale δ=1 essentially represents a domain that contains only one phase in a random chessboard, and this is where elementary bounds can be taken. Considering, for instance, the complex shear modulus, the ensemble averages ⟨μ1∗⟩ and ⟨J1∗⟩ are the weighted properties based on the relative statistical occurrences of both phases. In micromechanics ⟨μ1∗⟩ and ⟨J1∗⟩ correspond to the Voigt and Reuss bounds (e.g. [3]), which, for two-phase composites, take the forms Voigt bound: \ ⟨μ1∗⟩=vμa∗+(1−v)μb∗3.29and Reuss bound: \ ⟨J1∗⟩=vJa∗+(1−v)Jb∗=vμa∗+1−vμb∗,3.30with v being the volume fraction of the phase a and a, b denoting both component phases.

(c) Scaling function

With reference to Ranganathan & Ostoja-Starzewski [29], the scaling function quantifies the discrepancies between ensemble-averaged properties at any given mesoscale, showing that the smaller the discrepancy, the closer the material is to the RVE. This is now generalized to a viscoelastic random material setting. Begin with the constitutive equation in terms of the convolution integrals for isotropic viscoelasticity: σkl(t)=3∫0tk(t−τ)(13ϵ˙pp(τ)δkl)dτ+2∫0tμ(t−τ)(ϵ˙kl−13ϵ˙ppδkl)dτ,3.31which, in terms of complex moduli, reads σklA=3k∗(γ)(13ϵppAδkl)+2μ∗(γ)(ϵklA−13ϵppAδkl).3.32Keeping in mind that, under the boundary conditions (3.3) and (3.4), all the stress components are in phase and all the strain components are also in phase (albeit with a phase shift relative to stresses), the viscoelastic material on any mesoscale is characterized uniquely by two pairs of complex moduli: {μδ∗(γ), kδ∗(γ)} and {Jδ∗(γ),Lδ∗(γ)}. Thus, the ensemble-averaged complex modulus ⟨Cδ∗⟩ and compliance ⟨Sδ∗⟩ are ⟨Cδ∗⟩=2⟨μδ∗(γ)⟩K+3⟨kδ∗(γ)⟩Jand⟨Sδ∗⟩=12⟨Jδ∗(γ)⟩K+13⟨Lδ∗(γ)⟩J,3.33where J and K have been defined in (2.19). A contraction of ⟨Cδ∗⟩ and ⟨Sδ∗⟩ in (3.33) leads to a scaling function f* f∗(δ)=5⟨μδ∗(γ)⟩⟨Jδ∗(γ)⟩+⟨kδ∗(γ)⟩⟨Lδ∗(γ)⟩−6,3.34which effectively generalizes the scaling function of random elastic materials [29]. For random viscoelastic materials, the four moduli appearing in (3.34) are frequency dependent and also complex, rendering the scaling function f*(δ) complex rather than real.

Focusing on the special case of the plane stress state (σ₃₁=σ₃₂=σ₃₃=0), generalizing from the linear elastic case [32], the scaling function for viscoelasticity becomes f∗(δ)=2⟨μ2D∗,δ(γ)⟩⟨J 2D∗,δ(γ)⟩+⟨k2D∗,δ(γ)⟩⟨L2D∗,δ(γ)⟩−3,3.35where μ2D∗=σ¯12A/(2∗ϵ¯12A) and k2D∗=(σ11+σ22¯)A/(2(ϵ11+ϵ22¯)A). In polar coordinates, the scaling function is expressed as f∗(δ)=2|⟨μ2D∗,δ(γ)⟩||⟨J 2D∗,δ(γ)⟩|ei(ϕμ+ϕJ)+|⟨k 2D∗,δ(γ)⟩| |⟨L 2D∗,δ(γ)⟩| ei(ϕk+ϕL)−3.3.36

The general algorithm for determining the scaling function is given below. For demonstration, we take the example of general isotropic viscoelasticity, for which the scaling function follows (3.34). Only one loading frequency γ is presented in the algorithm, and for other frequencies the same procedure has to be repeated.

(a) Computational mechanics procedure

In order to quantitatively demonstrate the mesoscale hierarchies and the scaling function as well as their frequency dependencies, we run computational mechanics on the planar random chessboard of figure 1b. The BVPs are solved using the finite-element method (FEM) with a commercial solver Abaqus. For each realization and any given frequency, the two-dimensional bulk modulus k2D,δ∗ and two-dimensional shear modulus μ2D,δ∗ are calculated separately using (3.12) and (3.13). For the shear behaviour, the prescribed strain and stress tensors are set to be ϵ110=ϵ220=0, ϵ120≠0 for KUBC and σ110=σ220=0, σ120≠0 for SUBC. The mesoscale shear moduli are μ2D,δ∗=σ¯12Aϵ120andJ2D,δ∗=ϵ¯12Aσ120.4.1

Similarly, we have ϵ110=ϵ220=ϵ0≠0, ϵ120=ϵ210=0 in KUBC and σ110=σ220=σ0≠0, σ120=σ210=0 in SUBC for bulk behaviour. The mesoscale bulk moduli follow from k2D,δ∗=(σ11+σ22¯A)4ϵ0andL2D,δ∗=2(ϵ11+ϵ22¯A)2σ0.4.2

Simulations are performed over a set of realizations at each mesoscale and also over a range of loading frequencies. These realizations, mesoscales and frequencies are specified below. Upon ensemble averaging, the frequency-dependent apparent properties at mesoscale δ are derived: from 4.1:  ⟨μ2D,δ∗(γ)⟩;⟨J2D,δ∗(γ)⟩4.3and from 4.2:  ⟨k2D,δ∗(γ)⟩;⟨L 2D,δ∗(γ)⟩.4.4

In order to facilitate a comparison, the material properties are kept the same as in the time domain investigation [9], where one phase is taken as linear viscoelastic while the second one is elastic. The viscoelastic behaviour is implemented using the Prony series as (2.20) with only one term considered for simplicity, which, in terms of mechanical analogues, corresponds to a generalized Maxwell (or Zener) model. Values of all the relevant coefficients are listed in table 1, with the elastic phase being twice as stiff in Young’s modulus but having the same Poisson ratio as the viscoelastic phase. The bulk and shear behaviours of the viscoelastic phase are set to be different, with the bulk modulus less sensitive to time (k¯1<g¯1), which is consistent with many results on polymers [41].

In the finite-element simulations, all the cells of the chessboard are assumed to be perfectly bonded and simulations are performed at δ=2,4,8,16 over the frequency range from 0.05 to 50 Hz. At each mesoscale, both stochastic BVPs are handled by a Monte Carlo sampling of the random medium Bδ, with the number of realizations generated chosen such that the ensemble-averaged behaviour stabilizes to within a small error. For the moderate material contrast taken here, the numbers of realizations for δ=4,8,16 are, respectively, 100, 60, 40. The Monte Carlo approach is not necessary for δ=2, where the entire ensemble of 16 realizations of B2 can fully be accounted for. Each BVP is solved quasi-statically in Abaqus (using the key word *VISCO) and only the steady-state stress and strain outputs are used to determine mesoscale properties for that realization. The stress and strain data in the domain are first averaged and then fitted to a harmonic function to extract amplitude and phase angle.

To better resolve the stress and strain distributions at the interface of two phases, quadratic elements are used instead of bilinear elements; the mesh density becomes more important in ensuring the accuracy of the solution obtained using FEM in the case of strong mismatch in the phases’ properties [33]. For our microstructure of a random chessboard at volume fraction 0.5, a comprehensive mesh sensitivity analysis was conducted in Zhang & Ostoja-Starzewski [9], where it was found that, for the chosen constituent properties, the mesh density of 4×4 elements per one square of the chessboard provides a stable solution and sufficient accuracy.

(b) Hierarchies of mesoscale bounds

Upon Monte Carlo sampling, the averaged mesoscale properties of two-dimensional shear modulus and two-dimensional bulk modulus for all the mesoscales and both boundary conditions are plotted in figures 3 and 4. The frequency dependence of the viscoelastic phase is determined analytically through (2.21) and (2.22) and the Voigt and Reuss bounds are also given by formulae (3.29) and (3.30).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Hierarchies of mesoscale bounds on the two-dimensional complex shear modulus: (a) absolute value and (b) phase angle.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Hierarchies of mesoscale bounds on the two-dimensional complex bulk modulus: (a) absolute value and (b) phase angle.

As shown in figures 3 and 4, the ensemble mesoscale properties, in terms of both magnitude and phase angle, exhibit a clear frequency dependence, which means that the effect of viscoelasticity has been blended into the composites. The numerically derived mesoscale properties fall between those of the component viscoelastic (VE) phase and elastic (E) phase, and also the traditional Voigt bound and Reuss bound. There is a consistent trend for both (KUBC and SUBC based) ensemble mesoscale properties to converge to one another, from above and from below, as the mesoscale δ increases, with ever lower discrepancy, thus providing converging bounds on the effective (RVE) properties of the composite. For the particular system, the mesoscale bounds at δ=16 provide very close estimates for both bulk and shear moduli of the random chessboard at volume fraction 0.5.

The hierarchies on the magnitude of complex shear modulus |〈μ*〉| in three dimensions carry over to two dimensions. The two-dimensional bulk modulus is related to that in three dimensions through 1/k_2D=1/3μ+4/9k [42]. The hierarchies in the real and imaginary parts of the three-dimensional complex modulus k* and μ* are preserved. To sum up, the conclusions in two dimensions are analogous to those in three dimensions. The theoretical hierarchies (3.27) and (3.28) in terms of the magnitude of two-dimensional complex shear and bulk moduli are confirmed in the numerical results of figures 3a and 4a. Additionally, new hierarchies in phase angles are brought out in figures 3b and 4b, which are −⟨ϕJ1∗⟩≤⋯≤−⟨ϕJδ′∗⟩≤⟨−ϕJδ∗⟩≤⋯≤−⟨ϕJ∞∗⟩=⟨ϕμ∞∗⟩≤⋯≤⟨ϕμδ∗⟩≤⟨ϕμδ′∗⟩≤⋯≤⟨ϕμ1∗⟩,∀δ′<δ4.5and −⟨ϕL1∗⟩≤⋯≤−⟨ϕLδ′∗⟩≤⟨−ϕLδ∗⟩≤⋯≤−⟨ϕL∞∗⟩=⟨ϕk∞∗⟩≤⋯≤⟨ϕkδ∗⟩≤⟨ϕkδ′∗⟩≤⋯≤⟨ϕk1∗⟩,∀δ′<δ.4.6

If the mesoscale properties are plotted as a function of mesoscale δ for multiple frequencies, the rate of convergence in mesoscale bounds is found to be dependent on the loading frequency. In fact, as figure 5 shows, the discrepancy between mesoscale bounds on |μ*| is larger for the same δ at lower frequencies. This is because the viscoelastic material is softer than the elastic phase and, at lower frequencies, this introduces larger discrepancies in the material properties. Thus, a larger domain is required to homogenize the random material to RVE. Hence, in homogenization of the viscoelastic composite, the actual size of RVE varies with loading frequency and depends on the actual discrepancies between the constituents at that frequency. For comparison, the mesoscale bounds for the corresponding elastic problem, which contains the same microstructure but with only elastic properties (g₁,k₁,τ₁=0 in table 1), are also included. The elastic problem can be regarded as a limiting case at an extremely high loading frequency in which the loss modulus of (2.21) goes to 0.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Hierarchies of mesoscale bounds for several frequencies γ.

(c) Normalized scaling function

Finally, we would like to determine the form of scaling function for plane stress viscoelasticity (3.35). For the random chessboard, it should be a function of δ, frequency γ, and material properties of both component phases (k1∗,μ1∗ and k2∗,μ2∗), f∗(δ,k1∗,μ1∗,k2∗,μ2∗,γ)=2⟨μ2D,δ∗(γ)⟩⟨J 2D,δ∗(γ)⟩+⟨k2D∗,δ(γ)⟩⟨L 2D,δ∗(γ)⟩−3.4.7

Now, note that the scaling function f*(δ) should go to 0 as δ→∞ since the KUBC and SUBC controlled moduli are identical for the RVE, i.e. limδ→∞Re( f∗)=limδ→∞Im( f∗)=0or| f∗|=0.4.8Also, since at δ=1, the KUBC and SUBC controlled moduli are just the Voigt and Reuss bounds μ2D∗V(γ):=⟨μ2D,1∗(γ)⟩,k2D∗V(γ):=⟨k 2D,1∗(γ)⟩andJ2D∗R(γ):=⟨J 2D,1∗(γ)⟩,L 2D∗R(γ):=⟨J 2D,1∗(γ)⟩,}4.9we have f∗(δ=1,γ)=2μ2D∗V(γ)J2D∗R(γ)+k2D∗V(γ)L 2D∗R(γ)−3.4.10Since f*(δ,γ) attains a maximum at δ=1 and then monotonically decreases as δ increases, we introduce a normalized scaling function g*(δ,γ): g∗(δ,γ)=f∗(δ,γ)| f∗(1,γ)|,4.11such that 0≤|g*(δ,γ)|≤1, |g*(1,γ)|=1 and |g∗(∞,γ)|=0.

As discussed earlier, the viscoelastic material exhibits different properties at different frequencies. Hence, the simulation of a composite at any different frequency can be regarded as a combination of different component phases. In this analysis, we consider seven loading frequencies spread over a wide range from 0.05 to 50 Hz. The analytic values of the modulus and phase angle in two dimensions at those frequencies are listed in table 2. Since the elastic problem is also a special case of viscoelasticity and in order to test the applicability range of g*(δ,γ), we add one more scenario where the composite is made of two elastic phases at large contrast (E₁/E₂=8), i.e. mat A and mat C (E=240 MPa, ν=0.3) in table 2.