Effective bending modulus of thin ply fibre composites with uniform fibre spacing

Cauchy homogenisation methods are based on the assumption that materials show size independent behaviour. For example, the effective bending modulus of a beam is constant for all cross-sectional areas and is equal to the tensile modulus. However, as demonstrated by numerous authors, the stiffness of heterogenous materials under non uniform loading (e.g. bending) can be influenced by size effects, when the microstructural size scale is comparable to the macroscopic one. In other words, an effective elastic modulus can underestimate or overestimate the bending stiffness depending on the specimen size. In this work, the effect of the thickness dimension on the effective bending modulus of composite plies is investigated. The composite is represented by cylindrical fibres uniformly distributed in a uniform matrix. The cases of fibres stiffer and more compliant than the matrix are considered. The results show that the effective bending modulus depends on fibre volume fraction ratio, mismatch of fibre/matrix material properties and the number of fibres through the thickness. The implications for the analysis of bending of thin ply composites are discussed. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Classical homogenisation techniques are based on the assumption that effective material properties are independent of sample size and applied boundary conditions. However, the properties of heterogenous materials can be size dependent under nonuniform loading, if the influence of the microstructural length scale is sufficiently large to affect the overall macroscopic behaviour ( Wheel et al., 2015 ). Recent success in the manufacture of thin ply fibre reinforced polymer composites, i.e. plies as thin as 20 μm in Innovations and Testbeds (2018) , has motivated interest in the investigation of size effects in composite plies. Experimental evidence of size dependent bending behaviour is reported in the literature for materials with fibrous or granular structures, such as cortical bone ( Yang and Lakes, 1982 ), crystalline marble ( Vardoulakis et al., 1998 ) and porcelain ( Ece and e Nakagawa, 2002 ), cellular materials, such as metallic foams ( Brezny and Green, 1990;Andrews et al., 2001 ) and honeycomb panels ( Yasui,20 0 0 ), and composite materials, such as fibrereinforced concrete ( Nguyen et al., 2013 ). Materials with inclu- * Corresponding author. sions more compliant than the matrix were shown experimentally to have increasing stiffness with decreasing size ( Beveridge et al., 2013 ). However, contradictory size effects are also reported in the literature for cortical bone, showing that the size effect may be controlled not just by the stiffness of the inclusions but also by their surface area to volume ratio ( Choi et al., 1990 ).
Generalized continuum theorems can be used to incorporate size dependent material behaviour by including additional degrees of freedom, e.g. micropolar theory, or by considering higher derivatives of displacements, e.g. strain gradient theories, Eringen (1972) , Eringen and Edelen (1972) ; Eringen (2012) . Size effects in heterogenous materials can also be formulated through the use of a Cosserat characteristic length to derive material properties in generalised continua. Bigoni and Drugan (2007) formulated Cosserat effective moduli for composite materials for the case of inclusions more compliant than the matrix. In this case, the effective modulus increases with decreasing size, a phenomena termed 'positive size effect' in Wheel et al. (2015) . However, it is shown in Bigoni and Drugan (2007) that Cosserat effective properties could not be defined for the case of inclusions stiffer than the matrix, significantly limiting its applicability. In Gao and Mahmoud (2014) , Gao (2015) the Cosserat theory was modified to include additional size dependent surface effects to allow for this case, resulting in  ( Gortsas et al., 2018 ) used the boundary element method to determine the length scale parameter to be used in a strain gradient theory model for a composite beam under transverse bending (as opposed to longi-tudinal bending examined in our work). The results showed that, as for the Cosserat theory, the length scale parameter for a gradient elastic material cannot be defined for the case of fibres stiffer than the matrix. In Wheel et al. (2015) , Euler-Bernoulli beam theory was used to formulate the stiffness expression for a composite laminate beam consisting of layers (plies) with different elastic moduli loaded in three point bending. Their study showed that the nature of the size effect (termed size softening or size stiffening) is governed by the relative stiffness of the outer layer-a stiff outer layer yields a positive size effect (stiffening), a more com pliant outer layer gives the opposite result (softening). Furthermore, a characteristic Cosserat length was defined for a micropolar beam, following ( Beveridge et al., 2013 ), which was only applicable for the case of a stiffer outer layer. The work of Frame et al. (2017) extends the work of Wheel et al. (2015) , applying the general approach developed in Wheel et al. (2015) to the case of cortical bone.
Our work follows a similar approach to that in Wheel et al. (2015) . The study in Wheel et al. (2015) examined the size effect for multiple composite plies of different stiffness. Our focus is at the microstructural ply level incorporating fibre/matrix interactions. Fibres as represented as circular cylinders uniformly distributed in a homogenous matrix. Both cases of compliant and stiff fibres, relative to the matrix, are considered. By the use of Euler-Bernoulli beam theory, the effective bending modulus is determined, which is shown to be size dependent. Specifically, the effective bending modulus is a function of fibre volume fraction, mismatch of fibre/matrix modulus and the number of fibres through the thickness (or equivalently the ply thickness). The results have been validated through 3D FE analysis. The effective bending stiffness is used to define explicitly the characteristic length used in the formulation of a micopolar beam stiffness. The results are consistent with the work of Bigoni and Drugan (2007) in which they showed, using strain energy minimisation, that the Cosserat characteristic length for composite materials depends on fibre radius, volume fraction ratio and mismatch of material properties. The novelty of this work is the identification and quantification of the size effect for a composite ply under bending in the fibre direction, taking into account the fibre/matrix geometry and material mismatch at the microstructural level. This is the first paper to specifically address this issue using a combination of analytical and finite element approaches. While this paper focuses at the ply level, we expect that the size effect identified in this work will also affect the bending behaviour of a multi-ply laminate of the type studied in Wheel et al. (2015) .
The paper is laid out as follows: Section 2 derives the equation for the equivalent bending modulus using Euler-Bernoulli beam theory. Section 3 provides FE validation of the results of Section 2 . Section 3 derives the characteristic length for a generalised continuum based on the equivalent bending modulus. Conclusions are provided in Section 4 . Appendix A and Appendix B provide additional detail on the derivation of the equations.

Euler-Bernoulli beam model of thin ply composite with uniform fibre distribution
According to Euler-Bernoulli beam theory, there is a linear relationship between bending moment, M , and curvature, κ, for a beam, where D is the bending stiffness and κ is the second derivative of deflection, v , with respect to position along the beam, x , For a homogeneous cross section with elastic modulus E and second moment of area I with respect to the neutral axis, D = EI.
In this work, we consider a composite ply under pure bending. The bending stiffness for the composite cross-section, D , is determined by summing the product of the elastic modulus of each material (fibre and matrix) with the second moment of area about the neutral axis of the cross section. The effective bending modulus, E b , is then defined as the ratio of D to the homogenised second moment of area about the neutral axis, defined as I . We consider a cross section of N fibre-matrix unit cells repeated through the ply thickness, as in Fig. 1 . Cylindrical fibres are uniformly distributed and arranged symmetrically about the central axis, which is the neutral axis for the bending problem. Two classes of problem, having odd or even number of fibres, N , are shown separately in Fig. 1 (a) and (b), respectively. The plate is under pure bending in the longitudinal (fibre) direction, as illustrated in Fig. 1 (c). The elastic moduli are designated E f and E m for the fibre and matrix, respectively.
In order to calculate the bending rigidity, D , of the cross section, the second moment of area with respect to the neutral axis for each fibre and matrix region should be calculated separately.
The bending stiffness is then, when the number of fibres, N , is odd and when the number of fibres, N , is even. In Eqs. (3) and (4) , I n x | j is the second moment of area of the n th material region, with respect to the neutral axis ( x -axis), for material j ( j = f, m for fibre or matrix, respectively). The parallel axis theorem is used to calculate the second moment of area with respect to the neutral axis, which is parallel to the centroidal axis ( c -axis) of each material region and I c | j is the second moment of area with respect to the c-axis for material j (where j = f, m for fibre or matrix, respectively), with cross sectional area A j for j = f, m, if N odd and 1 ≤ n ≤ N + 1 2 if N even and 1 ≤ n ≤ N 2 (5) The effective bending modulus, E b , is calculated as the ratio of bending stiffness of the heterogenous medium, D , calculated using Eqs. (3) or (4) , to the second moment of area, I , of the equivalent homogenised cross section, where t is the ply thickness = Nd, as in Fig. 1 . Combining Eqs.
(3) to (7) (see Appendix A ) we obtain, where α is the dimensionless geometry parameter, α = r/d, and β is the dimensionless material ratio, β = E f /E m . Note that for the square packing arrangement in Fig. 1 , the maximum value of α is 0.5. In Eq. (8) E is the tensile modulus calculated by the rule of mixtures for a composite, noting that volume fraction, v f = πα 2 .
Eq. (8) shows that the effective bending modulus is given by the effective tensile modulus, which is size independent, plus an additional contribution, which is size dependent, due to the nature of the heterogenous material, i.e. through the dependence on α, β and N . Eq. (8) shows that for a thick ply ( N → ∞ ), the second term is negligible, so E b = E , which is consistent with assumptions of classical homogenisation techniques. Also, if the medium is  comprised of a single material ( β = 1 ), no size effect is produced, In order to investigate size softening or stiffening of a composite ply, two cases are considered. In the first case we assume that the fibres are stiffer than the matrix ( β > 1). From geometrical considerations, 0 ≤ α ≤ 0.5. Therefore, the second term in Eq. (8) is always negative so E b < E . In other words, size softening occurs for β > 1 and E therefore overestimates the effective bending modulus. For the second case with fibres more compliant than the matrix ( β < 1), the second term in Eq. (8) is always positive so E b > E . In other words, size stiffening occurs and E underestimates the effective bending modulus. These results are con-sistent with the work of Wheel et al. (2015) , in which the bending stiffness of a composite beam with K plies of two materials ( E 1 and E 2 ) and thicknesses ( t 1 and t 2 ) has been investigated. In Wheel et al. (2015) , the effective bending modulus for the laminate, E bl , is given by, and the function f was obtained from Euler-Bernoulli beam theory. In this work, the size effect at the individual ply level is investigated, explicitly accounting for fibre/matrix geometry, which is not addressed in Wheel et al. (2015) .

Criterion to define a thin composite ply
In the previous section, we have shown that a heterogenous material can show size dependent behaviour. In order to quantify the size effect, a modulus variation parameter, R , is defined as the relative difference between the bending elastic modulus and the tensile modulus (based on the rule of mixtures), i.e., Thus, R represents the level of heterogeneity-the greater the magnitude of R , the more heterogeneous is the medium. As discussed previously, for fibres stiffer than the matrix can determine for a given material mismatch, β, a value of geometry ratio, α, which maximises R , i.e., We define α 0 to be the value of α satisfying Eq. (12) ,  The largest value of modulus variation, R m , for a given material mismatch ratio, β, is then obtained, Eq. (13) shows that α 0 depends only on material mismatch ratio, β, and is independent of ply thickness, here quantified by number of fibres, N . Furthermore, while α 0 is a local minimum for fibres stiffer than the matrix, it is a local maximum for fibres more compliant than matrix.
While a size effect exists at any size scale, through Eq. (8) , it is important to determine when such an effect is significant (i.e. measurable). For this purpose, we consider that a size effect ex-ists when | R | > R , with R appropriately selected for the problem of interest. We also define a limiting thickness, defined by the number of fibres, N , for a given α and β for which | R | > R . From Eq. (12) this is given by, Thus, a thin ply is a ply with number of fibres N ≤ N satisfying | R | ≥ R . Furthermore, N reaches its maximum value N max for a given material mismatch ratio, β, at geometry ratio α 0 .
In the subsequent sections we examine separately the size effect under bending, quantified through the above equations for the  two cases: β > 1, fibres stiffer than the matrix, and β < 1 fibres more compliant than matrix. In this study R = 0 . 05 .

Size effect in composite plies for fibres stiffer than the matrix
As discussed in Section 2.1 , for fibres stiffer than the matrix, β > 1, the effective bending modulus decreases with decreasing ply thickness (anti size effect). In other words, the elastic modulus calculated by the rule of mixture overestimates the effective bending modulus, and R < 0. The results for the relevant range of α, β and N are plotted in Fig. 2 . Fig. 2 shows that by increasing the material mismatch ratio, β, for the relevant range 0.1 ≤ α ≤ 0.5, the magnitude of R increases, so the size effect is enhanced, E b E .
Therefore, for a constant ply thickness, N , the medium with the stiffest fibres will see the greatest size effect. Note that this does not imply that increasing the fibre stiffness reduces the bending modulus, as the homogenised modulus, E , increases with increasing fibre stiffness.
In Fig. 3 the dependence of R on geometry ratio, α, is shown.
As discussed in Section 2.1 , it is showed that for each value of β there is a value of α for which the magnitude of R is a maximum. Comparing Fig. 3 (a) and (b) also confirms that this value is independent of N . Fig. 4 , shows a schematic diagram of the variation of R with α, where smaller values of α correspond to matrix dominated conditions and larger values of α correspond to fibre dominated conditions. The highest heterogeneity level (largest magnitude of R ) occurs for α = α 0 between these extreme values. Fig. 5 (a) shows that α 0 decreases with increasing β, which means that for larger values of β, the maximum value of modulus variation, R m , occurs at a smaller value of α. As β increases R m becomes almost independent of β and depends only on N . This can be deduced by taking the limit of Eq. (14) ,   β (for α = α 0 ). For example, for β ≥ 40 (a typical value for a glass fibre polymer composite) significant size effects are observed for N ≤ 4.

Size effect in composite plies with fibres more compliant than the matrix
As discussed in Section 2.1 , for fibres more compliant than the matrix ( β < 1), size stiffening occurs, with the effective bending modulus increasing with decreasing size. In other words the elastic modulus calculated by the rule of mixture underestimates the effective bending modulus ( R > 0). Fig. 7 shows values of R for a wide range of α, β and N values. It shows that for a constant value of N , the modulus variation R is larger for smaller values of material mismatch ratio, β, and reaches its maximum at β = 0 , which corresponds to cylindrical voids in a matrix. Thus, if we have two plies with the same thickness, a ply with the smallest value of β shows a greater size effect. This trend is illustrated directly in Fig. 8 where R is plotted as a function of β. In some cases, though not all, there is a local maximum for a given β value at which the greatest heterogeneity occurs, R m . Fig. 9 shows a schematic diagram of the dependence of the modulus ratio R on geometry ratio α for a given value of β. Two possibilities occur-for the curve labelled 2 there is a local maximum corresponding to α 0 calculated from Eq. (13) . For the curve labelled 1, the solution to Eq. 13 lies outside the physical range, α ≥ 0.5 and therefore α 0 is taken equal to 0.5. Eq. (13) shows that values of β ≤ 0.15 correspond to curve 1 in Fig. 9 ; values of β > 0.15 correspond to curve 2 in Fig. 9 as shown in Fig. 8 . Fig. 10 (a) shows the values of α 0 as a function of material mismatch ratio, β. Values of β < 0.15 correspond to case 1 in Fig. 9 , i.e. the value of α 0 lies outside the physical range of α. Fig. 10 (b) shows that the largest value of modulus variation R m 10% for the cases considered is for N = 3 and β = 0 . As for the case of β > 1, Fig. 11 (a) determines the number of fibres, N , for which significant size effects are observed for the range of α and β values examined. Fig. 11 (b) shows that the thickest ply to have significant size effect ( R 5%) has N max = 4 fibres, corresponding to α 0 = 0 . 5 for β = 0 .

Finite-element validation of results
In this section, the results derived in Section 2 have been validated through full 3 D linear elastic FE simulations of a ply under bending load. The FE model is illustrated in Fig. 13 . Fig. 13 (a) shows the geometry of the model. The heterogeneous composite beam is modelled with a representative volume element (RVE) of width d , thickness t = Nd and half length L /2. Symmetric boundary conditions are applied on the back surface and pure bending about the x 3 axis is applied as a linear traction in the x 1 direction on the front surface. Periodic boundary conditions are applied to the left and right surfaces of the model. In order to apply Euler Bernoulli beam theory, consistent with the approximations used in the derivation of Eq. 8 and subsequent results in Section 2 , the length of the beam, L , should be significantly greater than the ply thickness, t ( L / t > 5) ( Kosmatka, 1995 ) and independent of end effects. Therefore the value L = 6 t is chosen for the model. A cross section of a typical FE model is illustrated in Fig. 13 (b) for a model containing three fibres. The commercial FE code, Abaqus, with linear reduced integration hexahedral elements (C3D8R) was employed ( Hibbitt et al., 2018 ). Linear geometry has been assumed, consistent with the Euler-Bernoulli assumption in Section 2 . In the model shown in Fig. 13 (b) there are a total number of 622,459 nodes (1,867,377 degrees of freedom) with 597,240 3D elements. Meshing was done such that the nodes on the left and right faces in Fig. 13 (a) correspond to each other, en- abling periodic boundary conditions on these surfaces to be implemented, as shown in Eq. 17 . Node 2, indicated in Fig. 13 (a), is fully restrained.
u i le f t = u i right − u i node 1 for i = 1 , 2  Fig. 14 (a) shows that close to the surface where the tractions are applied the strain distribution does not correspond to pure bending, due to edge effects. However, sufficiently far from the surface, X 1 > t , the strains are linear in X 2 and do not depend on coordinate X 1 . Thus the Euler-Bernoulli condition of plane sections remaining plane is satisfied, even in the case of a mismatched material. It may, however, be noted that the longitudinal stress ( σ 11 ) is not uniform at the fibre/matrix boundary due to the material mismatch, even far from the surface where the tractions are applied.
In order to calculate the effective bending modulus from the FE model, corresponding to the values obtained in Section 2 , the displacement v = u 2 as a function of X 1 on the mid-plane of the model ( X 2 = t/ 2 ) is obtained from the FE analysis. The curvature, κ, is then determined from Eq. (2) . The total moment, M , through the RVE cross section is obtained from the ap-  A comparison between the FE model and the analytical solutions is provided in Table 1 for N = 3 . The values chosen are those which give the most significant size effect, as determined from Section 2 for N = 3 . In general the FE predictions and the analytical solution are in close agreement. The largest effect for cases with β > 1 is obtained for β = 100 and occurs for both the analytical and FE results. For the case α = 0 . 4 , β = 50 there is some discrepancy between the FE and analytical result. This is thought to be due to numerical errors in the FE mesh. For the cases β < 1 the FE predictions are in good agreement with the analytical results, with the FE model predicting a somewhat larger size effect than the analytical solution.

Equivalent thin ply composite model for a generalized continuum
The bending behaviour of a heterogeneous composite ply has been investigated in Section 2 using Euler-Bernoulli beam theory, with the size effect directly determined through the formulation of the effective bending modulus E b . An equivalent way to capture the same behaviour is to consider a micropolar homogenised medium in which the body is subjected to additional micro-moments ( Yang et al., 2002;Park and Gao, 2006 ). We can thus rewrite the effective bending modulus determined from Eq. (18) to be consistent with the formulation for a micopolar beam with a characteristic length ˆ l , following the approach of Wheel et al. (2015) , Beveridge et al. (2013) , i.e.
where t is ply thickness and in our work t = Nd. Thus we obtain from Eq. (18) (see Appendix B ), Eq. (19) shows that the characteristic length, ˆ l is an explicit function of the microstructure, through the geometry ratio α and material mismatch ratio β. We can rewrite, Eq. (19) in the general

Conclusions
In this work, the size effect resulting from the bending of a composite ply with cylindrical fibres uniformly distributed in a matrix has been investigated through the use of Euler-Bernoulli beam theory. It has been shown that the effective bending modulus depends on fibre/matrix material mismatch, fibre volume fraction and ply thickness. It has been shown that composites with fibres stiffer than the matrix have a negative size effect (effec-tive bending stiffness decreases as the size decreases) so the elastic modulus calculated by the rule of mixtures will overestimate the bending modulus for thin plies. Conversely, size stiffening is predicted when the fibres are compliant than the matrix, so the effective modulus calculated by the rule of mixtures underestimates the bending modulus in this case. These results have been confirmed by 3D linear elastic FE modelling. It is found that significant size effects are observed when the number of fibres through the ply thickness, N ≤ 4. In other words, size ef- fects are significant when the ply thickness is less than 4 times the fibre spacing. Typically, in a fibre reinforced composite material, fibre spacing, d = 15 μm. Thus, we expect size effects to be significant when the ply thickness is less than 60 μm. The results can be interpreted in terms of a bending characteristic length for a micropolar beam, with the characteristic length determined directly from the microstructural length scale. However, for a micropolar formulation this characteristic length is only defined for the case when the fibres are more compliant than the matrix.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. where I c | f is the second moment of area of the fibre about the centroidal axis, as shown in Fig. A.16 (a), and A f is the cross sectional area of the fibre, of radius, r . Similarly, the second moment of area and cross sectional area of the matrix as illustrated in Fig. A.16 (b) region are defined as, where d is the fibre spacing. The bending stiffness, D , of a cross section of width d , containing N uniformly distributed fibres through the thickness, as illustrated in Fig. A.17 , is then given by where E m and E f are the Youngs moduli of the fibre and matrix, respectively and N is the number of fibres. We consider separately the cases of odd ( Fig. A.17 (a)) and even ( Fig. A.17 (b)) numbers of fibres.

Odd number of fibres
The second moment of area about the neutral axis of the cross section (the x -axis in Fig. A.17 ) is given by where n is the current fibre/marix region (ranges from 2 to (N + 1) / 2 ) and the subscript j refers to matrix ( j = m ) or fibre ( j = f ). The second moment of area for n = 1 is given by Eq. (A.1) or Eq. (A.3) .

Appendix B. Characteristic length of a micropolar beam
The equations given in Section 4 are derived in this part of appendix. According to Beveridge et al. (2013) , the stiffness, k , of a micropolar beam with width, d , thickness t and length L under three point bending is given by Eq. (B.1) , Where l is the characteristic length of the beam. In terms of bending stiffness, D , The effective bending modulus E b is then calculated, where I is the second moment of area of the equivalent homogenised cross-section. From Appendix A , where E = E m πα 2 (β − 1) + 1 .

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.ijsolstr.2020.04.004 .