A micromechanical approach for the micropolar modeling of heterogeneous periodic media

Computational homogenization is adopted to assess the homogenized two-dimensional response of periodic composite materials where the typical microstructural dimension is not negligible with respect to the structural sizes. A micropolar homogenization is, therefore, considered coupling a Cosserat medium at the macro-level with a Cauchy medium at the micro-level, where a repetitive Unit Cell (UC) is selected. A third order polynomial map is used to apply deformation modes on the repetitive UC consistent with the macro-level strain components. Hence, the perturbation displacement field arising in the heterogeneous medium is characterized. Thus, a newly defined micromechanical approach, based on the decomposition of the perturbation fields in terms of functions which depend on the macroscopic strain components, is adopted. Then, to estimate the effective micropolar constitutive response, the well known identification procedure based on the Hill-Mandel macro-homogeneity condition is exploited. Numerical examples for a specific composite with cubic symmetry are shown. The influence of the selection of the UC is analyzed and some critical issues are outlined.


INTRODUCTION
he use of composite materials in various fields of engineering, both for standard and innovative applications, has been widely researched.A thorough understanding of the mechanical behavior of existing materials is a fundamental step towards the design of new composites, characterized by increasingly high performances.Various approaches, marked out by different formulations and modeling the materials at different scales, have been proposed to deal with the constitutive response of composite materials.This study focuses on homogenization techniques, a very effective tool to obtain accurate results with low computational efforts.The actual heterogeneous medium is analyzed at two different scales: the macro-scale, where an equivalent homogenized medium is considered, characterized by overall effective mechanical properties, and the micro-scale, where detailed information about the texture, geometry and constitutive laws of the constituents are available.Different continuum models can be adopted, at the two levels.The classical Cauchy continuum provides an appropriate description of the actual heterogeneous response in the case of small microscopic length compared to the macro-scale structural length [1,2,3,4].On the contrary, when strong strain and stress gradients at the macro-level occur, or when the microscopic length of the constituents is comparable to the wavelength of variation of the strain and stress mean fields at the macro-level, some intrinsic limits emerge.This is due to the fact that the Cauchy theory does not account for length T scales.By adopting generalized continua this limit is overcome.Many authors have focused on coupling different continuum models at the two scales.In most cases, at the microscopic level, the classical Cauchy continuum is adopted, especially because nonlinear constitutive relationships are well-established in this framework.Among various generalized continua (second-gradient, couple stress, micropolar or multifield), a micropolar Cosserat continuum at the macro-level and a Cauchy continuum at the micro-level are used here to study the homogenized response of periodic composite materials.The computational homogenization technique adopted effectively predicts the macroscopic behavior of composite materials [5,6].Since composite materials, characterized by regular textures are analyzed, a Unit Cell (UC) is selected at the micro-level.Consistently with the strain-driven approach, the two levels are linked through a kinematic map based on a third order polynomial expansion, previously derived by the authors in [7], from an original idea developed in [1].The displacement field at the micro-level is represented as the superposition of the kinematic map and an unknown perturbation field, due to the heterogeneous nature of the material.For classical first order computational homogenization the perturbation fields are periodic [8], but this is not true when higher order polynomial terms are considered, as underlined in [9,10].In this study, the problem of determining the displacement perturbation fields, with particular reference to its influence on the structural response evaluation, is investigated.To this end, the three following techniques are adopted.The first is based on the solution of the Boundary Value Problem (BVP) by applying periodic boundary conditions on the UC.The second is the 3 steps homogenization, presented by the authors in [11], and based on the decomposition of the perturbation fields in terms of functions which depend on the macroscopic strain components.Finally the BVP is solved by applying special boundary conditions on the UC, as derived in [7].Furthermore, the identification of the homogenized linear elastic 2D Cosserat constitutive parameters is performed, by using the Hill-Mandel technique, based on the generalized macro-homogeneity condition presented in [11].As known, this technique has inherent limitations, leading to physically inconsistent results [9,12,13].For example, higher order constitutive components are identified, also when a homogeneous elastic material at the micro-level is considered.Despite the drawbacks, this technique has been widely used, at least when asymptotic techniques [14] cannot be applied, as in the case of coupling micropolar and classical continua.The influence of the selection of the UC is analyzed and some key issues are outlined.By considering two different UCs, selected for representing the composite texture, it emerges that the constitutive response of the homogenized medium depends on the choice of the cell.Indeed, while the elastic Cauchy coefficients are irrespective of the UC selected, this does not occur for the bending and skew-symmetric shear Cosserat coefficients, at least with regard to computational homogenization.This fact is also confirmed by the results obtained from the structural applications.Two numerical tests are presented to highlight the main aspects of the presented micropolar computational homogenization technique and to emphasize the differences obtained using the three procedures to describe the perturbation fields.

MICROPOLAR HOMOGENIZATION
he computational homogenization technique used here adopts a 2D micropolar continuum at the macro-level and a classical 2D Cauchy continuum at the micro-level.At typical macro-level material point, the displacement vector is defined, where 1 U and 2 U are the translational degrees of freedom and  is the rotational degree.The micropolar strain vector is characterized by six components as follows: where 1 E , 2 E and 12  are the axial and symmetric shear strains, 1 K and 2 K are the curvature components, while  is the skew-symmetric shear component.The compatibility equations can be written in compact form as: where the compatibility operator is defined as According to the strain driven approach, the macroscopic strain components, evaluated at , X are used as input variables for the microscopic level.The kinematic map, expressed in function of the vector , E is imposed on the UC, properly defining a BVP.At the micro-level a repetitive rectangular UC is selected, whose size is and its centre is located at the macroscopic point , X characterized by the displacement field of the UC domain  .The displacement field, resulting in the UC after solving the BVP, can be represented as the superposition of the assigned field    u X,x and a perturbation field    : The strain vector at the microscopic level is derived by applying the kinematic operator defined for the 2D Cauchy problem and, in expanded form, it results as: , with and = 0 with ,i  indicating the partial derivative with respect to .i x According to Eq. ( 5), the strain can be written as: The third order polynomial map, proposed in [5 ,6] and modified in [7], is used.Different material symmetries can be considered ranging between the isotropic and the orthotropic case.In the considered orthotropic case, the kinematic map can be written in compact form as: where 1 2 1 = / e e  , 1 e and 2 e are the values of Young's modulus of the equivalent homogenized orthotropic material, 12  the Poisson ratio and: The dependence of the kinematic map on the effective elastic coefficients is a consequence of the enforcement of the balance equations in the UC, ensuring that the kinematic map is the solution of the BVP in the UC made of the homogenized material.
The quantities 1  , 2  and 3  account for the effects of the perturbation parts of the displacement field in determining the average macroscopic strains and are presented in detail in [11].
In this work linear elastic isotropic behavior is assumed for the constituents at the micro-level.A straightforward extension to the case of non-linear material behavior can be, however, easily made.

HETEROGENEOUS MATERIAL: PERTURBATION FIELD IN THE UC DOMAIN
he characterization of the perturbation field    u X,x , arising in the UC when a heterogeneous medium is taken into account, is discussed.Depending on the choice of the boundary conditions imposed at the micro-level, very different results can be obtained in terms of displacement fields solution of the BVP.
It is well established that, in the case of first order homogenization,    u X,x is a periodic field.In this instance, the Periodic Boundary Conditions (PBCs) are suitable to correctly reproduce the unknown perturbation field.The extension to the case in which higher order polynomial terms are considered in the kinematic map is not trivial and it is not possible to a-priori assume the periodicity of    u X,x , as remarked in [7,9,11].In this section, three approaches are introduced to characterize the perturbation field.
In the first approach, the classical periodicity conditions are considered.The second technique assumes the decomposition of the perturbation field in different contributions, related to the first, second and third order gradients of the kinematic map as proposed in [11].The last approach is based on the enforcement of proper boundary conditions (BCs) on the UC, as described in [7], resulting from the analysis of the actual perturbation field distribution in the RVE undergoing remote fully displacement BCs.Finally, the comparison between the numerical results obtained using the adopted procedures for a paradigmatic example of a two-phase composite material, characterized by cubic symmetry, is presented.

Procedure A: Periodic BCS (PBCs) in the UC
The first procedure to characterize the field    u X,x involves the solution of the BVP in the UC under standard PBCs.Corresponding points on opposite sides of the UC are constrained to undergo the same perturbation displacement.The following conditions are imposed on the sides of the UC: where the dependence on X has been omitted for the sake of brevity.In the presence of nonvanishing components 1 K , 2 K and  , this assumption leads to unrealistic distributions of the perturbation field.

T Procedure B: 3 Step homogenization
This procedure has been proposed by the authors in [11] as an extension to the 2D micropolar computational homogenization of the homogenization of second gradient continua via the asymptotic approach [12,15].Here, the basic idea is recalled and the main steps are addressed.
It is assumed that the total perturbation displacement vector    u X,x can be expressed as the sum of three fields evaluated in sequence.Initially, only the first order terms of the kinematic map, multiplying the components 1 E , 2 E and 12  in (8), are activated; subsequently, the effects of the quadratic terms related to 1 K and 2 K are considered and, finally, the third order term associated with  is taken into account.Therefore, the vector results as: When only the linear terms of the kinematic map are considered, the case of the first order homogenization is recovered.
Here, it is assumed that   1 r X,x  is evaluated as the product of unknown functions times the independent components of the first gradient of the kinematic map, as: where   When the presence of the curvatures 1 K and 2 K is also considered, with = 0  , the perturbation field can be expressed as: where now     being periodic functions evaluated in the UC when 1 0 K  and 2 0 K  .Finally, when the component  is also taken into account, it results that: The field  

Procedure C: analysis of the perturbation field in the RVE
The characterization of the perturbation field  u is performed in [7] by evaluating its actual distribution.To this end, a RVE obtained as assemblage of a large number of UCs for a selected two-phase periodic composite material is considered and remote fully displacement BCs are prescribed.In particular, the micropolar deformation modes are imposed on the boundary, according to the kinematic map in Eq. ( 8), and the RVE response is evaluated by finite element method.Hence, the distribution of the perturbation field arising in the central UC of the RVE is taken as the benchmark.Thus, the suitable BCs to impose on the single UC are derived, in order to reproduce, with a satisfactory level of accuracy, the actual distribution of the perturbation field.Some selected two-phase composite materials, characterized by material symmetries, ranging from cubic to orthotropic, are analyzed.In all the cases considered, similar distributions of the perturbation displacement fields on the UC boundary emerge.Differently from the case of the first order homogenization procedure, where periodic BCs are suitably adopted, in the analyzed cases more complex BCs have to be considered, which are different for the two components of  u .In Fig. 1 the derived BCs are summarized.In the first row, the applied Cosserat macroscopic deformation components are reported; in the second row the BCs for the component  1 u along the horizontal and vertical edges of the UC are schematically reported, while in the third row those for the displacement component 2 u  are shown.The symbol "p" indicates periodic BCs , see Eq. 12; "s" refers to skew-periodic BCs, that is: Finally, "0" indicates zero perturbation displacement BCs.
Figure 1: Boundary conditions to impose on the UC to evaluate the perturbation fields, derived in [7].

Comparison between Procedures A, B and C
The approaches presented were compared by carrying out some numerical tests.This subsection is devoted to a qualitative discussion of the results obtained that are reported in detail in [11].
Reference is made to a specific two-phase composite material, characterized by cubic symmetry, whose texture is made from a soft matrix with stiff square inclusions, both isotropic and regularly spaced.The volume fraction , f defined as the ratio between the area of the inclusions and the total area of the UC, is set equal to 36%.As reference solution the distribution of the perturbation field arising in the central UC of a sufficiently large RVE, undergoing remote fully displacement boundary conditions, is considered.The three procedures lead to the same correct results, if the classical macroscopic strain components 1 , E 2 E and 12  are activated.For the components 1 K , 2 K and ,  different considerations apply.For both components 1 K and  Procedure B provides the best estimate of the perturbation fields.In the case of 1 K only the vertical component of the perturbation field, evaluated along the horizontal edges of the UC, slightly differs from the referential solution, while in the case of  a very good approximation is guaranteed anywhere on the edges.
To be noted is that the solution corresponding to the curvature 2 K can be obtained by rotating that evaluated for 1 K .
Adopting Procedure A, the field  u differs qualitatively from the actual solution and it can be concluded that this procedure leads to grossly erroneous results.
Finally, the results obtained by applying Procedure C are close to those evaluated with Procedure B, although the vertical perturbation components along the horizontal lines in presence of 1 K and  are worse approximated.

IDENTIFICATION PROCEDURE: HILL-MANDEL MACRO-HOMOGENEITY CONDITION
he identification procedure adopted in this study is based on the generalized Hill-Mandel macro-homogeneity condition.The virtual work evaluated at the macroscopic Cosserat point is set equal to the average virtual work of the heterogeneous Cauchy medium in the UC.Thus, the following expression holds: where Σ is the micropolar stress vector evaluated at the macroscopic point, while σ is the Cauchy stress vector at the typical point of the UC.After solving the BVP on the UC and determining the microscopic stress and strain fields, σ and ε , the homogenized Cosserat elastic constitutive matrix C can be derived by using Eq. ( 22).Considering a two-phase composite material with a regular arrangement of the inclusions, characterized by orthotropic texture, the homogenized Cosserat elastic constitutive matrix, expressed in a reference frame aligned with the principal axes of the material, results as:  The results, presented in [11], are critically discussed to highlight the dependence of the identification equivalent coefficients on the centering of the UC.Indeed, differently from the Cauchy coefficients, which are proved to be irrespective of the choice of UC, for the bending and skew-symmetric shear micropolar coefficients this is not straightforward, at least in the framework of computational homogenization.
In the following the focus is on the determination of 44 C , 55 C and 66 C , governing the bending and skew-symmetric shear behavior of the micropolar equivalent medium.

T Aiming at estimating the coefficient 44
C , the macroscopic strain component 1 1 K  is applied to the UC, with all the other components set equal to zero.Note that, due to the cubic symmetry of the composite material, the case 2 1 K  leads to results which are the rotated results of the case 1 1 K  , so that 55 44 .

C C 
For the two selected UCs, Procedures A, B and C are applied.It emerges that Procedures B and C lead to homogenized constitutive parameters that differ by about 8% for the same UC, while Procedure A provides results that differ by an order of magnitude.Moreover, the obtained results depend on the choice of the UC.Indeed, the values of 44 C computed for the two UCs differ by 14%, for Procedure B and about 27% for Procedure C.
These results can be explained by taking into account the explicit expression of the internal work at the right-hand side of Eq. ( 22).It emerges, in fact, that the relative position of the two (stiffer and softer) constituents, with respect to the UC center, strongly influences the value of the identified parameter.This problem is well known [2,8,11] and is related to the definition of the higher-order or couple stresses as the volume average of the product of microscopic stresses and microscopic coordinates over the UC.Similar considerations apply when the component 1   is considered, while all the other macro-level strain components are set equal to zero.Again, different results are obtained for the two UCs and for the adopted methods.Further numerical tests are performed to investigate on the influence of the size of the RVE [9].Various square RVEs are considered, taking into account assemblages of 3 3  , 5 5  , 7 7  , 9 9  ,… 15 15  UCs, subjected to the BCs shown in Fig. 1, which correspond to the Procedure C. The average internal work is evaluated over the entire RVE domains.It emerges that, by introducing proper scaling factor depending on the size L of the square RVE, in both cases of 44 C and 66 C , as the RVE size increases, the different UCs converge to the same quantity, from above and from below, respectively.

Rectangular wall under vertical loading
rectangular wall made from a periodic composite material is studied.In Fig. 3 the geometry is reported together with loading and boundary conditions.For each case, it is assumed that the heterogeneous material is obtained by adopting the UC1 in Fig. 2-a or UC2 in Fig. 2-b.In Fig. 3 the arrangement considered for the case of 16 12  UCs, by adopting both UC1 (left side) and UC2 (right side), is shown.The numerical simulations are performed considering the response of the two heterogeneous materials characterized by the arrangements shown in Fig. 3, compared with the response of the homogenized micropolar media.The structural A stiffness is evaluated by dividing the total base vertical reaction by the maximum vertical displacement measured at the midspan of the top edge.In Tab. 1 the values of the stiffness, obtained considering UC1 for the three assemblages ( 16 12  , 8 6  and 4 3  ), are shown.In the first column the results obtained with the micromechanical model and representing the reference solution are reported.In the second column the stiffness values, computed by adopting a standard first order computational homogenization (Cauchy), are shown, normalized with respect to the reference solution. UCs and by about 29% for 4 3  UCs.In Tab. 2 the same results reported in Tab. 1 are shown, when UC2 is taken into account.Also in this case, the results obtained via Procedure B are in very good agreement with the micromechanical model.Considering the results in both Tab. 1 and 2, it emerges that the micropolar effects become more evident as the number of UCs decreases, since the ratio between the microstructural size, directly related to the dimension of the inclusions, and a typical structural dimension, increases.As expected, the Cauchy model gives the same results in all the considered cases and is suitable to correctly estimate the structural response as the above ratio increases.Moreover, it is noteworthy that the position of the inclusions in the heterogeneous medium significantly influences the response.Indeed, especially for the 4 3  UCs significantly different results are obtained by adopting the two arrangements.Finally, to investigate the influence of the aspect ratio on the global elastic response, three different geometries are considered, corresponding to B/H 1.3  , B/H 1.6  and B/H 2  .The first geometry is the same presented above corresponding to 16 12  UCs.In Tab. 3 and 4 the results of the structural stiffness (adopting the same normalization as before), as the aspect ratio changes, are reported for UC1 and UC2, respectively.It can be remarked that Procedure B confirms to be the most suitable to reproduce the actual behavior of the heterogeneous medium.Moreover, the slenderer the wall ( B/H 1.3  ) the closer the values obtained with the first order technique are to the response of the micromechanichal models.The micropolar effects are, thus, more pronounced for the squat wall ( B/H 2  ).Also in this case, Procedure B gives results which best match those obtained with the micromechanical model.

Simple shear test of a composite strip
As a second example, a displacement driven shear test on a 2D structure, made from the same composite medium considered in the previous example, is presented.A strip with B/H 0.1  is considered assuming plane strain condition.In Fig. 4 the schematic of the geometry and boundary conditions is shown for two possible strips, corresponding to assemblages of UC1 and UC2 in Fig. 2-a and e 2-b, respectively.The strips are fixed at both bottom and top edges and a horizontal displacement d=H/100 is prescribed at the top.In Fig. 5-a the displacement horizontal component, evaluated along the vertical symmetry axis of each strip along the vertical axis of the strip is shown.The responses of the micromechanical models are represented in solid line (heter1) and in dash-dot line (heter2) for the strips reported in Fig. 4-a and 4-b, respectively.The homogenized constitutive coefficients adopted in the micropolar models are evaluated using Procedure B. The squares represent the response of the homogenized micropolar model adopting the UC1 in Fig. 2-a (CosB_UC1) and the circles represent the response with UC2 in Fig. 2-b (CosB_UC2).Finally, the dotted line refers to the response obtained by the standard first order computational homogenization (Cauchy), able to reproduce a linear variation of the horizontal displacement component along the height of the strip.No relevant differences between the results obtained with the two micromechanical models arise and both the homogenized micropolar models can satisfactorily follow the expected displacement distribution.In Fig. 5-b the rotation evaluated along the vertical symmetry axis of each strip versus the vertical axis is reported.Owing to the symmetry of this displacement component with respect to a horizontal axis located at the mid height, only one half of the strip height is depicted.Due to the different natures of the compared models, the rigid rotation  , i.e. the skew- symmetric part of the displacement gradient, is reported for micromechanical and Cauchy models, while the Cosserat rotation  is shown for the micropolar homogenized model.Line styles and symbols have the same meaning as in Fig. 5

CONCLUSIONS
he micropolar computational homogenization are considered.The perturbation fields in the presence of higher order polynomial boundary conditions is analyzed adopting three different procedures, which lead to different results.The first method is based on the imposition of periodic boundary conditions, classically adopted in the standard first order homogenization, also in the presence of higher order terms of the polynomial map and provides qualitatively incorrect results.The second procedure ( 3Step Homogenization) adopted, although the most complex, produces the best results.Finally, the methodology characterizing the perturbation fields on the basis of the direct observation of the large heterogeneous medium behavior gives results slightly differing from the 3 Step Homogenization, but is much simpler.Thus, this appears as the best compromise between accuracy and efficiency in correctly reproducing the actual trends of perturbation fields, when a single UC is analyzed.Moreover, the identification of the homogenized linear elastic constitutive parameters adopting the classical Hill-Mandel procedure is addressed.Considering different UCs, referred to the same composite material, it emerged that the constitutive response of the homogenized medium depends on the choice of the cell, adopting all the three procedures exploited for reproducing the perturbation fields.Indeed, while the elastic Cauchy coefficients are irrespective of the UC selected, for the bending and skew-symmetric shear micropolar coefficients this does not occur, at least not in the framework of computational homogenization.

T
Two structural examples are presented.The aim is to discuss the key aspects of the adopted micropolar formulation and to stress the differences obtained using the three presented procedures in terms of the global response.It would appear relevant to carry out further developments, focusing on the formulation of the kinematic map linking the macro-and micro-levels, as well as on the improvement of the identification procedure.The purpose is to better clarify some ongoing issues and to solve inherent limitations of the procedure, as for example the contradictory result consisting in nonvanishing micropolar effects in the presence of homogeneous materials.


the ratio between the dimensions of the UC and

2 r 1 r 2 γ 1 γ
X,x is the only unknown field.Following the same procedure as for the first term   is expressed as the product of unknown functions and   X, x , i.e. the nonvanishing components of the first gradient of   X, x .Then, the unknown field   2 r X,x  is represented in the form:

2 γ
the product of unknown functions times the relevant and nonvanishing components of the first gradient of   X, x , and it results:

6  and 4 3 
First, the capability of the homogenized Cosserat model to account for size effects is analyzed.The following dimensionless parameters B/H 1.3  , b/H 0cases are analyzed, considering the wall made from the periodic repetition of UCs: 16 12  , 8 UCs.

Figure 3 :
Figure 3: Schematic of the rectangular wall: geometry, loading and boundary conditions.
Finally, in the last three columns, Cos A, Cos B and Cos C refer to the responses obtained using the micropolar homogenized model whose homogenized elastic coefficients 44 C , 55 C and 66 C are derived by means of Procedure A, Procedure B and Procedure C, respectively.These are also normalized with respect to the reference solution.The comparison of the values collected in Tab. 1 highlights that Cos B provides the best estimation of the structural stiffness, while Cos A gives a response overestimating by out 23% the actual stiffness for 16 12

Table 1 :
Structural stiffness in the case of different assemblages of UC1: Heter = micromechanical model; Cauchy= homogenized first order model; Cos A = homogenized Cosserat with Procedure A; Cos B = homogenized Cosserat with Procedure B; Cos C = homogenized Cosserat with Procedure C.

Table 2 :
Structural stiffness in the case of different assemblages of UC2: Heter = micromechanical model; Cauchy= homogenized first order model; Cos A = homogenized Cosserat with Procedure A; Cos B = homogenized Cosserat with Procedure B; Cos C = homogenized Cosserat with Procedure C.

Table 3 :
Structural stiffness for different ratios of B/H in the case of UC1: Heter = micromechanical model; Cauchy= homogenized first order model; Cos A = homogenized Cosserat with Procedure A; Cos B = homogenized Cosserat with Procedure B; Cos C = homogenized Cosserat with Procedure C.

Table 4 Structural
stiffness for different ratios of B/H in the case of UC2: Heter = micromechanical model; Cauchy= homogenized first order model; Cos A = homogenized Cosserat with Procedure A; Cos B = homogenized Cosserat with Procedure B; Cos C = homogenized Cosserat with Procedure C.