International Journal of Solids and Structures 238 (2022) 111386
Available online 22 December 2021
0020-7683/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Contents lists available at ScienceDirect
International Journal of Solids and Structures
journal homepage: www.elsevier.com/locate/ijsolstr
Verification of asymptotic homogenization method developed for periodic
architected materials in strain gradient continuum
Hua Yanga, B. Emek Abalib,∗, Wolfgang H. Müllera, Salma Barbourac, Jia Li c
aTechnische Universität Berlin, Institute of Mechanics, Einsteinufer 5, 10587, Berlin, Germany
bUppsala University, Division of Applied Mechanics, Box 534, 751 21 Uppsala, Sweden
cSorbonne Paris North University, Laboratoire des Sciences des Procédés et des Matériaux, CNRS3407, 93430 Villetaneuse, France
ARTICLE INFO
Keywords:
Strain gradient elasticity
Asymptotic homogenization method
Finite element method
Constitutive parameters identification
ABSTRACT
Strain gradient theory is an accurate model for capturing the size effect and localization phenomena. However,
the challenge in identification of corresponding constitutive parameters limits the practical application of the
theory. We present and utilize asymptotic homogenization herein. All parameters in rank four, five, and six
tensors are determined with the demonstrated computational approach. Examples for epoxy carbon fiber
composite, metal matrix composite, and aluminum foam illustrate the effectiveness and versatility of the
proposed method. The influences of volume fraction of matrix, the stack of RVEs, and the varying unit cell
lengths on the identified parameters are investigated. The homogenization computational tool is applicable
to a wide class materials and makes use of open-source codes in FEniCS. We make all of the codes publicly
available in order to encourage a transparent scientific exchange.
1. Introduction
Composite materials have been widely used in engineering prac-
tice. Due to the heterogeneous nature of composites, the mechanical
properties of such materials are dependent on their substructures, for
example, the material properties of matrix and reinforcements, the
shape of inclusions, or the volume fraction of matrix, etc. An accurate
determination of effective properties of these heterogeneous media
plays an important role in the design and analysis of composites. Ex-
periments could be conceived to evaluate their effective properties, but
it is also possible to compute effective material parameters by means
of homogenization methods (Boutin,1996;Dirrenberger et al.,2019),
which reduce demands for experiments and enable to comprehend
microstructure influence on the macroscale in any complex geometries.
Homogenization techniques (Arabnejad and Pasini,2013;Chen
et al.,2020;Yvonnet et al.,2020;Jakabčin and Seppecher,2020)
allow to represent a heterogeneous elastic material, at the microscale,
as an equivalent homogeneous elastic material at the macroscale.
Although of primary importance, the conventional homogenization
fails to describe the mechanical response when the heterogeneity of
the material is of the same order of the macroscale. This inaccuracy
is due to the fact that the conventional homogenization methods are
based on a separation of scales, given by 𝜖=𝑙∕𝐿,𝑙 ≪ 𝐿. Here,
𝑙represents the typical length scale characteristic of the microstruc-
tural heterogeneity and 𝐿stands for the macroscopic length scale. If
∗Corresponding author.
the microstructure consists of relatively small heterogeneity, or the
macroscopic length is infinitely large, classical homogenization gives
an adequate estimate of the average macroscopic properties (Hollister
and Kikuchi,1992). However, if the size of heterogeneity is of the same
order of magnitude as that of the macroscopic problem, conventional
homogenization technique fails. For example, the size effect occurs
when the length scale of the macroscopic heterogeneous materials (𝐿)
approaches the length scale of the underlying heterogeneity (𝑙). An up-
scaling of Cauchy theory indicates that additional terms are necessary
in the constitutive equations in order to predict the size effect observed
in experiments (Müller,2020).
In order to incorporate these additional terms, different homoge-
nization techniques are proposed in the literature, for example in the
framework of generalized mechanics (Mindlin and Eshel,1968;Erin-
gen,1999;Altenbach and Forest,2016) such as micropolar theory (Ku-
mar and McDowell,2004;Dos Reis and Ganghoffer,2012), couple
stress (Skrzat and Eremeyev,2020), strain gradient theory (Yvon-
net et al.,2020;Kouznetsova et al.,2002;Goda and Ganghoffer,
2016;Abdoul-Anziz et al.,2019;Weeger,2021;Forest and Trinh,
2011), and micromorphic continuum (Rokoš et al.,2019). The task
of obtaining homogenized constitutive equations for generalized con-
tinua is challenging and a number of debates are active in the liter-
ature (Kumar and McDowell,2004;Dos Reis and Ganghoffer,2012;
Liu and Su,2009;Eremeyev,2016;Ganghoffer and Reda,2021).
https://doi.org/10.1016/j.ijsolstr.2021.111386
Received 2 June 2021; Received in revised form 1 December 2021; Accepted 2 December 2021
International Journal of Solids and Structures 238 (2022) 111386
2
H. Yang et al.
Many methods have been proposed to construct strain gradient con-
tinua by means of asymptotic homogenization approaches (Bacigalupo
et al.,2018;Boutin et al.,2017), multi-scale computational approaches
(Kouznetsova et al.,2004), dynamic methods (Rosi et al.,2018;Rosi,
2020), and several other identification techniques (dell’Isola et al.,
2019a,b;Misra and Poorsolhjouy,2015;Alibert and Della Corte,2019;
Rahali et al.,2015). Asymptotic homogenization method improves
descriptions by exploiting higher order terms and considering their
role in macroscopic behaviors. In Li (2011), Li and Zhang (2013),
Barboura and Li (2018), Yang et al. (2020) and Abali and Barchiesi
(2021), an asymptotic homogenization based solution has been utilized
to determine parameters of composite materials. Two issues were
addressed therein. One is that the identified strain gradient parameters
are all zero when structures are homogeneous. The other one is that
these parameters are independent of stack of RVEs.
In this paper, we briefly recall the homogenization method de-
scribed in Yang et al. (2020), which is based on the formal analysis
in Barboura and Li (2018). A complete computational methodology
determining all parameters in 2D and 3D has been achieved recently
in Abali and Barchiesi (2021). We basically use the same procedure and
verify the computational implementation by several numerical sanity
checks. In this way, we reveal an important limitation that remains
undetected during a formal analysis. Homogenization begins with two
materials of different properties. In the case of one material with a
nearly zero stiffness, the difference in properties may cause numerical
problems in the implementation. For example, a structure with voids
is a benchmark case for this issue. Indeed, we propose a change in the
formulation in order to make the numerical implementation robust and
the methodology more general. We apply the method to determine 2D
and 3D composite materials effective parameters as well as verify the
results by additional simulations. For determining parameters, we use
epoxy carbon fiber composite material, SiC/Al metal matrix composite,
and aluminum foam. As a benchmark, we choose aluminum foam.
The content of this paper is structured as follows: In Section 2, the
underlying method is explained in order to clearly present the addition
proposed herein. In Section 3, the details of a complete numerical
implementation are demonstrated. In Section 4, effective parameters
for 2D and 3D composite materials including epoxy-carbon fiber com-
posites, metal matrix composite, and aluminum foam are identified as
well as the aforementioned two challenges have been exploited for
checking the robustness of the implementation. In Section 5we discuss
the positive definiteness and in Section 6we verify the parameters by
using a strain gradient simulation. The homogenization computational
tool is developed based on open-source codes in FEniCS. It allows for
all kinds of 2D or 3D composite materials constructed by periodic
microstructures. The codes are made publicly available in Abali (2020)
in order to enable a transparent scientific exchange.
2. Homogenization method
We start from an assertion that the deformation energy for the
domain representing RVE, 𝛺𝑃, at the microscale is equal to the energy
for the RVE at the macroscale, namely
∫𝛺𝑃
𝑤md𝑉=∫𝛺𝑃
𝑤Md𝑉 . (1)
The superscripts ‘‘m’’ and ‘‘M’’ are used to denote microscopic and
macroscopic quantities, respectively. At the microscale, detailed mi-
crostructures are present in the RVE. At the macroscale, the same
domain is modeled by a homogeneous ‘‘metamaterial’’. We emphasize
that an RVE may be different from a unit cell. A unit cell is the
simplest repeating unit of heterogeneity. Spatial repetition of unit cells
composes an RVE. At the microscale, the first order theory is used,
as a consequence, we need to have a second order theory at the
macroscale (Mandadapu et al.,2021). Now by starting with a linear
strain measure, in the case of a linear material model, we obtain a
quadratic deformation energy,
∫𝛺𝑃
1
2𝐶m
𝑖𝑗𝑘𝑙𝑢m
𝑖,𝑗 𝑢m
𝑘,𝑙 d𝑉=∫𝛺𝑃(1
2𝐶M
𝑖𝑗𝑘𝑙𝑢M
𝑖,𝑗 𝑢M
𝑘,𝑙 +𝐺M
𝑖𝑗𝑘𝑙𝑚𝑢M
𝑖,𝑗 𝑢M
𝑘,𝑙𝑚
+1
2𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛𝑢M
𝑖,𝑗𝑘𝑢M
𝑙,𝑚𝑛 )d𝑉 .
(2)
Displacement fields at micro- and macroscales are indicated by ‘‘m’’ and
‘‘M’’, respectively. 𝐶m
𝑖𝑗𝑘𝑙 is given in each material point of the RVE. We
begin with the known microscale and search for its corresponding ho-
mogenized effective parameters. The effective coefficients, 𝐶M
𝑖𝑗𝑘𝑙,𝐺M
𝑖𝑗𝑘𝑙𝑚,
and 𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛 are the unknowns that we are searching for. We emphasize
that the quadratic energy and symmetric strain measure lead to the
minor symmetries 𝐶m
𝑖𝑗𝑘𝑙 =𝐶m
𝑗𝑖𝑘𝑙 =𝐶m
𝑖𝑗𝑙𝑘, 𝐶M
𝑖𝑗𝑘𝑙 =𝐶M
𝑗𝑖𝑘𝑙 =𝐶M
𝑖𝑗𝑙𝑘,𝐺M
𝑖𝑗𝑘𝑙𝑚 =
𝐺M
𝑗𝑖𝑘𝑙𝑚 =𝐺M
𝑖𝑗𝑙𝑘𝑚,𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛 =𝐷M
𝑗𝑖𝑘𝑙𝑚𝑛 =𝐷M
𝑖𝑗𝑘𝑚𝑙𝑛 and major symmetries
𝐶m
𝑖𝑗𝑘𝑙 =𝐶m
𝑘𝑙𝑖𝑗 ,𝐶M
𝑖𝑗𝑘𝑙 =𝐶M
𝑘𝑙𝑖𝑗 ,𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛 =𝐷M
𝑙𝑚𝑛𝑖𝑗𝑘 of the classical and strain
gradient stiffness tensors. In what follows, the connections between the
microscopic material parameters and macroscopic ones are established.
Let us investigate the macroscopic case for an RVE, 𝛺𝑃. Firstly, the
geometric center of the RVE is defined as c
𝑿=1
𝑉∫𝛺𝑃𝑿d𝑉. A Taylor
expansion of the macroscopic displacement around the center of the
RVE is written as
𝑢M
𝑖(𝑿) = 𝑢M
𝑖|||c
𝑿
+𝑢M
𝑖,𝑗 |||c
𝑿
(𝑋𝑗−
c
𝑋𝑗) + 1
2𝑢M
𝑖,𝑗𝑘|||c
𝑿
(𝑋𝑗−
c
𝑋𝑗)(𝑋𝑘−
c
𝑋𝑘),
𝑢M
𝑖,𝑙(𝑿) = 𝑢M
𝑖,𝑗 |||c
𝑿
𝛿𝑗𝑙 +1
2𝑢M
𝑖,𝑗𝑘|||c
𝑿
(𝛿𝑗𝑙(𝑋𝑘−
c
𝑋𝑘)+(𝑋𝑗−
c
𝑋𝑗)𝛿𝑘𝑙),
=𝑢M
𝑖,𝑙|||c
𝑿
+𝑢M
𝑖,𝑙𝑘|||c
𝑿
(𝑋𝑘−
c
𝑋𝑘),
𝑢M
𝑖,𝑙𝑚(𝑿) = 𝑢M
𝑖,𝑙𝑘|||c
𝑿
𝛿𝑘𝑚 =𝑢M
𝑖,𝑙𝑚|||c
𝑿
.
(3)
Then by using the spatial averaging, we have
⟨𝑢M
𝑖,𝑗 ⟩=1
𝑉∫𝛺𝑃
𝑢M
𝑖,𝑗 d𝑉=𝑢M
𝑖,𝑗 |||c
𝑿
,
⟨𝑢M
𝑖,𝑗𝑘⟩=1
𝑉∫𝛺𝑃
𝑢M
𝑖,𝑗𝑘 d𝑉=𝑢M
𝑖,𝑗𝑘|||c
𝑿
.
(4)
Therefore, the macroscopic energy of an RVE reads as follows (the de-
tailed derivation can be found in Yang et al.,2020), as the macroscopic
stiffness tensors are constant in space,
∫𝛺𝑃(1
2𝐶M
𝑖𝑗𝑘𝑙𝑢M
𝑖,𝑗 𝑢M
𝑘,𝑙 +𝐺M
𝑖𝑗𝑘𝑙𝑚𝑢M
𝑖,𝑗 𝑢M
𝑘,𝑙𝑚 +1
2𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛𝑢M
𝑖,𝑗𝑘𝑢M
𝑙,𝑚𝑛)d𝑉 ,
=𝑉
2𝐶M
𝑖𝑗𝑙𝑚⟨𝑢M
𝑖,𝑗 ⟩⟨𝑢M
𝑙,𝑚⟩+𝑉 𝐺M
𝑖𝑗𝑘𝑙𝑚⟨𝑢M
𝑖,𝑗 ⟩⟨𝑢M
𝑘,𝑙𝑚⟩
+𝑉
2(𝐶M
𝑖𝑗𝑙𝑚
𝐼𝑘𝑛 +𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛)⟨𝑢M
𝑖,𝑗𝑘⟩⟨𝑢M
𝑙,𝑚𝑛⟩,
(5)
𝐼𝑘𝑛 =1
𝑉∫𝛺𝑃
(𝑋𝑘−
c
𝑋𝑘)(𝑋𝑛−
c
𝑋𝑛) d𝑉 . (6)
At the microscale, the asymptotic homogenization method is used to
approximate the deformation energy for the RVE. We introduce a small
parameter 𝜖, which is defined as 𝜖=𝑙
𝐿, where 𝑙is the characteristic
length of the microstructure, 𝐿is the length of the macroscopic struc-
ture as shown in Fig. 1. We remark that 𝜖is the so-called homothetic
ratio, which shows the scaling law for strain gradient moduli. This
property will be illustrated later. A local coordinate is then introduced
as
𝑦𝑗=1
𝜖(𝑋𝑗−
c
𝑋𝑗),(7)
which is used to describe the local fluctuations caused by microscopic
heterogeneity. Variable 𝑿is associated with the macroscopic scale. The
displacement field for the RVE at the microscale is thus approximated
with regard to 𝜖as
𝒖m(𝑿) = 0
𝒖(𝑿,𝒚) + 𝜖1
𝒖(𝑿,𝒚) + 𝜖22
𝒖(𝑿,𝒚) + ⋯.(8)
International Journal of Solids and Structures 238 (2022) 111386
3
H. Yang et al.
Fig. 1. The heterogeneous continuum and its equivalent homogenized continuum.
For a linear elastostatics problem, we propose to write the governing
equations within the RVE, as follows:
(𝐶m
𝑖𝑗𝑘𝑙𝑢m
𝑘,𝑙),𝑗 +𝜌m𝑓𝑖= 0 ,(9)
where 𝜌m𝑓𝑖are volume forces, 𝜌mis the mass density at the microscale,
hence it is a function in 𝑿effected by the heterogeneous structure. We
stress that this interpretation of using mass density at the microscale
has an important remedy to the generality of the computational im-
plementation. This change in the formulation is for the first time in
the literature and numerically beneficial in the case of voids. Since
voids possess numerically zero mass density, their contribution to the
homogenization is weakened by exploiting this amendment in the
formulation. By substituting Eq. (8) to Eq. (9) and gathering terms
having the same order in 𝜖leads to the following equations:
•in the order of 𝜖−2
𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙
𝜕0
𝑢𝑘
𝜕𝑦𝑙)= 0 ; (10)
•in the order of 𝜖−1
(𝐶m
𝑖𝑗𝑘𝑙
𝜕0
𝑢𝑘
𝜕𝑦𝑙),𝑗 +𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙
0
𝑢𝑘,𝑙)+𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙
𝜕1
𝑢𝑘
𝜕𝑦𝑙)= 0 ; (11)
•in the order of 𝜖0
(𝐶m
𝑖𝑗𝑘𝑙
0
𝑢𝑘,𝑙),𝑗 +(𝐶m
𝑖𝑗𝑘𝑙
𝜕1
𝑢𝑘
𝜕𝑦𝑙),𝑗 +𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙
1
𝑢𝑘,𝑙)
+𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙
𝜕2
𝑢𝑘
𝜕𝑦𝑙)+𝜌m𝑓𝑖= 0 .
(12)
The only possible solution of Eq. (10) is to restrict 0
𝑢𝑖(𝑿)as
0
𝑢𝑖=0
𝑢𝑖(𝑿).(13)
Because 0
𝑢𝑖(𝑿)is only dependent on the macroscopic coordinates, from
Eq. (8), by a coefficient comparison, we obtain that it may be cho-
sen as the macroscopic displacement 0
𝑢𝑖(𝑿) = 𝑢M
𝑖(𝑿). After substitut-
ing Eq. (13) into Eq. (11), and introducing 𝜑𝑎𝑏𝑐 =𝜑𝑎𝑏𝑐 (𝒚)which is
𝒚-periodic with zero average value ∫𝛺𝑃𝜑𝑎𝑏𝑖 d𝑉= 0, we obtain
𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙(𝜕𝜑𝑎𝑏𝑘
𝜕𝑦𝑙
+𝛿𝑎𝑘𝛿𝑏𝑙))= 0 .(14)
Consequently, the solution of Eq. (11) is given as
1
𝑢𝑖=𝜑𝑎𝑏𝑖𝑢M
𝑎,𝑏(𝑿) +
1
𝑢𝑖(𝑿),(15)
where 1
𝑢𝑖=
1
𝑢𝑖(𝑿)are integration constants.
By recalling the governing equation at the macroscale with an
analogous suggestion to use a macroscale mass density, we have
𝐶M
𝑖𝑗𝑘𝑙𝑢M
𝑘,𝑙𝑗 −𝐷M
𝑖𝑗𝑘𝑙𝑚𝑛𝑢M
𝑙,𝑚𝑛𝑘𝑗 +𝜌M𝑓𝑖= 0 ,(16)
with 𝜌M=1
𝑉∫𝛺𝑃𝜌md𝑉and the usual axiom that body forces are
scale independent such that 𝑓𝑖remains the same at the micro- and
macroscales. By neglecting the fourth order term in Eq. (16), we obtain
𝑓𝑖= −
𝐶M
𝑖𝑗𝑘𝑙𝑢M
𝑘,𝑙𝑗
𝜌M.(17)
By plugging Eq. (13), Eq. (15) (with 1
𝑢𝑖(𝑿)=0), and Eq. (17) into
Eq. (12) and introducing 𝜓𝑎𝑏𝑐𝑖 which is 𝒚-periodic with zero average
∫𝛺𝑃𝜓𝑎𝑏𝑐𝑖 d𝑉= 0, the solution of 2
𝑢𝑖may be given as:
2
𝑢𝑖=𝜓𝑎𝑏𝑐𝑖𝑢M
𝑎,𝑏𝑐 (𝑿) +
2
𝑢𝑖(𝑿),(18)
where 2
𝑢𝑖(𝑿)are integration constants in 𝒚. The fourth order tensor 𝜓𝑎𝑏𝑐𝑑
must satisfy
𝜕
𝜕𝑦𝑗(𝐶m
𝑖𝑗𝑘𝑙(𝜕𝜓𝑎𝑏𝑐𝑘
𝜕𝑦𝑙
+𝜑𝑎𝑏𝑘𝛿𝑙𝑐 ))+𝐶m
𝑖𝑐𝑘𝑙(𝜕𝜑𝑎𝑏𝑘
𝜕𝑦𝑙
+𝛿𝑘𝑎𝛿𝑙𝑏)−𝜌m
𝜌M𝐶M
𝑖𝑐𝑎𝑏 = 0 .
(19)
We emphasize that the last term is a source term simply applying the
loading to the system at the microscale by considering mass densities.
By neglecting the mass density ratio, one applies a source term even in
the case of voids that may lead to numerically inconsistent results for
𝝍parameters. The microscale displacement field is rewritten as
𝑢m
𝑖(𝑿,𝒚) = 𝑢M
𝑖(𝑿) + 𝜖𝜑𝑎𝑏𝑖(𝒚)𝑢M
𝑎,𝑏(𝑿) + 𝜖2𝜓𝑎𝑏𝑐𝑖(𝒚)𝑢M
𝑎,𝑏𝑐 (𝑿) + ⋯.(20)
By using Eq. (20) and the latter on the left-hand side of Eq. (2) the
microscopic energy becomes
∫𝛺𝑃
1
2𝐶m
𝑖𝑗𝑘𝑙𝑢m
𝑖,𝑗 𝑢m
𝑘,𝑙 d𝑉=
𝑉
2(
𝐶𝑎𝑏𝑐𝑑 ⟨𝑢M
𝑎,𝑏⟩⟨𝑢M
𝑐,𝑑 ⟩+
𝐺𝑎𝑏𝑐𝑑𝑒⟨𝑢M
𝑎,𝑏⟩⟨𝑢M
𝑐,𝑑𝑒⟩+
𝐷𝑎𝑏𝑐𝑑𝑒𝑓 ⟨𝑢M
𝑎,𝑏𝑐 ⟩⟨𝑢M
𝑑,𝑒𝑓 ⟩),
(21)
International Journal of Solids and Structures 238 (2022) 111386
4
H. Yang et al.
with
𝐶𝑎𝑏𝑐𝑑 =1
𝑉∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝐿𝑎𝑏𝑖𝑗 𝐿𝑐𝑑𝑘𝑙 d𝑉 ,
𝐺𝑎𝑏𝑐𝑑𝑒 =2𝜖
𝑉∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝐿𝑎𝑏𝑖𝑗 𝑀𝑐𝑑𝑒𝑘𝑙 d𝑉
𝐷𝑎𝑏𝑐𝑑𝑒𝑓 =𝜖2
𝑉∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝑀𝑎𝑏𝑐𝑖𝑗 𝑀𝑑𝑒𝑓𝑘𝑙 d𝑉 ,
(22)
The appearance of 𝜖2is due to the fact that Eq. (26) is expressed in the
local coordinate 𝒚(The fifth order tensor 𝑴is only related to 𝒚).
𝐿𝑎𝑏𝑖𝑗 =𝛿𝑖𝑎𝛿𝑗𝑏 +𝜕𝜑𝑎𝑏𝑖
𝜕𝑦𝑗
,
𝑀𝑎𝑏𝑐𝑖𝑗 =𝑦𝑐(𝛿𝑖𝑎𝛿𝑗𝑏 +𝜕𝜑𝑎𝑏𝑖
𝜕𝑦𝑗)+(𝜑𝑎𝑏𝑖𝛿𝑗𝑐 +𝜕𝜓𝑎𝑏𝑐𝑖
𝜕𝑦𝑗).
(23)
Based on Eq. (2) the effective parameters are calculated by
𝐶M
𝑎𝑏𝑐𝑑 =1
𝑉∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝐿𝑎𝑏𝑖𝑗 𝐿𝑐𝑑𝑘𝑙 d𝑉 , (24)
𝐺M
𝑎𝑏𝑐𝑑𝑒 =𝜖
𝑉∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝐿𝑎𝑏𝑖𝑗 𝑀𝑐𝑑𝑒𝑘𝑙 d𝑉 , (25)
𝐷M
𝑎𝑏𝑐𝑑𝑒𝑓 =𝜖2
𝑉(∫𝛺𝑃
𝐶m
𝑖𝑗𝑘𝑙𝑀𝑎𝑏𝑐𝑖𝑗 𝑀𝑑𝑒𝑓𝑘𝑙 d𝑉−𝐶M
𝑎𝑏𝑑𝑒 ∫𝛺𝑃
𝑦𝑐𝑦𝑓d𝑉).(26)
It should be remarked that the Eq. (24) coincides with the well
known asymptotic homogenization method. The classical stiffness ten-
sor is scale independent. However, as observed from the Eq. (26), strain
gradient stiffness parameters depend on 𝜖2. Indeed, these parameters
emerge related to the substructure and vanish as 𝜖= 0 meaning that the
substructure diminishes. We stress that this distinction is of importance
and comes out of the proposed methodology quite naturally. As obvious
in Eq. (7), the homothetic ratio, 𝜖, acts as a multiplier between the
macroscopic length scale (in global coordinates, 𝑿) and microscopic
length scale (in local coordinates, 𝒚). In this way, we acquire different
𝑮Mand 𝑫Mcoefficients for the same RVE in larger structures without
repeating the calculations. The role of 𝜖will be further illustrated by
using numerical examples.
3. Numerical implementation
In order to identify effective parameters, Eqs. (24) and (26) need to
be resolved, which requires 𝝋and 𝝍. The tensors 𝝋and 𝝍are the solu-
tions of Eqs. (14) and (19), which are solved numerically by the finite
element method. As shown in Fig. 2, six cases 𝜑11𝑖, 𝜑22𝑖, 𝜑33𝑖, 𝜑23𝑖, 𝜑13𝑖,
𝜑12𝑖in total in 3D need to be computed under periodic boundary
conditions. After using integration by parts, considering the constraints
of zero average for 𝝋, the following weak form for 𝜑𝑎𝑏𝑘 is generated
∫𝛺𝑃(𝐶m
𝑖𝑗𝑘𝑙(𝜕𝜑𝑎𝑏𝑘
𝜕𝑦𝑙
+𝛿𝑎𝑘𝛿𝑏𝑙))𝜕δ𝜑𝑎𝑏𝑖
𝜕𝑦𝑗
d𝑉+ δ ∫𝜆𝑎𝑏𝑖𝜑𝑎𝑏𝑖 d𝑉= 0 (27)
and then immediately we have
∫𝛺𝑃(𝐶m
𝑖𝑗𝑘𝑙(𝜕𝜑𝑎𝑏𝑘
𝜕𝑦𝑙
+𝛿𝑎𝑘𝛿𝑏𝑙))𝜕δ𝜑𝑎𝑏𝑖
𝜕𝑦𝑗
d𝑉
+∫𝛺𝑃
𝜆𝑎𝑏𝑖δ𝜑𝑎𝑏𝑖 d𝑉+∫𝛺𝑃
δ𝜆𝑎𝑏𝑖𝜑𝑎𝑏𝑖 d𝑉= 0 ,
(28)
where over underlined indices, no summation convention is applied. All
fields with a variational delta, δ, denote a corresponding test function
such that 𝝋and 𝝀are unknowns. For each case of 𝜑11𝑖, 𝜑22𝑖, 𝜑33𝑖, 𝜑23𝑖,
𝜑13𝑖, 𝜑12𝑖, a corresponding Lagrange multiplier 𝜆11𝑖, 𝜆22𝑖, 𝜆33𝑖, 𝜆23𝑖, 𝜆13𝑖,
𝜆12𝑖, is employed in order to enforce the zero average constrains of
Table 1
Material properties used for 2D epoxy-carbon fiber composite. 𝐸Young’s modulus, 𝜈
Poisson’s ratio, and 𝜌mass density.
Type 𝐸in GPa 𝜈 𝜌 in kg/m3
Matrix (Epoxy) 17.3 0.35 1780
Inclusion (Carbon fiber) 35.9 0.30 1650
Table 2
Voigt notation used for 2D strain tensors.
𝐴123
𝑖𝑗 11 22 12
Table 3
Voigt kind-notation used for 2D strain-gradient tensors.
𝜃123456
𝑖𝑗𝑘 111 112 221 222 121 122
𝝋(Bleyer,2018). Likewise, the weak form for calculating 𝜓𝑎𝑏𝑐𝑖 reads
∫𝛺𝑃((𝐶m
𝑖𝑗𝑘𝑙(𝜕𝜓𝑎𝑏𝑐𝑘
𝜕𝑦𝑙
+𝜑𝑎𝑏𝑘𝛿𝑙𝑐 ))𝜕δ𝜓𝑎𝑏𝑐𝑖
𝜕𝑦𝑗
−𝐶m
𝑖𝑐𝑘𝑙(𝜕𝜑𝑎𝑏𝑘
𝜕𝑦𝑙
+𝛿𝑘𝑎𝛿𝑙𝑏)δ𝜓𝑎𝑏𝑐𝑖 +𝜌m
𝜌M𝐶M
𝑖𝑐𝑎𝑏δ𝜓𝑎𝑏𝑐𝑖 )d𝑉
+∫𝛺𝑃
𝜆𝑎𝑏𝑐𝑖δ𝜓𝑎𝑏𝑐𝑖 d𝑉+∫𝛺𝑃
δ𝜆𝑎𝑏𝑐𝑖𝜓𝑎𝑏𝑐𝑖 d𝑉= 0 .
(29)
There are 18 weak forms in 3D to be solved for 𝜓111𝑖,𝜓112𝑖, ...𝜓123𝑖..
The weak forms have been solved by the FEniCS platform. CAD
model and mesh files are created by using an open-source software
SALOME (Abali,2017). Triangle for surfaces elements and tetrahedron
for volume elements are used to discretize the system by using the
algorithms from NetGen Mesh Generator. An RVE needs to fulfill pe-
riodic boundary conditions such that the corresponding edges (in 2D)
or surfaces (in 3D) are matching for nodes to be defined as the same
degree of freedom in order to enforce the periodic boundary conditions
as shown in Fig. 3.
4. Numerical examples
The proposed homogenization method provides a unified analysis
for general 2D and 3D composites. It may be used to homogenize fiber
reinforced composites, particulate composites, and porous materials. In
order to show the predictive capability of the proposed method, four
examples are demonstrated in the following.
4.1. 2D epoxy-carbon fiber composite
A 2 dimensional carbon fibers reinforced epoxy composite structure
is investigated. The material properties1for both constituents (matrix
and inclusion) are shown in Table 1. The size of the unit cell is 1 mm.
The fiber is of circular shape, its radius is 0.45 mm, thus, the volume
fraction of matrix is 36.4%.
Voigt notations as presented in Tables 2,3are used to represent
rank four, five, six tensors as matrices (analogous to Voigt’s notation).
A bottleneck in the homogenization may be the missing convergence
criteria. We propose a simple yet effective approach by using the mate-
rial symmetry class of the analyzed microstructure. Owing to the cubic
material symmetry, we know that 𝐶1111 =𝐶2222,𝐷111111 =𝐷222222. The
convergence analysis is conducted as shown in Table 4, by comparing
the ratios 𝐶1111∕𝐶2222, and 𝐷111111∕𝐷222222. When they tend to be 1, the
computation is converged. The solutions for 𝝋and 𝝍are presented in
1Values of material properties are taken from matweb.com.
International Journal of Solids and Structures 238 (2022) 111386
5
H. Yang et al.
Fig. 2. The flowchart of the numerical implementation.
Fig. 3. Periodic boundary conditions applied in FEM. Left: Right/bottom (red) and
left/top (green) edges have the same mesh. Right: Only corresponding surfaces are
shown, they have the same mesh. Same mesh is necessary for implementing periodic
boundary conditions. (For interpretation of the references to color in this figure legend,
the reader is referred to the web version of this article.)
Table 4
Convergence analysis. With the increasing of degrees of freedom, the ratios 𝐶1111∕𝐶2222 ,
and 𝐷111111∕𝐷222222 reach 1.
DOFs 𝐶1111 GPa 𝐶2222 GPa 𝐶1111∕𝐶2222 𝐷111111 N𝐷222222 N𝐷111111∕𝐷222222
1342 38.6 38.7 99.7% 510.4 496.0 103.0%
22362 38.9 38.9 100.0% 505.3 506.1 100.0%
90226 39.0 39.0 100.0% 506.4 505.8 100.0%
Fig. 4. It is observed as expected that these fluctuations are all periodic.
Furthermore, due to the fact that the material is cubic, rotating 𝜑22,
𝜓111,𝜓221,𝜓122 by 90◦gives the same shapes as 𝜑11,𝜓222,𝜓112,
𝜓121.
Fig. 4. Solutions for 𝝋and 𝝍. Color distribution showing the cubic symmetry
resulted local fluctuation in 𝝋and 𝝍fields. Color bars are omitted since we analyze
qualitatively.
The identified effective classical and stain gradient stiffness tensors
are shown as follows:
𝐶M
𝐴𝐵 =⎛⎜⎜⎝
39.0 18.0 0.0
18.0 39.0 0.0
0.0 0.0 10.0⎞⎟⎟⎠
GPa,
International Journal of Solids and Structures 238 (2022) 111386
6
H. Yang et al.
Table 5
Voigt notation used for 2D strain-gradient tensors proposed in Auffray et al. (2013).
𝛼123456
𝑖𝑗𝑘 111 221 122 222 112 121
𝐷M
𝜃𝛾 =⎛⎜⎜⎜⎜⎜⎜⎝
506.4 181.9 0.0 −0.0 0.0 −182.2
181.9 −299.4 0.0 0.0 0.0 −176.2
0.0 0.0 181.2 −175.4 −183.0 0.0
0.0 0.0 −175.4 −298.5 181.2 0.0
0.0 0.0 −183.0 181.2 505.8 0.0
−182.2 −176.2 0.0 0.0 0.0 181.0
⎞⎟⎟⎟⎟⎟⎟⎠
N.
It is found that there are three independent parameters in the stiffness
tensor and six independent parameters in the strain gradient stiffness
tensor. This observation is consistent with Auffray et al. (2015,2009)
for cubic materials. Albeit we circumvent of showing, the implemen-
tation successfully computes all parameters of 𝑮as (numerical) zeros,
as expected from the cubic material symmetry as well. By using the
Voigt notation similar to the approach as in Auffray et al. (2015,
2009,2013) in Table 5, the strain gradient stiffness matrix is made to
be block-diagonal; each diagonal block matrix includes only non-zero
parameters, and each diagonal block matrix is invariant under every
cyclic permutation of 𝑿axis, 𝒀axis, and 𝒁axis (Auffray et al.,2013).
Therefore, the Voigt notation proposed in Auffray et al. (2013) will be
used throughout the paper.
𝐷M
𝛼𝛽 =⎛⎜⎜⎜⎜⎜⎜⎝
506.4 181.9 −182.2 0.0 0.0 0.0
181.9 −299.4 −176.2 0.0 0.0 0.0
−182.2 −176.2 181.0 0.0 0.0 0.0
0.0 0.0 0.0 505.8 181.2 −183.0
0.0 0.0 0.0 181.2 −298.5 −175.4
0.0 0.0 0.0 −183.0 −175.4 181.2
⎞⎟⎟⎟⎟⎟⎟⎠
N.
4.2. Interpretation of the homothetic ratio
When determining the strain gradient moduli, physical relevance
of the so-called homothetic ratio, 𝜖, is needed for determining the
correct value. Two coordinate systems are scaled to each other by this
homothetic ratio. Let us consider specific cases as shown in Fig. 5. In
Fig. 5(a), the macroscopic length is 𝐿= 4 mm and the microscopic
length is 𝑙= 1 mm with 𝜖=𝑙
𝐿=1
4. RVE is of length 𝐿=1 mm in global
coordinates 𝑿, but it is measured as 𝑙=4 mm in local coordinate 𝒚.
Since Eq. (26) are expressed in the local coordinate 𝒚, the parameters
in 𝐷M
𝑎𝑏𝑐𝑑𝑒𝑓 are calculated in the local coordinate. Thus, the length of the
computational domain in Eq. (26) is 4 mm. Likewise, in Fig. 5(b), the
length of integration domain is 2 times larger than that in Fig. 5(a).
However 𝜖=1
8is half of the former one. This leads to the equal values
for strain gradient moduli. Consequently, in the last section, 𝜖can be
chosen as, for example, 1
4or 1
8, as long as the corresponding length of
integration domain is chosen accordingly.
Indeed, a scaling rule occurs for the strain gradient moduli. For
example, in Fig. 5(c), the length of RVE is half of that in Fig. 5(a).
The same macroscopic length equals calculated integrals in Eq. (26).
The differences of the obtained strain gradient parameters originate
from the 𝜖2as presented in Eq. (26). The strain gradient parameters
for Fig. 5(a) are 4 times larger than those for Fig. 5(c). This scaling
factor is calculated as the ratio between 𝜖2, also equal to the square of
ratio of the unit cell lengths. Therefore, herein, we conclude that the
strain gradient moduli are indeed not related to the macroscopic length
but the microscopic length. This interpretation is indeed in coincidence
with the well-known size effect in the literature. We emphasize that the
substructure affects the values in 𝐶M
𝑖𝑗𝑘𝑙, but not its ratio with respect to
the macroscale. Therefore, for different substructures, 𝐶M
𝑖𝑗𝑘𝑙 needs to be
recalculated. For the same substructure but different homothetic ratios,
they remain the same.
Table 6
Voigt notation used for 3D strain tensors.
𝐴1 2 3 4 5 6
𝑖𝑗 11 22 33 23 13 12
4.3. 3D cases
In the followings, we consider 3D cases. The effective parameters
in the classical stiffness tensor and strain gradient stiffness tensor
for a carbon fibers reinforced epoxy composite, a (hard) spherical
particles reinforced (soft) matrix, a metal matrix composite, and an
aluminum foam will be investigated. The used Voigt notations for these
3 dimensional cases are displayed in Tables 6 and 7.
4.3.1. 3D fiber reinforced composite
Carbon fiber is modeled by using a cylindrical inclusion in 3D. In
order to compare and validate the results, the same material properties
shown in Table 1 are used for inclusion and matrix. The radius of
the cylinder is of 0.45 mm so that the volume fraction of matrix
reads 36.4%, which are both equal to the example shown in 2D. The
calculated parameters are shown as follows and as in Box I.
𝐶M
𝐴𝐵 =⎛⎜⎜⎜⎜⎜⎜⎝
38.6 17.9 18.0 0.0 0.0 0.0
17.9 38.6 18.0 0.0 0.0 0.0
18.0 18.0 40.1 0.0 0.0 0.0
0.0 0.0 0.0 10.2 0.0 0.0
0.0 0.0 0.0 0.0 10.2 0.0
0.0 0.0 0.0 0.0 0.0 9.7
⎞⎟⎟⎟⎟⎟⎟⎠
GPa,
We stress that the algorithm computes 𝑮as well, but as expected
from the centro-symmetry in the substructure, all coefficients of 𝑮van-
ish. The unidirectional laminate carbon reinforced epoxy composite is
a transverse isotropic material. There are five independent parameters
in the classical stiffness tensor, as shown below
𝐶M
𝐴𝐵 =⎛⎜⎜⎜⎜⎜⎜⎝
𝑐1𝑐1− 2𝑐5𝑐20.0 0.0 0.0
𝑐1− 2𝑐5𝑐1𝑐20.0 0.0 0.0
𝑐2𝑐2𝑐30.0 0.0 0.0
0.0 0.0 0.0𝑐40.0 0.0
0.0 0.0 0.0 0.0𝑐40.0
0.0 0.0 0.0 0.0 0.0𝑐5
⎞⎟⎟⎟⎟⎟⎟⎠
.
We stress that the computed parameters are satisfying this condition
within a tolerance of ±6.7%. After investigating the strain gradient
stiffness tensor, we find the relations between higher order parameters
as shown in Fig. 6.
Excluding the parameters of the same value, there are 28 parameters
in 𝑫for this transverse isotropic material which is given in Box II. We
emphasize that some of the 28 parameters could be linearly dependent
that leads to a reduction of independent coefficients. Moreover, the
corresponding parameters in 2D and 3D stiffness tensors are equal
within a tolerance of ±4.4%, for example, 𝐶1111 or 𝐷111111 in the 2D
stiffness tensors are equal to those in the 3D tensors. This verifies
the calculated results. In order to further examine the homogeniza-
tion method, computations for different volume fraction of matrix are
conducted as presented in Fig. 7.
The results are shown in Fig. 8. It is observed that with the increas-
ing of the volume fraction of matrix, absolute values of most of effective
parameters decrease. This is due to the fact that matrix (epoxy) is softer
than inclusion (carbon). It should be emphasized that when the volume
fraction of matrix is 1, namely the material is purely homogeneous, the
higher order parameters vanish as expected.
Further investigations are carried out for RVEs by varying their
sizes (1mm ×1 mm ×1 mm and 2mm ×2 mm ×2 mm as well
as 3mm ×3 mm ×3 mm) as shown in Figs. 9 and 10. It is found
that all coefficients remain constant, which indicates that the obtained
parameters are independent of the repetition of RVEs.
International Journal of Solids and Structures 238 (2022) 111386
7
H. Yang et al.
Fig. 5. Visualization regarding the meaning of the homothetic ratio 𝜖.
Table 7
Voigt notation used for 3D strain-gradient tensors.
𝛼1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
𝑖𝑗𝑘 111 221 122 331 133 222 112 121 332 233 333 113 131 223 232 231 132 123
Fig. 6. The structure of strain gradient stiffness tensor for transverse isotropic materials. It is found that the first two 5 ×5 matrices in the diagonal are equal, for example,
𝐷111111 =𝐷222222. In the third 5 ×5 matrix in the diagonal, it is also observed that 𝐷333113 =𝐷333223,𝐷333131 =𝐷333232 ,𝐷113113 =𝐷223223,𝐷113131 =𝐷223232 ,𝐷113232 =𝐷131223,
𝐷131131 =𝐷232232. In the 3 ×3 matrix, 𝐷231231 =𝐷132132 ,𝐷231123 =𝐷132123.
International Journal of Solids and Structures 238 (2022) 111386
8
H. Yang et al.
𝐷M
𝛼𝛽 =
⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝
506.2 180.1 −178.8 213.5 17.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
180.1 −297.1 −168.8 −11.4 −93.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
−178.8 −168.8 180.3 −100.6 −64.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
213.5 −11.4 −100.6 −321.4 −284.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
17.3 −93.9 −64.7 −284.1 55.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 506.5 180.2 −178.8 213.6 16.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 180.2 −297.4 −169.0 −11.5 −94.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 −178.8 −169.0 180.2 −100.6 −64.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 213.6 −11.5 −100.6 −322.0 −283.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 16.9 −94.0 −64.7 −283.8 55.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 164.1 4.0 −207.8 4.0 −207.7 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 5.9 39.6 −4.9 −47.3 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −207.8 39.6 181.9 −47.3 −126.9 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 −4.9 −47.3 6.2 39.3 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −207.7 −47.3 −126.9 39.3 182.1 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −124.2 −143.1 −67.6
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −143.1 −124.5 −67.6
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −67.6 −67.6 22.3
⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠
N.
Box I.
𝐷M
𝛼𝛽 =
⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝
𝑑1𝑑2𝑑3𝑑4𝑑50 0 0 0 0 0 0 0 0 0 0 0 0
𝑑6𝑑7𝑑8𝑑90 0 0 0 0 0 0 0 0 0 0 0 0
𝑑10 𝑑11 𝑑12 0 0 0 0 0 0 0 0 0 0 0 0 0
𝑑13 𝑑14 0 0 0 0 0 0 0 0 0 0 0 0 0
𝑑15 0 0 0 0 0 0 0 0 0 0 0 0 0
𝑑1𝑑2𝑑3𝑑4𝑑50 0 0 0 0 0 0 0
𝑑6𝑑7𝑑8𝑑90 0 0 0 0 0 0 0
𝑑10 𝑑11 𝑑12 0 0 0 0 0 0 0 0
𝑑13 𝑑14 0 0 0 0 0 0 0 0
𝑑15 0 0 0 0 0 0 0 0
𝑑16 𝑑17 𝑑18 𝑑17 𝑑18 0 0 0
𝑑19 𝑑20 𝑑21 𝑑22 0 0 0
𝑑23 𝑑22 𝑑24 0 0 0
𝑑19 𝑑20 0 0 0
𝑑23 0 0 0
𝑑25 𝑑27 𝑑28
𝑑25 𝑑28
Sym. 𝑑26
⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠
.
Box II.
Fig. 7. Different volume fraction of matrix. 𝑙= 1 mm,𝑟1= 0.45 mm,𝑟2= 0.35 mm.
Effective parameters are studied for unit cells with varying sizes as
displayed in Fig. 11. The smaller unit cells are generated by hypothet-
ically scaling the larger one. Therefore, the volume fraction of matrix
is identical in these cases. It is found in Fig. 12 that the parameters in
the classical stiffness tensor remain the same, but the ones in the strain
gradient stiffness tensor vary by changing the unit cell lengths. This
fact is because of 𝐶M
𝑖𝑗𝑘𝑙 being invariant regarding the microstructural
size. However the effective strain gradient ones are sensitive to the
homothetic ratio 𝜖. These higher order parameters follow a scaling rule.
For example, the parameters can be obtained for the unit cell size of
0.5 mm × 0.5 mm × 0.5 mm by multiplying a scaling factor with the
effective parameters of the unit cell size of 1mm×1mm×1mm. The
International Journal of Solids and Structures 238 (2022) 111386
9
H. Yang et al.
Fig. 8. Effective material parameters with the changing of volume fraction of matrix. It should be noted that when the material is purely homogeneous (volume fraction of matrix
is 1), all higher order parameters vanish.
Fig. 9. RVEs constructed by 1 unit cell, 8 unit cells, 27 unit cells. 𝑙=1 mm, the radius of fiber is 0.45 mm.
scaling factor is the square of homothetic ratio 𝜖2, which is numerically
equal to the square of ratio of the unit cell lengths herein.
4.3.2. SiC/Al Metal Matrix Composite (MMC)
Aluminum-based MMCs have gained interest in engineering over
the past three decades. The insertion of a ceramic material into an
aluminum matrix leads to high stiffness and toughness of the compos-
ite material. In this section, the effective properties of SiC/Al metal
matrix composite are investigated. RVE models have been created for
three-dimensional spherical particles embedded into the metal matrix.
Filler is used as a reinforcement. Their volume ratios within the MMC
vary from 0% to 38.2% by volume. The material parameters taken
from Böhm et al. (2002) are compiled in Table 8.
Table 8
Material properties used for SiC/Al metal matrix composite material. 𝐸Young’s
modulus, 𝜈Poisson’s ratio, and 𝜌density.
Type 𝐸in GPa 𝜈 𝜌 in kg/m3
Matrix (Al2618-T4) 70 0.3 2900
Inclusion (SiC) 450 0.17 3100
The identified parameters for 61.8% volume ratio of matrix are
found as follows and as in Box III.
𝐶M
𝐴𝐵 =⎛⎜⎜⎜⎜⎜⎜⎝
163.3 50.5 50.5 0.0 0.0 0.0
50.5 163.5 50.5 0.0 0.0 0.0
50.5 50.5 163.6 0.0 0.0 0.0
0.0 0.0 0.0 46.4 0.0 0.0
0.0 0.0 0.0 0.0 46.3 0.0
0.0 0.0 0.0 0.0 0.0 46.3
⎞⎟⎟⎟⎟⎟⎟⎠
GPa,
International Journal of Solids and Structures 238 (2022) 111386
10
H. Yang et al.
Fig. 10. Effective material parameters with the repetition of RVEs (1mm ×1mm ×1 mm, 2mm ×2 mm ×2 mm, 3 mm ×3mm ×3 mm).
Fig. 11. Unit cells with the changing sizes. 𝑙=1 mm The volume fraction of matrix
are kept equal.
Effected by the cubic material symmetry of the RVE, there are three
independent parameters in the classical stiffness tensor,
𝐶M
𝐴𝐵 =⎛⎜⎜⎜⎜⎜⎜⎝
𝑐1𝑐2𝑐20 0 0
𝑐2𝑐1𝑐20 0 0
𝑐2𝑐2𝑐10 0 0
0 0 0 𝑐30 0
0 0 0 0 𝑐30
0 0 0 0 0 𝑐3
⎞⎟⎟⎟⎟⎟⎟⎠
.
As shown in Fig. 13, excluding the parameters of the same value,
there are 11 parameters found in the strain gradient stiffness tensor,
which is given in Box IV.
Please note, among these 11 parameters, some of them might be
linearly dependent. Further investigations are conducted for different
volume fraction of matrix, different sizes of selected RVE, and different
sizes of unit cells as indicated in Figs. 14,16,18. Results are displayed
in Figs. 15,17,19. It is observed that the higher order parameters are
zero when materials are homogeneous; they are independent of the
stack of RVEs, and they are sensitive to microstructural sizes as well
as following the scaling rule.
4.3.3. Aluminum foam
Aluminum foam is a highly porous metallic material with a cellular
substructure. The RVEs of the aluminum foam are modeled by using
a cubic inclusion, which are literally voids embedded in a matrix
made of aluminum. In order to avoid numerical problems, a small
number is assigned to the Young’s modulus of voids. Indeed, this
benchmark case is challenging to obtain consistently by using other
procedures in the literature, where the mass density in the microscale
is exchanged with the volume averaged mass density. The parameters
do not show a monotonous convergence. Herein, the change in the
formulation solves the problem by considering a distinction between
mass densities leading to correct results in 𝝍and thus in 𝑮and 𝑫pa-
rameters. The material properties used for aluminum foam is found in
Table 9.
International Journal of Solids and Structures 238 (2022) 111386
11
H. Yang et al.
Fig. 12. Effective material parameters with the changing lengths of unit cells, we emphasize that the substructure remains the same.
Fig. 13. The structure of the strain gradient stiffness tensor for cubic materials. Three 5 ×5 matrices in the diagonal are equal, For example, 𝐷111111 =𝐷222222. In the first 5 ×5
matrix, it is also observed that 𝐷111221 =𝐷111331,𝐷111122 =𝐷111133 ,𝐷221221 =𝐷331331,𝐷122122 =𝐷133133 ,𝐷221122 =𝐷331133,𝐷221133 =𝐷122331. In the 3 ×3 matrix, 𝐷231231 =𝐷132132 =𝐷123123,
𝐷231132 =𝐷231123 =𝐷132123.
Fig. 14. Changing volume fraction of matrix.
International Journal of Solids and Structures 238 (2022) 111386
12
H. Yang et al.
𝐷M
𝛼𝛽 =
⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝
7120.6 1075.0 −844.4 1077.8 −836.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1075.0 −2517.5 −788.5 −48.3 −275.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
−844.4 −788.5 1914.3 −276.5 −347.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1077.9 −48.3 −276.5 −2515.1 −789.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
−836.6 −275.8 −347.3 −789.7 1915.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 7130.6 1070.6 −850.1 1076.5 −840.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1070.6 −2504.4 −781.2 −48.1 −276.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 −850.1 −781.2 1915.4 −275.8 −347.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1076.5 −48.1 −275.8 −2513.8 −790.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 −840.1 −276.2 −347.6 −790.8 1917.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7141.0 1072.0 −848.8 1072.0 −849.8 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1072.0 −2490.1 −783.5 −47.5 −276.5 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −848.8 −783.5 1915.3 −275.4 −350.4 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1072.0 −47.5 −275.4 −2506.2 −786.2 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −849.8 −276.5 −350.4 −786.2 1915.4 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −596.3 −711.0 −709.6
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −711.0 −593.9 −708.1
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 −709.6 −708.1 −589.7
⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠
N.
Box III.
𝐷M
𝛼𝛽 =
⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝
𝑑1𝑑2𝑑3𝑑2𝑑30 0 0 0 0 0 0 0 0 0 0 0 0
𝑑4𝑑5𝑑6𝑑70 0 0 0 0 0 0 0 0 0 0 0 0
𝑑8𝑑7𝑑90 0 0 0 0 0 0 0 0 0 0 0 0
𝑑4𝑑50 0 0 0 0 0 0 0 0 0 0 0 0
𝑑80 0 0 0 0 0 0 0 0 0 0 0 0
𝑑1𝑑2𝑑3𝑑2𝑑30 0 0 0 0 0 0 0
𝑑4𝑑5𝑑6𝑑70 0 0 0 0 0 0 0
𝑑8𝑑7𝑑90 0 0 0 0 0 0 0
𝑑4𝑑50 0 0 0 0 0 0 0
𝑑80 0 0 0 0 0 0 0
𝑑1𝑑2𝑑3𝑑2𝑑30 0 0
𝑑4𝑑5𝑑6𝑑70 0 0
𝑑8𝑑7𝑑90 0 0
𝑑4𝑑50 0 0
𝑑80 0 0
𝑑10 𝑑11 𝑑11
𝑑10 𝑑11
Sym. 𝑑10
⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠
.
Box IV.
Table 9
Material properties used for aluminum foam. 𝐸Young’s modulus, 𝜈Poisson’s ratio, and
𝜌mass density.
Type 𝐸in GPa 𝜈 𝜌 in kg/m3
Matrix (Aluminum) 70 0.3 2700
Inclusion (Voids) 10−10 0.0 0.0
The identified parameters are found as follows and as given in
Box V.
𝐶M
𝐴𝐵 =⎛⎜⎜⎜⎜⎜⎜⎝
15.1 3.0 3.0 0.0 0.0 0.0
3.0 15.1 3.0 0.0 0.0 0.0
3.0 3.0 15.1 0.0 0.0 0.0
0.0 0.0 0.0 2.9 0.0 0.0
0.0 0.0 0.0 0.0 2.9 0.0
0.0 0.0 0.0 0.0 0.0 2.9
⎞⎟⎟⎟⎟⎟⎟⎠
GPa,
Three independent parameters and eleven parameters are observed
in the classical stiffness tensor and the strain gradient stiffness tensor,
respectively. This is consistent to the cubic material symmetry as
mentioned before. Investigations on the different volume fraction of
matrix, repetition of RVEs, changing sizes of unit cells are conducted
as displayed in Figs. 20,22,24. Corresponding outcomes are presented
in Figs. 21,23,25.
5. Remark on positive definiteness
As observed from the previous sections, negative values appear in
strain gradient stiffness tensors. This fact may raise concerns regarding
the positive definiteness of the strain energy function. In Mindlin and
Eshel (1968), dell’Isola et al. (2009), Nazarenko et al. (2020) and
Eremeyev et al. (2020), the issue of positive definiteness of the strain
energy function for strain gradient materials is addressed and bounds
on material parameters are provided. The bounds on strain gradient
constants consider the continuum to be purely local, which means that
the strain energy function is convex with respect to every material
point (Mindlin and Eshel,1968). However, when homogenizing the mi-
crostructures of composite materials with an equivalent strain gradient
International Journal of Solids and Structures 238 (2022) 111386
13
H. Yang et al.
Fig. 15. Effective material parameters with the changing of volume fraction of matrix.
Fig. 16. RVEs constructed by 1 unit cell, 8 unit cells, 27 unit cells.
𝐷M
𝛼𝛽 =
⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝
1130.3 185.4 288.8 184.8 288.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
185.4 1080.6 114.9 328.0 74.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
288.8 114.9 −42.8 74.5 160.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
184.8 328.0 74.5 1080.3 114.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
288.6 74.6 160.7 114.9 −42.6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1139.1 187.6 290.6 186.9 290.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 187.6 1081.0 114.7 328.4 75.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 290.6 114.7 −42.8 74.9 161.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 186.9 328.4 74.9 1080.4 115.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 290.2 75.0 161.0 115.4 −42.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1171.8 194.3 296.5 194.5 296.9 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 194.3 1082.3 115.6 329.8 76.3 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 296.5 115.6 −42.1 76.3 162.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 194.5 329.8 76.3 1082.1 115.8 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 296.9 76.3 162.0 115.8 −42.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 406.8 19.6 19.9
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 19.6 406.7 19.8
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 19.9 19.8 406.9
⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠
N.
Box V.
International Journal of Solids and Structures 238 (2022) 111386
14
H. Yang et al.
Fig. 17. Effective material parameters with the changing RVE sizes (1 ×1×1, 2 ×2×2, 3 ×3×3).
Fig. 18. Unit cells with changing lengths.
Fig. 19. Effective material parameters with the changing lengths of unit cells.
International Journal of Solids and Structures 238 (2022) 111386
15
H. Yang et al.
Fig. 20. Different volume fraction of matrix for the aluminum foam.
Fig. 21. Effective material parameters with changing of volume fraction of matrix.
Fig. 22. RVEs constructed by 1 unit cell, 8 unit cells, 27 unit cells.
continuum, we have a limited non-locality. The non-locality originates
from the energy equivalence as shown in Eq. (1). We emphasize that
the 𝜖is a finite number, 𝜖 < 1but not necessarily 𝜖 ≪ 1, which means
that the studied composite material has a finite macroscopic and micro-
scopic sizes. Therefore, the strain energy function averaged over this
microstructure size should be positive definite and not the pointwise
local strain energy function. Thus coefficients in the strain gradient
stiffness tensor could be negative as long as the strain energy density
function integrated over the periodic unit cell is positive definite. This
interpretation is aligned with in Kumar and McDowell (2004), Barboura
and Li (2018) and Li and Zhang (2013). Additionally, this condition is
always fulfilled if the microstructure energy density is positive definite,
since Eq. (1) is enforced.
6. Verification of the homogenized strain gradient models
In order to assess the homogenized strain gradient continuum model
developed in this paper, finite element computations are conducted
to evaluate the performance of the proposed model. To this end, a
cantilever beam bending problem is selected as presented in Fig. 26.
The beam is made out of aluminum foam with periodically aligned
International Journal of Solids and Structures 238 (2022) 111386
16
H. Yang et al.
Fig. 23. Effective material parameters with changing RVE sizes (1 ×1×1, 2 ×2×2, 3 ×3×3).
Fig. 24. Unit cells with changing lengths.
microstructures of 1 mm ×1 mm ×1 mm. As mentioned above, the
inclusions of the microstructures are voids. We use the determined
parameters with matrix volume fraction of 0.271 as shown in Fig. 21.
The length, width, and height of the beam are assigned to be 50 mm,
2 mm, and 2 mm. The left surface, 𝑋1= 0, of the beam is clamped
(𝑢1=𝑢2=𝑢3= 0). A traction is applied at the right surface 𝑋1=𝐿.
We conduct three simulations. A Direct Numerical Simulation
(DNS), where the microstructure is modeled in detail. This result is
accepted as correct. A homogenization simulation with second order
(strain gradient) theory, where 𝑫and 𝑮tensors are employed. A
homogenized simulation with first order theory, in other words, 𝑫
and 𝑮are set to zero (first order theory is used). As 𝐶1continuity
is required for the numerical implementation of the strain gradient
computations, the isogeometric analysis is used in the simulations for
the homogenized models. The codes developed and verified in Yang
et al. (2021) are used herein. The weak form for linear elastic strain
gradient materials is presented as
∫𝛺(𝜎𝑖𝑗 𝛿𝑢𝑖,𝑗 +𝜏𝑖𝑗𝑘𝛿𝑢𝑖,𝑗𝑘)d𝑉=∫𝜕𝛺
𝑡𝑖𝛿𝑢𝑖d𝐴, (30)
where 𝜎𝑖𝑗 and 𝜏𝑖𝑗𝑘 are the stress tensor and hyperstress tensor defined
by
𝜎𝑖𝑗 =𝜕𝑤M
𝜕𝑢𝑖,𝑗
, 𝜏𝑖𝑗𝑘 =𝜕𝑤M
𝜕𝑢𝑖,𝑗𝑘
.(31)
with 𝑤Mthe strain energy density. The body forces, double traction,
and the so-called wedge forces are all set to be zero, therefore the
pertinent terms in the weak form Eq. (30) are neglected. The traction is
applied incrementally from (0,0,0) to (0,0,0.001 MPa). The calculated
results of total displacement are shown in Fig. 27. Strain gradient
results match accurately DNS results. However, a significant deviation
from DNS is observed, if one uses first order theory. DNS shows that
the beam ‘‘act’’ stiffer than first order theory suggests, this experimen-
tally well known fact is called size effect. The difference vanishes as
homothetic ratio approaches zero, in other words, the same foam in a
larger beam shows no size effect. Such phenomenon is also observed
in the numerical investigation in Abali et al. (2015,2017) In order to
assess the models further, more investigations are made as displayed in
Figs. 28 and 29.
The elapsed time for DNS is 771.5 s by using a computer (In-
tel(R) Core(TM) i7-8565U CPU). The elapsed time for the homogenized
Cauchy model and the strain gradient model is 15.0 s and 59.7 s as
shown in Fig. 28(a). It is evident that by using the homogenization
techniques, computational efficiency is greatly improved. Additional
difficulty is the challenge for the meshing algorithm to construct a
high quality mesh for the microstructure. Element quality will not
be ensured around sharp edges leading to inconsistencies as well as
strain concentrations making the model error-prone. On the other
hand, a homogeneous structure faces only macroscopic sharp contours
such that a mesh convergence is feasible to minimize the numerical
errors. Indeed, the computational efficiency of the first order theory
is significantly larger than the second order theory, even in the same
type of mesh. But for the chosen homothetic ratio herein, we obtain an
inadequate result from the first order theory as indicated in Figs. 28(b)
and 29.
7. Conclusions
Asymptotic homogenization method has been employed to homog-
enize composite material into effective homogeneous strain gradient
continua. Main conclusions are summarized as follows:
•Purely computational analysis determines all the parameters in
the strain gradient theory. In particular, the parameters in the
rank five tensor and the rank six tensor.
•Numerical examples for 2D and 3D, stiff and soft inclusions, cubic
and transverse material symmetry cases have been conducted.
International Journal of Solids and Structures 238 (2022) 111386
17
H. Yang et al.
Fig. 25. Effective material parameters with changing lengths of unit cells.
Fig. 26. Schematic of a cantilever beam bending problem. The length, width, and height of the beam is 𝐿,𝑊,𝐻. The inclusions (voids) are presented as red and matrix is
indicated as blue (with less opacity for the sake of visualization). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of
this article.)
Fig. 27. Comparisons of total displacement among the heterogeneous Cauchy continuum, homogenized strain gradient continuum, and the homogenized Cauchy continuum in the
case of cantilever beam bending. Scaling factor 5000.
•In both 2D and 3D numerical examples, the effective strain gradi-
ent parameters vanish when materials are purely homogeneous,
they are independent of repetitions of RVEs and sensitive to
microstructural sizes.
•Without assuming a specific symmetry group, in the case of
cubic symmetry, all expected relations have been captured by the
proposed formalism.
•Physical meaning of the homothetic ratio 𝜖is interpreted, a so-
called scaling rule for effective strain gradient parameters has
been discussed. The method is valid when 𝜖is a finite value. 𝜖 < 1
is required but not necessarily 𝜖 ≪ 1.
•An evaluation of the performance of the determined strain gra-
dient parameters is done in 3D for the first time. It is found that
including the strain gradient terms in the homogenized model will
International Journal of Solids and Structures 238 (2022) 111386
18
H. Yang et al.
Fig. 28. Comparisons of the elapsed time, the relative errors of the strain energy and maximum displacement for the first order theory and second order theory results.
Fig. 29. Comparisons of the calculated strain energy and maximum displacement among DNS, first order theory, and second order theory results.
improve the accuracy of the prediction of the response of alu-
minum foams compared to classical (first order) homogenization.
The homogenization tool is applicable to any composite materials
with a periodic substructure at the microscale. Such multiscale are
nowadays possible to manufacture by 3D printers. Therefore, effec-
tive parameters determination is of interest for a possible topology
optimization. Further investigations will focus on the following aspects:
•To validate the identified parameters not only in statics but also in
vibration responses, buckling critical loads (Khakalo et al.,2018),
and wave propagation (Rosi and Auffray,2016;Eremeyev et al.,
2019).
•To apply the homogenization method to the analysis of 3D com-
posite materials with finite thickness. This may be achieved by
modeling the full thickness unit cell model and relieving the
out-of-plane periodicity of the unit cell (Nasution et al.,2014).
•To explore the possibility of studying more sophisticated meta-
materials such as the so-called pantographic structures (dell’Isola
et al.,2019a,b) or the biomimetic spinodoids metamaterials
(Portela et al.,2020) by the homogenization method.
•To extend the homogenization method to non-linear regime (For-
est,2020;ElNady et al.,2016) and multiphysics fields.
Declaration of competing interest
The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.
References
Abali, B.E., 2017. Computational Reality. In: Advanced Structured Materials, vol. 55,
Springer Nature, Singapore.
Abali, B.E., 2020. Supply code for computations. http://bilenemek.abali.org/.
Abali, B.E., Barchiesi, E., 2021. Additive manufacturing introduced substructure
and computational determination of metamaterials parameters by means of the
asymptotic homogenization. Contin. Mech. Thermodyn. 33 (4), 993–1009.
Abali, B.E., Müller, W.H., dell’Isola, F., 2017. Theory and computation of higher
gradient elasticity theories based on action principles. Arch. Appl. Mech. 87 (9),
1495–1510.
Abali, B.E., Müller, W.H., Eremeyev, V.A., 2015. Strain gradient elasticity with
geometric nonlinearities and its computational evaluation. Mech. Adv. Mater. Mod.
Process. 1 (1), 1–11.
Abdoul-Anziz, H., Seppecher, P., Bellis, C., 2019. Homogenization of frame lattices
leading to second gradient models coupling classical strain and strain-gradient
terms. Math. Mech. Solids 24 (12), 3976–3999.
Alibert, J.-J., Della Corte, A., 2019. Homogenization of nonlinear inextensible
pantographic structures by 𝛤-convergence. Math. Mech. Complex Syst. 7 (1), 1–24.
Altenbach, H., Forest, S., 2016. Generalized Continua as Models for Classical and
Advanced Materials. Springer.
Arabnejad, S., Pasini, D., 2013. Mechanical properties of lattice materials via asymptotic
homogenization and comparison with alternative homogenization methods. Int. J.
Mech. Sci. 77, 249–262.
Auffray, N., Bouchet, R., Brechet, Y., 2009. Derivation of anisotropic matrix for
bi-dimensional strain-gradient elasticity behavior. Int. J. Solids Struct. 46 (2),
440–454.
Auffray, N., Dirrenberger, J., Rosi, G., 2015. A complete description of bi-dimensional
anisotropic strain-gradient elasticity. Int. J. Solids Struct. 69, 195–206.
Auffray, N., Le Quang, H., He, Q.-C., 2013. Matrix representations for 3D strain-gradient
elasticity. J. Mech. Phys. Solids 61 (5), 1202–1223.
Bacigalupo, A., Paggi, M., Dal Corso, F., Bigoni, D., 2018. Identification of higher-order
continua equivalent to a Cauchy elastic composite. Mech. Res. Commun. 93, 11–22.
Barboura, S., Li, J., 2018. Establishment of strain gradient constitutive relations by
using asymptotic analysis and the finite element method for complex periodic
microstructures. Int. J. Solids Struct. 136, 60–76.
Bleyer, J., 2018. Numerical Tours of Computational Mechanics with FEniCS. Zenodo,
http://dx.doi.org/10.5281/zenodo.1287832.
Böhm, H.J., Eckschlager, A., Han, W., 2002. Multi-inclusion unit cell models for metal
matrix composites with randomly oriented discontinuous reinforcements. Comput.
Mater. Sci. 25 (1–2), 42–53.
Boutin, C., 1996. Microstructural effects in elastic composites. Int. J. Solids Struct. 33
(7), 1023–1051.
International Journal of Solids and Structures 238 (2022) 111386
19
H. Yang et al.
Boutin, C., dell’Isola, F., Giorgio, I., Placidi, L., 2017. Linear pantographic sheets:
Asymptotic micro-macro models identification. Math. Mech. Complex Syst. 5 (2),
127–162.
Chen, Q., Chatzigeorgiou, G., Meraghni, F., 2020. Extended mean-field homogeniza-
tion of viscoelastic-viscoplastic polymer composites undergoing hybrid progressive
degradation induced by interface debonding and matrix ductile damage. Int. J.
Solids Struct. 210, 1–17.
dell’Isola, F., Sciarra, G., Vidoli, S., 2009. Generalized Hooke’s law for isotropic second
gradient materials. Proc. R. Soc. A 465 (2107), 2177–2196.
dell’Isola, F., Seppecher, P., Alibert, J.J., Lekszycki, T., Grygoruk, R., Pawlikowski, M.,
Steigmann, D., Giorgio, I., Andreaus, U., Turco, E., Gołaszewski, M., Rizzi, N.,
Boutin, C., Eremeyev, V.A., Misra, A., Placidi, L., Barchiesi, E., Greco, L.,
Cuomo, M., Cazzani, A., Corte, A.D., Battista, A., Scerrato, D., Eremeeva, I.Z.,
Rahali, Y., cois Ganghoffer, J.-F., Müller, W., Ganzosch, G., Spagnuolo, M., Pfaff, A.,
Barcz, K., Hoschke, K., Neggers, J., Hild, F., 2019b. Pantographic metamaterials:
an example of mathematically driven design and of its technological challenges.
Contin. Mech. Thermodyn. 31 (4), 851–884.
dell’Isola, F., Seppecher, P., Spagnuolo, M., Barchiesi, E., Hild, F., Lekszycki, T., Gior-
gio, I., Placidi, L., Andreaus, U., Cuomo, M., Eugster, S.R., Pfaff, A., Hoschke, K.,
Langkemper, R., Turco, E., Sarikaya, R., Misra, A., Angelo, M.D., D’Annibale, F.,
Bouterf, A., Pinelli, X., Misra, A., Desmorat, B., Pawlikowski, M., Dupuy, C.,
Scerrato, D., Peyre, P., Laudato, M., Manzari, L., Göransson, P., Hesch, C., Hesch, S.,
Franciosi, P., Dirrenberger, J., Maurin, F., Vangelatos, Z., Grigoropoulos, C.,
Melissinaki, V., Farsari, M., Müller, W., Abali, B.E., Liebold, C., Ganzosch, G.,
Harrison, P., Drobnicki, R., Igumnov, L., Alzahrani, F., Hayat, T., 2019a. Advances
in pantographic structures: design, manufacturing, models, experiments and image
analyses. Contin. Mech. Thermodyn. 31 (4), 1231–1282.
Dirrenberger, J., Forest, S., Jeulin, D., 2019. Computational homogenization of architec-
tured materials. In: Estrin, Y., Bréchet, Y., Dunlop, J., Fratzl, P. (Eds.), Architectured
Materials in Nature and Engineering. In: Springer Series in Materials Science,
Springer, pp. 89–139.
Dos Reis, F., Ganghoffer, J., 2012. Construction of micropolar continua from the
asymptotic homogenization of beam lattices. Comput. Struct. 112, 354–363.
ElNady, K., Goda, I., Ganghoffer, J.-F., 2016. Computation of the effective nonlinear
mechanical response of lattice materials considering geometrical nonlinearities.
Comput. Mech. 58 (6), 957–979.
Eremeyev, V.A., 2016. On effective properties of materials at the nano-and microscales
considering surface effects. Acta Mech. 227 (1), 29–42.
Eremeyev, V.A., Lurie, S.A., Solyaev, Y.O., dell’Isola, F., 2020. On the well posedness
of static boundary value problem within the linear dilatational strain gradient
elasticity. Z. Angew. Math. Phys. 71 (6), 1–16.
Eremeyev, V.A., Rosi, G., Naili, S., 2019. Comparison of anti-plane surface waves in
strain-gradient materials and materials with surface stresses. Math. Mech. Solids
24 (8), 2526–2535.
Eringen, A.C., 1999. Theory of micropolar elasticity. In: Microcontinuum Field Theories.
Springer, pp. 101–248.
Forest, S., 2020. Continuum thermomechanics of nonlinear micromorphic, strain and
stress gradient media. Phil. Trans. R. Soc. A 378 (2170), 20190169.
Forest, S., Trinh, D.K., 2011. Generalized continua and non-homogeneous boundary
conditions in homogenisation methods. ZAMM-J. Appl. Math. Mech./Z. Angew.
Math. Mech. 91 (2), 90–109.
Ganghoffer, J., Reda, H., 2021. A variational approach of homogenization of het-
erogeneous materials towards second gradient continua. Mech. Mater. 158,
103743.
Goda, I., Ganghoffer, J.-F., 2016. Construction of first and second order grade
anisotropic continuum media for 3D porous and textile composite structures.
Compos. Struct. 141, 292–327.
Hollister, S.J., Kikuchi, N., 1992. A comparison of homogenization and standard
mechanics analyses for periodic porous composites. Comput. Mech. 10 (2), 73–95.
Jakabčin, L., Seppecher, P., 2020. On periodic homogenization of highly contrasted
elastic structures. J. Mech. Phys. Solids 144, 104104.
Khakalo, S., Balobanov, V., Niiranen, J., 2018. Modelling size-dependent bending,
buckling and vibrations of 2D triangular lattices by strain gradient elasticity models:
applications to sandwich beams and auxetics. Internat. J. Engrg. Sci. 127, 33–52.
Kouznetsova, V., Geers, M.G., Brekelmans, W.M., 2002. Multi-scale constitutive
modelling of heterogeneous materials with a gradient-enhanced computational
homogenization scheme. Internat. J. Numer. Methods Engrg. 54 (8), 1235–1260.
Kouznetsova, V., Geers, M.G., Brekelmans, W., 2004. Multi-scale second-order compu-
tational homogenization of multi-phase materials: a nested finite element solution
strategy. Comput. Methods Appl. Mech. Engrg. 193 (48–51), 5525–5550.
Kumar, R.S., McDowell, D.L., 2004. Generalized continuum modeling of 2-D periodic
cellular solids. Int. J. Solids Struct. 41 (26), 7399–7422.
Li, J., 2011. A micromechanics-based strain gradient damage model for fracture
prediction of brittle materials–Part I: Homogenization methodology and constitutive
relations. Int. J. Solids Struct. 48 (24), 3336–3345.
Li, J., Zhang, X.-B., 2013. A numerical approach for the establishment of strain gradient
constitutive relations in periodic heterogeneous materials. Eur. J. Mech. A Solids
41, 70–85.
Liu, S., Su, W., 2009. Effective couple-stress continuum model of cellular solids and
size effects analysis. Int. J. Solids Struct. 46 (14–15), 2787–2799.
Mandadapu, K.K., Abali, B.E., Papadopoulos, P., 2021. On the polar nature and
invariance properties of a thermomechanical theory for continuum-on-continuum
homogenization. Math. Mech. Solids 26 (11), 1581–1598.
Mindlin, R.D., Eshel, N., 1968. On first strain-gradient theories in linear elasticity. Int.
J. Solids Struct. 4 (1), 109–124.
Misra, A., Poorsolhjouy, P., 2015. Identification of higher-order elastic constants for
grain assemblies based upon granular micromechanics. Math. Mech. Complex Syst.
3 (3), 285–308.
Müller, W.H., 2020. The experimental evidence for higher gradient theories. In:
Bertram, A., Forest, S. (Eds.), Mechanics of Strain Gradient Materials. In: CISM
International Centre for Mechanical Sciences, vol. 600, Springer, pp. 1–18.
Nasution, M.R.E., Watanabe, N., Kondo, A., Yudhanto, A., 2014. A novel asymptotic
expansion homogenization analysis for 3-D composite with relieved periodicity in
the thickness direction. Compos. Sci. Technol. 97, 63–73.
Nazarenko, L., Glüge, R., Altenbach, H., 2020. Positive definiteness in coupled strain
gradient elasticity. Contin. Mech. Thermodyn. 1–13.
Portela, C.M., Vidyasagar, A., Krödel, S., Weissenbach, T., Yee, D.W., Greer, J.R.,
Kochmann, D.M., 2020. Extreme mechanical resilience of self-assembled
nanolabyrinthine materials. Proc. Natl. Acad. Sci. 117 (11), 5686–5693.
Rahali, Y., Giorgio, I., Ganghoffer, J., dell’Isola, F., 2015. Homogenization à la
Piola produces second gradient continuum models for linear pantographic lattices.
Internat. J. Engrg. Sci. 97, 148–172.
Rokoš, O., Ameen, M.M., Peerlings, R.H., Geers, M.G., 2019. Micromorphic compu-
tational homogenization for mechanical metamaterials with patterning fluctuation
fields. J. Mech. Phys. Solids 123, 119–137.
Rosi, G., 2020. Waves and generalized continua. In: Altenbach, H., Öchsner, A. (Eds.),
Encyclopedia of Continuum Mechanics. Springer, pp. 2756–2765.
Rosi, G., Auffray, N., 2016. Anisotropic and dispersive wave propagation within
strain-gradient framework. Wave Motion 63, 120–134.
Rosi, G., Placidi, L., Auffray, N., 2018. On the validity range of strain-gradient
elasticity: a mixed static-dynamic identification procedure. Eur. J. Mech. A Solids
69, 179–191.
Skrzat, A., Eremeyev, V.A., 2020. On the effective properties of foams in the framework
of the couple stress theory. Contin. Mech. Thermodyn. 1–23.
Weeger, O., 2021. Numerical homogenization of second gradient, linear elastic
constitutive models for cubic 3D beam-lattice metamaterials. Int. J. Solids Struct..
Yang, H., Abali, B.E., Timofeev, D., Müller, W.H., 2020. Determination of metamaterial
parameters by means of a homogenization approach based on asymptotic analysis.
Contin. Mech. Thermodyn. 32 (5), 1251–1270.
Yang, H., Timofeev, D., Abali, B.E., Li, B., Müller, W.H., 2021. Verification of strain
gradient elasticity computation by analytical solutions. ZAMM-J. Appl. Math.
Mech./Z. Angew. Math. Mech. 101 (12), e202100023.
Yvonnet, J., Auffray, N., Monchiet, V., 2020. Computational second-order homogeniza-
tion of materials with effective anisotropic strain-gradient behavior. Int. J. Solids
Struct. 191, 434–448.
