Homogenized out-of-plane shear response of three-scale ﬁber-reinforced composites

In the present work we embrace a three scales asymptotic homogenization approach to investigate the effective behavior of hierarchical linear elastic composites reinforced by cylindrical, uniaxially aligned ﬁbers and possessing a periodic structure at each hierarchical level of organization. We present our novel results assuming isotropy of the constituents and focusing on the effective out-of-plane shear modulus, which is computed exploiting the solution of the arising anti-plane problems. The latter are solved semi-analytically by means of complex variables and successfully benchmarked against the results obtained by ﬁnite elements. Our ﬁndings can pave the way for multiscale modeling of complex hierarchical materials (such as bone and tendons) at a negligible computational cost.


Introduction
Multiscale composites organized across two or more length scales are often encountered in nature, as well as in artificial materials designed to optimize specific properties (see, e.g., [10]). Relevant applications involving hierarchical systems include, but are not limited to, rocks and fracture [11], biomechanics and nanomedicine [28], poroelasticity [21], and the bone tissue [26].
Multiscale modeling via suitable homogenization techniques is of crucial importance in the analysis of composite materials. On one hand, it can be exploited to obtain information concerning the coarser scale behavior of a system at hand (where experimental measurements are often easier to provide) on the basis of given microstructural properties of the finer scales. On the other hand, it can provide a deeper understanding on how to modify the microstructural arrangement of the finer scales constituents to achieve the optimal design of artificial constructs (e.g. biomimetic materials).
The most widely exploited techniques dealing with this issue in terms of mechanical properties rely on either the aver-(a) (b) (c) Fig. 1 Cross sections of the periodic composite: a whole material, b ε 1 -structural level. c ε 2 -structural level age field techniques, such as Eshelby-based approaches or representative volume elements formulations, as well as the asymptotic homogenization technique, see, e.g., the review [9] for a comparison. The former have been much more extensively applied to complex, hierarchical systems, although also the computational potential and theoretical reliability of the asymptotic homogenization technique has been recently investigated (see, e.g., [4,19,20]). In particular, the reiterated homogenization technique [1,3,12,30] has been considered in the theoretical study of heterogeneous composites presenting multiple length scales. From the application point of view, an extensively discussed hierarchical hard tissue has been bone and several hierarchical modeling via homogenization and continuum micromechanics methods have been used in the analysis of its mechanical properties [7,14,29].
In the present work we apply the three-scale asymptotic homogenization developed in [25] (which is therein applied to hierarchical layered materials) to investigate the effective behavior of fibrous hierarchical composites. The linear elastic constitutive framework and the homogenization steps developed in [25] are briefly illustrated in Sects. 2 and 3, respectively. We apply the latter to find the effective out-of-plane shear of a hierarchical linear elastic reinforced composite assuming a square-like arrangement of uniaxially aligned cylindrical fibers at all hierarchical levels of organization. We solve the resulting anti-plane cell problems (which are depicted in Sect. 4) for isotropic and piecewise constant materials properties using results of our previous works (see e.g. [5,6]). Our novel results (which match those computed numerically via finite elements) are presented and discussed in Sect. 5 in terms of the out-of-plane shear versus the fibers' volume fraction. We conclude the manuscript (cf. Sect. 6) highlighting possible applications to realistic scenarios of interest, such as hierarchical modeling of fibril bundles and/or fusing mineral crystals as analyzed for example in [22] in the context of aged bone modeling.

The linear elastic problem
Let us denote by ⊂ R 3 a periodic composite possessing two hierarchical levels of organization (Fig. 1). The domain at the first structural level (referred to as ε 1 ) is assumed to be reinforced by uniaxially aligned fibers, that is = 2 denotes the ensemble of fibers and ε 1 1 is the host medium at the ε 1 level. The interface between ε 1 1 and ε 1 2 is denoted by ε 1 . At the second structural level (referred to as ε 2 ), each fiber i ε 1 2 is supposed to be as well reinforced by aligned fibers oriented in the same direction of the composite fiber i ε 1 2 . Then, i 2 represents the ensemble of fibers and ε 2 1 is the host medium at the ε 2 -structural level. The interface between ε 2 1 and ε 2 2 is denoted by ε 2 . At this stage, we assume that all constituents of the hierarchical composite behave as linear elastic materials with constitutive relationship for the stress tensor given by, where ξ (u) is the elastic strain tensor and u(x) is the elastic displacement at x = (x 1 , x 2 , x 3 ) with components u 1 (x), u 2 (x) and u 3 (x), being the components of u(x) in an orthonormal cartesian vector basis. The fourth rank tensor C with components C i jkl (i, j, k, l = 1, 2, 3) is the stiffness tensor, which is assumed to be phase-wise smooth and positive definite and satisfies the standard symmetries, i.e., componentwise respectively, where ω is a second order symmetric tensor and > 0 is a constant. Then, ignoring inertia and volume forces, the problem in reads where n is the outward unit vector normal to the surface ∂ and u * and S * are the prescribed displacement and traction on ∂ = ∂ d ∪ ∂ n (with ∂ d ∩ ∂ n = ∅), respectively. Moreover, continuity conditions for displacements and tractions are imposed on both ε 1 and ε 2 , i.e.
where n η = (n 3 ) represent the outward unit vectors to the surfaces ε 1 and ε 2 , respectively. The operator • denotes the jump across the interface between two constituents.

Three scales asymptotic homogenization
We consider three different scales, namely d 1 , d 2 and L, which characterize the different structural sizes and we assume that they are well-separated, i.e.
Using relation (3), two formally independent variables are introduced. Additionally, each field and material property is considered to be η and ζ periodic [15]. As a consequence of the performed spatial scale decoupling (4) and using the chain rule, it holds that where ∇ j indicates that the derivative is performed with respect to j = x, η, ζ . We now assume that the elastic displacement u ε can be represented as a power series in terms of the small parameters ε 1 and ε 2 , namely whereũ (0) is defined as Since the quantities involved vary on the η and ζ scales, the following cell average operators are defined where |Y | and |Z | represent the periodic cell volumes. Here and subsequently (unless necessary), the variable dependence is dropped out for convenience.

Homogenization technique
The homogenization technique, as depicted in [25], is sketched below. First, we substitute expansion (6) into (1) and (2), and then, use the chain rule and equate in powers of ε 2 "freezing" the small parameter ε 1 . This allows to find the effective elastic properties at the ε 2 -structural level and use the results as the inputs for the problems arising at the ε 1 -structural level when equating in powers of ε 1 . Finally, the effective coefficients of the problem are obtained. The procedure is detailed as follows: Step 1 Substituting expansion (6) into (1); using result (5) and considering terms in powers of ε 2 . with where g is a vector. Since the right hand side of Eq. (8) is zero, the solvability condition is satisfied [2]. Then, i.e. the homogeneity of (8) together with (9)-(10), leads to a periodic ζ -constant solution.
(ii) To O(ε 2 ) and using the fact that ξ ζ kl (ũ 0 ) = 0, By the ζ -periodicity of C ε and the solvability condition, Eq. (11) has a ζ -periodic solution which is unique up to an additive constant. In particular, since the problem is linear andũ 0 does not depend on ζ ,ũ 1 can be written as wherẽ The third rank tensorχ is ηand ζ -periodic, and satisfies the local problem with where G is a third rank tensor. Equations (13)- (15) constitute the ε 2 -cell problem. The condition χ ζ = 0 is imposed in order to guarantee uniqueness. (iii) To O(ε 2 2 ), applying the average operator • ζ and taking into account the ζ periodicity of the involved functions, wherě is the effective stiffness tensor at the ε 1 -structural level.
Again the solvability condition is applied to (18)- (20), to obtain (ii) To O(ε 1 ), using the fact that ξ η kl (u 0 ) = 0 and applying the ζ average operator, By the η-periodicity ofČ ε and the solvability condition, the above equation has a η-periodic solution which is unique up to an additive constant. In particular, since the problem is linear and u 0 does not depend on η, where the third rank tensor χ is η-periodic and solution of Equations (22)-(24) represent the ε 1 -cell problem. The condition χ η = 0 is imposed for guarantee uniqueness. (iii) To O(ε 2 1 ), applying the average operator • η and taking into account the η periodicity of the functions involved. The homogenized problem becomes is the effective stiffness tensor.

Out-of-plane shear mechanical response
Since we are considering uniaxially aligned fibers, the threedimensional cell problems can be equivalently formulated as two-dimensional problems over the cross-section of the cell. That is, by symmetry the auxiliary third rank tensorsχ and χ do not depend on ζ 3 and η 3 , respectively. Therefore, we can refer to the spatial coordinates vectors as η = (η 1 , η 2 ), ζ = (ζ 1 , ζ 2 ) (keeping the same notation for the sake of simplicity), and the third components of the normal vectors n ζ and n η reduce to zero, see also the appendix reported in [19].
We assume that C ε is piecewise constant, such that, the parametric dependence of the ε 2 -cell problem on the variables η and x is lost andχ depends only on ζ . As a consequence,Č ε is also piecewise constant (as it is averaged on ζ ), and therefore χ is depending only on η, so that finally alsoĈ is piecewise constant.
Applying the above mentioned dimensional reduction, and the fact that C ε andČ ε are piecewise constants, the differential problems (13)-(15), (22)-(24) can be rewritten as and respectively, with j, q = 1, 2 and i, k, l, p = 1, 2, 3. In (26) and (27),Ỹ γ andZ γ (γ = 1, 2), denote the two-dimensional unit cells at the ε 1 -and ε 2 -structural levels, respectively, with restrictions to the matrix (γ = 1) and fiber (γ = 2). The interface curves between the constituentsỸ 1 ,Ỹ 2 andZ 1 ,Z 2 are denoted by˜ η and˜ ζ , respectively. Moreover, C γ,η i jkl and C γ,ζ i jkl are the stiffness tensors of the constituent γ at the ε 1and ε 2 -structural levels, respectively. In particular, the stiffness tensor of the composite fiber at the ε 1 structural level is C 2,η i jkl =Č i jkl . We further assume that the constituents at the ε 2 level are isotropic, therefore the fiber phase at the ε 1 -structural level can be at most monoclinic (i.e. 13 independent elastic coefficients), see, e.g. [15,16]. In our case, it is however tetragonal symmetric (6 independent elastic coefficients), as we are dealing with square-symmetric arrangement of cylindrical fibers (the cell's cross section corresponds to a square embedding a circle, see also [15,27]). In particular, the following relationships hold for both structural levels: and, as both C ε andČ ε are in particular orthotropic Here we focus on the shear mechanical response related to the third (out-of-plane) component of the elastic displacement. For a generic monoclinic material, the relevant constitutive relations are [15] where σ 13 and σ 23 are the out-of-plane shear stresses and ξ 13 and ξ 23 are the shear strains. In our case, the resulting effective stiffnessĈ is still expected to be at most monoclinic (as the individual phases are tetragonal and the fibers are uniaxially aligned also on the ε 1 -structural level).
We then aim at obtaining the effective elastic coefficients related to relationships (30) and (31), to characterize the outof-plane shear mechanical response. This can be done by performing a standard decomposition of the problems (26) and (27) into in-plane and antiplane problems. Exploiting (28) and (29), and fixing i, p = 3 in (26) and (27), the resulting nontrivial anti-plane cell problems can be expressed in terms of the doubly periodic functionsχ γ 33α =χ γ α (ζ ) for ζ ∈Z γ and χ γ 33α = χ γ α (η) for η ∈Ỹ γ as follows [15,27] (L ζ α ) and with α = 1, 2. Once solved the cell problem (32), the relevant effective coefficients at the ε 1 structural level are computed by Solving (33), the effective coefficients are given bŷ Since we assume square-symmetry also on the ε 1 structural level, the resulting effective stiffness tensor is actually expected to exhibit tetragonal symmetry, so that in particular as also shown by the results that follow. The effective shear mechanical response is then completely characterized by the effective out-of-plane shear modulusĈ 1313 .

Results and discussion
The theory of analytical functions in [13] can be applied to solve the cell problems (32) and (33), where doubly periodic harmonic functions need to be found (See "Appendix A"). In order to compute the effective properties at both structural levels, it is necessary to truncate the systems (46) and (47) into appropriate orders. Then, whereā 1 (ā 1 ) denotes the conjugate of the complex coefficient a 1 (ã 1 ) and "i" is the imaginary unit. Using formulas (41)-(44) the effective coefficients at both structural levels can be found. For the sake of exemplifying, we choose the following material properties C 1,ζ 1313 = 1, C 2,ζ 1313 = 10 and C 1,η 1313 = 5, and made a parametric study by varying the volume fraction of the fiber at both structural levels. Particularly, the effective coefficients are computed for the same fiber volume fraction at each structural level, i.e. φ = φ ζ = φ η . Then, we proceed in the following way: 1. Given the properties C 1,ζ 1313 and C 2,ζ 1313 and for a fixed volume fraction ofZ 2 , denoted by φ ζ , the cell problem at the ε 2 -structural level (32) is solved. To find the coefficient a l , the infinite linear system (46) is truncated at a certain orderK and solved. The matrix W ζ k is written in terms of φ ζ . 2. To compute the effective coefficientsČ 1313 ,Č 1323 ,Č 2313 andČ 2323 of the composite fiber ε 1 2 , the solutionã l of the linear system is substituted in (41) and (42). 3. Now, the cell problem at the ε 1 -structural level (33) is solved. In particular, we only need to fix C 1,η 1313 since C 2,η 1313 =Č 1313 was found in the previous step. At this point the volume fraction of the fiberỸ 2 must be fixed. One way to proceed is to make φ η = φ ζ . 4. Finally, to find a l , and consequently the effective coeffi-cientsĈ 1313 ,Ĉ 1323 ,Ĉ 2313 andĈ 2323 , the infinite linear system (47) is truncated at a certain order K and solved. In this case, the matrix W η k is equal to W ζ k , which is already computed. Figure 2 presents the values obtained for the effective coefficientsČ 1313 (left) andĈ 1313 (right) as a function of the fiber volume fraction. Specifically, the truncation order of the systems was fixed toK = K = 3, thus, decreasing the computational cost. As a consequence of constituent's isotropic behavior, we obtainČ 1313 =Č 2323 ,Ĉ 1313 =Ĉ 2323 andČ 2313 =Č 1323 =Ĉ 2313 =Ĉ 1323 = 0. A comparison between the results by the three scales asymptotic homogenization method and the finite element method using FreeFEM++ is also shown in Fig. 2. Specifically, we approximate the involved functions by piecewise linear continuous finite elements. As observed, the curves overlap even when the truncation orders of the infinite linear systems are very small. Interestingly, for the particular set of parameters fixed, the effective behavior ofĈ 1313 (in the context of elasticitŷ C 1313 is the effective shear modulusμ) behaves as a quadratic curve in function of the fiber volume fraction. This behavior is a consequence of the fact that at the ε 2 -structural level the fibers are stiffer than the host medium, and at the ε 1 -structural level the matrix becomes the stiffer constituent. In several works similar patterns have been observed. For instance, the enhancement of the effective out-of-plane Young modulus was obtained in [24] for a periodic bilaminate, and in [31], for a particular type of composite material.

Concluding remarks
In the present work, the effective behavior of hierarchical linear elastic composites at different length scales was studied. The approach permits the study of problems where several length scales are present (the two-scale asymptotic homogenization only deals with two different length scales). The novelty of this work relies in the solution, for the first time, of the effective elastic shear stiffness considering a nested arrangement of cylindrical and uniaxially aligned fibers. The corresponding cell problems are solved analytically using complex variables and the results are compared with a finite element approach. Since analytical formulas are found for computing the effective coefficients the computational cost is very low.
Future developments of this work will concern the analysis of the in-plane problems [16,27] to fully characterize the effective mechanical response of hierarchical threescale materials. Moreover, the extension of the present model to interface debonding conditions (see, e.g. [23]), as well as geometric heterogeneities (for example generalizing the non-macroscopically uniform homogenization approach developed in [8,17,18]) will provide more general and realistic results applicable to real-world hierarchical composites. The current method represents a first step towards computationally feasible multiscale modeling of complex hierarchical materials. For example, this approach could be exploited to model musculoskeletal mineralized tissues (such as bone and tendons) which are organized across several spatial scales. In this context, the three scales would represent the arragement on the basic constituents (i.e. collagen and mineral crystals), the resulting mineralized collagen fibril, and finally their packing into the mineralized collagen fibril bundle. Since the model is developed for fiber-reinforced composites, it can then serve as an approximation for both fibers embedded in a matrix, and for fusing inclusions of reinforcing material, such as for example the hydroxyapatite mineral crystals which are found in the aged bone tissue (see, e.g. [22]). equations are obtained in order to find the complex coefficientsã m and a m , where I is the infinite identity matrix, 0 is the infinite zero matrix,Ã 1 = ( 1ã1 , 1ã3 , . . .),Ã 2 = ( 2ã1 , 2ã3 , . . .), a l = 1ãl + 2ãl I (I 2 = −1), A 1 = ( 1 a 1 , 1 a 3 , . . .), A 2 = ( 2 a 1 , 2 a 3 , . . .), a l = 1 a l + 2 a l i, k W j = k W j 1 + k W j 2 i (j = ζ , η) with where R ζ and R η represent the fiber radius ofZ 2 andỸ 2 , respectively. Besides,