Size effects of mechanical metamaterials: a computational study based on a second-order asymptotic homogenization method

In this paper, size effects exhibited by mechanical metamaterials have been studied. When the sizescale of the metamaterials is reduced, stiffening or softening responses are observed in experiments. In order to capture both the stiffening and softening size effects fully, a second-order asymptotic homogenization method based on strain gradient theory is used. By this method, the metamaterials are homogenized and become effective strain gradient continua. The effective metamaterial parameters including the classical and strain gradient stiffness tensors are calculated. Comparisons between a detailed finite element analysis and the effective strain gradient continua model have been made for metamaterials under different boundary conditions, different aspect ratios, different unit cells (closed or open cells) and different topologies. It shows that both stiffening and softening size effects can be captured by using the effective strain gradient continua models.


Introduction
Mechanical metamaterials are of growing interest due to their extraordinary mechanical properties for engineering applications [13,16,30,51]. Mechanical metamaterials are designed with well-defined microstructures in order to achieve controllable properties. Consequently, it is not only the composition of materials but also the microstructures that govern the behavior of metamaterial [39]. Most of the existing designs of mechanical metamaterials are achieved either by utilizing origami and kirigami deformation mechanisms [17,52] or by tailoring the topology of the lattice cells and by arranging them periodically at the microscale [64]. For example, so-called pantographic microstructures [25][26][27]53,56,65] are constituted by two parallel arrays of beams connected by small internal cylinders (pivots) and display novel functionalities. Their mechanical properties are reported to be usually size dependent [4]. In other words, remarkable size effects are observed when metamaterials are tested. For example, in the experiments of [4,45] the samples become stiffer with decreasing sample size. But smaller is not always stiffer. For example, a softening effect is observed in the studies of [59].
In order to capture the size-dependent behaviors completely, a direct finite element computation with a detailed mesh and an appropriate constitutive relation connecting stress and strains is always possible [62] and usually provides good accuracy. However, due to the complexity of the microstructures, it carries too much computational burden even for modern computer technology. The homogenization method [5,12,[20][21][22]32,34,36,37,42,63] is thus used as an alternative way to model metamaterials as an equivalent effective H. Yang (B) · W. H. Müller Institute of Mechanics, Chair of Continuum Mechanics and Constitutive Theory, Technische Universität Berlin, Einsteinufer 5, 10587 Berlin, Germany E-mail: hua.yang@campus.tu-berlin.de homogeneous continuum. The homogenization modeling method is also motivated by the desire for an increasing understanding of the micromechanical mechanism behind the peculiar behaviors of metamaterials [28]. Despite the efficiency of the homogenization method when applied in practice, it shows a limitation. Metamaterials were homogenized into a Cauchy continuum in past studies [35]. However, the size-dependent behaviors could not be captured by such method. This is because the classical Cauchy theory lacks an internal length scale of the microstructure. Hence, it is unable to predict the size effects.
Consequently, homogenization in the framework of generalized mechanics such as micropolar theory [41] and strain gradient theory [15,19,47] was employed. The problem of getting homogenized constitutive equations for generalized continua is rather difficult, and many debates are active in the literature. Kumar and McDowell [41] derived the effective material parameters based on an energy equivalence assumption by means of micropolar elasticity. Dos Reis and Ganghoffer [28] used the asymptotic method and homogenized metamaterials to micropolar media. An effective couple-stress continuum model was built by Liu and Su [46]. The effective moduli are obtained based on a selected representative volume element (RVE) under certain boundary conditions. Many efforts have been made to develop strain gradient continua. Dell'Isola and colleagues developed second gradient homogenization techniques for the so-called pantographic structures which has been extensively studied, for example, in [3,9,11,[23][24][25][26][27]55]. Boutin et al. [18] obtained high-order continuum model by means of a micro/macro-identification procedure based on the asymptotic homogenization method of discrete media. Misra and Poorsolhjouy [48] identified higher-order elastic constants for grain assemblies based on granular micromechanics. Abdoul-Anziz and Seppecher [2] determined the effective behavior of 1D, 2D and 3D periodic structures through a Γ -convergence theorem. By means of homogenization, Li [43] established a strain gradient constitutive law. In [8,44], a fast Fourier transform (FFT) and an FEM approach have been used to deliver the effective classical stiffness tensor as well as strain gradient stiffness tensor. In [67], the effective metamaterial parameters are derived based on quadratic boundary conditions. Effective moduli were presented for metamaterials of different symmetries. In [61], based on an assumption of strain energy equivalence at micro-and macroscale of an RVE, a computational solution was provided for determining parameters of metamaterials. It is also shown that the strain gradient stiffness tensor determined in [61] vanishes when the structure is purely homogeneous, and it is insensitive to a repetition of the basic cells. Many authors compare direct FEM computations to equivalent generalized homogeneous continua. For example, in a recent study [66] size effects of lattice structures were examined by comparing direct FEM computations with an equivalent micropolar elastic model. It was found that size effects are not captured in the micropolar model.
The aim of the present paper is to provide quantitative comparative studies for metamaterials constructed by repetitive planar lattice structures between direct FEM computations, which are understood as correct solutions, and equivalent effective strain gradient continua. These simulations will examine size effects directly for metamaterials under different boundary conditions, of different aspect ratios, constructed by different unit cells and by different topologies. A second-order asymptotic analysis (presented in Sect. 2) is used to model the metamaterials by equivalent effective strain gradient continua. In this paper, it is intended to investigate the mechanism of size effects and to answer the following question: Do the metamaterial parameters derived by using the second-order asymptotic analysis allow to capture size effects accurately including both stiffening and softening effects?

Second-order asymptotic homogenization method based on strain gradient theory
Consider a two-dimensional periodic metamaterial Ω, which is constructed by a finite number of representative volume elements (RVEs) Ω m P at the microscale, namely: We remark that the RVE at the microscale Ω m P represents the detailed microstructure, such as the fibers and the matrices in a composite material. The same RVE at the macroscale Ω M P is modeled homogeneously. The equivalence of the strain energies for an RVE at the micro-Ω m P and at the macroscale Ω M P is postulated. The definitions of the micro-("m") and of the macroscale ("M") energy density, w m and w M , respectively, read where first-order theory is used at the microscale and second-order theory at the macroscale. C m i jkl is defined in every material point of the heterogeneous structure. C M i jkl and D M i jklmn are the homogenized classical stiffness tensor and strain gradient stiffness tensor of the corresponding equivalent continua. The use of a displacement gradient in Eq. (2) will bring simplicity. Indeed, it can be replaced by the linear strain tensor. A discussion of objectivity of the strain energy can be found in [61]. In what follows, two different cases are demonstrated: a heterogeneous case at the microscale with known material properties versus a homogeneous case at the macroscale with sought effective parameters.

The macroscale energy for an RVE
Consider the macroscale case for an RVE Ω M P . As it is customary in spatial averaging, we define the geometric center c X of the RVE. Then, we approximate the macroscale displacement by a Taylor expansion around the value at the geometric center by truncating after quadratic terms (in order to account for the strain gradient effect) and then calculate displacement gradients of this approximation: where according to [61] The macroscale energy of an RVE reads as follows (the detailed derivation can be found in [61]), as the macroscale stiffness tensors are constant in space, Therefore, the macroscale energy of an RVE is expressed in terms of the gradient of macroscopic deformation.
In what follows, it will be presented by using asymptotic homogenization method that the microscale energy can be expressed in terms of the gradient of macroscopic deformation. This will also lead to connections between the parameters.

The microscale energy for an RVE
Following the asymptotic homogenization method presented in [61], we reformulate the left-hand side of Eq.
(2). The asymptotic homogenization method decomposes length scales by making use of global coordinates, X, depicting the global variations, and by using local coordinates, y, presenting the local fluctuation. The local coordinates are introduced as: where is a homothetic ratio connecting global and local coordinates. Following [61], the displacement field of an RVE, Ω m P , at global coordinates, X, is approximated by an asymptotic expansion with regard to , Please note that the governing equation within the RVE must be fulfilled where f i is body forces. Substituting Eq. (8) to Eq. (9) and gathering terms having the same order in lead to the following equations: -in the order of 0 Indeed, the solution in the order of −2 can be given as 0 u i (X) depends only on the macroscopic coordinates. It is assumed to be the known macroscopic displacement Likewise, the solutions for the order of −1 and 0 can be obtained [61]. Therefore, the microscale displacement field can be rewritten as: where ϕ abi ( y) and ψ abci ( y) are assumed to be y-periodic. They are the solutions of the following two partial differential equations [61] We refer to [61] for a derivation of Eq. (15) and Eq. (16). By using Eq. (14) and the latter on the left-hand side of Eq. (2), the microscale energy becomes withC Fig. 1 The lattice structure and its equivalent homogenized continuum Note that the rank 5 tensor vanishes for centrosymmetric materials,Ḡ = 0.

The equivalence of the macro-and microscale strain energy for an RVE
Based on Eq. (2), the effective parameters including the classical Cauchy stiffness tensor as well as the strain gradient ones can be calculated by: after computing ϕ and ψ from Eqs. (15) and (16) in an RVE. Indeed, Eqs. (15) and (16) are solved numerically by means of the FEM [61].

Computational study
Size effects of metamaterials constructed by lattice-based microstructures are computationally examined. A quantitative assessment of the capability to capture size effects by means of the second-order homogenization method is conducted in this section. Three types of computations as indicated in Fig. 1 are conducted in this section. The first simulation is a direct computation with the detailed mesh of the metamaterials. The metamaterials are "made" from an isotropic materials with Young's modulus E = 100MPa and Poisson's ratio ν = 0.3. The second and third ones are the corresponding homogenized classical continua as well as homogenized strain gradient continua with the same sizes, boundary conditions as in the first simulations. The stiffness tensor C M i jkl is needed for the second simulation. Both C M i jkl and D M i jklmn are required for the third one. These two tensors are delivered by using the method that was outlined above.
The first and second simulations are conducted by means of standard FEM. A discretization with a quadratic Lagrange shape function was carried out. The implementation of strain gradient elasticity was achieved based on isogeometric analysis (IGA) in order to ensure an at least C 1 continuity [38,40,54]. A verification of the used IGA codes with analytical solutions can be found in [60]. All of the computations are done in the The relative errors are defined in Eq. (22). They serve for evaluation of the accuracy of different approaches: Figure 1 presents an example of the square lattice and its homogenized continua. In order to understand the mechanism of size effects better, a number of factors (boundary conditions, aspect ratios of lattice structures, different unit cells, different microstructural topologies) are investigated in what follows.

Examination of size effects for different boundary conditions
Square lattice structures (width = length = 4 mm) as shown in Fig. 2 are investigated here for different boundary conditions. The left sides of the structures are fixed, and on the right sides, prescribed displacements are imposed, resulting in extension, shear and torsion tests. The material parameters are displayed in Table 1, where V f denotes the volume ratio (the volume of ligaments over the whole cell), t is the thickness of the ligament, and l is the length of the cell. The square lattice structure is of cubic symmetry. There are 3 and 6 independent parameters in the classical stiffness tensor and strain gradient stiffness tensor [6,7], respectively, as shown in the following: The numerical results are compiled in Figs. 3 and 4. The following observations and interpretations can be made:     -With increasing aspect ratio, the microstructural effects become more important and the size effects become more pronounced. -Strain gradient continua show good agreement with direct computation, and the size effects can be grasped well even for the large aspect ratio value of 50.

Examination of size effects for different topologies of microstructures
The three topologies of lattices shown in Fig. 10, which have been intensively studied in [41,46,66], are used. The volume ratio of the lattices is set to be equal leading to a difference of the thickness of ligaments. The material parameters are shown in Table 3. Two kinds of numerical experiments (either applying prescribed displacement or tractions) were performed.

Numerical tests under prescribed displacements
Torsion tests for 4×4 lattice structures under prescribed displacements (torsion boundary conditions) as shown in Fig. 10 were conducted.
From the numerical experiments, we see in Figs. 11 and12: -Stiffening size effects are found for all three lattices. Strain gradient continua are able to grasp the stiffening effects accurately in all cases. -The stiffening size effects gradually decay from the square lattice, the mixed triangle A lattice, to the mixed triangle B lattice.

Numerical tests for prescribed tractions
As displayed in Fig. 13, lattices (50 mm × 1 mm) with different microstructural topologies were investigated. The left sides are clamped, and to the right sides, tractions are applied. The results for total displacements, total strain energies and the relative energies can be found in Figs. 14, 15 and 16. The following observations can be made: -The direct computations for lattices show smaller total displacements and larger total strain energies than the Cauchy continuum. This stiffening size effect is observed for all three lattices. -A good match between strain gradient continua and direct computations is visible for both total displacements and total strain energies. Strain gradient continua are able to grasp stiffening effects accurately in all of the cases. -The topologies of the microstructures have a big influence on the size effects.

Discussions
The results show that equivalent homogeneous strain gradient continua can be used as a replacement for metamaterials with lattice microstructures. Stiffening or softening size effects of metamaterials are both captured by strain gradient continua. The reason for the stiffening size effects is mainly due to -Beam bending effects From the extension and shear tests conducted in Sect. 3.1, it is concluded that the bending energies are negligible. Hence, no size effects are detected. In the case of torsion tests, beam bending energies cannot be ignored, resulting in a stiffening size effect. This size effect grows with increasing aspect ratios of the lattices as shown in Sect. 3.2. Size effects shown in different topologies of microstructures are originated from beam bending, as well. As shown in Sect. 3.4, the size effects are more distinct in square lattice than in mixed triangle A and B lattices due to the fact that the square lattice is a bending dominated lattice and the others are stretch dominated [5]. The reason for the softening size effects is mainly due to -Zigzagging edge arrangements contributing less stiffness than the interior.
The Cauchy theory does not have the capability to take the microstructural information into account. However, the differences of geometries and topologies of microstructures and edge arrangement can be described by strain gradient theory with additional constitutive parameters.
Another issue is the positive definiteness of the strain energy for the homogenized strain gradient continua since some negative values are observed in the strain gradient stiffness tensors. Indeed, in this paper, the assumption of equivalence of strain energies for one RVE is used. This means at the macroscale that the homogenized strain energies of the continua are equal to that at the microscale, where the constituents of the RVE are linear elastic ensuring a positive strain energy. For the FEM implementation, if the size of the FEM element is larger than or equal to that of the RVE, an outcome of positive energy for the FEM computation is ensured. This was demonstrated also in [8,41,44]. Fortunately, in the present work, isogeometric analysis was used, and it is possible to specify relatively large sized elements: So convergence of computations is assured by setting a higher polynomial degree of NURBS. Note that choosing a larger size of elements than RVE in FEM is a sufficient but not a necessary condition to ensure the positive definite of the effective material parameters.

Conclusions
A second-order asymptotic homogenization was used to model metamaterials of lattice-based microstructures as effective homogeneous strain gradient continua. Size effects of metamaterials under different boundary conditions, different aspect ratios, different unit cells (closed or open cells), different topologies were examined. It was shown that both stiffening and softening size effects can be fully captured by effective strain gradient continua. Stiffening size effects are attributed to beam bending, and softening size effects are found due to the zigzag edge arrangements. It should be remarked that the presented homogenized model can produce the symmetry properties considered in [6,7,29,50].
Finally, it is possible to envision for the near future the following extensions of this work: -The extension to 3D cases will be conducted in order to be able to predict out-of-plane deformations [11,23,33] and to study 3D metamaterials. -The application of the homogenization method to the analysis of pantographic (bi-pantographic) structures will be explored [9,10,26,27]. -A comparative study will be undertaken between the presented computational homogenization method and that shown in [23,31]. Indeed, in [31], the constitutive parameters of macro-homogenized model are identified by a fitting procedure to get the best possible agreement with micromodel. The homogenization is made "globally" and conjecturing the macro-continuum model "a priori." On the other hand, in the current work, the structure of the homogenized description and its parameters are explicitly deduced from an RVE. It will be interesting to investigate differences of the identified parameters by means of these two methods.  -Hencky-type discrete models are the fundamental building block for continuum models [57,58]. Physical systems studied in the present paper can be also modeled by the so-called Hencky-type models. Future investigations will also address this point. -Comparisons with the results presented in [2,18,29,48] will be beneficial to highlight the physical meaning of the identified higher-order constitutive parameters.