New insights on homogenization for hexagonal-shaped composites as Cosserat continua

In this work, particle composite materials with different kind of microstructures are analyzed. Such materials are described as made of rigid particles and elastic interfaces. Rigid particles of arbitrary hexagonal shape are considered and their geometry is described by a limited set of parameters. Three different textures are analyzed and static analyses are performed for a comparison among the solutions of discrete, micropolar (Cosserat) and classical models. In particular, the displacements of the discrete model are compared to the displacement fields of equivalent micropolar and classical continua realized through a homogenization technique, starting from the representative elementary volume detected with a numeric approach. The performed analyses show the effectiveness of adopting the micropolar continuum theory for describing such materials.


Introduction
Composite materials can be investigated by directly describing their constituents in a micromechanical discrete model or by homogenizing them as equivalent continua. In the former case, for example to study heterogeneous materials such as polycrystalline materials, jointed rock systems, and block masonry with periodic microstructures, the model involves a system with a large numbers degrees of freedom and as a consequence the computational costs is very high [1][2][3][4]; such cost increases by reducing the material scale [5]. In the latter case, by using homogenization techniques, the numerical analyses are faster and more efficient but it is fundamental to define equivalent continua that properly takes into account the influence of the microstructure on the macro-scale behaviour with particular reference to shape, size and texture of elements [6] and among them focus on those techniques developed for deriving scale-dependent models [7][8][9][10][11][12]. When classical kinematics is enriched with extra degrees of freedom for instance, homogenization procedures have been shown to provide more reliable models than in the case of classical, local, continua, even in the case of higher order continua [10,[13][14][15].
The selection of proper non-local continuum theories which retain memory of the fine internal structure of the material, shape, size and texture, is then a fundamental point in the macroscopic description of composite materials. Among various homogenization techniques, the development of methods and conceptual guidelines for obtaining continuous field descriptions by linking discrete and continuum solid mechanics have been successfully employed material science. In this work, a micromechanical approach based on the scale dependent discrete-continuum homogenization procedure defined in [10,35] for particle composite materials, addressed to derive the macroscopic mechanical properties of a micropolar continuum, is adopted. The micropolar theory, beside the standard displacement, introduces the rotation of the material point, termed microrotation, to be distinguished with the macrorotation (local rigid rotation). When the relative rotation, i.e. the difference between micro and macrorotation, is null the micropolar continuum, under some other dynamical costraints, behaves as a couple stress continuum [13]. As the effects of the relative rotation are expected to be prominent for strongly anisotropic particle composites, as widely investigated in [14,15] for masonrylike materials, we consider here the full (unconstrained) micropolar theory.
In order to highlight the advantages of micropolar elasticity description, the numerical solutions of equivalent micropolar and classical continua, obtained by putting some internal constraint to the micropolar model, are compared to the one of a discrete model assumed as benchmark solution. The discrete numerical solution has been obtained using the ABAQUS software in which the single blocks are rigid and elasticity is concentrated at the interfaces through elastic linear springs (joints) while for the continuum, classical and micropolar, solution, an in-house FEM code has been performed [37,54]. Moreover, an investigation is proposed in order to detect the representative volume element for hexagonal shaped composites. This geometry is typical for some polycrystals with thin interfaces such as Aluminia (Al 2 O 3 ), Zirconia (ZrO 2 ), Zinc Ozide (ZnO) or Tungsten-Carbide (WC). Considering a generic tile (e.g., orientation of the interfaces and internal angles of the hexagonal geometry) different material symmetries can be detected [55]. The hexagonal shaped composites are investigated starting from the previous works [37,39] and in particular considering three microstructure geometries [40].
In this paper, at first Cosserat theoretical aspects are briefly presented (Sect. 2) and then the geometries of microstructured composites are be introduced, in particular by studying the unite cell to detect the RVE necessary for the homogenization procedure (Sect. 3). Later the numerical implementation of the continuum model is shown (Sect. 4), the discrete model assumed as benchmark solution is properly described and finally a static analyses comparison among the discrete system, the micropolar and the classical continuum for a 2D panel is reported (Sect. 6).

Cosserat theoretical framework
The Cosserat model can be considered equivalent to an assembly of rigid particles undergoing homogeneous displacements and rotations which interact with each other via forces and couples. Therefore, it is suitable for simulating the mechanical response of granular, masonry-like and composite materials with either periodic or random microstructure [10,33,36].
For the present study a 2D rectangular panel with three different hexagonal microstructures is considered. In a linearized 2D framework there are two displacements and one rotation component, so that the generalized displacements u T ¼ u 1 u 2 x ½ applies. The strain displacement vector is e T ¼ e 11 e 22 e 12 e 21 k 1 k 2 ½ , where e 11 ; e 22 ; e 12 ; e 21 are the in-plane normal and shear strains and k 1 , k 2 are the micropolar curvatures. Note that the strain components are not symmetric e 12 6 ¼ e 21 . Moreover, the microrotation x is different from the local rigid rotation (macrorotation) h, defined as the skew-symmetric part of the gradient of and the difference between the two rotations h À x defines the strain measure of the relative rotation, corresponding to the skew-symmetric part of the strain. When the relative rotation equals zero h ¼ x then e 12 ¼ e 21 ¼ ou 2 ox 1 À ou 1 ox 2 as in the classical continuum, the micropolar continuum becomes a continuum with constrained rotations, that under other dynamical constrains, is a couple stress continuum [13,56]. The stress vector is represented as: r T ¼ r 11 r 22 r 12 r 21 l 1 l 2 ½ where r ij for i; j ¼ 1; 2 represent the classical normal and shear stress components and l 1 , l 2 are the microcouples. The stress components are not reciprocal, r 12 6 ¼ r 21 and the couple stress components l 1 , l 2 have to be introduced in order to satisfy the moment equilibrium of the micropolar body.
The kinematic compatibility relations can be written as: where: From the virtual work principle, equilibrium equations in terms of stresses and microcouples can be carried out as dU þ dV ¼ 0, being dU the internal work: and dV the potential of external loads, that in 2D domain takes the form: where f and p are the vectors of body forces and boundary tractions, respectively. Balance domain equations are then given by: and boundary tractions as: wherep are respectively, contact forces and couples applied at the boundary. The micropolar anisotropic constitutive equations take the form: where: Due to hyperelasticity, the constitutive matrix is symmetric (C 2 Sym).
Accounting for the constitutive equations, the internal work can be write: In Sect. 4 the principle of virtual work will be adopted for the numerical implementation.

Microstructure geometry
The geometry of the microstructure is realized considering a six edges parallelogram and it is possible to give different input settings to set the shape: the main parameters are the angles a 1 , a 2 and a 3 and the relative length ( Fig. 1 þ 1Þ that defines the ratio between the length of the base AE (or BD) and the l 1 (perpendicular distance between AE and BD), l 1 , l 2 are defined as: [40,57].
Keeping constant the lengths, and varying the angles a 2 and a 3 it is possible to obtain different geometries; in particular for a 2 ¼ 0 and a 3 ¼ 0 the lateral edges of hexagonal shape are vertical therefore the parallelogram acquires a rectangular shape. In this paper, for the patterns studied, the angles a 2 and a 3 assume the values: a 2 ¼ a 3 ¼ 30 for the regular hexagonal shape, a 2 ¼ a 3 ¼ À20 for the hourglass shape and a 2 ¼ Àa 3 ¼ À30 for the skew shape. Finally, by varying the relative scale the three geometries reported in Fig. 2 are obtained by setting the scale factor s: where s ¼ 1 for large blocks, s ¼ 0:5 for medium blocks and s ¼ 0:25 for small blocks. The ratios among the parameter l 2 of microstructure shape and the vertical side L y of panel analyzed in Sect. 5 are:

Representative volume element
The homogenization procedure adopted for deriving the material constants of the micropolar continuum is the one defined in [10], where at the fine level the material is described as a system of rigid blocks, here named elements, interacting through elastic contact elements, here referred as joints (e.g. springs of normal stiffness k 11 ¼ 785 N/m and shear stiffness k 22 ¼ 392:5 N/m), the interaction is represented by forces and couples. A kinematic correspondence map between discrete and continuum fields is assumed (generalization of the Cauchy-Born rule) and an energy equivalence criterion is adopted. To apply the homogenization procedure for periodic assemblies it is necessary to firstly define the Representative Volume Element (RVE), as the elementary volume element (unit cell) made of the minimal number of elements and joints sufficient to completely define the behavior of the discontinuous and heterogeneous material. As the micropolar continuum is intrinsically For the following analyses it is convenient to represent the constitutive matrix (8) in the form: It is possible to show that only the constitutive components of B and D contain internal length scales, and this also justify the classification as 'implicitly' non local material [29]. In order to define the RVE, here by exploiting the results of numerical simulations performed comparing micropolar and discrete solutions (Sect. 3.1.2), with

7-Blocks unit cell
In the first case, the unit cell is constituted of one central block and six joints, therefore seven elements are considered (Fig. 3). The common edges interact by the medium point of the joint and the calculations of material constants are referred to the gravity center of the cell.
The constitutive matrix for regular hexagonal shape is:  Finally, the constitutive matrix for skew shape is:

4-Blocks unit cell
In this case the unit cell considered is constituted of a system of four blocks with five joints, centered at the central joint (Fig. 4). The constitutive matrices are different from the cases above. The new matrices are: For the regular hexagons: It can be useful to observe that the constitutive matrices for the regular hexagon shapes are not orthotetragonal, the terms A 1111 6 ¼ A 2222 , A 1212 6 ¼ A 2121 as in the 7-blocks regular hexagon case because the joints are not in a symmetric configuration, invariant with respect to p 4 , therefore in this case the symmetry is lost homogenizing from the discrete system to the continuous system in an equivalent energetic approach [10]. The micropolar term D is scale dependent with the same ratio as in the above case, approximately D s¼1 reg % 4 D s¼0:5 reg % 16 D s¼0:25 reg . The constitutive matrix for the hourglass shape is: The constitutive matrix for the skew shape finally is:   It is useful to note that the Cauchy constitutive matrix does not change with the scale variation, differently from the Cosserat matrix.
The question of the RVE definition is still open and the answer can depend on the type of microstructure and on the consideration that in an energy equivalent discrete-continuum transition the material symmetry class of the discrete assembly must be preserved [10].
Here a numerical criterion is adopted. In order to select the RVE for the homogenization procedure here adopted, the external load work of the static analysis case studied in Sect. 5 has been calculated for the discrete assembly and for the micropolar model of both 7-blocks and 4-blocks unit cells. The comparison refers to the static analysis of a rectangular panel where the microstructure factor scale is 0.5 and the results are reported in Table 1. The analyses relative to the the scale 1 and scale 0.25 are omitted because they give analogously consequences. Basing on these results, the 7-blocks unit cell of Subsect. 3.1.1 has been selected as the RVE for the numerical simulations performed; furthermore the Cauchy constitutive matrices are calculated from the matrices showed in Subsect. 3.1.1.

Numerical implementation
In the next subsections a numerical implementation of the micropolar theory is resumed and the discrete models adopted for the simulations are described.

Continuum model
The present finite element framework is based on the previous works [15,38]. The finite element enforces an approximation through nodal kinematic parameters as: The displacement vector for the present problem is given by: d eT ¼ u 1 1 ::: u 4 1 u 1 2 ::: u 4 2 x 1 ::: the finite element has 12 degrees of freedom. The matrix of the shape functions takes the form [37,38] :  where N is the vector of the linear Lagrangian shape functions. The internal work becomes [39]: where B ¼ D N , thus the element stiffness matrix is: The potential of loads becomes: where F e and P e are volume and surface force vectors respectively. Linear shape functions are considered for the present implementation with reduced integration for the micro-rotation x [40]. To perform reduced integration the strain vector has to be reordered by separating strain terms which are fully integrated and the ones for which reduced integration is applied.
Once the problem is solved in terms of displacements other quantities such as stresses and relative rotation have to be post computed [40].

Discrete model
The discrete solution has been obtained using the FEM software ABAQUS; because of the complexity to build the model through the program interface, Python script codes have been prepared. Static analyses with three different discrete models are presented as a benchmark for the numerical study in order to show the goodness of the adopted homogenization technique [10].

Spring model
In the first model three linear elastic springs, longitudinal, transversal and rotational, are assigned at the center of each joints. For the edges parallel to axis x or y, springs are represented in the joint local reference system. The scripts are structured in the following principal steps: 1. The geometry of the hexagonal shape is defined; 2. The material properties are defined; 3. The parts of the model are created; 4. The boundary blocks of the panel are cut to define the rectangular geometry panel; 5. Edges are partitioned to define the middle point for each joint; 6. Rotational and translational springs are added for all internal joints.
The single block must behave as a rigid body, therefore the material is assigned as an elastic isotropic behavior with a high value of the Young modulus with respect to the joint stiffness in particular E ¼ 100 MPa. The normal stiffness is represented by the k 11 term in reference to the normal direction to the contact area and the tangential stiffness by the term k 22 along the contact area plane in the local system (Fig. 5). Energetic equivalence is used to carry out the rotational stiffness [39]: Fig. 5 Example of a single tile for the hourglass shape with elastic springs and local reference systems where d is the generic edge length of the single tile. Finally the boundary conditions, loads and the mesh of model are manually input.

Contact interaction
The elements, the joints and the texture of the discrete model are realized through Python script. However, in the present implementation the interactions between rigid blocks is considered as distributed contact with the ABAQUS interaction option ''surface interaction''. The pressure-overclosure relationship are prescribed by using a linear elastic law: the surfaces transmit contact pressure when the overclosure between them, measured in the contact (normal) direction, is greater than zero. The linear pressureoverclosure relationship is identical to a tabular relationship with two data points, where the first point is located at the origin, it is necessary to specify the slope of the pressure-overclosure relationship, k. For the the tangential direction it is set the friction formulation as ''frictionless'' to consider the case with no tangential stiffness, in this way ABAQUS assumes that surfaces in contact slide freely without friction. It is considered this particular case for the difficulty to define a linear elastic tangential behavior as done for the normal direction. Various numerical tests have been done and the results are very close to the springs model, obviously with only a normal and rotational springs; the quality of the test indicates the right calculation for the rotational stiffness.

Elastic Interfaces
The third approach considers linear elastic interfaces among the blocks. Two parts in ABAQUS are created: one for the rigid blocks and one for the elastic interfaces. The mechanical properties of the materials were defined for the rigid block as orthotropic behavior with a high Young modulus and for the elastic interfaces engineering constants are set in order to fictitiously model a frictionless elastic interface with G 13 ¼ G 23 ¼ G 12 ¼ 0. The Young modulus was calculated as EA=l ¼ k n s=l, therefore E ¼ k n s=A where: E is the Young modulus, A is the interface contact area, A ¼ td; with d the edge length (because in this case the contact is developed along the whole surfaces) and t is the panel thickness. k n is the normal spring stiffness and s thickness of the layer. Thus, assuming a unitary thickness t ¼ 1, the modulus becomes: The Young modulus E increases with the decrease of the scale. Subsequently from the blocks and interfaces parts two instances were created: then they were bonded by the ''tie constraints'', that give them the same displacements at the contact interface. For explanatory purposes, using the three discrete systems previously described, in Fig. 6 a comparison of the vertical displacements for a panel clamped at the bottom with a vertical concentrated force at the center of the upper edge is reported. From the vertical displacement contour plot a greater similarly between the spring and contact models can be detected, in particular for the panel region far from the applied load. The discrete model with elastic springs has been used for a comparison with the equivalent micropolar continuum model: this case is closer to the theoretical discrete model presented in [13]. Moreover, the springs model is faster with respect to the model with contact interactions for which the computational cost increases exponentially with the scale reduction. In addition, contact interfaces leads to a non linear model on the contrary springs interfaces are linear.

Simulations
In this section, the behavior of 2D rectangular panels with different microstructural shape and different scales is investigated in Fig. 7, the panel is clamped at the base and a vertical load is applied at the top side center, for a footprint length equals a ¼ L x =4, L x being the base length. The resultant of the load is a force of intensity P ¼ 10 N, the stiffness terms are k 11 ¼ 785 N/m and k 22 ¼ 392:5 N/m and the length of panel sides are L x ¼ 6:6 m, L y ¼ 7:7 m for the regular hexagonal shape, L x ¼ 5 m, L y ¼ 7:7 m for the hourglass shape and L x ¼ 5:85 m, L y ¼ 7:7 m for the skew shape.
In the next sections, a displacement field comparison among the discrete spring model and the micropolar and classical continuum models is reported.

Regular hexagonal shape
Horizontal components of displacement are shown (Fig. 8): the contour plot of the micropolar continuum is in line with the discrete system. There are not evident differences among the Cosserat and Cauchy model. This is expected because of the orthotetragonal configuration of regular hexagons, that implies corresponding results among all the models [14].
For vertical displacements the micropolar continuum contour plot is in line with the discrete system behavior, once again because of the orthotetragonal configuration no evident differences among the Cosserat and Cauchy model are revealed and their differences vanish when the microstructure scale factor become smaller (Fig. 9).

Hourglass shape
Horizontal displacements are shown (Fig. 10): the contour plot of the micropolar continua is in line with the discrete system. At the top of the panel a smaller maximum displacement increase can be noted with the scale reduction for the discrete system.   Also for the vertical displacements the micropolar continuum contour plot is in line with the discrete system behavior (Fig. 11). The effect of negative Poisson ratio is visible in this configuration which displays concentrated vertical displacements which are close to zero far from the applied load, this is due to the horizontal contraction of the material. A major displacement distribution can be noted for the micropolar model respect to the classical model: moreover, for the hourglass shape a major percolation [58] is visible for the Cosserat model respect to the Cauchy continuum and this behavior is closer to the discrete system plot.

Skew shape
The micropolar continuum model reproduces accurately the discrete system behavior, once again for the horizontal displacements there are not obvious differences among the three continuum models. For the discrete system, with the scale reduction a smaller displacement magnitude increase can be noted at the panel top (Fig. 12).
Also for the vertical displacements field the micropolar continuum gives a faithful representation for the discrete system contour plot and once again the Cosserat model reveal a major distribution of the level curves in comparison with the Cauchy model (Fig. 13). Lastly as already shown for the hourglass shape, in the skew shape contour plot the percolation effect is more evident for the micropolar model.

Conclusion
The present work investigates the static behavior of composite materials with three types of hexagonal microstructures as equivalent micropolar media. The comparison among the discrete model and the continuum models highlighted that assemblies of regular hexagons have an orthotetragonal behavior, therefore there is only a weak dependence on the scale, as in the case of materials with elements of small size, and it has been shown that their homogenized behavior is close to the behavior of classical elastic bodies. For the hourglass and skew configuration, the scale dependence is more evident both for the discrete system and for the micropolar continuum, which properly accounts for the non-symmetries as well as size effects exhibited by the discrete system, as consequence, this last should be preferred in this cases. It is then highlighted the scale dependence of the discrete system and how the Cosserat model takes into account it, differently for the Cauchy model, and how to consider the asymmetric behavior of anisotropic materials is fundamental.
Moreover, the static analyses comparison shows that the homogenization technique can be generalized for hexagonal geometry microstructure, both for concave shapes and convex shapes [10,13,35]. For the study of particle composites the possibility to use a proven homogenization technique involves in a faster and efficient computational modeling. Furthermore three discrete models have been used to study composite materials described as materials made of rigid particles and elastic interfaces. The comparison between them show the right assumption to calculate the interface rotational stiffness. The elastic spring model is a linear model more faster for the numerical simulation compared to other two and it well represent the mathematical theoretical model adopted for the discrete system representation. Finally, other geometries can be investigated in the future comparing the discrete system and the continuum model and as well as the dynamic case can be studied.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.