Material Symmetries in Homogenized Hexagonal-Shaped Composites as Cosserat Continua

: In this work, material symmetries in homogenized composites are analyzed. Composite materials are described as materials made of rigid particles and elastic interfaces. Rigid particles of arbitrary hexagonal shape are considered and their geometry described by a limited set of parameters. The purpose of this study is to analyze different geometrical conﬁgurations of the assemblies corresponding to various material symmetries such as orthotetragonal, auxetic and chiral. The problem is investigated through a homogenization technique which is able to carry out constitutive parameters using a principle of energetic equivalence. The constitutive law of the homogenized continuum has been derived within the framework of Cosserat elasticity, wherein the continuum has additional degrees of freedom with respect to classical elasticity. A panel composed of material with various symmetries, corresponding to some particular hexagonal geometries deﬁned, is analyzed under the effect of localized loads. The results obtained show the difference of the micropolar response for the considered material symmetries, which depends on the non-symmetries of the strain and stress tensor as well as on the additional kinematical and work-conjugated statical descriptors. This work underlines the importance of resorting to the Cosserat theory when analyzing anisotropic materials. Author Contributions: Conceptualization, N.F. and P.T.; methodology, P.T.; software, N.F.; validation, N.F., P.T.; formal analysis, N.F.; investigation, N.F.; resources, N.F., P.T. and R.L.; data curation, N.F.; writing–original preparation, N.F.; writing–review and editing, N.F., P.T. and R.L.; visualization, P.T. and R.L.; supervision, P.T.;


Introduction
In the study of composite materials, symmetries play a very important role in identify symmetry planes and peculiar material behaviors. It is well-known that composite materials can be investigated by directly analyzing their constituents in a micromechanical model or by homogenizing them as an equivalent continuum. For composites made of particles and matrix, for instance, the analysis of the single constituents might require high computational cost [1,2]. On the contrary, homogenization is faster and more efficient, but the effectiveness of such a way has to be validated and modeled using proper continuum theories [3][4][5][6]. In particular, the homogenization of complex interaction effects in composite materials needs internal scale parameters which are not negligible compared to the structural length scale.
The paper [7] describes a historical overview of different homogenization techniques which have been extended to non-classical continua. As it is well known in the literature, the classical continuum, by lacking in internal lengths, leads to ill-posed problems. This calls for the need of non-local descriptions which include internal length parameters in different ways; or by adding extra degrees of freedom, thus obtaining the so-called implicit/weak non-local descriptions, or by adding extra

Framework of Cosserat Theory
In the present work, a two-dimensional framework is considered for simplicity. The Cosserat kinematics take into account an additional degree of freedom of particle rotation ω other than classical displacement components u 1 and u 2 .
The kinematic compatibility equations in matrix form are where ε ij , i, j = 1, 2 and k i , i = 1, 2 are the strain components: axial, angular and curvature ordered in the strain vector {ε} T = {ε 11 ε 22 ε 12 ε 21 k 1 k 2 }. Note that the angular strain components are not reciprocal, ε 12 = ε 21 . In general in the micropolar model the microrotation, ω, is different from the local rigid rotation (macrorotation), θ, defined as the skew-symmetric part of the gradient of displacement, and the difference between the two rotations, θ − ω, defines the strain measure of the relative rotation that corresponds to the skew-symmetric part of the strain. When the relative rotation equals zero, ∂y + ∂u 2 ∂x , the micropolar continuum becomes a continuum with constrained rotations [31,59]. In the following, we focus on θ − ω as a peculiar strain measure of the micropolarity of the model under study.
The work conjugate stress vector to {ε} is written in the form {σ} T = {σ 11 σ 22 σ 12 σ 21 µ 1 µ 2 }, where σ ij for i, j = 1, 2 represent the normal and shear stress components, and µ i , i = 1, 2 the microcouples. The stress components are not reciprocal, σ 12 = σ 21 , and the couple stress components, µ 1 , µ 2 , have to be introduced in order to satisfy the moment equilibrium of the micropolar body. The Principle of Virtual Works can be written as where { f } and {p} are the vectors of body forces and boundary tractions, respectively. Balance equations and correspondent boundary conditions can be carried from Equation (2).
For the sake of conciseness, boundary terms are not reported. For further details and reading, strong form equations can be found in [50,60,61].
The micropolar constitutive equations take the compact form and extended matrix form Due to hyperelasticity, the constitutive matrix C in Equation (5) is symmetrical. This matrix is obtained via an energetic homogenization technique presented in [47]. Such homogenization starts from the definition of a Representative Volume Element (RVE) in which orientation and location of elastic interfaces are defined according to a hexagonal geometry selected and described in Section 4. Further details on constitutive elastic terms are given in the following sections.

Finite Element Implementation
A standard finite element implementation is considered for the numerical solution of the present Cosserat problem. The numerical framework follows the previous contributions of the authors [50,60,61]. An in-house finite element MATLAB code has been developed as an extension of a classical two-dimensional elastic continuum as given in [62].
A linear finite element model is considered with quadrilateral (Q4) elements. Linear Lagrange interpolation functions are used; thus, the approximation takes the form {u} = N{d e }, where {u} is the vector of degrees of freedom of the model, N is the matrix of shape functions (of size 3 × 12) and {d e } the vector of nodal displacements.
By inserting the finite element approximation in the strain definitions (1) and in the variational form of the equilibrium (2) the stiffness matrix and load vectors can be derived as where [B] is the matrix which includes the derivatives of the shape functions. Note that selective Gauss integration is performed on the shear terms due to different differentiation orders among displacements and rotational approximation. In the present work, body forces are neglected { f } and this is why they do not appear in the force vector definition above.
The present finite element approximation is based on displacements; therefore, derivative quantities are post-computed according to a 2 × 2 Gauss-Legendre grid, and inter-element averaging is done in the nodes without any extra weight [60].

Numerical Applications
The hexagonal pattern of Figure 1 has been defined according to the parametrization [58]. The single tile is defined by the relative length l r which identifies the unit parallelogram. The selected RVE is highlighted within the centroids of the tiles and outward unit normal vectors at the block interfaces used for computing the constitutive matrix according to the procedure presented by [47]. The relative length defines the ratio between AB ≡ DE and BD ≡ AE. The single tile is defined by a rectangle and two isosceles triangles with the base attached to BD and AE edges of the rectangle. The nodal coordinates of the single tile are given by: , l 1 = 1 − l 2 , l 4 = 0.5l 2 tan α 2 , l 5 = 0.5l 2 tan α 3 . Due to this definition, the geometric range of parameters α 2 and α 3 is Thus, by keeping the relative length l r constant and varying α 2 and α 3 , all possible configurations (not self-intersecting hexagons) are depicted in Figure 2. Note that the main diagonal represents the case wherein α 2 = α 3 and the horizontal and vertical lines represent α 3 = 0 and α 2 = 0, respectively. The two orthogonal lines divide the space in four areas: 1. The upper-right area has α 2 > 0 and α 3 > 0 where any combination is theoretically possible in the given domain range. In this area, regular hexagons are present as well as diamond shaped ones (according to Figure 3). 2. The lower-left area has α 2 < 0 and α 3 < 0 where not all shapes are achievable due to some self-intersections of hexagons. In this area, hourglass shaped hexagons are present (symmetric and not-symmetric ones) according to Figure 3. 3. The upper-left and lower-right areas have α 2 < 0 and α 3 > 0 and α 2 > 0 and α 3 < 0, respectively.
These shapes might be characterized by an asymmetric shape (according to Figure 3).
The constitutive matrices of several hexagonal configurations are considered in the following. The homogenization procedure follows the approach described by [47], where the adopted spring stiffness values at the elastic joint interfaces are K N = 785. Energetic equivalence is used to carry out rotational stiffness as: where d is the current interface length between two rigid particles in contact, for which the interactive stiffness is computed. Note that the definition of k t is a particular choice which depends on the Poisson coefficient of the heterogeneous medium. Considering all the possible geometric configurations illustrated in Figure 2 the equivalent mechanical properties of the equivalent micropolar material can be carried out according to the constitutive matrix (5) as in [50,61]. The nonzero components of the classical constitutive matrix A are depicted in Figure 4. Please note that for all possible hexagonal combinations, A 1112 = A 1121 = A 2212 = A 2221 = 0. Similar trends are observed for A 1111 and D 11 , as well as the couples A 1122 and A 1221 , finally for A 2222 and A 1212 as shown in Figure 4. Note that stiffnesses on the main diagonal are related (as well as off-diagonal terms which are related to Poisson effect and shear asymmetry). A 1111 and A 2121 increase as α 2 or α 3 increase on the contrary A 2222 and A 1212 have an opposite behavior (they increase when the two angles decrease). Off-diagonal terms A 1122 and A 1221 are positive when α 2 = α 3 > 0; they have negative values otherwise. Positive or negative values are also observed when α 2 = −α 3 > 0 and vice versa. The latter occur for asymmetric configurations of the tiles according to the nomenclature introduced in [60].
Coupling between stresses/curvatures and microcouples/strains are observed when chiral material behavior is present. As presented in [60], chiral behavior occurs when asymmetric tiles are considered for α 2 = −α 3 . As a matter of fact, B ijk components are zero for α 2 = α 3 as depicted in  Thus the constitutive relations for the present case are always in the form Given the aforementioned presentation of the constitutive configurations, the results for the following material configurations are shown: And their reference RVE is depicted in Figure 6. The constitutive matrices of the configurations selected and depicted in Figure 6 are listed below. For the rectangular configuration, it is The constitutive matrix for the hourglass configuration ( Figure 6b) is As it was already observed in [60], an auxetic (negative Poisson ratio) nature of the material is observed when hourglass tiles are considered.
The constitutive matrix for the diamond configuration (Figure 3c) is The diamond configuration results to be similar (but with higher stiffness coefficients for A 1111 and A 2121 due to the elongated shape) with respect to the regular one (see below or [60]) with positive Poisson contraction and uncoupled condition between classical and micro-couple quantities.
The constitutive matrix for the regular configuration ( Figure 3d) is the regular configuration results to be like an isotropic material as described already in [60]. The constitutive matrix for the skew configuration ( Figure 3e) is The skew configuration, as indicated above, has a chiral behavior due to the coupling of stresses and curvatures and microcouples and strains. Large values of coupling stiffness are present due to the large angles for |α 2 | and |α 3 |.
The constitutive matrix for the tip configuration ( Figure 3f) is Note that the present configuration has all the constitutive components of the matrix in Equation (8) not zero. In other words, a Poisson effect is present as well as a coupling between stresses/curvatures and microcouples/strains. Given the present hexagonal geometry with elastic interfaces which interact only through normal stiffnesses, this is the most general configuration possible. Table 1 lists in compact form all the nonzero constitutive stiffnesses according to the geometry selected. It is clear from this table that tiles with a vertical symmetry (rectangular, hourglass, diamond, regular) do not correspond to a chiral material as B ijk = 0. In the following, for all the six configurations depicted in Figure 6, an elastic problem is solved via the finite element method. A square panel of side L = 4 is considered and subjected to a top pressure on a limited area a = L/4 with a resulting equivalent concentrated force of P = 10 3 pointing downwards. The panel is clamped at the bottom. The problem is symmetrical for geometry and loads so only half of the domain is investigated. A regular squared finite element mesh of 16 × 32 Q4 elements is considered and depicted in Figure 7.
Horizontal displacement u 1 is depicted in Figure 8. It is evident that rectangular and regular hexagons correspond to orthotropic and orthotetragonal behaviors which in this case are very similar; the only difference between these two configuration is that orthotropic configuration has no Poisson effect. The hourglass configuration, which is auxetic (Poisson negative), has a stronger contraction of the material (larger horizontal displacement) close to the applied load as it is expected from auxetic materials. The other configurations that have all horizontally elongated tiles display very small horizontal motions, either positive or negative.
Vertical displacement u 2 is depicted in Figure 9. Rectangular and regular hexagons also have a similar behavior (except for the Poisson effect). The effect of negative Poisson ratio is visible in the hourglass configuration which displays concentrated vertical displacements and is close to zero displacements far from the applied load. This is due to the horizontal contraction of the material. Diamond and tip configurations, one not-chiral and one chiral, respectively, show similar behavior due to the horizontally elongated tiles. On the contrary, the smallest vertical displacements are measured for the skew configuration due to the engagement of tiles.
The horizontal stress σ 11 is shown in Figure 10 does not reflect the map given by the horizontal displacement ( Figure 8) due to the coupling between classical and micropolar quantities in the skew and tip configurations. In fact, the diamond configuration which registers small horizontal displacements has very localized horizontal stresses. Rectangular, auxetic and regular configurations have behaviors expected and already shown in [60]. Note that tip configuration, to which corresponds a small horizontal displacement, has a high and localized horizontal stress wherein the maximum value of the stress does not lay upon the symmetry axis as in all the other configurations. This translation of the maximum horizontal stress is due to the coupling between the horizontal stress σ 11 and the curvature k 2 due to B 112 = 0. Vertical stress σ 22 is shown in Figure 11. Stress percolation due to micropolar effect is shown by rectangular, auxetic, regular and skew configurations. Diamond and tip configurations display highly concentrated vertical stresses with stress fields that are negligible and close to the bottom boundary condition. The more evident stress percolation is observed in the hourglass configuration due to the higher D 22 /D 11 ratio with respect to the others. The diamond configuration which has the smallest D 22 /D 11 ratio observes a higher vertical stress concentration close to the applied load.
Similar to what has been observed in the previous works [50,60], the relative rotation, represented in Figure 12 is affected by not symmetric classical stiffnesses A 1212 = A 2121 , chiral behavior B ijk = 0 and micropolar stiffnesses D 11 , D 22 . As an example, the rectangular configuration, which has an orthotropic behavior, shows a higher relative rotation than the skew and tip configurations, which are chiral. The same can be said if the rectangular configuration is compared to the regular one, which is orthotetragonal, and results in having a smaller relative rotation. As for the vertical stress case, once again, the hourglass configuration has the more evident relative rotation and on the contrary, the diamond configuration has almost zero relative rotation.
In the following unsymmetrical shear strains, ε 12 and ε 21 are discussed and displayed in Figures 13 and 14. Classical behavior should be retrieved if ε 12 = ε 21 ; in other words, no micropolar effect should be observed. ε 12 and ε 21 include Although with some (in some cases small) variations, a micropolar effect is always observed for all configurations due to the fact that D ij = 0 in all cases and also A 1212 = A 2121 and B ijk = 0 for some cases. As expected, the configurations that observed a relative rotation close to zero in most of the solid (diamond, skew and tip) have ε 12 that is a classical continuum. At the same time, the highest micropolar effect is observed in the auxetic behavior which had the largest variation of the relative rotation contour plot.

Conclusions
This paper investigated the static behavior of materials with hexagonal microstructures. A homogenization approach is applied to selected lattice assemblies in order to obtain a micropolar constitutive relation which results to be a function of two geometric parameters. Up to twelve nonzero elastic constants have been activated with the present procedure and the particular behavior of six material configurations (termed rectangular, auxetic, diamond, regular, skew and tip) have been presented. The selected patterns correspond to the so-called orthotropic, orthotetragonal, auxetic and chiral material symmetries. In particular, chirality occurs when asymmetric tiles are selected and an auxetic solid is derived when rehentrant corners are used.
Different levels of nonlocality, that in the case of micropolar media are of the 'implicit' type, that are related to the geometric configuration, have been observed to be due to the bending moduli ratio, D 22 /D 11 ratio, and to chiral terms (B ijk = 0). However, chiral materials did not show the larger degree of micropolarity given by the auxetic configuration. None of the configurations correspond to classical continua, because no configuration has zero micropolar terms; however, the so-called diamond and tip configurations result in the ones with largest contour areas of zero relative rotation, that is, areas with strain that is quite symmetric, and then behave in a way that is similar to a classic material. Future developments of this work will be oriented to wider investigations of the effect of the relative rotations with a reliable measure of nonlocality.