Determination of Metamaterial Parameters by Means of a Homogenization Approach Based on Asymptotic Analysis

TL;DR

This paper presents a homogenization method based on asymptotic analysis to determine effective parameters of metamaterials, linking micro- and macro-scale properties through energy equivalence and finite element analysis.

Contribution

It introduces a homogenization approach that accurately derives macro-scale parameters from micro-scale models, ensuring consistency with classical schemes and size independence.

Findings

Successfully determines macro parameters using FEM boundary value problems.

Ensures additional parameters vanish for homogeneous materials.

Parameters depend on intrinsic structure size, not RVE choice.

Abstract

Owing to additive manufacturing techniques, a structure at millimeter length scale (macroscale) can be produced by using a lattice substructure at micrometer length scale (microscale). Such a system is called a metamaterial at the macroscale as the mechanical characteristics deviate from the characteristics at the microscale. As a remedy, metamaterial is modeled by using additional parameters; we intend to determine them. A homogenization approach based on the asymptotic analysis establishes a connection between these different characteristics at micro- and macroscales. A linear elastic first order theory at the microscale is related to a linear elastic second order theory at the macroscale. Relation for parameters at the macroscale is derived by using the equivalence of energy at macro- and microscales within a so-called Representative Volume Element (RVE). Determination of parameters are succeeded by solving a boundary value problem with the Finite Element Method (FEM). The proposed approach guarantees that the additional parameters vanish if the material is purely homogeneous, in other words, it is fully compatible with conventional homogenization schemes based on spatial averaging techniques. Moreover, the proposed approach is reliable as it ensures that such resolved additional parameters are not sensitive to choices of RVE consisting in the repetition of smaller RVEs but depend upon the intrinsic size of the structure.

Tables

Table 1. Table 1: Voigt ’s notation used for 2D strain tensors.

	$I$	1	2	3
	$ij$	11	22	12

Table 2. Table 2: Voigt ’s notation used for 2D strain-gradient tensors.

	$I$	1	2	3	4	5	6
	$ijk$	111	221	122	222	112	121

Table 3. Table 3: Material properties used in lattice structures.

Type	$E$ in MPa	$\nu$
Matrix	100.0	0.3
Inclusion	${10}^{-30}$	${10}^{-30}$

Table 4. Table 4: Parameters determined for the square lattice.

$C_{1111}$ in MPa	$C_{1122}$ in MPa	$C_{1212}$ in MPa
11.177	0.555	0.060
$D_{111111}$ in N	$D_{111221}$ in N	$D_{111122}$ in N
0.005379	0.042197	-0.047860
$D_{221221}inN$	$D_{221122}$ in N	$D_{122122}$ in N
1.597997	0.076341	0.033462

Table 5. Table 5: Identified non-zero strain gradient stiffness parameters for the square lattice structure in N.

Type	$D_{111111}$ in N	$D_{111221}$ in N	$D_{111122}$ in N	$D_{221221}$ in N	$D_{221122}$ in N	$D_{122122}$ in N
1$\times$1	0.005379	0.042197	-0.047860	1.597997	0.076341	0.033462
0.5$\times$0.5	0.001344	0.010549	-0.011965	0.399499	0.019085	0.008365
0.2$\times$0.2	0.000215	0.001688	-0.001914	0.063919	0.003054	0.001385

Equations

$$\cup {\Omega }^P=\Omega ,{\Omega }^P\cap {\Omega }^Q=\varnothing ,P,Q=1,2,3,...M,P\ne Q.$$

$$\begin{gathered}\int_{\Omega^{P}}{w^{\text{m}}}\,\mathrm{d}V=\int_{\Omega^{P}}{w^{\text{M}}}\,\mathrm{d}V\ ,\\ \int_{\Omega^{P}}\frac{1}{2}C^{\text{m}}_{ijkl}u^{\text{m}}_{i,j}u^{\text{m}}_{k,l}\,\mathrm{d}V=\int_{\Omega^{P}}\frac{1}{2}\big{(}C^{\text{M}}_{ijkl}u^{\text{M}}_{i,j}u_{k,l}^{\text{M}}+D_{ijklmn}^{\text{M}}u_{i,jk}^{\text{M}}u_{l,mn}^{\text{M}}\big{)}\,\mathrm{d}V\ ,\end{gathered}$$

$$\begin{array}{c}{u}_i^{\text{m}}=x_i^{\text{m}}-X_i,\\ u_i^{\text{M}}=x_i^{\text{M}}-X_i,\end{array}$$

$$\begin{array}{c}\stackrel{\text{c}}{X}=\frac{1}{V}{\int }_{\Omega }\mathbf{X}\mathrm{d}V,\end{array}$$

$$u_i^{\text{M}}(\mathbf{X})$$

$$u_{i,l}^{\text{M}}(\mathbf{X})$$

$$\displaystyle=u^{\text{M}}_{i,l}\Big{|}_{\overset{\text{c}}{\bm{X}}}+u^{\text{M}}_{i,lk}\Big{|}_{\overset{\text{c}}{\bm{X}}}(X_{k}-\overset{\text{c}}{X}_{k})\ ,$$

$$u_{i,lm}^{\text{M}}(\mathbf{X})$$

$$⟨u_{i,j}^{\text{M}}⟩$$

$$⟨u_{i,jk}^{\text{M}}⟩$$

$$\begin{gathered}\langle{u}_{i,j}^{\text{M}}\rangle=u^{\text{M}}_{i,j}\Big{|}_{\overset{\text{c}}{\bm{X}}}\ ,\ \ \langle{u}_{i,jk}^{\text{M}}\rangle=u^{\text{M}}_{i,jk}\Big{|}_{\overset{\text{c}}{\bm{X}}}\ .\end{gathered}$$

$$u_i^{\text{M}}(\mathbf{X})$$

$$u_{i,j}^{\text{M}}(\mathbf{X})$$

$$u_{i,jk}^{\text{M}}(\mathbf{X})$$

$$\begin{gathered}\int_{\Omega^{P}}\frac{1}{2}\big{(}C^{\text{M}}_{ijlm}u^{\text{M}}_{i,j}u^{\text{M}}_{l,m}+D^{\text{M}}_{ijklmn}u^{\text{M}}_{i,jk}u^{\text{M}}_{l,mn}\big{)}\,\mathrm{d}V=\frac{1}{2}C^{\text{M}}_{ijlm}\int_{\Omega^{P}}u^{\text{M}}_{i,j}u^{\text{M}}_{l,m}\,\mathrm{d}V+\\ \frac{1}{2}D^{\text{M}}_{ijklmn}\int_{\Omega^{P}}u^{\text{M}}_{i,jk}u^{\text{M}}_{l,mn}\,\mathrm{d}V=\frac{1}{2}C^{\text{M}}_{ijlm}\int_{\Omega^{P}}\Big{(}\langle{u}_{i,j}^{\text{M}}\rangle+\langle{u}_{i,jk}^{\text{M}}\rangle(X_{k}-\overset{\text{c}}{X}_{k})\Big{)}\times\\ \times\Big{(}\langle{u}_{l,m}^{\text{M}}\rangle+\langle{u}_{l,mn}^{\text{M}}\rangle(X_{n}-\overset{\text{c}}{X}_{n})\Big{)}\,\mathrm{d}V+\frac{1}{2}D^{\text{M}}_{ijklmn}\int_{\Omega^{P}}\langle{u}_{i,jk}^{\text{M}}\rangle\langle{u}_{l,mn}^{\text{M}}\rangle\,\mathrm{d}V=\\ =\frac{1}{2}V\Big{(}C^{\text{M}}_{ijlm}\langle{u}_{i,j}^{\text{M}}\rangle\langle{u}_{l,m}^{\text{M}}\rangle+(C^{\text{M}}_{ijlm}\bar{I}_{kn}+D^{\text{M}}_{ijklmn})\langle{u}_{i,jk}^{\text{M}}\rangle\langle{u}_{l,mn}^{\text{M}}\rangle\Big{)}\ ,\\ \end{gathered}$$

$$\begin{array}{c}{\bar{I}}_{kn}=\frac{1}{V}{\int }_{{\Omega }^P}(X_k-{\stackrel{\text{c}}{X}}_k)(X_n-{\stackrel{\text{c}}{X}}_n)\mathrm{d}V.\end{array}$$

$$y_j=\frac{1}{ϵ}(X_j-{\stackrel{\text{c}}{X}}_j),$$

$${\mathbf{u}}^{\text{m}}(\mathbf{X})=\stackrel{0}{u}(\mathbf{X},\mathbf{y})+ϵ\stackrel{1}{u}(\mathbf{X},\mathbf{y})+{ϵ}^2\stackrel{2}{u}(\mathbf{X},\mathbf{y})+\dots ,$$

$$\begin{gathered}\big{(}C_{ijkl}^{\text{m}}u^{\text{m}}_{k,l}\big{)}_{,j}+f_{i}=0\qquad\forall\bm{X}\in\Omega^{P}\ ,\end{gathered}$$

$$\begin{gathered}u^{\text{m}}_{i,j}=\Big{(}\overset{0}{u}_{i}(\bm{X},\bm{y})+\epsilon\overset{1}{u}_{i}(\bm{X},\bm{y})+\epsilon^{2}\overset{2}{u}_{i}(\bm{X},\bm{y})+\dots\Big{)}_{,j}=\\ =\overset{0}{u}_{i,j}+\frac{\partial\overset{0}{u}_{i}}{\partial y_{k}}\frac{\delta_{kj}}{\epsilon}+\epsilon\overset{1}{u}_{i,j}+\epsilon\frac{\partial\overset{1}{u}_{i}}{\partial y_{k}}\frac{\delta_{kj}}{\epsilon}+\epsilon^{2}\overset{2}{u}_{i,j}+\epsilon^{2}\frac{\partial\overset{2}{u}_{i}}{\partial y_{k}}\frac{\delta_{kj}}{\epsilon}+\dots{}\end{gathered}$$

$$\begin{gathered}\bigg{(}C_{ijkl}^{\text{m}}\Big{(}\overset{0}{u}_{k,l}+{\frac{1}{\epsilon}}\frac{\partial\overset{0}{u}_{k}}{\partial y_{l}}+\epsilon\overset{1}{u}_{k,l}+\frac{\partial\overset{1}{u}_{k}}{\partial y_{l}}+\epsilon^{2}\overset{2}{u}_{k,l}+\epsilon\frac{\partial\overset{2}{u}_{k}}{\partial y_{l}}\Big{)}\bigg{)}_{,j}+\\ +\frac{\partial}{\partial y_{j}}\bigg{(}C_{ijkl}^{\text{m}}\Big{(}{\frac{1}{\epsilon}}\overset{0}{u}_{k,l}+{\frac{1}{\epsilon^{2}}}\frac{\partial\overset{0}{u}_{k}}{\partial y_{l}}+\overset{1}{u}_{k,l}+{\frac{1}{\epsilon}}\frac{\partial\overset{1}{u}_{k}}{\partial y_{l}}+\epsilon\overset{2}{u}_{k,l}+\frac{\partial\overset{2}{u}_{k}}{\partial y_{l}}\Big{)}\bigg{)}+f_{i}=0\end{gathered}$$

$$\frac{\partial}{\partial y_{j}}\Big{(}C_{ijkl}^{\text{m}}\frac{\partial\overset{0}{u}_{k}}{\partial y_{l}}\Big{)}=0\ ;$$

$$\Big{(}C_{ijkl}^{\text{m}}\frac{\partial\overset{0}{u}_{k}}{\partial y_{l}}\Big{)}_{,j}+\frac{\partial}{\partial y_{j}}\big{(}C_{ijkl}^{\text{m}}\overset{0}{u}_{k,l}\big{)}+\frac{\partial}{\partial y_{j}}\Big{(}C_{ijkl}^{\text{m}}\frac{\partial\overset{1}{u}_{k}}{\partial y_{l}}\Big{)}=0\ ;$$

$$\begin{gathered}\big{(}C_{ijkl}^{\text{m}}\overset{0}{u}_{k,l}\big{)}_{,j}+\Big{(}C_{ijkl}^{\text{m}}\frac{\partial\overset{1}{u}_{k}}{\partial y_{l}}\Big{)}_{,j}+\frac{\partial}{\partial y_{j}}\big{(}C_{ijkl}^{\text{m}}\overset{1}{u}_{k,l}\big{)}+\frac{\partial}{\partial y_{j}}\Big{(}C_{ijkl}^{\text{m}}\frac{\partial\overset{2}{u}_{k}}{\partial y_{l}}\Big{)}+f_{i}=0\ .\end{gathered}$$

$$u_i^{\text{m}}(\mathbf{X},\mathbf{y})={\stackrel{0}{u}}_i(\mathbf{X})+ϵ{\phi }_{abi}(\mathbf{y}){\stackrel{0}{u}}_{a,b}(\mathbf{X})+{ϵ}^2{\psi }_{abci}(\mathbf{y}){\stackrel{0}{u}}_{a,bc}(\mathbf{X})+\dots .$$

$$\begin{gathered}\frac{\partial}{\partial y_{j}}\bigg{(}C_{ijkl}^{\text{m}}\Big{(}\frac{\partial\varphi_{abk}}{\partial y_{l}}+\delta_{ak}\delta_{bl}\Big{)}\bigg{)}=0\ ,\end{gathered}$$

$$\begin{gathered}\frac{\partial}{\partial y_{j}}\bigg{(}C_{ijkl}^{\text{m}}\Big{(}\frac{\partial\psi_{abck}}{\partial y_{l}}+\varphi_{abk}\delta_{lc}\Big{)}\bigg{)}+C_{ickl}^{\text{m}}\Big{(}\frac{\partial\varphi_{abk}}{\partial y_{l}}+\delta_{ka}\delta_{lb}\Big{)}-{C}_{icab}^{\text{M}}=0\ .\end{gathered}$$

$$u_i^{\text{m}}(\mathbf{X},\mathbf{y})=u_i(\mathbf{X}{)}^{\text{M}}+ϵ{\phi }_{abi}(\mathbf{y})u_{a,b}^{\text{M}}(\mathbf{X})+{ϵ}^2{\psi }_{abci}(\mathbf{y})u_{a,bc}^{\text{M}}(\mathbf{X})+\dots .$$

$$u_{i,j}^{\text{m}}$$

$$=u_{i,j}^{\text{M}}+\frac{\partial {\phi }_{abi}}{\partial y_j}{u}_{a,b}^{\text{M}}+ϵ{\phi }_{abi}{u}_{a,bj}^{\text{M}}+ϵ\frac{\partial {\psi }_{abci}}{\partial y_j}{u}_{a,bc}^{\text{M}}+{ϵ}^2{\psi }_{abci}{u}_{a,jbc}^{\text{M}}+\dots$$

$$u_{i,j}^{\text{m}}$$

$$\displaystyle=\Big{(}\delta_{ia}\delta_{jb}+\frac{\partial\varphi_{abi}}{\partial y_{j}}\Big{)}\Big{(}\langle{u}_{a,b}^{\text{M}}\rangle+\epsilon y_{c}\langle u_{a,bc}^{\text{M}}\rangle\Big{)}+\epsilon\langle u_{a,bc}^{\text{M}}\rangle\Big{(}\varphi_{abi}\delta_{jc}+\frac{\partial\psi_{abci}}{\partial y_{j}}\Big{)}+\dots$$

$$=L_{abij}⟨u_{a,b}^{\text{M}}⟩+ϵM_{abcij}⟨u_{a,bc}^{\text{M}}⟩+\dots$$

$$L_{abij}$$

$$M_{abcij}$$

Full text

Keywords: Metamaterial, Homogenization, Strain gradient theory, Elasticity, Asymptotic analysis

1 Introduction

Periodic lattice type structures involving large number of repetitive substructures continue to attract the interest of many researchers because of their fascinating properties like relatively low manufacturing costs, high specific stiffness, etc. [46, 47, 19, 62, 48, 64]. Mechanical respons of such a structure depends not only on the material, but also on the morphology of its substructure [61, 60]. Hence, metamaterial is used for such a substructure depending material.

In order to design and fabricate metamaterials for engineering applications, the accurate and efficient prediction of their mechanical performances is of importance [88, 27, 32, 92, 51]. Indeed, standard numerical techniques like the Finite Element Method (FEM) are utilized to achieve a modeling of such a structure with every detail of its substructure [94, 95]. However, it requires the mesh size to be at least one order smaller than the substructure geometric size leading to very high computational costs. Hence, homogenization techniques are developed to upscale the mechanical response at the microscale—the presence of the substructure leads to a composite material, which can be seen as a heterogeneous material—to the macroscale by defining an appropriate constitutive equation. Especially in composite materials, with fibers embedded in a matrix building a periodic substructure, micro- and macroscale behaviors are modeled by the same linear elastic model, also called a Cauchy continuum. The homogenization of such periodic structures toward an equivalent Cauchy continuum has been investigated thoroughly [45, 56, 20, 67, 96, 65].

Many approaches in the literature assume that there exists an Representative Volume Element (RVE) with periodic boundary conditions that precisely captures the deformation behavior of the whole geometry. Such an approach utilizes the energy equivalence of the RVE at both macroscale and microscale as also used in [52]. The effective properties of such homogenized continua is in good agreement with experiments [87] under the condition that $L\gg l$, where $L$ represents the macroscopic length scale, i.e., (mean value of ) the geometric dimensions of the whole structure, and $l$ represents the length scale of the microscale, namely, the geometric dimensions of the substructure. The quantity $l$ will be used as the “length scale” of a basic cell of the structure, as demonstrated in Fig. 1. Note that the concept of a basic cell is different from an RVE. It is evident that a basic cell can be regarded as an RVE, and stacking or gathering several basic cells can construct an RVE as well. Classical homogenization encounters limitations [58, 10] when $L$ is of a comparable order of scale with respect to $l$.

Size effects fail to be captured by a standard homogenization having the same order theory at both scales. A feasible approach is to use first order at the microscale and second order at the macroscale leading to additional parameters at the macroscale to be determined, we refer to various formulations of a second order theory in [63, 89, 6, 2, 3, 33, 37, 80, 78, 75, 85, 8, 74, 77, 81, 7, 1]. Higher order theories are called generalized mechanics and homogenization in the framework of generalized mechanics is a challenging task endeavored by many scientists, among others by [40, 71, 49, 31, 79, 55, 17]. Mostly, it is agreed that homogenization of an RVE by involving so-called higher gradient terms of the macroscopic field is a natural way to include size effect [14, 58, 41, 44, 43, 42]. Using gamma-convergence, homogenization results have been obtained in [5, 70, 4]. A remarkable class of structures being described at the macroscale by using a second gradient elasticity theory in pantographic structures [18, 86, 83], which have received a notable follow-up in the literature [30, 34, 76, 93, 91], also from a mathematically rigorous standpoint, regarding fundamental issues such as well-posedness [36].

A possibly promising homogenization technique is the asymptotic analysis, which has been used to obtain homogenized material parameters in [84]. This method decomposes variables to their global variations and local fluctuations. Such a decomposition is used to generate closed form equations to determine constitutive parameters as applied in one-dimensional problems, for example in the analysis of composites [21, 16], while 2D problems [13, 15, 22, 72, 29] have been investigated numerically. FEM is employed in [68] demonstrating that higher-order terms start dominating as the difference between parameters of composite materials increases. A second-order asymptotic and computational homogenization technique is proposed by [13] that solves the boundary value problems generated by the asymptotic homogenization with a quadratic ansatz. However, there are still two main issues which are not well discussed when trying to homogenize structures in the framework of generalized mechanics [90]:

• •

The first one concerns the compatibility such that parameters of the strain gradient stiffness tensor should vanish when the structure is purely homogeneous.

• •

The second one is about reliability such that the strain gradient stiffness tensor has to be insensitive to the repetition of the basic cell.

A successful attempt is made in [57, 58] establishing a connection between microscale parameters and macroscale parameters (by using the strain gradient theory) by proposing a “correction” term rendering the strain gradient stiffness tensor satisfying compatibility and reliability requirements. Different numerical solution methods are used for this approach, Fast Fourier Technique (FFT) is employed in [59] and FEM is exploited in [15]. We follow their methodology and propose an alternative derivation for this “correction” term in Section 2 and try to do it in a pedagogical way. Furthermore, we apply and validate the method for simple yet general 2D metamaterials in Sections 3 and 4 by using the FEM. In order to demonstrate its versatility, computations of the square lattice are performed in Section 5. The computations are performed with the aid of open-source codes developed by the FEniCS project [2]. The proposed method delivers all metamaterial parameters in 2D by using linear elastic material model at the microscale after a computational procedure as investigated in the following.

2 Connection of micro- and macroscale parameters

Consider a continuum body occupying a domain $\Omega$ in two-dimensional space, $\Omega\in\mathbb{R}^{2}$. The metamaterial embodies an RVE, $\Omega^{P}$, where periodically aligned RVEs constitute metamaterials domain,

[TABLE]

RVE at the microscale represents the detailed substructure like fibers and matrix in a composite material. The same RVE at the macroscale is modeled by a homogeneous metamaterial and we assume that their stored energy values are equivalent although the definitions at both scales differ. We use a first order theory for defining the energy (volume) density of an RVE at the microscale, $w^{\text{m}}$, whereas we utilize a second order theory at the macroscale for the energy density, $w^{\text{M}}$, leading to

[TABLE]

where and henceforth we apply Einstein summation convention over repeated indices and use a comma notation for space derivatives in $\bm{X}$. Moreover, all fields are expressed in Cartesian coordinates. The microscale stiffness tensor, $C^{\text{m}}_{ijkl}$, is a function in space. Consider a lattice substructure. Even if the trusses are of homogeneous material, voids between trusses generates a heterogeneous substructure at the microscale such that microscale stiffness tensor depends on space coordinates and possesses either the value of truss material or zero due to voids. In contrary, macrocale material tensors, $C^{\text{M}}_{ijkl}$ and $D^{\text{M}}_{ijklmn}$, are constant in space as they are generated by the homogenization procedure to be explained in the following. The continuum body at the reference frame has particles at coordinates $X_{i}$, where they move to $x_{i}$ under a mechanical loading. The displacement is the deviation from the reference frame and we emphasize that the microscale displacement, $u^{\text{m}}_{i}$, is different than the macroscale displacement, $u^{\text{M}}_{i}$,

[TABLE]

as the current positions of particles differ. This difference between ${x_{i}^{\text{m}}}$ and ${x_{i}^{\text{M}}}$ is illustrated in Fig. 1. For demonstrating the microscale deformation, the substructure is visualized as well. For simplicity, a well-known example is used, namely, composite materials with the red inclusion (fibers) embedded in the blue material (matrix). For the homogenized case, the same particle moves to $x_{i}^{\text{M}}$ expressed at the macro-scale without the substructure. We emphasize that micro- and macroscales are both expressed in the same coordinate system. Two different cases are examined, a heterogeneous case at the microscale with known material properties versus a homogeneous case at the macroscale with sought parameters. In order to identify the material parameters, the strain energy expressions for macro- and microscales are derived in what follows.

2.1 The macroscale energy for an RVE

Consider the macroscale case for an RVE, $\Omega^{P}$. As usual in spatial averaging, we define the geometric center $\overset{\text{c}}{\bm{X}}$ of the RVE,

[TABLE]

and 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 calculate displacement gradients of this approximation

[TABLE]

According to Eq. (5), spatial averaging the gradient terms of the displacement field reads

[TABLE]

As we acquire $\bar{I}_{k}=0$ from Eq. (4),

[TABLE]

After inserting Eq. (7) into Eq. (5), we obtain

[TABLE]

Now, by using the latter Eq. (8) on the right-hand side of Eq. (2), macroscale energy of an RVE is acquired as macroscale tensors are constant in space

[TABLE]

where

[TABLE]

Consequently, the macroscale energy of an RVE is expressed in terms of the gradient of macroscopic deformation. In what follows, it will be shown, by employing asymptotic homogenization analysis, that the microscale energy can be formulated in terms of the gradient of macroscopic deformation as well leading to connections between parameters.

2.2 The microscale energy for an RVE

By following the asymptotic homogenization method as in [28], we reformulate the left-hand side of Eq. (2). The asymptotic homogenization method separates length scales by using global coordinates, $\bm{X}$, for describing the global variation of the displacement, and by using local coordinates, $\bm{y}$, for describing the local fluctuation of the displacement. We refer to [69] and [35, Appendix B] for a more detailed investigation of the multiscale asymptotic analysis applied herein. We introduce the local coordinates,

[TABLE]

where $\epsilon$ is a homothetic ratio scaling global and local coordinates. We stress that the dimensions of an RVE in local coordinates can be arbitrarily chosen by varying $\epsilon$.

For example, as depicted in Fig. 2, the size of an RVE is given by $l$ in global coordinates, whereas it is denoted by $w$ in local coordinates. If we choose $l=0.001$ mm, as measured in global coordinates, $\bm{X}$, then it can be homothetically scaled to any dimension, such as $w=0.001$ mm or $w=1000$ mm in local coordinates, $\bm{y}$, by setting the homothetic ratio to $\epsilon=1.0$ or $\epsilon=10^{-6}$, in such a way that the size of the RVE is kept constant in the global coordinates. We remark that the homothetic ratio is used to describe the relationship for the sizes of an RVE between global and local coordinates; but the ratio between macroscale and microscale remains the same, $L/l=\text{const}$. We assume that the displacement field is a smooth function at the macroscopic level and $\bm{y}$-periodic in local coordinates resulting in vanishing mean local fluctuations within each RVE. Hence, the decomposition of the microscale displacement is additive into macroscale displacement and local fluctuations defined in different scales—they are independent.

By following [59], the displacement field of an RVE, $\Omega^{P}$, at global coordinates $\bm{X}$ is expanded by using asymptotic series with homothetic ratio $\epsilon$, where the corresponding coefficients in general do depend on global coordinates, $\bm{X}$, as well as local coordinates, $\bm{y}$, which are related by Eq. (11),

[TABLE]

where $\overset{n}{\bm{u}}(\bm{X},\bm{y})$ ($n$ = 0, 1, 2, …) are assumed to be $\bm{y}$-periodic. We shall see later that the first term $\overset{0}{\bm{u}}(\bm{X},\bm{y})$ is independent of $\bm{y}$. We apply now the elasticity problem in statics as it needs to be fulfilled within the RVE

[TABLE]

where the body force, $\bm{f}$, is a given function. By inserting Eq. (12) as well as using the chain rule with the aid of the relation in Eq. (11), we obtain

[TABLE]

Using the latter in Eq. (13)

[TABLE]

and then gathering terms having the same order in $\epsilon$ leads to the following terms:

• •

in the order of $\epsilon^{-2}$

[TABLE]

• •

in the order of $\epsilon^{-1}$

[TABLE]

• •

in the order of $\epsilon^{0}$

[TABLE]

By solving these partial differential equations, Eq. (12) can be rewritten as

[TABLE]

in which $\varphi_{abi}(\bm{y})$ and $\psi_{abci}(\bm{y})$ are both $\bm{y}$-periodic and they are the solutions of the following two partial differential equations:

[TABLE]

[TABLE]

A note should be made that the choice of the indices of the third order tensor $\bm{\varphi}$ and fourth order tensor $\bm{\psi}$ differs from those in [59, 15]. Since $\bm{\varphi}$ and $\bm{\psi}$ are expressed in the Cartesian coordinates, we choose to use lower indices like $\varphi_{abk}$ and $\psi_{abck}$ in the context, and they are mathematically and physically exactly identical to those in [59, 15]. We refer to Appendix for a derivation of Eq. (19), Eq. (20) and Eq. (21). Since the first term $\overset{0}{u}_{i}(\bm{X})$ depends only on the macroscopic coordinates, $\bm{X}$, it is assumed to be the known macroscopic displacement $\overset{0}{u}_{i}(\bm{X})=u^{\text{M}}_{i}(\bm{X})$ such that Eq. (12) provides

[TABLE]

We aim at defining the energy at microscale, thus, we need the gradient of the microscale displacement,

[TABLE]

with the same accuracy, i.e., after neglecting higher than second gradients and inserting Eq. (8) with the aid of Eq. (11)

[TABLE]

where

[TABLE]

By using the latter on the left-hand side of Eq. (2), microscale energy becomes

[TABLE]

with

[TABLE]

As we have assumed centro-symmetric materials, the rank 5 tensor vanishes, $\bar{\bm{G}}=0$. Immediately we observe by comparing with Eq. (9),

[TABLE]

where

[TABLE]

Therefore, we have generated an algorithm delivering effective parameters:

[TABLE]

after computing $\bm{\varphi}$ and $\bm{\psi}$ in an RVE.

3 Numerical solution of strain gradient homogenization problems

The overarching aim is to obtain classical stiffness tensors $C^{\text{M}}_{ijlm}$ and strain gradient stiffness tensors $D^{\text{M}}_{ijklmn}$, for their determination we need to solve Eqs. (20), (21). We restrict the analysis to 2D case, for the sake of simplicity such that all indices belong to $\{1,2\}$. Within the RVE, which is the computational domain, $\Omega$, Eqs. (20), (21) are solved by using the Galerkin procedure in the FEM with continuous shape functions. All boundaries are given as periodic boundary conditions, in other words, the values of $\phi_{abi}$, $\psi_{abci}$ are given by Dirichlet boundary conditions.

Indeed, the solutions of Eqs. (20), (21) are determined for specific $a$, $b$ indices (classical coefficients) as well as $a$, $b$, $c$ indices (strain gradient coefficients). Consider the case that $a=1$ and $b=1$ leading to the weak form of Eqs. (20) after multiplying by an arbitrary test function vanishing on Dirichlet boundaries and integrating by parts

[TABLE]

out of which we determine $\varphi_{abi}$ after solving for $ab={11,12,21,22}$. By knowing $\bm{\varphi}$, for example in the case of $a=1$, $b=1$, and $c=2$, we then solve

[TABLE]

The result in $\bm{\psi}$ after solving for $abc=\{111,112,121,122,211,212,221,222\}$. By inserting $\bm{\varphi}$ and $\bm{\psi}$ in Eq. (25) and then applying Eq. (30), we determine $\bm{C}^{\text{M}}$ and $\bm{D}^{\text{M}}$.

We have used the open-source software FEniCS for establishing the computations. The CAD models of the RVE have been created on the open source platform SALOME 7.6, and FEM discretizations of the CAD models were realized by a mesh generator NetGen built in SALOME 7.6. Applying the periodic conditions and building matrices are done via Python. We stress that the generated mesh has to possess perfectly matching vertices on opposite (periodic) boundaries for consistency. By NetGen, this has been automatically fulfilled by mapping the meshes between periodic surfaces. The mesh is transferred to the FEniCS and the numerical solution of weak forms have been obtained by using the iterative solver gmres with preconditioner jacobi with relative tolerance $10^{-5}$ and absolute tolerance $10^{-10}$ to ensure the accuracy of the calculations.

4 Identification of the classical and strain gradient stiffness tensors

In order to demonstrate the approach, the classical and strain gradient stiffness tensors are identified for a specific cases. Firstly, it is examined for its consistency by computing $\bm{C}^{\text{M}}$ and $\bm{D}^{\text{M}}$ in the case of a homogeneous material. As expected, the approach delivers zero for $\bm{D}^{\text{M}}$ (within the numerical tolerance). Concretely, the implementation delivers for $\bm{D}^{\text{M}}$ components $10^{-6}$ N or smaller for a material with Young’s modulus of 100 MPa and Poisson’s ratio of 0.3. This is consistent with the interpretation that for a homogeneous material all corresponding strain gradient material parameters vanish.

Then a simple geometry, the so-called square lattice structure in 2D is investigated. The square lattice structure has been widely used in engineering practice [9], as shown in Fig. 3 where gray lines build up a truss like structure. This inner structure is expected to deliver $D_{4}$ invariant material symmetry group [11, 12, 73].

As the microscale material parameters for the lattice structure, isotropic material properties are used with Young’s modulus, $E$, and Poisson’s ratio, $\nu$, as follows:

[TABLE]

Voigt’s notation is used for representing the tensors, for convenience, we refer to Table 1 and Table 2 for the chosen convention based on the work by [12].

4.1 Parameter determination for the square lattice structure

In the case of the square lattice structure, we consider the case that the material parameters of the inclusion are much smaller than those of the matrix, simply stated it is a additively manufactured truss like structure with rods out of a polymer and voids called inclusions. By choosing material properties as compiled in Table 3

and the volume fraction of the inclusion as $81\%$, we select different RVEs and determine the parameters. The RVEs are generated by repeating the corresponding basic cell, and the size of basic cell is kept constant; namely the RVEs are constituted by one cell, four cells, and nine cells, as depicted in Fig. 3.

The results for

[TABLE]

are compiled in Table 4.

In order to investigate how the size of the basic cell affects classical and strain gradient stiffness tensors, different sizes of basic cells ($0.2\times 0.2$, $0.5\times 0.5$ ) are selected and corresponding results are compared with those obtained with the basic cells size $1\times 1$, see Fig. 4 for the basic cells.

Due to the fact that these three structures share the same topology, the same material properties, as well as the same inclusion volume fraction, the corresponding classical stiffness tensors are identical. However, this is not the case for strain gradient stiffness tensors, as compiled in Table 5.

All non-vanishing parameters approach zero as the size of basic cells is decreasing. We remark that this fact is intuitively correct. Indeed, when the size of basic cells vanishes, the material becomes homogeneous resulting in vanishing $\bm{D}^{\text{M}}$. This computation also illustrates the role of the homothetic ratio $\epsilon$. To this end, let us consider the parameter $D_{221221}$ as shown in Table 5. In the case of a basic cell 1$\times$1, such a parameter is 4 times larger than that computed for the case of a basic cell 0.5$\times$0.5, and it is 25 times larger than that computed for the case of a basic cell 0.2$\times$0.2. The magnification factors (4 or 25) are equal to the square of homothetic ratios of these three basic cells as directly seen in Eq. (27).

5 Computational validation of determined parameters

In order to verify and validate the numerical values of the determined parameters, we perform three distinct computations: a computation at the microscale by incorporating the inner structure, a computation only with the determined classical stiffness tensor at the macroscale by using the homogenized structure and another computation with both determined classical stiffness tensor as well as strain gradient tensor at the macroscale by using the homogenized structure.

As suggested in [66, 82, 39], the problem of strain gradient elasticity is solved by using a weak form that, in the linear setting, leads to the $H^{2}$ norm about the trial solutions as well as test functions. Hence, the corresponding finite-dimensional approximations are guaranteed to lie in a function space which is at least of $C^{1}$ continuity. In order to obtain this property, the isogeometric FEM is employed with non-uniform rational Bezier splines (NURBS) based shape functions. The isogeometric FEM is able to ensure $C^{n}$ continuity in one single patch, which is appropriate for 2D simple geometries as being the case here. A detailed discussion of the NURBS basis and isogeometric FEM as well as the weak formulation of strain gradient elasticity can be found in [53, 54, 23, 38, 24, 25, 26, 50]. The deformation energy which quantitatively describes the overall deformation behavior of the structures are used to compare the results.

The boundary conditions for the simulations are shown in Fig. 5.

The left side of the structure is clamped and on the right side of the structure a rotation is prescribed along the center of the right edge. Two different types of computations are performed in the following subsections. In Section 5.1, the computations are done for the lattice structures with different macro sizes but the same sizes of basic cell, and in Section 5.2 we conduct computations for lattice structures with the same macro sizes but with different sizes of internal basic cells. The total volume remains the same in this case due to the fact that the ratio of the cell wall length to thickness of the basic cell is held constant with a ratio 1 to 10.

5.1 Computations for square lattices with the same basic cell sizes and varied macro sizes

In this section, the computations for the square lattice with the same sizes of the basic cell but with different macro sizes as shown in Fig. 6 are done. The size of the basic cell is 1 mm $\times$ 1 mm, and the selected lattices are of the macro sizes

• •

2 mm $\times$ 2 mm,

• •

4 mm $\times$ 4 mm,

• •

6 mm $\times$ 6 mm,

• •

10 mm $\times$ 10 mm.

The results of the simulations are shown in Fig. 7, where the vertical axis stands for the strain energy of the structures (in mJ) and the horizontal axis stands for the prescribed rotation (in rad). The black solid lines in Fig. 7 represents the results from the microscale. We understand this solution as the correct one. The blue dashed line represents the computations of the homogenized structure by using the classical stiffness tensor. The yellow dashed line represents the simulations for the homogenized structure taking the strain gradient effect into account.

The blue lines show a smaller strain energy with regard to the microscale due to the absence of the higher order strain gradient energy. We remark that, keeping the sizes of the basic cells unchanged, with increasing macro sizes of the structures, namely $L/l$ becoming larger and larger, the computations under classical elasticity theory approach to that of microscale. We may say that in a large macro scale $L/l>10$, the classical elasticity is adequate to guarantee the accuracy of the computation. However, when the macroscopic length scale is of the same order of its sizes of internal substructures, the strain gradient effect becomes significant. This phenomenon is also known as size effect.

5.2 Computations for square lattices with varied basic cell sizes and the same macro sizes

In this section, in order to further verify the identified parameters for square lattices with different basic cell sizes but the same macro sizes, computations are conducted herein. Three square lattices are selected as shown in Fig. 8. These three lattices possess the same macro sizes 4 mm $\times$ 4 mm, and the basic cell sizes of them are 1 mm $\times$ 1 mm, 0.5 mm $\times$ 0.5 mm and 0.2 mm $\times$ 0.2 mm for the left, the middle and the right lattice in Fig. 8, respectively, which divides the macro domain into 16, 64 and 400 basic cells.

The computations are shown in Fig. 9 and Fig. 10.

Fig. 9 indicates that with the increasing of basic cell sizes, the strain energy at the rotation of 0.2 radian shows an increasing trend shown under microscale. This scale-dependent (depends on $L/l$) phenomenon is also know as size effect. The computations for macroscale of elasticity are identical for these three cases due to the fact that ratio of the cell wall length to thickness of the basic cell is fixed with a ratio 1 to 10. The computation for macroscale of elasticity are independent of the scale ratio and they show a significant error compared with the microscale where the scale ratio ($L/l$) is getting smaller, and the size effect could not be ignored. When the scale ratio ($L/l$) is getting larger, which means the decreasing of the basic cell sizes, the computation for macroscale of elasticity are gradually approaching to the microscale. In such a case, for example $L/l>20$, the size effect can be ignored. It can be also observed from Fig. 10 (a), (b), (c) that the computations with strain gradient show a good match with the microscale quantitatively, which means by taking the strain gradient stiffness tensor into account, the size effect of the lattice structure is fully resolved.

6 Conclusions

A homogenization approach based on the asymptotic analysis has been exploited for developing a methodology to determine parameters in a metamaterial. Specifically, the strain gradient theory is used at the macroscale. The expressions of classical stiffness tensor and strain gradient stiffness tensor have been derived and the FEM has been successfully used to solve the partial differential equations generated from the homogenization procedure. The so-called square lattice structure has been investigated, and their material parameters are explicitly computed. The proposed approach guarantees that the parameters of strain gradient stiffness tensors vanish as the material becomes homogeneous. Moreover, it ensures that strain gradient related parameters are independent on the repetition of RVE, but dependent on the intrinsic size of the material. In order to validate the parameters determined by this methodology, additional numerical computations of the square lattice with different sizes have been performed. The numerical results show that the size effect of the lattice can be accurately captured by using the strain gradient theory with the parameters determined by the methodology applied herein. We stress that this methodology simply allows for any metamaterial made of a substructure with an RVE.

Appendix: Asymptotic solution for the displacement field

The asymptotic solution for an RVE are derived. Namely the solutions of Eq. (16), Eq. (17) and Eq. (18) are shown.

We start with Eq. (16). As $C_{ijkl}^{\text{m}}$ is a function of $\bm{y}$, the only possible general solution of Eq. (16) is to restrict $\overset{0}{u}_{i}(\bm{X})$ since it is $\bm{y}$-periodic and has a bounded gradient. The solution in the order of $\epsilon^{-2}$ can be given as

[TABLE]

Note that $\overset{0}{u}_{i}(\bm{X})$ depends only on the macroscopic coordinates; it is assumed to be the known macroscopic displacement $\overset{0}{u}_{i}(\bm{X})=u^{\text{M}}_{i}(\bm{X})$. By substituting Eq. (35) into Eq. (17), by introducing $\varphi_{abc}=\varphi_{abc}(\bm{y})$, for the inverse operation, we obtain

[TABLE]

Then the general solution of Eq. (17) can be given as

[TABLE]

where $\overset{1}{\bar{u}}_{i}=\overset{1}{\bar{u}}_{i}(\bm{X})$ are integration constants in $\bm{y}$.

Substituting Eq. (35) and Eq. (37) (with $\overset{1}{\bar{u}}_{i}(\bm{X})=0$) into Eq. (18) leads to

[TABLE]

Please note that the body force $\bm{f}$ keeps unchanged in micro- and macro-scales . We recall the governing equation in the macro-scale which reads [3]

[TABLE]

By neglecting the fourth order term in Eq. (39) and using $\overset{0}{u}_{i}(\bm{X})=u^{\text{M}}_{i}(\bm{X})$, we obtain

[TABLE]

Substituting Eq. (40) into Eq. (38) leads to

[TABLE]

As $\overset{0}{u}_{a,bc}$ is constant in $\bm{y}$, we can introduce $\psi_{abci}$ depending on $\bm{y}$ and decompose as follows:

[TABLE]

where $\psi_{abcd}=\psi_{abcd}(\bm{y})$ and $\overset{2}{\bar{u}}_{i}(\bm{X})$ are integration constants in $\bm{y}$. By substituting Eq. (42) (with $\overset{2}{\bar{u}}_{i}(\bm{X})=0$) into Eq. (41), it is found that the tensor $\psi_{abcd}$ must fulfill the following equation

[TABLE]

such that Eq. (12) provides

[TABLE]

Acknowledgements

We express our gratitude to Emilio Barchiesi, Ivan Giorgio, and Francesco dell’Isola for valuable discussions. We also thank David Kamensky for the help of implementation of isogeometric FEM in FEniCS.