<figure><img src="image_1.png" /></figure>

HAL
open science

# Multi-scale modeling and simulation of high-contrast periodic composite materials: second-order gradient theory

Alioune Nacro, Philippe Karamian, Sophie Lemaître

► To cite this version:

Alioune Nacro, Philippe Karamian, Sophie Lemaître. Multi-scale modeling and simulation of high-contrast periodic composite materials: second-order gradient theory. 2025. hal-05108135v2

HAL Id: hal-05108135

https://hal.science/hal-05108135v2

Preprint submitted on 10 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Alioune Nacro, Philippe Karamian-Surville and Sophie Lemaitre

Multi-scale modeling and simulation of high-contrast periodic composite materials: second-order gradient theory

Volume 353 (2025), p. 815-862

Online since: 10 July 2025

https://doi.org/10.5802/crmeca.307

<figure><img src="image_3.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# Multi-scale modeling and simulation of high-contrast periodic composite materials: second-order gradient theory

## Modélisation et simulation multi-échelle des matériaux composites périodiques à fort contraste : théorie du gradient d'ordre deux

Alioune Nacro<sup>a</sup>, Philippe Karamian-Surville<sup>®, * , a</sup> and Sophie Lemaitre<sup>a</sup>

**Abstract.** Within the framework of second-order gradient theory and the multi-scale modeling approach for periodic composite materials, a key challenge lies in determining the four tensors: $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$, which appear in the asymptotic expansion of the energy. This paper presents the numerical evaluation and simulation of these tensors across various geometric configurations, employing a modified Green's kernel-accelerated scheme to efficiently handle high-contrast cases. Extensive 3D simulations are conducted for multiple distinct geometries, accompanied by a detailed analysis of tensor values, computational efficiency, and classification. The results provide valuable insights into the behavior of these tensors across different morphological structures, contributing to a deeper understanding of advanced composite material modeling.

**Résumé.** Dans le cadre de la théorie du gradient d'ordre deux et de l'approche de modélisation multi-échelle des matériaux composites périodiques, un défi majeur réside dans la détermination des quatre tenseurs : $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, et $D^{0,0}$ qui apparaissent dans le développement asymptotique de l'énergie. Cet article présente l'évaluation numérique et la simulation de ces tenseurs pour différentes configurations géométriques, en utilisant un schéma modifié accéléré par le noyau de Green pour traiter efficacement les cas à fort contraste. De nombreuses simulations 3D sont menées pour plusieurs géométries distinctes, accompagnées d'une analyse détaillée des valeurs des tenseurs, de l'efficacité numérique et d'une classification. Les résultats fournissent des informations précieuses sur le comportement de ces tenseurs selon différentes structures morphologiques, contribuant ainsi à une meilleure compréhension de la modélisation avancée des matériaux composites.

**Keywords.** Second-order gradient, high contrast, periodic composite materials, multi-scale modeling, simulation, FFT, FEM.

**Mots-clés.** Gradient d'ordre deux, fort contraste, matériaux composites périodiques, modélisation multi-échelle, simulation, FFT, FEM.

Manuscript received 20 March 2025, revised and accepted 7 June 2025.

<footnote><sup>a</sup> Normandy University, UNICAEN, UMR 6139, CNRS, Nicolas Oresme Mathematics Laboratory, 14000 Caen, France</footnote>

<footnote>E-mails: alioune.nacro@unicaen.fr, philippe.karamian@unicaen.fr, sophie.lemaitre@unicaen.fr</footnote>

<footnote>* Corresponding author</footnote>
# 1. Introduction

As is well established, composite materials are defined as materials composed of at least two distinct constituents. They are widely utilized across various industrial sectors, including automotive, aerospace, construction, and sports. Over the past few decades, significant research has been conducted within the framework of first-order gradient theory, which encompasses two primary approaches. The first, the mean-field approach, has been extensively employed to estimate the effective properties of composite materials [1–9], utilizing analytical models such as the Mori–Tanaka method [10–12]. In contrast, the more accurate full-field approach relies on high-resolution numerical simulations, including the finite element method (FEM) and fast Fourier transform (FFT)-based methods. Both approaches aim to precisely characterize the effective properties of composite materials within the framework of first-order gradient theory.

In this study, the full-field approach is adopted for evaluating the effective properties of composite materials [13–16], employing the fast Fourier transform (FFT) method, with the finite element method (FEM) used for validation. Among the available numerical techniques, FFT-based methods have gained significant popularity over the past few decades due to their numerous advantages. Originally introduced by Moulinec and Suquet [16,17] in the context of computational homogenization [18–20] in the 1990s, these methods offer a computationally efficient alternative to traditional numerical approaches within a remarkably short computational time.

In particular, in this study, we re-examine the bibliographic foundations of the models presented to enhance the clarity of this paper. The original notations introduced by the respective authors have been retained to prevent any confusion with those associated with the model developed in this work. Non-local models naturally extend from the pioneering contributions of Kröner [21], Krumhansl [22], Eringen and Edelen [23], and Peerlings [24–27]. The theory of non-locality postulates that the response of a point within a material depends not only on its own deformation [28,29] but also on the deformation of neighboring points, potentially over long distances.

We classify non-local models into two main families: gradient-based models and those formulated using integral approaches. These so-called enriched models provide a more comprehensive description of material behavior by incorporating higher-order constitutive equations beyond those of the classical framework. However, the introduction of additional boundary conditions presents substantial implementation challenges, increasing the complexity of their numerical treatment.

In this paper, we focus on phenomenological approaches based on second-order gradient models, originally developed in the 1960s. Specifically, we examine the second-gradient model introduced by Toupin [30,31], Mindlin [32], and Green and Rivlin [28]. These authors demonstrated that the third-order tensor exhibits symmetry with respect to its first two indices and that the elastic energy is a quadratic function of $\epsilon(u^0)$ and $\partial^2 u^0$, expressed as:

$$W = \frac{1}{2} \epsilon(u^0) : \mathbb{L}^{0,0} : \epsilon(u^0) + \epsilon(\epsilon(u^0) : \mathbb{L}^{0,1} : \partial^2 u^0) + \epsilon^2(\partial^2 u^0 : \mathbb{L}^{1,1} : \partial^2 u^0),$$

where $\epsilon$ denotes the ratio between the characteristic length of the inclusion and that of the material.

According to the literature, similar to the model proposed by Toupin, Mindlin, Green, and Rivlin, we are interested in particular in the tensors $\mathbb{L}^{0,0}$ and $\mathbb{L}^{0,1}$ defined as follows using our notations:

$$\mathbb{L}^{0,0} = A^{0,0} + 2\epsilon B^{0,1},$$

$$\mathbb{L}^{0,1} = (C^{0,0} + D^{0,0}),$$

relationships examined in greater detail later in this article.
This article focuses on the application of second-order gradient theory to composite and periodic materials for the determination of these four key tensors: $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$ (refer to Appendix A for notation conventions) allowing a more precise and detailed analysis of material properties. One of the motivations of this work is the role of these second-order tensors in damage mechanics. Actually, several studies have proposed approaches for incorporating second-order gradient terms into damage and fatigue models for composite materials, including works by Pham [33,34], Marigo [35], and Benallal [36]. We will be explored damage and fatigue in a forthcoming study. Firstly, the aim of this study is to be able to calculate tensors numerically.

A comprehensive review of the literature [32,37–43] on multiscale modeling and simulation of high-contrast periodic composite materials, when integrated with second-order gradient theory, provides an in-depth exploration of the topic. This article examines various methodologies employed for tensor computation and discusses different tensor calculus techniques within the framework of second-order gradient theory, highlighting the associated challenges and opportunities.

Furthermore, it addresses the limitations of current research and proposes potential avenues for applying these techniques to enhance the understanding of composite behavior, thereby offering valuable insights for material design. The overarching objective is to model these composites as equivalent homogeneous materials, facilitating more accurate simulations within the framework of linear elasticity theory.

This article distinguishes itself through a multiscale approach that excludes energy theory [32,37–41] in determining the tensors $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$ [42,44] within the framework of second-order gradient theory for periodic composite materials. We demonstrate that the tensors $A^{0,0}$ and $B^{0,1}$ are identical. Furthermore, we introduce a numerical tool specifically designed for these computations, based on solving the Lippmann–Schwinger equation. The solution is obtained using the fast Fourier transform (FFT) in conjunction with an accelerated scheme algorithm and a modified Green’s kernel to effectively handle high-contrast cases [45,46]. Throughout this paper, *contrast* refers to the ratio between the Young’s modulus and Poisson’s coefficient of the inclusion and those of the matrix.

The article is structured as follows: the first section focuses on the theoretical aspects related to the calculation of the tensors $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$ (refer to Appendix A for the notation conventions), while the second section is dedicated to numerical simulations. The theoretical section is subdivided into two parts: the first part reviews the principles of homogenization and the multiscale approach, and the second part details the process of calculating the aforementioned tensors. The numerical section is divided into in two subsections. The first subsection presents a comparison between two numerical methods, namely the finite element method (FEM) and the fast Fourier transform (FFT), to validate the results. This comparison is based on different inclusion shapes (spherical and cylindrical) with a resolution of 128³. Subsequently, the FFT method is chosen for further simulations, 256³. In the next subsection, we presents the simulations for a representative volume element (RVE) with a resolution of 256³, divided into two parts: the first is devoted to inclusions immersed within the representative volume element (RVE), while the second addresses elongated inclusions that touch the edges of the RVE. The simulation includes various inclusion configurations: spherical, cylindrical, spherical/cylindrical, as well as sets composed of three cylindrical inclusions. The results related to the tensor $C^{0,0}$ reveal a classification of symmetry groups analogous to those of the tensors $A^{0,0}$ and $B^{0,1}$ [47–49].

Finally, a synthesis of the numerical results is presented for the various geometries: spherical, cylindrical, spherical/cylindrical, as well as sets composed of three cylindrical inclusions, followed by a comparative analysis. The conclusion emphasizes the importance of second-order gradient theory, particularly for high-contrast composites with long inclusions, such as long fibers.
## 2. Homogenization process and multiscale approach

The multi-scale method is founded on the asymptotic expansion techniques introduced by Bensoussan et al. [50,51], Sanchez-Palencia [52,53], and Mindlin [32]. Within the framework of the double-scale method, Sanchez-Palencia [54,55] delineates the periodic homogenization process for elasticity problems, which enables the determination of the effective tensor $C^{\text{hom}}$ through asymptotic expansion. The following provides a brief overview of the fundamental principles underlying this method.

### 2.1. Homogenization process

Let $\Omega$ be a three-dimensional periodic domain with period $P$, comprising various heterogeneities, as depicted in Figure 1. In this framework, the parameter $\varepsilon$ is defined as $\varepsilon = \frac{dl}{L}$, where $dl$ denotes the diameter of the inclusion, and $L$ represents a characteristic length of the material. Within this domain, the macroscopic coordinates of a point are given by $x = (x_1, x_2, x_3)$. The corresponding coordinates at the microscopic scale, within the representative volume element (RVE), are expressed as $y = \frac{x}{\varepsilon}$.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The principle of the homogenization technique.</figcaption></figure>

Given the external forces acting on the representative volume element (RVE), the mechanical equilibrium equations are expressed as follows:

$$ egin{cases} 
\operatorname{div}(\boldsymbol{\sigma}^{\varepsilon}) + f = 0 & \text{in } \Omega, \\ 
\boldsymbol{\sigma}^{\varepsilon} = C^{\varepsilon}(y)\operatorname{grad}(u^{\varepsilon}), & \\ 
u_{\Gamma_D}^{\varepsilon} = 0, & 
\end{cases} 
\quad (1) $$

where $\boldsymbol{\sigma}^{\varepsilon}$ denotes the stress tensor, $C^{\varepsilon}(y)$ represents the material stiffness tensor, and $u^{\varepsilon}$ corresponds to the displacement field within the domain $\Omega$. The boundary condition $u_{\Gamma_D}^{\varepsilon} = 0$ enforces a prescribed displacement on the Dirichlet boundary $\Gamma_D$.

In the first gradient homogenization theory [1,54,56,57], the scaling parameter $\varepsilon$ is assumed to approach zero, leading to the derivation of the first-order effective tensor $A^{0,0}$, commonly denoted as $C^{\text{hom}}$. Consequently, the governing equations in (1) reduce to:

$$ egin{cases} 
\operatorname{div}(\boldsymbol{\sigma}^0) + f = 0 & \text{in } \Omega, \\ 
\boldsymbol{\sigma}^0 = C^{\text{hom}}\operatorname{grad}(u^0), & \\ 
u_{\Gamma_D}^0 = 0. & 
\end{cases} 
\quad (2) $$
The homogenized tensor $C^{\text{hom}}$ is well established in classical homogenization theory and follows directly from the first-order gradient framework. Therefore, in this study, we extend the analysis to the second gradient theory, focusing on the determination of the higher-order effective tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$.

## 2.2. Asymptotic expansion

The asymptotic expansion method [32,50,52,54] follows a systematic sequence of steps. First, the displacement field $u^\epsilon$ is determined, and its gradient is computed. Next, this gradient is utilized to evaluate the associated stress field $\boldsymbol{\sigma}^\epsilon$. Finally, the divergence of $\boldsymbol{\sigma}^\epsilon$ is computed to obtain the governing equilibrium equations.

The asymptotic expansion of $u^\epsilon$ is traditionally formulated following the works of Bensoussan et al. [50,51], Sanchez-Palencia [52], and Mindlin [32]. The displacement field $u^\epsilon$ is expressed as:

$$u^{\epsilon} = u^{0}(x, y) + \epsilon u^{1}(x, y) + \epsilon^{2} u^{2}(x, y) + \epsilon^{3} u^{3}(x, y) + \dots \quad (3)$$

To introduce scale separation, we define the microscopic coordinate $y$ as $y = \frac{x}{\epsilon}$, where $\epsilon$ is a small parameter. The material stiffness matrix $C(y)$ is assumed to be positive definite and periodic.

To account for both macroscopic and microscopic scales, we introduce the “gradient” and “divergence” operators, denoted respectively as “grad” and “div”. These operators are decomposed as follows:

$$\left\{ \begin{array}{l} \operatorname{grad}(\bullet) = e_x \bullet + \epsilon^{-1} e_y \bullet, \\ \operatorname{div}(\bullet) = \operatorname{div}_x(\bullet) + \epsilon^{-1} \operatorname{div}_y(\bullet). \end{array} \right. \qquad (4)$$

Thus, the gradient of the displacement field expands as:

$$\operatorname{grad}(u^{\epsilon}) = \epsilon^{-1}(e_y u^0) + \epsilon^0(e_x u^0 + e_y u^1) + \epsilon^1(e_x u^1 + e_y u^2) + \epsilon^2(e_x u^2) + \dots \quad (5)$$

Accordingly, the stress field $\boldsymbol{\sigma}^\epsilon$ is expressed as:

$$\boldsymbol{\sigma}^{\epsilon} = \epsilon^{-1}[\mathbb{C}(y)(e_y u^0)] + \epsilon^0[\mathbb{C}(y)(e_x u^0 + e_y u^1)] \\ + \epsilon^1[\mathbb{C}(y)(e_x u^1 + e_y u^2)] + \epsilon^2[\mathbb{C}(y)(e_x u^2 + e_y u^3)] + \dots \quad (6)$$

Next, by substituting the asymptotic expansion of $\boldsymbol{\sigma}^\epsilon$ into the equilibrium equation (1), the balance equation transforms into:

$$\operatorname{div}_x(\boldsymbol{\sigma}^{\epsilon}) + \epsilon^{-1} \operatorname{div}_y(\boldsymbol{\sigma}^{\epsilon}) + f = 0. \qquad (7)$$

By systematically considering the scaling order of the parameter $\epsilon$, we obtain:

$$\begin{aligned} & \epsilon^{-2}[\operatorname{div}_y(\mathbb{C}(y)e_y u^0)] \\ & + \epsilon^{-1}[\operatorname{div}_x(\mathbb{C}(y)e_y u^0) + \operatorname{div}_y(\mathbb{C}(y)(e_x u^0 + e_y u^1))] \\ & + \epsilon^0[\operatorname{div}_x(\mathbb{C}(y)(e_x u^0 + e_y u^1)) + \operatorname{div}_y(\mathbb{C}(y)(e_x u^1 + e_y u^2))] + f \\ & + \epsilon^1[\operatorname{div}_x(\mathbb{C}(y)(e_x u^1 + e_y u^2)) + \operatorname{div}_y(\mathbb{C}(y)(e_x u^2 + e_y u^3))] + O(\epsilon) = 0. \end{aligned} \quad (8)$$

Each color represents a distinct order of magnitude with respect to the small parameter $\epsilon$, which characterizes the scale of rapid variations in the displacement field $u^\epsilon$ [32,42,50,52,54]. The different powers of $\epsilon$ highlight the contributions of multiple scales to the overall governing equation.
### 3. Effective determination of the tensors $\mathbb{A}^{0,0}$, $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$

This section focuses on the effective determination of the tensors $\mathbb{A}^{0,0}$, $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$, transitioning from a theoretical framework to an algorithmic approach based on the analysis of equation (8). This analysis results in four independent subproblems that highlight the role of the tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$.

In the following, the function space $H_{\text{per}}^{1}(P)$ is the periodic Sobolev space. Additionally, $\mathbb{I}$ denotes the identity tensor.

#### 3.1. Problem at order $\varepsilon^{-2}$

By identifying terms of order $\varepsilon^{-2}$, equation (8) simplifies to:

$$ \mathrm{div}_y(\mathbb{C}(y)e_y u^0) = 0. \qquad (9) $$

The corresponding variational formulation requires finding $u^0 \in H_{\text{per}}^{1}(P)$ such that:

$$ \int_P \mathbb{C}(y)(e_y u^0)(e_y u^*) dP = 0, \quad \forall u^* \in H_{\text{per}}^{1}(P). \qquad (10) $$

Since $\mathbb{C}(y)$ is positive definite, the Lax-Milgram theorem ensures the existence of a unique solution $u^0 \in H_{\text{per}}^{1}(P)$. Given that $u^0$ is independent of $y$, it follows that $u^0$ is solely a function of $x$.

#### 3.2. Problem at order $\varepsilon^{-1}$

Identifying terms of order $\varepsilon^{-1}$ in equation (8) leads to:

$$ \mathrm{div}_y(\mathbb{C}(y)(e_x u^0 + e_y u^1)) = 0. \qquad (11) $$

The corresponding variational formulation requires finding $u^1 \in H_{\text{per}}^{1}(P)$ such that:

$$ \int_P \mathbb{C}(y)(e_x u^0 + e_y u^1)(e_y u^*) dP = 0, \quad \forall u^* \in H_{\text{per}}^{1}(P). \qquad (12) $$

Since $\mathbb{C}(y)$ is positive definite, equation (12) must hold for all $u^* \in H_{\text{per}}^{1}(P)$, implying that $e_x u^0 + e_y u^1 = 0$.

Owing to the separation of scales, the displacement field $u^1(x, \frac{x}{\varepsilon})$ can be expressed as:

$$ u^1\left(x, \frac{x}{\varepsilon}\right) = U^1(x) + \chi^0(y) : e_x(u^0(x)), \qquad (13) $$

where $U^1(x)$ represents an integration constant.

The function $\chi^0(y)$ is the unknown corrector field, which satisfies the equation:

$$ \mathrm{div}_y(\mathbb{C}(y)(\mathbb{I} + e_y(\chi^0))) = 0. \qquad (14) $$

#### 3.3. Problem at order $\varepsilon^0$

By identifying terms of order $\varepsilon^0$, equation (8) simplifies to:

$$ \mathrm{div}_x(\mathbb{C}(y)(e_x u^0 + e_y u^1)) + \mathrm{div}_y(\mathbb{C}(y)(e_x u^1 + e_y u^2)) + f = 0. \qquad (15) $$

In the following, we assume the volumetric forces are negligible.

The variational formulation of this problem requires finding $u^2 \in H_{\text{per}}^{1}(P)$ such that:

$$ \int_P \mathbb{C}(y)(e_x u^0 + e_y u^1) \cdot (e_x u^*) dP + \int_P \mathbb{C}(y)(e_x u^1 + e_y u^2) \cdot (e_y u^*) dP = 0, \quad \forall u^* \in H_{\text{per}}^{1}(P). \qquad (16) $$
By exploiting periodicity, the equation simplifies to:

$$ \int_P \mathcal{C}(y)(e_x u^0 + e_y u^1) \cdot (e_x u^*) \, dP = 0, \quad \forall u^* \in H_{\text{per}}^1(P). \qquad (17) $$

Since $\mathcal{C}(y)$ is positive definite, equation (17) must hold for all $u^* \in H_{\text{per}}^1(P)$, ensuring the consistency of the formulation.

Within the framework of first-order gradient theory, the problem at order $\epsilon^0$ leads to the determination of the effective homogenized elasticity tensor, denoted as $\mathcal{C}^{\text{hom}}$ or $\mathcal{A}^{0,0}$. Its explicit expression is given by:

$$ \mathcal{A}^{0,0} = \mathcal{C}^{\text{hom}} = \int_P (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : (\mathbb{I} + e_y(\chi^0(y))) \, dP. \qquad (18) $$

The homogenized elasticity tensor $\mathcal{C}^{\text{hom}}$, also denoted as $\mathcal{A}^{0,0}$ (see Appendix A for exponent notation), is well established in the literature [1,37,41,44,54]. Various computational techniques, such as Fast Fourier Transform (FFT) and the Finite Element Method (FEM), have been developed for its efficient numerical evaluation. The following section will focus on the detailed derivation of the tensor $\mathcal{A}^{0,0}$, along with its index notation.

### 3.3.1. Order of $\mathcal{A}^{0,0}$

To rigorously determine the order of the tensor $\mathcal{A}^{0,0}$, we analyze the tensor products and contractions involved in its definition:

$$ \mathcal{A}^{0,0} = \mathcal{C}^{\text{hom}} = \int_P (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : (\mathbb{I} + e_y(\chi^0(y))) \, dP. \qquad (19) $$

The order determination follows from the structural analysis of each term:

**Component orders:**

* $\mathbb{I}$ is the fourth-order identity tensor.

* $\chi^0(y)$ is a third-order corrector tensor.

* The gradient operation applied to $\chi^0(y)$, namely $e_y(\chi^0(y))$, increases its order by one, resulting in a fourth-order tensor.

* $\mathcal{C}(y)$ is the fourth-order elasticity (rigidity) tensor.

**Addition of $\mathbb{I}$ and $e_y(\chi^0(y))$:**

Since both tensors are of order 4, their sum remains a fourth-order tensor:

$$ (\mathbb{I} + e_y(\chi^0(y))) \quad (\text{order 4}). $$

**First contraction:**

The first contraction is given by:

$$ (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y). $$

Contracting two fourth-order tensors along two indices reduces the total order by $2 \times p$, where $p = 2$ for a standard two-index contraction. Thus, the resulting tensor remains of order:

$$ 4 + 4 - 2 \times 2 = 4. $$

**Second contraction:**

The second contraction follows the same principle, reducing the order in the same manner:

$$ 4 + 4 - 2 \times 2 = 4. $$

After performing all contractions, the integral expression for $\mathcal{A}^{0,0}$ results in a fourth-order tensor. Since integration is applied element-wise and does not alter the tensor order, we conclude that $\mathcal{A}^{0,0}$ is a fourth-order tensor.
The subsequent section introduces the notation of $\mathcal{A}^{0,0}$ using index representation for further clarity and mathematical precision.

### 3.3.2. Index notation of the tensor $\mathcal{A}^{0,0}$

In this section, we express the tensor $\mathcal{A}^{0,0}$ using its index notation to explicitly illustrate its order and structure:

$$ \mathcal{A}_{mnpq}^{0,0} = \mathcal{C}_{mnpq}^{\text{hom}} = \int_P \mathcal{C}_{ijkh}(y) (\delta_{im}\delta_{jn} + e_{yij}(\chi_{mn}^0(y))(\delta_{kp}\delta_{hq} + e_{ykh}(\chi_{pq}^0(y)))) \, dP. \quad (20) $$

By explicitly incorporating index notation, the integral expression for $\mathcal{A}_{mnpq}^{0,0}$ confirms that it is a fourth-order tensor. Consequently, the order of $\mathcal{A}_{mnpq}^{0,0}$ is determined to be 4, with physical units of $(N/m^2)$, consistent with elasticity tensors.

The Kronecker delta $\delta_{ij}$ is a fundamental mathematical function of two integer variables, defined as:

$$ \delta_{ij} = \begin{cases} 1, & \text{if } i = j, \\ 0, & \text{otherwise.} \end{cases} \quad (21) $$

The next step, which constitutes the novel contribution of this work, is to address the problem at order $\epsilon^1$. Within the framework of second-order gradient theory, this analysis leads to the determination of the tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$, which will be systematically derived and analyzed in the subsequent sections.

### 3.4. Problem at order $\epsilon^1$

By identifying terms at order $\epsilon^1$, equation (8) takes the form:

$$ \operatorname{div}_x(\mathbb{C}(y)(e_x u^1 + e_y u^2)) + \operatorname{div}_y(\mathbb{C}(y)(e_x u^2 + e_y u^3)) = 0. \quad (22) $$

The variational formulation of this problem leads to finding $u^2 \in H_{\text{per}}^1(P)$ such that:

$$ \int_P \mathbb{C}(y)(e_x u^1 + e_y u^2) : (e_x u^*) \, dP + \mathbb{C}(y)(e_x u^2 + e_y u^3) : (e_y u^*) \, dP = 0 \quad \forall u^* \in H_{\text{per}}^1(P). \quad (23) $$

Due to the periodicity conditions, this equation simplifies to:

$$ \int_P \mathbb{C}(y)(e_x u^1 + e_y u^2) : (e_x u^*) \, dP = 0 \quad \forall u^* \in H_{\text{per}}^1(P). \quad (24) $$

Since $\mathbb{C}(y)$ is a positive-definite tensor, equation (24) must hold for all $u^* \in H_{\text{per}}^1(P)$. This directly implies that:

$$ e_x u^1 + e_y u^2 = 0. $$

Due to the separability of the scaling, the displacement field $u^2(x, \frac{x}{\epsilon})$ can be expressed as:

$$ u^2\left(x, \frac{x}{\epsilon}\right) = U^2(x) + \chi^0(y) : E^1(x) + \chi^1(y) : e_x(e_x u^0(x)). \quad (25) $$

where:

* $U^2(x)$ is an integration constant;

* $E^1(x) = (e_x(U^1(x)))$ represents the macroscopic strain field;

* $\chi^0(y)$ and $\chi^1(y)$ are the so-called corrector functions, determined as solutions to:

$$ \operatorname{div}_x(\mathbb{C}(y)(e_x u^1 + e_y u^2)) = 0. $$
Substituting the expressions for $(e_x u^1 + e_y u^2)$ and $(e_x u^0 + e_y u^1)$ into equation (24), we obtain:

$$ 
\begin{aligned}
& \int_P e_x(u^0) : (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : (\mathbb{I} + e_y(\chi^0(y))) : E^1(x) \, dP \\
& + \int_P e_x(u^0) : (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : \chi^0(y) : e_x(e_x(u^0)) \, dP \\
& + \int_P e_x(u^0) : (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : e_y(\chi^1(y)) : e_x(e_x(u^0)) \, dP = 0.
\end{aligned} 
\quad (26) $$

Within the framework of second-order gradient theory, the problem at order $\epsilon^1$ requires the determination of three fundamental tensors: $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$. These tensors characterize the homogenized response at this order and are derived through rigorous mathematical formulation. Their explicit expressions are provided below (see Appendix A for details on the exponent notation):

$$ \mathbb{B}^{0,1} = \int_P (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : (\mathbb{I} + e_y(\chi^0(y))) \, dP \quad (27) $$

$$ \mathbb{C}^{0,0} = \int_P (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : \chi^0(y) \, dP \quad (28) $$

$$ \mathbb{D}^{0,0} = \int_P (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) : e_y(\chi^1(y)) \, dP \quad (29) $$

The multi-scale expansion framework employed in this study has facilitated the derivation of tensor expressions analogous to those documented in the existing literature [37,40,44], obtained via the energy method. This substantiates the validity of our methodological approach but also underscores its originality. The results obtained reinforce the effectiveness of multi-scale tensor analysis, corroborating previous findings while introducing novel insights into the field.

The subsequent section focuses on establishing the order of the tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$, along with their index notation.

### 3.4.1. Order of the tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$

If $E$ and $F$ are tensors of orders $m$ and $n$, respectively, then the contraction of order $p$ between $E$ and $F$ (denoted as $E : F$, understood as a double contraction), where the contraction is performed $p$ times, results in a tensor of order $m+n-2p$. For this operation to be well-defined, it is necessary that $m \ge p$ and $n \ge p$.

To determine the order of $\mathbb{B}^{0,1}$, we analyze the tensor products and contractions in the equation (27):

**Component orders:**

* $\mathbb{I}$ is the identity tensor of order 4.

* $\chi^0(y)$ is a corrector tensor of order 3.

* The gradient of a rank-3 tensor, $e_y(\chi^0(y))$, increases the order by 1, making $e_y(\chi^0(y))$ a tensor of order 4.

* $\mathcal{C}(y)$ is the fourth-order rigidity tensor.

**Addition of $\mathbb{I}$ and $e_y(\chi^0(y))$:**

Since the addition of tensors does not change the highest order present, $\mathbb{I} + e_y(\chi^0(y))$ remains a tensor of order 4.

**First contraction operation:**

$$ (\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y) $$

Contracting a tensor of order 4 $(\mathbb{I} + e_y(\chi^0(y)))$ with another tensor of order 4 $\mathcal{C}(y)$ reduces the total order according to the contraction rule: $4+4-2 \times 2 = 4$. Thus, $(\mathbb{I} + e_y(\chi^0(y))) : \mathcal{C}(y)$ remains a tensor of order 4.
**Second contraction operation:**

Contracting the resulting tensor of order 4 with another ($\mathbb{I} + e_y(\chi^0(y))$) (also of order 4) reduces the order again: $4 + 4 - 2 \times 2 = 4$.

Considering the successive contractions, the integral expression for $\mathbb{B}^{0,1}$ results in a tensor of order 4 before integration. Since integration is an element-wise operation, it does not alter the tensor order. Consequently, $\mathbb{B}^{0,1}$ remains a tensor of order 4.

To determine the order of $\mathbb{C}^{0,0}$, we analyze the tensor products and contractions in the equation (28):

**Component orders:**

* $\mathbb{I}$ is the identity tensor of order 4.

* $\chi^0(y)$ is a corrector tensor of order 3.

* The gradient of a rank-3 tensor, $e_y(\chi^0(y))$, increases the order by 1, making $e_y(\chi^0(y))$ a tensor of order 4.

* $\mathbb{C}(y)$ is the rigidity tensor of order 4.

**Addition of $\mathbb{I}$ and $e_y(\chi^0(y))$:**

Since addition does not change the highest order term, $\mathbb{I} + e_y(\chi^0(y))$ remains a tensor of order 4.

**First contraction operation:**

$$ (\mathbb{I} + e_y(\chi^0(y))) : \mathbb{C}(y) $$

Contracting a tensor of order 4 ($\mathbb{I} + e_y(\chi^0(y))$) with another tensor of order 4 $\mathbb{C}(y)$ reduces the total order: $4 + 4 - 2 \times 2 = 4$. Thus, $(\mathbb{I} + e_y(\chi^0(y))) : \mathbb{C}(y)$ remains a tensor of order 4.

**Second contraction operation:**

Contracting the resulting tensor of order 4 with $\chi^0(y)$ of order 3 reduces the order further: $4 + 3 - 2 \times 2 = 3$.

After successive contractions, the integral expression for $\mathbb{C}^{0,0}$ results in a tensor of order 3 before integration. Since integration is an element-wise operation, it does not alter the tensor order. Consequently, $\mathbb{C}^{0,0}$ remains a tensor of order 3.

To determine the order of $\mathbb{D}^{0,0}$, we analyze the tensor products and contractions in the equation (29):

**Component orders:**

* $\mathbb{I}$ is the identity tensor of order 4.

* $\chi^0(y)$ is a corrector tensor of order 3.

* The gradient of a rank-3 tensor, $e_y(\chi^0(y))$, increases the order by 1, making $e_y(\chi^0(y))$ a tensor of order 4.

* $\chi^1(y)$ is identified as a tensor of order 4.

* The gradient of a rank-4 tensor, $e_y(\chi^1(y))$, increases the order by 1, making $e_y(\chi^1(y))$ a tensor of order 5.

* $\mathbb{C}(y)$ is the rigidity tensor of order 4.

**Addition of $\mathbb{I}$ and $e_y(\chi^0(y))$:**

Since addition does not change the highest order term, $\mathbb{I} + e_y(\chi^0(y))$ remains a tensor of order 4.

**First contraction operation:**

$$ (\mathbb{I} + e_y(\chi^0(y))) : \mathbb{C}(y) $$

Contracting a tensor of order 4 ($\mathbb{I} + e_y(\chi^0(y))$) with another tensor of order 4 $\mathbb{C}(y)$ reduces the total order: $4 + 4 - 2 \times 2 = 4$. Thus, $(\mathbb{I} + e_y(\chi^0(y))) : \mathbb{C}(y)$ remains a tensor of order 4.

**Second contraction operation:**

Contracting the resulting tensor of order 4 with $e_y(\chi^1(y))$ of order 5 reduces the order further: $4 + 5 - 2 \times 2 = 5$.
After successive contractions, the integral expression for $D^{0,0}$ results in a tensor of order 5 before integration. Since integration is an element-wise operation, it does not alter the tensor order. Consequently, $D^{0,0}$ remains a tensor of order 5.

The following section addresses the index notation for the tensors $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$.

### 3.4.2. Notation of $B^{0,1}$, $C^{0,0}$ and $D^{0,0}$ tensors with indices

In this section, the tensor $B^{0,1}$ is represented with indices to explicitly illustrate its order, using the Kronecker delta (21).

$$B_{mnpq}^{0,1} = \int_P C_{ijkh}(y)(\delta_{im}\delta_{jn} + e_{yij}(\chi_{mn}^0(y)))(\delta_{kp}\delta_{hq} + e_{ykh}(\chi_{pq}^0(y))) dP. \quad (30)$$

Considering the index notation, the integral expression of $B_{mnpq}^{0,1}$ results in a tensor of order 4. Thus, the order of the tensor $B_{mnpq}^{0,1}$ is determined to be 4, with the units (N/m²).

In the next section, we represent the tensor $C^{0,0}$ with its indices to clearly define its order.

$$C_{mnp}^{0,0} = \int_P C_{ijkh}(y)(\delta_{im}\delta_{jn} + e_{yij}(\chi_{mn}^0(y)))(\chi_p^{0kh}(y)) dP. \quad (31)$$

Taking into account the index notation, the integral expression for $C_{mnp}^{0,0}$ results in a tensor of order 3. Thus, the order of the tensor $C_{mnp}^{0,0}$ is determined to be 3, with the units (N/m).

In the existing literature, the $C^{0,0}$ tensor is typically expressed as a 5th-order tensor ($C_{(ij)mnp}^{0,0}$), which corresponds to augmenting the 3rd-order tensor ($C_{mnp}^{0,0}$) by adding two indices (ij) that represent stresses in different directions.

In the following section, the tensor $D^{0,0}$ is represented with its indices to clearly illustrate its order, using the Kronecker delta (21).

$$D_{mnpqr}^{0,0} = \int_P C_{ijkh}(y)(\delta_{im}\delta_{jn} + e_{yij}(\chi_{mn}^0(y))(e_{ykh}(\chi_p^{1qr}(y))) dP. \quad (32)$$

Considering the index notation, the integral expression for $D_{mnpqr}^{0,0}$ results in a tensor of order 5. Thus, the order of the tensor $D_{mnpqr}^{0,0}$ is determined to be 5, with the units (N/m).

The next section of this article delves into the demonstration of the symmetry properties of the tensors $A^{0,0}$, $B^{0,1}$, and $C^{0,0}$.

## 4. Symmetry properties of tensors $A^{0,0}$, $B^{0,1}$, and $C^{0,0}$

The work of Mindlin [32,57] and Toupin [30,31] demonstrated that the third-order tensor exhibits symmetry with respect to its first two indices. Consequently, the tensor $C_{ijk}^{0,0}$ is symmetric with respect to the indices i and j. The physical interpretation of this tensor was provided by Mindlin [32,57] and Polizzotto [58], who described it as a double stress tensor, denoted $C_{ijk}^{0,0}$, representing a double force oriented in the $x_k$ direction, applied to a plane with normal $x_i$, and possessing a lever arm in the $x_j$ direction.

This interpretation leads to an alternative rule—complementary to the aforementioned symmetry—whereby the first two indices of the tensor $C_{ijk}^{0,0}$ may be interchanged. In other words, the second index j denotes the direction of the lever arm, while the first index i corresponds to the normal of the loaded plane. This principle, commonly referred to as “Mindlin’s rule”, holds within the framework of gradient distortion elasticity (GDE) theory, wherein the dual stresses exhibit symmetry with respect to the first two indices—precisely the case in our study.

Mindlin’s findings [32,57] were subsequently extended by Exadakylos and Vardoulakis [59], as well as Lazar [60]. Consequently, the representation of the tensor $C_{ijk}^{0,0}$ is provided in Appendix D, further justifying its symmetry. In the following section, we explore the symmetry group configurations of the tensors $A^{0,0}$, $B^{0,1}$, and $C^{0,0}$ through a generalized approach.
Group symmetry theory [47–49] has established that an isotropic tensor $A^{0,0}$ or $B^{0,1}$ is characterized by two independent coefficients, with the third coefficient being uniquely determined by the first two. Similarly, the tensor $C^{0,0}$ is also defined by two independent coefficients.

Furthermore, group symmetry theory [47–49] has established that a tetragonal tensor $A^{0,0}$ or $B^{0,1}$ is characterized by six independent coefficients. Similarly, the tensor $C^{0,0}$ is defined by four independent coefficients.

Finally, group symmetry theory [47–49] has established that an orthotropic tensor $A^{0,0}$ or $B^{0,1}$ is characterized by nine independent coefficients. In contrast, the tensor $C^{0,0}$ is defined by six independent coefficients.

## 5. Connection and interpretation of higher-order stiffness tensors obtained with existing literature in gradient elasticity frameworks theory

To clarify the connection between the elastic energy expressed in terms of strain and its spatial derivatives, and the emergence of higher-order stiffness tensors, we begin by recalling the asymptotic expansion of the total elastic energy $W^\epsilon$ in powers of the characteristic microstructural length scale $\epsilon$:

$$W^{\epsilon} = W^{0} + \epsilon W^{1} + O(\epsilon^{2}). \quad (33)$$

Following the asymptotic homogenization approach detailed in [44], the leading-order energy $W^0$ takes the standard form of classical elasticity:

$$W^0 = \frac{1}{2} \int_P e_x(u^0) : A^{(0,0)} : e_x(u^0) dP, \quad (34)$$

where $e_x(u^0)$ is the symmetric strain tensor associated with the macroscopic displacement field $u^0$, and $A^{(0,0)}$ is the classical fourth-order stiffness tensor.

The first-order correction $W^1$ incorporates coupling between strain and its gradient, and is expressed as:

$$W^1 = \int_P e_x(u^0) : \mathbb{B}^{(0,1)} : e_x(u^0) dP + \int_P e_x(u^0) : \mathbb{C}^{(0,0)} : e_x(e_x(u^0)) dP + \int_P e_x(u^0) : \mathbb{D}^{(0,0)} : e_x(e_x(u^0)) dP, \quad (35)$$

where $\mathbb{B}^{(0,1)}$ is a fourth-order tensor, $\mathbb{C}^{(0,0)}$ and $\mathbb{D}^{(0,0)}$ are fifth-order tensors to capture higher-order effects.

Combining the expressions for $W^0$ and $W^1$, we obtain:

$$W^{\epsilon} = \frac{1}{2} \int_P e_x(u^0) (A^{(0,0)} + 2\epsilon B^{(0,1)}) e_x(u^0) dP + \epsilon \int_P e_x(u^0) (C^{(0,0)} + D^{(0,0)}) e_x(e_x(u^0)) dP + O(\epsilon^2). \quad (36)$$

This structure aligns with the second-order gradient elasticity models originally formulated by Toupin [30,31], Mindlin [32], and Green and Rivlin [28], as well as with the modern reformulation in [61], where the elastic energy is written as:

$$W^{\epsilon} = \frac{1}{2} e_x(u^0) : \mathbb{L}^{0,0} : e_x(u^0) + \epsilon e_x(u^0) : \mathbb{L}^{0,1} : e_x(e_x(u^0)) + \epsilon^2 e_x(e_x(u^0)) : \mathbb{L}^{1,1} : e_x(e_x(u^0)) + O(\epsilon^2), \quad (37)$$

where the tensorial moduli $\mathbb{L}^{0,0}$, $\mathbb{L}^{0,1}$, and $\mathbb{L}^{1,1}$ generalize the material response to include gradient effects.
By identifying terms order-by-order in $\epsilon$, the following equivalences are established:

$$ \mathbb{L}^{0,0} = \mathcal{A}^{(0,0)} + 2\epsilon \mathbb{B}^{(0,1)}, \quad \mathbb{L}^{0,1} = \mathcal{C}^{(0,0)} + \mathbb{D}^{(0,0)}. $$

Thus, the higher-order tensors $\mathcal{C}^{(0,0)}$ and $\mathbb{D}^{(0,0)}$ contribute additively to the effective second-gradient tensor $\mathbb{L}^{0,1}$, which plays the same role as the so-called $\mathbb{S}$ tensor in [1, (22)], provided that the minor symmetry $C_{ijk\ell m} = C_{ikj\ell m}$ is enforced. This mapping is summarized in Table 1, which provides a one-to-one correspondence between the present notation and that found in the literature.

<table><caption>**Table 1.** Correspondence of tensor notations and their physical roles.</caption><thead><tr><th>Tensor notation</th><th>Tensor order and role</th></tr></thead><tbody><tr><td>A<sup>(0,0)</sup></td><td>Fourth-order classical stiffness tensor</td></tr><tr><td>B<sup>(0,1)</sup></td><td>Fourth-order perturbation tensor (contribution to L<sup>0,0</sup>)</td></tr><tr><td>C<sup>(0,0)</sup></td><td>Fifth-order tensor</td></tr><tr><td>D<sup>(0,0)</sup></td><td>Fifth-order gradient tensor</td></tr><tr><td>L<sup>0,0</sup></td><td>Fourth-order stiffness tensor (sum of A<sup>(0,0)</sup> and 2εB<sup>(0,1)</sup>)</td></tr><tr><td>L<sup>0,1</sup></td><td>Effective “second-gradient” coupling tensor (sum of C<sup>(0,0)</sup> and D<sup>(0,0)</sup>)</td></tr></tbody></table>

# 6. Numerical approach and simulations

## 6.1. *Numerical approach*

This section presents results from two in-house computational codes. The first, based on the finite element method (FEM) [62], analyzes the mesh and computes the tensors $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$. Despite its computational cost, this approach serves as a benchmark for validation. The second code, based on the fast Fourier transform (FFT) method, incorporates a modified Green's kernel to handle high-contrast cases [45,46,63]. This approach significantly enhances computational efficiency, enabling faster evaluation of the tensors $A^{0,0}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$. Further details of the algorithms are provided in the subsequent section (see Appendix C).

In the next section, we present a graph depicting the logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and their respective values when the contrast is set to 1. In this graph, the logarithm of the contrast is plotted along the x-axis, while the logarithm of the coefficient ratio is represented on the y-axis.

Additionally, we provide a graph illustrating the logarithm of the absolute value of the ratio between the maximum coefficient of the tensor $D^{0,0}$ and its maximum value at a contrast of 1. Here, the x-axis represents the logarithm of the contrast, while the y-axis displays the logarithm of the absolute value of the coefficient ratio.

Furthermore, it is important to note that contrast is defined as the ratio between the properties of the inclusion and those of the matrix. All graphical representations are generated for contrast values ranging from 0.001 to 1000.

Finally, the occurrence of negative values in logarithmic calculations arises from the fact that the initial values are very small.

### 6.1.1. Validation FEM vs FFT: spherical inclusion in resolution of 128³

In this section, we present the results by comparing the coefficient $C_{111}$ of the $C^{0,0}$ tensor obtained using both the finite element method (FEM) and the fast Fourier transform (FFT) method for a volume fraction of $f_{sp} = 25\%$. The analysis considers a spherical inclusion inside
a representative volume element (RVE) of unit cubic shape, defined over the domain [0, 1]³. The center of the spherical inclusion is located at coordinates ($p_x = 0.5, p_y = 0.5, p_z = 0.5$) in three-dimensional space, with a radius of $r = 0.39$.

A resolution of $128^3$ was employed exclusively for this section with contrast values ranging from 0.001 to 1000, where contrast is defined as the ratio of the Young's modulus and Poisson's coefficient of the inclusion's properties to those of the matrix. The finite element method (FEM), employing a discretization technique, was also used to compute the results. Initially, a pixelated representation at a resolution of $128^3$ was generated. These pixels were then used to construct a tetrahedral mesh, which facilitated the computation of the tensors. For the representation of the tensor $C^{0,0}$, the coefficient $C_{111}$ was selected, noting that in this specific case, $C_{111} = C_{222} = C_{333}$.

The results obtained using the finite element method (FEM) were compared with those derived from the fast Fourier transform (FFT) method, and the relative error was computed, using the FEM curve as the reference to quantify discrepancies between the two approaches. Additionally, computation times were recorded to assess the performance of both methods.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast (FEM vs FFT) for a single spherical inclusion in a representative volume element (RVE) at a resolution of $128^3$.</figcaption></figure>

<table><caption>**Table 2.** Simulation time FEM vs FFT with a single spherical inclusion, a resolution of $128^3$.</caption><thead><tr><td>Contrast</td><td>Time with FEM (s)</td><td>Time with FFT (s)</td></tr></thead><tbody><tr><td>0.001</td><td>80350</td><td>450</td></tr><tr><td>0.01</td><td>59822</td><td>359</td></tr><tr><td>0.1</td><td>57030</td><td>169</td></tr><tr><td>1</td><td>53040</td><td>141</td></tr><tr><td>10</td><td>56199</td><td>159</td></tr><tr><td>32</td><td>58300</td><td>200</td></tr><tr><td>100</td><td>60725</td><td>400</td></tr><tr><td>180</td><td>64820</td><td>405</td></tr><tr><td>320</td><td>80645</td><td>430</td></tr><tr><td>570</td><td>82340</td><td>470</td></tr><tr><td>1000</td><td>104747</td><td>500</td></tr></tbody></table>
To summarize, it has been observed that the finite element method (FEM) generally produces more accurate results compared to those obtained using the fast Fourier transform (FFT) for the $A^{0,0}$ and $B^{0,1}$ tensors. For the $C^{0,0}$ tensor (see Figure 2a), the FEM results yield a curve that outperforms the one obtained with the FFT. However, the relative error between the two methods, as shown in Figure 2b, remains below 3%, which is generally deemed acceptable. Despite this, the FEM approach is computationally expensive (see Table 2), which has led us to prioritize the FFT method for ongoing work, while intermittently using FEM for validation purposes.

### 6.1.2. Validation FEM vs FFT: cylindrical inclusion in resolution of 128³

In this section, we present the results of comparing specific coefficients of the $C^{0,0}$ tensor in the case of a cylindrical inclusion with a volume fraction of $f_{cy} = 5.71\%$, embedded within a cubic-shaped representative volume element (RVE) with dimensions $[0, 1]^3$. The center of the cylinder is positioned at the coordinates $(p_x = 0.5, p_y = 0.5, p_z = 0.5)$ in three-dimensional space. The cylinder has a radius of $r = 0.135$, and its axis, aligned along the y-axis, is characterized by the directional vector $(l_x = 0, l_y = 1, l_z = 0)$. The height of the cylinder is $h = 1$.

A resolution of 128³ is also employed exclusively in this section, with contrast values ranging from 0.001 to 1000. Calculations have been carried out in the same way as in the previous case.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast (FEM vs FFT) for a single cylindrical inclusion in a representative volume element (RVE) at a resolution of 128³.</figcaption></figure>
In the representation of the tensor $C^{0,0}$, the coefficients $C_{111}$ and $C_{222}$ are selected, with $C_{111} = C_{333}$. However, $C_{222}$ differs due to the directional orientation of the inclusion, which is aligned along the y-axis. The finite element method (FEM) results are compared with those obtained using the fast Fourier transform (FFT), and the relative error is calculated, with the FEM curve serving as the reference for assessing the differences between the two methods. Additionally, computation times are recorded to evaluate the performance of both approaches.

Negative values in logarithmic calculations arise from initially small values.

In conclusion, it is observed that the finite element method (FEM) tends to yield superior results compared to those obtained from the fast Fourier transform (FFT) for the $A^{0,0}$ and $B^{0,1}$ tensors. For the tensor $C^{0,0}$ (see Figures 3a and 3c), the results from the finite element method (FEM) exhibit a curve that is superior to that obtained using the fast Fourier transform (FFT). However, the relative error between the two methods, as shown in Figures 3b and 3d, remains below 7% for a resolution of $128^3$. While this error may be considered relatively high, it is deemed acceptable, particularly since it decreases with increasing resolution. Consequently, we plan to continue our work using a resolution of $256^3$.

Despite its advantages, the finite element method (FEM) is computationally intensive (see Table 3), prompting us to prioritize the fast Fourier transform (FFT) for further analyses at the $256^3$ resolution, while intermittently employing the finite element method (FEM) to validate our results.

In the following section, numerical simulations are performed using our in-house developed code for spherical and cylindrical inclusions. These simulations enable the determination of the homogenized stiffness tensor ($C^{hom}$ or $A^{0,0}$) based on first-order gradient theory, as well as the tensors $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$ derived from second-order gradient theory. Several distinct cases are examined to illustrate the numerical values of these tensors.

<table><caption>**Table 3.** Simulation time FEM vs FFT with a single cylindrical inclusion, resolution of $128^3$.</caption><thead><tr><td>Contrast</td><td>Time with FEM (s)</td><td>Time with FFT (s)</td></tr></thead><tbody><tr><td>0.001</td><td>51240</td><td>401</td></tr><tr><td>0.01</td><td>32020</td><td>350</td></tr><tr><td>0.1</td><td>29844</td><td>201</td></tr><tr><td>1</td><td>28120</td><td>130</td></tr><tr><td>10</td><td>29932</td><td>141</td></tr><tr><td>32</td><td>31825</td><td>191</td></tr><tr><td>100</td><td>34628</td><td>307</td></tr><tr><td>180</td><td>36514</td><td>321</td></tr><tr><td>320</td><td>44710</td><td>349</td></tr><tr><td>570</td><td>48840</td><td>379</td></tr><tr><td>1000</td><td>53676</td><td>400</td></tr></tbody></table>

## 6.2. Numerical simulations

Several simulations were conducted at a resolution of $256^3$ voxels, on a representative volume element (RVE) composed of inclusion(s) within a unit cubic domain $[0, 1]^3$, with contrast levels ranging from 0.001 to 1000, where contrast is defined as the ratio between the Young's modulus and Poisson's ratio of the inclusion and those of the matrix. These simulations were categorized into two groups: the first investigated inclusions embedded within a representative volume
element (RVE), while the second focused on elongated inclusions in contact with the edges of the RVE. This approach enables a comprehensive analysis of the effects of different inclusion geometries.

Initially, our study focused on a spherical inclusion, which inherently exhibits no preferred orientation. Subsequent simulations examined short and long cylindrical inclusions, analogous to short and long fibers, respectively. Additionally, a combination of these various types of inclusions was analyzed to assess their overall impact, culminating in a study of cylindrical inclusions oriented along the three principal axes to investigate potential compensatory effects.

A volume fraction of 5.71% was deliberately chosen to facilitate the combination of different cases. These simulations enabled the determination of the homogenized tensor $A^{0,0}$ or $C^{\text{hom}}$ using first-order gradient theory, as well as the tensors $\mathbb{B}^{0,1}$, $C^{0,0}$, and $\mathbb{D}^{0,0}$ using second-order gradient theory.

It is important to emphasize that, in all configurations analyzed, the tensor $A^{0,0}$ is identical to the tensor $\mathbb{B}^{0,1}$ with remarkable precision, exhibiting an accuracy of $10^{-5}$ beyond the decimal point. This demonstrates an exceptionally precise correspondence between the two tensors. Although the computational approaches for $A^{0,0}$ and $\mathbb{B}^{0,1}$ differ, both yield identical results, thereby confirming the robustness and reliability of our numerical framework.

In this study, the tensor $C^{0,0}$ is represented in matrix form as two distinct blocks (see Appendix D): an upper block and a lower zero block. This representation facilitates the integration of $C^{0,0}$ into the symmetry group of the tensors $A^{0,0}$ and $\mathbb{B}^{0,1}$.

Additionally, we analyze the logarithm of the ratio between the coefficients of the $C^{0,0}$ tensor and the corresponding coefficient value when the contrast is set to 1. In this representation, the logarithm of the contrast is plotted along the x-axis, while the logarithm of the coefficient ratio is displayed on the y-axis.

Furthermore, the tensor $\mathbb{D}_{mnpqr}^{0,0}$ has been represented (see Appendix E) as a curve by plotting the logarithm of the absolute value of the ratio between the maximum coefficient of $\mathbb{D}_{mnpqr}^{0,0}$ and its maximum value at a contrast of 1, as a function of the logarithm of the contrast. This approach enables the observation of the tensor's variation with respect to contrast. In other words, the $L_\infty$ norm was chosen to characterize the $\mathbb{D}_{mnpqr}^{0,0}$ tensor, defined as: $||\mathbb{D}_{mnpqr}^{0,0}||_\infty = \max_{m,n,p,q,r} |\mathbb{D}_{mnpqr}^{0,0}|$.

The next phase of our study explores an alternative geometric configuration, focusing on inclusions embedded within a representative volume element (RVE). This approach facilitates a detailed analysis of the influence of inclusion shape on the properties of the investigated tensors. The primary objective is to gain a deeper understanding of the impact of geometry on the overall mechanical behavior of the material.

### 6.2.1. Inclusions immersed in a representative volume element (RVE)

In this section, we have conducted simulations incorporating various types of inclusions within a representative volume element (RVE). The studied inclusion shapes include spheres, and cylinders, as well as combinations of spheres with cylinders. The analysis concludes with an investigation of cylindrical inclusions oriented along three different directions.

The next phase of our study focuses on a distinct geometric configuration, specifically examining a spherical inclusion. This approach enables a detailed assessment of the influence of inclusion shape on the properties of the studied tensors. The primary objective is to gain a deeper understanding of the impact of geometry on the overall behavior of the material.

**Spherical inclusion.** Two volume fractions were considered: $f_{sp} = 5.71\%$ with a sphere radius of $r = 0.239$, and $f_{sp} = 25\%$ with a radius of $r = 0.39$. The sphere is centrally positioned at coordinates $(p_x = 0.5, p_y = 0.5, p_z = 0.5)$ within the three-dimensional space.

The two graphs presented in Figure 5 illustrate the variation of the tensor $C^{0,0}$ as a function of the logarithm of contrast for different components. The analysis of the $C^{0,0}$ tensor coefficients
<figure><img src="image_7.png" /></figure>

Case of a single sphere:

• Contrast: 0.001 up to 1000.

• Resolution: 256³.

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Illustration of a single sphere inside a 3D domain.</figcaption></figure>

reveals intriguing behaviors, highlighting the significance of isotropic symmetry in spherical inclusion materials. It is evident that the coefficients of the $C^{0,0}$ tensor increase as the volume fraction increases, as shown in Figure 5.

For the coefficients $C_{111} = C_{222} = C_{333}$ (see Figure 5a), we observe that they vary from -0.2 to a maximum of 0.2 for $f_{sp} = 5.71\%$ and from -0.6 to a maximum of 0.6 for $f_{sp} = 25\%$ with the logarithm of contrast. However, they remain identical, which confirms the absence of a preferential direction, a characteristic feature of spherical inclusions and isotropic symmetry. These variations, observed across all contrast levels, confirm the material's isotropic nature and the absence of directional bias.

Similarly, the other coefficients $C_{112} = C_{113} = C_{221} = C_{223} = C_{331} = C_{332}$ (see Figure 5b) remain constant across varying contrast levels, illustrating the uniformity of the material's properties. These results are particularly significant for material design, as they demonstrate that the properties remain uniform despite variations in contrast—an essential aspect for applications requiring directional consistency. Furthermore, since the coefficients have magnitudes on the order of unity, this further substantiates the isotropic symmetry, as illustrated in Figure 6.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single sphere.</figcaption></figure>
Figure 6 presents the configuration of the tensor $C^{0,0}$, where elements of the same color (red and blue) correspond to identical values. The red squares, positioned along the diagonal of the upper block, have equal values due to the absence of a preferred direction. The blue squares, which are also uniform, represent the other values in the upper block. At the bottom, the “0” symbols indicate null values or inactive elements.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** $C^{0,0}$ tensor configuration with a single sphere.</figcaption></figure>

Figure 7 shows that variations in contrast induce only a slight variation in $||D^{0,0}||_\infty$. Additionally, it is observed that the coefficients of the $D^{0,0}$ tensor increase as the volume fraction increases. This suggests that contrast does not significantly influence the behavior of this particular component of the tensor when dealing with a spherical inclusion.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single sphere.</figcaption></figure>

**Cylindrical inclusion with cylinder height of 0.5.** In this section, the cylindrical inclusion, with a volume fraction $f_{cy} = 5.71\%$, is centrally positioned at coordinates $(p_x = 0.5, p_y = 0.5, p_z = 0.5)$ in three-dimensional space. It has a radius of $r = 0.165$ and a height of $h = 0.5$. The cylinder's orientation along the y-axis is defined by the directional vector $(l_x = 0, l_y = 0.25, l_z = 0)$. Simulations are carried out at contrast levels ranging from 0.001 to 1000, where contrast is defined as the ratio of inclusion to matrix, with a resolution of $256^3$ voxels. This study focuses on the analysis of the tensors $A^{0,0}$ or $C^{hom}$, $B^{0,1}$, $C^{0,0}$, and $D^{0,0}$ under varying contrast.
<figure><img src="image_13.png" /></figure>

Case of a single cylinder following y:

• $f_{cy}$: 5.71%.

• Contrast: 0.001 up to 1000.

• Resolution: 256³.

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient at a contrast of 1, plotted as a function of the logarithm of contrast for a single cylindrical inclusion with a height of 0.5.</figcaption></figure>

We present three curves in Figure 9 that represent the coefficients of the $C^{0,0}$ tensor as a function of the logarithm of contrast for a material with a cylindrical inclusion of height 0.5. These graphs highlight the distinct behaviors of the various coefficients of the $C^{0,0}$ tensor under different contrast conditions.
The first curve (see Figure 9a) represents the behavior of the coefficients $C_{222} \neq (C_{111} = C_{333})$, showing a slight increase with increasing contrast. While $C_{222}$ increases and reaches values higher than $C_{111}$ and $C_{333}$, the other two coefficients show a more moderate increase. This pronounced difference underscores the influence of the inclusion's orientation along the y-axis, corresponding to the $C_{222}$ coefficient, suggesting strong directional dependence in the material under high-contrast conditions.

The second curve (see Figure 9b) highlights the behavior of the coefficients $C_{113} = C_{331}$, while the third curve (shown in Figure 9c) depicts the evolution of the coefficients $C_{112} = C_{221} = C_{223} = C_{332}$. These coefficients show a consistent trend, increasing with contrast, which reflects greater anisotropy in the material at higher contrast levels.

<figure><img src="image_15.png" /><figcaption>**Figure 10.** $C^{0,0}$ tensor configuration with a single cylinder with a cylinder height of 0.5.</figcaption></figure>

Figure 10 illustrates the $C^{0,0}$ tensor, obtained from the results shown in Figure 9. In this configuration, elements of the same color represent identical values, indicating symmetry within the tensor. The upper block uses colors such as red, yellow, blue, and green to denote equal values among corresponding tensor elements, while the yellow color represents a unique value due to the orientation of the inclusion along the y-axis. The lower block of the tensor consists entirely of zeros, signifying that all elements in this section are null. This configuration, with a non-zero upper block and a lower block with zeros, highlights the structured symmetry of the $C^{0,0}$ tensor.

In summary, Figure 9 establishes a clear link between the length of inclusions in a composite material and its mechanical properties, as illustrated by the coefficients of the $C^{0,0}$ tensor. This tensor, which is essential for describing the elastic properties of the material, shows a marked variation in its coefficients depending on the length of the inclusion. The elongation of the inclusion significantly affects the mechanical properties along its main axis, enabling the material to withstand increased stress or deformation. The coefficients of the $C^{0,0}$ tensor, expressed in newtons per meter (N/m), quantify the force or tension per unit length, which is typical for characterizing the elasticity of materials. Consequently, a proportional increase in the length of the inclusion leads to an increase in these coefficients, indicating an enhanced ability of the material to resist or transmit forces along the axis of the inclusion.

Finally, we turn our attention to the $D^{0,0}_{mnpqr}$ tensor. As illustrated in Figure 11, the curve reveals that, for a cylindrical inclusion with a height of 0.5, the variations in contrast significantly influence $||D^{0,0}||_\infty$. The tensor's behavior demonstrates a clear increasing trend with rising contrast, suggesting that contrast is a pivotal factor in determining $||D^{0,0}||_\infty$. Based on these findings, it is evident that the second-order gradient theory must be considered when dealing with high contrast in cylindrical inclusions.

The next phase of our study shifts focus to a different geometry, specifically investigating a mixture with the previous inclusions: single spherical and cylindrical inclusions. This approach
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single cylinder with a cylinder height of 0.5.</figcaption></figure>

facilitates the analysis of how the shape of the inclusion influences the properties of the studied tensors. The objective is to gain a deeper understanding of the impact of inclusion geometry on the material's behavior.

**Spherical and cylindrical inclusions with cylinder height of 0.5.** In this section, we perform a numerical simulation of a representative volume element (RVE) comprising both spherical and cylindrical inclusions ($f_{sp} = f_{cy} = 5.71\%$), with a total volume fraction of 11.42%. The sphere is centered at coordinates ($p_x = 0.25, p_y = 0.5, p_z = 0.5$) with a radius of $r = 0.239$. The cylindrical inclusion is positioned at ($p_x = 0.75, p_y = 0.5, p_z = 0.5$) and has a radius of $r = 0.165$ and a height of $h = 0.5$, oriented along the y-axis with a directional vector ($l_x = 0, l_y = 0.25, l_z = 0$).

<figure><img src="image_17.png" /><figcaption>**Figure 12.** Illustration of a single sphere and cylinder with a cylinder height of 0.5 inside a 3D domain.</figcaption></figure>

**Case of a single sphere and cylinder following y:**

* Total volume fraction: 11.42%.

→ With a $f_{sp}$ of 5.71%.

→ With a $f_{cy}$ of 5.71%, the same as the case of a single cylindrical inclusion with $h = 0.5$.

* Contrast: 0.001 up to 1000.

* Resolution: $256^3$.

<figure><img src="image_18.png" /></figure>

Figure 13a demonstrates the behavior of the coefficients $C_{111} \neq C_{222} \neq C_{333}$, showing that $C_{222}$ does not effectively capture the directional properties of the cylinder due to the equal height of both the cylinder and the sphere. This outcome aligns with expectations, as the side-by-side arrangement of the sphere and cylinder leads to a directional dominance along the x-axis. Furthermore, the coefficients vary between -0.5 and 0.5, indicating that they are of the order of unity, which underscores the relevance of the $C^{0,0}$ tensor in characterizing the material properties.
Figure 13b illustrates the variation of coefficients $C_{113} = C_{331}$, while Figure 13c demonstrates the variation of the coefficients $C_{112} = C_{221}$ and $C_{223} = C_{332}$. These coefficients fluctuate between -0.4 and 0.4, confirming that they are of the order of unity. This fluctuation highlights the significant role of the $C^{0,0}$ tensor in revealing how material responses become increasingly anisotropic under high contrast conditions.

<figure><img src="image_19.png" /><figcaption>**Figure 13.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast for a single sphere and cylinder, with the cylinder having a height of 0.5.</figcaption></figure>

Figure 14 presents the configuration of the $C^{0,0}$ tensor derived from the results of Figure 13, where elements of the same color correspond to identical values, indicating symmetry within the tensor. The upper block is color-coded, with red, yellow, black, blue, green, and cyan representing equal values among corresponding tensor elements. However, the red, yellow, and black colors represent unique values that differ from the others. The lower block of the tensor is entirely composed of zeros, indicating that all elements in this section are null. This configuration, with a non-zero upper block and a lower block of zeros, highlights the structured symmetry of the $C^{0,0}$ tensor.

In summary, as contrast increases, the coefficients of the $C^{0,0}$ tensor also increase. We observe that the direction of the short cylinder, which has the same size as the sphere, becomes indistinguishable, aligning with our expectations. Additionally, the side-by-side positioning of the sphere and cylinder contributes to the directional predominance along the x-axis. Moreover,
the cumulative effect of the tensor coefficients for the sphere and cylinder indicates synergistic interactions within the composite. In this configuration, the $C^{0,0}$ tensor exhibits orthotropic behavior. These results are pivotal for understanding the material's behavior under varying contrast conditions.

Figure 15 indicates that, in the case of single spherical and cylindrical inclusions with a cylinder height of 0.5, variations in contrast significantly affect $||D^{0,0}||_\infty$. The tensor behavior shows an increasing trend with increasing contrast, suggesting that contrast plays a crucial role in determining $||D^{0,0}||_\infty$.

<figure><img src="image_20.png" /><figcaption>**Figure 14.** $C^{0,0}$ tensor configuration with a single sphere and cylinder with a cylinder height of 0.5.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 15.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single sphere and cylinder with a cylinder height of 0.5.</figcaption></figure>

The next phase of our study focuses on a particular geometry, involving the analysis of three cylindrical inclusions with a height of 0.75.

**Three cylindrical inclusions with cylinder height of 0.75.** In this section, we perform a numerical simulation of a representative volume element (RVE) within a cubic space $[0, 1]^3$, containing three cylindrical inclusions, each contributing to a total volume fraction of 12.94%. The centers
of all three cylinders are positioned at the coordinates ($p_x = 0.5, p_y = 0.5, p_z = 0.5$). Each cylinder has a radius of $r = 0.156$ and a height of $h = 0.75$. The orientations of the cylinders are aligned along the x, y, and z axes, respectively, with directional vectors ($l_x = 0.375, l_y = 0, l_z = 0$) for the first, ($l_x = 0, l_y = 0.375, l_z = 0$) for the second, and ($l_x = 0, l_y = 0, l_z = 0.375$) for the third.

<figure><img src="image_22.png" /></figure>

**Case of three cylindrical inclusions positioned along the three axes, each with a volume fraction of 5.71%:**

* Total volume fraction: 12.94%.

* Contrast: 0.001 up to 1000.

* Resolution: 256³.

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Illustration of three cylinder with cylinder height of 0.75 inside a 3D domain.</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 17.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a three cylinder with cylinder height of 0.75.</figcaption></figure>

The two graphs presented in Figure 17 illustrate the variation of the tensor $C^{0,0}$. For the coefficients $C_{111} = C_{222} = C_{333}$ (see Figure 17a), we observe that they vary from -0.6 to a maximum of 0.6 with the logarithm of contrast, yet remain identical. This behavior indicates the absence of a preferential direction, which is characteristic of configurations with three cylindrical inclusions oriented along the x, y, and z axes, thereby demonstrating isotropic symmetry.

For the coefficients $C_{112} = C_{113} = C_{221} = C_{223} = C_{331} = C_{332}$, as shown in Figure 17b, we observe an increase with increasing contrast, ranging from -0.8 to 0.2. These coefficients, being of the order of unity, demonstrate a significant material response to variations in contrast.

In summary, the graphs presented in Figure 17 establish a clear relationship between the inclusion length in a composite material and its mechanical properties, as illustrated by the coefficients of the tensor $C^{0,0}$. This tensor, crucial for describing the elastic properties of the material,
shows marked variations in its coefficients as a function of inclusion length. Consequently, a proportional increase in the inclusion length leads to a rise in these coefficients, indicating the material's improved ability to resist or transmit forces along the inclusion's axis. Additionally, there is a cumulative effect in the coefficients of the $C^{0,0}$ tensor for each cylinder, with the overall response reflecting the combined influence of all three cylindrical inclusions.

Figure 18 displays the configuration of the tensor $C^{0,0}$, where elements of the same color (red and blue) are assigned identical values. The red squares, arranged along the diagonal of the upper block, exhibit identical values, reflecting the absence of a preferred directional influence. This uniformity results from the strategic placement of three cylindrical inclusions, each oriented differently: the first along the x-axis, the second along the y-axis, and the third along the z-axis. The blue squares, also identical to each other, represent the other values in the upper block. At the bottom, the “0” symbols indicate null values or inactive elements, representing areas unaffected by the inclusions.

Finally, in Figure 19, the curve indicates that, for the case of three cylindrical inclusions with a cylinder height of 0.75, variations in contrast affect $||D^{0,0}||_\infty$. The behavior of the tensor shows an upward trend with increasing contrast, highlighting the critical role of contrast in determining $||D^{0,0}||_\infty$.

<figure><img src="image_26.png" /><figcaption>**Figure 18.** $C^{0,0}$ tensor configuration with a three cylinder with cylinder height of 0.75.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 19.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with three cylinders with a cylinder height of 0.75.</figcaption></figure>
This analysis reveals interesting behaviors that underscore the importance of isotropic symmetry in materials with three cylindrical inclusions, each oriented differently: the first along the x-axis, the second along the y-axis, and the third along the z-axis. These variations highlight how orientation affects the tensor properties, contributing to a deeper understanding of isotropic responses in these materials.

**Summary and comparison of first numerical results.** The geometries of inclusions in composites show varied responses under contrast-related variations. Spherical inclusions retain isotropic properties, while cylindrical inclusions exhibit significant anisotropy. The coefficients of the $C^{0,0}$ tensor are close to unity and increase with contrast, even at low volume fractions. When combined, the anisotropy persists and the cumulative effect of the $C^{0,0}$ tensor coefficients becomes more pronounced, particularly in scenarios where the volume fraction of the inclusion is high. These results highlight the crucial role of inclusion geometry in the mechanical response of materials, particularly in high-contrast composites.

In conclusion, second-order gradient theory is essential to accurately model these composites, particularly when the inclusions are submerged, resulting in a significant volume fraction. A comparative analysis of the $C_{111}$ coefficient in the $C^{0,0}$ tensor is presented for detailed study, in Figure 20, highlighting the importance of this advanced modelling approach.

The results presented in Figure 20 demonstrate that an increase in the volume fraction leads to a corresponding rise in the coefficients of the $C^{0,0}$ tensor. This effect becomes more pronounced at higher contrast levels, highlighting the substantial influence of contrast on the material's behavior. Furthermore, when multiple inclusions are combined, a cumulative effect on the coefficients of the $C^{0,0}$ tensor is observed. These findings emphasize the critical importance of incorporating second-order gradient theory in the modeling of high-contrast composite materials.

Figure 20 establishes a clear correlation between the length of inclusions in a composite material and its mechanical properties, as reflected by the variations in the coefficients of the $C^{0,0}$ tensor. The elongation of the inclusions significantly influences the material's mechanical properties along the direction of extension, enhancing its capacity to endure greater stress or strain along the major axis.

<figure><img src="image_28.png" /><figcaption>**Figure 20.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of the contrast.</figcaption></figure>
The next phase of this study shifts focus to a new inclusion geometry, involving elongated inclusions that extend to the edges of the representative volume element (RVE). This configuration facilitates the exploration of how the shape of inclusions influences the properties of the tensors under investigation. The goal of this approach is to deepen our understanding of how geometric configurations affect the material's overall behavior.

### 6.2.2. Elongated inclusions touching the edges of the representative volume element (RVE)

In this section, we present simulations that incorporate various types of inclusions within a representative volume element (RVE). The inclusion shapes analyzed include spheres, cylinders, as well as combinations of spheres with cylinders. The study concludes with an investigation into cylindrical inclusions oriented along three different directions.

The next phase of this study shifts focus to a different geometry, specifically involving cylindrical inclusions. The goal of this approach is to explore the influence of inclusion shape on the properties of the tensors under investigation. Ultimately, this will deepen our understanding of how geometry impacts the overall behavior of the material.

**Cylindrical inclusion with cylinder height of 1.** In this section, we perform a numerical simulation of a representative volume element (RVE) comprising a cylindrical inclusion embedded within a cubic domain $[0, 1]^3$. The cylinder occupies a volume fraction of $f_{cy} = 5.71\%$, and is centrally located at the coordinates $(p_x = 0.5, p_y = 0.5, p_z = 0.5)$ within the three-dimensional space. The inclusion has a radius of $r = 0.135$ and is oriented along the y-axis, defined by the directional vector $(l_x = 0, l_y = 0.5, l_z = 0)$. Its height is $h = 1$.

<figure><img src="image_29.png" /><figcaption>**Figure 21.** Illustration of a single cylinder with a cylinder height of 1 inside a 3D domain.</figcaption></figure>

We present three curves in Figure 22 that depict the coefficients of the $C^{0,0}$ tensor as a function of the logarithm of contrast for a material containing a cylindrical inclusion.

The first curve (see Figure 22a) shows the behavior of the coefficients $C_{222}$ compared to $C_{111}$ and $C_{333}$, emphasizing how each coefficient responds to increasing contrast. While $C_{222}$ exhibits the steepest increase after a contrast value of 100, reaching significantly higher values than $C_{111}$ and $C_{333}$, the other two coefficients show more moderate increases. This pronounced difference highlights the influence of the inclusion's orientation along the y-axis (associated with $C_{222}$), suggesting that the material exhibits a strong directional dependence under high-contrast conditions. The rapid increase of $C_{222}$ implies a significant enhancement in stiffness in the y-direction.

The second curve (see Figure 22b) illustrates the variation of the coefficients $C_{113} = C_{331}$, while the third curve (see Figure 22c) depicts the behavior of the coefficients $C_{112} = C_{221} = C_{223} = C_{332}$. All of these coefficients display similar behavior: they remain nearly constant for contrast values below 100. However, beyond this threshold, a significant increase is observed, reflecting a greater anisotropy of the material at higher contrast values. This suggests that while the material behaves
isotropically at low contrast, it becomes progressively more anisotropic as the contrast increases. The $C^{0,0}$ tensor thus plays a crucial role under high-contrast conditions.

Figure 23 shows the configuration of the $C^{0,0}$ tensor derived from the results in Figure 22, where elements of the same color are associated with identical values, reflecting the symmetry within the tensor. The upper block is color-coded: red, yellow, blue, and orange regions correspond to equal values, illustrating the symmetry between the respective tensor elements. Notably, the yellow elements represent a distinct value, which differs from the others due to the orientation of the inclusion along the y-axis. The lower block of the tensor is entirely composed of zeros, signifying that all elements in this section are null. This configuration, with a non-zero upper block and a lower block filled with zeros, emphasizes the structured symmetry of the $C^{0,0}$ tensor.

<figure><img src="image_30.png" /><figcaption>**Figure 22.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single cylinder with a cylinder height of 1.</figcaption></figure>

In summary, the results shown in Figure 22 reveal that at low contrast, the material exhibits relatively uniform properties in all directions, as indicated by the nearly constant behavior of the tensor coefficients. However, as the contrast increases, a significant anisotropy emerges across all coefficients, indicating that the material becomes increasingly directional. This highlights the importance of the $C^{0,0}$ tensor under high contrast conditions.
In the subsequent section, we turn our attention to the $D_{mnpqr}^{0,0}$ tensor, a second-order tensor. In Figure 24, the curve indicates that, for a cylindrical inclusion, variations in contrast significantly affect $||D^{0,0}||_\infty$. The tensor behavior exhibits an upward trend with increasing contrast, suggesting that contrast plays a critical role in determining $||D^{0,0}||_\infty$. Based on these results, it is evident that second-order gradient theory must be considered for high-contrast scenarios involving cylindrical inclusions, as the values of the tensors become significantly large.

<figure><img src="image_31.png" /><figcaption>**Figure 23.** $C^{0,0}$ tensor configuration with a single cylinder with a cylinder height of 1.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 24.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single cylinder with a cylinder height of 1.</figcaption></figure>

**Spherical and cylindrical inclusions with cylinder height of 1.** In this section, we perform a numerical simulation of a representative volume element (RVE) with a unit cubic shape [0, 1]³, incorporating both spherical and cylindrical inclusions. The total volume fraction of the inclusions is $f_{sp} + f_{cy} = 5.71\% + 5.71\% = 11.42\%$. The spherical inclusion is positioned at coordinates $(p_x = 0.25, p_y = 0.5, p_z = 0.5)$ with a radius of $r = 0.239$. The cylindrical inclusion is located at $(p_x = 0.75, p_y = 0.5, p_z = 0.5)$ with a radius of $r = 0.135$. Its axis is oriented along the y-axis, defined by the directional vector $(l_x = 0, l_y = 0.5, l_z = 0)$, and it has a height of $h = 1$.

Three curves (see Figure 26) represent the coefficients of the $C^{0,0}$ tensor as a function of the logarithm of contrast for a material containing spherical and cylindrical inclusions.

Figure 26a illustrates the behavior of the coefficients $C_{111} \neq C_{222} \neq C_{333}$, where $C_{222}$ exhibits a particularly strong increase beyond a contrast of 100, reaching significantly higher values than
<figure><img src="image_33.png" /></figure>

### Case of a single sphere and cylinder following y:

* Total volume fraction: 11.42%.

* With a $f_{sp}$ of 5.71%.

* With a $f_{cy}$ of 5.71%, the same as the case of a single cylindrical inclusion for $h = 1$.

* Contrast: 0.001 up to 1000.

* Resolution: $256^3$.

<figure><img src="image_34.png" /><figcaption>**Figure 25.** Illustration of a single sphere and cylinder with a cylinder height of 1 inside a 3D domain.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 26.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single sphere and cylinder with a cylinder height of 1.</figcaption></figure>

$C_{111}$ and $C_{333}$. This pronounced increase in $C_{222}$ indicates a strong directional dependence along the y-axis, suggesting that the material exhibits enhanced stiffness or mechanical resistance in this direction at high contrasts.
Figure 26b illustrates the behavior of the coefficients $C_{113} = C_{331}$ as a function of the logarithm of contrast. Both coefficients follow an identical trend, remaining nearly constant up to a contrast of 100, after which they increase sharply. This sharp increase, particularly after a contrast of 100, suggests that the material becomes highly sensitive to directional variations at high contrast. The nearly identical behavior of $C_{113} = C_{331}$ implies a symmetry in the material's response along different directions under these conditions.

Figure 26c shows the variation of the coefficients ($C_{112} = C_{221} \neq C_{223} = C_{332}$), which remain relatively low up to a contrast of 100. Beyond this threshold, the coefficients increase significantly by approximately a factor of 8, indicating that the material's response becomes anisotropic under high contrast conditions.

<figure><img src="image_36.png" /><figcaption>**Figure 27.** $C^{0,0}$ tensor configuration with a single sphere and cylinder with a cylinder height of 1.</figcaption></figure>

Figure 27 presents the configuration of the $C^{0,0}$ tensor, obtained from the results of Figure 26, where elements of the same color correspond to identical values, indicating symmetry within the tensor. The upper block is color-coded: red, yellow, black, blue, green, and cyan areas represent equal values, demonstrating symmetry between the corresponding tensor elements. The red, yellow, and black colors, however, represent unique values that differ from the others. The lower block of the tensor is entirely composed of zeros, indicating that all elements in this part are null. This configuration, with a non-zero upper block and a lower block filled with zeros, highlights the structured symmetry of the $C^{0,0}$ tensor.

In summary, at low contrast (contrast < 100), the coefficients of all figures remain relatively stable and close to zero or unity, indicating isotropic behavior and minimal directional dependence. At high contrast (contrast > 100), significant changes occur in the coefficients, with some showing sharp increases. These trends highlight the material's transition to anisotropic behavior at high contrast. These results are crucial for understanding material behavior under different contrast conditions.

Finally, we examine the $D_{mnpqr}^{0,0}$ tensor, a second-order tensor, using the $L_\infty$ norm to represent $D_{mnpqr}^{0,0}$ as in the previous case. In Figure 28, the curve indicates that, in the case of a material with both spherical and cylindrical inclusions and a cylindrical height of 1, variations in contrast significantly affect $||D_{mnpqr}^{0,0}||_\infty$. The behavior of the tensor shows an increasing trend as the contrast increases, suggesting that contrast plays a crucial role in determining $||D_{mnpqr}^{0,0}||_\infty$.

The next step of our study focuses on a different geometry, involving the study of three cylindrical inclusions with a cylinder height of 1. This approach allows for the analysis of the influence of the inclusion shape on the properties of the tensors studied. The objective is to better understand the impact of geometry on the material's behavior.
<figure><img src="image_37.png" /><figcaption>**Figure 28.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $C^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a single sphere and cylinder with a cylinder height of 1.</figcaption></figure>

**Three cylindrical inclusions with cylinder height of 1.** In this section, we perform a numerical simulation of a representative volume element (RVE) with cubic dimensions [0, 1]$^3$, containing three cylindrical inclusions that together account for a total volume fraction of 14.36%. The centers of the cylinders are positioned at coordinates ($p_x$ = 0.5, $p_y$ = 0.5, $p_z$ = 0.5) in three-dimensional space. Each inclusion has a radius of $r$ = 0.135 and a height of $h$ = 1. The first cylinder is oriented along the x-axis, with the directional vector ($l_x$ = 0.5, $l_y$ = 0, $l_z$ = 0), the second along the y-axis ($l_x$ = 0, $l_y$ = 0.5, $l_z$ = 0), and the third along the z-axis ($l_x$ = 0, $l_y$ = 0, $l_z$ = 0.5).

<figure><img src="image_38.png" /><figcaption>**Figure 29.** Illustration of three cylinder with a cylinder height of 1 inside a 3D domain</figcaption></figure>

The two graphs presented (see Figure 30) illustrate the variation of the $C^{0,0}$ tensor as a function of the logarithm of contrast for different components. The analysis of the $C^{0,0}$ tensor reveals noteworthy behaviors that underscore the significance of isotropic symmetry in the material consisting of three cylindrical inclusions, with the first oriented along the x-axis, the second along the y-axis, and the third along the z-axis.

For the coefficients $C_{111} = C_{222} = C_{333}$ (see Figure 30a), it is observed that they range from -0.5 to a maximum value of 8 with the logarithm of contrast, yet remain identical across all three components. This indicates the absence of a preferential direction, characteristic of the isotropic
symmetry of the three cylindrical inclusions arranged along the x, y, and z axes. The consistent variation from -0.5 to 8 across both low and high contrast regimes further affirms this isotropic behavior.

Similarly, the other coefficients $C_{112} = C_{113} = C_{221} = C_{223} = C_{331} = C_{332}$ (see Figure 30b) remain constant irrespective of the contrast, demonstrating the uniformity of the material's properties. These findings are of particular relevance to material design, as they show that the material retains its uniform properties despite variations in contrast, which is essential for applications requiring isotropy. Furthermore, a significant increase in the coefficient values of the $C^{0,0}$ tensor is observed beyond a contrast of 100, indicating a marked directional and morphological dependence under high-contrast conditions.

It appears that in the configuration comprising three cylindrical inclusions oriented along the x, y, and z axes, the tensor maintains its isotropic nature but exhibits a substantial increase in coefficient values with increasing contrast, particularly at higher contrast levels. Additionally, the negative values observed in the logarithmic region can be attributed to the initially low values of the coefficients; since the logarithm of a small number is negative, this results in negative values in the graph.

Figure 31 illustrates the configuration of the $C^{0,0}$ tensor, where elements sharing the same color (red and blue) are assigned identical values. The red circles, arranged along the diagonal of the upper block, represent identical values due to the lack of a preferential direction, as the cylinders are positioned along the x, y, and z axes. The blue circles, also identical to each other, correspond to other values in the upper block. The lower block is filled with “0” symbols, signifying null values or inactive elements.

In summary, Figures 30a and 30b confirm a clear relationship between the length of inclusions in a composite material and its mechanical properties, as evidenced by the coefficients of the $C^{0,0}$ tensor. This tensor, which plays a critical role in describing the material's elastic properties, demonstrates a pronounced variation in its coefficients as a function of inclusion length. The elongation of the inclusion significantly influences the material's mechanical properties in the direction of extension, enabling the material to withstand greater stress or deformation along its primary axis.

<figure><img src="image_39.png" /><figcaption>**Figure 30.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast with a three cylinder with a cylinder height of 1.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>**Figure 31.** $C^{0,0}$ tensor configuration with a three cylinder with a cylinder height of 1.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Figure 32.** Logarithm of the absolute value of the ratio between the maximum of the coefficients of the tensor $D^{0,0}$ and the maximum of the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast for a system with three cylindrical inclusions, each having a cylinder height of 1.</figcaption></figure>

Finally, the norm of $D^{0,0}_{mmpqr}$ tensor provides a clear and effective method for quantifying the magnitude of its components. In Figure 32, the curve illustrates that, in the case of three cylindrical inclusions with a height of 1, variations in contrast have a significant impact on $||D^{0,0}||_\infty$. The tensor exhibits an increasing trend as the contrast rises, emphasizing the crucial role of contrast in determining $||D^{0,0}||_\infty$.

In the continuation of this study, a comprehensive analysis of the three previously investigated cases is conducted to draw specific conclusions for each. This analysis is subsequently followed by a general summary, which includes a graphical comparison of the coefficients of the tensor $C^{0,0}$, emphasizing the differences observed across the various cases examined.

**Summary and comparison of last numerical results.** The spherical inclusion is demonstrated to preserve isotropic properties with minimal variations as a function of contrast. In contrast, the cylindrical inclusion exhibits an increasing degree of anisotropy, particularly at high contrast levels, even with a small volume fraction. When both inclusions are combined, the anisotropy persists and is nearly identical to that observed for the cylindrical inclusion alone, with only slight variations. This indicates that the presence of a spherical inclusion does not significantly affect the material's overall behavior when combined with a cylindrical inclusion. These findings
underscore the critical role of inclusion geometry in influencing the mechanical response of materials, especially in high-contrast composites. In conclusion, it is evident that second-order gradient theory is essential for accurately modeling high-contrast composites with cylindrical inclusions, such as long fibers. The remainder of the paper provides a comparative analysis of the $C_{111}$ coefficient of the $C^{0,0}$ tensor, enabling a detailed comparison of the various cases studied.

<figure><img src="image_42.png" /><figcaption>**Figure 33.** Logarithm of the ratio between the coefficients of the tensor $C^{0,0}$ and the coefficient for a contrast of 1, plotted as a function of the logarithm of contrast.</figcaption></figure>

The results obtained (see Figure 33) reveal that spherical inclusions, regardless of their volume fraction, as well as short cylindrical inclusions (see Figure 20), exhibit negligible variations in their coefficients across all contrast levels. In contrast, long cylindrical inclusions display a significant increase in the coefficients at high contrast levels, following the relatively low values observed at lower contrasts. When a spherical inclusion is combined with a long cylindrical inclusion, the variations in the coefficients remain moderate. These findings suggest that spherical inclusions have a minimal influence on the tensor variations, whereas the presence of long cylindrical inclusions leads to a marked increase in the coefficients of the $C^{0,0}$ tensor. Consequently, this confirms the importance of second-order gradient theory in accurately modeling high-contrast composites, particularly those incorporating long cylindrical inclusions, such as long fibers.

Figure 33 clearly illustrates the relationship between the inclusion length in composite materials and their mechanical properties, as evidenced by the coefficients of the $C^{0,0}$ tensor. This tensor, which is crucial for characterizing the elastic behavior of the material, shows substantial variations in its coefficients as a function of the inclusion length. The elongation of the inclusion significantly influences the mechanical properties along the direction of extension, thereby enhancing the material's ability to withstand higher stress or deformation along its principal axis. In the subsequent section, a table summarizing the simulation times for all cases is presented.

### 6.2.3. *Simulation time*

In this final subsection, we aim to highlight the efforts and advancements made in utilizing the MPI and OpenMP libraries for parallel programming, along with the results obtained at a contrast level of 1000. We successfully integrated the MPI/OpenMP libraries into our numerical
simulation framework, achieving a notable performance improvement. Specifically, the time reduction observed was a factor of 4, as summarized in Table 4.

The numerical simulation is executed on the Myria cluster, which is optimized for shared-mode operation. Jobs are submitted in batch through the Slurm job scheduler. Comprehensive details regarding the setup and usage can be found in the “User Guide” section, as cited in [64–68].

<table><caption>**Table 4.** Simulation time.</caption><thead><tr><td>RVE</td><td>Time with MPI/OpenMP (s)</td><td>Time without MPI/OpenMP (s)</td><td>Time gain</td></tr></thead><tbody><tr><td>Spherical inclusion for f<sub>sp</sub> = 5.71%</td><td>769.04</td><td>4228.67</td><td>5.49</td></tr><tr><td>Spherical inclusion for f<sub>sp</sub> = 25%</td><td>904</td><td>4523</td><td>5</td></tr><tr><td>Cylindrical inclusion for h = 0.5</td><td>858</td><td>4547</td><td>5.30</td></tr><tr><td>Cylindrical inclusion for h = 1</td><td>998.34</td><td>4857.35</td><td>5</td></tr><tr><td>Spherical/cylindrical inclusions for h = 0.5</td><td>312.23</td><td>1300.26</td><td>4.16</td></tr><tr><td>Spherical/cylindrical inclusions for h = 1</td><td>349.89</td><td>1400.87</td><td>4</td></tr><tr><td>Three cylindrical inclusions for h = 0.75</td><td>725.55</td><td>4152.85</td><td>5.72</td></tr><tr><td>Three cylindrical inclusions for h = 1</td><td>887.16</td><td>4452.62</td><td>5.02</td></tr></tbody></table>

In this section, we illustrate the homogenization approach based on higher-order gradient theory through numerical results obtained for different three-dimensional cell configurations.

### 6.2.4. Numerical homogenization of periodic microstructures: impact of basic and unit cell configuration on the computation of effective tensors

To illustrate the homogenization approach based on higher-order gradient theory, we present numerical results obtained from several three-dimensional cells [69–72]: a basic unit cell (see Figures 34a and 36a) and supercells with dimensions 2 × 2 × 2 (see Figures 34b and 36b), 3 × 3 × 3 (see Figures 34c and 36c), and 4 × 4 × 4 (see Figures 34d and 36d). Two types of inclusions are considered: spherical and cylindrical. We begin by analyzing the results corresponding to spherical inclusions, followed by those associated with cylindrical inclusions. These simulations allow us to study the mechanical effects induced by first- and second-order stiffness tensors, as well as their influence. The analyses are performed for contrast ratios equal to 10 and 1000, with a volume fraction of 16% for spherical inclusions and 3.84% for cylindrical inclusions.

**Case of a spherical inclusion.** To illustrate the effect of spherical inclusions on the effective mechanical properties of a composite, various three-dimensional periodic representative volume elements (RVEs) of increasing sizes—from the basic unit cell (SP1) to super-cells 2×2×2, 3×3×3, and 4×4×4 (SP2 to SP4)—have been studied, as shown in Figure 34. The classical first-order homogenized tensor $A^{0,0}$ and second-order tensor $B^{0,1}$ have been computed for two contrast ratios $E_i/E_m = 10$ and $E_i/E_m = 1000$, with a volume fraction of 16% for the spherical inclusion, as reported in Tables 5 and 6. It has been observed that the main components $A_{111}$, $A_{2222}$, $A_{3333}$ and $B_{111}$, $B_{2222}$, $B_{3333}$ have remained identical, indicating isotropic effective behavior, and have been found to be nearly constant regardless of the cell size. This stability has suggested that even the basic cell SP1 has been sufficient to represent the periodic microstructure in this case.
However, a slight increase in the components of the tensors $A^{0,0}$ and $B^{0,1}$ have been observed with high contrast (from 1.53533 to 1.54146), indicating a greater sensitivity to gradient effects in high-contrast materials. Finally, in Figure 35, the evolution of the components $C_{111} = C_{222} = C_{333}$ of the tensor $C^{0,0}$, obtained via second-gradient homogenization, has been plotted as a function of the number of cells. The results have shown a rapid convergence, with constant values for all cells. Moreover, a significant increase in effective stiffness has been highlighted at high contrast, confirming that the gradient model has captured non-local mechanical effects induced by the microstructure.

<figure><img src="image_43.png" /><figcaption>**Figure 34.** Basic cells and unit cell in periodic microstructures with a spherical inclusion in 3D.</figcaption></figure>

<table><caption>**Table 5.** Components of tensor $A^{0,0}$ with a volume fraction of 16%.</caption><thead><tr><th></th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 10</th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 1000</th></tr><tr><th></th><th>A<sub>1111</sub></th><th>A<sub>2222</sub></th><th>A<sub>3333</sub></th><th>A<sub>1212</sub></th><th>A<sub>1111</sub></th><th>A<sub>2222</sub></th><th>A<sub>3333</sub></th><th>A<sub>1212</sub></th></tr></thead><tbody><tr><th>SP1</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP1</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP2</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP2</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP3</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP3</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP4</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP4</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr></tbody></table>

In the following section, we focus on the case of a cylindrical inclusion within a three-dimensional framework. This configuration allows us to analyze how the geometric shape of the inclusion influences the effective mechanical properties of the composite material. Using a series of periodic cells similar to those employed in the spherical inclusion case, we evaluate both the classical and higher-order homogenized tensors to compare the mechanical effects induced by a more pronounced geometric anisotropy. This study also aims to highlight the sensitivity of the gradient-based model to variations in the shape and orientation of the inclusion.

**Case of a cylindrical inclusion.** In this section, a microstructure containing a cylindrical inclusion has been considered, modeled using three-dimensional periodic cells of increasing sizes (CY1 to CY4), as illustrated in Figure 36. The classical homogenized tensor $A^{0,0}$, presented in Table 7, has been computed, and it has revealed an effective anisotropy induced by the elongated shape of the cylinder: the components $A_{111}$ and $A_{222}$ have been found to be lower than $A_{333}$, indicating greater stiffness along the cylinder's axis. This trend has been further amplified under high contrast conditions ($E_i/E_m = 1000$), where the overall stiffness values have increased while maintaining the anisotropic structure. The tensor $B^{0,1}$, which accounts for gradient effects and is presented in Table 8, has exhibited values very close to those of $A^{0,0}$, and has preserved the
<figure><img src="image_44.png" /><figcaption>**Figure 35.** Coefficients of the $C^{0,0}$ tensor as a function of the number of cells.</figcaption></figure>

<table><caption>**Table 6.** Components of tensor $\mathbb{B}^{0,1}$ with a volume fraction of 16%.</caption><thead><tr><th></th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 10</th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 1000</th></tr><tr><th></th><th>B<sub>1111</sub></th><th>B<sub>2222</sub></th><th>B<sub>3333</sub></th><th>B<sub>1212</sub></th><th>B<sub>1111</sub></th><th>B<sub>2222</sub></th><th>B<sub>3333</sub></th><th>B<sub>1212</sub></th></tr></thead><tbody><tr><th>SP1</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP1</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP2</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP2</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP3</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP3</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr><tr><th>SP4</th><td>1.53533</td><td>1.53533</td><td>1.53533</td><td>0.38539</td><td>SP4</td><td>1.54146</td><td>1.54146</td><td>1.54146</td></tr></tbody></table>

same anisotropic pattern with strong stability across all cell sizes. This behavior has suggested that even the basic unit cell CY1 has been sufficient to capture the primary effects of this geometry. Finally, Figure 37 has shown the evolution of the components of the tensor $C^{0,0}$, obtained through second-gradient homogenization. It has been observed that $C_{333}$ has been significantly higher than $C_{111}$ and $C_{222}$, confirming increased stiffness along the longitudinal axis of the cylinder. Moreover, the consistency of the results from CY2 onward has confirmed a rapid convergence of the effective tensors. These observations have highlighted the strong influence of inclusion geometry on the homogenized mechanical response and have demonstrated the ability of the gradient-based model to capture non-local effects induced by microstructural anisotropy.

## 7. Conclusion

This paper investigates second-order gradient theories and multi-scale modeling of high-contrast periodic composite materials. In this work, we have demonstrated the ability to efficiently, reliably, and rapidly evaluate and classify these materials for various levels of contrast, resolution, and geometric configurations. Although more complex geometries can be analyzed, their exploration falls beyond the scope of this study.

This paper focuses on the numerical evaluation and simulation of the four tensors $A^{0,0}$, $\mathbb{B}^{0,1}$, $C^{0,0}$, and $\mathbb{D}^{0,0}$ across different geometric configurations. The results obtained provide crucial
<figure><img src="image_45.png" /><figcaption>**Figure 36.** Basic cells and unit cell in periodic microstructures with a cylindrical inclusion in 3D.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 37.** Coefficients of the $C^{0,0}$ tensor as a function of the number of cells.</figcaption></figure>

<table><caption>**Table 7.** Components of tensor $A^{0,0}$ with a volume fraction of 3.84%.</caption><thead><tr><th></th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 10</th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 1000</th></tr><tr><th></th><th>A<sub>1111</sub></th><th>A<sub>2222</sub></th><th>A<sub>3333</sub></th><th>A<sub>1212</sub></th><th>A<sub>1111</sub></th><th>A<sub>2222</sub></th><th>A<sub>3333</sub></th><th>A<sub>1212</sub></th></tr></thead><tbody><tr><th>SP1</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP1</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP2</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP2</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP3</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP3</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP4</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP4</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr></tbody></table>

<table><caption>**Table 8.** Components of tensor $B^{0,1}$ with a volume fraction of 3.84%.</caption><thead><tr><th></th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 10</th><th colspan="4">E<sub>i</sub>/E<sub>m</sub> = 1000</th></tr><tr><th></th><th>B<sub>1111</sub></th><th>B<sub>2222</sub></th><th>B<sub>3333</sub></th><th>B<sub>1212</sub></th><th>B<sub>1111</sub></th><th>B<sub>2222</sub></th><th>B<sub>3333</sub></th><th>B<sub>1212</sub></th></tr></thead><tbody><tr><th>SP1</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP1</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP2</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP2</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP3</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP3</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr><tr><th>SP4</th><td>1.58368</td><td>1.58368</td><td>1.62709</td><td>0.39842</td><td>SP4</td><td>1.60240</td><td>1.60240</td><td>1.68561</td><td>0.40410</td></tr></tbody></table>
insights into the behavior of these tensors across various shapes, offering a deeper understanding of advanced composite materials.

Three-dimensional simulations utilizing the fast Fourier transform (FFT) and the finite element method (FEM) demonstrate the computational efficiency and accuracy of the proposed models in approximating the real properties of composites. The influence of inclusion geometry on the mechanical response of the material is highlighted through numerical simulations. Specifically, spherical inclusions, as well as sets of three cylindrical inclusions oriented along the three principal directions (x, y, z), result in isotropic behavior, whereas cylindrical inclusions introduce pronounced anisotropy and enhance directional stiffness under high-contrast conditions. The geometry of the inclusions was analyzed using three-dimensional simulation software.

The combination of spherical/cylindrical inclusions induces anisotropy. The results reveal a cumulative effect for inclusions embedded within a representative volume element (RVE). In the case of elongated inclusions extending to the RVE boundaries, the mechanical response closely resembles that of cylindrical inclusions alone, emphasizing the negligible influence of spherical inclusions in this context. These findings confirm that inclusion geometry plays a pivotal role in defining the mechanical properties of composite materials and that second-order gradient theory is essential for accurately modeling fiber-reinforced composites, which are extensively utilized in industry.

The in-house computational code developed in C++ for processing tensors associated with this theory has demonstrated its efficiency in enabling parallelized simulations and facilitating the investigation of various morphological and statistical aspects of composite materials. Future enhancements will focus on adapting this code for GPU-based processing, enabling large-scale simulations and the generation of extensive databases, ultimately contributing to the optimization of composite materials for specific applications.

## Appendix A. Representation of tensor exponents

The asymptotic expansion of the energy $W^\epsilon$ takes the form of a series in powers of $\epsilon$, as follows:

$$W^{\epsilon} = W^{0} + \epsilon W^{1} + O(\epsilon). \qquad (38)$$

The expression for the energy $W^0$ (see reference [44]) is given by:

$$W^{0} = \frac{1}{2} \int_{P} E^{0}(x) : A^{(0,0)} : E^{0}(x) dP. \qquad (39)$$

The first exponent of the tensor $A^{(0,0)}$ represents the order of the term on the left-hand side. Since the term is $E^0(x)$, we retain 0. Similarly, the second exponent represents the order of the term on the right-hand side, which is also $E^0(x)$, so we retain 0. Hence, we obtain the tensor $A^{(0,0)}$.

The expression for the energy $W^1$ (see reference [44]) is as follows:

$$W^{1} = \int_{P} E^{0}(x) : \mathbb{B}^{(0,1)} : E^{1}(x) dP + \int_{P} E^{0}(x) : \mathbb{C}^{(0,0)} : \nabla E^{0}(x) dP + \int_{P} E^{0}(x) : \mathbb{D}^{(0,0)} : \nabla E^{0}(x) dP, \quad (40)$$

where:

* $E^0(x) = e_x(u^0) = \epsilon_x(u^0)$;

* $E^1(x) = e_x(U^1) = \epsilon_x(U^1)$;

* $\nabla E^0(x) = e_x(e_x(u^0)) = \epsilon_x(\epsilon_x(u^0))$.

The first exponent of the tensor $\mathbb{B}^{(0,1)}$ represents the order of the term on the left-hand side, which is $E^0(x)$, so we retain 0. The second exponent indicates the order of the term on the right-hand side, which is $E^1(x)$, so we retain 1. Therefore, we obtain the tensor $\mathbb{B}^{(0,1)}$.
The first exponent of the tensor $\mathbb{C}^{(0,0)}$ represents the order of the term on the left-hand side, which is $E^0(x)$, so we retain 0. The second exponent indicates the order of the term on the right-hand side, which is $\nabla E^0(x)$, so we retain 0. Therefore, we obtain the tensor $\mathbb{C}^{(0,0)}$.

The first exponent of the tensor $\mathbb{D}^{(0,0)}$ represents the order of the term on the left-hand side, which is $E^0(x)$, so we retain 0. The second exponent indicates the order of the term on the right-hand side, which is $\nabla E^0(x)$, so we retain 0. Hence, we obtain the tensor $\mathbb{D}^{(0,0)}$.

## Appendix B. Behavioral law

The expression for the energy $W^\epsilon$ is as follows:

$$ W^{\epsilon} = \frac{1}{2} \epsilon(u^0) \mathbb{A}^{(0,0)} \epsilon(u^0) + \epsilon(\epsilon(u^0) \mathbb{B}^{(0,1)} \epsilon(u^0) + \epsilon(u^0) \mathbb{C}^{(0,0)} \partial^2 u^0 + \epsilon(u^0) \mathbb{D}^{(0,0)} \partial^2 u^0) + O(\epsilon). \quad (41) $$

To obtain the stress, $\sigma$ represents the stress tensor (or internal stress tensor) in a material. It describes the relationship between the deformation of the material and the deformation energy density, as follows:

$$ \boldsymbol{\sigma} = \frac{\partial W^{\epsilon}}{\partial \epsilon(u^0)}. \quad (42) $$

Thus, $\sigma$ measures the material's response (the stress) as a function of its strain $\epsilon(u^0)$.

$$ \boldsymbol{\sigma} = (\mathbb{A}^{(0,0)} + 2\epsilon \mathbb{B}^{(0,1)})\epsilon(u^0) + \epsilon(\mathbb{C}^{(0,0)} + \mathbb{D}^{(0,0)})\partial^2 u^0 + O(\epsilon). \quad (43) $$

Thus, $\tau$ is obtained by differentiating the energy $W^\epsilon$ with respect to the second gradient of displacement $\partial^2 u^0$. Physically, this means that $\tau$ is associated with higher-order strain gradient effects, often referred to as hyperstress:

$$ \tau = \frac{\partial W^{\epsilon}}{\partial \partial^2 u^0}, \quad (44) $$

$$ \tau = \epsilon(\mathbb{C}^{(0,0)} + \mathbb{D}^{(0,0)}) + O(\epsilon). \quad (45) $$

With $\sigma$ representing the stress in the composite material and $\tau$ denoting the double stresses, the physical interpretation was provided by Polizzotto [58].

## Appendix C. Algorithm

The algorithm used to compute the homogenized tensor $\mathbb{A}^{0,0}$ is based on the modified Green's accelerated kernel scheme, which is specifically designed to handle high-contrast cases [45,46]. The tensors $\mathbb{B}^{0,1}$, $\mathbb{C}^{0,0}$, and $\mathbb{D}^{0,0}$ are derived through the evaluation of both the displacement gradient and the displacement. The expression for the modified Green's operator is as follows:

The tensor $\widehat{\Gamma}^0_{ijkl}(r)$ is defined as follows:

$$ \widehat{\Gamma}^0_{ijkl}(r) = \frac{(\lambda^0 + 2\mu^0)(r_i r_l^* \delta_{jk})_{\text{sym}} + \lambda^0 [(r_i r_l^* s_{jk})_{\text{sym}} - \text{Re}(r_i r_j^*) \text{Re}(r_k r_l^*)] - \mu^0 r_i r_j r_k^* r_l^*}{\mu^0 [2(\lambda^0 + \mu^0) - \lambda^0 |r_1^2 + r_2^2 + r_3^2|^2]} \quad (46) $$

where:

$$ r_i = \frac{k_i}{|k|}, \quad r_i^* = \frac{k_i^*}{|k^*|} $$

with $k$ and $k^*$ representing the “discrete” gradient and divergence operators, respectively. The tensor $s$ is defined as:

$$ s_{jj} = 4 \operatorname{Im}(r_i r_k^*)^2, \quad s_{jk} = -4 \operatorname{Im}(r_k r_j^*) \operatorname{Im}(r_k r_i^*), \quad \text{for } i \neq j \neq k \neq i, \quad (47) $$

where $\operatorname{Im}(-)$ denotes the imaginary part of the complex quantity.

Additionally, $F_i$ is defined as the inverse of the stiffness matrix.
**Algorithm 1 Elasticity Algorithm Based on an Accelerated Scheme**

1: **Initialize**: $\epsilon^0(\mathbf{x}) \leftarrow E$

2: Set convergence criterion: $\text{acc} \leftarrow 10^{-6}$.

3: **while** no convergence **do**

4: **Convergence Test:**

5: Compute $\boldsymbol{\sigma}(\mathbf{x}) \leftarrow \mathbb{C}(\mathbf{x})\boldsymbol{\epsilon}^n(\mathbf{x})$

6: $\hat{\boldsymbol{\sigma}}(\xi) \leftarrow \text{FFT}(\boldsymbol{\sigma})$

7: $\epsilon_{eq} \leftarrow \frac{\sqrt{\langle\langle \|\xi\hat{\boldsymbol{\sigma}}(\xi)^n\|\rangle\rangle}}{\|\hat{\boldsymbol{\sigma}}(0)\|}$

8: **if** $\epsilon_{eq} < \text{acc}$ **then**

9: **Stop:** Convergence achieved.

10: **end if**

11: Compute $\boldsymbol{\tau}^n(\mathbf{x}) \leftarrow (\mathbb{C}(\mathbf{x}) + \mathbb{C}^0)\boldsymbol{\epsilon}^n(\mathbf{x})$

12: $\hat{\boldsymbol{\tau}}^n(\xi) \leftarrow \text{FFT}(\boldsymbol{\tau}^n)$

13: $\hat{\boldsymbol{\epsilon}}_{\text{comp}}^n(\xi) \leftarrow \hat{\Gamma}^0(\xi)\hat{\boldsymbol{\tau}}^n(\xi)$ for $\xi \neq 0$

14: $\hat{\boldsymbol{\epsilon}}_{\text{comp}}^n(0) \leftarrow E$

15: $\boldsymbol{\epsilon}_{\text{comp}}^n(\mathbf{x}) \leftarrow \text{FFT}^{-1}(\hat{\boldsymbol{\epsilon}}_{\text{comp}}^n)$

16: $\epsilon_{\text{comp}} \leftarrow \frac{\sqrt{\langle\langle \|\boldsymbol{\epsilon}^n - \boldsymbol{\epsilon}_{\text{comp}}^n\|\rangle\rangle^2}}{\|E\|}$

17: $\boldsymbol{\epsilon}^{n+1}(\mathbf{x}) \leftarrow \boldsymbol{\epsilon}^n(\mathbf{x}) - 2(\mathbb{C}(\mathbf{x}) - \mathbb{C}^0)^{-1}\mathbb{C}^0(\boldsymbol{\epsilon}_{\text{comp}}^n(\mathbf{x}) - \boldsymbol{\epsilon}^n(\mathbf{x}))$

18: $\hat{\boldsymbol{\epsilon}}_{\alpha\beta}^n(\chi^0) \leftarrow F_i \cdot \hat{\boldsymbol{\sigma}}_{\alpha\beta}(\xi)$

19: $\boldsymbol{\epsilon}_{\alpha\beta}^n(\chi^0) \leftarrow \text{FFT}^{-1}(\hat{\boldsymbol{\epsilon}}_{\alpha\beta}^n(\chi^0))$

▷ Deformation

20: $\hat{\chi}_{\alpha\beta}^{0n} \leftarrow \frac{\hat{\epsilon}_{\alpha\beta}^{n}(\chi^0)}{i\xi_{\alpha}}$

21: $\chi_{\alpha\beta}^{0n} \leftarrow \text{FFT}^{-1}(\hat{\chi}_{\alpha\beta}^{0n})$

▷ Displacement

22: $\frac{\partial}{\partial y_l} \mathbb{C}_{nijkl}(\boldsymbol{e}_{y_{nj}}(\hat{\boldsymbol{\chi}}_{pq}^{1r}) + \delta_{jr}\hat{\boldsymbol{\chi}}_{pq}^{0n}) = -\mathbb{C}_{n j k r}(\delta_{n p} \delta_{j q} + e_{y_{ij}}(\hat{\boldsymbol{\chi}}_{pq}^{0n}))$

▷ Force $e_y(\chi^1)$

23: $-i\xi_l \mathbb{C}_{nijkl}(\boldsymbol{e}_{y_{nj}}(\hat{\boldsymbol{\chi}}_{pq}^{1r}) + \delta_{jr}\hat{\boldsymbol{\chi}}_{pq}^{0n}) = -\mathbb{C}_{n j k r}(\delta_{n p} \delta_{j q} + e_{y_{ij}}(\hat{\boldsymbol{\chi}}_{pq}^{0n}))$

24: $\boldsymbol{e}_{y_{nj}}(\hat{\boldsymbol{\chi}}_{pq}^{1r}) = \frac{1}{-i\xi_l \mathbb{C}_{nijkl}}(-\mathbb{C}_{n j k r}(\delta_{n p} \delta_{j q} + e_{y_{ij}}(\hat{\boldsymbol{\chi}}_{pq}^{0n})) + i\xi_l \delta_{j r} \hat{\boldsymbol{\chi}}_{pq}^{0n})$

25: $\boldsymbol{e}_{y_{nj}}(\chi_{pq}^{1r}) \leftarrow \text{FFT}^{-1} \boldsymbol{e}_{y_{nj}}(\hat{\boldsymbol{\chi}}_{pq}^{1r})$

26: **end while**

At convergence, the average stress $\langle\langle \boldsymbol{\sigma}(\mathbf{x}) \rangle\rangle$ is calculated, and $\mathbb{A}^{0,0}$ or $\mathbb{C}^{\text{hom}}$ is obtained using the following relation:

$$ \langle\langle \boldsymbol{\sigma}(\mathbf{x}) \rangle\rangle = \mathbb{A}^{0,0} : \langle\langle \varepsilon(\mathbf{x}) \rangle\rangle $$

The tensor $\mathbb{B}^{0,1}$ is determined from the deformation calculated using equation (30), considering the specific indices. Each coefficient of the tensor $\mathbb{B}^{0,1}$ is calculated individually, resulting in a total of 36 coefficients. For example, the process of calculating the coefficient $\mathbb{B}_{1111}$ of the tensor $\mathbb{B}^{0,1}$ is illustrated below.

$$ \mathbb{B}_{1111} = C_{1111}(1 + e_{11}(\chi_{11}^0)) + C_{1122}e_{22}(\chi_{11}^0) \\ + C_{1133}e_{33}(\chi_{11}^0) + C_{1123}e_{23}(\chi_{11}^0) + C_{1131}e_{31}(\chi_{11}^0) + C_{1112}e_{12}(\chi_{11}^0) $$

The tensor $\mathbb{C}^{0,0}$ is determined from the deformation and displacement obtained using equation (31), taking into account the specific indices. Each coefficient of the tensor $\mathbb{C}^{0,0}$ is calculated separately, resulting in a total of 18 coefficients. For example, the process of calculating the coefficient $C_{111}$ of the tensor $\mathbb{C}^{0,0}$ is illustrated below.
The coefficient $C_{111}$ is calculated as follows:

$$
\begin{equation}
\begin{aligned}
C_{111} ={}& C_{1111} \chi_1^{01} + C_{1122} \chi_1^{02} + C_{1133} \chi_1^{03} + C_{1123} \chi_1^{02} + C_{1131} \chi_1^{03} + C_{1112} \chi_1^{01} \\
& + C_{1111} e_{11}(\chi_{11}^0) \chi_1^{01} + C_{1122} e_{11}(\chi_{12}^0) \chi_1^{01} + C_{1133} e_{11}(\chi_{13}^0) \chi_1^{01} + C_{1123} e_{11}(\chi_{21}^0) \chi_1^{01} \\
& + C_{1131} e_{11}(\chi_{31}^0) \chi_1^{01} + C_{1112} e_{11}(\chi_{12}^0) \chi_1^{01} + C_{2211} e_{11}(\chi_{11}^0) \chi_1^{02} + C_{2222} e_{11}(\chi_{22}^0) \chi_1^{02} \\
& + C_{2233} e_{11}(\chi_{33}^0) \chi_1^{02} + C_{2223} e_{11}(\chi_{23}^0) \chi_1^{02} + C_{2231} e_{11}(\chi_{31}^0) \chi_1^{02} + C_{2212} e_{11}(\chi_{12}^0) \chi_1^{02} \\
& + C_{3311} e_{11}(\chi_{11}^0) \chi_1^{03} + C_{3322} e_{11}(\chi_{22}^0) \chi_1^{03} + C_{3333} e_{11}(\chi_{33}^0) \chi_1^{03} + C_{3323} e_{11}(\chi_{23}^0) \chi_1^{03} \\
& + C_{3331} e_{11}(\chi_{31}^0) \chi_1^{03} + C_{3312} e_{11}(\chi_{12}^0) \chi_1^{03} + C_{2311} e_{11}(\chi_{11}^0) \chi_1^{02} + C_{2322} e_{11}(\chi_{22}^0) \chi_1^{02} \\
& + C_{2333} e_{11}(\chi_{33}^0) \chi_1^{02} + C_{2323} e_{11}(\chi_{23}^0) \chi_1^{02} + C_{2331} e_{11}(\chi_{31}^0) \chi_1^{02} + C_{2312} e_{11}(\chi_{12}^0) \chi_1^{02} \\
& + C_{3111} e_{11}(\chi_{11}^0) \chi_1^{03} + C_{3122} e_{11}(\chi_{22}^0) \chi_1^{03} + C_{3133} e_{11}(\chi_{33}^0) \chi_1^{03} + C_{3123} e_{11}(\chi_{23}^0) \chi_1^{03} \\
& + C_{3131} e_{11}(\chi_{31}^0) \chi_1^{03} + C_{3112} e_{11}(\chi_{12}^0) \chi_1^{03} + C_{1211} e_{11}(\chi_{11}^0) \chi_1^{02} + C_{1222} e_{11}(\chi_{22}^0) \chi_1^{02} \\
& + C_{1233} e_{11}(\chi_{33}^0) \chi_1^{02} + C_{1223} e_{11}(\chi_{23}^0) \chi_1^{02} + C_{1231} e_{11}(\chi_{31}^0) \chi_1^{02} + C_{1212} e_{11}(\chi_{12}^0) \chi_1^{02}.
\end{aligned}
\tag{48}
\end{equation}
$$

The tensor $D^{0,0}$ is determined from the deformation and volumetric forces, as obtained through equation (32), while considering the specific indices. Each coefficient of the tensor $D^{0,0}$ is calculated individually, resulting in a total of 243 coefficients. As an example, the procedure for calculating the coefficient $D_{1111}$ of the tensor $D^{0,0}$ is demonstrated below:

$$
\begin{equation}
\begin{aligned}
D_{1111} ={}& C_{1111}(1 + e_{11}(\chi_{11}^0))(e_{11}(\chi_1^{11})) + C_{2211}(e_{22}(\chi_{11}^0))(e_{11}(\chi_1^{11})) + C_{3311}(e_{33}(\chi_{11}^0))(e_{11}(\chi_1^{11})) \\
& + C_{2311}(e_{23}(\chi_{11}^0))(e_{11}(\chi_1^{11})) + C_{3111}(e_{22}(\chi_{11}^0))(e_{11}(\chi_1^{11})) + C_{1211}(e_{12}(\chi_{11}^0))(e_{11}(\chi_1^{11})).
\end{aligned}
\tag{49}
\end{equation}
$$

## Appendix D. $C^{0,0}$ tensor representation

$$
C^{0,0} = \begin{bmatrix} C_{111} & C_{112} & C_{113} \\ C_{221} & C_{222} & C_{223} \\ C_{331} & C_{332} & C_{333} \end{bmatrix} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}
$$

## Appendix E. $D^{0,0}$ tensor representation

The tensor $D_{mnpqr}^{0,0}$ has been able to be represented in matrix form, where $m$ has been assigned the values 1, 2, and 3.

$$
D_{1npqr}^{0,0} = \begin{bmatrix}
D_{1111} & D_{1121} & D_{1131} & D_{11112} & D_{11212} & D_{11312} & D_{11113} & D_{11213} & D_{11313} \\
D_{1112} & D_{1122} & D_{1132} & D_{11122} & D_{11222} & D_{11322} & D_{11113} & D_{11223} & D_{11323} \\
D_{1113} & D_{1123} & D_{1133} & D_{11113} & D_{11232} & D_{11332} & D_{11113} & D_{11233} & D_{11333} \\
D_{1211} & D_{1221} & D_{1231} & D_{12112} & D_{12212} & D_{12312} & D_{12113} & D_{12213} & D_{12313} \\
D_{1212} & D_{1222} & D_{1232} & D_{12122} & D_{12222} & D_{12322} & D_{12123} & D_{12223} & D_{12323} \\
D_{1213} & D_{1223} & D_{1233} & D_{12132} & D_{12232} & D_{12332} & D_{12133} & D_{12233} & D_{12333} \\
D_{1311} & D_{1321} & D_{1331} & D_{13112} & D_{13212} & D_{13312} & D_{13113} & D_{13213} & D_{13313} \\
D_{1312} & D_{1322} & D_{1332} & D_{13122} & D_{13222} & D_{13322} & D_{13123} & D_{13223} & D_{13323} \\
D_{1313} & D_{1323} & D_{1333} & D_{13132} & D_{13232} & D_{13332} & D_{13133} & D_{13233} & D_{13333}
\end{bmatrix}
$$
$$D_{2npqr}^{0,0} = \begin{bmatrix} D_{21111} & D_{21211} & D_{21311} & D_{21112} & D_{21212} & D_{21312} & D_{21113} & D_{21213} & D_{21313} \\ D_{21121} & D_{21221} & D_{21321} & D_{21122} & D_{21222} & D_{21322} & D_{21123} & D_{21223} & D_{21323} \\ D_{21131} & D_{21231} & D_{21331} & D_{21132} & D_{21232} & D_{21332} & D_{21133} & D_{21233} & D_{21333} \\ D_{22111} & D_{22211} & D_{22311} & D_{22112} & D_{22212} & D_{22312} & D_{22113} & D_{22213} & D_{22313} \\ D_{22121} & D_{22221} & D_{22321} & D_{22122} & D_{22222} & D_{22322} & D_{22123} & D_{22223} & D_{22323} \\ D_{22131} & D_{22231} & D_{22331} & D_{22132} & D_{22232} & D_{22332} & D_{22133} & D_{22233} & D_{22333} \\ D_{23111} & D_{23211} & D_{23311} & D_{23112} & D_{23212} & D_{23312} & D_{23113} & D_{23213} & D_{23313} \\ D_{23121} & D_{23221} & D_{23321} & D_{23122} & D_{23222} & D_{23322} & D_{23123} & D_{23223} & D_{23323} \\ D_{23131} & D_{23231} & D_{23331} & D_{23132} & D_{23232} & D_{23332} & D_{23133} & D_{23233} & D_{23333} \end{bmatrix}$$

$$D_{3npqr}^{0,0} = \begin{bmatrix} D_{31111} & D_{31211} & D_{31311} & D_{31112} & D_{31212} & D_{31312} & D_{31113} & D_{31213} & D_{31313} \\ D_{31121} & D_{31221} & D_{31321} & D_{31122} & D_{31222} & D_{31322} & D_{31123} & D_{31223} & D_{31323} \\ D_{31131} & D_{31231} & D_{31331} & D_{31132} & D_{31232} & D_{31332} & D_{31133} & D_{31233} & D_{31333} \\ D_{32111} & D_{32211} & D_{32311} & D_{32112} & D_{32212} & D_{32312} & D_{32113} & D_{32213} & D_{32313} \\ D_{32121} & D_{32221} & D_{32321} & D_{32122} & D_{32222} & D_{32322} & D_{32123} & D_{32223} & D_{32323} \\ D_{32131} & D_{32231} & D_{32331} & D_{32132} & D_{32232} & D_{32332} & D_{32133} & D_{32233} & D_{32333} \\ D_{33111} & D_{33211} & D_{33311} & D_{33112} & D_{33212} & D_{33312} & D_{33113} & D_{33213} & D_{33313} \\ D_{33121} & D_{33221} & D_{33321} & D_{33122} & D_{33222} & D_{33322} & D_{33123} & D_{33223} & D_{33323} \\ D_{33131} & D_{33231} & D_{33331} & D_{33132} & D_{33232} & D_{33332} & D_{33133} & D_{33233} & D_{33333} \end{bmatrix}$$

## Acknowledgments

Most of the computations detailed in this paper have been carried out at the cluster of the Center of Informatics Resources of Higher Normandy (CRIANN – Centre Régional Informatique et d'Applications Numériques de Normandie).

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## References

[1] S. Lemaitre, Modélisation des matériaux composites multiphasiques à microstructures complexes. Étude des propriétés effectives par des méthodes d'homogénéisation, PhD thesis, Université de CAEN Normandie (France), 2017.

[2] A. C. Eringen, "Nonlocal polar elastic continua", *Int. J. Eng. Sci.* **10** (1972), no. 1, pp. 1-16.

[3] M. Lazar and G. Po, "On Mindlin's isotropic strain gradient elasticity: Green tensors, regularization, and operator-split", *J. Micromech. Mol. Phys.* **3** (2018), no. 03n04, article no. 1840008.

[4] M. Lazar, E. Agiasofitou and T. Böhlke, "Mathematical modeling of the elastic properties of cubic crystals at small scales based on the Toupin-Mindlin anisotropic first strain gradient elasticity", *Contin. Mech. Thermodyn.* **34** (2022), no. 1, pp. 107-136.

[5] G. Po, N. C. Admal and M. Lazar, "The Green tensor of Mindlin's anisotropic first strain gradient elasticity", *Mater. Theory* **3** (2019), pp. 1-16.

[6] P. Cartraud, S. Bourgeois and M. Sharafaty-Zangeneh, "Homogénéisation d'ordre supérieur de composites élastiques", 2005. 7e Colloque National en Calcul des Structures.

[7] J. L. Auriault and D. Caillerie, "Quelques remarques sur les méthodes d'homogénéisation", *Rev. Franc. Géotech.* (1989), no. 49, pp. 43-50.

[8] B. Durand, Modélisation numérique et caractérisation expérimentale de matériaux architecturés au comportement élastique de second gradient, PhD thesis, École des Ponts ParisTech (France), 2022.

[9] J.-C. Michel, H. Moulinec and P. Suquet, "Effective properties of composite materials with periodic microstructure: a computational approach", *Comput. Methods Appl. Mech. Eng.* **172** (1999), no. 1-4, pp. 109-143.
[10] B. Gommers, I. Verpoest and P. Van Houtte, "The Mori–Tanaka method applied to textile composite materials", *Acta Mater.* **46** (1998), no. 6, pp. 2223–2235.

[11] R. Christensen, H. Schantz and J. Shapiro, "On the range of validity of the Mori–Tanaka method", *J. Mech. Phys. Solids* **40** (1992), no. 1, pp. 69–73.

[12] T. Mori and K. Tanaka, "Average stress in matrix and average elastic energy of materials with misfitting inclusions", *Acta Metall.* **21** (1973), no. 5, pp. 571–574.

[13] T.-H. Tran, V. Monchiet and G. Bonnet, "A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media", *Int. J. Solids Struct.* **49** (2012), no. 5, pp. 783–792.

[14] M. Schneider, "Convergence of FFT-based homogenization for strongly heterogeneous media", *Math. Methods Appl. Sci.* **38** (2015), no. 13, pp. 2761–2778.

[15] H. Moulinec, P. Suquet and G. W. Milton, "On the convergence test of FFT-based methods", in *ECCOMAS Congress 2016 VII European Congress on Computational Methods in Applied Sciences and Engineering*, 2016.

[16] H. Moulinec and P. Suquet, "A fast numerical method for computing the linear and nonlinear mechanical properties of composites", *C. R. Méc.* **318** (1994), pp. 1417–1423.

[17] H. Moulinec and P. Suquet, "A numerical method for computing the overall response of nonlinear composites with complex microstructure", *Comput. Methods Appl. Mech. Eng.* **157** (1998), no. 1-2, pp. 69–94.

[18] M. Schneider, "A review of nonlinear FFT-based computational homogenization methods", *Acta Mech.* **232** (2021), no. 6, pp. 2051–2100.

[19] S. Lucarini, M. V. Upadhyay and J. Segurado, "FFT based approaches in micromechanics: fundamentals, methods and applications", *Model. Simul. Mat. Sci. Eng.* **30** (2021), no. 2, article no. 023002 (98 pages).

[20] M. Bornert, T. Bretheau and P. Gilormini, *Homogénéisation en mécanique des matériaux*, Hermès Science, 2001, p. 250.

[21] E. Kröner, "Elasticity theory of materials with long range cohesive forces", *Int. J. Solids Struct.* **3** (1967), no. 5, pp. 731–742.

[22] J. A. Krumhansl, "Some considerations of the relation between solid state physics and generalized continuum mechanics", in *Mechanics of Generalized Continua: Proceedings of the IUTAM-Symposium on the Generalized Cosserat Continuum and the Continuum Theory of Dislocations with Applications, Freudenstadt and Stuttgart (Germany) 1967*, Springer, 1968, pp. 298–311.

[23] A. C. Eringen and D. G. B. Edelen, "On nonlocal elasticity", *Int. J. Eng. Sci.* **10** (1972), no. 3, pp. 233–248.

[24] R. H. J. Peerlings, *Enhanced damage modelling for fracture and fatigue*, PhD thesis, Eindhoven University of Technology (The Netherlands), 1999.

[25] R. H. J. Peerlings, R. de Borst, W. A. M. Brekelmans, J. H. P. de Vree and I. Spee, "Some observations on localisation in non-local and gradient damage models", *Eur. J. Mech. A Solids* **15** (1996), pp. 937–953.

[26] R. H. J. Peerlings, R. de Borst, W. A. M. Brekelmans and J. H. P. de Vree, "Gradient enhanced damage for quasi-brittle materials", *Int. J. Numer. Methods Eng.* **39** (1996), no. 19, pp. 3391–3403.

[27] R. H. J. Peerlings, M. G. D. Geers, R. de Borst and W. A. M. Brekelmans, "A critical comparison of nonlocal and gradient-enhanced softening continua", *Int. J. Solids Struct.* **38** (2001), no. 44-45, pp. 7723–7746.

[28] A. E. Green and R. S. Rivlin, "Multipolar continuum mechanics", *Arch. Ration. Mech. Anal.* **17** (1964), pp. 113–147.

[29] H. Askes and L. J. Sluys, "Explicit and implicit gradient series in damage mechanics", *Eur. J. Mech. A Solids* **21** (2002), no. 3, pp. 379–390.

[30] R. A. Toupin, "Elastic materials with couple-stresses", *Arch. Ration. Mech. Anal.* **11** (1962), no. 1, pp. 385–414.

[31] R. A. Toupin, "Theories of elasticity with couple-stress", *Arch. Ration. Mech. Anal.* **17** (1964), no. 2, pp. 85–112.

[32] R. D. Mindlin, "Second gradient of strain and surface-tension in linear elasticity", *Int. J. Solids Struct.* **1** (1965), no. 4, pp. 417–438.

[33] K. Pham, *Construction et analyse de modèles d'endommagement à gradient*, PhD thesis, Université Pierre et Marie Curie–Paris VI (France), 2010.

[34] K. Pham, H. Amor, J.-J. Marigo and C. Maurini, "Gradient damage models and their use to approximate brittle fracture", *Int. J. Damage Mech.* **20** (2011), no. 4, pp. 618–652.

[35] K. Pham and J.-J. Marigo, "From the onset of damage to rupture: construction of responses with damage localization for a general class of gradient damage models", *Contin. Mech. Thermodyn.* **25** (2013), pp. 147–171.

[36] A. Benallal and J.-J. Marigo, "Bifurcation and stability issues in gradient theories with softening", *Model. Simul. Mat. Sci. Eng.* **15** (2006), no. 1, S283–S295.

[37] H. A. Anziz, *Homogénéisation de composites élastiques périodiques à fort contraste: conception de métamatériaux de second gradient*, PhD thesis, Université de Toulon (France), 2018.

[38] C. Boutin, J. Soubestre and M. Dietz, "Milieux à double gradient pour matériaux renforcés. Modélisation théorique et expérimentale", in *CFM 2011-20ème Congrès Français de Mécanique*, AFM, 2011, 6 pages.
[39] R. D. Mindlin and N. N. Eshel, "On first strain-gradient theories in linear elasticity", *Int. J. Solids Struct.* **4** (1968), no. 1, pp. 109–124.

[40] D.-K. Trinh, *Méthodes d'homogénéisation d'ordre supérieur pour les matériaux architecturés*, PhD thesis, École Nationale Supérieure des Mines de Paris (France), 2011.

[41] T.-H. Tran, *Comportement homogénéisé des matériaux composites: prise en compte de la taille des éléments microstructuraux et des gradients de la déformation*, PhD thesis, Université Paris-Est (France), 2013.

[42] S. B. E. Salah, A. Nait-Ali, M. Gueguen and C. Nadot-Martin, "Non-local modeling with asymptotic expansion homogenization of random materials", *Mech. Mater.* **147** (2020), article no. 103459 (18 pages).

[43] H. Moulinec and P. Suquet, "Comparison of FFT-based methods for computing the response of composites with highly contrasted mechanical properties", *Phys. B: Condens. Matter* **338** (2003), no. 1-4, pp. 58–60.

[44] S. B. E. Salah, *Modélisation non-locale et stochastique de matériaux à fort gradient de propriétés par développement asymptotique*, PhD thesis, ISAE-ENSMA Ecole Nationale Supérieure de Mécanique et d'Aérotechnique, Poitiers (France), 2019.

[45] F. Willot, "Fourier-based schemes for computing the mechanical response of composites with accurate local fields", *C. R. Méc.* **343** (2015), no. 3, pp. 232–245.

[46] A. Finel, "A fast and robust discrete FFT-based solver for computational homogenization", 2024. Online at https://arxiv.org/abs/2405.11168.

[47] M. François, "Tenseurs en mécanique", 2018. Lecture notes "M2 Mécanique et fiabilité des structures".

[48] M. François, G. Geymonat and Y. Berthaud, "Determination of the symmetries of an experimentally determined stiffness tensor: application to acoustic measurements", *Int. J. Solids Struct.* **35** (1998), no. 31-32, pp. 4091–4106.

[49] S. Forte and M. Vianello, "Symmetry classes and harmonic decomposition for photoelasticity tensors", *Int. J. Eng. Sci.* **35** (1997), no. 14, pp. 1317–1326.

[50] A. Bensoussan, J.-L. Lions and G. Papanicolaou, *Asymptotic analysis for periodic structures*, AMS Chelsea Publishing, Providence, RI, 2011, pp. xii+398. Corrected reprint of the 1978 original.

[51] J.-L. Lions, "Asymptotic expansions in perforated media with a periodic structure", *Rocky Mt. J. Math.* **10** (1980), no. 1, pp. 125–140.

[52] E. Sanchez-Palencia, *Non-homogeneous media and vibration theory*, Lecture Notes in Physics, Springer, 1980, p. 406.

[53] E. Sanchez-Palencia, "Comportements local et macroscopique d'un type de milieux physiques heterogenes", *Int. J. Eng. Sci.* **12** (1974), no. 4, pp. 331–351.

[54] J. Sanchez-Hubert and E. Sanchez-Palencia, *Introduction aux méthodes asymptotiques et à l'homogénéisation*, Masson, 1992, p. 280.

[55] B. Miara and E. Sanchez-Palencia, "Asymptotic analysis of linearly elastic shells", *Asymptotic Anal.* **12** (1996), no. 1, pp. 41–54.

[56] R. A. Toupin, "Elastic materials with couple-stresses", *Arch. Ration. Mech. Anal.* **11** (1962), no. 1, pp. 385–414.

[57] R. D. Mindlin, "Micro-structure in linear elasticity", *Arch. Ration. Mech. Anal.* **16** (1964), pp. 51–78.

[58] C. Polizzotto, "A note on the higher order strain and stress tensors within deformation gradient elasticity theories: Physical interpretations and comparisons", *Int. J. Solids Struct.* **90** (2016), pp. 116–121.

[59] G. E. Exadaktylos and I. Vardoulakis, "Microstructure in linear elasticity and scale effects: a reconsideration of basic rock mechanics and rock fracture mechanics", *Tectonophysics* **335** (2001), no. 1-2, pp. 81–109.

[60] M. Lazar and G. A. Maugin, "Nonsingular stress and strain fields of dislocations and disclinations in first strain gradient elasticity", *Int. J. Eng. Sci.* **43** (2005), no. 13-14, pp. 1157–1184.

[61] T.-H. Tran, V. Monchiet and G. Bonnet, "A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media", *Int. J. Solids Struct.* **49** (2012), no. 5, pp. 783–792.

[62] O. Débordes, "Homogenization computations in the elastic or plastic collapse range applications to unidirectional composites and perforated sheets", in *Proceedings of the 4th International Symposium Innovative Numerical Methods in Engineering*, Springer, 1986, pp. 453–458.

[63] A. E. Green and R. S. Rivlin, "Multipolar continuum mechanics", *Arch. Ration. Mech. Anal.* **17** (1964), pp. 113–147.

[64] CRIANN, *Centre Régional Informatique et d'Applications Numériques de Normandie (CRIANN)*. Online at https://www.criann.fr (accessed on October 1, 2024).

[65] CRIANN, *Guide d'utilisation du calculateur Myria*. Online at https://services.criann.fr/services/calcul/cluster-myria/utilisation (accessed on October 1, 2024).

[66] CRIANN, *Guide d'utilisation du cluster Austral*. Online at https://services.criann.fr/services/hpc/cluster-austral/guide (accessed on October 1, 2024).

[67] CRIANN, *Modules*. Online at https://services.criann.fr/services/hpc/cluster-austral/guide/modules (accessed on October 1, 2024).
[68] CRIANN, *Description de la configuration*. Online at https://services.criann.fr/services/hpc/cluster-austral/architecture (accessed on October 1, 2024).

[69] S. Barboura and J. Li, "Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures", *Int. J. Solids Struct.* **136** (2018), pp. 60–76.

[70] H. Yang, B. E. Abali, W. H. Müller, S. Barboura and J. Li, "Verification of asymptotic homogenization method developed for periodic architected materials in strain gradient continuum", *Int. J. Solids Struct.* **238** (2022), article no. 111386 (19 pages).

[71] V. Monchiet, N. Auffray and J. Yvonnet, "Strain-gradient homogenization: a bridge between the asymptotic expansion and quadratic boundary condition methods", *Mech. Mater.* **143** (2020), article no. 103309 (22 pages).

[72] W. J. Drugan and J. R. Willis, "A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites", *J. Mech. Phys. Solids* **44** (1996), no. 4, pp. 497–524.