New RVE concept and FFT methods in micromechanics of composites subjected to body force with compact support

Abstract

We consider static linear elastic composite materials (CMs) with periodic structure. The core of the proposed methodology is the generation of a novel dataset using specially designed body force fields with compact support (BFCS), enabling a new RVE concept that reduces the infinite periodic medium to a finite domain without boundary artifacts. This functionally reduced RVE is used for translated averaging of DNS results, efficiently computed via a newly developed fast Fourier transform (FFT)-based solver for BFCS loading. The resulting dataset captures localized field responses and is used to train machine learning (ML) and neural networks (NN) models to learn effective nonlocal surrogate operators. These operators accurately predict macroscopic responses while reflecting microstructural features and nonlocal interactions. By accounting for field localization while simultaneously eliminating influences from finite sample size and boundary effects, it provides a physically grounded and data-driven framework for constructing accurate surrogate models for the homogenization of complex materials.

Keywords

Microstructures inhomogeneous material non-local methods multiscale modeling fast Fourier transform

1. Introduction

Periodic composites, due to their inherent microstructural regularity, are well-suited for multiscale homogenization frameworks [1 –3]. In asymptotic homogenization, originally developed by Babuška and further formalized in Bakhvalov and Panasenko [4] and Fish [1], the response of a periodic medium is approximated by separating scales under the assumption that the unit cell (UC) is much smaller than the overall material size that leads to homogenized coefficients via unit-cell problems. Conversely, computational homogenization methods resolve the microstructural field equations numerically, enabling the modeling of nonlocal [5] and inelastic behavior [6 –8]. A prominent framework is the FE² scheme, where the macroscopic finite element problem is coupled to microscale RVEs embedded at each Gauss integration point [9 –11]. Each RVE solves a boundary value problem consistent with macroscopic deformation modes, typically under periodic boundary conditions (PBCs). While FE² provides a rigorous scale-bridging strategy, its computational cost is significant due to the nested finite element discretizations, often requiring high-performance computing or model order reduction techniques for tractability in large-scale simulations.

Fast Fourier transform (FFT)-based numerical homogenization techniques, first introduced by Moulinec and Suquet [12,13], offer a computationally efficient alternative to finite element methods (FEMs), reducing complexity from $O (N^{2})$ to $O (N \log N)$ via spectral solvers that leverage the discrete Fourier transform (DFT). The original scheme relies on the solution of the periodic Lippmann–Schwinger (L-S) integral equation, reformulated in Fourier space where convolution with the Green’s operator becomes an element-wise multiplication, enabling efficient iterative solvers. FFT methods are generally grouped into three technical categories: (1) LS-type solvers: These include polarization-based iterative schemes and Krylov-subspace methods (e.g., conjugate gradient (CG), biconjugate gradient (BiCG)) for solving the linearized form of the L-S equation [14 –18]. Recent developments have incorporated machine learning (ML)—e.g., FFT-neural network (NN) hybrids for multiscale prediction in complex woven composites [19]. (2) Fourier–Galerkin methods: Derived from a variational (weak-form) framework, these approaches utilize trigonometric polynomial basis functions with consistent projection operators to enhance convergence and accuracy, particularly for materials with discontinuities or high contrast [20,21]. (3) Displacement-based FFT solvers: These formulations treat the displacement fluctuation field as the primary unknown and solve equilibrium equations directly in Fourier space, avoiding the need for a reference medium and facilitating the use of preconditioned Krylov solvers [22]. Advanced FFT-based solvers are now used in highly nonlinear regimes (e.g., elastoplasticity [23]) and in the context of crystal plasticity [24]. Many implementations are open-source and exploit domain decomposition and parallel FFT libraries to scale across distributed memory architectures [25], achieving optimal $O (N \log N)$ performance in large-scale simulations.

The representative volume element (RVE) is fundamental for predicting the effective behavior of heterogeneous materials. It must be large enough to capture the essential microstructural features while ensuring that macroscopic responses are independent of boundary conditions (BCs) and representative of the bulk material. Hill’s classical definition [26] requires macroscopically uniform BCs and effective properties described by homogenized moduli. RVE size selection involves achieving scale separation, where the microstructural length scale a satisfies $a ≪ Λ ≪ L$ , with Λ being the applied field scale and L the macroscale domain length. The minimal domain size where effective properties stabilize is taken as the RVE. Related to this is the concept of the statistically equivalent RVE (SERVE), which uses micro-computed tomography (micro-CT) imaging and statistical analysis to construct realistic computational domains (for detailed methodologies, see [7,27 –30]).

When scale separation is violated, statistical homogeneity breaks down, leading to nonlocal coupling between stress and strain fields. This interaction is governed by a tensorial kernel, requiring the use of effective nonlocal operators—either integral or higher-order differential forms—instead of classical effective moduli [26]. These operators capture the influence of distant points in the material. Nonlocal models fall into two categories: strongly nonlocal (such as strain-based or displacement-based methods like peridynamics) and weakly nonlocal (e.g., strain- or stress-gradient theories). This shift necessitates redefining the RVE to correspond with the chosen nonlocal operator, applicable to both random [31 –33] and periodic [34 –37] composites. These generalized RVEs are defined as the minimal domain size at which a predefined effective (nonlocal) operator—unlike the classical effective moduli of Hill [26]—stabilizes. The RVE becomes even more critical when accounting for combined nonlocal effects from BCs, intrinsic material behavior, and phase interactions.

The evolution of effective nonlocal operator theory has been greatly accelerated by the integration of ML and NN techniques, introducing new levels of flexibility and modeling power. Initial efforts by Silling [38] and You et al. [39,40] demonstrated how DNS data could be used to construct surrogate integral operators for complex materials. More recently, nonlocal neural operators have emerged as tools for learning mappings between function spaces [41,42]. Numerous neural operator architectures have been developed, including DeepONet, PCA-Net, Graph Neural Operators, FNO, and LNO, each tailored to different aspects of operator learning. Comparative reviews can be found in the literature [42 –45]. For nonlocal mechanics, the Peridynamic Neural Operator (PNO) [46] and its heterogeneous extension HeteroPNO [47] provide physics-aware modeling of peridynamic interactions. Physics-Informed Neural Networks (PINNs) further enhance NN models by embedding governing equations directly as soft constraints [44,48,49], ensuring consistency with physical laws. When neural operators are combined with PINNs [43,50,51], the resulting frameworks can accurately model complex nonlinear, heterogeneous, and nonlocal material behavior with strong generalization capabilities.

While ML and NN methods have advanced material modeling significantly, they frequently miss key micromechanical considerations—such as scale separation, boundary effects, and the RVE—which are essential for accurate predictions in both linear and nonlinear regimes. To bridge this gap, the proposed approach constructs new types of compressed effective datasets for complex microstructures (random or periodic), based on an innovative RVE framework. This novel RVE concept is independent of the constitutive behavior of individual phases and the analytical form of surrogate operators. Instead, it leverages field concentration factors within each phase to characterize the microstructure effectively. The resulting effective datasets, enriched by this micromechanically informed RVE, should be compatible with any ML or NN architecture for predicting nonlocal surrogate operators. This innovative RVE concept ensures the accuracy of predictions by removing potential issues related to size scale, boundary layers, and edge effects.

The proposed approach is composed of several key components, some of which are fundamentally novel and mark a significant departure from traditional techniques. At the core of the methodology is the generation of a new type of effective dataset, derived through the application of specifically designed body force fields with compact support (BFCS). These body forces are spatially localized, allowing for precise control of excitation within the material sample. This feature enables the formulation of a new RVE concept, in which the infinite periodic microstructure is effectively reduced to a finite computational domain, without sacrificing the fidelity of microstructural response or introducing boundary-related artifacts. This RVE reduction is not merely geometric, but functional: the effective dataset is built through translated averaging of direct numerical simulations (DNSs) performed within this finite domain. To enable efficient simulation, a newly developed FFT-based solver tailored for CMs is employed. This FFT solver is adapted to handle the response of materials subjected to body forces with compact support, leading to fast, scalable, and accurate computation of the local fields. Once this enhanced effective dataset is constructed—containing information that is rich in both spatial resolution and microstructural mechanics—it serves as a training ground for ML and NN models. These models are trained to learn the effective surrogate nonlocal operators that can predict the material response under arbitrary macroscopic loading conditions. Importantly, the surrogate operator constructed via this framework reflects both the nonlocal interactions inherent in the material and the fine-scale features encoded in the effective dataset. Therefore, the proposed methodology effectively integrates rigorous micromechanical principles with advanced data-driven modeling techniques. By accounting for field localization while simultaneously eliminating influences from finite sample size and boundary effects, it provides a physically grounded and broadly applicable framework for constructing accurate surrogate models for the homogenization of complex materials.

The structure of the paper is organized as follows. Section 2 introduces both the classical and modified forms of the L-S equations. Section 3 presents the traditional and newly proposed RVE concepts, along with the methodology for selecting the corresponding effective datasets. A brief overview of the FFT framework, adapted to support later developments, is provided in Section 4. In Section 5, new FFT-based methods are introduced for analyzing composite materials (CMs) under BFCS. Finally, Section 6 details how the newly generated effective datasets are integrated into existing ML and NN approaches for constructing surrogate models.

2. Modified Lippmann–Schwinger equation

We consider a linear elastic body occupying an open simply connected bounded domain $w \subset R^{d}$ with a smooth boundary $Γ_{0}$ and with an indicator function W and space dimensionality d ( $d = 2$ and $d = 3$ for 2D and 3D problems, respectively). The domain w with the boundary $Γ^{0}$ contains a homogeneous matrix $v^{(0)}$ and a periodic field $X = (v_{i})$ of heterogeneity $v_{i}$ with the centers $x_{i}$ , indicator functions $V_{i}$ and bounded by the closed smooth surfaces $Γ_{i}$ $(i = 1, 2, \dots)$ . It is presumed that the heterogeneities can be grouped into phases $v^{(q)} (q = 1, 2, \dots, N)$ with identical mechanical and geometrical properties. We first consider the local basic equations of thermoelasticity of composites

\begin{matrix} \nabla \cdot σ (x) & = & - b (x), \end{matrix}

(1)

\begin{matrix} σ (x) & = & L (x) ε (x), or ε (x) = M (x) σ (x), \end{matrix}

(2)

\begin{matrix} ε (x) & = & \nabla^{s} u, \nabla \times ε (x) \times \nabla = 0, \end{matrix}

(3)

where ⊗ and × are the tensor and vector products, respectively, $\nabla^{s}$ is the symmetric gradient operator, $\nabla^{s} u : = [\nabla \otimes u + {(\nabla \otimes u)}^{⊤}] ∕ 2$ , and ${(.)}^{⊤}$ denotes the matrix transposition; b is the body force. L(x) and ${M(x) \equiv L(x)}^{- 1}$ are the known phase stiffness and compliance tensors. In particular, for isotropic constituents, the local stiffness tensor $L (x)$ is given in terms of the local bulk modulus $k (x)$ and the local shear modulus $μ (x)$ and:

L (x) = (dk, 2 μ) \equiv dk (x) N_{1} + 2 μ (x) N_{2},

$N_{1} = δ \otimes δ ∕ d, N_{2} = {I - N}_{1}$ $(d = 2 or 3$ ) whereas δ and I are the unit second-order and fourth-order tensors. For all material tensors g (e.g., $L, M)$ , the notation $g_{1} (x) \equiv g (x) - g^{(0)} = g_{1}^{(m)} (x)$ $(x \in v^{(m)}, m = 0, 1, \dots, N$ ) is used. The upper index ^(m) indicates the components, and the lower index i shows the individual heterogeneity; $v^{(0)} = w ∖ v$ , $v \equiv \cup v^{(k)} \equiv \cup v_{i}, V (x) = \sum V^{(k)} = \sum V_{i} (x)$ , and $V^{(k)} (x)$ and $V_{i} (x)$ are the indicator functions of $v^{(k)}$ and $v_{i}$ , respectively, equals 1 at $x \in v^{(m)}$ and 0 otherwise, $(m = 0, k; k = 1, 2, \dots, N; i = 1, 2, \dots)$ .

Substituting equations (2) and (3₁) into equation (1) leads to the equilibrium equation (1) being expressed in the form:

ℒ (u) (x) + b (x) = 0, ℒ (u) (x) : = \nabla [L \nabla u (x)],

(4)

with $ℒ (u) (x)$ representing a second-order elliptic differential operator.

The body force density $b (x)$ is assumed to have compact support, be self-equilibrated called BFCS, and vanishes outside a specified loading region $ℬ^{b} : = b (0, B^{b})$ :

\int b (x) d x = 0, b (y) \equiv 0 for y \notin b (0, B^{b}) : = {y | | y | \leq B^{b}},

(5)

where $b (0, B^{b})$ denotes a ball of radius $B^{b}$ centered at the origin $x = 0$ .

A linear-elastic reference material is introduced, characterized as homogeneous and isotropic with a stiffness tensor $L^{(0)}$ . Consider the governing equation for an infinite homogeneous medium occupying $ℝ^{d}$ ( $d = 1, 2, 3$ ), subject to the body force density $b (x)$ defined in equation (5):

ℒ^{(0)} (u^{b (0)}) (x) + b (x) = 0,

(6)

where $ℒ^{(0)}$ denotes the elliptic operator associated with the homogeneous stiffness tensor $L^{(0)}$ .

The BFCS $b (x)$ (5) is assumed to be self-equilibrated and periodic with respect to a body force unit cell (BFUC) $Ω_{00}^{b} \supset ℬ^{b}$ . For notational simplicity in describing periodic BFCS, we restrict our attention to the 2D case, where the entire domain is represented as a union of square UCs, $w = \cup Ω_{ij}^{b}$ ( $i, j = 0, \pm 1, \pm 2, \dots$ ), with the grid of centers $Λ^{b} = {x^{b Λ}}$ . Let $Ω_{00}^{b}$ denote a representative BFUC, bounded by corner points $x_{kl}^{bc}$ ( $k, l = \pm 1$ ), and with boundary $Γ^{b 0} = \cup Γ_{ij}^{b 0}$ . Each segment $Γ_{ij}^{b 0}$ separates the central cell $Ω_{00}^{b}$ from its adjacent neighbor $Ω_{ij}^{b}$ , where the indices $(i, j)$ satisfy $i = 0, \pm 1$ and $j = \pm (1 - | i |)$ (see Figure 19.1 in Buryachenko [52]). The representative BFUC $Ω_{00}^{b}$ undergoes deformation in the same repetitive manner as all neighboring cells, ensuring periodicity of the body force:

b (x - χ) = b (x), χ \in Λ^{b} .

(7)

This force distribution induces a displacement field given by ( $x \in Ω_{00}^{b}$ )

u^{b (0)} (x) \equiv - {(ℒ^{(0)})}^{- 1} b .

(8)

Alternatively, this displacement can be represented via the periodic Green operator $G^{(0)} (x)$ corresponding to the Navier equation (4) for the homogeneous reference tensor $L^{(0)}$ :

u^{b (0)} (x) = \int G^{(0)} (x - y) b (y) d y .

(9)

Loosely speaking, the Green operator $G^{(0)} (x)$ may be interpreted as the inverse of the reference stiffness operator, characterizing the response of the infinite medium to localized force distributions.

We now proceed to the consideration of a CM whose local stiffness is given by $L (x) = L^{(0)} + L_{1} (x)$ , where $L^{(0)}$ is the stiffness tensor of a homogeneous reference medium and $L_{1} (x)$ represents the spatial variations.

The total displacement field $u (x)$ can then be decomposed as:

u (x) = u^{b (0)} (x) + u_{1} (x),

(10)

where $u^{b (0)} (x)$ solves the homogeneous problem with the reference modulus:

\begin{matrix} ℒ^{(0)} (u^{b (0)}) (x) & = & - b (x), \end{matrix}

(11)

\begin{matrix} ℒ (u) (x) & = & ℒ^{(0)} (u^{b (0)}) . \end{matrix}

(12)

Substituting these into the equilibrium equation (2) yields:

\nabla L^{(0)} \nabla u_{1} (x) = - \nabla L_{1} (x) \nabla u (x),

(13)

which leads to an implicit integral representation of the strain field using a Green operator called a modified Lippmann–Schwinger (ML-S) equation:

ε (x) = ε^{b (0)} (x) + \int U^{(0)} (x - y) τ (y) d y,

(14)

where $ε^{b (0)} (x) = \nabla^{s} u^{b (0)} (x)$ is the symmetric gradient of the reference displacement, and $τ (y) : = L_{1} (y) ε (y)$ is the polarization tensor. The kernel $U^{(0)}$ is the second derivative of the reference Green tensor $G^{(0)}$ , i.e., $G^{(0)}$ : $U_{ijkl}^{(0)} (x) = \nabla_{j} \nabla_{l} G_{(ij) (kl)}^{(0)}$ which vanishes as $| x | \to \infty$ . Here, the symmetrization in the lower indices (denoted by parentheses) ensures the tensor satisfies symmetry requirements of elasticity. It should be noted that the solutions to equations (11) and (12), given in the form of the volume integral representations (8) and (14), respectively, can also be obtained by alternative methods, such as the FEM.

It should be noted that the ML-S equation (14) represents a linear, locally elastic counterpart of the Additive General Integral Equations (AGIEs) proposed for CMs with arbitrary microstructures—random (statistically homogeneous or inhomogeneous), periodic, or deterministic (neither random nor periodic). These AGIEs are formulated to accommodate a broad spectrum of material behaviors, including linear and nonlinear responses, coupled or uncoupled multiphysical effects, and local, weakly nonlocal (strain/stress gradient), or strongly nonlocal (peridynamic) phase properties.

The solution to the integral equation (14) can be expressed in terms of a Neumann series expansion, which iteratively approximates the solution by a convergent sequence of successive operator applications. This representation is valid under conditions ensuring the contraction property of the associated integral operator—typically guaranteed when the contrast between material properties is sufficiently small or the norm of the perturbation operator is less than one. Formally, the Neumann series takes the form:

ε = [\sum_{k = 0} {(U^{(0)} * L_{1})}^{k}] ε^{b (0)} .

(15)

The strain Eshelby–Green tensor $U^{(0)} (x)$ is known explicitly in the Fourier domain for an infinite, isotropic reference medium with Lamé coefficients $λ^{(0)}$ and $μ^{(0)}$ (see, e.g., [12,53]):

\begin{matrix} U^{(0)} & = & \frac{1}{4 μ^{(0)}} U_{1} + \frac{λ^{(0)} + μ^{(0)}}{μ^{(0)} (λ^{(0)} + 2 μ^{(0)})} U_{2}, \end{matrix}

(16)

\begin{matrix} U_{1 | ijkl} (ζ) & = & - | ζ |^{- 2} (δ_{ki} ζ_{h} ζ_{j} + δ_{hi} ζ_{k} ζ_{j} + δ_{kj} ζ_{h} ζ_{i} + δ_{hj} ζ_{k} ζ_{i}), \end{matrix}

(17)

\begin{matrix} U_{2 | ijkl} (ζ) & = & | ζ |^{- 4} ζ_{i} ζ_{j} ζ_{k} ζ_{h}, \end{matrix}

(18)

It is also worth noting that for any square-integrable field $τ (x)$ , the convolution $(U^{(0)} * τ) (0) = 0$ , ensuring consistency with the equilibrium condition.

Equation (14) is supplemented by the PBCs at the BFUC $Ω_{00}^{b}$ :

ε (x) #, σ (x) \cdot n (x) = - #, x \in \partial Ω_{00}^{b}

(19)

where the strain field $ε (x)$ is assumed to be periodic (denoted by $ε #$ ), while the corresponding traction vector $σ (x) \cdot n (x)$ is anti-periodic (denoted by $σ \cdot n = - #$ ). These conditions ensure compatibility and equilibrium across the boundaries of adjacent periodic cells. The PBC (19) at $b (x) \equiv 0$ are equivalent to the homogeneous remote BCs (also called the kinematic uniform boundary conditions (KUBC) and static uniform boundary conditions (SUBC), respectively) with some symmetric constant tensors either $ε^{w_{Γ}}$ or $σ^{w_{Γ}}$

\begin{matrix} u (y) & = & ε^{w_{Γ}} y, \forall y \in Γ_{0 u} = Γ_{0}, \\ t (y) & = & σ^{w_{Γ}} n (y), \forall y \in Γ_{0 σ} = Γ_{0}, \end{matrix}

(20)

correspond to the analyses of the equations for either strain or stresses, respectively, which are formally similar to each other.

It is further assumed that the CM $w = \cup Ω_{ij}$ ( $i, j = 0, \pm 1, \pm 2, \dots)$ exhibits periodic microstructure, with a UC of geometry $Ω_{ij}$ and a grid of centers $Λ = {x_{ij}}$ of $Ω_{ij}$ , and that the body force $b (x)$ is defined over an enlarged domain $Ω_{00}^{b}$ such that representative UC $Ω_{00} \subset Ω_{00}^{b}$ :

L (x - χ) = L (x), χ \in Λ .

(21)

Under these assumptions, the ML-S equation (14), along with the PBC (19), remains valid over $Ω_{00}^{b}$ . However, the standard PBCs on the smaller reference unite cell (UC) $Ω_{00}$ :

ε (x) #, σ (x) \cdot n (x) = - #, x \in \partial Ω_{00}

(22)

are generally violated due to the larger support of the forcing term $b (x)$ , which extends beyond $Ω_{00}$ . As a result, equilibrium between adjacent UCs $Ω_{ij}$ is not automatically satisfied, and the problem must be formulated on the extended domain $Ω_{ij}^{b}$ to ensure consistency.

We consider two distinct problems. The first involves equation (4) subjected to BFCS (5) on an isolated domain $Ω_{00}^{b}$ with a free edge. The second involves equation (4) in an infinite periodic medium subjected to the periodic body force (7) with PBCs (19) (or equivalently, remote BCs (20)). At the end of Section 5, we establish the conditions under which the solutions to these problems coincide within the isolated domain $ε (x) \in Ω_{00}^{b}$ and the representative domain $Ω_{00}^{b}$ of the periodic medium.

We now consider the classical L-S problem in the absence of body forces, i.e., $b \equiv 0$ ( $x \in R^{d}$ ), with PBCs applied on the boundary of the UC $Ω_{00}$ , as defined in equation (22). In this setting, we adopt an alternative decomposition of the displacement field, replacing equation (10) with:

u (x) = ⟨ ε ⟩ \cdot x + u_{1} (x),

(23)

where $⟨ ε ⟩$ is a prescribed constant macroscopic strain, and $u_{1} (x)$ is a periodic fluctuation field. Substituting this decomposition into the governing equations yields the same differential identity as in equation (13), albeit with a different fluctuation field $u_{1} (x)$ . This leads directly to the classical L-S (see [12,13]) integral equation:

ε (x) = ⟨ ε ⟩ + \int U^{(0)} (x - y) τ (y) d y .

(24)

Here, $U^{(0)} (x)$ refers to the strain-based Green operator previously introduced in equation (15) (see [52 –54]). It is worth noting that some works use $Γ^{(0)} = - U^{(0)}$ instead of $U^{(0)}$ (see [12,13]). It is noteworthy that the ML-S equation (14) is simpler than the classical L-S equation (24) where its free term $⟨ ε ⟩$ (24) is defined by the averaged solution. In contrast, the free term $ε^{b} (x)$ in equation (14) is a deterministic function that must be estimated in advance using equation (9).

In this formulation (24), the total strain field $ε (x)$ is periodic, and the macroscopic strain $⟨ ε ⟩$ serves as the driving term. The solution to the integral equation (24) can be represented using a Neumann series

ε = [\sum_{k = 0} {(U^{(0)} * L_{1})}^{k}] ⟨ ε ⟩,

(25)

the properties of which are detailed in Milton [55]. The operator representations on the right-hand sides of equations (15) and (25) are identical and, consequently, are subject to the same conditions of convergence.

The method (25) is equivalent to an explicit operator recurrence procedure, initiated from the expression $ε^{1} = ⟨ ε ⟩$ :

ε^{[i + 1]} = ⟨ ε ⟩ + U^{(0)} * (L_{1} ε^{[i]}) .

(26)

Michel et al. [56] observed that the strain Green’s operator satisfies the identity:

U^{(0)} * (L^{(0)} ε) = ⟨ ε ⟩ - ε .

(27)

Substituting equation (27) into equation (26) yields an equivalent recurrence relation in the form:

ε^{[i + 1]} = ε^{[i]} + U^{(0)} * (L ε^{[i]}) with ε^{[1]} = ⟨ ε ⟩ .

(28)

We observe that the recurrence counterpart (28) of the classical L-S equation (24) is reduced to its equivalent recursive representation (28) by applying equation (27), under the condition that $L_{1} (x) \equiv L^{(0)}$ .

We note that the ML-S equation (14) can be viewed as a generalization of the classical form equation (24), in which the constant macroscopic strain $⟨ ε ⟩$ is replaced by a non-uniform strain field $ε^{b (0)} (x)$ that possesses compact support, mirroring the localized nature of the body force $b (x)$ introduced in equation (5).

3. Classical and new RVE concepts

The concept of the RVE, originally introduced by Hill [26], has a long and complex history, often marked by debate and reinterpretation (for a detailed account, see Ostoja-Starzewski et al. [30]). To faithfully convey the essence of Hill’s original definition and its foundational importance, we cite directly from Hill’s [26] seminal work, which provides a rigorous basis for the RVE concept (at $b (x) \equiv 0).$

Definition 2.1. RVE (a) is structurally entirely typical of the whole mixture on average, and (b) contains a sufficient number of inclusions for the apparent overall moduli to be effectively independent of the surface values of traction and displacement, so long as these values are “macroscopically uniform.” The contribution of this surface layer to any average can be negligible by taking the sample large enough.

The primary role of the RVE is to provide an estimation of the effective moduli $L^{*}$ , which in turn relies on the evaluation of phase-averaged field quantities. For CMs with periodic microstructures, the RVE coincides with the Unit Cell (RVE ≡UC), and the application of remote homogeneous BCs (20) naturally reduces to PBCs (22). Although Hill [26] originally formulated the RVE concept in the context of linearly locally elastic CMs, his Definition 2.1 notably omits the term “constitutive law,” allowing for a natural extension of the concept to nonlinear and nonlocal material models.

A representative response can be obtained via DNS of microstructural volume elements (MVEs), either synthetically generated or experimentally extracted (e.g., by micro-CT). The homogenized overall properties $L_{KUBC}^{A}$ and $M_{SUBC}^{A}$ are computed under kinematically (20₁) and statically (20₂) uniform BCs, respectively. While generally different, their discrepancy

L_{KUBC}^{A} - {(M_{SUBC}^{A})}^{- 1} \to 0,

(29)

decreases with increasing MVE size, enabling estimation of the RVE size and the effective moduli $L^{*}$ . The RVE concept assumes sufficiently large—ideally infinite—material domains (see, e.g., [30,52]). It should be emphasized that, although increasing the number of realizations in this intuitive EVE framework can achieve statistical convergence in estimating $L^{*}$ (i.e., the mean sample response), it does not inherently remove edge-related scale effects. In practice, intuitive RVEs (“sufficiently large sample”) are typically constructed using two main approaches: (1) simulated random inclusion fields and (2) image-based models derived from micro-CT scans of real materials (see, e.g., [57 –59]). A pictorial interpretation of the RVE concept is illustrated in Figure 1, which shows a micro-CT scan of a heterogeneous material subjected to remote homogeneous loading, as defined by equation (20₁).

Figure 1.

CT image with remote BC (20).

Under the application of PBC (22) (or (20₁)), the estimation of effective moduli $L^{*}$ (the “upscaling” procedure) and the determination of field concentration factors $A (x)$ for points $x \in v_{i}$ (the “downscaling” procedure) are expressed as:

\begin{matrix} L^{*} = L^{(0)} + R^{*}, ⟨ τ ⟩ = R^{*} ⟨ ε ⟩, ⟨ ε ⟩_{i} (x) = A^{*} (x) ⟨ ε ⟩ \end{matrix}

(30)

highlighting the intrinsic coupling between upscaling and downscaling—two complementary aspects of the same homogenization framework.

Under the PBC as defined in equation (22), it follows from equations (30) that both the effective moduli $L^{*}$ and the effective strain field ${⟨ ε ⟩}_{i} (x)$ are independent of the particular choice of UC $Ω_{ij}$ . However, this invariance breaks down in the presence of a body force applied to the special BFUC $Ω_{00}^{b}$ , where the PBC at $Ω_{00}$ (22) is no longer satisfied. In this case, the resulting strain field $ε (x)$ , expressed in the coordinate system associated with the body force $b (x)$ (with $x \in Ω_{00}^{b}$ ), acquires a dependence on the specific UC $Ω_{ij}^{b}$ containing the point x.

To eliminate this dependency and restore generality, we analyze the behavior of periodic CMs under a fixed body force while allowing the microstructure (i.e., the grid of UC centers $x_{m} \in Λ$ ) to undergo a rigid translation. For concreteness, we restrict the discussion to the linear elasticity model defined by equations (2), and the microstructural strain solution is denoted by $ε_{0} (x)$ . When the grid of UC centers undergoes a rigid translation by a vector χ (i.e., $Λ_{0} \to Λ_{χ}$ ), the material stiffness $L (x, χ)$ (2), the phase indicator function $V_{i} (x, χ)$ , and the resulting strain field solution $ε (x, χ)$ governed by equations (1)–(3) will also transform accordingly

\begin{matrix} L (x, χ) & = & L_{0} (x - χ), V_{i} (x, χ) = V_{i 0} (x - χ), \end{matrix}

(31)

\begin{matrix} b (x, χ) & = & b_{0} (x), u (x, χ) ≢ u_{0} (x - χ) . \end{matrix}

(32)

The inequality in equation (32₂) holds because the body force field $b (x)$ is fixed.

Accordingly, for every translation vector $χ \in V_{x}$ , the solution $u (x, χ)$ derived from equations (32) enables the definition of an effective (macroscopic) displacement field over the domain w

\begin{matrix} ⟨ u ⟩ (x) & = & \frac{1}{{\bar{V}}_{x}} \int_{V_{x}} u (x, χ) d χ, ⟨ u ⟩^{l (1)} (x) = \frac{1}{{\bar{V}}_{x}} \int_{V_{x}} u (x, χ) V_{i} (x, χ) d χ, \end{matrix}

(33)

\begin{matrix} ⟨ σ ⟩ (x) & = & \frac{1}{{\bar{V}}_{x}} \int_{V_{x}} σ (x, χ) d χ, ⟨ σ ⟩^{(1)} (x) = \frac{1}{{\bar{V}}_{x}} \int_{V_{x}} σ (x, χ) V_{i} (x, χ) d χ . \end{matrix}

(34)

It is interesting that equations (33₁) and (34₁) formally looks as averaging over the moving-window $V_{x}$ although the operations (33₁) and (34₁) are conceptually different and obtained by averaging over the number of the displacement realizations $u (x, χ)$ produced by a parallel transform of $L (x, χ)$ (32) rather than by averaging of one realization $u_{0} (x)$ over the moving-window. The averages defined in equations (33) and (34) are more accurately described as ensemble averages, taken over the statistical set of all translated configurations of a given periodic microstructure, where the translations χ are uniformly distributed within the periodicity cell $V_{x}$ . This translation-based averaging procedure is general and applies to periodic CMs with arbitrary constitutive behavior and under any inhomogeneous loading conditions. In so doing, the solution ${\cdot} (x, χ)}$ for each $χ \in V_{x}$ can be obtained by any suitable numerical method (e.g., finite element analysis (FEA) or FFT, see Section 5). A particular instance of such an averaging scheme—specifically, in the form of equation (34₁)—was introduced in the context of asymptotic homogenization by Smyshlyaev and Cherednichenko [37] and later discussed by Ameen et al. [34]. Notably, Smyshlyaev and Cherednichenko [37] attribute this idea to a personal communication with J.R. Willis, in the context of periodic media subjected to periodic loading. Interestingly, the statistical averages in (33₂), and (34₂) can be interpreted as a probabilistic reformulation of the classical student problem: “What is the probability that a randomly dropped coin (representing an inclusion $v_{i}$ ) lands on a fixed point $x \in R^{d}$ ?” In the specific cases of equations(33₂) and (34₂), the “coin” (i.e., inclusion $v_{i}$ ) belongs to a periodically structured grid $Λ_{χ}$ undergoing random translations, further emphasizing the ensemble nature of the averaging. An illustration of the averaging scheme described by equation (34) for the 1D case is shown in Figure 2. In this example, a fixed macroscopic point $X \in b (0, B^{b})$ corresponds to a set of local points $x_{j} = X - χ$ located within the translated UC grids $Λ_{χ} = Λ_{j}$ for $j = 1, \dots, 6$ . The estimation of the average stress field $⟨ σ ⟩ (X)$ is carried out by summing the stress values $σ (X, χ) = σ (x_{j})$ across all realizations of the translated grid $Λ_{j}$ . In contrast, the average stress within the inclusions, denoted by ${⟨ σ ⟩}^{(1)} (X)$ , is computed by summing only over those realizations where the point X falls inside an inclusion—specifically, the realizations corresponding to $j = 1, 2, 3, 6$ . Notably, due to nonlocal effects, it is possible to observe non-zero average stress values $⟨ σ ⟩ (Y) \neq 0$ at points $Y \notin b (0, B^{b})$ , even though the body force vanishes at those points, i.e., $b (Y) = 0$ . This nonlocal average $⟨ σ ⟩ (Y)$ is again computed by summing $σ (X, χ) = σ (x_{j})$ over all six realizations ( $j = 1, \dots, 6$ ), while the inclusion-specific average ${⟨ σ ⟩}^{(1)} (Y)$ is based only on realizations $j = 1, 2, 6$ , where Y lies within an inclusion. The use of capital letters for the coordinates X and Y highlights their macroscopic nature.

Figure 2.

Scheme of translated averaging.

After the averaging (33) and (34) over the grid translation $Λ χ$ , we can introduce an effective dataset $D^{p}$ for periodic-structured CMs as follows:

\begin{matrix} D^{p} & = & {D_{k}^{p}}_{k = 1}^{N}, D_{k}^{p} = {⟨ u_{k} ⟩ (b_{k}, x), ⟨ σ_{k} ⟩ (b_{k}, x), \\ ⟨ ε_{k} ⟩^{(1)} (b_{k}, x), ⟨ σ_{k} ⟩^{(1)} (b_{k}, x), b_{k} (x)}, \end{matrix}

(35)

where each set of effective parameters is computed for a given body force realization $b_{k}$ using the averaging procedures defined in equations (33) and (34) during the offline stage. The coordinates $x \in ℝ^{d}$ represent macroscopic spatial variables; hence, the effective dataset $D^{p}$ provides macroscopic-level information only. Although each effective sub-dataset $D_{χ k}^{p}$ , corresponding to a specific grid $Λ χ$ , is generated through DNS, the complete effective dataset $D_{k}^{p}$ is obtained using the translated averaging procedures defined in equations (33) and (34). This averaging technique is a central component of computational analytical micromechanics (CAM), as defined in Buryachenko [5,60]). Furthermore, the effective datasets $D^{r}$ for random-structure CMs, also derived by Buryachenko [60] using CAM, are formally equivalent to the periodic effective dataset $D^{p}$ given in equation (35). As a result, the estimation procedure for $D^{p}$ in equation (35) is likewise referred to as the CAM method. We use the displacement field $u (x)$ in equation (35) (see also equation (37)), only because the source publications for Section 6 are given as $(u (x), b (x))$ rather than the more typical $(ε (x), b (x))$ ; absent this data restriction, one could equally replace $u (x)$ with the strain field $ε (x)$ .

The representation for the dataset $D^{p}$ (35) was obtained for either a finite size body force UC $Ω_{00}^{b}$ or a full size $Ω_{00}^{b} = R^{d}$ with finite size compact support $ℬ^{b}$ (5) (i.e., $B^{b} < \infty$ ). Reformulation and generalization of the classical definition by Hill [26] 2.1 enables one to formulate a flexible definition sufficient for our current interests with self-equilibrated body force $b (x)$ (5).

Definition 2.2. RVE is structurally entirely typical of the whole CM area, which is sufficient for all apparent effective dataset $D_{k}^{p}$ ( $k = 1, \dots, N$ ) (35) to be effectively stabilized outside $x \notin RVE$ (i.e., the strains and stresses vanish in $x \in \bar{RVE} : = Ω_{00}^{b} ∖ RVE$ ) in the infinite periodic structure CMs.

A critical issue in the analysis of periodic CMs is the appropriate choice of PBC (22) at the interface of UCs, which provide a connection to field distributions, which are the neighboring UCs. PBC (22) is strictly valid under homogeneous remote loading conditions (20)and zero body force $b (x) \equiv 0$ (5). However, in the more general case, where a nonzero body force field is present, the standard PBC given by equation (22) becomes inaccurate and may not reflect the true mechanical response of the material. In such scenarios, it is important to note that imposing PBC at the UC level is not strictly necessary, especially when DNS is employed. The RVE size is treated as a learned (or adjustable) parameter, selected so as to satisfy the tolerance criterion

\begin{matrix} | σ (x) |, | ε (x) | < tol for \forall x \in \bar{RVR} : = R^{d} ∖ RVE \end{matrix}

(36)

ensuring that the geometrical and mechanical representativeness of the RVE is sufficient for accurate homogenization. Moreover, owing to the compact support of the applied body forces ${b_{k} (x)}_{k = 1}^{N}$ , the original problem posed on the infinite domain $ℝ^{d}$ for periodic CM is effectively reduced to a finite domain corresponding to the RVE. This reduction significantly simplifies the analysis. As a result, complications typically associated with finite sample size or boundary-induced (edge) effects are inherently avoided.

Definition 2.2 is conceptually distinct from Definition 2.1. The latter pertains to an idealized, asymptotically large (infinite) domain, where effective material parameters are obtained through a formal limiting process for a “sufficiently large” sample. In contrast, Definition 2.2 introduces a more practical concept—a finite-size RVE—serving as an initial approximation that is amenable to computational or experimental implementation and subsequent refinement. The key distinction, however, lies in the fundamentally different remote BCs (22) and BFCS (6) employed in Definitions 2.1 and 2.2, respectively. This definition considers a heterogeneous medium periodically structured occupying the entire space $ℝ^{d}$ . Similarly, no reference is made to an explicit “effective moduli” (or “effective nonlocal operator,” see for references [5]). Rather, attention is focused on the exterior domain $\bar{RVE} = ℝ^{d} ∖ RVE$ , where effective dataset ${D_{k}^{p} (x)}_{k = 1}^{N}$ (35) reaches a stabilized state. Stabilization of the effective dataset $D_{k}^{p}$ (35) along with the proper selection of the RVE implies that all effective parameters within the annular region bounded by $| x | = B^{RVE}$ and $| x | = B^{RVE} + B^{b} ∕ 2$ remain constant within a prescribed tolerance (36). When this condition is met, the external region $| x | > B^{RVE} + B^{b} ∕ 2$ can be excluded from the simulation, allowing the infinite medium to be accurately modeled using a finite-sized sample. In this way, a correctly chosen RVE eliminates edge effects—commonly encountered as boundary-layer artifacts, especially when $⟨ ε ⟩ (x)$ remains nonzero near the domain boundary. This phenomenon is discussed, for example, on Buryachenko [54]. Conversely, if the RVE is not appropriately selected (i.e., $B^{RVE}$ is insufficiently large), the subsequent application of the effective dataset $D$ (35) (as discussed in Section 6) will result in numerical errors originating from both finite sample size and residual edge effects.

Figure 3 illustrates the spatial arrangement $ℬ^{b} \subset RVE \subset Ω_{00}^{b}$ for CMs with two types of microstructures: periodic (Figure 3(a)) and deterministic (Figure 3(b)). In Figure 3(a), we depict a representative periodic configuration X of inclusions, characterized by the set of centers $Λ_{χ}$ . Although periodicity underlies the inclusion layout, it is not directly used in simulations. Instead, datasets $D_{χ k}^{pj}$ are generated by applying various realizations of BFCS loading (2.6), and evaluating the local response within the BFUC $Ω_{00}^{b}$ using a suitable numerical method (e.g., FEA of FFT-based solvers from Section 5). Importantly, stress and strain fields are assumed to vanish in the exterior of the RVE ( $\bar{RVE}$ ), allowing the treatment of a finite set of inclusions $v_{i} \subset Ω_{00}^{b}$ without explicitly enforcing periodicity. The periodic structure $Λ_{χ}$ is leveraged only during the translation averaging process (see equations (33) and (34)). In a similar manner, Figure 3(b) represents deterministic microstructures $X^{j}$ ( $j = 1, 2, \dots$ ), which are neither periodic nor stochastic and may, for instance, correspond to distinct samples obtained from computed tomography (CT) imaging. As in the periodic case, effective datasets $D_{k}^{d j}$ are computed within $Ω_{00}^{b}$ for different deterministic configurations $X^{j}$ and fixed body force fields $b_{k} (x)$ . These yield an ensemble-based effective dataset defined (as (35)) by

\begin{matrix} D^{d} & = & {D_{k}^{d}}_{k = 1}^{N}, D_{k}^{d} = {⟨ u_{k} ⟩ (b_{k}, x), ⟨ σ_{k} ⟩ (b_{k}, x), \\ ⟨ ε_{k} ⟩^{(1)} (b_{k}, x), ⟨ σ_{k} ⟩^{(1)} (b_{k}, x), b_{k} (x)}, \end{matrix}

(37)

where statistical averaging $⟨ \cdot ⟩$ is performed over the set of deterministic configurations $X^{j}$ , capturing the effective material response to each BFCS $ℬ^{b}$ (5) realization. No specific assumptions are made about the geometry of the regions $ℬ^{b}$ or the RVE itself (cf. Definition 2.2) nor about their relative sizes. Spherical shapes shown in Figure 3 are adopted for clarity and visualization purposes only, without implying any essential geometric constraint.

Figure 3.

Schemes of $ℬ^{b} \subset RVE \subset Ω_{00}^{b}$ for CM with periodic (a) and deterministic (b) structures.

A particularly significant contribution to computational micromechanics is presented in the recent study by Silling et al. [61], which proposes a novel framework for analyzing 2D random composites. The authors examine a composite microstructure generated via Monte Carlo simulations within a finite square domain w, where each phase is modeled using nonlocal elasticity theory (peridynamics). The studied sample, containing approximately 900 circular inclusions (see Figure 4(a)) would conventionally be considered an ideal candidate for implementing the classical RVE concept, as defined in Definition 2.1. In a notable departure from this traditional approach, the authors introduce a self-equilibrated BFCS, $b (x)$ , as described in equation (5). This loading scheme is consistent with the methodologies developed in Buryachenko [62,63] and fulfills the criteria outlined in the generalized RVE concept of Definition 2.2. Specifically, the distance between the RVE boundary and the outer boundary of the computational domain satisfies $dist (\partial RVE, ∂w) \approx 10 a = 100 λ$ , where λ is the characteristic spacing of the microstructural lattice. The effectiveness of this approach is substantiated by the nearly vanishing strain field $ε (x)$ observed near the external boundary ∂w, as shown in the violaceous region in Figure 4(b). This outer region corresponds to $\bar{RVE} = w ∖ RVE$ , where strain localization is negligible.

Figure 4.

(a) Simulated structure in w. (b) DNS of strains $ε (x)$ in w.

Although the explicit term “RVE” is not used in Silling et al. [61], the core conceptual structure clearly aligns with the generalized RVE framework, as represented by the central multicolored (non-violaceous) domain, which functions as the effective RVE. To the best of the author’s knowledge, Figure 4(b) provides the clearest illustration of the new RVE concept. We now consider numerical parametric analyses modeling the simultaneous interaction of three types of nonlocal effects arising from inhomogeneous applied fields $b (x)$ , intrinsic material nonlocality (characterized by the peridynamic horizon $l_{δ}$ ), and interactions among inclusions. These analyses also serve to validate the main concepts introduced in the present study, namely the AGIEs (14), the new RVE concept (see Definition 2.2), and the effective dataset (35). Specifically, a numerical interpretation of Figure 3 was presented in Buryachenko [62,63] for a randomly inhomogeneous bar subjected to a self-equilibrated body force (5), where $b (x) = 0$ for $| x | > B^{b}$ and $b (x) = - b (- x)$ for $| x | \leq B^{b}$ ; $c^{b} = 0$ (see Figure 5). In Figure 5, which shows the effective displacement $⟨ u ⟩ (X)$ vs. $x ∕ B^{b}$ , curve 4 corresponds to the case $l_{δ} = 0$ (the classical engineering approach), while curves 1–3 are obtained for fixed parameters $l_{δ} ∕ a = 1$ and $c^{(1)} = 0.5$ . The ratio $B^{b} ∕ a$ takes the values $0.25, 0.375, 0.5$ , corresponding to curves 1, 2, and 3, respectively.

Figure 5.

Average displacement $⟨ u ⟩ (x)$ vs. $x ∕ a^{b}$ .

Owing to nonlocal effects, the extent of long-range interactions is not limited solely by the condition $| x | \leq B^{b}$ , but instead extends over a larger domain $| x | \leq a^{l - r}$ $(a^{l - r} \approx 3 B^{b})$ . Consequently, the original problem defined on the infinite domain $x \in R^{1}$ is effectively reduced to a finite computational domain. This reduction implies that, for the considered scale ratios $a ∕ B^{b} ∕ l_{δ}$ , the displacement field must stabilize within the domain $| x | \leq a^{l - r} \approx 3 B^{b}$ , such that $⟨ u ⟩ (X) \approx$ const. (or, equivalently, the strain vanishes as indicated in equation (79)), $| x | > a^{l - r}$ .

A qualitative examination of Figure 5 provides convincing evidence supporting the validation of the three principal concepts introduced in this work: the AGIE, the new RVE concept, and the effective dataset $D^{p}$ . In this regard, Buryachenko [62] previously proposed a peridynamic, linear bond-based formulation—representing a particular case—of the AGIE in $R^{d}$ , although the term “AGIE” was not explicitly used at that time. Likewise, the new RVE concept had already been introduced implicitly in Buryachenko [62,63], even though the term “RVE” was not explicitly employed. Consequently, the domain $| x | \leq a^{l - r}$ effectively serves as an RVE for the effective displacement field $⟨ u ⟩ (x)$ , with its size governed by the scale ratios $a ∕ B^{b} ∕ l_{δ}$ .

Furthermore, the averaged displacement field $⟨ u ⟩ (x)$ shown in Figure 5 represents neither more nor less than a specific component of the newly defined effective dataset $D_{k}^{r}$ , evaluated for the particular loading $b (x)$ illustrated in the upper-right corner of Figure 5. At the time of its original introduction, this quantity was referred to as a “compressed generated dataset.” Importantly, the considered one-dimensional (1D) examples possess a distinct micromechanical priority over their potential two-dimensional (2D) and three-dimensional (3D) numerical counterparts. These 1D settings completely eliminate auxiliary side effects and ad hoc assumptions—such as hypotheses H2 and H3 (see Chapter 9 in Buryachenko [52])—as well as the influence of probability density functions of inclusion locations, which often obscure the clear identification of the individual roles played by the three key achievements: the AGIE, the new RVE concept, and the effective dataset $D$ . Both quantitative and qualitative validation of these concepts is thus clearly and independently demonstrated in Figure 5 and in the results reported in Buryachenko [62,63].

At the opposite end of the numerical validation spectrum lies the colored 2D example shown in Figure 4, which considers only a single deterministic microstructure and a single set of scale ratios $a ∕ l_{δ} ∕ B^{b} ∕ | Ω_{00} |$ . Owing to this highly specific configuration, it becomes difficult—even at a qualitative level—to disentangle the individual contributions of the three new concepts. Strictly speaking, only the new RVE concept is qualitatively manifested in Figure 4, as the results reflect a complex superposition of multiple interacting physical mechanisms, an issue that was not addressed in Silling et al. [61].

The ability of the RVE concept in Definition 2.2 to unify previously established methods and modeling approaches is of particular importance. This unifying capability not only confirms the internal consistency of the generalized definition but also demonstrates its flexibility and wide applicability. Such integrative potential is rare in micromechanics and represents a major step forward in establishing a more comprehensive and adaptable modeling paradigm for complex composite systems.

Due to the compact support of the applied body forces ${b_{k} (x)}_{k = 1}^{N}$ , the problem originally posed over the infinite domain $ℝ^{d}$ for both periodic and deterministic structure CMs can, in practice, be reduced to a finite computational domain corresponding to the RVE. This localization greatly simplifies the analysis by removing the necessity of simulating the full-space domain. As a result, common issues related to limited sample size or edge effects are naturally avoided. This type of problem can be tackled using a variety of computational mechanics methods, including the FEM, BEM, and others. Notably, Moulinec and Suquet [12,13] introduced an alternative and highly efficient technique based on FFTs—commonly referred to as “FFT-based methods”—to solve computational homogenization problems of the type given in L-S equation (24) with PBC (22) at UC $Ω_{00}$ . Their approach, known as the basic scheme, represented a major advancement at the time, as prior micromechanical simulations had been dominated by FEM-based techniques. Thanks to the advantages discussed in Section 1, FFT-based micromechanical methods now serve as competitive alternatives to FEM in a broad range of applications, including linear and nonlinear homogenization, local and nonlocal constitutive modeling, dislocation mechanics, and multiscale simulations. However, the solution of the ML-S equation (14) with PBC (19) at BFUC $Ω_{00}^{b}$ remains less explored. In this case, due to the periodicity requirement, the body-force UCs $Ω_{ij}^{b}$ (with $i, j = 0, \pm 1, \pm 2, \dots$ ) are assumed to be finite in extent, i.e., $Ω_{00}^{b} \neq ℝ^{d}$ . Nevertheless, due to the vanishing of the fields in the exterior region $\bar{RVE}$ , extending the domain $Ω_{00}^{b}$ does not affect the solution within the RVE—even in the limit where $Ω_{00}^{b} \to ℝ^{d}$ . This approach enables solving equation (14) using the same FFT-based method for both periodic structures ( $Ω_{00} \neq ℝ^{d}$ ) and deterministic structures ( $Ω_{00}^{b} = ℝ^{d}$ ). In the next section, we focus on FFT-based methods for solving equation (14) under the assumption that the body force $b (x)$ (5) has compact support.

4. DFT

In this section, we provide a brief overview of the notation and key properties of the DFT, tailored for direct use in the subsequent analysis.

We consider a UC $Ω_{00}^{b} = \prod_{α = 1}^{d} [0, l_{α}^{Ω}]$ , discretized using a uniform grid consisting of $N_{1} \times \dots \times N_{d}$ nodes along each spatial direction $α = 1, \dots, d$

\begin{matrix} x_{N}^{k} = k \cdot h \equiv {[k_{1} h_{1}, \dots, k_{d} h_{d}]}^{⊤}, k \in Z_{N}^{d} : = {m \in Z^{d} | 0 \leq m_{α} \leq N_{α} - 1} . \end{matrix}

(38)

hereafter, $ℤ_{N}^{d}$ denotes the set of d-tuples of integers corresponding to the grid indices, where each $N_{α}$ ( $α = 1, \dots, d$ ) is assumed to be even. We define the vector $N : = {(N_{1}, \dots, N_{d})}^{⊤}$ , with total number of grid points given by $| N | = \prod_{α = 1}^{d} N_{α}$ . The grid spacing in the α-th direction is $h_{α} = l_{α}^{Ω} ∕ N_{α}$ , and $e_{α}$ denotes the unit vector in the same direction. A function $f : R^{d} \to R^{d}$ is said to be $Ω_{00}^{b}$ -periodic if it satisfies:

\begin{matrix} f (x + \sum_{α} l_{α}^{Ω} k_{α} e_{α}) = f (x), \end{matrix}

(39)

for $x \in Ω_{00}$ and $k \in Z^{d}$ . In this section, we adopt the notation $Ω_{00}^{b} = \prod_{α} [0, l_{α}^{Ω}]$ , which is commonly used in the context of the DFT and its associated indexing. This contrasts with the earlier convention $Ω_{00}^{b} = \prod_{α} [- l_{α}^{Ω}, l_{α}^{Ω}]$ employed in Section 3.

If the integrable function f satisfies the periodicity condition (39), then its Fourier transform has a discrete frequency spectrum. The DFT provides a means to transition between the spatial domain $Ω_{00}^{b}$ and the corresponding frequency domain $F_{d}$ . The DFT $\hat{f} = F_{d} (f)$ and its inverse (iDFT) $f = F_{d}^{- 1} (\hat{f})$ for discrete periodic functions are defined as follows (see, e.g., [64], [65]):

\begin{matrix} {\hat{f}}_{N} (k) & = & \sum_{n \in Z_{N}^{d}} f_{N} (n) ω_{N}^{- kn}, \end{matrix}

(40)

\begin{matrix} f_{N} (n) & = & \frac{1}{| N |} \sum_{k \in Z_{N}^{d}} {\hat{f}}_{N} (k) ω_{N}^{kn} \end{matrix}

(41)

where $k, n \in ℤ_{N}^{d}$ represent multi-index vectors associated with the N-point discretization, and $f_{N} (n)$ and ${\hat{f}}_{N} (k)$ are the corresponding sampled sequences (which may be complex-valued). The Fourier kernel $ω_{N}^{k \cdot n}$ is given by:

\begin{matrix} ω_{N}^{kn} = \exp (2 πi \sum_{α} \frac{k_{α} n_{α}}{N_{α}}), \end{matrix}

(42)

with $i = \sqrt{- 1}$ and α indexing the spatial directions. To emphasize the link between discrete and continuous Fourier representations, we interpret $f_{N} (n)$ and ${\hat{f}}_{N} (k)$ as sampled versions of the continuous functions $f (x)$ and $\hat{f} (ζ)$ . This interpretation allows the introduction of explicit spatial sampling steps $Δ h$ and frequency intervals $Δ ζ$ , such that values may be written as $f (x_{N}^{n})$ and $\hat{f} (ζ_{N}^{k})$ with $ζ_{N}^{k} = k \cdot Δ ζ$ . This approach has the advantage of associating each sample with its actual position in physical or frequency space, rather than a purely integer-based index. For discussions on various DFT indexing conventions, their relationships, and reasons for choosing one over another, see, for example, [66]. Nonetheless, when formulating DFT-related theorems, the sampled coordinates are typically abstracted to integer indices n and k for simplicity (see [67]).

The key properties of the DFT that are particularly relevant for the subsequent analysis include:

\begin{matrix} {\hat{f}}_{N} (k + N) & = & {\hat{f}}_{N} (k), \end{matrix}

(43)

\begin{matrix} ℱ_{d} (a f + b g) & = & a \hat{f} + b \hat{g}, \end{matrix}

(44)

\begin{matrix} ℱ_{d} (\nabla f) & = & i k \hat{f}, \end{matrix}

(45)

\begin{matrix} ℱ_{d}^{- 1} ℱ_{d} (f) & = & f, \end{matrix}

(46)

\begin{matrix} ℱ_{d} {(f ⊙ g)}_{N} (k) & = & {\hat{f}}_{N} (k) \cdot {\hat{g}}_{N} (k), \end{matrix}

(47)

where

\begin{matrix} (f ⊙ g) (x_{N}^{n}) : = \sum_{j} f (x_{N}^{n− j}) \cdot g (x_{N}^{j}) . \end{matrix}

(48)

The operation described is a circular (or periodic) convolution. One important property of circular convolution is that it is commutative, meaning the order of the sequences does not affect the result:

\begin{matrix} (f ⊙ g) (x_{N}^{n}) = (g ⊙ f) (x_{N}^{n}) . \end{matrix}

(49)

Furthermore, scalar multiplication distributes over convolution, so a constant factor can be applied to either operand without changing the outcome:

\begin{matrix} a (f ⊙ g) (x_{N}^{n}) = ((a f) ⊙ g) (x_{N}^{n}) = (f ⊙ (a g)) (x_{N}^{n}) . \end{matrix}

(50)

Among the most fundamental properties of the DFT is the convolution theorem, which transforms the convolution operation in physical (real) space into a simple pointwise multiplication in the Fourier domain—greatly simplifying numerical computation. In the 1D case, the circular convolution of two N-periodic sequences f and g results in a new N-length sequence s defined as:

\begin{matrix} s (n) = f (n) ⊙ g (n) : = \sum_{m = 0}^{N - 1} f (m) g ((n - m) \mod N), \end{matrix}

(51)

where the modulo operation ensures periodicity by wrapping indices around the finite domain: if $m = m_{0} + lN$ with $m_{0} \in 0, 1, \dots, N - 1$ and $l \in ℤ$ , then $m mod N = m_{0}$ .

Since the DFT input and output (equations (40) and (41), respectively) are arrays of $N = | N |$ elements (real or complex numbers), they are similar to N–element vectors, and the DFT can be represented as the product of such a vector with an $N \times N$ matrix. Thus, the DFT (equation (40)) is a linear transformation of the discrete function $f_{N}$ and can be written in matrix form as

\begin{matrix} {\hat{f}}_{N} & = & F \cdot f_{N}, \end{matrix}

(52)

where F is a square Vandermonde matrix constructed from $ω_{N}^{- 1}$ (equation (42)), fully dense and invertible, with

\begin{matrix} F^{- 1} = \frac{1}{| N |} f (- ω^{- 1}) . \end{matrix}

In general, the DFT is a computationally costly algorithm, requiring $(O (N^{2}))$ operations, as each of the N Fourier components depends on NNN input values. However, this cost can be drastically reduced to $O (N \log_{2} N)$ by exploiting the symmetries of the transformation matrix F (equation (52)), as achieved in the so-called FFT algorithm, first proposed by Cooley and Tukey [68] and further developed in subsequent works (see, e.g., [24]). It should be noted that without the FFT, the DFT would be impractical for numerical methods; it is therefore the FFT, rather than the DFT itself, that gives its name to the numerical approaches briefly outlined in this section.

5. FFT methods for CMs subjected to BFCS loading

The original scheme by [12] of FFT method, based on fixed-point (Picard) iterations, can, in a retrospective sense, be considered like an ignition spark for a wide range of FFT-based methods for CMs. So, the implicit L-S equation (53) (see also equation (24)) arising in linear elastic micromechanics can be transformed into a simple multiplication operation in the Fourier space, as shown in equation (54)

\begin{matrix} ε (x) & = & ε^{w_{Γ}} + U^{(0)} * τ (x), \end{matrix}

(53)

\begin{matrix} \hat{ε} (k) & = & {\hat{U}}^{(0)} (k) \hat{τ} (k) (k \neq 0), \hat{ε} (0) = ε^{w_{Γ}} . \end{matrix}

(54)

Equation (53) is formulated under the PBCs given by (22) at UC $Ω_{00}$ . If the PBC are replaced with the BFCS loading (5) with the corresponding at PBC (19) at BFUC $Ω_{00}^{b}$ , the classical formulation given by equations (53) and (54) transforms accordingly into a new pair of ML-S equations that reflect this loading framework

\begin{matrix} ε (x) & = & ε^{b (0)} (x) + U^{(0)} * τ (x), \end{matrix}

(55)

\begin{matrix} \hat{ε} (k) & = & \hat{ε^{b (0)}} (k) + {\hat{U}}^{(0)} (k) \hat{τ} (k), (k \neq 0), \end{matrix}

(56)

and $\hat{ε} (0) = \hat{ε^{b (0)}} (0)$ . Under BFCS loading (5), the field periodicity within the UC $Ω_{00}^{b}$ , which holds in equation (53), is lost. As a result, $Ω_{00}^{b}$ is replaced by a larger mesocell $Ω_{00}^{b}$ that fully contains the RVE, i.e., $RVE \subset Ω_{00}^{b}$ (see Figure 3). In other words, equations (53) and (55) are formulated over different domains, $Ω_{00}$ and $Ω_{00}^{b}$ , respectively, with $Ω_{00}$ being a subset of $Ω_{00}^{b}$ (i.e., $Ω_{00} \subset Ω_{00}^{b}$ ). The size of the periodically distributed mesocells $Ω_{00}^{b}$ acts as a postprocessing learning parameter, chosen so that the strain field $ε (y)$ vanishes in the boundary layer region $y \in \bar{RVE} : = Ω_{00}^{b} ∖ RVE$ . This outer region imposes vanishing PBCs (22). Specifically, the vanishing of the strain field $ε (y)$ in $y \in Ω_{00}^{b} ∖$ RVE enables $ε (x)$ ( $x \in Ω_{00}^{b}$ ) to be regarded as periodic within an extended medium where the mesocell $Ω_{00}^{b}$ serves as the new periodicity cell. The resulting equations (55) and (56) can be solved iteratively using the same Picard (fixed-point) method as used for (53) and (54). This approach yields estimates of the effective dataset $D_{k}^{p}$ (35) for CMs with deterministic structure—a class for which FFT methods have not previously been applied. In this context, each specific translation $χ \in V_{x}$ (33) is treated as a deterministic structure. Consequently, the effective dataset $D_{k}^{p}$ (35) is obtained by averaging over these translations, using the statistical formulation given in (33) and (34).

Just as the Picard iteration scheme for solving equations (53) and (54) under PBCs (22) formed the basis for FFT-based homogenization methods for periodic composites [23 –25, 69], the analogous scheme for equations (55) and (56) under BFCS loading (5) may similarly drive the development of a new generation of FFT approaches. These would extend existing methods to handle both deterministic and periodic structures within a unified framework, accommodating non-periodic loading in $Ω_{00}$ while preserving the efficiency of Fourier-based solvers for $Ω_{00}^{b}$ . Deterministic microstructures could be embedded in mesocells with BFCS loading, enabling FFT analysis beyond classical PBC assumptions. For periodic structures, statistical averaging over translations (via effective dataset $D_{k}^{p}$ (35)) naturally integrates with this scheme. As with the original FFT method by Moulinec and Suquet [12], the proposed extension promises significant advances in modeling composites with high contrast, nonlinearity, or nonlocal effects.

So, by the use of the DFT’s properties (46) and (47), the discretized ML-S equation (55) can be presented in the form (see [16] for similar manipulation of L-S equation (53)

\begin{matrix} (ℐ - ℬ^{U}) e = e^{b (0)}, ℬ^{U} = ℱ_{d}^{- 1} \hat{U} ℱ_{d} ℒ_{1} \end{matrix}

(57)

where all the matrices exhibit a block-diagonal structure, e.g., $\hat{U} = [δ_{nm} {\hat{U}}_{αβ}^{(0) nm}],$ $ℐ = [δ_{nm} δ_{αβ}]$ , $ℒ_{1} = [δ_{nm} L_{1 | αβ}^{nm}]$ , ( $n, m \in Z_{N}^{d}$ , $α, β = 1, \dots, d$ ) whereas $e = (ε_{α}^{n})$ and $e^{b (0)} = ({(ε^{b (0)})}_{α}^{n})$ ( $n \in Z_{N}^{d}$ , $α = 1, \dots, d$ ). The cost of multiplication by $ℬ^{U}$ is governed by the forward DFT and iDFT, both efficiently performed in $O (| N | \log | N |)$ operations using the FFT. This makes system (57) well suited to iterative schemes. In particular, the original FFT-based basic scheme proposed in Moulinec and Suquet [12] expresses the solution of (57) as a Neumann series of the matrix inverse ${(ℐ - ℬ^{U})}^{- 1}$ :

\begin{matrix} e^{[j]} = \sum_{p = 0}^{j} {(ℬ^{U})}^{p} e^{b (0)} . \end{matrix}

(58)

which is similar in form to the fixed-point iterative method (Picard’s iterations) ( $ε^{[0]} (x) = ε^{b (0)}$ )

Algorithm 1. Modified basic scheme
$\begin{matrix} Data: L^{(0)}, L (x), tol, ε^{b (0)} (x) \\ Result: ε (x) \\ ε^{[0]} = ε^{b (0)} (x) \\ while: Δ_{σ} > tol do \\ τ^{[j]} (x) = L_{1} (x) ε^{[j]} (x) & (59) \end{matrix}$
$\begin{matrix} {\hat{τ}}^{[j]} (k) = ℱ_{d} (τ^{[j]} (x)), & (60) \end{matrix}$
$\begin{matrix} {\hat{ε}}^{[j + 1]} (k) = {\hat{U}}^{(0)} {(k)}^{L} {\hat{τ}}^{[j]} (k), & (61) \end{matrix}$
$\begin{matrix} ε^{[j + 1]} (x) = ℱ_{d}^{- 1} ({\hat{ε}}^{[j + 1]} (k)) + ε^{b (0)} (x), & (62) \end{matrix}$
end

where $k \neq 0$ and ${\hat{ε}}^{[j + 1]} (0) = {\hat{ε}}^{b (0)} (0)$ in equation (61). The convergence behavior of the modified basic scheme applied to equations (55) and (56) is the same as the original basic scheme used for equations (53) and (54). A common stopping criterion, first introduced by Moulinec and Suquet [13], is based on the relative equilibrium residual:

\begin{matrix} Δ_{σ} : = \frac{| \nabla \cdot σ |_{L_{2}}}{| Ω_{00}^{b} ⟨ σ ⟩ |} = \frac{[\sum_{ξ} | ξ \cdot σ |^{2}]^{1 ∕ 2}}{| \hat{σ} (0) |} < tol, \end{matrix}

(63)

where $∥ • ∥_{L_{2}}$ represent the $L_{2}$ -norm of the vector field and ||•|| is the Frobenius norm of the tensor. By Parseval’s theorem, this criterion is efficiently evaluated in Fourier space, making it well-suited for FFT-based iterative solvers.

Thus, the ML-S equation (55) can be solved in Fourier space similarly to a pointwise product, by using the convolution property (47), which reduces it to equation (61). In this form, the polarization τ (61) depends on the strain field $ε (x)$ and local stiffness $L (x)$ in real space $x \in R^{d}$ , requiring both forward (60) and inverse (62) DFTs. Brisard and Dormieu [15] (see also [23]) interpreted the discretization by Moulinec and Suquet [12] as similar to an under-integrated conforming Galerkin discretization of the Hashin–Shtrikman variational principle, with voxel-wise constant strain ansatz functions.

An alternative approach to enhance the convergence of FFT-based solvers, beyond the basic scheme, was introduced by Zeman et al. [16] and Vondřejc et al. [70] in the context of vector field homogenization for electrostatics. In this formulation, the integral L-S equation (53) is discretized using the trigonometric collocation method, wherein trigonometric polynomials serve as basis functions for interpolating field quantities in real space. The projection of the continuous equation onto this discrete function space yields a linear system of equations, where the nodal values of the strain field are treated as the primary unknowns. This discrete system is structurally similar to that arising in Galerkin-type methods and can be efficiently solved using Krylov-subspace algorithms, such as the CG or BiCG methods by Zeman et al. [16]. The resulting discrete counterpart of the modified L-S equation (55) takes the form:

\begin{matrix} ε (x) - ℱ_{d}^{- 1} {{\hat{U}}^{(0)} (k) : ℱ_{d} [L_{1} (x) : ε (x)]} = ε^{b (0)} (x), \end{matrix}

(64)

here, $ε (x)$ denotes the discrete (nodal) representation of the strain field. The left-hand side of equation (64) defines a linear operator acting on ε, which can be directly utilized within a Krylov-subspace solver. The structure of this linear operator is similar to that required by iterative methods such as the CG algorithm. The complete procedure for the Krylov-based scheme is summarized in Algorithm 2.

Algorithm 2. Modified Krylov-based scheme
$\begin{matrix} Data: L^{(0)}, L (x), tol, ε^{b (0)} (x) \\ Result: ε (x) \\ ε = 0 \\ Solve A (ε (x)) = ε^{b (0)} (x) by conjugate gradients with tol \\ A (ε (x)) = ε (x) - ℱ_{d}^{- 1} {{\hat{U}}^{(0)} (k) : ℱ_{d} [L_{1} (x) : ε (x)]} & (65) \end{matrix}$

Unlike traditional strain-based FFT schemes, displacement-based approaches solve directly for the displacement field, reducing memory usage and yielding a determinate system. This avoids the rank-deficiency issues of strain-based formulations and enables the use of preconditioners and alternative iterative solvers.

An alternative displacement-based FFT formulation for the ML-S equation (55)—extending a similar approach developed for equation (53) (see [22,23,25])—focuses on directly solving for the displacement field instead of the strain. By substituting the displacement decomposition (10) into the constitutive relation (2), the equilibrium equation (1) can be reformulated accordingly

\begin{matrix} \nabla \cdot [L (x) : \nabla^{s} u_{1} (x)] = - \nabla \cdot [L (x) : ε^{b (0)}] . \end{matrix}

(66)

This approach eliminates the need for a reference medium and supports both standard and staggered discretizations via discrete differential operators. The key idea is to derive a fully determined system in Fourier space by expressing the problem in terms of displacement fluctuations at each grid point. By removing symmetries associated with the real Fourier transform and excluding the zero-frequency component, the formulation yields a determinate linear system suitable for preconditioning. Rather than introducing a reference medium—as in equation (13)—to express the microstructural dependence $L (x)$ through an eigenstrain, equation (66) is directly transformed into Fourier space to compute spatial derivatives. Thus, applying the Fourier transform to the linear momentum conservation equation (66) yields

\begin{matrix} \hat{d} : ℱ_{d} [^{L} L (x) : ℱ_{d}^{- 1} (\hat{s} \cdot {\hat{u}}_{1})] = - \hat{d} : ℱ_{d} [L (x) : ε^{b (0)}], \end{matrix}

(67)

Using the Fourier derivative property (45), the operators in Fourier space are defined as

\begin{matrix} \hat{d} (k) : = {\hat{d}}_{pqr} (k) = i k_{r} δ_{pq}, \hat{s} (k) : = {\hat{s}}_{pqr} (k) = \frac{1}{2} (i k_{q} δ_{pr} + i k_{p} δ_{qr}), \end{matrix}

(68)

which correspond to the divergence and symmetric gradient operators, respectively, with dependence on the frequency vector k. Thus, both the differential form of the equilibrium equation (67) and its integral counterpart (55) reduce, under the DFT, to algebraic operations in Fourier space. When the Fourier transform is discretized as a DFT, equation (68) becomes a linear system over complex variables with $u_{1}$ as the primary unknown in Fourier space. This system can be efficiently solved using either direct or iterative solvers, potentially with preconditioning in Fourier space (see [22]).

It is important to emphasize that the linear system formulated in equation (67) is not expressed in terms of a traditional matrix of coefficients but rather through the action of a linear operator. This system can be formally written as

\begin{matrix} A (u_{1}) = b^{A}, \end{matrix}

(69)

where $A$ denotes a linear operator acting on the vector $u_{1}$ , and $b^{A}$ represents the known right-hand side derived from equation (67). Given the operator-based structure of equations (69), classical matrix-based solvers are not directly applicable. Among the various iterative methods suitable for operator-defined systems, the CG method is particularly effective and is employed in this study. As discussed by Barrett et al. [71] (see also [15]), the key requirement for implementing the CG method in this context is the ability to compute the action of the operator $A$ on a given input vector $u_{1}$ , without explicitly assembling a global stiffness matrix. This matrix-free formulation is highly advantageous in terms of memory efficiency and scalability. The evaluation of the operator action $A (u_{1})$ involves a structured sequence of operations that leverage the FFT framework. Naturally, the FFT is efficiently employed to perform the required DFTs and their inverses (iDFTs). For completeness, we provide below the pseudocode for the standard (unpreconditioned) CG algorithm. The method advances by iteratively constructing sequences of vectors: the approximations to the displacement fluctuation field $u_{1}^{[j]}$ (line (71) in Algorithm 3), the corresponding residuals $r^{[j]}$ (line (71)), and the conjugate (search) directions $p^{[j]}$ used to update both iterates and residuals. Each iteration involves two key scalar products (lines (70) and (73)) that yield step size parameters ensuring the satisfaction of orthogonality and conjugacy conditions between residuals and search directions. It is worth noting that the primary computational expense per iteration of the CG algorithm stems from the evaluation of the operator $A (u_{1})$ . This operation entails the same number of forward and inverse FFTs as a single iteration of the modified basic scheme. Consequently, the computational cost of applying $A$ is largely governed by the DFT operations $ℱ_{d}$ and $ℱ_{d}^{- 1}$ , both of which scale as $O (| N | \log | N |)$ due to the efficiency of FFT algorithms.

Algorithm 3. Modified CD scheme
$\begin{matrix} Data: L (x), tol, ε^{b (0)} (x) \\ Result: ε (x) \\ u_{1}^{[0]} = b^{A}, r^{[0]} = b^{A} - A (u_{1}^{[0]}), p^{[0]} = r^{[0]} \\ for j = 0, 1, 2, \dots, until convergence do \\ α^{[j]} : = (r^{[j] ⊤} r^{[j]}) ∕ (p^{[j] ⊤} A^{⊤} (p^{[j]})) & (70) \end{matrix}$
$\begin{matrix} u_{1}^{[j + 1]} : = u_{1}^{[j]} + α^{[j]} p^{[j]} & (71) \end{matrix}$
$\begin{matrix} r^{[j + 1]} : = r^{[j]} - α^{[j]} A (p^{[j]}) & (72) \end{matrix}$
$\begin{matrix} β^{[j]} : = (r^{[j + 1] ⊤} r^{[j + 1]}) ∕ (r^{[j] ⊤} r^{[j]}) & (73) \end{matrix}$
$\begin{matrix} p^{[j + 1]} : = r^{[j + 1]} + β^{[j]} p^{[j]} \\ end for \\ ε (x) = u_{1} (x) + ε^{b (0)} (x) . & (74) \end{matrix}$

In the classical setting of symmetric positive-definite (SPD) systems, these conditions ensure that each new iterate minimizes the error in an energy norm, and convergence is achieved in at most n steps for an n-dimensional problem. Interestingly, Vondřejc et al. [70] have shown that the CG algorithm remains effective within the local micromechanical framework even when the system matrix—such as in equation (57)—is not strictly symmetric. This empirical robustness opens the door for broader applications of CG in FFT-based schemes, including those where standard SPD assumptions do not hold. This insight is particularly relevant in the context of micromechanical formulations where the governing operators often arise from variational or weak forms that yield non-symmetric but still well-conditioned systems. Furthermore, within each CG iteration, a few inner iterations are performed to solve the ML-S formulation (see equations (55) and (56)). These steps closely follow the structure of the classical basic scheme of the L-S equation (53) and are embedded within the Krylov solver framework, similar to the nested iterative strategy outlined by Lahellec et al. [72]. In the context of the L-S equation (53) (and also of the ML-S equation (55)), numerical experiments confirm the robustness of this approach with respect to variations in internal parameters such as mesh resolution and material heterogeneity. Notably, the method demonstrates markedly improved convergence rates for problems characterized by high-contrast material coefficients—while maintaining a low per-iteration computational overhead due to the efficient use of FFT-based operator evaluations.

This establishes a formal correspondence between the DFTs of the classical L-S equation (54) and its modified counterpart (56), along with their respective equilibrium formulations (cf. [25] and Algorithm 3). This structural similarity enables the extension of numerous FFT-based algorithms—originally developed over the past three decades for micromechanics of periodic media under PBC (22), starting from the pioneering work by Moulinec and Suquet [12] (see also the comprehensive reviews [23,25])—to analogous formulations that accommodate BFCS loading (5) under the different PBCs (19). Since the fields vanish in the surrounding buffer zone $\bar{RVE} : = Ω_{00}^{b} ∖ RVE$ , see equation (36), the internal field distributions $σ (x)$ and $ε (x)$ within the domain $x \in RVE$ remain unaffected by further enlargement of the BFUC $Ω_{00}^{b}$ . Consequently, the FFT-based solution of the ML-S equation (55) is identical for CMs with either periodic ( $Ω_{00} \neq ℝ^{d}$ ) or deterministic ( $Ω_{00} = ℝ^{d}$ ) microstructures. In contrast, composite images obtained from micro-CT or scanning electron microscopy (SEM) typically contain several hundred or even thousands of inclusions (see [73,74]). Importantly, such CT-derived images can be interpreted as observational snapshots (or “windows of observation”) of deterministic composite structures. The observation window w is subjected to BCs (either (20₁), (20₂), or (22)). To enable numerical simulation using FEM or FFT methods, periodization of w is usually performed by rearranging inclusions near the boundary ∂w, which inevitably introduces boundary layer effects. In contrast, within the new proposed RVE framework, a similar rearrangement is applied only in regions $\bar{RVE}$ where the fields vanish. As a result, boundary layer effects are inherently and permanently eliminated.

Furthermore, upon projection of the equilibrium equations to a discrete functional space, the governing equations reduce to a linear system in which the unknown is the discrete, nodal representation of the displacement fluctuation field $u_{1}$ . This system can be efficiently solved using iterative Krylov subspace methods, such as the CG or BiCG algorithms, which are well-suited to large-scale, sparse systems. Specifically, the left-hand side of equation (69), $A (u_{1})$ , represents the application of a linear operator (rather than an explicit matrix) to the field $u_{1} (x)$ and thus lends itself naturally to matrix-free implementations. This operator can be directly utilized in Krylov solvers, including those provided by high-performance computing libraries such as PETSc [75] (see also [76]). Importantly, in the context of the L-S equation (53) (the ML-S equation (55) can be similarly considered), the studies by Zeman et al. [16] and Vondřejc et al. [70] have demonstrated that, in contrast to traditional acceleration techniques for FFT-based schemes (e.g., [77] and [14] methods), this Krylov-based method maintains the same per-iteration computational complexity as the original basic scheme. . That is, each iteration incurs a cost comparable to that of a single application of the ML-S operator in the basic FFT scheme, primarily driven by the required FFT and inverse FFT operations. Furthermore, the convergence behavior of the Krylov-based approach is shown to be largely independent of the choice of the auxiliary reference medium $L^{(0)}$ , eliminating the need for fine-tuning this parameter. The method offers significantly improved convergence rates over the modified basic scheme (59)–(61), particularly for moderate stiffness contrasts. However, its performance still degrades in scenarios involving high contrast in material stiffness. In such cases, the condition number of the system matrix associated with equation (55) becomes excessively large, rendering the system ill-posed and causing the iterative solver to converge slowly or stagnate. Despite this, the algorithm remains robust and efficient for a broad class of heterogeneous media and stands out for its simplicity, generality, and ability to leverage existing high-performance Krylov solvers.

In broad terms, FFT methods can be grouped into two categories. The first category originates from the L-S equations (53) and (54), while the second is derived from the ML-S equations (55) and (56). These groups differ only in their free terms, $ε^{w Γ}$ and $ε^{b (0)} (x)$ , respectively, while the kernels of the corresponding operators are identical. Consequently, the convergence proofs for the methods corresponding to the L-S equations (53) and (54) and the ML-S equations (55) and (56) coincide and are omitted for brevity. The most apparent distinction lies in the respective representative computational domains, UC $Ω_{00}$ and BFUC $Ω_{00}^{b}$ . A more essential difference concerns the intended outcomes of solving equations (53) and (54) versus equations (55) and (56). Solving (53) and (54) provides a complete solution, yielding both DNS results and the effective moduli $L^{*}$ (77). By contrast, solving (55) and (56) also serves the additional purpose of determining the RVE. If the RVE is contained within $Ω_{00}^{b}$ , the resulting DNS is then employed to construct $D^{p}$ . If the RVE is not fully contained within RVE $⊄ Ω_{00}^{b}$ , its size $Ω_{00}^{b}$ must be increased, and the solution of equations (55) and (56) repeated. It is also important to clearly distinguish between the BFCS $b (x)$ in (5) and the periodic $b (x)$ in (7). The body force $b (x)$ has compact support $ℬ^{b}$ when the identified RVE lies within $Ω_{00}^{b}$ . In this case, the solution $ε (x)$ in the isolated domain $Ω_{00}^{b}$ (associated with the BFCS $b (x)$ from (5)) is identical to the solution in the periodic medium with the representative BFUC $Ω_{00}^{b}$ (associated with the periodic body force $b (x)$ from (7)). Thus, if and only if RVE $\subset Ω_{00}^{b}$ (since the condition $ℬ^{b} \subset Ω_{00}^{b}$ is not sufficient), the notions of the BFCS $b (x)$ (5) and periodic body force $b (x)$ (7) coincide. This equivalence enables the generalization of the FFM methods to solve the ML-S equation (14)—rather than L-S equation (24) (see [12 –25]).

6. Effective elastic moduli and surrogate operators

Let us consider periodic CM with inclusion centers $Λ_{00}$ and PBC (22) at UC $Ω_{00}$ . By the use of the Gauss–Ostrogradsky theorem, we define the overall macrostress ${σ} = ⟨ σ ⟩^{Ω}$ , and the overall macrostrain ${ε} = ⟨ ε ⟩^{Ω}$ of the UC $Ω_{00}$

\begin{matrix} ⟨ σ ⟩^{Ω} & : = & | Ω_{00} |^{- 1} \int_{Ω_{00}} σ (x) d x = | Ω_{00} |^{- 1} \int_{Γ^{0}} t (s) \overset{S}{\otimes} s d s, \end{matrix}

(75)

\begin{matrix} ⟨ ε ⟩^{Ω} & : = & | Ω_{00} |^{- 1} \int_{Ω_{00}} ε (x) d x = | Ω_{00} |^{- 1} \int_{Γ^{0}} u (s) \overset{S}{\otimes} n (s) d s, \end{matrix}

(76)

in terms of the traction $t (s) : = σ (s) n (s)$ and the displacement $u (s)$ on the geometrical boundary of the UC $s \in Γ^{0}$ with the outward normal unit vectors $n (s)$ on $Γ^{0}$ . The effective stiffness $L^{*}$ is estimated as a proportionality factor between the UC’s averages of the stresses $⟨ σ ⟩^{Ω}$ and strains $⟨ ε ⟩^{Ω}$

⟨ σ ⟩^{Ω} = L^{*} ⟨ ε ⟩^{Ω} .

(77)

Estimating the involved macrovariables is performed by the micro-to-macro transition. Evaluations of effective moduli (77) by the FFT methods are well-developed directions in micromechanics (see Section 1 for references).

We now focus on BFCS, as described by equation (5). The foundational work of Silling and co-authors [38 –40] introduced ML approaches to develop surrogate nonlocal operators for CMs, using effective datasets generated via DNS. Specifically, their dataset structure is given by

\begin{matrix} D^{DNS} = {D_{k}^{DNS}}_{k = 1}^{N}, D_{k}^{DNS} = {u (b_{k}, x), b_{k} (x)}, \end{matrix}

(78)

where each realization corresponds to a different body force loading $b_{k} (x)$ . Their work primarily focused on 1D heterogeneous bars under wave-like loading at the boundary and oscillatory body forces. In contrast, the current study (see [62,63]) extends this idea to more general periodic and deterministic microstructures by replacing the DNS dataset $D^{DNS}$ (78) with datasets $D^{p}$ (35) and $D^{d}$ (37), corresponding to periodic and deterministic configurations, respectively. Here, the BFCS $b (x)$ serves both as a means of loading and as an input to infer nonlocal constitutive behavior. The datasets are compressed into surrogate forms ${\tilde{D}}^{p}$ and ${\tilde{D}}^{d}$ to reduce computational cost while preserving essential micromechanical detail. Each of these surrogate datasets is used to learn a corresponding nonlocal operator $ℒ_{γ}$ , characterized by a convolution-type integral operator:

\begin{matrix} ℒ_{γ} [⟨ u_{k} ⟩] (x) & = & Γ (x), \\ ℒ_{γ} [⟨ u_{k} ⟩] (x) & = & \int K_{γ} (| x - y |) (⟨ u_{k} ⟩ (y) - ⟨ u_{k} ⟩ (x)) d y, \end{matrix}

(79)

where $γ = b, σ, ϵi, σi$ indexes one of four different operator types—associated with either displacement, stress, or their localized forms—and $Γ_{k} (x)$ represents the corresponding averaged response: $- b_{k} (x)$ , $⟨ σ ⟩ (x)$ , ${⟨ ε ⟩}_{i} (x)$ , or ${⟨ σ ⟩}_{i} (x)$ . To determine the optimal kernel $K_{γ}^{*}$ for each case, the following minimization problem is solved:

K_{γ}^{*} = \arg \min_{K_{γ}} \sum_{k = 1}^{N} | | ℒ_{γ} [⟨ u_{k} (b_{k}) ⟩] (x) - Γ_{k} (x) ∥_{l_{2}}^{2} + ℛ (K_{γ}) .

(80)

where $R (K_{γ})$ is a regularization term (e.g., Tikhonov regularization) to stabilize the inverse problem. The kernel $K_{γ}$ is parameterized using Bernstein polynomial bases, and the optimization is performed using the Adam algorithm [78]. This formulation generalizes and strengthens earlier works [39,40,79 –81], providing a framework that combines physically motivated loading, rigorous micromechanical averaging, and data-driven learning of nonlocal operators across a wide class of complex materials.

The methods by You et al. [39] and You et al. [40] rely on uncompressed DNS datasets (78), which store full microscale displacement fields for each BFCS $b_{k} (x)$ , resulting in large data volumes. In contrast, the compressed effective datasets $D^{I}$ ( $I = p, d$ ) —whether periodic or deterministic—require no full-field DNS and instead use micromechanics-based averaging to extract effective quantities more efficiently. For a linearized homogeneous peridynamic medium under remote homogeneous BCs ((20₁) or (20₂)), classical peridynamics yields local moduli directly from the constitutive relation [82]. Similarly, for a surrogate homogeneous medium under the same BCs, the effective stiffness tensor is given by

\begin{matrix} L^{*} & = & \int K_{σ} (| x - y |) (y - x) d y, \end{matrix}

(81)

providing a compact expression for the homogenized elastic response. An analogy of the strain concentration (30) can be expressed as

⟨ ε ⟩_{i} (x) = \int K_{ϵi} (| x - y |) (y - x) d y ⟨ ε ⟩ .

(82)

The surrogate operators in equations (79) and (80) are fixed and limited to modeling linear responses. To overcome this limitation, nonlocal neural operators have been proposed to learn general mappings between function spaces, offering greater flexibility and adaptability [41,42]. Traditional artificial neural networks (ANNs), such as fully connected neural networks (FCNNs), define local nonlinear operators. For example, an L-layer FCNN $Ψ (x)$ : $R^{d_{x}} \to R^{d_{u}}$ maps input x to output u through successive transformations:

z^{l} (x) = A^{a} (w^{l} z^{l - 1} (x) + b^{l}), u (x) = w^{L} z^{L - 1} (x) + b^{L},

(83)

where $A^{a}$ is a nonlinear activation function (e.g., ReLU or tanh), and the learnable parameters are $ϑ = {w^{l}, b^{l}}_{l = 1}^{L}$ . Importantly, this type of operator is local, since the output at a point x depends only on the input at that same point. In contrast, nonlocal neural operators generalize this by incorporating integral terms to model spatial interactions and long-range dependencies. For instance, a typical nonlocal layer takes the form:

z^{l} (x) = A^{a} (w^{l} z^{l - 1} (x) + b^{l} + K^{l} (z^{l - 1}) (x)),

(84)

where $K^{l}$ represents an integral operator involving a learnable kernel $K^{l}$ , enabling the model to aggregate information from across the spatial domain. A wide range of architectures has been developed based on this idea, including Deep Operator Networks (DeepONet), PCA-Net, Graph Neural Operators, Fourier Neural Operators (FNO), and Laplace Neural Operators (LNO). These methods differ in how they implement the nonlocal interactions and the structure of the kernel $K^{l}$ . Comparative reviews and benchmarks of these neural operator frameworks are available in the literature [42,43,45,83].

PINNs [48,49,83 –88] incorporate governing physical laws—such as equation (2)—directly into the training process by embedding the residuals of these equations into the loss function. This ensures that the NN solutions remain consistent with the underlying physics. When PINNs are combined with neural operator frameworks [43,50,51], the resulting models can effectively learn complex material behavior, including nonlinear responses, microstructural heterogeneity, and nonlocal interactions, while maintaining strong generalization performance. Nevertheless, these approaches are typically limited to problems defined on finite computational domains and are not inherently suited for direct modeling of infinite media.

The procedure for generating surrogate models is outlined in Figure 6. Block 1 (DNS, see $D^{DNS}$ , equation (78)) was used by [79], [46], [39], [80], and [81]. In our CAMNN, we instead obtain compressed datasets $D^{p}$ (35) or $D^{d}$ (37), respectively (Block 2). These serve as inputs to Block 3 (Optimization), replacing the more data-intensive $D^{DNS}$ . This substitution significantly reduces effective dataset size and improves efficiency, requiring only minor adjustments between Blocks 2 and 3. Block 3 in Figure 6 then yields either a single surrogate operator (e.g., $K_{b}^{*}$ or $G$ , Block 4, see [38 –40]) or a set of surrogate models (e.g., equation (79) in Block 5). The only modification to the proposed approach is the replacement of Block 1 with Block 2 of the new compressed effective dataset; no adjustments to the existing Block 3 are required.

Figure 6.

The scheme of obtaining of surrogate model set.

Periodic structure CM with the PBC (22) are totally dominated in analyses by the most popular ANN methods such as are NO [42 –44]; PNO [46,47]; PINN [43,48 –50,78,89 –92]; and EINN [93 –95]. A key advantage of the BFCS loading (4) lies in its ability to generate effective datasets $D^{p}$ and $D^{d}$ while avoiding issues related to sample size, boundary layers, and edge effects. In addition, for periodic CMs, the CAM framework relies on solving problems (33) and (34) on the RVE, which may include several UCs, without requiring any specific PBCs (19) at UC interfaces are lost. Consequently, the direct application of both asymptotic homogenization [1,4] and computational homogenization methods [6 –11]—which inherently assume PBCs (19)—becomes questionable. Therefore, the revised RVE concept (see Definition 2.2), formulated under the general BFCS loading (5), is fundamental to the CAMNN approach—and to its extension to the ML and NN techniques discussed in Section 6—when applied to both periodic and deterministic composite structures.

Significant progress has been achieved in developing effective operators for both random [31 –33] and periodic [34 –37] structures (see also [52] for further references). Yet, all of these methods are confined to estimating predefined operators—most notably the ubiquitous fourth-order differential operator. If one dares to step outside this mold, for instance, by replacing a strain-type model with a displacement-type strongly nonlocal model, the entire micromechanical problem must be rebuilt and solved from the ground up. Moreover, selecting an alternative micromechanical method in this process poses an additional, nontrivial challenge. The present approach breaks free from this constraint. Here, all micromechanical investigation culminates in the determination of a new effective dataset $D^{p}$ (or $D^{d}$ ). From that moment, the computational model becomes entirely agnostic to both the composite’s microstructure and the numerical method (FEA, FFT, or otherwise) used to effective obtain the effective dataset. This effective dataset $D^{p}$ (or $D^{d}$ ) can then be seamlessly approximated by ML and NN techniques using any predetermined—or entirely “a priori” undefined—surrogate model. No further micromechanical computation is ever required. In essence, once the effective dataset is built, the heavy machinery of micromechanics can be switched off—the future analysis is instant.

7. Conclusion

To clarify the essence of the proposed approach, we briefly highlight the following key aspects: the novelty of the problem formulation, the presented solution, and the potential directions for its further development—both theoretical and practical. To the best of the author’s knowledge, the BFCS loading (5) has not previously been applied in micromechanics—whether for random or periodic structures—likely due to an initial underestimation of its practical significance. However, the principal motivation for introducing BFCS loading lies in the fundamentally new opportunity it provides: namely, its role as a training parameter for the estimation of unspecified surrogate nonlocal operators (see Section 6). This is achieved through several steps. First, the BFCS loading (5) is employed to define a new RVE concept (see Definition 2.2), via a novel general translation-based averaging procedure (33) and (34). The required DNS for each specific grid $Λ_{χ}$ in this averaging process can be carried out using any numerical method, particularly FFT-based techniques (see Section 5). The effective implementation of the new RVE concept relies critically on the accurate construction of the effective datasets $D^{p}$ (35) and $D^{d}$ (37). Only once these datasets have been precisely established can they be incorporated into ML and NN frameworks (see Section 6). Notably, the practical relevance of the ML-S equation (14) under compact support loading conditions remains limited—unless enhanced through ML and NN techniques. Without incorporating into ML and NN components, the BFCS (5) equation offers minimal practical utility, with the exception of certain special cases, such as the laser heating problem discussed in Buryachenko [96]. This explains why the ML-S equation (14)—despite being conceptually simpler than the original L-S formulation (24)—has historically received little attention. Only through a fully integrated framework—combining BFCS loading (5) with the new RVE concept, refined effective datasets $D^{r}$ or $D^{d}$ , and ML and NN tools—can the approach yield practical and predictive results. Symbolically, the proposed approach can be represented by the following sequence:

BFCF \to AGIE \to RVE \to D_{χk}^{p} \to D^{p} \to ML & NN .

(85)

The AGIE, the new RVE concept, and the effective dataset $D^{p}$ (obtained by modified FFT method, see Section 5), whose validity is established in Section 3, along with their coherent integration into the ML and NN framework, represent fundamentally new theoretical advances. To the best of the author’s knowledge, these concepts have not appeared in the literature by other authors—even at the idea level—either separately or in combination, in their present form. This unified methodology (85) (see [96] for details) enables the prediction of a broad class of “a priori” undefined surrogate operators—unlike the predefined effective nonlocal operators examined in [31 –37]. These operators encompass both macroscopic effective (nonlocal) properties and local concentration fields while removing dependence on sample size, BCs, and edge effects. These capabilities mark a substantial advancement and underscore the transformative potential of the ML-S equation when used in a data-driven context.

In standard practice, PINNs enforce physics by minimizing residuals of governing equations, BCs, and initial conditions for a prescribed material model, making them essentially model-driven rather than data-driven [44,48,49]. In contrast, in the proposed approach (85), these physical constraints are embedded a priori in the effective dataset $D^{p}$ , rendering explicit enforcement of governing- and boundary-condition losses unnecessary at the optimization stage. Thus, while both approaches are physics-informed, here the physics is enforced before training rather than during it. This ensures the highest attainable physical fidelity with controlled numerical accuracy and avoids boundary-layer effects, making the direct replacement $D^{DNS} \to D^{p}$ in PNO and HeteroPNO more natural than embedding $D^{p}$ into a PINN framework.

Regarding DNS on the specific grid of inclusions $Λ_{χ}$ , this paper extends FFT-based methods to the field of micromechanics for CMs with either periodic or deterministic structures, subjected to BFCS loading (5). The advantages of the FFT approach in the field of computational homogenization are (1) its very efficient numerical response, (2) the reduced memory allocation needs, (3) the possibility of using 2D/3D, and (4) the periodicity of the fields, which does not require the additional cost to impose periodicity in the FEM. The FFT approach is proposed by the use of the most popular tools and concepts exploited in the local elasticity of CMs (the modified L-S (14) equation-based approaches) and adapted to the case of BFCS loading (5).

Over the past three decades, following the seminal work of Moulinec and Suquet [12], numerous high-performance FFT-based algorithms have been developed to enhance convergence and accuracy in the analysis of microstructures with arbitrary phase contrast under finite deformations. These advancements have enabled the efficient simulation of nonlinear material behavior—such as plasticity, viscoplasticity, damage, fracture, and fatigue—on standard desktop computers, eliminating the need for high-performance computing clusters or supercomputers (see reviews in [25], [23], [24]). These FFT solvers have also been extended to address multiphysics problems involving coupled mechanical, thermal, electrical, magnetic, and pyroelectric effects. For instance, they have been employed in chemo-thermo-mechanical modeling of batteries, including phase-field approaches to damage, where multiple fields interact through coupled partial differential equations. A promising future direction would be to generalize these FFT methods—originally developed for PBC as in equation (19)—to accelerated solvers capable of handling BFCS-type loading conditions (equation(5)), such as those discussed in Section 5, with similarly improved convergence properties. It should be emphasized that the convergence behavior of the FFT schemes for the L-S equation (53) and for the ML-S equation (55) is identical, since the operator $U^{(0)}$ is the same in both equations; they differ only in their free terms.

An emerging direction for improving predictive capabilities involves the use of ML and NN techniques to construct surrogate operators. However, despite their power, these methods often overlook fundamental micromechanical principles such as scale effects, boundary layers, and the concept of the RVE. To address this limitation, the proposed CAMNN approach generates fundamentally new, compressed effective datasets for both periodic and deterministic composite microstructures. The novelty of CAMNN lies in replacing the dataset $D^{DNS}$ (78) with either of the effective datasets $D^{p}$ (35) or $D^{d}$ (37), while the Block 3 Optimization in Figure 6 remains unchanged. The newly introduced RVE concept (Definition 2.2) represents a significant departure from classical definitions. It is revolutionary in that it does not depend on the constitutive behavior of the individual phases or the specific form of the surrogate operator being predicted. Instead, it is grounded in the behavior of field concentration factors within the composite phases. This abstraction makes the definition more flexible and broadly applicable across various modeling approaches. Crucially, the generated effective datasets embed this generalized RVE concept as a core component, enabling seamless integration into existing ML and NN frameworks for the prediction of nonlocal surrogate operators. Incorporating this concept has the potential to significantly enhance both the reliability and generalizability of ML/NN-based models, particularly in complex systems where micromechanical accuracy is essential. By systematically eliminating size effects, boundary layers, and edge effects, the CAM approach ensures more robust and physically consistent predictions.

In summarizing, structuring, and generalizing the proposed approach, we highlight that it seamlessly integrates (three ingredients) a constitutive law agnostic RVE definition under BFCS loading, an FFT solver adapted to that loading, and ML/NN-based surrogates trained on compact, physics informed datasets. We begin with a novel RVE concept that is defined through BFCS loading (5)—his is entirely independent of any constitutive model, be it local, gradient-enhanced, peridynamic, or multiphysics (see classification of constitutive laws in Maugin [97]). This concept hinges on using field concentration factors within each phase (denoted in our notation as $D_{χ k}^{p}$ (35) and $D_{k}^{d_{j}}$ (37)) to characterize the microstructure. Crucially, it does not rely on a sample needing to be “large enough” or on specific phase behavior, making it universally applicable. Next, we drive this RVE into a numerical scheme based on the ML-S equations (55) and (56) solved via FFT (the second ingredient). What makes our method stand out is replacing the usual homogeneous BCs (20₁) or (20₂) with BFCS loading (5). This transforms the standard FFT solver for L-S equation (53) and (54) —originally developed for linear, high-contrast matrices or polycrystals—into a more general tool that handles BFCS-driven, possibly nonlinear or nonlocal composites. Finally, rather than training ML models on full-field DNS data $D^{DNS}$ (78), we train on the reduced effective datasets generated via the new RVE and FFT steps. By feeding in $D^{p}$ (35) or $D^{d}$ (37)—as opposed to massive DNS outputs—our CAMNN method builds surrogate operators with much smaller, noise-free inputs, while the optimization framework remains the same (the third ingredient, see Block 3 Optimization in Figure 6). Altogether, the innovation lies in how the BFCS-driven RVE + FFT yields compact, constitutive-independent microstructure descriptors that feed directly into ML/NN models (a set of surrogate operators, see Section 6). This creates a modular, robust pipeline that generalizes across local/nonlocal, linear/nonlinear, and multi-physics CMs.

Each of the three principal components of the proposed methodology addresses a remarkably broad spectrum of problems (see, e.g., the end of Section 5). These problems, despite their complexity or specificity, can often be effectively treated through relatively minor modifications or strategic reformulations of well-established problem classes. In the present work, the author has intentionally focused on a conceptual and schematic exposition of these reductions, highlighting the fundamental ideas together with their rigorous and explicitly demonstrated validation (see Figure 5 and the accompanying discussion), which provide a solid foundation for extending classical approaches. Readers who are interested in the technical nuances, computational strategies, and practical implementations of the proposed framework are encouraged to pursue these aspects independently. A comprehensive treatment of such applications—while undoubtedly important and promising—falls outside the scope of the current theoretical investigation and will be addressed in future work or in applied extensions of the present study.

Footnotes

Acknowledgements

The author acknowledges Dr Stewart A. Silling for the helpful comments, fruitful personal discussions, and encouragement suggestions. The author acknowledges the reviewers for the encouraging comments that initiated a significant correction of the manuscript. Permission to reproduce Figure 4(a) and (b) from Silling et al. [] was granted by Springer Nature (License Number 6044880317237).

ORCID iD

Valeriy A Buryachenko

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Fish

. Practical multiscaling. Chichester: John Wiley & Sons, 2014.

Ghosh

. Micromechanical analysis and multi-scale modeling using the Voronoi cell finite element method. Computational mechanics and applied analysis. Boca Raton, FL: CRC Press, 2011.

Zohdi

Wriggers

. Introduction to computational micromechanics. Berlin: Springer, 2008.

Bakhvalov

Panasenko

. Homogenisation: averaging processes in periodic media. Moscow: Nauka, 1984.

Buryachenko

. Peridynamic micromechanics of composites: a review. J Peridyn Nonlocal Model 2024; 6: 531–601.

Kouznetsova

Brekelmans

WAM

Baaijens

FPT

. An approach to micro-macro modeling of heterogeneous materials. Comput Mech 2001; 27: 37–48.

Matouš

Geers

MGD

Kouznetsova

, et al. A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. J Comput Phys 2017; 330: 192–220.

Terada

Kikuchi

. A class of general algorithms for multi-scale analyses of heterogeneous media. Comput Method Appl M 2001; 190: 5247–5464.

Geers

MGD

Kouznetsova

Brekelmans

WAM

. Multi-scale computational homogenization: trends and challenges. J Comput Appl Math 2010; 234: 2175–2182.

10.

Kanouté

Boso

Chaboche

, et al. Multiscale methods for composites: a review. Arch Comput Method E 2009; 16: 31–75.

11.

Raju

Tay

Tan

VBC

. A review of the FE2 method for composites. Multiscale Multidiscip Model Exp Des 2021; 4: 1–24.

12.

Moulinec

Suquet

. Fast numerical method for computing the linear and nonlinear properties of composites. C R Acad Sci Paris 1994; 318: 1417–1423.

13.

Moulinec

Suquet

. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Method Appl M 1998; 157: 69–94.

14.

Michel

Moulinec

Suquet

. A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int J Numer Meth Eng 2001; 52: 139–160.

15.

Brisard

Dormieux

. FFT-based methods for the mechanics of composites: a general variational framework. Comput Mater Sci 2010; 49: 663–671.

16.

Zeman

Vondřejc

Novák

, et al. Accelerating an FFT-based solver for numerical homogenization of periodic media by conjugate gradients. J Comput Phys 2010; 229: 8065–8071.

17.

Wicht

Schneider

Böhlke

. Anderson-accelerated polarization schemes for FFT-based computational homogenization. Int J Numer Meth Eng 2021; 122: 2287–2311.

18.

Wang

Fang

, et al. A novel FFT framework with coupled non-local elastic-plastic damage model for the thermomechanical failure analysis of UD-CF/PEEK composites. Compos Sci Technol 2024; 251: 110540.

19.

Wang

, et al. Artificial neural network-based homogenization model for predicting multiscale thermo-mechanical properties of woven composites. Int J Solids Struct 2024; 301: 112965.

20.

de Geus

TWJ

Vondřejc

Zeman

, et al. Finite strain FFT-based non-linear solvers made simple. Comput Method Appl M 2017; 318: 412–430.

21.

Zeman

de Geus

TWJ

Vondřejc

, et al. A finite element perspective on nonlinear FFT-based micromechanical simulations. Int J Numer Meth Eng 2017; 111: 903–926.

22.

Lucarini

Segurado

. DBFFT: a displacement-based FFT approach for non-linear homogenization of the mechanical behavior. Int J Eng Sci 2019; 144: 103131.

23.

Schneider

. A review of nonlinear FFT-based computational homogenization methods. Acta Mech 2021; 232: 2051–2100.

24.

Segurado

Lebensohn

Llorca

. Computational homogenization of polycrystals. Adv Appl Mech 2018; 51: 1–114.

25.

Lucarini

Upadhyay

Segurado

. FFT-based approaches in micromechanics: fundamentals, methods, and applications. Model Simul Mater Sc 2022; 30: 023002.

26.

Hill

. Elastic properties of reinforced solids: some theoretical principles. J Mech Phys Solids 1963; 11: 357–372.

27.

Bargmann

Klusemann

Markmann

, et al. Generation of 3D representative volume elements for heterogeneous materials: a review. Prog Mater Sci 2018; 96: 322–384.

28.

Kanit

Forest

Galliet

, et al. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 2003; 40: 3647–3679.

29.

Moumen

Kanit

Imad

. Numerical evaluation of the representative volume element for random composites. Eur J Mech A Solids 2021; 86: 104181.

30.

Ostoja-Starzewski

Kale

Karimi

, et al. Scaling to RVE in random media. Adv Appl Mech 2016; 49: 111–211.

31.

Drugan

. Micromechanics-based variational estimations for a higher-order nonlocal constitutive equation and optimal choice of effective moduli for elastic composites. J Mech Phys Solids 2000; 48: 1359–1387.

32.

Drugan

. Two exact micromechanics-based nonlocal constitutive equations for random linear elastic composite materials. J Mech Phys Solids 2003; 51: 1745–1772.

33.

Drugan

Willis

. A micromechanics-based nonlocal constitutive equation and estimates of representative volume elements for elastic composites. J Mech Phys Solids 1996; 44: 497–524.

34.

Ameen

Peerlings

RHJ

Geers

MGD

. A quantitative assessment of the scale separation limits of classical and higher-order asymptotic homogenization. Eur J Mech A Solids 2018; 71: 89–100.

35.

Kouznetsova

Geers

Brekelmans

. Size of a representative volume element in a second-order computational homogenization framework. Int J Multiscale Com 2004; 2: 575–598.

36.

Kouznetsova

Geers

Brekelmans

. Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy. Comput Method Appl M 2004; 193: 5525–5550.

37.

Smyshlyaev

Cherednichenko

. A rigorous derivation of strain gradient effects in the overall behavior of periodic heterogeneous media. J Mech Phys Solids 2000; 48: 1325–1357.

38.

Silling

. Propagation of a stress pulse in a heterogeneous elastic bar. Technical Report SAND2020-8197, Sandia National Laboratories, Albuquerque, NM, 5 December 2020.

39.

You

Silling

, et al. Data-driven learning of nonlocal models: from high-fidelity simulations to constitutive laws. arXiv [preprint] 2020. DOI: 10.48550/arXiv.2012.04157.

40.

You

Silling

, et al. Nonlocal operator learning for homogenized models: from high-fidelity simulations to constitutive laws. J Peridyn Nonlocal Model 2024; 6: 709–724.

41.

Kovachki

Azizzadenesheli

, et al. Neural operator: graph kernel network for partial differential equations. arXiv [preprint] 2003. DOI: 10.48550/arXiv.2003.03485.

42.

Lanthaler

Stuart

. Nonlocal and nonlinearity imply universality in operator learning. arXiv [preprint] 2024. DOI: 10.48550/arXiv.2304.13221.

43.

Goswami

Bora

, et al. Physics-informed neural operators. arXiv [preprint] 2022. DOI: 10.48550/arXiv.2111.03794.

44.

Chao

. Physics-Informed Neural Networks (PINN) for computational solid mechanics: numerical frameworks and applications. Thin-Walled Struct 2024; 205: 112495.

45.

Kumara

Yadav

. Deep learning algorithms for solving differential equations: a survey. Journal of Experimental and Theoretical Artificial Intelligence 2025; 37: 609–648.

46.

Jafarzadeh

Silling

Liu

, et al. Peridynamic neural operators: a data-driven nonlocal constitutive model for complex material responses. arXiv [preprint] 2024. DOI: 10.48550/arXiv.2401.06070.

47.

Jafarzadeh

Silling

Zhang

, et al. Heterogeneous peridynamic neural operators: discover biotissue constitutive law and microstructure from digital image correlation measurements. arXiv [preprint] 2024. DOI: 10.48550/arXiv.2403.18597.

48.

Raissi

Perdikaris

Karniadakis

. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys 2019; 378: 686–707.

49.

Karniadakis

Kevrekidis

, et al. Physics-informed machine learning. Nat Rev Phys 2021; 3: 422–440.

50.

Faroughi

Pawar

Fernandes

, et al. Physics-guided, physics-informed, and physics-encoded neural networks and operators in scientific computing: fluid and solid mechanics. J Comput Inf Sci 2024; 24: 040802.

51.

Wang

Yin

. Interpretable physics-encoded finite element network to handle concentration features and multi-material heterogeneity in hyperelasticity. Comput Method Appl M 2024; 431: 117268.

52.

Buryachenko

. Local and nonlocal micromechanics of heterogeneous materials. New York: Springer, 2022, p. 226.

53.

Mura

. Micromechanics of defects in solids. 2nd ed. Mechanics of Elastic and Inelastic Solids. Berlin: Springer, 1987.

54.

Buryachenko

. Micromechanics of heterogeneous materials. New York: Springer, 2007, p. 129.

55.

Milton

. The theory of composites. Cambridge: Cambridge University Press, 2022.

56.

Michel

Moulinec

Suquet

. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Method Appl M 1999; 172: 109–143.

57.

König

Carvajal-Gonzalz

Downs

, et al. Modelling and analysis of 3-D arrangements of particles by point process with examples of application to biological data obtained by confocal scanning light microscopy. J Microsc 1991; 161: 405–433.

58.

Ohser

Mücklich

. Statistical analysis of microstructures in material science. Chichester: John Wiley & Sons, 2000.

59.

Torquato

. Statistical description of microstructures. Annu Rev Mater Res 2002; 32: 77–111.

60.

Buryachenko

. Critical analyses of RVE concepts in local and peridynamic micromechanics of composites. J Peridyn Nonlocal Model 2024; 7: 9.

61.

Silling

Jafarzadeh

. Peridynamic models for random media found by coarse graining. J Peridyn Nonlocal Model 2024; 6: 654–683.

62.

Buryachenko

. Effective nonlocal behavior of peridynamic random structure composites subjected to body forces with compact support and related prospective problems. Math Mech Solids 2023; 28: 1401–1436.

63.

Buryachenko

. Effective displacements of peridynamic heterogeneous bar loaded by body force with compact support. Journal of Multiscale Computational Engineering 2023; 21: 27–42.

64.

Amidror

. Mastering the discrete Fourier transform in one, two or several dimensions. London: Springer-Verlag, 2013.

65.

Marks

. Handbook of Fourier analysis and its applications. New York: Oxford University Press, 2009.

66.

Briggs

Henson

. The DFT: an Owner’s manual for the discrete Fourier transform. Philadelphia, PA: SIAM, 1995.

67.

Brigham

. The fast Fourier transform and its applications. Hoboken, NJ: Prentice-Hall, 1988.

68.

Cooley

Tukey

. An algorithm for the machine calculation of complex Fourier series. Math Comput 1965; 19: 297–301.

69.

Buryachenko

. Fast Fourier transform in peridynamic micromechanics of composites. Math Mech Solids 2024; 29: 1844–1878.

70.

Vondřejc

Zeman

Marek

. Analysis of a fast Fourier transform-based method for modeling of heterogeneous materials. In: Lirkov

Margenov

Waśniewski

(eds) Large-scale scientific computing. Berlin: Springer, 2012, pp. 515–522.

71.

Barrett

Berry

Chan

, et al. Templates for the solution of linear systems: building blocks for iterative methods. 2nd ed. Philadelphia, PA: SIAM, 1994.

72.

Lahellec

Michel

Moulinec

, et al. Analysis of inhomogeneous materials at large strains using fast Fourier transforms. In: Miehe

(ed.) IUTAM symposium on computational mechanics of solid materials at large strains. Berlin: Springer, 2003, pp. 247–258.

73.

Bellens

Guerrero

Vandewalle

, et al. Machine learning in industrial X-ray computed tomography—a review. CIRP J Manuf Sci Technol 2024; 51: 324–341.

74.

Echlin

Lenthe

Pollock

. Three-dimensional sampling of material structure for property modeling and design. Integr Mater Manuf Innov 2014; 3: 1–14.

75.

Balay

Abhyankar

Adams

, et al. PETSc users manual 3.7. Technical report, Argonne National Laboratory, Argonne, IL, April 2016.

76.

Plimpton

Kohlmeyer

Coffman

, et al. fftMPI, a library for performing 2d and 3d FFTs in parallel. Technical report, Sandia National Laboratories, Albuquerque, NM, 2018.

77.

Eyre

Milton

. A fast numerical scheme for computing the response of composites using grid refinement. Eur Phys J Appl Phys 1999; 6: 41–47.

78.

Kingma

. Adam: a method for stochastic optimization. arXiv 2014. DOI: 10.48550/arXiv.1412.6980.

79.

Fan

D’Elia

, et al. Bayesian nonlocal operator regression: a data-driven learning framework of nonlocal models with uncertainty quantification. J Eng Mech 2023; 149: 04023049.

80.

You

Trask

, et al. Data-driven learning of robust nonlocal physics from high-fidelity synthetic data. Comput Method Appl M 2021; 374: 113553.

81.

You

Zhang

Ross

, et al. Learning deep implicit Fourier neural operators (IFNOs) with applications to heterogeneous material modeling. Comput Method Appl M 2022; 398: 115296.

82.

Silling

Zimmermann

Abeyaratne

. Deformation of a peridynamic bar. J Elasticity 2003; 73: 173–190.

83.

Daryakenari

Shen

, et al. State-space models are accurate and efficient neural operators for dynamical systems. arXiv 2024. DOI: 10.48550/arXiv.2409.03231.

84.

Cuomo

Schiano

Cola

, et al. Machine learning through physics-informed neural networks: where we are and what’s next. J Sci Comput 2022; 92: 88.

85.

Haghighata

Raissi

Moure

, et al. A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Comput Method Appl M 2021; 379: 113741.

86.

Harandi

Moeineddin

Kaliske

, et al. Mixed formulation of physics-informed neural networks for thermo-mechanically coupled systems and heterogeneous domains. Int J Numer Meth Eng 2024; 125: e7388.

87.

Kim

Lee

. A review of physics informed neural networks for multiscale analysis and inverse problems. Multiscale Sci Eng 2024; 6: 1–11.

88.

Ren

Lyu

. Mixed form-based physics-informed neural networks for performance evaluation of two-phase random materials. Eng Appl Artif Intel 2024; 127: 107250.

89.

Eghbalpoor

Sheidaei

. A peridynamic-informed deep learning model for brittle damage prediction. Theor Appl Fract Mec 2024; 131: 104457.

90.

Haghighat

Bekar

Madenci

, et al. A nonlocal physics-informed deep learning framework using the peridynamic differential operator. Comput Method Appl M 2021; 385: 114012.

91.

Ning

Cai

Dong

, et al. A peridynamic-informed neural network for continuum elastic displacement characterization. Comput Method Appl M 2023; 407: 115909.

92.

Paszke

Gross

Massa

, et al. Pytorch: an imperative style, high-performance deep learning library. Adv Neur In 2019; 32: 8024–8035.

93.

Zhou

. A nonlocal energy-informed neural network for peridynamic correspondence material models. Eng Anal Boundary Elem 2024; 160: 273–297.

94.

Zhou

. A nonlocal energy-informed neural network based on peridynamics for elastic solids with discontinuities. Comput Mech 2024; 73: 233–255.

95.

Zhou

. Transfer learning enhanced nonlocal energy-informed neural network for quasi-static fracture in rock-like materials. Comput Method Appl M 2024; 430: 117226.

96.

Buryachenko

. Unified micromechanics theory of composites. arXiv 2025. DOI: 10.48550/arXiv.2503.14529.

97.

Maugin

. Non-classical continuum mechanics: a dictionary. Singapore: Springer, 2017.